Luca Bortolussi

Diego Latella

Mieke Massink

Stochastic Process Algebra and Stability Analysis of Collective Systems

Keywords: Fluid flow, process algebra, crowd dynamics, self-organisation

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

A key factor to allow modern cities to reach or maintain a good and sustainable quality of life for their increasingly numerous inhabitants is the development of systems that are relying on a much more decentralised and distributed design that is adapting itself to dynamically changing circumstances [START_REF] Naphade | Smarter cities and their innovation challenges[END_REF]. Examples of such future systems are electricity grids that can cope with many local electricity producers and consumers, and a decentralised organisation of transportation and information. Such large scale collective adaptive systems rely on the continuous feedback between vast numbers of participants of different kinds, and, as is wellknown, can be expected to show complex dynamic and emergent behaviour or perhaps even exploit such behaviour [START_REF] Frei | Advances in complexity engineering[END_REF]. Also the formal analysis of such systems poses many new research challenges.

Process algebras have been specifically designed for the compositional and high-level modelling and analysis of distributed concurrent systems. Recently some of them, in particular PEPA [START_REF] Hillston | A compositional approach to performance modelling[END_REF] and Bio-PEPA [START_REF] Ciocchetta | Bio-PEPA: A framework for the modelling and analysis of biological systems[END_REF], have been provided with a fluid semantics based on ordinary differential equations (ODE) [START_REF] Hillston | Fluid flow approximation of PEPA models[END_REF] providing a scalable approach to the analysis of agent coordination in large collective Work partially supported by the project "FRA-UniTS".

systems. In this paper we further exploit this development to study stability aspects of collective adaptive dynamic systems in a symbolic and numeric way. Stability analysis provides important information about the predictability of dynamic systems and their sensitivity to parameter values. Numeric stability analysis is of particular interest for the analysis of distributed adaptive strategies that are applied in asymmetric situations. Such situations occur naturally in many real-world, natural or designed, systems. An analytical approach to stability analysis in the presence of asymmetry is infeasible in most cases. To the best of our knowledge this is the first time that a method is proposed for stability analysis that starts from a stochastic process algebra specification of agent coordination in a collective dynamic system. The use of a process algebra greatly facilitates the modelling of variants of a system at a high level. Designers can then focus on the coordination strategies instead of having to manipulate the underlying, possibly large set of non-compositional, ODEs in ad-hoc ways.

We illustrate the approach and related tool chain by analysing a new variant of a collective model of spontaneous self-organisation of drinking parties in the squares of cities in Spain, also known as "El Botellón" [START_REF] Rowe | El Botellón: Modeling the movement of crowds in a city[END_REF]. In this variant the parameter of the system is the level of socialisation, i.e. the average number of friends people have, instead of the the probability to find a partner to chat with 1 . The model shows a number of surprising behavioural aspects. However, the main contribution of the paper is to illustrate and explore what stochastic process algebra can offer to provide high-level models of coordination in large scale collective systems, in combination with well-established numeric and symbolic analysis techniques for a systematic stability analysis of such systems. In particular in case of models showing various forms of asymmetry.

Related work can be found in several directions. In [START_REF] Massink | Modelling non-linear crowd dynamics in Bio-PEPA[END_REF] the original model of Rowe and Gomez was analysed with Bio-PEPA. Both empirical and analytical justification was provided for the good correspondence found between stochastic simulation results and the Bio-PEPA based fluid flow approximation. Moreover, a comparison of analytical results with those obtained numerically via Bio-PEPA was provided. For this reason the focus was on symmetric models that can be handled analytically. It does not include a systematic study of the stability aspects of the model. In [START_REF] Bortolussi | Revisiting the limit behaviour of "El Botellon[END_REF] a variant based on the socialisation level has been analysed in an analytical and numerical way. The study addresses a stability analysis of a symmetric model limited to three squares, but it contains no stability analysis of asymmetric models. Recently, the use of Bio-PEPA has also been explored in the field of swarm robotics, where it was used to model a swarm decision making strategy [START_REF] Massink | Analysing robot swarm decision-making with Bio-PEPA[END_REF][START_REF] Massink | Fluid analysis of foraging ants[END_REF]. This collective decision-making strategy has been used as a benchmark for the application of stochastic process algebra, and in particular Bio-PEPA with locations, in this new field of application. It was shown that important aspects of swarm robotics can be addressed such as cooperation and space-time characteristics, but also emergent behaviour.

In the following we first introduce the crowd model and Bio-PEPA followed by a description of the tool pipeline and its application.

A Socialisation Level Based Crowd Dynamics Model

Rowe and Gomez [START_REF] Rowe | El Botellón: Modeling the movement of crowds in a city[END_REF] analyse the movement of crowds between four city squares, connected in a ring by streets, using an agent based approach. The movement of individual people is simulated by agents following a simple set of rules. At every step each agent tries to find a partner to chat with. If this succeeds it stays where it is; else, it moves randomly to an adjacent square. It is assumed that the probability of the latter is (1 -c) pi-1 when this agent is at square i, and p i > 0 is the number of agents currently at square i. The parameter c (representing the chat probability, 0 ≤ c ≤ 1) is the probability that an agent finds a partner to talk to and thus remains in the square. Rowe and Gomez showed that the emergence of crowds in their model is directly linked to a critical threshold value of the chat probability. If the value is below the threshold, the population remains evenly distributed over the squares, while walking around. If the value is above the threshold, the population tends to gather into one or a few squares. However, the probability to meet a friend in a crowded city is in general not the same as the probability to find a friend when it is less crowded. People tend to have a fixed number of friends given a city population, and the larger the number of people walking around, the more of them will turn out not to be one of your friends. This consideration leads to an alternative crowd model, introduced in [START_REF] Bortolussi | Revisiting the limit behaviour of "El Botellon[END_REF], in which the chat probability is defined as c = s/N , where N is the size of the population and s is the level of socialisation of the population, i.e. the average number of friends that people have. Using this alternative definition of c and an n × n routing matrix Q for n squares, i.e. Q i,j is the probability that a person moves to square j, given that she decided to leave square i, the ODE for population level N of this model is:

