Borzoo Bonakdarpour
email: borzoo@cs.uwaterloo.ca

Marius Bozga
email: marius.bozga@imag.fr

Jean Quilbeuf
email: jean.quilbeuf@imag.fr

Automated Distributed Implementation of Component-based Models with Priorities *

Keywords: C.2.4 [Computer-Communication Networks]: Distributed Systems[Distributed applications], D.1.3 [Programming Techniques]: Concurrent Programming-Distributed programming, Parallel programming, D.2.13 [Software Engineering]: Reusable Software-Reuse models, D.4.7 [Operating Systems]: Organization and Design-Real-time and embedded systems, F.3.1 [Logics Specifying and Verifying and Reasoning about Programs-Logic of programs, I.2.2 [Artificial Intelligence ]: Automatic Programming-Program transformation Theory, Design, Languages, Reliability, Performance Component-based modeling, Automated transformation, Distributed systems, BIP, Correctness-by-construction, Committee coordination, Conflict resolution

In this paper, we introduce a novel model-based approach for constructing correct distributed implementation of component-based models constrained by priorities. We argue that model-based methods are especially of interest in the context of distributed embedded system due to their inherent complexity. Our three-phase method's input is a model specified in terms of a set of behavioural components that interact through a set of high-level synchronization primitives (e.g., rendezvous and broadcasts) and priority rules for scheduling purposes. Our technique, first, transforms the input model into a model that has no priorities. Then, it transforms the deprioritized model into another model that resolves distributed conflicts by incorporating a solution to the committee coordination problem. Finally, it generates distributed code using asynchronous point-topoint send/receive primitives. All transformations preserve the properties of their input model by ensuring observational equivalence. The transformations are implemented and our experiments validate their effectiveness.

INTRODUCTION

Correct design and implementation of computing systems has been an ongoing research topic in the past three decades. This problem is significantly more challenging in the context of distributed systems due to a number of factors such as non-determinism, non-atomic execution of processes, race conditions, and occurrence of faults. Correctness of distributed implementations is of significant importance in the context of embedded applications, as such applications are often employed in safety-critical systems. Model-based development of embedded distributed applications aims at increasing their integrity by using explicit models employed in clearly defined transformation steps leading to correct-byconstruction artifacts. This approach is beneficial, as one can ensure functional correctness of the system by dealing with a high-level formally specified model that abstracts implementation details and then derives a correct implementation through a series of transformations that terminates when an actual executable code is obtained.

In this paper, we focus on the BIP framework [START_REF] Basu | Modeling heterogeneous real-time components in BIP[END_REF] as our formal modelling language. BIP (Behaviour, Interaction, Priority) is based on a semantic model encompassing composition of heterogeneous components. The behaviour of components is described as an automaton or Petri net extended by data and functions given in C++. BIP uses a diverse set of composition operators for obtaining composite components from a set of components. The operators are parametrized by a set of interactions between the composed components. Finally, priorities are used to specify different scheduling mechanisms 1 . Transforming a BIP model into a distributed implementation involves addressing three fundamental issues:

1. Concurrency: Components and interactions should be able to run concurrently while respecting the sequential semantics of the high-level model.

Conflict resolution:

Interactions that share a common component can potentially conflict with each other.

Enforcing priorities:

When two interactions can execute simultaneously, the one with higher priority must be executed.

These issues introduce challenging problems in a distributed setting. The conflict resolution issue can be addressed by incorporating solutions to the committee coordination problem [START_REF] Chandy | Parallel program design: a foundation[END_REF] for implementing multiparty interactions. For example, Bagrodia [START_REF] Bagrodia | A distributed algorithm to implement n-party rendevouz[END_REF] proposes different solutions with different degrees of parallelism. The most distributed solution is based on the drinking philosophers problem [START_REF] Chandy | The drinking philosophers problem[END_REF], and has inspired the approaches by Pérez et al. [START_REF] Pérez | An order-based algorithm for implementing multiparty synchronization[END_REF] and Parrow et al. [START_REF] Parrow | Multiway synchronizaton verified with coupled simulation[END_REF]. In the context of BIP, a transformation addressing all the three challenges through employing centralized scheduler is proposed in [START_REF] Basu | Distributed semantics and implementation for systems with interaction and priority[END_REF]. Moreover, in [START_REF] Bonakdarpour | Automated conflict-free distributed implementation of component-based models[END_REF][START_REF] Bonakdarpour | From high-level component-based models to distributed implementations[END_REF], we propose transformations that address the concurrency issue by breaking the atomicity of interactions and conflict resolution by embedding a solution to the committee coordination problem in a distributed fashion. On the contrary, designing transformations that enforce priorities between interactions in a distributed setting remains unaddressed in spite of the vital role specifying priorities plays in designing systems.

Motivation

Priorities are widely used in system design, as a way of scheduling events. Below, we present examples of how applying priorities can guide a system to satisfy certain properties:

• Ensuring safety. Safety properties are normally of the form "nothing bad happens during the system execution". In the context of concurrent and distributed computing, such bad things are often due to existence of a set of processes competing over a resource. Priorities can be used to resolve such race conditions. For instance, one way to prevent two processes to enter a critical section simultaneously is to give explicit priority to one process. Dynamic priorities can then be used to ensure non-starvation.

• Improving performance.

In distributed systems, it is often the case that certain resources have higher demands. For example, in group mutual exclusion [START_REF] Jurdzinski | Small progress measures for solving parity games[END_REF], as Mittal and Mohan argue [START_REF] Mittal | A priority-based distributed group mutual exclusion algorithm when group access is non-uniform[END_REF], in many commonly considered systems, group access requests are non-uniform. Hence, in order to improve the performance, it is reasonable to devise algorithms that give priority to groups that require resources with higher demand. A concrete example of group mutual exclusion is the well-known readers/writers problem. In most cases, we give priority to readers to improve the performance.

• Reducing non-determinism. Non-determinism in distributed and concurrent computing is one of the sources of obtaining a diverse set of behaviours. In many scenarios and in particular, in embedded applications, it is desirable to guide the system to behave in a certain predictable fashion.

The main challenge in ensuring priorities in a distributed setting is their correct implementation. This is due to the fact that components need to obtain a reliable knowledge about enabledness of interactions, so that only the interaction with highest priority is executed. In [START_REF] Basu | Priority scheduling of distributed systems based on model checking[END_REF], the authors propose a model checking approach that determines whether actions of a given Petri net can be executed without violating priority rules. However, the downside of this approach is (1) it has scaling issues, as it uses model checking, and (2) in most cases the local knowledge of processes is shown to be insufficient to decide whether or not an action can be executed. Other approaches include applying customized algorithms to implement priority rules for specific problems in distributed computing (e.g., [START_REF] Mittal | A priority-based distributed group mutual exclusion algorithm when group access is non-uniform[END_REF]).

These examples demonstrate the demand for developing methods that automatically construct a correct distributed implementation by starting from a high-level model along with a set of priority rules. This way, all implementation issues are dealt with by transformation algorithms and designers only need to make minimal effort to develop models.

Contributions

Our contributions in this paper are as follows:

• We propose a transformation that, given a high-level BIP model with priorities, generates a BIP model without priorities, that behaves equivalently. This corresponds to the first step in Figure 1.

• We show the correctness of this transformation by proving that the initial and transformed models are observationally equivalent.

• We apply the transformation introduced in [START_REF] Bonakdarpour | From high-level component-based models to distributed implementations[END_REF] to derive a distributed model, where multiparty interactions are implemented in terms of asynchronous point-to-point send/receive primitives. This corresponds to the second step in Figure 1. From this distributed model, we generate distributed code, as explained in [START_REF] Bonakdarpour | Automated conflict-free distributed implementation of component-based models[END_REF][START_REF] Bonakdarpour | From high-level component-based models to distributed implementations[END_REF], which completes the design flow from the initial BIP model with priorities to a correct distributed implementation.

• Finally, we validate the effectiveness of our approach by modelling a jukebox application in BIP and conducting experiments on the generated distributed code.

The jukebox application incorporates priorities to manage demands on reading discs and our experiments show that the overhead of our transformations has minimal effect on the benefit of using priorities.

Organization.

The rest of the paper is organized as follows. In Section 2, we present the basic semantics model of BIP. Then, in Section 3, we describe our transformation for deriving a model that has no priorities. Our approach for deriving a distributed model and code is presented in Section 4. We discuss our case study and experimental results in Section 5. Finally, we conclude in Section 6.

BASIC SEMANTIC MODELS OF BIP

In this section, we present operational global state semantics of BIP [START_REF] Basu | Modeling heterogeneous real-time components in BIP[END_REF]. BIP is a component framework for constructing systems by superposing three layers of modelling: Behaviour, Interaction, and Priority. Atomic Components. We define atomic components as transition systems extended with a set of ports and a set of variables. Each transition is guarded by a predicate on the variables, triggers an update function, and is labelled by a port. The ports are used for communication among different components and each port is associated with a subset of variables of the component.

Definition 1 (Atomic Component

). An atomic component B is a labelled transition system represented by a tuple (Q, X, P, T ) where:

• Q is a set of control states.
• X is a set of variables.

• P is a set of communication ports. Each port is a pair (p, Xp) where p is a label and Xp ⊆ X is the set of variables bound to p. By abuse of notation, we denote a port (p, Xp) by p.

• T is a set of transitions of the form τ = (q, p, g, f, q ′ ) where q, q ′ ∈ Q are control states, p ∈ P is a port, g is the guard of τ and f is the update function of τ . g is a predicate defined over the variables in X and f is a function that computes new values for X according to the previous ones.

We denote X the set of valuations on X, and Q × X the set of local states. Let (q, v) and (q ′ , v ′ ) be two states in Q × X, p be a port in P , and v ′′ p be a valuation in Xp of Xp. We write (q, v)

p(v ′′ p ) -→ (q ′ , v ′ ), iff τ = (q, p, g, f, q ′ ) ∈ T , g(v) is true, and v ′ = f (v[Xp ← v ′′ p ]
), (i.e., v ′ is obtained by applying f after updating variables Xp associated to p by the values v ′′ p ). When the communication port is irrelevant, we simply write (q, v) → (q ′ , v ′ ). Similarly, (q, v) p → means that there exists a transition τ = (q, p, g, f, q ′ ) such that g(v) is true; i.e., p is enabled in state (q, v).

Figure 2(a) shows an atomic component B, where Q = {s}, X = {n}, P = {(p, {n})}, and T = {(s, p, g, f, s)}.

Here g is always true and f is the identity function.

Interactions. For a model built from a set of n atomic components {Bi = (Qi, Xi, Pi, Ti)} n i=1 , we assume that their respective sets of ports and variables are pairwise disjoint; i.e., for any two i = j in {1..n}, we require that Pi ∩ Pj = ∅ and Xi ∩ Xj = ∅. Thus, we define the set P = n i=1 Pi of all ports in the model as well as the set X = n i=1 Xi of all variables. An interaction a is a triple (Pa, Ga, Fa), where Pa ⊆ P is a set of ports, Ga is a guard, and Fa is an update function, both defined on the variables associated by the ports in Pa (i.e., p∈Pa Xp). By Pa = {pi}i∈I , we mean that for all i ∈ I, pi ∈ Pi, where I ⊆ {1..n}. We denote by F i a the projection of Fa on Xp i .

Priorities. Given a set γ of interactions, a priority between two interactions specifies which one is preferred over the other. We define such priorities through a partial order π ⊆ γ × γ. We write aπb if (a, b) ∈ π, which means that a has less priority than b.

Definition 2 (Composite Component). A composite component (or simply component) is defined by a set of components, composed by a set of interactions γ and a priority partial order π ⊆ γ ×γ. We denote B def = πγ(B1, . . . , Bn) the component obtained by composing components B1, • • • , Bn using the interactions γ and the priorities π.

Note that if the system does not contain any priority, we may omit π.

Definition 3 (Composite Component Semantics).

The behaviour of a composite component without priority γ(B1, • • • , Bn), where Bi = (Qi, Xi, Pi, Ti) and →i is the transition relation between states of Bi, is a transition system (Q, γ, X, →γ), where Q = Π n i=1 Qi, X = n i=1 Xi and →γ is the least set of transitions satisfying the rule:

a = ({pi}i∈I , Ga, Fa) ∈ γ Ga(v1, . . . , vn) ∀i ∈ I. (qi, vi) = (q ′ i , v ′ i ) ∀i ∈ I. (qi, vi) p i (v ′′ p i ) -→ i (q ′ i , v ′ i ), v ′′ p i = F i a (v1, . . . , vn) ((q1, v1), . . . , (qn, vn)) a →γ ((q ′ 1 , v ′ 1 ), . . . , (q ′ n , v ′ n ))
We denote (q, v) the state of γ(B1,

• • • , Bn) that correspond to the states (q1, v1), • • • , (qn, vn) of the components B1, • • • , Bn.
We define the behaviour of the composite component B = πγ(B1, . . . , Bn) as the transition system (Q, γ, X, →π) where →π is the least set of transitions satisfying the rule: Intuitively, the first inference rule specifies that a composite component B = γ(B1, . . . , Bn) can execute an interaction a ∈ γ, iff (1) for each port pi ∈ Pa, the corresponding atomic component Bi can execute a transition labelled by pi, and (2) the guard Ga of the interaction evaluates to true in the current state. Execution of the interaction modifies components' variables by first applying update function Fa to associated variables and then function fi inside each component. The states of components that do not participate in the interaction stay unchanged. The second inference rule simply filters out transitions which are not maximal with respect to priorities. A transition is executed only if no other one with higher priority is enabled. It may be desirable to always execute interaction a when possible. This can be done by adding the two priority rules bπa and cπa. We denote the obtained component by πγ(B1, . . . , B4). We will use this example to illustrate the transformations presented in this paper.

(q, v) a →γ (q ′ , v ′ ) ∀a ′ ∈ γ. aπa ′ =⇒ (q, v) a ′ →γ (q, v) a →π (q ′ , v ′ ) s p p n (a) An atomic component s 1 p 1 p 1 n 1 s 2 p 2 p 2 n 2 s 3 p 3 p 3 n 3 s 4 p 4 p 4 n 4 a = p 1 p 2 ⌈n 1 > n 2 ⌋ sw(n 1 , n 2 ) b = p 2 p 3 ⌈n 2 > n 3 ⌋ sw(n 2 , n 3 ) c = p 3 p 4 ⌈n 3 > n 4 ⌋ sw(n 3 , n 4 ) (b) A BIP composite component that
We now introduce the notion of conflicting interactions. Intuitively, two interactions a1 and a2 are weakly conflicting iff they share a common component.

Definition 4 (Weak Conflict). Two interactions a1 and a2 are weakly conflicting (denoted a1 ⊕a2) iff there exist two ports p and q in some component B such that p ∈ Pa 1 and q ∈ Pa 2 .

This kind of conflict is called weak because it is weaker than the definition of conflict in [START_REF] Bonakdarpour | From high-level component-based models to distributed implementations[END_REF], that we call here strong conflict. Two interactions are strongly conflicting iff they share a common port or there is a couple of ports with one member in each interaction such that these two ports label two conflicting transitions of the same component. Clearly, strong conflict implies weak conflict but the converse is not true.

DEPRIORITIZING A BIP MODEL

In this section, we describe our approach to transform a BIP model B into an equivalent model without priorities, denoted B. Intuitively, our transformation proceeds as follows:

1. First, it replaces atomic components in B by functionally equivalent send/receive atomic components, where atomicity of transitions and interactions is broken. This first transformation, already used in [START_REF] Basu | Distributed semantics and implementation for systems with interaction and priority[END_REF][START_REF] Bonakdarpour | Automated conflict-free distributed implementation of component-based models[END_REF][START_REF] Bonakdarpour | From high-level component-based models to distributed implementations[END_REF] separates the synchronization from computation on component transitions and enables the concurrent execution of atomic components.

2. Secondly, it inserts manager components for handling interactions. These managers detect enabledness of interactions and schedule them for execution according to priority rules. Managers interact with each other through multi-party interactions in order to maintain a consistent view on the state of the system.

Breaking Atomicity

The transformation of atomic components splits each transition into two consecutive steps: (i) an offer that publishes the current state of the component, and (ii) a notification that triggers the update function. The intuition behind this transformation is that the offer transition corresponds to sending information about component's intention to interact with the other components. The notification transition receives the response from the scheduler, once some interaction has been completed. Local update functions can then be executed concurrently and independently by components upon notification reception.

The offer transition publishes its enabled ports through a special port named o. Enabled ports are encoded through a list of Boolean variables. After the computation of the local function, this list is updated to the ports that are enabled at the next control state. Notification transitions are triggered by corresponding ports from the original atomic component.

Definition 5 (Transformed atomic components).

Let B = (Q, X, P, T ) be an atomic component. The corresponding transformed atomic component is

B ⊥ = (Q ⊥ , X ⊥ , P ⊥ , T ⊥ ), such that: • Q ⊥ = Q ∪ {⊥s |s ∈ Q}.
• X ⊥ = X ∪ {xp}p∈P , where each xp is a Boolean variable indicating whether port p is enabled.

• P ⊥ = P ∪ {o}, where o is the offer port. All variables in X ⊥ are associated to o (i.e., Xo = X ⊥ ). • For each transition τ = (q, p, g, f, q ′ ) ∈ T , we include the following two transitions in T ⊥ :

1. offer τ q o = (⊥q, o, go, fo, q) where go is true, fo is the identity function, and 2. notification τ q p = (q, p, gp, fp, ⊥ q ′ ) where gp is true and fp applies fτ on X and for each port r ∈ P , it sets xr to true if τ ′ = (q ′ , r, g ′ , f ′ , q ′′ ) ∈ T for some q ′′ and g ′ is true. Otherwise, xr is set to false.

In Definition 5, states {⊥s |s ∈ Q} from where the component sends offers, are called busy or unstable states. States Q, from where the component is waiting to receive a notification, are called stable states.

Figure 3 shows the transformed version of the atomic component shown in Figure 2(a). Initially, the component is in busy state ⊥s and the value of xp is true; i.e., the component is willing to interact on port p. Then, it sends an offer through port o containing the current values of xp and n and reaches stable state s. The reception of a notification corresponds to the p-labelled transition that brings back the component to the initial busy state.

Interaction Managers

The set of managers are introduced to execute interactions according to the global semantics of the original BIP model described in Section 2. To this end, a manager component for an interaction a has to (i) detect enabledness of a by listening to offers sent by atomic components, (ii) trigger the execution of a, (iii) notifies atomic components as well as the other conflicting managers, whenever the interaction is executed.

Let us observe that if two interactions are weakly conflicting, then executing one can change the status of the other. For instance, let a and b be two interactions, such that a ⊕ b; i.e., they share some component B. Obviously, executing a triggers a transition in component B. This transition can result in changing the status of interaction b. That is, until component B completes its local execution and sends a new offer, the status of enabled ports and values of variables in B can change.

Definition 6 (Interaction Manager). Let a ∈ γ be an interaction, where Pa = {pi}i∈I . The interaction manager Ma is an atomic component Ma = (Q, X, P, T ) defined as follows:

• The set of control states is Q = {undef , en, dis, exc}.
Intuitively, in state undef (undefined), the manager does not have enough information to decide whether or not interaction a is enabled. This is normally because some offers have not been received yet. In states en (enabled) and dis (disabled), the manager knows that a is enabled or disabled, respectively. In state exc (executing), the interaction a is being executed.

• The set of variables is X = {b a i }i∈I ∪ {x a p i }∪X a p i p i ∈a . For every component Bi, the manager holds a Boolean variable b a i which is true iff component Bi is in a stable state, that is, waiting for a notification. For every port pi ∈ a, the manager holds respectively, a Boolean x a p i which indicates the status of the port (i.e., enabled or disabled) and variables X a p i that is, data associated to the port pi.

• The set of ports P and their associated variables is presented in Table 1.

• The set of transitions T and their intuitive meaning is presented in Table 2.

undef dis dis a exc en n a ⊕ a o i b a i := true ⌈(b a 1 ∧ x a p 1 )∧ (b a 2 ∧ x a p 2 )∧ (n a 1 > n a 2 )⌋ ι ⊕ a ⌈(b a 1 ∧ b a 2 )∧ (¬x a p 1 ∨¬x a p 2 ∨n a 1 ≤ n a 2 )⌋ ι ⊕ a ⊕dis a start a b a 1 := f alse b a 2 := f alse sw(n a 1 , n a 2 ) x a p 1 n a 1 o a 1 x a p 2 n a 2 o a 2 n a b a 1 b a 2 ⊕ a ⊕dis a dis a start a
∧ x a p 1 ) ∧ (b a 2 ∧ x a p 2 ) ∧ (n a 1 > n a
2 ) switches from undefined to enabled state. The two first conjuncts ensures that (1) B1 and B2 are in stable state, and (2) p1 and p2 are enabled. The latter conjunct corresponds to the guard of interaction a in Figure 2(b). Likewise, the transition from undef to dis allows reaching the state where the interaction a is disabled. The update function associated to τstart sets b a 1 and b a 2 to false and then swaps the variables n a 1 and n a 2 . Both n a 1 and n a 2 are associated to the notification port na, so their new values are sent back to the component.

