Sylvain Corlay
email: sylvain.corlay@gmail.com.

A fast nearest neighbor search algorithm based on vector quantization

Keywords: vector quantization, fast nearest neighbor search, Voronoi diagram, Delaunay triangulation, principal component analysis

In this article, we propose a new fast nearest neighbor search algorithm, based on vector quantization. Like many other branch and bound search algorithms [1, 10], a preprocessing recursively partitions the data set into disjointed subsets until the number of points in each part is small enough. In doing so, a search-tree data structure is built. This preliminary recursive data-set partition is based on the vector quantization of the empirical distribution of the initial data-set.

Unlike previously cited methods, this kind of partitions does not a priori allow to eliminate several brother nodes in the search tree with a single test. To overcome this difficulty, we propose an algorithm to reduce the number of tested brother nodes to a minimal list that we call "friend Voronoi cells". The complete description of the method requires a deeper insight into the properties of Delaunay triangulations and Voronoi diagrams.

Introduction

The problem of nearest neighbor search, also known as the post office problem [START_REF] Knuth | Art of Computer Programming[END_REF] has been widely investigated in the area of computational geometry. It is encountered for many applications, as pattern recognition and vector quantization.

The post-office problem has been solved near optimally for the case of low dimensions. Algorithms differ on their practical efficiency on real data sets. For large dimensions, most solutions have a complexity that grows exponentially with the dimension, or require a bigger query time than the obvious brute force algorithm. In fact, it has been noticed that, if n is the size of the data set and d is the dimensionality, the best choice becomes linear search when d > K log(n) for some positive constant K which depends on the chosen algorithm. This effect is known as the curse of dimensionality.

As concerns the application to (Voronoi) vector quantization, nearest neighbor projections are recognized to represent the critical part of most codebook optimization algorithms. In this case, the big amount of nearest neighbor searches we have to do shows that a preprocessing of the data-set will be profitable if it reduces the average query time. Still, in some particular cases, the codebook is chosen so that nearest neighbor search is performed easily, (as when dealing with product quantization). Moreover, non-Voronoi quantization methods can also be designed in order to simplify the projection procedure while preserving some important properties of optimal quantizers, as the stationarity in the quadratic case.

Let us also point out that a field recently emerged under the name of dual quantization [START_REF] Pagès | Intrinsic stationarity for vector quantization: Foundation of dual quantization[END_REF][START_REF] Pagès | Sharp rate for the dual quantization problem[END_REF]. In this context, the nearest neighbor search, i.e. the location of a point in a Voronoi partition, is replaced by the analogous procedure in the Delaunay triangulation. This localization procedure in Delaunay triangulations have been widely investigated in the practical viewpoint in terms of reduction of its computational complexity. We refer to Devillers, Pion and Teillaud for a review on this subject [START_REF] Devillers | Walking in a triangulation[END_REF].

Many nearest neighbor search algorithms rely on a recursive partitioning of the data-set resulting in a search-tree data structure [START_REF] Louis | Multidimensional binary search trees used for associative searching[END_REF][START_REF] Mcnames | A fast nearest-neighbor algorithm based on a principal axis search tree[END_REF]. The method proposed by McNames in [START_REF] Mcnames | A fast nearest-neighbor algorithm based on a principal axis search tree[END_REF] improved the classical Kd-tree algorithm [START_REF] Louis | Multidimensional binary search trees used for associative searching[END_REF] by taking advantage of the shape of the data-set thanks to principal component analysis. The "principal axis tree" algorithm performs much faster than the classical Kd-tree when the coordinates of the data-set are correlated and it seems to take better the growth of dimensionality.

In our case, the proposed algorithm uses vector quantization as a clustering method to perform this recursive partitioning and to take advantage of the geometry of the data-set. It is classical background that when dealing with empirical distributions, the quadratic vector quantization problem is equivalent to the reduction of the intraclass inertia of the related partition, and the specification of the classical Lloyd algorithm to this case turns out to be the k-means clustering algorithm.

We will see that one draw-back of this kind of partition is that, as other tree-based search algorithms, after determining the closest neighbor of a query in a leaf-node of the tree, the procedure has to move up to the parent node and determine whether brother nodes have to be explored or not. Unlike Kd-tree and "principal axis tree", our so-called "quantization tree" can't eliminate several brother nodes by with a single test. This is the motivation for the development of our friend node algorithm.

The paper is organized as follows. Section 1 is devoted to classical definitions and notations related to vector quantization. The link with the classification problem is pointed out. Section 2 recalls in mind some definitions of computational geometry which will be useful in the sequel. As both the fields of vector quantization and algorithmic geometry deal with the notion of Voronoi diagram, we apply ourselves to distinguish the corresponding definitions and notations. Section 3 makes a brief presentation of both the Kd-tree [START_REF] Louis | Multidimensional binary search trees used for associative searching[END_REF] and "principal axis tree" [START_REF] Mcnames | A fast nearest-neighbor algorithm based on a principal axis search tree[END_REF] algorithms. We deal with some optimizations that will be applicable with our quantization tree. Section 4 presents the "crude" quantization tree, i.e. without using any friend node algorithm. It is presented as the natural counterpart these two branch and bound algorithms with a quantization based partition of the data-set. Section 5 presents the friend node algorithm which was discussed above. Finally, the last section provides some performance comparisons between the different algorithms on various data-sets.

Vector quantization and Voronoi tessellations

We consider (Ω, A, P) a probability space and E a (real) finite dimensional Euclidean space. The principle of a random variable X taking its values in E is to approach X by a random variable Y taking a finite number of values in E.

Definition 1 (quantizer). In this surrounding, the discrete random variable Y is a quantizer of X. (1)

If X ∈ L p ,
Definition 2 (Voronoi partition). Consider N ∈ N * , Γ = {γ 1 , • • • , γ N } ⊂ E and let C = {C 1 , • • • , C N } be a Borel partition of E. C is a Voronoi partition associated with Γ if ∀i ∈ {1, • • • , N }, C i ⊂ {ξ ∈ E, |ξ -γ i | = min j∈{1,••• ,N } |ξ -γ j |}. If C = {C 1 , • • • , C N } is a Voronoi partition associated with Γ = {γ 1 , • • • , γ N }, it is clear that ∀i ∈ {1, • • • , N }, γ i ∈ C i . C i is called Voronoi slab associated with γ i in C and γ i is the center of the slab C i . We denote C i = slab C (γ i ). For every a ∈ Γ, W (a|Γ) is the closed subset of E defined by W (a|Γ) = ß y ∈ E, |y -a| = min γ∈Γ |y -γ| ™ . Definition 3 (Nearest neighbor projection). Consider Γ ⊂ E a finite subset of E. A nearest neighbor projection onto Γ is an application Proj Γ that satisfies ∀x ∈ E, x -Proj Γ (x) = min γ∈Γ |x -γ|.
To be more precise, if Proj Γ is a measurable nearest neighbor projection onto Γ, there exists a Voronoi

partition C = {C 1 , • • • , C N } associated to Γ such that Proj Γ = N i=1 γ i 1 Ci .
Proposition 1.1. Let X be an E-valued L p random variable, and Y taking its values in the settled point set

Γ = {y 1 , • • • , y N } ⊂ E where N ∈ N. Set " X Γ the random variable defined by " X Γ := Proj Γ (X)
where Proj Γ is a nearest neighbor projection on Γ, called a Voronoi Γ-quantizer of X. Then we clearly have X -"

X Γ ≤ |X -Y | a.s.. Hence X -" X Γ p ≤ X -Y p .
A consequence of this proposition is that solving the minimization problem (1) amounts to solving the simpler minimization problem

min { X -Proj Γ (X) p , Γ ⊂ E, card(Γ) ≤ N } . (2)
The quantity X -Proj Γ (X) p is called the mean L p -quantization error. When this minimum is reached, we refer to L p -optimal quantization. The problem of the existence of a minimum have been investigated for decades on its numerical and theoretical aspects in the finite dimensional case [START_REF] Pagès | A space quantization method for numerical integration[END_REF][START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF]. For every N ≥ 1, the L p -quantization error is Lipschitz-continuous and reaches a minimum. An N -tuple that achieves the minimum has pairwise distinct components, as soon as card(supp(P X )) ≥ N . This result stands in the general case of a random variable valued in a reflexive Banach space [START_REF] Luschgy | Functional quantization of Gaussian processes[END_REF]. If card(X(Ω)) is infinite, this minimum strictly decreases to 0 as N goes to infinity. The asymptotic rate of convergence, in the case of non singular distributions is ruled by the Zador theorem [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF]. A non-asymptotic upper bound for the quantization error is also available [START_REF] Luschgy | Functional quantization rate and mean regularity of processes with an application to Lévy processes[END_REF]. We now focus on the quadratic case (p = 2). For a L 2 random variable X, we now denote C N (X) the set of L 2 -optimal quantizers of X of level N and e N (X) the minimal quadratic distortion that can be achieved when approximating X by a quantizer of level

N . A quantizer Y of X is stationary (or self-consistent) if Y = E[X|Y ].
Proposition 1.2 (Stationarity of L 2 -optimal quantizers). A (quadratic) optimal quantizer is stationary.

The stationarity is a particularity of the quadratic case. In other L p cases, a similar property involving the notion of p-center occurs. A proof is available in [START_REF] Graf | Optimal quantizers for Radon random vectors in a Banach space[END_REF].

Definition 4 (Centroidal projection). Let C = {C 1 , • • • , C N } be a Borel partition of E. Let us define for 1 ≤ i ≤ N , G i = ß E[X|X ∈ C i ] if P[X ∈ C i ] = 0,
Var(X) = E î X -Proj C,X (X) 2 ó :=(1) + E î Proj C,X (X) -E[X] 2 ó :=(2)
.

The variance of the probability distribution X decomposes itself as the sum of the intraclass inertia (1) and the interclass inertia (2).