dx i dt = -x i (1 -s/N ) N xi-1 + j x j (1 -s/N ) N xj -1 Q j,i (1)
where x i denotes the fraction of the population that is in location i. Here we assume that Q is symmetric, i.e. Q i,j = Q j,i , so Q is a stochastic and symmetric matrix. It is further assumed that Q is irreducible (this is not a limitation since otherwise the city can be split into its connected components.) The above ODE is also the fluid flow interpretation of a Bio-PEPA model of this scenario that will be presented in Section 4 with the only difference that the latter is defined directly on the population sizes itself and not on their fractions. The basic symmetric model has also an interesting fluid limit, i.e. an ODE model resulting from letting N → ∞. For N going to ∞ one obtains (see [START_REF] Bortolussi | Revisiting the limit behaviour of "El Botellon[END_REF] for an analytical derivation):

dx i dt = -x i e -sxi + j x j e -sxj Q j,i (2)
This is a non linear ODE, that has all its solutions in the unit simplex

∆ n = {x ∈ R n : x i ≥ 0 and i x i = 1} if the initial condition is in ∆ n .
3 Bio-PEPA and Fluid Flow Analysis

Before presenting a process algebra model of the crowd scenario we briefly recall Bio-PEPA [START_REF] Ciocchetta | Bio-PEPA: A framework for the modelling and analysis of biological systems[END_REF], a language that has originally been developed for the modelling and analysis of biochemical systems. The main components of a Bio-PEPA system are the "species" components, describing the behaviour of individual entities of a species, and the model component, describing the interactions between the various species. The initial amounts of each type of entity or species are given in the model component. The syntax of the Bio-PEPA components is defined as:

S ::= (α, κ) op S | S + S | C with op = ↓ | ↑ | ⊕ | | and P ::= P £ ¡ L P | S(x)
where S is a species component and P is a model component. In the prefix term (α, κ) op S, κ specifies the multiples of an entity of species S involved in an occurring action α2 . The prefix combinator "op" represents the role of S in the action, or conversely the impact that the action has on the species. Specifically, ↓ indicates a reduction of the population involved in the action, ↑ indicates an increase as a result of the action. The operators ⊕, and play a role in an action without leading to increments or decrements in the involved populations and have a defined meaning in the biochemical context. Except from , that will play a role in the shorthand notation introduced below, we will not need these operators in this paper. The operator "+" expresses the choice between possible actions, and the constant C is defined by an equation C=S. The process P £ ¡ L Q denotes synchronisation between components P and Q, the set L determines those actions on which the components P and Q are forced to synchronise, with £ ¡ * denoting a synchronisation on all common actions. In S(x), the parameter x ∈ IR represents the initial amount of the species. The frequency with which an action occurs is defined by its (functional) rate. This rate is the parameter of a negative exponential distribution. Its value may be a function of the size of the populations involved in the interaction.

Bio-PEPA comes with a notion of discrete locations that may contain species. A Bio-PEPA system with locations consists of a set of species components, also called sequential processes, a model component, and a context (locations, functional rates, parameters, etc.). The prefix term (α, κ) op S@l is used to specify that the action is performed by S in location l. The notation α[I → J] S is a shorthand for the pair of interactions (α, 1)↓S@I and (α, 1)↑S@J that synchronise on action α. This shorthand is very convenient when modelling agents migrating from one location to another as we will see in the next section. Bio-PEPA is given an operational semantics, based on Continuous Time Markov Chains (CTMCs), and a fluid semantics, based on ordinary differential equations (ODE) [START_REF] Ciocchetta | Bio-PEPA: A framework for the modelling and analysis of biological systems[END_REF]. The Bio-PEPA language is supported by a suite of software tools which automatically process Bio-PEPA models and generate internal representations suitable for different types of analysis [START_REF] Ciocchetta | Bio-PEPA: A framework for the modelling and analysis of biological systems[END_REF][START_REF] Ciocchetta | The Bio-PEPA Tool Suite[END_REF]. These tools include mappings from Bio-PEPA to differential equations (supporting a fluid flow approximation), stochastic simulation models [START_REF] Gillepie | Exact stochastic simulation of coupled chemical reactions[END_REF], and PRISM models [START_REF] Kwiatkowska | PRISM: Probabilistic model checking for performance and reliability analysis[END_REF].

Bio-PEPA Crowd Model

In this section we define a Bio-PEPA specification of the crowds scenario presented in Section 2. Let us consider a small ring topology with 4 city squares in a 2 × 2 grid, denoting them by 00, 01, 10 and 11, allowing bi-directional movement between squares. In Bio-PEPA the city squares are modelled as locations called sq00, sq01, sq10 and sq11. Parameter c defines the chat-probability and parameter d the degree or number of streets connected to a square. In the symmetric topology d = 2 for each square. The chat-probability is defined as the fraction of the socialisation factor s w.r.t. the total population N , i.e. c = s/N . The actions modelling agents moving from square X to square Y are denoted by fXtY . The associated functional rate for f00t01 with P @sq00 denoting the population in sq00 at time t is defined as: f00t01 = (P @sq00 * (1 -c) (P @sq00-1) )/d; the other rates are defined similarly. The behaviour of a typical agent moving between squares is modelled by sequential component P . For example, f00t01 [sq00 → sq01] P models that an agent present in sq00 moves to sq01 according to the functional rate defined for the action f00t01 .

P = f00t01 [sq00 → sq01] P + f01t00 [sq01 → sq00] P + f00t10 [sq00 → sq10] P + f10t00 [sq10 → sq00] P + f01t11 [sq01 → sq11] P + f11t01 [sq11 → sq01] P + f01t11 [sq10 → sq11] P + f11t10 [sq11 → sq10] P ;
Finally, the model component defines the initial conditions of the system, i.e. in which squares the agents are located initially, and the relative synchronisation pattern. If, initially, there are 1000 agents in sq00 this is expressed by P @sq00[1000]. The fact that moving agents need to synchronise follows from the definition of the shorthand operator →.

((P @sq00[1000] £ ¡ * P @sq01[0]) £ ¡ * (P @sq10[0])) £ ¡ * (P @sq11[0])
The total number of agents P @sq00+P @sq10+P @sq01+P @sq11 is invariant and amounts to 1000 in this specific case. The occupancy measure, i.e. the fraction of the population in sq00 can be defined as P sq00 = P @sq00/N and similarly for the other squares. The fluid semantics of Bio-PEPA leads to an ODE that is very similar to Eq. ( 1) except that it is defined on the actual population sizes rather than their fractions:

dP @sq i dt = -P @sq i (1 -s/N ) P @sqi-1 + j P @sq j (1 -s/N ) P @sqj -1 Q j,i (3)
Using the Bio-PEPA plugin tool suite a first insight in the behaviour of the above model for different values of the socialisation factor s can be obtained using e.g. stochastic simulation [START_REF] Gillepie | Exact stochastic simulation of coupled chemical reactions[END_REF] or one of the built-in ODE solvers. For example, for s = 5 an interesting so called 'metastable' behaviour can be observed in ODE trajectories such as the one shown in Fig. 2(a), where the fractions of the population present in squares sq00 through sq11, are denoted by x i , for i ∈ {1, .., 4}. Until time 150 just over 40% of the population is in each of the nonadjacent squares x 1 and x 3 , and slightly less than 10% in each of the remaining squares. Then suddenly this situation changes and square x 1 gets almost all of the population. This is just one example of the typical kind of behaviours that may occur in non-linear systems such as these.

To get a more complete overview of potential emergent behaviour of a collective dynamic system it may be useful to construct a bifurcation diagram of the system. This is a diagram that shows for each value of a selected parameter of the model its stationary points for that value. For each stationary point it also shows whether it is stable or unstable, i.e. whether a system would remain in a state forever once it is reached, or whether it could still move on from there reaching other states. The selected parameter of interest in our case is the socialisation value s. Fig. 1 shows the bifurcation diagram for square x 1 of the crowds model with four squares and for s ranging from 2.75 to 6. For example, for s = 3 we see that the model has one stable stationary point with value 0.25. This means that for s = 3 the system ends up in a stable state in which 25% of the population is expected to be in square x 1 .

In fact, the model has a stationary point at 0.25 for all values of s considered, i.e. the vector of four squares x sym = 1 4 , 1 4 , 1 4 , 1 4 is always a stationary point, but it is not stable for all values of s; its stability changes at s = 4. The stability of this equilibrium has been analysed in an analytical way in [START_REF] Bortolussi | Revisiting the limit behaviour of "El Botellon[END_REF] for a model with an arbitrary number n of squares, but with symmetric routing matrix Q. There it was proven that the stability status of this stationary point changes from stable into unstable when s becomes equal to n. It is stable for models in which s < n and unstable for those in which s > n. However, in [START_REF] Bortolussi | Revisiting the limit behaviour of "El Botellon[END_REF] no analytic results have been given for the other stationary points, due to the difficulty in identifying them in general. Stability of all fixed points has been discussed in [START_REF] Bortolussi | Revisiting the limit behaviour of "El Botellon[END_REF] only for a symmetric, fully connected, model with three squares. Such difficulties are also common for models that show irregular structure or other forms of asymmetry. When investigating the behaviour of the model in more detail using stochastic simulation (not shown) it turns out that in the model with 4 squares and s < 3.25 the population is migrating between squares in such a way that it is evenly distributed over the squares. This means that at any point in time one would see approximately 25% of the population in each of the squares. However, this situation changes for s > 3.25. Stochastic simulation shows that in that case most of the population eventually gathers at random in one of the squares.

In the following we propose an analysis pipeline, starting from a Bio-PEPA specification, that can be used to compute the stationary points numerically instead of analytically for various values of s, such that a bifurcation diagram of the kind shown in Fig. 1 can be generated. As we can see in the figure, the change of stability of x sym is correctly identified by the numerical method, which also predicts other stationary points, corresponding to three different configurations:

1. Most agents in a single square. This is the top branch in the diagram. These are stable stationary points, but emerge for s greater than 3.5, approximately.

For instance, for s = 5, in this case 98% of agents are in a single square. 2. Agents evenly split between two opposite squares. This is an unstable stationary point, emerging for s > 4, which is in fact a saddle node3 . For s = 5, this configuration corresponds to 43% of agents in two opposite squares (for instance, sq00 and sq11), and 7% in the remaining ones. Notice that this information can be deduced by observing the numeric values of fixed points. 3. Agents evenly split between three adjacent squares. This is again an unstable equilibrium of the system. for s = 5, we have 29% of agents in each of the three adjacent squares and 0.13% in the remaining one. Interestingly, configuration (2) above, despite being a saddle node, hence unstable, has a quite strong stable manifold. The effect is that solutions starting nearby the stable manifold (e.g. close to the plane x 1 = x 3 ) remain for some time close to it, before escaping to one of the stable fixed points. This gives rise to the typical metastable behaviour that is shown in the ODE trajectory in Fig. 2(a). As could be expected, the higher the socialisation coefficient s, the more intense is the attractive behaviour of the stable manifold, hence the longer solutions of the ODE of the model, that start from initial values in the neighbourhood of such a point, remain attracted to it. This can be seen in Figure 2(b), where the time spent nearby the stable manifold is shown as a function of s. In Section 6 we will apply the analysis pipeline proposed in the next section to several variants of the crowd model. Each of these variants is characterised by a different form of asymmetry e.g. caused by squares with different attractivity, or by asymmetry in topology.

Numeric/Symbolic Stability Analysis

Before discussing the tool pipeline, we briefly sketch the methods we are using in Octave/Matlab to investigate the stability behaviour of fixed points and to generate bifurcation diagrams for different asymmetric variants of the crowd model that will be presented in next section. The idea, as discussed previously, is to combine the features of stochastic process algebras as modelling languages, Bio-PEPA specifically, with the numerical algorithms and analysis routines of a platform like Matlab [18] or Octave [1]. Assume that we have obtained an ODE function in the input format of Matlab or Octave, i.e. an m-file. Such an ODE function models the changes in the behaviour of the system over time. The idea is then to use numerical routines to look for the set of zeros of this function, i.e. points for which the system is stationary, for a fixed set of parameters (e.g. in our case for the socialisation coefficient s). More precisely, we can use the Octave/Matlab method fsolve which looks for zeros of non-linear equations and which incorporates different solvers (in our case we used the Levenberg-Marquardt algorithm [START_REF] Nocedal | Numerical Optimization[END_REF]). As fsolve finds one zero, depending on the initial values, we run the algorithm many times to look for more zeros, starting from random initial coordinates (and from some predefined fixed points, like the unit vectors). In order to avoid considering all possible symmetric solutions of the problem, we post-process the set of zeros found in this way, to keep only one zero among those equivalent under symmetry. For instance, in a fully symmetric model with ring topology as the one in Section 4, there is a zero which puts almost all the population in one square, for each square. The post-processing then removes all but one of those zeros. Finally, for each such stationary point, we investigate its stability by computing numerically the Jacobian matrix and its eigenvalues, according to standard methods in dynamic systems theory [START_REF] Strogatz | Non-linear dynamics and chaos: with applications to physics[END_REF][START_REF] Khalil | Nonlinear systems[END_REF]. After discarding one zero eigenvalue (which is always present due to the conservation of the total number of agents), we check whether the real part of the remaining ones is positive or negative. If all eigenvalues have a negative real part (and are sufficiently far away from zero to account for numerical errors), we declare the point stable, otherwise we mark it unstable. If there is an eigenvalue too close to zero, we mark its status as unknown. In order to generate bifurcation diagrams, we perform this operation for a range of values of the parameter of interest, e.g. in this case the socialisation coefficient s.

Tool Pipeline. In order to link the above described procedure to Bio-PEPA, we export the model in the SBML format [4]. Such an export is already available via the Bio-PEPA plugin tool suite [START_REF] Ciocchetta | The Bio-PEPA Tool Suite[END_REF]. Both Matlab and Octave have a toolbox importing from SBML files [START_REF] Keating | SBMLToolbox: an SBML toolbox for MATLAB users[END_REF], generating an m-file computing the vector field (ODE) corresponding to the fluid semantics of Bio-PEPA [START_REF] Ciocchetta | Bio-PEPA: A framework for the modelling and analysis of biological systems[END_REF]. Once an m-file has been obtained, it can be used within the routines explained above by creating a function handle. However, following the above procedure, we have not yet obtained the limit model discussed in Section 2, as from Bio-PEPA we export the N -dependent model (see Sect. 4). Even if the two sets of ODE will produce very similar solutions when N is large enough, working with the limit model would be preferable, as this seems to reduce the numerical errors caused by the exponentiation. A limit model can be obtained by exploiting a computer algebra system such as the symbolic toolbox in Matlab [START_REF]Matlab symbolic toolbox[END_REF] or the open source software Maxima [START_REF]SBML: Systems Biology Markup Language[END_REF]. This requires that the Bio-PEPA specification is exported to an m-file containing a symbolic definition of the ODE equations in Matlab or an equivalent format suitable for use with Maxima. This is not very difficult to automatise. Once this operation is performed, the symbolic calculus routines can be exploited to compute the limit of the vector field and to compute the Jacobian matrix symbolically, increasing the precision of the method and speeding up the numerical analysis phase. The resulting functions can be either exported from Maxima to Matlab/Octave by a suitable script that generates an appropriate m-file, or by converting a symbolic function into a numerical one in Matlab.

Results

In this section we present some results for variants of the basic crowd model enriching it either by adding new behaviours or by breaking the symmetry between squares. Due to space limitations, we will mainly report on results for a model with four squares connected in a ring topology. However, we have also obtained interesting results for a larger a model with 9 squares disposed in a 3 × 3 grid-like topology. The modifications in the basic model are essentially of four types:

-Breaking the symmetry between squares (in a symmetric topology like the ring one) by assigning to each square an uneven attractiveness coefficient. Each agent then chooses the next square to go according to their relative attractiveness.

-Breaking the symmetry in the routing, by assuming that certain connections between squares can be crossed only in a given direction, i.e. by introducing one way pedestrian streets (which can be enforced for instance by the presence of police). -Breaking the symmetry in the topology, by having different degrees of connection between different squares, e.g. the 9 squares model (not shown). -Breaking the symmetry in the behaviour of agents, assuming that an agent first decides if she wants either to leave the square or look for somebody to chat to. In this second case, we assume that she leaves the square with the same probability as in the standard model.

In particular, in each of these cases, we will discuss the stability of stationary points as a function of the socialisation parameter s, or of other parameters, like attractiveness or the probability to leave. As discussed in Section 5, the results always consider the limit fluid model (see Eq.( 2)), for a population going to infinity. Due to convergence of vector fields, the results for the limit model can be expected to be the same as those of the ODE with explicit dependence on N , if N is reasonably large. In particular, we have seen basically no differences between results for the limit model and those for a model with a total population of N = 1000.

square model with attractive squares

The first kind of asymmetric model we consider is one in which an attractiveness coefficient is assigned to each square. Higher attractiveness is modelled by a higher value of the coefficient. When agents decide to leave a square, their decision to which of the adjacent squares to go is now proportional to the relative attractiveness of the adjacent squares. In particular, we consider a situation in which only one square (by convention, square sq00), has a higher attractiveness than the others. Hence, the attractiveness coefficient of square sq00 is equal to a ≥ 1, while that of the other squares is set to 1. In Bio-PEPA this is modeled by replacing the transition rates in the symmetric model of Sect. 4 by the following ones, assuming attractiveness coefficients a00 for square sq00, a01 for square sq01 and so on, e.g.: f00t01 = (P @sq00 * (1 -c) (P @sq00-1) ) * (a01/(a01 + a10));

In Fig. 3 we show a bifurcation diagram as a function of the socialisation factor s, for the fraction of people in the first square, x 1 , for a model with a = 3. As we can see, the solution in which the more attractive sq00 gets the larger amount of people is always stable. However, for s around 3.5, a new set of stationary points appears, with a stable and an unstable branch. The stable branch corresponds to situations in which most of the people stay in one of the two squares adjacent to the attractive one. This stable equilibrium is quite surprising, as one would expect that people always move towards the attractive square. Note that this happens for s < 4, i.e. the predicted threshold for s in the symmetric model.

The other stable equilibrium that emerges around s = 4.5 corresponds to the situation in which all agents are in the square opposite to the attractive one.

Hence, even in the presence of asymmetric attractiveness, we may still obtain a pattern of emergent behaviour in which all people are gathered at a random single square. So, contrary to intuition, for sufficiently large values of s, people may concentrate also in non-attractive squares. Furthermore, as expected, we loose the symmetric equilibrium in which agents are uniformly distributed over the squares and find no replacement for it.

4 square model with one-way streets In this model we assume that streets between squares are unidirectional rather than bi-directional. In particular, we assume that the ring can be traversed only clockwise, with no square having a higher attractiveness than others. This model is still ergodic, and moreover the four squares still behave symmetrically. However, the routing matrix is no longer symmetric. In the Bio-PEPA model this can be easily obtained by removing the related directions of movement and their corresponding rate definitions.

The pattern of the bifurcation diagram we obtain is very similar to the one for the symmetric model in Fig. 1. Indeed, the only difference is that now some eigenvalues of the Jacobian in steady states have imaginary values, suggesting that we may in fact have stable and unstable foci rather than single nodes [START_REF] Strogatz | Non-linear dynamics and chaos: with applications to physics[END_REF], i.e. convergence to the fixed point happens by damped oscillations. In any case, this effect is very weak, and cannot be observed at the resolution scale at which we plot trajectories.

4 square model with independent leaving probability In this model we change the behaviour of single agents. In particular, we assume that each agent first chooses if it wants to leave the square (with a "boredom" probability p) or look for another one to chat with. In the latter case, it behaves like in the original model. In Bio-PEPA this is modeled by replacing the transition rates in the model with attractiveness by: f00t01 = (1 -p) * (P @sq00 * (1 -c) (P @sq00-1) ) * (a01/(a01 + a10)) + p * (P @sq00); Note that in case all squares have an attractiveness coefficient equal to 1 we can study the effect of the boredom probability separately from the effect of attractiveness of squares. This is what is assumed in the following.

The model with independent leaving probability has a second parameter in addition to s, namely the boredom probability p. Intuitively, this probability should influence the overall behaviour of the model quite radically. If it is large enough (maybe in case of a meeting of eremites), then we do not expect people to gather in a single square, but rather to see them uniformly distributed over the four squares while moving between them. This indeed happens for p = 0.15 and for any value of s considered. On the other hand, for small values of p, the behaviour manifested by the system becomes much more complicated, as shown in Fig. 4. In this case, for p = 0.05, we can observe a frequent change in the stability status of fixed points, due to bifurcation events in which stationary points split, change stability status, and move towards other fixed points, generating a cascade of bifurcations.

The most interesting feature is the stability of the symmetric equilibrium. Initially it is stable, but then, as s increases, it undergoes a bifurcation event and becomes unstable (for s around 5 here), as in the original model. At the same time, a new stable branch emerges, corresponding to the emergent behaviour in which most people are concentrated in a single square. However, for s around 12, this stable behaviour undergoes a new bifurcation event and becomes unstable. Furthermore, as s increases further, this unstable equilibrium approaches the symmetric one and hits it around s = 16. When this happens, the unstable symmetric equilibrium becomes stable again, and as s increases even further, it becomes the only (stable) fixed point. This is a counter-intuitive behaviour of the model: as the socialisation factor increases, instead of obtaining a higher probability of people concentrating in a single square, exactly the opposite effect emerges.

The assumption that agents can leave a square with a small fixed probability, no matter whether friends are present or not, is quite reasonable. For example one may receive a text message from a friend in another square or have other things to do. The behaviour of this model shows that, when the socialisation coefficient is large enough (in reality, for a population of a few thousand people, we can expect a value of s well above 30), the chat probability mechanism is not able to fully explain the emergent behaviour of 'El Botellon'. Other mechanisms have to be taken into account and also their potential interference. Among these are most likely some asymmetry breaking phenomena such as the different degrees of attractiveness of the squares.

To obtain deeper insight in why we observe the effect shown in Fig. 4, we can compare the total exit rate of an agent from a square for the symmetric model and the model with boredom probabilities. In particular, we compare these rates as a function of s, for a fixed fraction α of the total population in the square. In the symmetric case, we obtain e -αs (see Eq.( 2)), which is a decreasing function of s, exponentially approaching zero for large values of s. Hence, for a sufficiently large value of s, the rate of leaving a square becomes very small: if there are many agents in the square (α close to 1), then we can expect nobody will move, as the exit rate from the crowded square will be much smaller than the exit rate of the other squares. On the other hand, for the model with boredom probability, the probability at which a single agent leaves a square is p + (1 -p)e -αs , which for large s converges to p. This means that for large s, the effect of the chat probability to remain in a square is negligible compared to the effect of the boredom probability with which agents leave a square. This, in turn, implies that the ODE we obtain for large s is essentially a set of linear ODE with a small non-linear perturbation term, hence they converge to the unique equilibrium of the linear system, which is the symmetric one in which the population is uniformly distributed over the squares.

In order to break this effect, we may want to find out whether the introduction of asymmetric attractiveness can counterbalance the disruptive effect of boredom probability on the emergent behaviour of the model. Fig. 6 (left) shows the bifurcation diagram for the model which sq00 is three times as attractive as the other squares, a = 3, and boredom probability with p = 0.05. We observe the same pattern as in the symmetric model: for large enough values of s, the system converges to the perturbed symmetric equilibrium (due to asymmetric attractiveness). An additional variation can be to consider a boredom probability which is inversely proportional to the attractiveness of squares, for instance equal to p/a 2 , in order to take into account the effect that it is more unlikely that people just leave a square where interesting events are going on. In this case, in the presence of the compensating effect of boredom probability for large s, we can enlarge the range of s for which we observe an emergent party (see Fig. 6 (right)). So, a combination of effects of the attractiveness on the choice of the next square and on the boredom probability can still be used to explain the emergent behaviour of 'El Botellon'.

Discussion and Further Work

The engineering of large scale collective dynamic systems is a relatively new domain of software engineering and requires effective and scalable formal analysis methods. We illustrated how the exploitation of the combination of process algebras, designed specifically for concurrent systems, and techniques typical of dynamic systems analysis, such as stability analysis, could provide valuable tools for such engineering approaches. In this paper we proposed a tool pipeline that, starting from a Bio-PEPA specification of non-linear asymmetric collective coordination, can produce bifurcation diagrams. These can be used to analyse the effect of potentially interfering coordination mechanisms on the stability properties of a system. We have illustrated the combined numeric and symbolic approach on a number of variants of a model of crowd dynamics that represents various kinds of asymmetry. The automatisation of the approach is feasible and part of future work. An issue that needs further investigation is the scalability of the method. Fluid ODEs are independent on the population size, but not on the number of states of the agents. Finding numerically all the zeros of a vector field can be challenging in large dimensions. One approach would be to parallelise the code and use more efficient zero finding numerical routines. An alternative and more promising direction is to exploit the formal nature of process algebras to reduce the agent's state space by using behavioural equivalences or abstract interpretation. Furthermore we plan to integrate the approach with the other already available analysis tools for Bio-PEPA, such as fluid flow analysis and stochastic simulation, via the Bio-PEPA tool suite.

Fig. 1 .

 1 Fig. 1. Bifurcation diagram for the symmetric 4 square model.

Fig. 2 .

 2 Fig. 2. (2(a)) A trajectory for s = 5 showing metastable behaviour. The initial conditions are x1 = 0.4 + δ, x2 = 0.1, x3 = 0.4 -δ, x4 = 0.1, with δ small. (2(b)) Time spent close to the metastable equilibrium, as a function of s. Initial conditions are as above, with δ = 10 -7 .

Fig. 3 .

 3 Fig. 3. Bifurcation diagram for the 4 square model with asymmetric attractiveness, as a function of s, holding a fixed to 3 .

Fig. 4 .

 4 Fig. 4. Bifurcation diagram for the 4 square model with boredom probability p = 0.05.

Fig. 5 .

 5 Fig. 5. Bifurcation diagram for the attractive square in 4 square model with boredom probability. for p = 0.05 and attractiveness equal to 3 (left), and attractiveness dependent p (right).

Studied in the original work by Row and Gomez[START_REF] Rowe | El Botellón: Modeling the movement of crowds in a city[END_REF] and in[START_REF] Massink | Modelling non-linear crowd dynamics in Bio-PEPA[END_REF]

The default value of κ is 1 in which case we simply write α instead of (α, κ).

Unstable stationary points come in different kinds. One of these kinds is called a 'saddle node', informally it is attracting from two opposite sides and repellent from the two other sides of the point.

Acknowledgments. This research has been partially funded by the EU-IP project ASCENS (nr. 257414), the EU-FET project QUANTICOL (nr. 600708) and the Italian MIUR-PRIN project CINA.


Tiago Cogumbreiro

Francisco Martins

Vasco Thudichum Vasconcelos

Coordinating phased activities while maintaining progress

In order to develop reliable applications for parallel machines, programming languages and systems need to provide for flexible parallel programming coordination techniques. Barriers, clocks and phasers constitute promising synchronisation mechanisms, but they exhibit intricate semantics and allow writing programs that can easily deadlock. We present an operational semantics and a type system for a fork/join programming model equipped with a flexible variant of phasers. Our proposal allows for a precise control over the maximum number of synchronisation steps each task can be ahead of others. A type system ensures that programs do not deadlock, even when they use multiple phasers.

Introduction

The key to develop scalable parallel applications lies in using coordination mechanisms at the "right" level of abstraction [START_REF] Cole | Photoshop scalability: Keeping it simple[END_REF]. Rather than re-inventing the wheel with ad hoc solutions [START_REF] Xiong | Ad hoc synchronization considered harmful[END_REF], programmers should resort to off-the-shelf coordination mechanisms present in programming languages and systems. Barriers, in their multiple forms [START_REF] Aditya | Semantics of barriers in a non-strict, implicitlyparallel language[END_REF][START_REF] Barnes | Barrier Synchronisation for occampi[END_REF][START_REF] Cavé | Habanero-Java: the new adventures of old X10[END_REF][START_REF] Charles | X10: an object-oriented approach to non-uniform cluster computing[END_REF][START_REF] Dagum | OpenMP: An Industry-Standard API for Shared-Memory Programming[END_REF][START_REF] Leijen | The design of a task parallel library[END_REF][START_REF] Reinders | Intel Threading Building Blocks: Outfitting C++ for Multi-core Processor Parallelism[END_REF] constitute one such coordination mechanism. A barrier allows multiple tasks to synchronise at a single point, in such a way that: a) before synchronisation no task has crossed the barrier, and b) after synchronisation all tasks have crossed the barrier.