Connecting Managers

The transformed atomic components and interaction managers are interconnected using three types of interactions: (i) offer interactions where components send their enabled ports to corresponding managers, (ii) notification interactions where managers notify components after execution of an interaction, and (iii) schedule interactions where priority rules are handled.

We now formally define the deprioritized model, by specifying how we connect the components defined so far. Let γ (i) denote the set of all interactions in γ that involve the component Bi.

Definition 7 (Deprioritized model). Given a model

B = πγ(B1, • • • , Bn), with γ = {a1 • • • am}, we define its de- prioritized version as B = γ(B ⊥ 1 , • • • , B ⊥ n , Ma 1 , • • • , Ma m ), where B ⊥
i is obtained from Bi as explained in definition 5, Ma j is obtained from aj as explained in definition 6, and γ contains the following interactions:

• Offer interactions.

For each i ∈ {1 • • • n}, γ contains the interaction offi, where

P of f i = {oi} ∪ a∈γ (i) {o a i }.
For each interaction a ∈ γ (i) , the update function F of f i sets the values of variables {x a p i } ∪ X a p i to the values of {xp} ∪ Xp associated to oi, where p is the of port Bi involved in a. Offer interactions have no guard and they only copy data from the sender component to the manager.

• Notifications interactions.

For each interaction a ∈ γ, where a = {pi}i∈I , γ contains the interaction nota, such that Pnot a = na ∪ {pi}i∈I . This interaction notifies each component which port has been selected. The update function Fnot a copy back data to each component Bi involved in a. That is, the values of Xp i (in Bi) are set to the values of X a p i (from Ma). • Schedule interactions. For each interaction a ∈ γ, γ contains the interaction ã:

Pã = {starta} ∪ {⊕c|c ⊕ a, c > a} ∪ {disc|c ⊕ a, aπc} ∪ {⊕disc|c ⊕ a, aπc}
This interaction has no guard. For each interaction c weakly conflicting with a, the update function Fã sets variable b c i of manager Mc to false through the port ⊕c if {a, c} ⊆ γ (i) . In other terms, the starta interaction informs the manager Mc that the components causing the weak conflict with a have moved and are not in their stable state anymore. This information maintains coherence between the b c i variable in each manager Mc and the actual state of component Bi. Figure 5 presents the deprioritized model from Figure 2(b). Please note that the port names have been shortened for space reasons (e.g. sa and da stand for starta and disa respectively. For offer and notification interactions, we interpret a triangle port as a send port (i.e., for sending offers) and bullet port as a receive port (i.e., for receiving offers). Note that offers and notifications only copy variables between components and managers.

o 1 p 1 B ⊥ 1 o 2 p 2 B ⊥ 2 o 3 p 3 B ⊥ 3 o 4 p 4 B ⊥ 4 M a n a o a 1 o a 2 ⊕ a ⊕ d a d a s a M b n b o b 2 o b 3 ⊕ b ⊕ d b d b s b M c n c o c 3 o c 4 ⊕ c ⊕ d c d c s c ã b c
If we assume priorities bπa and cπa for the model in Figure 2(b), we obtain the following schedule interactions: a has no higher priority interaction and is weakly conflicting with b, thus Pã = {starta, ⊕ b }. Executing ã will set the variable b b 2 to false in M b , since B2 will become busy. b has less priority than a and is weakly conflicting with both a and c, thus Pb = {start b , ⊕disa, ⊕c}. c has less priority than a and is weakly conflicting with b, thus Pc = {startc, disa, ⊕ b }.

Correctness

We now show that the above transformation preserves the semantics of the original BIP model. By preserving the original semantics, we mean ensuring observational equivalence between the original model and the transformed model. This is proved in Theorem 1.

Let

B = πγ(B1, • • • , Bn) be a BIP model and B = γ(B ⊥ 1 , • • • , B ⊥ n , Ma 1 , • • • , Ma m ) be its unprioritized version. We denote q = (q1, • • • , qn) a state of B and q = (q1, • • • , qn, s1, • • • , sm) a state of B. We show that B is observationally equivalent to B.
The observable actions of B are the interactions γ. The observable actions of B are only the schedule interactions, that is {ã|a ∈ γ}. The remaining interactions in B, namely offers offi and notifications nota, are unobservable and are denoted β. We denote q β → q′ if a β action brings the system from state q to state q′ . Proof. In any reachable state, if a manager reaches the state exc then the corresponding notification is enabled, since schedule interactions and boolean variables bi ensure that each component may receive only one notification after each offer. Similarly, if any component reaches an unstable state, then the corresponding offer is enabled.

Offer interactions are independent since they do not share any port nor change a common variable. Thus, the order of their execution does not change the final state.

Notification interactions (that correspond to interactions of the original model, augmented by a notification port) enabled from a reachable state are not conflicting since schedule interactions handle weak conflicts. Thus, notification interactions are independent and their order of execution does not change the final state. We can conclude that β → is confluent.

From proposition 1 and 2, for each reachable state q of B, there is a unique state denoted [q] such that q β * → [q] and

[q] β →.
We recall the definition of observational equivalence of two transition systems A = (QA, P ∪ {β}, →A) and B = (QB, P ∪ {β}, →B). It is based on the usual definition of weak bisimilarity [START_REF] Milner | Communication and concurrency[END_REF], where β-transitions are considered unobservable. The same definition is trivially extended for atomic and composite BIP components.

Definition 8 (Weak Simulation). A weak simulation over A and B, denoted

A ⊂ B, is a relation R ⊆ QA × QB, such that we have ∀(q, r) ∈ R, a ∈ P : q a →A q ′ =⇒ ∃r ′ : (q ′ , r ′ ) ∈ R ∧ r β * aβ * → B r ′ and ∀(q, r) ∈ R : q β →A q ′ =⇒ ∃r ′ : (q ′ , r ′ ) ∈ R ∧ r β * →B r ′
A weak bisimulation over A and B is a relation R such that R and R -1 are both weak simulations. We say that A and B are observationally equivalent and we write A ∼ B if for each state of A there is a weakly bisimilar state of B and conversely. We consider the correspondence between observable actions of B and B as follows. To each interaction a ∈ γ, where γ is the set of interactions of B, we associate the schedule interaction ã of B.

Theorem 1. B ∼ B.
Proof. We define the relation R between the states of B and the states of B as follows: the couple (q, q) is in the relation R if the states of atomic components

B ⊥ 1 , • • • , B ⊥ n in [q]
are the same as in q. Formally, we have (q, q) ∈ R if

[q] = (q1, • • • , qn, s1, • • • , sm) and q = (q1, • • • , qn).
We show that R is an observational equivalence by proving the next three assertions:

(i) If (q, q) ∈ R and q β → r then (r, q) ∈ R.

(ii) If (q, q) ∈ R and q ã → r then ∃r : q a → r and (r, r) ∈ R.

(iii) If (q, q) ∈ R and q a → r then ∃r : q β * ã -→ r and (r, r) ∈ R.

The point (i) comes from the definition of R.

(ii) If the interaction ã is enabled, then manager Ma is in state en, which implies that at equivalent state q:

• All ports of a are enabled and the guard Ga is true, since the guard of the τen transition is true

• No higher priority interaction is enabled since ã is enabled only when managers corresponding to such interactions are in state dis.

Thus we have q a → r, and the reader can easily check that (r, r) ∈ R.

(iii) From q we can reach [q] by using only β transitions. In state [q], since every atomic component has sent an offer, the state of each manager will be either en or dis, according to the status of the corresponding interaction at state q in B. Then since a is enabled at state q, Ma is in state en at state [q]. If there is any interaction b with higher priority than a, then it is disabled in state q, thus the manager M b is in state dis at state [q]. Thus ã is enabled at state [q] and we have q β * ã -→ r. Executing the notification interaction na and the offer interactions from components involved in a lead B in a state where atomic components have the same state as in r. Thus (r, r) ∈ R.

BUILDING A DISTRIBUTED MODEL: THE 3-TIER ARCHITECTURE

Once we construct a model with no priorities as prescribed in Section 3, one can apply the technique presented in [START_REF] Bonakdarpour | Automated conflict-free distributed implementation of component-based models[END_REF] to generate distributed code. We now briefly recap this technique. The code generation is accomplished in two steps. First, from a given BIP model, we generate another BIP model that only incorporates asynchronous message passing as interactions (denoted SR-BIP). Then, we transform the SR-BIP model into a set of executables -one per atomic component -that communicate using asynchronous message passing primitives such as MPI or TCP sockets. We only review the first step.

Distributed execution of interactions may introduce conflicts even if we do not consider priorities. Thus, our target SR-BIP model in a transformation should have the following three properties: (1) preserving the behaviour of each atomic component, (2) preserving the behaviour of interactions, and (3) resolving conflicts in a distributed manner. Moreover, we require that interactions in the target model are asynchronous message passing.

We design our target BIP model based on the three tasks identified above, where we incorporate one tier for each task. Since several distributed algorithms exist in the literature for conflict resolution, we design the tier corresponding to conflict resolution so that it provides appropriate interfaces with minimal restrictions. As a running example, we use the part of the model presented in Figure 5 the transformation explained in Subsection 3.1 to break atomicity), and