Proof:

Var(X) = E î X -Proj C,X (X) + Proj C,X (X) -E[X] 2 ó = E î X -Proj C,X (X) 2 ó =(1) + E î Proj C,X (X) -E[X] 2 ó =(2) + 2E X -Proj C,X (X), Proj C,X (X) -E[X] := (3)
.

Now (3) = 0 since Proj C,X (X) = E X Proj C,X (X) .

Backgrounds on theory of polytopes

Let E be a d dimensional vector space and E * its dual.

Definition 5 (k-flat).

A k-flat is a k-dimensional affine subspace E.

Definition 6 (convex polyhedron and convex polytope). A convex polyhedron is the intersection of a finite subset of closed halfspaces. If it is bounded, it is a convex polytope.

Definition 7 (cell).

A cell is the intersection of a finite set of flats and open halfspaces. And thus, equivalently, it is the relative interior of a convex polyhedron. If R ⊂ E, we denote cell(R) the relative interior of the convex hull of R.

Definition 8 (simplex). A simplex is cell(R) where R is a set of affinely independent points.

• A 2-dimensional simplex is the interior of a triangle.

• A 3-dimensional simplex is the interior of a tetrahedron.

Definition 9 (circumsphere). A circumsphere of a set R ⊂ E is a sphere S of E such that R ⊂ S. Proof: The proof is straightforward using the same approach as for the previous lemma. Lemma 2.3. If C is a non empty closed convex subset of E, distinct of E, then every point of the boundary ∂C belongs to a supporting hyperplane of C.

(x) = α. Consider v ∈ E such that ∀x ∈ E φ(x) = x|v . Consider a ∈ H ∩ C. We may assume that ∀x ∈ C φ(x) = x|v ≥ α.
Proof: a ∈ ∂C ⇒ ∀k ∈ N * , ∃x k ∈ B a, 1 k , x k / ∈ C. We denote y k = p C (x k ) the projection of x k on C, z k = x k -y k x k -y k .
Owing to the characterization of the projection on a closed convex subset, we have

∀z ∈ C, x k -p C (x k ), x k -z = |x k -p C (x k )| 2 - ≤0 x k -p C (x k ), z -p C (x k ) ≥ |x k -p C (x k )| 2 > 0 because x k / ∈ C.
Every vector z k lying on the unit sphere of E (which is compact), one can extract a subsequence of z φ(k) that converges to a vector v, with |v| = 1. As (x k ) 1≤k converges to a, by continuity of p C and of the scalar product, we have

∀z ∈ C, v, a -z = lim k→+∞ z φ(k) , x φ(k) -z ≥ 0.
In other words C is contained in the halfspace {z ∈ E, v, az ≥ 0}. Moreover, as a is in the corresponding hyperplane H, H is a supporting halfspace of C.

Definition 11 (face).

A face of a convex polyhedron P is the relative interior of the intersection of a hyperplane supporting P with the closure of P .

Proposition 2.4. Let P be a convex polyhedron, a face of P is a cell, and a face of a face of P is a face of P .

Definition 12 (k-face).

A k-face is a a face whose affine closure has dimension k.

Definition 13 (cell complex). A cell complex is a finite collection of pairwise disjoint cells so that the face of every cell is in the collection.

Definition 14 (opposite k-faces). Two distinct k-cells of a cell complex are opposite if they have a common (k -1)-face.

Definition 15 (triangulation). Let S be a finite point set of E. A triangulation T of S is a cell complex whose union is the convex hull of S and whose set of 0-cells is S.

Definition 15 is a non standard definition because cells are not required to be simplices. This formalism is due to Steven Fortune [START_REF] Fortune | Voronoi diagrams and Delaunay triangulations[END_REF].

Definition 16 (proper triangulation). A proper triangulation is a triangulation whose all cells are simplices.

Any triangulation can be completed to a proper triangulation by subdividing non simplicial cells.

Voronoi diagrams and Delaunay triangulations

Voronoi diagram

Let E be a d-dimensional Euclidean space, and S a finite subset of E. In the following, elements of S will be called sites.

Definition 17 (Voronoi cell). For a nonempty subset of S, R ⊂ S, the Voronoi cell of R, denoted V (R) is the set of all points in E that are equidistant from all sites in R, and closer to every site of R than to any site not in R.

Proposition 2.5.

• Clearly, is r ∈ S, V ({r}) is the set of all points strictly closer to r than to any other site. In particular, it is the interior of the Voronoi slab associated to r in S. (See the definition of a Voronoi slab in Section 1.)

• V (R) may be empty. Proposition 2.6 (Empty circumsphere property). For R ⊂ S, cell(R) is a Delaunay cell if and only if there is is a circumsphere of R that contains no site of S\R in its interior.

• Any point of E lies in V (R) for some R ⊂ S. Definition 18 (Voronoi diagram). The Voronoi diagram V is the collection of all nonempty Voronoi cells V (R) for R ⊂ S.
Proof: Such a circumsphere can be obtained with center an point in the Voronoi cell V (R). V is a cell complex that partitions E.

D is a triangulation of S.

3. V and D are linked with the following duality relation:

For R, R ′ ⊂ S, V (R) is a face of V (R ′ ) if and only if D(R ′ ) is a face of D(R).
4. V (R) is unbounded if and only if every site of R is on the boundary of the convex hull of S.

We refer to [START_REF] Fortune | Voronoi diagrams and Delaunay triangulations[END_REF] for a detailed proof.

Locality

Definition 21 (locally Delaunay). We consider two opposite d-cells cell(R) and cell(R ′ ) in a triangulation T with circumspheres C and C ′ . cell(R) and cell(R ′ ) are locally Delaunay if R ′ \R is outside of C. This is equivalent to R\R ′ outside of C ′ . A triangulation is locally Delaunay if every pair of opposite d-cells is locally Delaunay. Lemma 2.8 (Delaunay and locally Delaunay). A triangulation is Delaunay if and only if it is locally Delaunay.

We refer to [START_REF] Fortune | Voronoi diagrams and Delaunay triangulations[END_REF] for a detailed proof.

Definition 22 (General position). Let S be a nonempty finite set of sites in E. S is in general position if no d + 1 points of S are affinely dependent and if no d + 2 points of S lie on a common sphere.

Definition 23 (Incircle list). In the following, if S is a finite nonempty set of sites, D is a Delaunay triangulation of S and x ∈ E is a settle point, we call incircle list and denote ICL D (x) the set of d-cells of D whose circumsphere contains x.

If S is in general position, no Delaunay cell of S is degenerate. Every cell of the triangulation is a simplex and for any R ⊂ S, V (R) has dimension d + 1 -|R|.

Computing the Delaunay triangulation and the Voronoi diagram

Whereas the Voronoi diagram was defined before the Delaunay triangulation, it has been recognized that it is easier to devise algorithms in terms of Delaunay triangulation, especially because of the locality property 2.8.

A common data structure for Delaunay triangulations is a graph structure where each simplex is a "node". The node contains the indices of the d + 1 sites of the simplex and the pointers to the adjacent simplices. Null pointers are used when the simplices lie on the boundary of the triangulation. Cells of lower dimension are not directly represented in the graph structure. Another convenient convention is that the kth pointer stored in the node corresponds to the facet obtained by deleting the kth site in the node. Moreover the order is chosen so that the orientation of every simplex in the triangulation remains always positive.

Here, we present the principles of incremental algorithms for Delaunay triangulations. In this kind of algorithms, sites are added one by one, and the Delaunay triangulation is modified to include each new site. Many other algorithms have been designed for computing the Delaunay triangulation, especially in dimension 2. Moreover, computing the Delaunay triangulation of the Voronoi diagram in the onedimensional case simply amounts to sorting the data set. An advantage of incremental algorithms is that they are valid in any dimension. Moreover, for another purpose in the following, we will need a new algorithm (the friend node algorithm presented in Section 5) that requires a stage which is very similar to the insertion of a new point in the Delaunay triangulation. Hence we will focus here on incremental algorithms. Let S = (s 1 , • • • , s N ) be a nonempty finite set of sites of E of cardinal N . We define the sets

S k := (s 1 , • • • , s k ) for k ∈ {1, • • • , N }. Now,
for a settled i < N , let us consider D i the Delaunay triangulation of S i . We inspect the situation of s i+1 with respect to the Delaunay triangulation D i . From this analysis, the Delaunay triangulation will be modified locally to build a new Delaunay triangulation D i+1 of S i+1 . When all the sites of S will be processed, we will have the complete Delaunay triangulation D of S.

Three situations can occur, if S is in general position:

1. s i+1 lies in the interior convex hull of S i .

2. s i+1 does not lie in any circumsphere of any simplex of D i .

3. s i+1 lies outside of the convex hull of S i but belongs to a circumsphere of a simplex of D i .

(1) In the first situation, let denote S := ICL Di (s i+1 ) and F 1 , • • • , F p the external faces of S of any dimension k < d. We can show that the cell complex defined by

D i+1 := (D i \S) ∪ cell(F j , s i+1 ) j , 1 ≤ j ≤ p ∪ {s i+1 }
is the Delaunay triangulation associated to S i+1 . In a more general setting, we have the following property: Proposition 2.9 (star-shaped incircle list). Let S be a nonempty finite set of sites of E and x ∈ E that lies on the convex hull of S. Consider C the union of the d-cells of ICL D (x) and of all its faces. Then C is star-shaped from x, that is for any point p ∈ C, [x, p] ⊂ C.

(2) The second situation is the simplest. If F 1 , • • • , F p are the external faces of the triangulation D i (of any dimension k < d) that are visible from s i+1 . We can show that the cell complex defined by

D i+1 := D i ∪ cell(F j , s i+1 ) j , 1 ≤ j ≤ p ∪ {s i+1 }
is the Delaunay triangulation associated to S i+1 .