Programs that use a single barrier to coordinate all their tasks do not deadlock. If using a single barrier may reveal itself quite limited in practice, groups of tasks that use multiple barriers may easily deadlock. To address this issue, the X10 programming language [START_REF] Charles | X10: an object-oriented approach to non-uniform cluster computing[END_REF] proposes clocks, a deadlock-free coordination mechanism. Clocks later inspired primitives in other languages, such as Java [START_REF]Java Specification Request JSR-166[END_REF] (as of version 7), Habanero Java [START_REF] Shirako | Phasers: a unified deadlock-free construct for collective and point-to-point synchronization[END_REF] (HJ), and an extension to OpenMP [START_REF] Shirako | Unifying barrier and point-to-point synchronization in OpenMP with phasers[END_REF].

For some applications, the semantics of traditional barriers is overly inflexible. With this mind, Shirako et al. introduced phasers [START_REF] Shirako | Phasers: a unified deadlock-free construct for collective and point-to-point synchronization[END_REF], a primitive that allows some tasks to cross the barrier before synchronisation, thus relaxing condition a) above. Phasers allow for asymmetric parallel programs, including multiple producer/single consumer applications where, at each iteration, the consumer waits for the producers to cooperate in assembling an item. In a different direction, Albrecht et al. proposed a partial barrier construct [START_REF] Albrecht | Loose synchronization for large-scale networked systems[END_REF] that only requires some tasks to arrive at the barrier, thus relaxing condition b). This form of barriers can be used in applications that synchronise on a subset of all tasks, allowing, e.g., computation to progress quickly even when in presence of slow tasks.