• B ⊥ 1 • • • B ⊥
n are copied from the deprioritized model, since they have already been transformed by the deprioritization.

Interaction Protocol. This tier consists of a set of components each hosting a set of interactions from the deprioritized BIP model. Conflicts between interactions included in the same component are resolved by that component locally. For instance, interactions ã and b in Figure 5 are grouped into component IP1 in Figure 6. Thus, the conflict between ã and b is handled locally in IP1. To the contrary, the conflict between b and c has to be resolved using the third tier of our model. The interaction protocol also evaluates the guard of each interaction and executes the code associated with an interaction that is selected locally or by the upper tier. The interface between this tier and the component tier provides ports for receiving enabled ports from each component and notifying the components on permitted port for execution (ports nã, nb, nc).

Conflict Resolution Protocol.

This tier accommodates an algorithm that solves the committee coordination problem [START_REF] Chandy | Parallel program design: a foundation[END_REF] to resolve conflicts between interactions hosted in separate interaction protocol components. For instance, the external conflict between interactions b and c is resolved by the central component CRP in Figure 6. We emphasize that the structure of components in this tier solely depends upon the augmented committee coordination algorithm. Incorporating a centralized algorithm results in one component CRP as illustrated in Figure 6. Other algorithms, such as ones that use a circulating token [START_REF] Bagrodia | A distributed algorithm to implement n-party rendevouz[END_REF] or dining philosophers [START_REF] Bagrodia | Process synchronization: Design and performance evaluation of distributed algorithms[END_REF][START_REF] Chandy | Parallel program design: a foundation[END_REF] result in different structures in this tier and are discussed in detail in [START_REF] Bonakdarpour | From high-level component-based models to distributed implementations[END_REF]. The interface between this tier and the Interaction Protocol involves ports for receiving requests to reserve an interaction (labelled r) and responding by either success (labelled ok) or failure (labelled f ).

CASE STUDY

In this section, we use a jukebox example to illustrate our deprioritization transformation and conduct experiments to study the effectiveness of our method (see the models in Figures 7 and8). This model represents a system, where a set of readers (R1, . . . , R4) need to access the data located on discs (D1, . . . , D4). A reader may need any disc. Access to the disc is managed by jukebox components (J1, J2) that can load any disc to make it available for reading. Each pair (Di, J k ), i ∈ {1 • • • 4} and k ∈ {1, 2}, has two interactions: (1) a load i,k interaction for loading the disc in the jukebox and an unload i,k interaction for unloading it. Each reader Rj is connected to a jukebox through a readj interaction. During the test, we simulate execution of interactions by waiting a given amount of time. Namely, we wait 100ms for load /unload and 500ms for read.

Figure 8 presents the behaviour of atomic components and the data transfer on interactions. To ensure that all discs are eventually loaded, each jukebox keeps a list of discs to load, namely to load. Each time a disc is loaded, it is removed from the list by the load transition in the jukebox component. The guard of a load interaction prevents the disc to be loaded if it is not on the list. When the to load list becomes empty, it is reinitialized to the set of all discs on the unload interaction. The variable current contains the identity (i.e., 1 . . . 4) of the disc currently loaded in the jukebox, and is updated by the load interaction. In order to ensure that the reader gets the correct data, a guard on the {read , data} interaction holds, only if the disc in the jukebox (current) is the one to be read (to read ). Each reader has a sequence of 2 discs to read. The variable to read contains the id of the next disc to be read. It is initialized with the first value (not shown in the figure), and is updated after the first read. This model has finite runs: the execution terminates when all readers have read the two discs they needed.

We consider two versions of the model. The first model, denoted B ∅ , does not contain priorities. The second model, denoted Bπ, is the B ∅ restricted by two types of priorities:

• Priorities to enforce termination. We give priority to the read interactions over the unload interactions. Formally, it corresponds to the sets of priorities {unloadi,1 π readj 4}}, for each jukebox. This ensures that any enabled read interaction will be executed before the disc is unloaded. Since each disc is eventually loaded, each read interaction will take place and the execution terminates. Otherwise, sequences of load /unload interaction could occur forever. Note that here, assuming fairness ensures that the model eventually terminates.