(3) In the third situation, if we denote S = ICL Di (s i+1 ) the set of elements of D i whose circumsphere contains s i+1 and F 1 , • • • , F p are the external faces) of this set which are not visible from s i+1 and F p+1 , • • • , F p+q are the external faces of D i that are not faces of elements of S and that are visible from x i+1 . We can show that the cell complex defined by

D i+1 := (D i \S) ∪ cell(F j , s i+1 ) j , 1 ≤ j ≤ p ∪ {s i+1 }
is the Delaunay triangulation associated to S i+1 .

The first triangulation D d+1 is made of a simple simplex defined by the d + 1 first inserted points.

One important modification of the incremental algorithm consists in inserting sites in a random order. Its expected running time is better than the worst case running time for the incremental algorithm.

The worst case complexity of computing the Delaunay triangulation of n points in a

d dimensional Euclidean space E is O n log(n) + n ⌈ d 2 ⌉ .
On the practical implementation The first step is the Localization. It consists in finding whether the new site x is in the convex hull of S or not, and if it is the case, in what Delaunay cell of the triangulation T S x lies. A survey on localization methods is available in [START_REF] Devillers | Walking in a triangulation[END_REF]. When x is inside of the convex hull of S, the localization procedure return the index of the the Delaunay cell where it lies. This corresponds to the situation (1). When x is outside of this convex hull, the localization returns a Null pointer. This corresponds to the situations (2) or (3).

The second step consists in finding the list of the Delaunay cells whose circumsphere contains x (the incircle list). In the situation (1), this list contains at least the Delaunay cell where x is located. Owing to the Proposition 2.9, we know that the union of these Delaunay cells is star-shaped so that it can be determined locally by testing connected cells in the graph structure presented above.

The last step consists in deleting the Delaunay cells of the incircle list and connecting the new site to the external faces of the incircle list or the visible faces of the convex hull of S depending on the situation (1), ( 2) or (3).

Classical examples of fast nearest neighbor search algorithms in low dimensions

Given a set of n points, {x 1 , • • • , x n } ⊂ E, the nearest neighbor problem is to find the point that is closest to a query point q ∈ E. Many algorithms have been proposed to avoid the large computational cost of the obvious brute force algorithm. When one has to perform a big amount of nearest neighbor searches, a preprocessing of the data set will be profitable if it reduces the average query time.

The problem is optimally solved in the case of dimension 1, where the best algorithm is, as a preprocessing to sort the data set by the unique coordinate of its points. (Approximative cost of O(n ln(n))). The search algorithm consists of a simple binary search whose cost is ln(n) ln(2) + O(1). In the case of low dimensions, most fast search algorithms still have an approximative preprocessing cost of O(n log(n)) and an average search cost in O(log(n)) in low dimension. The criterion of choice among them relies on

• their effective speed on real data sets,

• the required memory,

• the sensitivity of the speed to the dimensionality.

A first obvious optimization called partial distance search (P.D.S.) consists of a simple modification of the brute force search: during the calculation of the distance, if the partial sum of square differences exceeds the distance to the nearest neighbor found so far, the calculation is aborted. This almost always speeds up the nearest neighbor search procedure.

The Kd-tree algorithm

The Kd-tree algorithm is the archetype of the branch-and-bound nearest neighbor search tree. It is very popular because of its simplicity.

Building the tree:

• Every point of the data set is associated to the root node.

• The data set is being sorted by its first coordinate. Then it is divided in two subsets of cardinal

n 2 + 1 or n 2 .
• Each subset is associated to a child node of the root node.

• The process is repeated on each child node recursively using the coordinate axis in a cyclic order, until there are less than two points in each node.

Searching in the tree: Let q be the query point.

• The search procedure begins by searching in what child node q is (depending of its first coordinate).

• This child node is then searched, and the process is repeated recursively until a terminal node is reached.

• A trivial nearest neighbor search is performed in the terminal node. (Partial Distance Search optimization can be used.)

• The procedure moves up to the parent of the terminal node.

• If the distance d 2 between q and the hyperplane that splits the data set is smaller than the distance d min to the nearest neighbor found so far, the other child node is searched.

• The procedure continues its way back to the root node.

Complexity: Except in one dimension where the search complexity is logarithmic (it amounts to a binary search), the worst case of the Kd-tree corresponds to the case where every node of the tree is explored. Then the worst case complexity is time exponential. The distances to every point is computed. The complexity of the preprocessing is O(d × n log(n)).

The principal axis tree algorithm

The Principal Axis Tree (PAT) is a generalization of the Kd-tree proposed by McNames in [START_REF] Mcnames | A fast nearest-neighbor algorithm based on a principal axis search tree[END_REF]. Instead of using a coordinate axis to sort the data set, its principal axis is used at each step. Moreover, the number of child node in the tree can be greater than 2 at each generation.

Building the tree:

• Every point of the data set is associated to the root node.

• The data set is being sorted by its projection on its principal axis. Then it is partitioned in n c subsets whose cardinality is ö

n nc ù + 1 or ö n nc ù . d2 dmin q p
Figure 2: K-d tree elimination condition: if the distance d 2 between the query point q and the brother node is smaller than the distance d min to the nearest neighbor found so far, say p, the brother node has to be explored.

• Each subset is associated to a child node of the root node.

• The process is repeated on each child node recursively until there are less than n c points in each node.

• At each step, the principal axis, and maximal and minimal values of subset's projection on the principal axis are kept in memory.

Optimizing the elimination condition: We refer here to Figure 3. We can improve the lower bound to the points that belong to child nodes of brother nodes. For any point q in region 1 and x in region 2, we have

d 2 (q, x) ≥ d 2 q2 + d 2 2x
. This result is then used again to get a lower bound to points in region 3, and 4 and so on.

d 2 2x ≥ d 2 23 ∀x ∈ Region 3, d 2 (q, x) ≥ d 2 q2 + d 2 23 + d 2 34 ∀x ∈ Region 4.
Searching in the tree: Let q be the query point.

• The search process begins by searching in which child node q is (by computing its projection on principal axis).

• This child node is then searched, and the process is repeated recursively until a terminal node is reached.

• A partial distance search is then performed in the terminal node.

• The procedure moves up to the parent of the terminal node.

• The elimination condition is checked to decide if brother nodes have to be searched or not.

• The procedure continues its way back to the root node.

Choice of parameter n c : For normal or uniform random data sets (and distribution of query points), best overall performances are obtained with n c = 7 (independently from dimensionality for d < 10).

(The same optimal value is obtained by McNames in [START_REF] Mcnames | A fast nearest-neighbor algorithm based on a principal axis search tree[END_REF].) In the case where the data set is an optimal quantizer of those distributions, best performance is obtained with n c = 13.

Complexity: Space storage is O(n). Except in the one-dimensional setting where the search complexity is logarithmic (it comes to a binary search), the worst case of the Kd-tree corresponds to the case where every node of the tree is explored. Then the worst case complexity is time exponential (2 n comparisons of coordinates). n distances are computed. The complexity of the preprocessing is O(d × n log(n)).

Algorithm performance: On a 5000 points Gaussian data set in R 2 , the depth of the tree is 4.

• 27 (partial) distances,

• 15 scalar products,

• 9 binary searches are performed in average. Why using this space partitioning ? The idea is that good empirical performance of PAT are due to the fact that it takes advantage of the shape of the data set. Yet obviously when both query point distribution and data sets lie on a smaller dimension (k < d) subspace of E, one retrieves the same complexity as when using the same algorithm on a k dimensional space. This intrinsic dimension is often less than the spatial dimension of the space. In a more general setting, PAT takes advantage of high correlations in the data set coordinates. However if one uses the same number of child nodes n c in Kd-tree and PAT tree, we see that

• Preprocessing time is longer for PAT than for Kd-tree.

• The first traversal of the tree to a terminal node is more costly (projections have to be computed).

But PAT is still faster because its geometrical partition of the space fits the data set in a more relevant way. To be precise, it happens less often than one has to search a brother node with PAT than with Kd-tree.

In [START_REF] Wim | An efficient branch and bound seach algorithm for computing k nearest neighbors in a multidimensional vector space[END_REF], the same space decomposition was proposed for the nearest neighbor search problem (but using the only 2 child node at each generation). They justify the use of this decomposition using a heuristic criterion, according to which the best possible decomposition of the data-set into two subsets for branch and bound nearest neighbor search is to split the data set with respect to its projection on the principal axis.

A new quantization based tree algorithm

As we have seen in previous sections, a good space decomposition that fits to the data distribution may lead to a faster branch and bound nearest neighbor search algorithm, if less brother nodes have to be explored. The traversal of the tree can be a little more expensive if it is compensated by the gain due to the fact that less nodes are explored.

Principal component analysis and optimal quantization are two types of projection of a probability distribution. Similar inertia decompositions hold in the quadratic case (Huyghens lemma).

PAT is based on a recursive space decomposition based on the principal component analysis of the underlying data set. The initial idea here is to design a branch and bound algorithm based on a recursive quantization of the empirical distribution of the underlying data set.

The crude quantization tree algorithm

Building the tree:

• Every point of the data set is associated to the root node.

• The data set is being partitioned into n c subsets corresponding to the Voronoi cells of an optimized quantizer of the empirical distribution of the data set.

• Each subset is associated to a child node of the root node.

• The process is repeated on each child node recursively until there are less than a certain number of points in each node.

Some other computations are done during the preprocessing that will be detailed further on.

Remark. One notices that the resulting search tree is not balanced and may have some longer branches.

Searching in the tree: Let q be the query point.

• By performing trivial nearest neighbor researches in the node's quantizer the search algorithm traverses the tree to a terminal node where a trivial partial distance search is performed.

• The procedure moves up to the parent of the terminal node.

• The elimination condition, (developed further on) is checked to decide whether brother nodes have to be searched or not.