Reasoning about such enriched constructs is usually far from trivial: the semantics is intricate and most languages lack a precise specification (including Java barriers and HJ phasers). In this work we propose a calculus for a fork/join programming model that unifies the various forms of barriers, including X10 clocks [START_REF] Charles | X10: an object-oriented approach to non-uniform cluster computing[END_REF], HJ phasers (both bounded [START_REF] Shirako | Phaser beams: Integrating stream parallelism with task parallelism[END_REF] and unbounded [START_REF] Shirako | Phasers: a unified deadlock-free construct for collective and point-to-point synchronization[END_REF]), and Java barriers [START_REF]Java Specification Request JSR-166[END_REF].

Our proposal not only subsumes those of X10, HJ and Java, but further increases the flexibility of the coordination mechanism. We allow tasks to be ahead of others up to a bounded number of synchronisation phases, and yet guarantee that well-typed programs are deadlock free. In contrast, HJ allows tasks to be ahead of others by an arbitrary number of phases, which is of limited interest, since fast tasks may eventually exhaust computational resources, such as buffer space.

To summarise, our contributions are:

a flexible barrier construct that allows for a precise control over the maximum number of phases each task can be ahead of others; an operational semantics for a fork/join programming model that captures barrier-like coordination patterns found in X10, HJ, and Java; a type system that ensures progress, hence the absence of deadlocks, in addition to the usual type preservation.