| i ∈ {1, • • • , 4}, j ∈ {1, 2}} and {unloadi,2 π readj | i ∈ {1, • • • , 4}, j ∈ {3,
• Priorities to speed up execution. By inspecting the discs requested by the readers, we know that some discs are more often needed than others. Thus, we give higher priority to the corresponding load interactions. Here, we give higher priority to Disc 1 in Jukebox 1 by adding the following set of priorities: {loadi,1 π load1,1 | i ∈ {2, 3, 4}}.

For both versions B ∅ and Bπ, we generate the corresponding deprioritized models B∅ and Bπ. In Table 3, we present the size -the number of atomic components and the number of interactions -of these different models, in the columns labelled "Orig.". We then apply the transformation provided in [START_REF] Bonakdarpour | From high-level component-based models to distributed implementations[END_REF] to the models B ∅ , B∅ , and Bπ to obtain a distributed version of each model including a centralized scheduler2 . The number of Send/Receive components and interactions contained in the distributed version of these models is given in the columns labelled "S/R" in Table 3. We simulate the execution of these models on two different platforms. The first one is centralized, where only one processor is available to execute all components. The second one is fully decentralized, where each atomic component has its own processor. We assume that executing a load, unload or read interaction completely blocks the processor. For each couple (model, execution platform), we measure the average time of terminating executions. The results are presented in Table 3.

As mentioned earlier, we applied our deprioritization transformation to model B ∅ although we can directly obtain a distributed model. By comparing the execution times of B ∅ and B∅ on the centralized platform, we observe that our deprioritization transformation does not introduce a significant overhead, even if it increases the number of components and interactions. More importantly, the distributed execution of B∅ is almost twice faster than B ∅ . This is due to the fact all time consuming computations in B ∅ are on interactions, which are all executed on the same processor (the one hosting the scheduler). When switching to B∅ , these interactions are executed by the manager components and, hence, run concurrently on different processors.

Furthermore, the model Bπ runs faster than B ∅ on a centralized platform. In this scenario, priorities enforce a better scheduling -we first load the discs that are often used and we do not perform an unload if a reader has something left to read -and thus reduce the total execution time. Again, switching to decentralized execution gives almost twice better results, as (time consuming) interactions are now running concurrently.

CONCLUSION

In this paper, we proposed an automated method to derive correct distributed implementation from high-level component-based models encompassing prioritized multiparty interactions. Our method consists of three steps: (1) one transformation to deprioritize the initial model, (2) a transformation from [START_REF] Bonakdarpour | Automated conflict-free distributed implementation of component-based models[END_REF][START_REF] Bonakdarpour | From high-level component-based models to distributed implementations[END_REF] that generates a distributed model from the deprioritized model by resolving interaction conflicts, and (3) a final transformation from the distributed model into C++ code. All steps preserve observational equivalence between the input and output models. We illustrated our approach using a non-trivial example and presented encouraging experimental results.

There exist several research directions for future work. First, more rigorous and deeper case studies and experiments are needed to completely understand the overheads introduced by our transformations. Since deprioritization is an independent step of our method and is isolated from conflict resolution (i.e., step two), one can study the overhead of each step separately. Another direction is to devise a committee coordination algorithm for conflict resolution that takes priority issues into account. This allows us to incorporate such an algorithm directly in our 3-tier model [START_REF] Bonakdarpour | From high-level component-based models to distributed implementations[END_REF]. This approach can potentially have less overhead than the one presented in this paper. Finally, one can speed-up distributed execution of models with priorities by detecting disabled interactions as early as possible. Such detection can benefit from knowledge-based methods (e.g., [START_REF] Basu | Priority scheduling of distributed systems based on model checking[END_REF]).

  sorts integers ni, obtained by gluing 4 atomic components using 3 interactions.

Figure 2 :

 2 Figure 2: Atomic and composite components in BIP

Figure 2 (

 2 b) illustrates a composite component γ(B1, • • • , B4), where each Bi is identical to component B in Figure 2(a). The set γ of interactions is {a, b, c}, where a = ({p1, p2}, n1 > n2, sw(n1, n2)) and function sw swaps the values of its arguments. Interactions b and c are defined in a similar fashion. Interaction a is enabled when ports p1 and p2 are enabled and the value of n1 (in B1) is greater than the value of n2 (in B2). Thus, the composite component B sorts variables n1 • • • n4, such that n1 contains the smallest and n4 contains largest values.

Figure 3 :

 3 Figure 3: Transformed version of one atomic component from Figure 2(b)

Figure 4 :

 4 Figure 4: The manager component for interaction a between components B1 and B2 in Figure 2(b).

Figure 5 :

 5 Figure 5: Deprioritized version of model from Figure 2(b).

→

  is terminating.Proof. Each β action involve at least a component. Each component can take part in at most 2 β actions, 1 notification and 1 offer, then no other β action is possible until an ã action is executed. Thus at most 2n consecutive β-steps are possible. Proposition 2. From any reachable state q of B, β → is confluent.

Figure 6 :

 6 Figure 6: Distributed version of the deprioritized model from Figure 5.

4 Figure 7 :

 47 Figure 7: BIP Model for the jukebox example.

Figure 8 :

 8 Figure 8: Behaviour of jukebox components and interactions.

Table 1 :

 1 Ports of a manager component

	port	variables	description
	o a i	{x a p i } ∪ X a p i	receives offers from atomic compo-nent Bi
	ι	∅	change status to enabled or disabled (internal port)
	starta	∅	triggers interaction execution
	na	{X a p i }	notifies atomic components upon execution
	disa	∅	signals disabled status to other man-agers
			gets notified about execution of
	⊕a	{b a i }	weakly conflicting interactions by
			other managers
	⊕disa	{b a i }	similar to port ⊕a, but for interac-tions with higher priority

Table 2 :

 2 Transitions of a manager component

	Transition	Guard / Function	Description
	undef	o a i -→ undef	-/ b a i := true	receive offer from Bi
	undef	ι -→ en	Ga ∧ ∀i ∈ I. (b a i ∧ x a p i ) / -	change state to enabled
	undef	ι -→ dis	(∀i ∈ I.b a i ) ∧ (¬Ga ∨ ∃i ∈ I.¬x a p i ) / -	change state to disabled
	en	starta → exc	-/ {b a i } := f alse; {X a p i }:=Fa({X a p i })	execute interaction, apply update function.
	exc		na → undef	-/ -	notifies atomic components on execution
	dis	disa → dis	-/ -	signals disabled state
	dis undef ⊕a → undef ⊕a → undef en → undef ⊕a	-/ -	gets notified about execution of a weakly conflicting interaction
	dis	⊕disa → undef	-/ -	gets notified about execution of a higher priority weakly conflicting interaction
	Figure 4 represents the manager for interaction a in Figure
	2(b). It contains the variables b a 1 and b a 2 since interaction a
	involves components B1 and B2. The manager contains two
	offer ports o a 1 and o a 2 . Port o a i , i ∈ {1, 2}, is associated with variables (1) x a p

i , which indicates the status of port pi in Bi, and (2) n a i , that are local copies of variables ni associated to ports pi in Figure 2(b). All these variables are refreshed upon receiving an offer through ports o a i . The transition from undef to en guarded by (b a 1

Table 3 :

 3 Model size and execution time (s) for different implementations of Figure7.

		Model Size # Atoms # Interactions	Execution time
		Orig. S/R Orig.	S/R	Cent. Decent.
	B ∅	10	11	20	28	15.2	11.0
	B∅	30	31	70	148	12.0	5.9
	Bπ	30	31	70	154	5.4	2.8

Although our focus is on BIP, all results in this paper can be applied to any model that is specified in terms of a set of components synchronized by broadcast and rendezvous interactions.

We cannot transform directly Bπ into such a distributed model since the transformation presented in[START_REF] Bonakdarpour | From high-level component-based models to distributed implementations[END_REF] does not take priorities into account.

* The research leading to these results has received funding from the European Community's Seventh Framework Programme [FP7/2007-2013] under grant agreement no 248776 (PRO3D) and no 257414 (ASCENS), from ARTEMIS JU grant agreement ARTEMIS-2009-1-100230 (SMECY) and from Canada ORF RE03-045, NSERC DG 357121-2008, and IS09-06-037 grants.


Ioannis Emiris

Angelos Mantzaflaris

Bernard Mourrain

Voronoï Diagrams of Algebraic Distance Fields

Keywords: Voronoï diagram, subdivision algorithm, bisector curve, anisotropic diagram, lower envelope

We design and implement an efficient and certified algorithm for the computation of Voronoï Diagrams (VD's) constrained to a given domain. Our framework is general and applicable to any VD-type where the distance field is given explicitly or implicitly by a polynomial, notably the anisotropic VD or VD's of non-punctual sites. We use the Bernstein form of polynomials and DeCasteljau's algorithm to subdivide the initial domain and isolate bisector, or domains that contain a Voronoï vertex. The efficiency of our algorithm is due to a filtering process, based on bounding the field over the subdivided domains. This allows to exclude functions (thus sites) that do not contribute locally to the lower envelope of the lifted diagram. The output is a polygonal description of each Voronoï cell, within any user-defined precision, isotopic to the exact VD. Correctness of the result is implied by the certified approximations of bisector branches, which are computed by existing methods for handling algebraic curves. First experiments with our C++ implementation, based on double precision arithmetic, demonstrate the adaptability of the algorithm.

Introduction

Voronoï Diagrams (VD's) have a surprising variety of applications, e.g. in path planning, computer vision and machine perception, meshing etc. It is a widely studied subject, and several algorithms exist for their computation.

The VD of a given set of geometric objects (sites) in the plane is the partition of the plane into regions (cells), where each site S is associated to the region consisting of all points for which S is the nearest site, compared to any other site of the (so called) generating set. If we take a set of points in the plane as generating set, and as distance between a site and an arbitrary point the length of their connecting line segment, we have the classic VD under the Euclidean metric. The distance induces a scalar distance field over the plane for every site in the generating set.

There are at least two ways to generalize the construction of VD's, in order to meet applications' needs; to allow non-punctual sites, e.g. circles, under the Euclidean distance field, or to attach distance fields other than the Euclidean to the sites, for example the anisotropic diagram (see Sect. 1.1). These generalizations lead to curved VD's, having algebraic bisectors, see [START_REF] Boissonnat | Curved Voronoï diagrams[END_REF] for an expository paper.

The bibliography on VD's is vast, and many alternative strategies have been proposed, but let us mention some of the latest developments. Boada et al. [START_REF] Boada | Approximations of 2d and 3d generalized voronoi diagrams[END_REF] use a subdivision approach to compute approximate VD's. The closest site to the corners of a subdivided domain is computed and used to deduce the Voronoï cell in which the box belongs to. Their output is proved to converge to the exact VD, but the result is not always topologically correct, since corner signs alone are not in one-to-one correspondence with the possible configurations of bisectors inside rectangular domains. A similar idea, using the centres of pixels as points of reference and graphics hardware is implemented in [START_REF] Hoff | Fast computation of generalized Voronoï diagrams using graphics hardware[END_REF]. Their discrete computation has to be carried out on every pixel of the screen, and is subject to resolution errors.

Setter et al. [START_REF] Setter | Constructing twodimensional Voronoï diagrams via divide-and-conquer of envelopes in space[END_REF] compute VD's using divide-and-conquer of lower envelopes and exact computations. They recursively build the diagram, by going down to two sites and then merge up a full VD. They provide good treatment of degenerate cases, and a complexity analysis for randomized inputs. They comment that computation speed may be limited due to extensive use of symbolic computations.

Divide-and-conquer techniques are also utilized by Aichholzer et al. [START_REF] Aichholzer | Divide & conquer algorithms for Voronoi diagrams revisited[END_REF], to propose a VD algorithm for general shapes. They exploit connections to medial axis and use a plane-sweep technique. Their method avoids computing redundant pieces of bisectors that will then be rejected.

Seong et al. [START_REF] Seong | Voronoï diagram computations for planar nurbs curves[END_REF] explore the VD of parametric NURBScurves. They approximate bisectors as algebraic curves in parameter-space and use real solving to find their intersections globally, and then trim away the unwanted parts.

Emiris et al. [START_REF] Emiris | Exact Delaunay graph of smooth convex pseudo-circles: general predicates, and implementation for ellipses[END_REF] present an efficient, certified algorithm for VD's, specialized to certain families of closed planar curves, notably ellipses. They use adapted numerical techniques to reduce the degree of predicates to be evaluated. Their algorithm becomes exact if they choose to use a special iterated resultant for the "InCircle" predicate.

Algebraic curves are the main objects of study when looking at VD's and their bisector sets. Meshing algebraic curves is regarded as a black-box operation in the present work. Several subdivision techniques are available for this task. An important highlight of certain approaches is that the output is guaranteed to be isotopic to the input curve. In [START_REF] Liang | Subdivision methods for the topology of 2d and 3d implicit curves[END_REF] they give algorithms for (local and global) isotopic computation of algebraic curves in 2D and 3D. They use univariate solving, enveloping techniques and fast regularity criteria. In [START_REF] Mantzaflaris | A subdivision approach to planar semi-algebraic sets[END_REF], these ideas are extended to semialgebraic sets, notably to the arrangement of several curves.

Our method computes VD's by means of the lower envelope of the distance fields. The latter exhibit an algebraic degree which is usually quite low, thus handling them in an algebraic fashion is computationally advantageous. The method is inspired by, and uses tools from geometric modeling and CAGD (Computer Aided Geometric Design). Polynomial curves constitute a major branch of the research in CAGD, and flexible curve representations have been developed in this frame. Therefore, we aim at bringing data representations and algorithmic experience from CAGD to this fundamental problem of computational geometry.

Our algorithm applies on the distance field induced by the sites and does not require the sites themselves in the input. It provides polygonal approximations of the Voronoï cells and is applicable to all inputs where the distance field of every site can be expressed by a polynomial. This expression may be explicit, but also implicit. In Sect. 1.1 we give a list of distance fields for commonly encountered VD's.

The implementation is done in C++, using double precision with controlled rounding modes to certify the subdivision process. The algorithm is naturally parallelizable, since it applies locally and independently on different domains of the plane. We did not focus on insertion of new sites, or point queries, yet these operations can easily be added, without significant changes. On the other hand, we focus on genericity of the framework, the correctness of the result and the ability to treat arbitrary diagrams, if certain reasonable algorithmic prerequisites are met (cf. Sect. 2).

The rest of the text is organized as follows. We complete our introduction by a list of common VD's, as well as some details on curve representations. We present the core of our algorithm in Sect. 2, with details on the filtering, subdivision and the recovery phase. Then we extend the method to implicitly given distance fields in Sect. 3, and conclude with some experiments in Sect. 4.

Voronoï diagrams & distance fields.

In this section we recall a non-exhaustive list of VD-types that we are interested in, and briefly introduce the main tools and representations that build up our algorithm.

For a point q = (x, y), the distance between q and the Voronoï site attached to f i , as viewed by site i, is given by dist i (q) = f i (q). Now one can define the Voronoï cell of i-th site as:

Vor(i) = {q ∈ R 2 : dist i (q) ≤ dist j (q), j = 1 . . . n}
Therefore the VD can be retrieved as the projection of the lower envelope of the distance fields [START_REF] Edelsbrunner | Voronoï diagrams and arrangements[END_REF], also known as the minimization diagram.

In the simple case of Euclidean VD the input is a list of (squared) distance functions (f 1 , . . . , f n ). The lower envelope changes if we work with squared fields, yet the VD (that is, its projection) remains intact under squaring, or under any other invertible and strictly increasing transformation. Different distance fields dist i (q), q ∈ R 2 , give rise to different types of VD's, including the: -VD of points p i = (x i , y i ) under the p -metric, p even:

dist i (q) = (x -x i ) p + (y -y i ) p . -Anisotropic diagram of points p i = (x i , y i ) with weights w i ∈ R [11]: dist i (q) = (q -p i ) T M i (q -p i ) -w i , with M i ∈ R 2×2 symmetric positive-definite matrices. -Power (or Laguerre) diagram of points with weights: dist i (q) = q -p i 2 -w 2 i . -Möbius diagram of points: dist i (q) = v i q -p i 2 -w i , with v i , w i ∈ R.
-Apollonius (or additively weighted) diagram of disks [START_REF] Emiris | The predicates of the apollonius diagram: Algorithmic analysis and implementation[END_REF] with centers p i and radii w i : dist i (q) = qp iw i . -Euclidean VD of ellipses [START_REF] Emiris | Exact Delaunay graph of smooth convex pseudo-circles: general predicates, and implementation for ellipses[END_REF], or Euclidean VD of general closed parametric curves [START_REF] Gravesen | Curves and surfaces represented by polynomial support functions[END_REF]. An algebraic function f can be represented over a bounded domain D = [a, b] × [c, d] by its Bernstein coefficients over D. A number of properties of this basis, e.g. convexity, variation diminishing, positivity etc, make it suitable for stable numerical computations. Also, DeCasteljau's algorithm can be applied to split the representation, i.e. produce the coefficients over sub-domains of D (cf. [START_REF] Farin | Curves and surfaces for CAGD[END_REF]). We use the notation f i | D for the restriction of f i (x, y) in D, that is, the Bernstein representation of f i over the domain D. For instance, for the anisotropic diagram, the f i 's are conics, thus every f i | D is defined over D by 9 coefficients.

Note that some VD types in the previous list do not have polynomial distance fields, nor do they admit squaring to eliminate radicals. In such cases, it is usually possible to express the distance field implicitly. For instance, if we want to use our framework to compute the Apollonius diagram, i.e. the Euclidean VD of circles with centers p i and radii w i , we need to treat the corresponding distance function dist i (q) = qp iw i , which may be transformed to the algebraic equation (dist i (q) + w i ) 2 = qp i 2 . Therefore, distance field values z := dist i (x, y) are given in implicit form by the cone

g i (x, y, z) := (x -x i ) 2 + (y -y i ) 2 -(z + w i ) 2 = 0, (1)
i.e., given (x 0 , y 0 ), values of dist i (x 0 , y 0 ) are computed by solving g i (x 0 , y 0 , z) = 0 for the smallest root > 0 w.r.t. z. This tri-variate function can be represented similarly over a cuboid, by tensoring 3 Bernstein bases. In Sect. 3, we shall adapt our algorithm to use implicit representations.

The algorithm

Our algorithm consists of two phases: The subdivision phase (Alg. 1), followed by the reconstruction phase (Alg. 2). Two algorithmic ingredients must be provided in order to apply the method on a specific VD type:

(i) Field bound: A way to bound the value range of a distance field over a given domain. (ii) Bisector tracking: A way to compute an approximation of bisector(s) or vertices in a given domain. These computations need to be carried out on axis-aligned boxes of the subdivision. Note that at this level we do not emphasize on efficiency, e.g. the bounds could be bad, or the approximation very loose. Nevertheless, it is expected that these computations converge to the actual distance value or bisector, when the size of the boxes becomes smaller.

The first box that is computed as soon as the algorithm is launched is the one corresponding to the initial domain D 0 , carrying all distance fields of the input, in Bernstein form.

In the subdivision phase we compute a graph of boxes that span the VD. To do so, the field bound is used in a filtering process in order to exclude most boxes and reach down to boxes intersecting the VD. These boxes contain Bernstein representations of sites that contribute to the part of the diagram inside the box.

In the reconstruction phase the boxes of the subdivision are traversed and the part of the VD that intersects the box is meshed locally using the second ingredient. These bisector segments are finally stitched together to recover the full VD.

Subdivision phase.

The idea of a space-subdivision scheme is typical in CAGD. A big domain of interest is divided into smaller ones until some (fast computed) conditions are fulfilled. Then every small item is treated independently.

We denote by

B(D) = {f i | D , f j | D , .
. . } the list of control grids that are held for a subdivided domain D. Also, by signature of a box sig(B(D)) we mean the site labels present in the list, i.e. sig(B(D)) = {i :

f i | D ∈ B(D)}.
The main loop in Alg. 1 starts by filtering the box B(D). This allows to remove large valued f i 's, i.e. those that correspond to sites that are "away" from the box. Then, we check if the remaining fields are at most three. In this case the box is not subdivided anymore. Otherwise the box is split in two sub-domains, along the longest of its sides. The new boxes are pushed in the stack and the next loop begins.

If B(D) is left with only one f i | D , then D ⊆ Vor(i), hence this box does not span the VD. If there are 2 active functions then the bisector locus f i | Df j | D = 0 possibly intersects the box. If there are 3 or more functions then there probably exists a Voronoï vertex in D. Boxes that contain degenerate Voronoï vertices, i.e. vertices that are equi-distant from more than 3 sites will always hold > 3 functions. Hence

Upper bounds and filtering. Not all bisectors are

Voronoï edges. Furthermore, only a specific piece of the bisector curve of two sites is contained in the VD, the one connecting two adjacent Voronoï vertices. For instance, in the classic VD of points, only a linear amount (w.r.t. the number of sites) of line segments out of the totality of (line) bisectors contribute to the VD. This implies that, locally, only a few sites contribute to the VD and therefore the majority of sites should be filtered out. This idea leads to the filtering process described herein.

The filter is based on the following fact: If the control grids of any two f i | D , f j | D do not intersect, then the one that is superior to the other, say the i-th, does not contribute to the lower envelope over D, or equivalently Vor(i) ∩ D = ∅. This is a direct consequence of the variation diminishing property of Bernstein representation [START_REF] Farin | Curves and surfaces for CAGD[END_REF].

Doing the previous test directly for a set of n sites requires O(n 2 ) time to check all pairs. But we can do better: we compute an upper bound U D on the lower envelope over D, and then compare every control grid f i | D against this bound, excluding those that are found to be over it. This leads to a linear time filter, w.r.t. n = #B(D).

Before computing U D , we need bounds on every f i | D . For this we use the minimum (resp. maximum) Bernstein coefficient of f i | D to bound every f i from below (resp. above). These extreme coefficients yield two supporting planes, parallel to xy-plane, that enclose the values f i (D). More sophisticated bounds exist, see [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF][START_REF] Peters | Sleves for planar spline curves[END_REF][START_REF] Reif | Best bounds on the approximation of polynomials and splines by their control structure[END_REF] and references therein. We choose to employ constant time bounds, since they are quite efficient in practice. Indeed, the control grid is known to converge with quadratic speed to the function it represents [START_REF] Reif | Best bounds on the approximation of polynomials and splines by their control structure[END_REF], thus a small number of refinements is needed to separate the control grids of two non-intersecting patches. Now, to compute U D , an upper bound on the lower envelope over D, consider (x, y) ∈ D and let L(x, y) be the value of the lower envelope at (x, y) ∈ D. We have:

L(x, y) = min i {f i (x, y)} ≤ min i {maxcoef(f i | D )}.
(

) 2
This directly proposes to take

U D = min{maxcoef(f i | D )}.
Having U D , we can check which distance fields are bounded over this limit, and filter them out. As soon as the process starts, there is a pre-processing of the graph G. In this step, all the boxes of the graph G that touch the boundary of D 0 are connected in a CCW-directed loop. This loop constrains the reconstruction phase in D 0 . First we discuss how every single box is treated, when encountered in the traversal. Every bisector in the box is approximated inside a domain D by line segments connecting its intersections with ∂D, or other bisectors. For getting the arrangement inside D, we use intersection points and the sign of derivatives (cf. [START_REF] Mantzaflaris | A subdivision approach to planar semi-algebraic sets[END_REF]). Computation of intersections of bisectors is carried out by a Bernstein solver [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF]. Then, signature information that is attached to bisectors is used to identify Voronoï vertices.

If there are no degenerate vertices, and an adequetely small tolerance, the reported diagram is isotopic to the exact one. In practice, subdivision around a degenerate vertex goes down to ε-precision. At this point the vertex is considered to be unique and is attached to all involved bisectors. Therefore, it is possible that two (or more) close vertices are reported as one degenerate vertex. The impact of this change on the global structure of the diagram is small: we may neglect some edges of the dual graph (ie. the generalization of Delaunay's graph, see for instance [START_REF] Emiris | Exact Delaunay graph of smooth convex pseudo-circles: general predicates, and implementation for ellipses[END_REF]), but no false edges are introduced.

The same holds even in the case of a tiny cell of size lower than the working tolerance that may be discarded, due to agglomeration of close vertices. The cell is reported empty, therefore some edges attached to the respective site are missing in the dual graph.

The defining polynomial f if j of a bisector (i, j) is negative over the (open) Voronoï cell of i and positive over the cell of j. The traversal of the spanning boxes uses this sign information: The signs of the bisector on the corners of D give us the correct orientation for tracking the cells corresponding to i and j.

Alg. 2 summarizes the reconstruction process. Starting from an arbitrary box, the boundary of the cell of site i is tracked by traversing the box-span G. The navigation is done based on the sign of the bisector f i -f j . When a vertex is reached, j is updated, so that we follow the next edge of V or(i), and the process continues until the recovery of the whole partition V .

The output of this algorithm is a polygonal approximation of each Voronoï cell. Fig. 3 In this section we sketch how the framework is extended to implicitly defined distance fields.

There are distance fields that cannot be made polynomial by a mere squaring, or by another suitable transformation. A nice polynomial formula even for Euclidean distance fields no longer exists when the sites are not points. In such cases it is natural to represent the function in implicit form, that is, as a polynomial g ∈ R[x, y, z], s.t. a triple satisfying g(x, y, z) = 0, with z > 0 and minimal among the solutions implies that (x, y) is at distance z from the site. Some diagrams that fall in this class are the Apollonius diagram [START_REF] Emiris | The predicates of the apollonius diagram: Algorithmic analysis and implementation[END_REF], the VD of ellipses [START_REF] Emiris | Exact Delaunay graph of smooth convex pseudo-circles: general predicates, and implementation for ellipses[END_REF], or the VD of closed curves given by support function representation [START_REF] Gravesen | Curves and surfaces represented by polynomial support functions[END_REF].

In order to have a representation of g i (x, y, z) over the initial domain D 0 ⊂ R 2 , first we compute an interval I 0 for the third variable. This must contain the z-values of zeros of g. For this we use an upper bound (e.g. Cauchy's bound) for the roots of the (univariate) interval polynomial g i (D, z). Then we feed the subdivision Alg. 1 with g i | D0×I0 , for every site i in the generating set.

The same ingredients, i.e. the "field bound" and "bisector tracking" (Sect. 2) are needed, in order to run the scheme. The subdivision and reconstruction is done the same way, on 3D boxes that enclose implicit space-curve branches. 3.2 Bisector tracking. Bisectors are defined by the (projection of the) intersection of two implicit surfaces, thus they are implicit spatial curves. One way to get the projection is using resultants, but this would be costly to do in every box, and would increase the degree of the polynomial to handle. Isotopic meshing of 3D curves, given in Bernstein form, has been implemented in [START_REF] Liang | Subdivision methods for the topology of 2d and 3d implicit curves[END_REF], and yields a suitable way to derive the projection. Fig. 4 shows a snapshot of the implicit method. Sites are segments on the real axis and distance fields are implicit cones. The box D 1 × I 1 contains 3 bi-variate Bernstein polynomials. Based on U D1 and the lower bound µ I1 (c), the f c | D1 is filtered out.

3.3

Deriving the implicit equation. The offset curve of a site is closely related to the distance function. The z-offset is the set of points at distance z from the site. Therefore, we need to compute an implicit equation containing the z-offset, having z as a parameter.

For the Apollonius diagram, the z-offset of a circle-site centered at p i and with radius w i is the circle of the same center and radius z + w i (see [START_REF] Aichholzer | Divide & conquer algorithms for Voronoi diagrams revisited[END_REF] ). Let n(t) : R → S 1 be a parametrization of the circle S 1 . A supported curve (cf. [START_REF] Gravesen | Curves and surfaces represented by polynomial support functions[END_REF]) is a curve given by a parametrization h(t) over the unit circle, which associates to every point n(t) of S 1 the point on the curve that has normal vector n(t). The offset curve at distance z has support h(t) + z. After elimination of t from (h(t) + z) • n(t) -(x, y) = 0 we arrive to the implicit form of the z-offset, ie. the distance field of the supported curve.

Anton et al. [START_REF] Anton | The offset to an algebraic curve and an application to conics[END_REF] compute, using resultants, the offset curve of radius z of a conic, e.g. an ellipse E. This is a polynomial of degree 4 w.r.t. z and of degree 8 w.r.t. x, y. For any specific point q = (x, y) outside of the ellipse, we have an irreducible univariate polynomial that has exactly one real positive solution, corresponding to dist(q, E). We note that this equation can also be derived if we consider a support function parametrization of the ellipse.

Experimentation

We run experiments on Euclidean, anisotropic and VD under the p metric, for p = 4, 8. For these tests, we set ε = 0.001 and D 0 = [-2, 2] 2 . Also, for the sake of a smooth result, we set a maximum admissible box-size, i.e. all boxes are subdivided until their longest side is at most = 0.05. Table 1 reports execution time (in secs) and number of boxes produced, for random inputs of 50 up to 800 sites in D 0 . The timings are of the same order for all VD types. This can be explained by the fact that the process discards the linear nature of the Euclidean VD. On the other hand, it is applied to the non-linear, anisotropic case, where bisectors are general conics, with the same performance. The reason is that the total time depends mostly on the number of boxes that are produced as output of Alg. 1. This adaptability to general topologies, or mixed diagrams (see Fig. 5), seems to be an advantage of our approach.

The experiments showed a stable method. For example, nearly touching bisectors, as in the cell on the top right of Fig. 1, were correctly separated. Near-degenerate Voronoï vertices were reported as a single vertex of high valency. We note that for computing bounds for univariate root-finding or for the filtering, the rounding mode in the floating point operations is controlled, so that we always enclose the exact values. Nevertheless, exact treatment of all degenerate cases would require manipulation of algebraic numbers, something that we refrain from, for the sake of efficiency.

Fig. 1 .

 1 Fig. 1. Anisotropic diagram. Subdivision box span by Alg. 1 is shown.

Algorithm 1 :

 1 Subdivision phaseInput: A set of distance fields f 1 , f 2 , . . . , fn ∈ R[x, y], a threshold ε > 0 and a rectangle domain D 0 . Output: A graph G of boxes that span the minimization diagram of f 1 , . . . , fn. Compute B(D 0 ) := {f 1 | D 0 , . . . , fn| D 0 }; Initialize stack Q and add B(D 0 ) on top of it; Initialize empty box graph G; while Q is not empty do Pop a box B(D) from Q; Apply filter (Sect. 2.2) on B(D); if |B(D)| ≤ 3 or |D| < ε then Add B(D) := {(f i -f j )| D : i < j ∈ sig(B(D))} to G ; else Split B(D) into B(D 1 ) and B(D 2 ) ; Update adjacency graph G with B(D 1 ), B(D 2 ) ; Push B(D 1 ) and B(D 2 ) into Q ; return G;they will be subdivided until reaching threshold size ε > 0.Whenever we arrive to a box with |B(D)| ≤ 3, or if the box reaches the threshold size ε, the box is mutated : This means that we form differences f i | Df j | D , i < j and we store them in the place of f i | D , i ∈ sig(B(D)). Thus the box is transformed into bisector box. This box will be further treated in the next phase.Three main operations constitute every iteration of the subdivision loop: First, an application of the filter of Sect. 2.2 on the subdivided boxes. Then, an execution of DeCasteljau's algorithm to split the representation of the distance functions. Finally, an update of the adjacency graph, i.e. connection of newly created boxes with their surrounding neighbors. These operations are O(n) for a box with |B(D)| = n. The overall time of Alg. 1 depends on the efficiency of the filter, since it affects the stopping condition |B(D)| ≤ 3.

Fig. 2

 2 illustrates the filtering process in 1D. Three pointsites a, b, c yield the red distance fields, that are squared Euclidean metrics. Domain D 1 holds all three f a , f b , f c , and the dotted line is the level of U D1 , equal to maxcoef(f b | D1 ). Since f c | D1 is over the dotted line, it will be filtered out of the list B(D 1 ). Similarly, maxcoef(f c | D2 ) defines U D2 . Consequently, f b | D2 is excluded and D 2 is dominated by site c, i.e. D 2 ⊂ Vor(c).

2. 3

 3 Cell reconstruction phase. The box-graph computed by Alg. 1 spans the VD, having edges between adjacent boxes. Traversing this graph, by navigating based on signature and sign of the bisector leads to the cells of the VD. In particular, we shall traverse all boundary-contours of the Voronoï cells in counter-clock wise (CCW) order, completing open cells with pieces of the boundary of the initial domain D 0 .

2 Fig. 2 .

 22 Fig. 2. process applied to domains D 1 and D 2 , in computing the VD of tree points a, b, c ∈ R.

1 Fig. 3 .Algorithm 2 :

 132 Fig. 3. Approximating the bisector of a and b. main defines the bisector point (in 1D) of sites a, b. Algorithm 2: Cell reconstruction phase Input: A graph G of boxes that span the minimization diagram of f 1 , . . . , fn. Output: A partition V of D 0 , defined by meshing the minimization diagram of f 1 , . . . , fn. Join border boxes of G in a CCW loop ; Initialize V [i] = ∅, i = 1, . . . , n ; foreach unvisited node (box) B(D) of G do foreach label i ∈ sig(B(D)) do foreach label j ∈ sig(B(D)), j = i do • Traverse G, starting from B(D), and track bisector (i, j) along boxes with label i ; • If some B (D) contains a vertex, set j := j , (i, j ) ∈ B (D) and continue, until B(D) is reached for the second time ; • Add tracked cell component to V [i]. return V;

Fig. 4 .

 4 Fig. 4. Instance of Apollonius diagram in 1D: sites a, b, c are line segments and distance fields are bi-variate (implicit) cones. 3.1 Field bound. For g| D×I = γ ijk B i dx (x) B j dy (y) B k dz (z) we set m k = min i,j {γ ijk } and M k = max i,j {γ ijk }. Let m g (z) :=

Acknowledgment

  This research received funding from [FP7/2007-2013], Marie Curie ITN SAGA, [PITN-GA-2008-214584].

Fig. 5 .

 5 Fig. 5. A mixed VD of anisotropic, and sites with p distance function for p = 2, 4 or 6. Offset curves (in red) reveal the nature of each site.