• The procedure continues its way back to the root node.

Consistency of the space decomposition:

Implementing only the way down to the terminal node (with n c = 7 in both principal axis tree and quantization tree), we naturally do not obtain always the index of the nearest neighbor. But we have noticed that the result is more often the right one with the quantization tree than with the principal axis tree.

For instance, in dimension 2, on a 5000 points Gaussian data set, on a million Gaussian query points, we notices:

• 56 percent of false results with PAT.

• 16 percent of false results with the quantization tree.

Similar results are obtained with other values of the parameters and other data set distributions. This empirical test makes us reasonably optimistic about the performance of a branch and bound tree based on this decomposition.

Still, the cost of the way through the search tree is more expansive with the quantization tree (as described above).

• For the "quantization tree", we have to perform trivial nearest neighbor search to find the right child node.

• For "principal axis tree", we only compute a projection and perform a binary search.

Moreover, it was proved in [START_REF] Tarpey | Principal points and self-consistent points of elliptical distributions[END_REF] that in the case of Gaussian distributions, the affine subspace spanned by stationary quantizers correspond to the first principal components of the considered Gaussian distribution. (This result, extended to the infinite dimensional case in [START_REF] Luschgy | Functional quantization of Gaussian processes[END_REF] allows to efficiently compute optimal quadratic quantizers of bi-measurable Gaussian processes.) Hence, in this case, this shows that the quantization tree with two branches at each generation is related to the principal axis tree.

First elimination condition If the center of the Voronoi cell corresponding to the current node is A, the first rough method to decide whether a brother node with center B has to be explored or not is compute the distance d 2 of the query point Q to the Leibniz halfspace H(B, A). Then the node corresponding to point B is explored if d 2 is smaller than the distance to the nearest neighbor found so far, d 1 . We

have d 2 = AB 2 -AQ cos α and QB 2 = QA 2 + AB 2 -2AQAB cos α so that ⇒ cos α = QA 2 +AB 2 -QB 2 2AQAB . This yields d 2 = QB 2 -QA 2 2AB
. Hence, the computation of the distance to the Leibniz halfspace requires one subtractions QA 2 -QB 2 , (QA 2 and QB 2 can be computed during the search in the quantizer in the parent node), and one multiplication by 1 2AB . ( 1 2AB can be computed during the preprocessing.) Then, it is clear that the nearest brother node correspond to the second nearest neighbor in the quantizer, and the second nearest to the third nearest neighbor, and so on. Hence, brother nodes have to be explored in the order defined by the distances of its centers the query point.

We can also use the same optimization of the lower bound proposed by McNames in [START_REF] Mcnames | A fast nearest-neighbor algorithm based on a principal axis search tree[END_REF] and presented in Section 3.2. Referring to Figure 4, the lower bounds d i are recursively incremented when exploring brother nodes. Performance of this first quantization tree algorithm. This first algorithm has been implemented and its empirical performances has been compared to the two previously exposed PAT and Kd-tree in terms of empirical performances. Intermediate performances between our implementations of Kd-tree and PAT were obtained in small dimensions. Although, as we will see further in empirical tests, it seems to take better the increase of dimensionality. The preprocessing time, that requires small quantizer computations is also more costly than both PAT and Kd-tree.

Optimizations for the quantization tree

To reduce the average query time, we are now proposing a new optimization procedure which reduces the number of brother nodes to be checked. Let us consider the Voronoi diagram plotted in Figure 5. In this figure, we obviously know that when the query point is in a cell A, its nearest neighbor cannot be in cell B, because cell B is "hidden" by closer cells. One has to give a precise mathematical sense to "hidden" in this sentence. However, in the quantization tree as it has been described, the distance of query point to H(a, b) has to be computed.

A first idea is to compute for each 1 ≤ i ≤ n c a list of "friends" among brother nodes in which the nearest neighbor can be when q is in cell i.

This list has to be large enough to ensure that it contains the nearest neighbor but as small as possible in order to reduce the computations of elimination conditions.

As concerns the choice of the parameter n c , we have to take in consideration that increasing n c makes the depth of the tree smaller but also makes the nearest neighbor search slower for each generation of the search tree.

How can we obtain a friend Voronoi slabs list? The first observation about obtaining such a friend list is that it is not a simple problem. Indeed, this list is a priori not reduced to slabs whose corresponding Voronoi cells are adjacent in the Voronoi diagram. Moreover, in some cases, the minimal friend list can be quiet large. So is the case for unbounded Voronoi slabs for example. Figure 6: In these cases, the nearest neighbor of the query point q may be p although p is not in an adjacent Voronoi cell.

5 Some optimizations for the quantization tree algorithm

In Section 2.1, basic definitions about Voronoi diagrams and Delaunay triangulations that are prerequisites to this section ahev been recalled.

Remark (Voronoi slabs and Voronoi cells). From their respective definitions, one can easily deduce the following properties:

• Let S ⊂ E be a finite set of sites, let C be an associated Voronoi partition and consider s ∈ S.

Then it is clear that V ({s}) =

• ¸ slab C (s).

• The points of the Voronoi cells V (R) with R ⊂ S and card R > 1 belong to the boundaries of Voronoi slabs.

• As a consequence, for s ∈ S, as the boundary V ({s}) is constituted with its faces of lower dimensions, previous remark yields V ({s}) = slab(s) and δ slab S (s) = ∂V S ({s}).

Notations: In the following of this section, if S ⊂ E is a finite set of sites in E, one will denote T S the Delaunay triangulation of S, DG S the Delaunay graph of S, V S its Voronoi diagram. For R ⊂ S, For every point p of the underlying codebook,

• Compute s = Proj S (p) and P I S (p).

• Then for every point s ′ ∈ P I S (p), insert s in the set of friends of node s ′ .

This method gives a first algorithm to compute friend list. Still, when the data set is large, it is very expensive because one has to deal with all the points of the data set.

In fact it is possible to compute an acceptable friend list thanks to the same result 5.3 without using the points of the underlying data set.

Fast friend node algorithm: In this section, another method to compute friend node lists is devised which does not need to deal with the complete underlying data set but only the underlying codebook.

When keeping the same notations, the principle of the method is to compute for every slab S (s), s ∈ S ,of the Voronoi partition C S , the set U P I S (s) := p∈slabS (s) P I S (p). It is the union of all the pseudo-insertions of points of slab S (s). If one is able to compute this set, the resulting friend nodes algorithm simply writes: For every point s ∈ S,

• Compute U P I S (s).

• Then for every point s ′ ∈ U P I S (s), insert s in the set of friends of node s ′ .