The paper is organised as follows. The next section addresses related work. Section 3 presents the syntax and the operational semantics of our language. Thereafter, we introduce the type system and the main results of our work. Section 5 concludes the paper while putting forward lines of future research.

Related work via an example

Figure 1a sketches a parallel breadth-first search algorithm to find an exit in a labyrinth. The algorithm uses two groups of tasks: one traverses the labyrinth (modelled by a graph) and another inspects the visited nodes for an exit. The sketch was originally implemented in OpenMP by Süß and Leopold [START_REF] Süß | Implementing irregular parallel algorithms with OpenMP[END_REF]. The algorithm proceeds iteratively, handling at each step (or phase) all nodes at the same depth. If a node is an exit, the algorithm terminates; otherwise it computes the node's neighbours and places them in a buffer to be processed in a later phase. The tasks synchronise at the end of each phase.

Figure 1b depicts two execution states for phases i and i + 1. At phase i, tasks t 1 , . . . , t N read nodes from buf 1 (that stores the nodes of depth level i), compute their descendants, and then in buf 2 (via function traverseNodes()). Tasks c 1 , . . . , c M read nodes from buf 1 and look for an exit node (function checkNodesForAnExit()). Both groups of tasks synchronise (at phaser c) and advance to phase i + 1. Notice that at phase i + 1 the tasks c 1 , . . . , c M continue to process buf 1, while the traversal group t 1 , • • • , t N is handling nodes at depth level i + 1 (buf 2). This situation is possible due to the bound K assigned to phaser c upon launching the traversal tasks (line 5). Bound K denotes the number of phases the traversal tasks may be ahead of the checker tasks. The number of buffers is equal to K + 1.

The synchronisation of tasks is challenging and involves phasers c and t. Phaser c is used to synchronise traversal tasks with checker tasks, while the traversal tasks further synchronise among themselves at phaser t. The bound c:K in line 5 reads as "phaser b has bound K in the spawned task." On the other hand, bound t:0 (also in line 5) means that the traversal activities must all be in the same phase. This synchronisation scheme ensures that the traversal tasks must process all nodes at the same depth and only after that advance to the next level. The traversal tasks set the pace for the checker tasks with phaser c. In their turn, the checker tasks (lines 12-16) use bound c:0, thus enforcing that they simultaneously process the same depth level and do not overtake the tasks in the traversal group.

Barriers, clocks, and phasers are insufficient for this sort of coordination. Barriers and clocks are too inflexible, since no task can cross the barrier before the others arrive. HJ phasers are too loose, since when not used as barriers, tasks run unconstrained and may overflow buffers. Phasers beams [START_REF] Shirako | Phaser beams: Integrating stream parallelism with task parallelism[END_REF] are a step towards this sort of control, but its semantics is only informally described and they do not guarantee deadlock-freedom. Our proposal permits different tasks to specify the maximum number of phases they can be ahead of the slowest.

In contrast, phaser beams define such a number on a per-phaser basis. To unify regular phasers and phaser beams, we also supply an operation that skips waiting on a phase.

Saraswat and Jagadeesan presented a calculus for the X10 language that includes barriers, fork/join, conditional atomic blocks, and hierarchical shared memory [START_REF] Saraswat | Concurrent clustered programming[END_REF]. The authors define a small-step operational semantics and claim that X10 programs without conditional atomic blocks do not deadlock, yet no deadlock-freedom theorem is formally proved. Lee and Palsberg present a calculus, called FX10, with two constructs from X10: fork/join and atomic blocks [START_REF] Lee | Featherweight X10: a core calculus for async-finish parallelism[END_REF]. FX10 is suited for inter-procedural analysis through type inference and includes a formal proof of a fragment of the deadlock theorem stated by Saraswat and Jagadeesan. The type system used to identify may-happen-parallelism is further explored in [START_REF] Agarwal | May-happen-in-parallel analysis of X10 programs[END_REF]. Other formal studies on fork/join semantics include [START_REF] Aditya | Semantics of barriers in a non-strict, implicitlyparallel language[END_REF][START_REF] Arvind | λs: an implicitly parallel λcalculus with letrec, synchronization and side-effects[END_REF]. Our previous work defines an operational semantics and a type system for a calculus with a fork/join programming model and clocks [START_REF] Martins | Types for X10 Clocks[END_REF], but does not include a deadlock-freedom theorem.

X10 clocks and HJ phasers are language-based approaches, whereas the Java barrier (also called phaser) is a library-based approach. Features that appear simultaneous in our calculus, in X10, and in HJ are: controlled barrier registration to avoid non-determinism, advancement on multiple barriers, and enforced barrier deregistration before activity termination to avoid barrier-related deadlocks. Features that appear in our calculus and in HJ alone include phaser visibility restricted to finish scopes to avoid deadlocks between phaser and finish; enforced deregistration before terminating a finish to avoid deadlocks between barriers. HJ and X10 automatically deregister from barriers when activities terminate. Instead we require explicit operations on phasers, in order to obtain a clearer operational semantics. Our type system provides enough information to guide a compiler into automatically inserting such operations, if desired. Finally, the ability to specify the maximum number of phases a task may be ahead of the slowest is unique to our language.

Syntax and operational semantics

Our language, inspired by X10 and HJ, uses activities to organise independent computations and features two coordination mechanisms to control concurrency: phaser and finish. For the sake of simplicity, the language focuses on task coordination, providing very little in the way of describing complex computations.

The syntax of the language is defined in Figure 2. It relies on a base set of variables, ranged over by x, and a set of natural numbers, ranged over by m and n. Bounds B map phaser names to the phaser bounds. We use the standard abbreviation t; e to denote the expression let x = t in e when variable x does not occur in expression e. A term t is transformed into an expression via let x = t in x. We describe the language constructs along with the presentation of the operational semantics. Figure 3 introduces the syntax of a machine state. The run-time system relies on two additional disjoint sets, H and L. Set H contains phaser names, ranged over by h. Activity names, l, are taken from set L. A state S of a computation comprises a shared phaser map Q and an activity map A. Each phaser map Q stores the available phasers, mapping phaser names to phasers. A phaser maps activity names to local views. A local view consists of the phase n the activity is in and an arrival status s set to ar when the activity arrives at the phaser. An activity map A maps activity names l to activities a. There are two kinds of activities: regular activities B; e consist of the bound for each phaser the activity is registered with, and an expression e; finish activities S; B; e additionally include a state S comprising the activities spawned within a finish instruction. Given any map, say X, we write dom X for the domain of X and range X for the co-domain of X. In a map X {x : u} we assume x does not occur in dom X.