The question is: how can we compute U P I S (s)? Remark. Let us recall that, thanks to Proposition 2.5, (δ slab S (s) = ∂V S ({s}).

Proof: Consider x ∈ slab S (s) such as s ′ ∈ P I S (x). Let us define x * , such that {x * } = [x, s ′ ] ∩ ∂V S (s).

• One has H(x * , s ′ ) ⊃ H(x, s ′ ). So V S ({s ′ }) ∩ H(x * , s ′ ) ⊃ V S ({s ′ }) ∩ H(x, s ′ ), hence V S ({s ′ }) ∩ H(x, s ′ ) = ∅ ⇒ V S ({s ′ }) ∩ H(x * , s ′ ) = ∅ that is equivalent to s ′ ∈ P I(x *
) thanks to the Lemma 5.3.

• Finally, ∀x ∈ slab S (s), ∀s ′ ∈ P I S (x), ∃x * ∈ δ slab S (s) such that s ′ ∈ P I S (x * ).

Remark. As there are not a finite number of sites on the boundaries, this does not give an effective method for computing U P I S (s) yet.

As seen in Section 2.1, computing the set P I S (x) corresponds almost to the same algorithm as the insertion procedure in an incremental triangulation algorithm, that is:

• Localization of x in the triangulation,

• Computation of the set ICL(x),

• U I S (x) is the set of points that belong to a cell of ICL(x) plus, if x is outside the convex hull of S, the points of the external faces of T S that are visible from x.

Lemma 5.6. Let S be a non empty finite set of sites in E. We consider the circumsphere C of Delaunay d-cell of the Delaunay triangulation T S . We denote c its center and r its radius. Let s be a site of S.

If V S ({s}) ⊂ C = ∅ then c + r |s-c| (s -c) ∈ V S ({s}).
The proof is straightforward. This leads to an algorithm to compute sets (U P I S (s)) s∈S .

• For every Delaunay d-cell D of T S -Compute the center c and radius r of its circumsphere.

-For every site s ∈ S that is not in D, compute p := c + r s-c |s-c| ∈ V S ({s}), and check if the site s is the nearest neighbor of p in S. If so is the case, then the points of the Delaunay d-cells D belong to U P I S (s).

• Then deal with unbounded Voronoi cells:

-For every external face F of the Delaunay triangulation, compute a normal vector u F directed toward the exterior of the convex hull of S.

-For two distinct external faces F 1 and F 2 of the Delaunay triangulation, if u F1 , u F2 > 0 then for every (s 1 , s 2 ) ∈ F 1 × F 2 , s 1 ∈ U P I S (s 2 ) and s 2 ∈ U P I S (s 1 ).

In Figure 8, we present some friend Voronoi lists in the 2-dimensional case.

Test with real data sets

To perform the following tests, the quantization tree algorithm and the friend-node optimization have been implemented in C++. Because of the additional feature related to computational geometry that we needed, as the pseudo-insertion procedure, we had to implement a Delaunay triangulation. All the figures presented in this article were generated with this implementation of the Voronoi diagram with which we performed the following tests.

Tests on Gaussian and uniform data sets

In Tables 9, 10 and 11, we report the execution time for 10 millions nearest neighbor queries on datasets of size 5000 generated with independent Gaussian pseudo random variables and with a uniform distribution on the hypercube. The best overall performances were obtained with n c = 35 children by node for the quantization tree. The tests were performed with an Intel Pentium Dual CPU at 2GHz. We noticed that in dimension d = 2 and d = 3, we had intermediate performances between the "principal axis tree" and the Kd-tree algorithms. In dimension 4, the performance of the "principal axis tree" and the "quantization tree" are close one to each other. Finally, it seems that the quantization tree has a better behaviour in dimensions greater than 5 where it significantly outperforms the two other implemented methods. 6.47s 10.33s 15.91s 32.62s Principal axis tree 0.74s 1.52s 2.81s 6.71s 16.53s 28.03s 47.53s Kd-tree 1.54s 2.82s 5.46s 10.64s 18.50s 31.60s 55.71s Remark (Computational cost or the preprocessing for the friend cell algorithm). An important fact that we have experienced is that, in higher dimensions, the friend cells list becomes bigger and there is no more competitive advantage in using it in dimension higher than 7 (when having less than 30 branches per generation in the quantization tree). Moreover, as it requires to compute Delaunay triangulations during the preprocessing, whose complexity exponentially increases with the dimension, the computational cost of the friend cell preprocessing makes it useless in higher dimensions.

  the quantization error is the L p norm of |X -Y |, where | • | denotes the Euclidean norm on E. The minimization of this error yields the following minimization problem min{ X -Y p , Y : Ω → E measurable , card(Y (Ω)) ≤ N }.

  If a does not belong to the boundary of C, ∃ε ≥ 0, B(a, ε) ⊂ C so for any λ > 0 small enough, aλv ∈ C and α ≤ φ(aλv) = a|vλ v 2 < a|v = α which yields a contradiction. Consequently a ∈ ∂C.

Corollary 2 . 2 .

 22 Every point of the boundary of a convex subset of E belongs to one of its supporting hyperplanes.

  Delaunay triangulation Definition 19 (Delaunay cell). If R ⊂ S, and V (R) is a non empty Voronoi cell, then the Delaunay cell D(R) is cell(R). Definition 20 (Delaunay triangulation). The Delaunay triangulation D of S is the collection of Delaunay cells D(R), where R varies over subsets of S with V (R) non empty.

Figure 1 :

 1 Figure 1: Voronoi diagram and Delaunay triangulation of a data set S of size 10. We have C ∈ V S ({s 1 , s 2 }). So C is the center of an empty circumsphere of {s 1 , s 2 }. The point C 123 is the center of the circumsphere of the Delaunay triangle {s 1 , s 2 , s 3 }.

4 Figure 3 :

 43 Figure 3: Elimination condition of the principal axis tree.

3 Figure 4 :

 34 Figure 4: Optimization of the elimination condition for the quantization tree d 2 ≥ d 2 1 + d 2 2 + d 2 3 .

Figure 5 :

 5 Figure 5: Cell B is "hidden" from cell A.

Figure 7 :

 7 Figure7: If the query point q lies on the dark grey region H(x, s) ∩ V S ({s}) its nearest neighbor may be x.

Lemma 5 . 5 .

 55 With the same notations, one has U P I S (s) = p∈δ slabS(s) P I S (p). In other words, we have to check points of the boundary δ slab S (s) of slab S (s).

Figure 8 :

 8 Figure 8: Examples of friend Voronoi cells in a two-dimensional Voronoi diagram in the case of a bounded Voronoi cell (left) and in the unbounded case (right). In both case, the dark gray region is the considered Voronoi cell and the light gray regions are the friend Voronoi cells.

d = 2 d

 2 = 3 d = 4 d = 5 d = 6 d = 7 d = 8 Quantization tree 1.76s

Figure 10 :

 10 Figure 10: Execution time of 10 millions random queries on a data set of 10000 points, generated with a Gaussian pseudo random generator.

Figure 11 :

 11 Figure 11: Execution time of 10 millions random queries on a data set of 5000 points, generated with a uniform pseudo random generator.

  The centroidal projection associated C and X is the application Proj C,X : x →

	0	in the other case,	the centroids associated with X and C.
				N	G i 1 Ci (x).
				i=1
	Lemma 1.3 (Huyghens, variance decomposition). Let X be a E-valued L 2 random variable, N ∈ N *
	and C = (C i ) 1≤i≤N a Borel partition of E. Consider Proj C,X = projection. Then one has,	N i=1	G i 1 Ci the associated centroidal

  Definition 10 (supporting halfspace). Let C be a convex subset of E. A hyperplane H supports C if H ∩ C = ∅ and C is contained into one of the closed halfspaces defined by H.

	Lemma 2.1. Let C	E be a convex subset of E. If H is a supporting hyperplane of C, then every
	point of H ∩ C is a frontier point of C.
	Proof: Let H be a supporting hyperplane of C of equation φ

  2.75s 5.35s 8.93s 15.99s 28.06s 52.31s Principal axis tree 1.21s 1.86s 4.49s 10.87s 20.14s 41.56s 82.30s Kd-tree 1.88s 3.71s 8.54s 17.13s 31.06s 60.67s 118.93s Figure 9: Execution time of 10 millions random queries on a data set of 5000 points, generated with a Gaussian pseudo random generator.

		d = 2 d = 3 d = 4	d = 5	d = 6	d = 7	d = 8
	Quantization tree 2.59s 3.87s	6.46s 11.90s 27.54s 45.78s 84.63s
	Principal axis tree 1.33s 2.44s	4.94s 12.78s 41.02s 62.33s 119.88s
	Kd-tree	2.82s 5.20s 11.32s 24.20s 47.51s 87.61s 164.52s

The author is very grateful to Gilles Pagès (LPMA University Paris VI) for his helpful remarks and comments, and to Johan Mabille (Natixis) for his advices concerning the practical implementation.

V S (R) will represent the Voronoi cell of R in S. If C S is a Voronoi partition associated to S, and s ∈ S, slab S (s) will denote the Voronoi slab associated to S is the Voronoi partition C. Lemma 5.3. Let S ⊂ E be a nonempty finite set of sites in E and x ∈ E\S. Consider s ∈ S, the following assertions are equivalent:

Proof:

) is obvious owing to Proposition 5.1.

The intermediate value theorem shows that there exists λ * such that ∆ • f (λ * ) = 0 and thus f (λ * ) ∈ V S∪{x} ({x, s}).

The first modification made in the quantization tree algorithm is to assume that the points of the quantizer at each generation are points of the underlying codebook Γ. (In order to fulfill this requirement, we project an optimal quantizer onto the codebook.)

by their nearest neighbor projection on S. Consider q ∈ E. If q ∈ slab S (s) and t = Proj Γ (s) then {t, s} ∈ DG S∪{t} .

Proof: This is a straightforward consequence of the previous lemma.

Notation: Let S be a set of sites in E. For a point t in E, we denote P I S (t) = s ∈ S, {s, t} ∈ DG S∪{t} . The notation P I stands for "Pseudo-Insertion".

From an algorithmic viewpoint, the Delaunay graph of S being computed, P I S (t) stands for the sets of points in S, that are connected to t when updating the Delaunay graph to take account of this new point.

Implementing a procedure that computes P I S (t) is very similar to the insertion procedure of point t in T S . First friend node algorithm: This leads to a first method to compute a friend list:


Giovanni Funchal

Matthieu Moy

Modeling of Time in Discrete-Event Simulation of Systems-on-Chip

Today's consumer electronics industry uses modeling and simulation to cope with the complexity and time-to-market challenges of designing high-tech devices. In such context, Transaction-Level Modeling (TLM) is a widely spread modeling approach often used in conjunction with the IEEE standard SystemC discrete-event simulator.

In this paper, we present a novel approach to modeling time that distinguishes between instantaneous actions and tasks with a duration. We argue that this distinction should be natural to the user. In addition, we show that it gives us important insight and better comprehension of what actions can overlap in time. We are able to exploit this distinction to parallelize the simulation, achieving an important speedup and exposing subtle software bugs related to parallelism.

We propose a set of primitives and discuss the design decisions, expressiveness and semantics in depth. We present a research simulator called jTLM that implements all these ideas.

I. Introduction

Most of the functionality of modern high-tech consumer electronic devices is often grouped into a single integrated circuit, which is called a system-on-chip (SoC). The design of such systems is very complex and faces issues such as the time-to-market pressure. To cope with these constraints, there is a trend towards using models of the hardware in discrete-event simulation frameworks.

We call virtual prototype a specifically designed model of the hardware intended for development of software before the real, physical hardware is available. Virtual prototypes range from very high-level, application simulators such as that included in the iPhone SDK [1] to ones more targeted at low-level software development (drivers, etc.).

A. Transaction-level modeling

Transaction-level modeling [START_REF] Ghenassia | Transaction-Level Modeling with SystemC: TLM Concepts and Applications for Embedded Systems[END_REF] is a widely used technique for designing virtual prototypes intended for early software development. This approach tries to provide the "right" abstraction level in the sense of keeping just This paper has been partially supported by the French ANR project HELP (ANR-09-SEGI-006).

enough details so as to maintain the behavior perceived from a software programmer's point-of-view in what concerns the functionality.

Conceptually, a TLM model is a set of components, which represent hardware blocks (typically: CPUs, DMAs, memories, timers). The behavior is described in concurrent processes inside each component. Components communicate through interconnections, which represent memory-mapped buses. An interconnection transports transactions, which are abstractions of data exchanges. In the industry, most TLM models are written in SystemC [START_REF]IEEE 1666-2005 Standard SystemC Language Reference Manual[END_REF], [START_REF]TLM-2.0.1 User Manual, Open SystemC Initiative[END_REF], which is a C++ simulation library containing a cooperative, discrete-event scheduler.

B. Motivation and contributions

Previous works [START_REF] Helmstetter | Test coverage for loose timing annotations[END_REF], [START_REF] Cornet | A method for the efficient development of timed and untimed transaction-level models of systems-on-chip[END_REF] have identified some bad modeling practices in real-world TLM/SystemC models. Our experience suggests that some of these bad practices may be caused by confusion of modeling concepts inherent to the TLM approach with their implementation in the SystemC language.

In this context, we have developed a custom simulator in an attempt to measure the extent to which common modeling issues could be eliminated by providing new primitives that best suit the programmers' intents. We introduced the resulting framework (jTLM ) in [START_REF] Funchal | jTLM: an experimentation framework for the simulation of transaction-level models of systems-on-chip[END_REF]. That paper focused on the study of concurrency and confronted two modes of simulation scheduling: cooperative and preemptive.

In this paper, we now concentrate on the modeling of time. Typical TLM simulators such as SystemC do not allow durations to be expressed. In other words, actions do not "take" time, they happen instantaneously with respect to the simulated time. Although durations can be emulated by an instantaneous action followed (or preceded) by a pause, we argue that this is not natural to the user.

We introduce a novel way of modeling tasks with a duration and show that it gives us important insight and better comprehension of what actions can overlap. The notion of tasks with duration should be more intuitive to the user, since it is closer to the way the concrete system actually behaves.

This can also be exploited during simulation in various ways. For instance, it would be relatively easy to implement a trace back-end to rich visualization tools such as SimVision [START_REF] Simvision | [END_REF], allowing them to display non-instantaneous tasks as boxes on a time diagram.

In addition, having a clear definition of when tasks overlap allows us to better parallelize the simulation, achieving an important speedup and exposing more subtle software bugs related to parallelism that could not be detected previously.

C. Structure of the paper

To summarize, the main contribution of this paper is an innovative notion of simulation time, for which we provide both a clear semantics and an implementation. This improves on existing approaches by allowing the modeling of durations directly (i.e. not as a workaround).

The rest of the paper is structured as follows: Section II briefly presents the jTLM simulator. We introduce the main contribution, the notion of tasks with a duration, first informally in Section III, then more detailed in Section V. Section IV provides a discussion of the motivation and consequences of our choices. Then, we present experimental results in Section VII and related work in Section VIII.

II. Background: Overview of jTLM

In this section, we shall present jTLM briefly through a running example. Let us consider the virtual prototype of a system with two processors, a shared RAM memory and a DMA controller.

The DMA controller is a common hardware accelerator for transferring blocks of memory without overhead to the CPU. The software can program the DMA through its slave port, by writing to special addresses that are routed to registers in the DMA instead of the memory.

The DMA then performs the memory transfer through its master port.

A. Architecture and instantiation

Figure 2 shows the architecture of the example. Components represent the structure, and are connected by directed ports. The components are instantiated and connected before the simulation starts. Master ports initiate transactions; slave ports treat the requests. The interconnection routes transactions according to a map that tells which addresses belong to which slave ports.

B. Describing concurrency

Concurrency inside each component is described using "behaviors", which are sequential action streams in the form of a piece of code. jTLM provides two ways of simulating concurrency: cooperative and preemptive, detailed in [START_REF] Funchal | jTLM: an experimentation framework for the simulation of transaction-level models of systems-on-chip[END_REF]. We will give more details about the type of actions that can be performed in a behavior in the next sections.

The software in each of the processors in the example is straightforwardly described by wrapping its code inside a behavior. Hardware (in particular state machines), can also be described using a behavior. In our example, part of the DMA that is responsible for executing the transfers is described in the DMA f sm behavior (List. 1).

C. Simulated time/ordering of actions

jTLM has a notion of simulated time that represents or approximates the time by which actions happen on the concrete system. Simulated time is in principle completely disconnected from the wall-clock time, i.e. the time taken by the simulation when running on an ordinary computer. Figure 3 illustrates this fact by plotting both axes perpendicularly.

Although simulation time is usually slower than wallclock time, it can be faster depending on the precision of the model. In the field of wireless sensor networks, for instance, the interesting property is usually battery life, and therefore the simulators should be much faster than real-time. Time in jTLM is discrete, which means that the simulated time "hops" from one integer value to another. The user defines the granularity of the smallest "hop" which we call time unit. Since the time unit can be arbitrarily small, this does not limit the expressive power of jTLM.

Simulated time influences the order in which actions are observed in the simulator, i.e. an action that happens at simulated time 5 must happen before all actions that happen at greater simulated time. Actions at the same instant are not ordered by simulated time. We shall see more details in Section V.

D. Basic simulator primitives

jTLM is implemented on top of Java, which means that everything from the Java language is directly available to the user. In addition, jTLM introduces new actions to Java in the form of simulator primitives.

This section recalls the basic simulator primitives, originally introduced in [START_REF] Funchal | jTLM: an experimentation framework for the simulation of transaction-level models of systems-on-chip[END_REF] and very similar to those present in SystemC. The reader might want to skip to the Section III where we present the new primitives related to the modeling of tasks with a duration.

• awaitTime pauses the caller for the amount of simulated time specified as a parameter; and

• awaitEvent pauses until another behavior calls signalEvent on the same event. These primitives are directly inspired from SystemC's wait and notify. Our running example puts these primitives into practice in the DMA f sm behavior (List. 1); and in the register write implementation of the DMA slave port (List. 2).

The master.read and master.write commands in the DMA f sm behavior initiate transactions at the DMA's master port. On the other hand, the DMA slave port receives transactions and implements them. Our example has four registers: the "from" address, the "to" address, the transfer size and the START_TRANS register. The START_TRANS register is special in the sense that writing any value to it starts the DMA. This is carried out using an event.

By default, simulation time does not progress unless the behavior is waiting (inside awaitTime or awaitEvent). Everything else is instantaneous. At this point, the awaitTime(4) (line 7) is the closest way to model the fact that the reads and writes take some time (and this is what would have been done in SystemC).

III. New primitives

A. Tasks with known duration: consumesTime consumesTime is a new primitive that allows the modeling of sequences of actions that take time (a task ). In contrast with awaitTime(T) that creates an interval of time during which no action is executed, consumesTime(T) creates an interval during which a piece of code is executed.

Fig. 4 illustrates the semantics of consumesTime(T) by comparing three possible ways of expressing the fact that the computation of a function f() takes 30 units of simulation time (starting at time t = 10). The computation of f() itself is represented in bold dashed red.

Fig. 4a represents the state-of-the-art prior to this paper. To model a duration, we must use an instantaneous task followed by a delay (f(); awaitTime(30);). In simulation, the whole computation is executed at time 10, and the simulation time is increased afterwards.

In Figures 4b and4d, the task is modeled using consumesTime(30) { f(); }. In this case, each of the actions in f() will be executed somewhere in the interval [START_REF] Schumacher | parSC: Synchronous Parallel SystemC Simulation on Multi-Core Host Architectures[END_REF]40]. However the exact "path" of time inside the task is not precisely specified, only its duration. We represent this fact by a "curved" line in the figure, although in simulation the actual path will be a staircase function (because simulated time is discrete). The minimum width of the steps (granularity) is user-controlled by the time unit.

Two situations are worth mentioning. First, if the computation of f() is slow, simulation time might have to be blocked from advancing too much (at t = 40 in Figure 4b), in order to guarantee that every action of f() lies in the interval [START_REF] Schumacher | parSC: Synchronous Parallel SystemC Simulation on Multi-Core Host Architectures[END_REF]40].

If, for instance, the task consumesTime(30) {f();} was followed by an instantaneous task g() (Fig. 4c), then g() would be executed 30 units of simulated time after the consumesTime starts, no matter how long f() takes to execute in wall-clock time.

Second, while consumesTime(30) { f(); } allows f() to be executed at any time in the interval [START_REF] Schumacher | parSC: Synchronous Parallel SystemC Simulation on Multi-Core Host Architectures[END_REF]40], it cannot force the computation to last until the end of this interval. In other words, it is also possible for the computation of f() to finish early (at t = 30 in Figure 4d), in which case the task will suspend while the rest of the platform drives the advance of the simulated time, until it reaches 40.

Actually, the simulator has the freedom to execute the actions of the task at any time within the interval, including the particular case where the computation is executed completely at the start of the task. In other words, the execution of f(); awaitTime(30); is included in the possible executions of consumesTime(30) {f();}. By using consumesTime instead of awaitTime, we can relax the semantics of the model where it would otherwise be over-specified.

We can apply this primitive to the DMA example presented in the previous section. In contrast with the List. 1 where we had to model the time taken by the reads and writes using instantaneous actions followed by a wait, we can now express the fact that the DMA transfer takes an amount of time proportional to the length of the transfer by using a dedicated primitive. One of the immediate advantages of this is the readability of the code, as shown in List. 3. The rest of the advantages are discussed in Section IV.

B. Loose timings

Loose timings [START_REF] Helmstetter | Test coverage for loose timing annotations[END_REF] consists in the addition of nondeterminism to the duration of tasks or awaitTime statements. This increases the faithfulness of the model when we have only partial or imprecise knowledge of the timings. List. 4 shows an application of loose timings in conjunction with consumesTime to provide a ±10% tolerance. In the simulator, this is implemented using randomization (possibly non-uniformly distributed).

C. Tasks with unknown duration: consumesUnknownTime consumesTime allowed modeling tasks of known duration. There are times however when we cannot quantify in advance the duration of an action. This can be either because we are very early in the design-flow and such timing information is not available, or because the duration depends on external actions.

For instance, in a high-level model of a complex processor with pipelines, caches, etc., it is very hard to quantify precisely the time taken by each iteration of the low-level software loop in List. 5. This loop implements

} }

Listing 6: Software loop in jTLM with consumesUn-knownTime part of the lock acquire of a mutex. The test_and_set instruction generates a transaction that atomically reads the address x and replaces it with value 1. If the value read was 0, then the lock was free and we have just acquired it, otherwise the lock was busy and we must try again.