Figure 4 introduces a small step reduction relation on states, S 1 → S 2 , capturing the non-deterministic choice of which activity to evaluate next. Auxiliary functions and predicates are in Figure 5. The phaser creation instruction, newPhaser n, evaluates to a fresh phaser name h; it also registers under h the current activity l with a local view composed of bound n and phase 0. All other phaser-related operations evaluate to (). Activities deregister from a phaser h via an expression drop h, thus removing the local view from the phaser. Instruc- Instruction awaitAll waits until activity l becomes unblocked. The current phase of a phaser is the natural number corresponding to the smallest local view among all activities registered in the phaser (Figure 5). An activity l is unblocked when it has arrived at all phasers and each bound allows progress (i.e., the bound is larger than the difference between the current phase and the phaser's phase). For each phaser activity l is registered with, rule R-await advances the phase and sets the arrival status back to un, via function commit(l, Q). HJ and X10 implicitly arrive at all non-arrived barriers before advancing.

h ∈ dom Q Q ; A {l : B; let x =
; let x = () in e } (R-arrive) unblocked(Q, l, B) Q ; A {l : B; let x = awaitAll in e } → commit(l, Q) ; A {l : B; let x = () in e } (R-await) Q ; A {l : B; let x = skipAll in e } → commit(l, Q) ; A {l : B; let x = () in e } (R-skip) l2 ∈ dom A ∪ dom B1 ∪ dom B2 Q ; A {l1 : B1; let x = async B2 e2 in e1 } → copy(l1, l2, dom B2, Q) ; A {l1 : B1; let x = () in e1 } {l2 : B2; e2 } (R-async) l2 ∈ dom A ∪ dom B Q; A {l1 : B; let x = finish e1 in e2 } → Q; A {l1 : ∅; {l2 : ∅; e1 } ; B; let x = () in e2 } (R-finish) S1 → S2 Q; A {l : S1; B; e } → Q; A {l : S2; B; e } (R
Expression skipAll simply advances the phase and does not wait for other activities. This operation can be used to let an activity use a phaser repeatedly without waiting for others. Spawning a new activity with rule R-async

unblocked(Q, l, ∅) unblocked(Q, l, B) Q(h) = P P (l) = n; ar m > n -phase(P ) unblocked(Q, l, B {h : m}) Phase commit partial function, commit(l, Q) = Q: commit(l, Q {h : P {l : n; ar }}) = commit(l, Q) {h : P {l : n + 1; un }} l / ∈ dom P commit(l, Q {h : P }) = commit(l, Q) {h : P } commit(l, ∅) = ∅ Local view copy partial function, copy(l1, l2, H, Q) = Q: h ∈ H P (l1) = r copy(l1, l2, H, Q {h : P }) = copy(l1, l2, H, Q) {h : P {l2 : r}} h / ∈ H copy(l1, l2, H, Q {h : P }) = copy(l1, l2, H, Q) {h : P } copy(l1, l2, H, ∅) = ∅
Halted state predicate, halted(S): halted( Q; {l1 : ∅; v1 , . . . , ln : ∅; vn } ) Fig. 5: Phaser-related functions and predicates.

evaluates expression e concurrently, by augmenting the activity map with a new activity B 2 ; e . For each phase in bounds B 2 , we copy the local views from the spawning activity l 1 to the spawned activity l 2 , as captured by function copy(l 1 , l 2 , H, Q), where H is a set of phaser names, defined in Figure 5. Spawned activities inherit the arrival statuses, for otherwise, depending on the order of reduction of the spawning and spawned activities, the spawned activity may or may not participate in the synchronisation of the current phase, inducing an undesirable non-determinism in the phaser semantics [START_REF] Martins | Types for X10 Clocks[END_REF].

Rule R-finish evaluates expressions of the form let x = finish e 1 in e 2 by suspending expression let x = () in e 2 and by evaluating e 1 in a newly created state. The finish activity (a triple as in Figure 3) holds a state (comprising an empty phaser map and an activity ∅; e 1 that is not registered on any phaser), the current bounds B, and the suspended expression let x = () in e 2 . Such an activity will then reduce, via rule R-activity, until halted (a predicate introduced in Figure 5). Then, the suspended activity resumes execution by means of rule R-join.

The evaluation of let-expressions is standard. Rule R-let replaces variable x by value v in continuation e. Nested let bindings are unfold with rule R-unfold.

We complete this section with a pair of examples leading to deadlocks, thus motivating the need for the type system in the next section. An activity that, before terminating, does not deregister from every phaser it is registered with will cause every activity that synchronises with that phaser to deadlock. Consider a program composed of two activities. Activity l1 creates a phaser x (line 2) and in the subsequent line spawns an activity l2 also registered with x; the latter activity does nothing (not even deregisters from phaser x). Activity l1 synchronises and deadlocks at line 5, forever waiting for activity l2 to arrive. An activity that forgets to arrive before awaiting other activities deadlocks itself and the remaining participants in the synchronisation. In the next program, activity l1 creates a phaser x and spawns activity l2 that simply arrives at x (line 4) and awaits for l1 (line 5). Activity l1 forgets to arrive at x but still awaits (line 6) and therefore deadlocks along with activity l2 (in line 5).

1 // activity l1 2 let x = newPhaser 0 in 3 async {x:0} ( // activity l2 4 arrive x; 5 awaitAll); // l2 deadlocks here 6 awaitAll // l1 deadlocks here Listing 1.2: An activity that forgets to arrive at a phaser.

The deadlocked state: {h : {l 1 : 0; un , l 2 : 0; ar }}; {l 1 : {h : 0}; let w = awaitAll in w , l 2 : {h : 0}; let z = awaitAll in z }

Type system and results

This section introduces our type system and its main results, namely type preservation and progress for typable states.

We rely on a set of type variables, ranged over by α. The syntax of types is defined by the grammar in Figure 6, and include those for the unit constant, unit, for the natural numbers, nat, and for phasers, α. We assign a different (singleton) type α to each phaser, in order to track how phasers are used throughout the program. The type system for our programming language is also defined in Figure 6. Typings Γ are maps from variables and phaser names to types. Arrival maps Φ map type variables (singleton phaser types) into arrival status. The re- lation for well-formed types Φ τ ensures that activities only make use of the phasers they are registered with. The typing rules for bounds Γ ; Φ 1 B : Φ 2 assign an arrival map Φ 2 to bounds B, under a context consisting of a typing Γ and an arrival map Φ 1 . Rule T-bound-cons ensures that B maps phasers into natural numbers, and also that it holds distinct phasers. The fact that phasers are associated to singleton types enables us to track aliasing in bounds.

The typing rules for values Γ ; Φ v : τ are straightforward. Rule T-val asserts that the value, either a variable or a phaser name, must be in typing Γ and that its type well formed.

For expressions we define a type and effect system Γ ; Φ 1 e : (τ, Φ 2 ) stating that expression e is of type τ and effect Φ 2 . The effects are important to track the changes in the arrival map, forced by the evaluation of an expression. The type of a newPhaser term is a new singleton type α that is also introduced in the effect. All remaining terms are of type unit. The effect a drop v term is the incoming arrival map from which α was removed, so that value v cannot be further used (cf. rule T-val). Rule T-arrive ensures that activities can arrive at phaser α only once, by requiring that the phaser's arrival status is un and by changing it into ar.

Terms awaitAll and skipAll mark the end of a phase: the rules check that all phasers have arrived and then reset the phasers to un. For example, in line 6 of Listing 1.2, activity l1 does not arrive at x, so we must type the term awaitAll under a typing {x : α} and an arrival map {α : un} which does not succeed according to rule T-await.

In rule T-async, the spawned activity e is checked against an arrival map Φ 2 for the phasers in B. Furthermore, e must deregister from all its phasers before terminating (hence the empty effect for e). The effect of the finish itself is the incoming phaser map, so that the spawning activity inherits the arrival status of the phasers. For example, we reject the program in Listing 1.1, since the spawned activity does not drop all its phasers before terminating. In fact, the empty task in line 3 must be typed under context {x : α} and phaser map {α : un}, but the arrival map is not empty for the unit term ().

Rule T-finish also forces e to deregister from all the phasers it has created, and therefore finish e has no effect on the arrival map Φ. In order to avoid deadlocks, we prevent e from accessing any existing phaser, thus eliminating (nested) dependencies between phasers and finish. The typing rule for let is standard.

The typing rules for states are introduced in Figure 7. We rely on a set of activity names Λ, a phase difference map ∆ mapping pairs of activity names to integer values (not necessarily natural numbers), and a phase difference tree map Σ mapping activity names to pairs composed of a phase difference map (for the root activity) and a phase difference tree map (for the children activities, if any). A state Q; A can be seen as a set of activity trees, the trees in A. Regular activities B; e are leaf nodes, whereas finish activities S; B; e are internal nodes whose children are the activities in S. When type checking a state, the topology of the phase difference tree Σ matches that of the activity tree.

We use ∆ P to check that the difference of phases between activities l i and l j are recorded in ∆(l i , l j ), for any pair l i , l j registered with P . Judgement Γ l Q : Φ collects the arrival statuses of every phaser activity l is registered with. Judgement ∆; Λ Q : Γ checks the phase difference and the registered activity names for each phaser in Q, while building a context Γ for Q.

We have two rules for activities. For regular activities, rule T-act ensures that the bounds B and the arrival map Φ mentions the same phasers, while ensuring that expression e deregisters from all phasers when before terminating (the effect of typing the expression is the empty arrival map). For finish activities, rule T-f-act ensures that both the state S and the finish continuation B; e are Phase difference for phasers, ∆ P :

∆(li, lj) = ni -nj ∀1 ≤ i, j ≤ k ∆ {l1 : n1; , . . . l k : n k ; } (T-dif)
Arrival map of a phaser map, Γ l Q : Φ: Typing rules for activity maps, Σ; Γ Q A:

Γ l Q : Φ Γ (h) = α P (l) = ; s Γ l (Q {h : P }) : (Φ {α : s}) Γ l Q : Φ l / ∈ dom P Γ l (Q { : P }) : Φ Γ l ∅ : ∅ (
Γ l Q : Φ ∆1; Σ1; Γ ; Φ a Σ; Γ Q A Σ {l : ∆1; Σ1 }; Γ Q A {l : a} ∅; Γ Q ∅ (T-am-cons,T-am-nil)
Typing rule for states, ∆; Σ S:

∆; dom A Q : Γ Σ; Γ Q A ∆; Σ Q; A (T-state)
Fig. 7: Typing rules for states well typed. To type check a state ∆; Σ Q; A , rule T-state uses the activity names in A and the phase difference ∆ to type check the phaser map Q; it also checks that the activity map A is well typed according to the phase difference tree Σ.

We complete this section by presenting the main results of the paper.

Lemma 1. If Σ {l : ∆ 1 ; Σ 1 }; Γ Q A {l : a}, then there exists Φ such that Γ l Q : Φ, ∆ 1 ; Σ 1 ; Γ ; Φ a, and
Σ; Γ Q A.
Theorem 1 (Subject reduction). If ∆ 1 ; Σ 1 S 1 and S 1 → S 2 , then there exists ∆ 2 and Σ 2 such that ∆ 2 ; Σ 2 S 2 .

Proof (Sketch). The proof follows by induction on the derivation of the reduction step. In each case we have to exhibit ∆ 2 and Σ 2 . For R-phaser we make use of a weakening lemma. Cases R-drop, R-arrive, R-await, R-skip, and R-async are similar. We prove that changes made in the phaser map after reduction have no effect on any activity besides the one under reduction. When the derivation of the reduction step ends with rule R-activity, we know that Σ 1 = Σ {l : ∆ 1 ; Σ 1 }; by induction it follows that ∆ 2 ; Σ 2 S 2 . We take ∆ 2 = ∆ 1 and Σ 2 = Σ {l : ∆ 2 ; Σ 2 }. We apply Lemma 1 to the hypothesis to, and the induction hypothesis to complete the proof. Case R-join, and R-unfold are similar to R-activity. For R-finish, we know that Σ 1 = Σ {l 1 : ∅; ∅ }; we take ∆ 2 = ∆ 1 and Σ 2 = Σ {l 1 : ∅; {l 2 : ∅; ∅ } }. We use a strengthening lemma: Γ ; ∅ e : (τ, ∅) and Γ only contains phaser names in its domain implies ∅; ∅ e : (τ, ∅). Strengthening is applied to the newly created state. For R-let we need a substitution lemma, strengthening and weakening.

For progress, we start by extracting a total order for the activities in a typable state.

Lemma 2. If ∆; Σ Q; then the relation {(l 1 , l 2 ) | P ∈ range Q, P (l 1 ) = n 1 ; , P (l 2 ) = n 2 ; , n 1 ≤ n 2 } is a total order.
Theorem 2 (Progress). If ∆; Σ S 1 then S 1 is either halted(S 1 ) or there is a state S 2 such that S 1 → S 2 .

Proof (Sketch). Activities can block for a number of reasons, easily deduced from the reduction rules in Figure 4. The cases for drop, arrive, skipAll, async and let are easily dismissed by a simple analysis of the typing derivation rules. For example, when the reduction step ends with rule R-async, we must show that ∆; Σ Q; A {l 1 : B 1 ; let x = async B 2 e 2 in e 1 } and l 2 ∈ dom(A, B 1 , B 2 ) implies that copy(l 1 , l 2 , dom B 2 , Q) is defined. We proceed by induction on the structure of Q. The interesting case is when Q is Q {h : P }, where we must show that l 1 ∈ dom P . From ∆; Σ S we know that Γ l1 Q : Φ. By showing that Γ l Q : Φ implies dom Q = dom Γ , we obtain h ∈ dom Γ . Then we show that h ∈ dom Γ and Γ l Q : Φ implies l ∈ dom Q(h).

Otherwise suppose that all activities in S 1 are blocked at awaitAll. Lemma 2 ensures that there is a total order on activity names. Since the order is total, there is one activity name that is smaller than all others; let it be l. From ∆; Σ Q; A {l : B; let x = awaitAll in e 1 } , we know that P (l) = m; ar , for some m. We also know that phase(P ) = min{n + 1 | n; ar ∈ range P } = 1 + min{n | n; ar ∈ range P } = 1 + m. Hence l is unblocked and we have attained a contradiction.

In absence of infinite computations, it follows from the above theorem that all typable states eventually reach a halted state.

Conclusion and further work

We presented a calculus and a type system for a fork/join programming model with a flexible phaser mechanism. We favour explicit operations, yielding phaser operators which are simpler than those of Habanero Java [START_REF] Shirako | Phasers: a unified deadlock-free construct for collective and point-to-point synchronization[END_REF]. Our proposal unifies the semantics of clocks [START_REF] Charles | X10: an object-oriented approach to non-uniform cluster computing[END_REF], regular phasers [START_REF] Shirako | Phasers: a unified deadlock-free construct for collective and point-to-point synchronization[END_REF], and phaser beams [START_REF] Shirako | Phaser beams: Integrating stream parallelism with task parallelism[END_REF], but goes further by allowing tasks to be ahead of others by a bounded number of phases.

HJ permits writing applications where the same task is registered with a phaser as a regular barrier and with another phaser that disregards all forms of synchronisation. For such cases, our skipAll operation is not enough. We can however introduce unbounded local views for tasks registered at phasers with an infinite bound (tasks that do not want to synchronise, only to influence others). An expression advance v, available only for activities registered with an infinite bound, would advance a single phaser. We believe that such an extension can be easily accommodated in our system.

To focus on the intricacies of synchronisation, we kept our language very simple. Extensions required for real world programming include provision for unbounded computations (in the form of recursion or loops) and for a mutable store (in the form of imperative variables). Loops can be introduced in our calculus while causing little interference with the model and results, as long as they preserve an invariant on the registered phasers at every iteration. In order to build circular buffers of phasers, HJ applications may create phasers within a loop. A possible workaround to accommodate such a feature in our language is to introduce a primitive that allocates an array of phasers.

Fig. 1 :

 1 Fig. 1: A parallel breadth-first search algorithm.

||Fig. 2 :Fig. 3 :

 23 Fig. 2: The syntax of expressions.

  newPhaser n in e } → Q {h : {l : 0; un }} ; A {l : B {h : n}; let x = h in e } (R-phaser) Q {h : (P {l : })} ; A {l : B {h : }; let x = drop h in e } → Q {h : P } ; A {l : B; let x = () in e } (R-drop) Q {h : (P {l : n; un )}} ; A {l : B; let x = arrive h in e } → Q {h : (P {l : n; ar )}} ; A {l : B

Fig. 4 :

 4 Fig. 4: Small step semantics for states, S → S.

  Phase function for local views, phase(r) = n: phase( n; un ) = n phase( n; ar ) = n + 1 Phase partial function for phasers, phase(P ) = n: phase(P ) = min{phase(r) | r ∈ range P } Unblocked predicate, unblocked(Q, l, B):

1 // activity l1 2

 2 let x = newPhaser 0 in 3 async {x:0} (); // forgets drop x 4 arrive x; 5 awaitAll // l1 deadlocks here Listing 1.1: An activity that forgets to deregister from a phaser. The deadlocked state: {h : {l 1 : 0; ar , l 2 : 0; un }}; {l 1 : {h : 0}; let w = awaitAll in w , l 2 : {h : 0}; () }, }

ΓFig. 6 :

 6 Fig. 6: The syntax of types and the typing rules.

  T-ar-cons,T-ar-skip,T-ar-nil)Typing context of a phaser map, ∆; Λ Q : Γ :α ∈ range Γ dom P ⊆ Λ ∆ P ∆; Λ Q : Γ ∆; Λ (Q {h : P }) : (Γ {h : α}) ∆; Λ ∅ : ∅ (T-pm-cons,T-pm-nil) Typing rules for activities, ∆; Σ; Γ ; Φ a: Γ ; Φ B : Φ Γ ; Φ e : ( , ∅) ∅; ∅; Γ ; Φ B; e ∆; Σ S ∅; ∅; Γ ; Φ B; e ∆; Σ; Γ ; Φ S; B; e (T-act,T-f-act)

The syntax of types:τ ::= unit | nat | αWell-formed types, Φ τ :Φ unit Φ nat Φ, α : s α (T-wf-u, T-wf-n, T-wf-p)Typing rules for bounds, Γ ; Φ B : Φ:Γ ; Φ1 B : Φ2 Γ ; Φ1 v1 : α Γ ; Φ1 v2 : nat Γ ; Φ1 (B {v1 : v2}) : (Φ2 {α : Φ1(α)}) Γ ; Φ ∅ : ∅ (T-bound-cons,T-bound-nil)Typing rules for values, Γ ; Φ v : τ :Γ ; Φ () : unit Γ ; Φ n : nat Γ (v) = τ Φ τ Γ ; Φ v : τ (T-unit, T-nat, T-val)Typing rules for terms and for expressions, Γ ; Φ t : (τ, Φ) and Γ ; Φ e : (τ, Φ):Γ ; Φ v : nat α ∈ dom Φ Γ ; Φ (newPhaser v) : (α, Φ {α : un}) (T-phaser) Γ ; Φ {α : s} v : α Γ ; Φ {α : s} (drop v) : (unit, Φ) Γ ; Φ {α : un} v : α Γ ; Φ {α : un} (arrive v) : (unit, Φ {α : ar}) (T-drop,T-arrive) Γ ; {α1 : ar, . . . , αn : ar} awaitAll : (unit, {α1 : un, . . . , αn : un}) (T-await) Γ ; {α1 : ar, . . . , αn : ar} skipAll : (unit, {α1 : un, . . . , αn : un}) (T-skip) Γ ; Φ1 B : Φ2 Γ ; Φ2 e : ( , ∅) Γ ; Φ1 (async B e) : (unit, Φ1) Γ ; ∅ e : (τ, ∅) Γ ; Φ (finish e) : (unit, Φ) (T-async,T-finish)

Acknowledgements. This work was partially supported by project PTDC/EIA-CCO/122547/2010. The first author would like to thank Vivek Sarkar for welcoming him at the Habanero group at Rice University, during the year of 2012. We are grateful to Jun Shirako and anonymous referees for their feedback on this paper and to Vivek Sarkar for discussions related to phasers.