It is a common practice to put a cpu_relax() instruction in the loop, in the case the lock is found busy. This is important for performance reasons: depending on the processor architecture, the implementation of this instruction might cause the processor to enter a lowconsumption mode or reduce its bus bandwidth usage.

The interesting property in the example is that the loop will not exit until it reads x = 0. However, the write to x is an unpredictable external action by a different behavior. We cannot know in advance when and if the variable x will be written. This means that there is no a-priori bound to the time between when the polling loop starts, and when the new value of x allows it to complete.

Of course, the intuition says that any update to x should become visible (i.e. the behavior reading x is given an opportunity to be executed) after a reasonable amount of time. This is indeed the case most of the time in real life, but modeling this requires some sort of fairness.

Having these considerations in mind, we have devised a primitive that would fulfill our intents precisely in the above cases. We call it consumesUnknownTime. Informally, this primitive guarantees that the time taken by a task is at least enough to satisfy any conditions needed to complete the task, but not necessarily bounded. If the conditions are never met, the task never finishes. In other words, while consumesTime fixes the duration of the task and lets the actions be scheduled at any time in the specified interval, consumesUnknownTime executes the actions contained in the task until they complete regardless of the time they take to execute.

In the preemptive mode of jTLM [START_REF] Funchal | jTLM: an experimentation framework for the simulation of transaction-level models of systems-on-chip[END_REF], the software loop of the previous example can be implemented in a straightforward way with the help of this new primitive (List. 6). The only thing that requires special attention is cpu_relax, which is redirected to Java Threads' library yield() method.

Some may argue that polling loops are to be avoided when possible, but real-life hardware systems sometimes have no alternative, and the model of the system must then use the same mechanism (e.g. if the actual software uses polling and is embedded directly in the TLM platform). The existence of consumesUnknownTime in jTLM is not meant to promote polling as a synchronization mechanism, but to allow modeling it efficiently when needed.

Consider another example: consumesUnknownTime {f();} g();

During the execution of the consumesUnknownTime statement, simulated time is not constrained, i.e. other behaviors are allowed to continue their execution, advance time, etc. However, g() will always be executed after f() completes, both in terms of simulation and of wallclock time. The current simulation time when g() starts its execution is defined as the simulation time when consumesUnknownTime {f();} completed.

As a closing remark to this section, notice that we do not provide an implementation of consumesUnknownTime in the cooperative mode of jTLM. Section V-D is dedicated to this issue.

IV. Discussion

One of the fundamental ideas of jTLM is to expose the parallelism of the simulated platform to the user, as opposed to hiding it as it is done in traditional cooperative discrete-event simulators. The system being simulated is parallel (since hardware is intrinsically parallel), and it is desirable to exploit this parallelism in the simulation. A detailed discussion on the topic can be found in [START_REF] Funchal | jTLM: an experimentation framework for the simulation of transaction-level models of systems-on-chip[END_REF].

This section discusses the improvements achieved by the new primitives proposed in this paper to two domains: exploiting the physical parallelism of the host machine for better parallelization, and finding more bugs early in the systems-on-chip design-flow.

A. Better simulator parallelization

Nowadays machines are usually multi-core, and simulating a concurrent system in a purely sequential way is clearly suboptimal. Still, the current industry standard in the domain, SystemC, cannot easily be made parallel since the standard requires "co-routine" semantics (i.e. cooperative multitasking thus no preemption). Previous attempts to parallelize the execution in a semanticspreserving way [START_REF] Bouzouzou | Accélération des simulations de systèmes-sur-puce au niveau transactionnel[END_REF] showed the difficulty of the problem. In practice, this would also be inefficient since most TLM programs have lots of shared variables that would need to be protected (the RAMs in particular).

Other approaches, like that in [START_REF] Schumacher | parSC: Synchronous Parallel SystemC Simulation on Multi-Core Host Architectures[END_REF] and in our previous work with jTLM [START_REF] Funchal | jTLM: an experimentation framework for the simulation of transaction-level models of systems-on-chip[END_REF], chose to break with co-routine semantics. These approaches are able to parallelize the simulation of unrelated actions occurring at the same simulation instant. Still, none of them is able to run two unrelated actions that happen at different simulation times in parallel.

Unfortunately, the situation where only one process runs at a given simulation time is quite common in TLM models because they do not use clocks. The presence of loose timings makes it even more common.

The new notion of tasks with a duration in jTLM allows to overcome this limitation. Back to our running example of Section II, suppose that CPU 0 executes a task A during interval [START_REF] Schumacher | parSC: Synchronous Parallel SystemC Simulation on Multi-Core Host Architectures[END_REF]20] and CPU 1 executes a task B during interval [15, 30]. There is an overlap between times 15 and 20. This means that if the code of task A is not finished when the simulation time reaches 15, then the code of task B can start in parallel with the code of A.

While the chance of having simultaneous, instantaneous actions is small, tasks with a duration have much greater chances to overlap. Hence, the parallelism is exploited far better, with a little help from the user.

B. Finding more bugs

In the context of virtual prototyping, it is important to distinguish bugs of the prototype, and bugs that not only appear on the prototype but also correspond to bugs in the real system. In the latter case, it is desirable to exhibit bugs as early as possible, avoiding surprises later in the design-flow.

Cooperative multitasking completely eliminates simulation parallelism. This avoids bugs of the prototype due to parallelism, simplifying the development, but may also hide bugs of the real system that could otherwise have been found. In particular, if the simulator does not expose the parallelism of the hardware to the software, then mis-synchronized software will run correctly on the simulator, but fail to run on the final chip. We have already showed in [START_REF] Funchal | jTLM: an experimentation framework for the simulation of transaction-level models of systems-on-chip[END_REF] how physical parallelism within a simulation instant could exhibit software bugs. However, previous works did not solve the problem of finding bugs that span across instants.

Returning to our example, suppose that one of the processors writes to a shared variable at time t = 10, and that another reads the same variable at t = 11 without any synchronization. This probably represents a bug in the real system: i.e. since there is no synchronization, nothing guarantees that the read will see the effect of the write. In fact, the read may even return an unpredictable value ("out-of-thin-air") because the write may have been buffered or still be in progress. This condition is known as a "data-race" [START_REF] Adve | Shared memory consistency models: A tutorial[END_REF] and is very difficult to debug.

Traditional simulators, even with parallelism within instants, would have to insert artificial synchronization at the end of instant t = 10 in order to advance simulation time and unblock any behaviors executing at t = 11. The synchronization implies that, in the simulator, the write at t = 10 is globally performed before t = 11. In other words, the bug described above is not visible on such simulators.

With the introduction of tasks with a duration, the above example can be modeled by a read and a write within overlapping tasks. Then there will be no synchronization between them, and the data-race is made visible in two ways: (1) the jTLM scheduler can choose the order of the read and write operations, possibly leading to different interleavings; (2) the semantics of the read and write are the ones of concurrent read/write in the underlying Java memory model [START_REF] Microsystems | Java Memory Model and Thread Specification[END_REF]. This means that the operations are not always atomic, and allows finding bugs that are not exposed by simple interleaving semantics, such as data-races.

While the (1) above could be implemented with some effort in SystemC, the (2) is really not straightforward, because of the complexity of modern memory models. A detailed discussion of this issue is clearly out of the scope of the present paper.

C. Non-reproducibility

Compared to cooperative simulators, the preemptive mode of jTLM has a drawback in terms of reproducibility due to the way it exploits concurrency to achieve parallelism. This is not a consequence of the notion of noninstantaneous tasks presented in this paper, but rather inherent to jTLM [START_REF] Funchal | jTLM: an experimentation framework for the simulation of transaction-level models of systems-on-chip[END_REF]. We still provide a cooperative mode for jTLM if the user needs reproducible executions.

V. Details of the semantics

By default, actions in jTLM are instantaneous. Each of the primitives we have presented so far controls the advance of simulated time from the point of view of the caller. The previous sections gave the intuition of the semantics of the primitives in the general case. We now discuss the details of the semantics in special circumstances.

A. awaitTime and instantaneous actions

If a behavior calls awaitTime [START_REF] Schumacher | parSC: Synchronous Parallel SystemC Simulation on Multi-Core Host Architectures[END_REF] at time 0, then the behavior cannot do any action for 10 units of time, and its next action (f() in Fig. 5) will only begin at time 10. If any other behavior has an action at t < 10, then this action (g()) will be simulated before f(). If there are several instantaneous actions at the same time (f() and h()), the simulated time cannot advance until they all have completed.

t = 0 Wall-clock time Simulation time await (10) f() await(5) g() await(5) h()
t = 5 t = 10

B. consumesTime and tasks with duration

When a behavior starts at time 0 a task with the duration of 20 time units (consumesTime(20)), all actions that are part of that task will be simulated somewhere between times 0 and 20, inclusive. This means that a task with duration 0, such as consumesTime(0){x()}, is equivalent to an instantaneous task x(); and that an empty task, such as consumesTime(10){}, is equivalent to a call to awaitTime [START_REF] Schumacher | parSC: Synchronous Parallel SystemC Simulation on Multi-Core Host Architectures[END_REF].

Figure 6 represents two tasks that overlap over the instant t = 20. Actions in the overlapping portion are not ordered by simulated time, but may be ordered by dependency (see next section). In the Fig. 6, f() and g() are overlapping and not ordered by dependency, and thus can be run in parallel during the simulation. For the time being, nesting calls to consumesTime, awaitTime or awaitEvent inside a consumesTime is forbidden. We provide some ideas on how to deal with this in the conclusion.

C. awaitEvent and signalEvent

When a behavior starts an awaitEvent, its next action will only begin after another behavior calls signalEvent on the same event. Simulated time can pass while the behavior is inside awaitEvent. Once the event is signaled, the behavior wakes up immediately and executes its next action at the same simulated time as the notification. The signalEvent/awaitEvent pair implies dependency, which means that in Figure 7, even though f() and g() occur at the same simulated time, there is no overlap.

D. consumesUnknownTime

The semantics of consumesUnknownTime is the hardest to understand because we cannot quantify its duration in advance. While instantaneous tasks completely prevent time from advancing, and fixed duration tasks (consumesTime) block time from advancing too much, consumesUnknownTime imposes no upper bound on the amount of time that can pass. However, the lower bound is such that any conditions needed to complete the code inside the task must be satisfied before it can finish.

A trivial implementation that is always valid consist in letting an infinite amount of time pass (equivalent to awaitTime(∞)) without executing the body of the task. Of course, we wish to avoid this kind of unrealistic executions. In the preemptive mode of jTLM, our implementation achieves this by relying on the fairness of the underlying operating system. However, there is no straightforward way to do the same thing in the cooperative mode. We therefore do not provide an implementation of cooperative consumesUnknownTime for the time being, but this will be the subject of future work.

VI. jTLM implementation

The algorithm behind jTLM involves relatively tricky synchronization at some points because of concurrency.

For brevity, we focus here only on the principle of the algorithm, omitting low-level details.

Let B be the set of all behaviors, E be the set of all events, and T be the time unit. The jTLM scheduler has the following data structures: AG : T → 2 B is a priority queue implementing an agenda of behaviors waiting or consuming time; and P B : E → 2 B maps a set of pending behaviors for each event. Additionally, the scheduler must keep track of the current simulated time (ct ∈ T) and integer counters for the number of behaviors running an instantaneous section of code (bi), awaiting an event (be) and consuming unknown time (bu). Each behavior is associated with a semaphore used to block and unblock it. We must also identify the calling behavior (cb ∈ B) which invokes the primitives.

The heart of the scheduler is in the timeElapse() method, which is called at the beginning of the simulation and every time bi reaches zero. This method sets the current time to the earliest deadline in AG and then either: (a) in the preemptive mode, wakes up all behaviors with deadline ct; or (b) in the cooperative mode, picks a random behavior among those with deadline ct and wakes it up. Here, "waking up" consists of incrementing bi and releasing a semaphore to unblock the behavior. The simulation ends if timeElapse() is called but AG = ∅ ∧ bu = 0. awaitTime(t) inserts a pair (ct+t, cb) into AG, and then decrements bi, checks if bi = 0 (in which case timeElapse() must be called), and blocks the behavior by acquiring a semaphore. awaitEvent(e) is similar, but increments be instead of adding an element to AG.

In the preemptive mode, signalEvent(e) simply wakes up all b ∈ P B(e), removing them from the set, and decrementing be by |P B(e)|. In the cooperative mode, we must only let one behavior run among the possible candidates in P B(e), cb and AG(ct). To achieve this, our (unoptimized) algorithm follows a similar procedure but inserts {ct} × P B(e) ∪ {cb} into AG and then calls timeElapse().

Although the semantics of consumesTime and awaitTime are very different, their implementation in preemptive mode are very similar. The only (subtle) difference is: before acquiring the semaphore, consumesTime first passes control to the user code. While it is running, the pair (ct+ t, cb) previously inserted into AG acts as a "reservation" ensuring that the simulated time does not advance more than it should. In the cooperative mode, consumesTime(t) randomly splits the time interval [0, t] into two calls to awaitTime before and after executing the user code. We use a technique based on Java's inline unnamed functors to implement this.

As discussed in Sec. V, we only implement consumesUnknownTime in the preemptive mode. This implementation is similar to that of consumesTime, but increments bu instead of adding an element to AG and decrements bu instead of acquiring the semaphore after running the user code.

VII. Benchmarks

This paper focuses mainly on modeling issues. Nevertheless, we provide some results obtained from experiments in this section. Our intention is to show to which extent the concurrency described with the help of the new primitives can be exploited in a parallel simulation.

For this, we need an "embarrassingly" parallel example, which does not need to be realistic but allows us to measure the maximum expected improvement. Our model is organized as one main processor and 100 hardware accelerators, each of which calculates 200 hexadecimal digits of π using the BBP formula [START_REF] Bailey | On the rapid computation of various polylogarithmic constants[END_REF]. The main processor assigns tasks to each accelerator in a loop, with a small delay between each assignment, then waits for them all to complete. We conducted our benchmarks on a server with 32 Intel Xeon processors.

In the first version, each accelerator is implemented with an instantaneous task followed by an awaitTime. Because of the small delay between each assignment, the instantaneous tasks happen at different simulation times, so there is no overlap. One task cannot start before the previous has completed and let the time elapse. The computation takes 220s.

In the second version, the accelerators are implemented using consumesTime to model the fact that each accelerator task has a duration. The simulator knows that there is an overlap between tasks and is able to exploit this by running the tasks in parallel. Now, the computation takes only 20s.

The results confirm our intuitive expectations and show the practical benefits of consumesTime for parallelization. In addition, we are convinced from our tests and our experience that these primitives are natural for the user and using them requires little effort.

VIII. Related work

A. SystemC

The SystemC [START_REF]IEEE 1666-2005 Standard SystemC Language Reference Manual[END_REF] modeling language is the current industry standard for developing transaction-level models. Strictly speaking, SystemC is a C++ library that includes a simulation scheduler and data-types specially designed for describing hardware structures such as wires and registers.

For describing concurrency, SystemC includes a notion of process. The SystemC standard requires processes to have co-routine semantics: they run in isolation and in a cooperative manner. Each process either runs to completion or calls at some point a primitive that yields control back to the scheduler. In other words, the scheduler is not able to preempt the running process. In addition, actions of a process are always instantaneous w.r.t. the simulated time.

In contrast, jTLM includes both a preemptive mode and a cooperative mode [START_REF] Funchal | jTLM: an experimentation framework for the simulation of transaction-level models of systems-on-chip[END_REF]. In the preemptive mode, jTLM can naturally run processes in parallel both inside of an instant and by exploiting the notion of overlap of tasks, as seen in Section IV.

B. SpecC

SpecC [14] is another language for system-level simulation of hardware systems. Unlike SystemC and jTLM, SpecC is implemented as a new language, not as a library. Its underlying concepts are the usual ones in discreteevent simulation. The semantics of SpecC do not impose a preemptive or cooperative simulation, and the user code is expected to work in either cases. However, the reference implementation is cooperative.

Simulated time in SpecC is managed in a similar way to that of SystemC, with instantaneous actions and waitfor statements.

An interesting feature of SpecC is the notion of timing constraints ( §2.4.9 of the LRM [14]). This feature allows the user to label arbitrary points in the execution, and specify constraints to the difference of simulation time between the execution of pairs of labels. This is a very general specification mechanism, which can be used to express constraints on independent tasks, but also on overlapping or nested time intervals.

However, unlike the consumesTime primitive of jTLM, this feature is only a specification mechanism: it does not influence the behavior of the simulation. Instead, it only allows one to check that the execution complies with the specification. The simulation itself is still derived from instantaneous actions spread in time with waitfor statements.

IX. Conclusion

This paper presented a novel approach to modeling time in discrete-event simulators of transaction-level models of systems-on-chip. Previous works such as [START_REF]IEEE 1666-2005 Standard SystemC Language Reference Manual[END_REF], [14] did not distinguish between instantaneous actions and tasks that take time.

We have proposed a notion of tasks with a duration that we think is both intuitive and natural to the user. In addition, we showed that this notion could be exploited:

• to enrich trace visualization tools, allowing them to display tasks as boxes on a time diagram;

• to derive a clear definition of overlapping tasks;

• to effortlessly achieve an important simulation speedup by enabling parallel execution of actions occurring at different simulation times;

• to expose bugs in simulation that may correspond to bugs in the real system, by removing the constraint that actions at different simulation times are necessarily synchronized.

Furthermore, we have described our implementation of a custom research simulator called jTLM [START_REF] Funchal | jTLM: an experimentation framework for the simulation of transaction-level models of systems-on-chip[END_REF]. This allowed us to confirm the benefits of our approach both in terms of ease of use and parallelization.

In our future work, we intend to study how these features could be incorporated in a more complex simulator such as SystemC. Another direction of future work concerns consumesUnknownTime. It is particularly difficult to avoid implementations of this primitive that only have trivial executions, i.e. where some behaviors never get a chance to run.

Finally, the last direction of work is related to the nesting of calls to consumesTime and awaitTime. In this sense, there are many pitfalls that should be avoided. For instance, the naive approach could lead to the situation in Figure 8. This figure shows a task with duration 50 and, after 40 units of time have passed, a nested task with duration 20 starts. There is not enough time left in the outer task to accommodate for the inner task. In this case, we would have to either abort or violate the duration of one of the tasks, neither of which is satisfactory. One way to get around this issue is to make inner calls suspend the outer calls for their duration as shown in Fig. 9. To implement this, one would need to keep track of the nesting level, in order to be able to decide when exiting a task if there is an outer task that needs to be resumed.

Figure 1 :

 1 Figure 1: Conceptual view of a typical TLM model

Figure 2 :

 2 Figure 2: Our running example, a virtual platform in jTLM

Figure 3 :

 3 Figure 3: Simulation time vs. wall-clock time

Figure 4 :

 4 Figure 4: Semantics of the new jTLM primitives

Figure 5 :

 5 Figure 5: Semantics of awaitTime

Figure 6 :

 6 Figure 6: Semantics of consumesTime

Figure 7 :

 7 Figure 7: Semantics of awaitEvent/signalEvent

Figure 8 :

 8 Figure 8: The naive semantics for nested calls to consumesTime

Figure 9 :

 9 Figure 9: Suspend semantics for nested calls to consumesTime

Acknowledgments

Many thanks to Nabila Abdessaied, Mohamed El Aissaoui and Rafael Velasquez, who contributed to the initial version of the jTLM implementation and Jérôme Cornet and Laurie Lugrin for the remarks.