Riccardo Rastelli
email: riccardo.rastelli@ucdconnect.ie

Pierre Latouche

Nial Friel

Choosing the number of groups in a latent stochastic block model for dynamic networks

Keywords: Stochastic Block Models, Dynamic Networks, Greedy Optimisation, Bayesian Inference, Integrated Completed Likelihood

Latent stochastic block models are flexible statistical models that are widely used in social network analysis. In recent years, efforts have been made to extend these models to temporal dynamic networks, whereby the connections between nodes are observed at a number of different times. In this paper we extend the original stochastic block model by using a Markovian property to describe the evolution of nodes' cluster memberships over time. We recast the problem of clustering the nodes of the network into a model-based context, and show that the integrated completed likelihood can be evaluated analytically for a number of likelihood models. Then, we propose a scalable greedy algorithm to maximise this quantity, thereby estimating both the optimal partition and the ideal number of groups in a single inferential framework. Finally we propose applications of our methodology to both real and artificial datasets.

Introduction

In the last few years, there has been an increasing amount of data stored characterising interactions between individuals or, more generally, units of interest. An interaction, represented by a triple (i, j, η), indicates that units i and j have a connection at a specific time point η. Interactions can for instance describe email exchanges between individuals or posts on social media. In Biology, units can correspond to genes and interactions to regulation events between the genes. A natural approach to model the set of all observed interactions is to rely on a dynamic graph where each unit is associated with a node and an edge (i, j) is present between two units i and j, at time η, if the interaction (i, j, η) is recorded.

A long series of methods have been proposed recently to cluster the nodes of dynamic networks in order to summarise the information hidden in such data. A large number of these consider the Stochastic Block Model (SBM) [START_REF] Wang | Stochastic blockmodels for directed graphs[END_REF][START_REF] Nowicki | Estimation and prediction for stochastic blockstructures[END_REF] as a starting point and propose extensions to the dynamic framework.

Moreover, they are usually discrete in time, that is, predefined time intervals are introduced and interactions during those time intervals are aggregated to obtain snapshots indexed by a discrete time variable t. In the case of a binary dynamic network, two nodes are connected in a snapshot if they have at least one interaction in the corresponding time interval. We recall that the SBM was originally introduced to cluster nodes in static networks where nodes and edges do not evolve through time. The model assumes that nodes are spread in latent clusters, which, in practice, have to be inferred from the data.

The probability for two nodes to connect is then only dependent on their respective clusters. Because no constraints are imposed on the connection probabilities, various types of clusters can be extracted from the data, which makes methodologies based on SBM applicable to networks with different connectivity structures [START_REF] Daudin | A mixture model for random graphs[END_REF]). The dynamic SBM-like model of [START_REF] Yang | Detecting communities and their evolutions in dynamic social networks -a Bayesian approach[END_REF] allows, for example, each node to switch its cluster at time t depending on its state at time t -1. A transition matrix is employed to characterise the switching probabilities. While clusters can change over time, fixed connection probabilities are used. Conversely, the model of [START_REF] Xu | Dynamic stochastic blockmodels for time-evolving social networks[END_REF] relies on evolving connection probabilities and temporal changes are described through a state space model. Therefore, the inference requires the use of optimisation tools such as the Kalman filter and the Rauch-Tung-Striebel smoother. The work of [START_REF] Yang | Detecting communities and their evolutions in dynamic social networks -a Bayesian approach[END_REF] was then extended by [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF] for the clustering of nodes in dynamic networks where edges are not necessarily binary. We emphasise that theoretical results are also provided in their paper to show that dynamic SBM-like models should not let both the clusters and connection parameters evolve through time without incurring into identifiability issues.

In the static framework, many extensions have been considered for the SBM model. Many of them have then been adapted to deal with dynamic networks. For example, the well-known mixed-membership SBM of [START_REF] Airoldi | Mixed membership stochastic blockmodels[END_REF] has been extended to a dynamic framework by several works including [START_REF] Xing | A state-space mixed membership blockmodel for dynamic network tomography[END_REF], [START_REF] Ho | Evolving cluster mixed-membership blockmodel for time-evolving networks[END_REF] and [START_REF] Kim | Nonparametric multi-group membership model for dynamic networks[END_REF]. Note that the latent position model of [START_REF] Hoff | Latent space approaches to social network analysis[END_REF], which is also popular in the network community, was also adapted by [START_REF] Sarkar | Dynamic social network analysis using latent space models[END_REF] and [START_REF] Friel | Interlocking directorates in Irish companies using a latent space model for bipartite networks[END_REF] to deal with dynamic interactions.

In this paper, our objective is to model the evolution of the cluster memberships over time by relying on a Markovian property. We take advantage of prior conjugacy to integrate out analytically most of the model parameters, with an approach similar to that of [START_REF] Mcdaid | Improved Bayesian inference for the stochastic block model with application to large networks[END_REF] and [START_REF] Côme | Model selection and clustering in stochastic block models based on the exact integrated complete data likelihood[END_REF]. This so-called collapsing allows one to obtain an analytical expression for the integrated completed data likelihood for a number of likelihood models. A greedy optimisation algorithm is then employed for inferential purposes. It allows estimation of both the number of clusters and the cluster memberships of the nodes to the clusters.

Section 2 introduces the model and the notation. A Bayesian hierarchical structure is then proposed in Section 3. Finally, the optimisation procedure is given in Section 5 and experiments are carried out in Sections 6, 7, and 8 to assess the proposed methodology.

Markovian Stochastic Block Model

The observed data consist of a sequence of network objects X = X (t) t∈T defined on the same set of nodes V = {1, . . . , N }. We focus our attention on dynamic networks, where T = {1, . . . , T } denotes the temporal span.

The random variable X (t) ij models the value exhibited by the edge from i to j at time t, ∀t ∈ T and ∀i, j ∈ V. We outline our methodology on directed networks, however it applies also to undirected network as we illustrate in the applications we consider. A typical scenario is that of a binary dynamic network, where:

x (t) ij =
1, if an edge from i to j appears at time t; 0, otherwise.

(1) Also, self-edges are not allowed, i.e. x

(t)

ii = 0, ∀t ∈ T and ∀i ∈ V. The random variable Z (t) i characterises the allocation of a node, whereby

Z (t) i = g
indicates that node i is allocated to group g at time t, for a certain g ∈ K = {1, . . . , K}.

The set

Z (t) = Z (t) i i∈V
therefore corresponds to a random hard clustering (i.e. a partitioning) of V, ∀t ∈ T . We also denote the full set of allocations with Z = Z (t)

t∈T

.

The total number of groups in this underlying structure is denoted by K, hence Z (t) i ∈ K, ∀t ∈ T and ∀i ∈ V. Note that K represents the total number of groups in Z, implying that, at each time frame t, the current number of non-empty groups K (t) may be any number in K.

The allocations characterise the connection profiles of the nodes of the network, in that nodes allocated to the same group have their edges drawn from the same probability distribution. The edges of the dynamic network are marginally dependent, however they are conditionally independent given the allocations:

p (X |Z) = t∈T p X (t) Z (t) , (2)
and ∀t ∈ T :

p X (t) Z (t) = i∈V j∈V:j =i p X (t) ij Z (t) i , Z (t) j
.

(3)

The distribution for the value of a single edge is determined by the allocations of the nodes, as follows:

p X (t) ij = x z (t) i = g, z (t) j = h, Θ = f (x; θ gh ) , (4)
where Θ and θ gh are collections of model parameters. In the binary network case, f corresponds to the mass probability function of a Bernoulli variable with success probability

θ gh ∈ [0, 1].
In such case the likelihood of the model may be written as:

L X = g∈K h∈K t∈T i∈V: z (t) i =g j∈V: j =i; z (t) j =h θ x (t) ij gh (1 -θ gh ) 1-x (t) ij .
(5)

Concerning the modelling of the temporal evolution of the network, a Markovian property on the nodes' allocations is adopted. The sequences

Z i = Z (t) i t∈T
are assumed to be independent Markov chains realised using the same K × K kernel matrix Π. Therefore the prior structure on the allocations factorises as follows:

p (Z|Π) = p Z (1) t∈T \{1} p Z (t) Z (t-1) , Π , (6)
where, ∀t ∈ T \ {1} and ∀i ∈ V:

p Z (t) i = h Z (t-1) i = g, Π = π gh ∈ [0, 1]. (7)
Note that the rows of Π must sum to 1, whereas no constraint is imposed on the columns of the same matrix. Combining ( 6) and ( 7), the prior on the allocations may be alternatively written as:

p (Z|Π) = p Z (1) g∈K h∈K π R gh gh , (8)
where R gh denotes the total number of switches from group g to h.

Bayesian hierarchical structure

We propose a hierarchical structure to further model the parameters θ and Π. In particular, we focus on the use of conjugate priors, since these allow one to integrate out (collapse) most of the model parameters.

Concerning the transition probabilities, we specify a Dirichlet distribution over the rows of Π:

(π g1 , . . . , π gK ) ∼ Dir (δ g1 , . . . , δ gK ) , ∀g ∈ K.

By analytically integrating out the parameters {π gh } g,h , the following compound distribution ("marginal prior") for the allocations arises:

p (Z|δ) = p Z (1) g∈K h∈K Γ (δ gh + R gh ) Γ h∈K δ gh + h∈K R gh Γ h∈K δ gh h∈K Γ (δ gh ) . ( 10
)
We assume a uninformative flat prior distribution for the hyperparameters {δ gh } g,h and thus fix all of them to 1.

A separate and independent Multinomial-Dirichlet model may be specified for Z (1) .

However, for these starting allocations, we opt for a more pragmatic and parsimonius approach: we use a Multinomial distribution with parameters {α g } g∈K , where:

α g ∝ t∈T \{1} i∈V 1 Z (t) i =g , ∀g ∈ K. ( 11
)
The distribution on the initial states so defined approximates the stationary distribution

associated to Π. Γ (a + b) Γ (a) Γ (b) Γ (a + η gh ) Γ (b + N gh -η gh ) Γ (a + b + N gh ) , ( 12
)
where a and b are the Beta hyperparameters, η gh counts the number of edges from group g to h:

η gh = t∈T i∈V: z (t) i =g j∈V: j =i; z (t) j =h X (t) ij (13)
N gh = t∈T i∈V: z (t) i =g j∈V: j =i; z (t) j =h 1. ( 14
)
Note that, similarly to the hyperparameter δ, we make a "symmetric" assumption on the hyperparameters a and b, making them identical over all of the groups. In the applications, these are both fixed to 1 yielding a uninformative flat prior distribution.

Figure 1 shows a graphical model summarising the hierarchical structure and the dependencies between the variables introduced.

δ φ π θ A similar modelling framework has been recently proposed by [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF]. In this paper, we take advantage of some of their results on model-identifiability and we propose a radically different estimation technique.

Z (t -1) Z (t) X (t -1) X (t)

Incomplete weighted networks

When dealing with weighted networks, graphs are often incomplete, in that not all the edges exhibit a value. The presence-absence of edge values can be modelled using a

Bernoulli random variable ρ ij , which is equal to one if the edge from i to j is present and is 0 otherwise. This adds an extra layer in the modelling, which is unrelated to the likelihood model specified on the present edge values. In fact, any likelihood structure can be used to model the edges that appear (ρ ij = 1) in the network.

Assuming a Beta prior on the presence-absence indicator, the marginal likelihood of the network is simply the product of the two compound distributions:

p (X |Z, a ρ , b ρ , φ) = p (ρ|Z, a ρ , b ρ ) p (X |Z, ρ, φ) ; (15)
where (a ρ , b ρ ) and φ denote the hyperparameters for ρ and θ, respectively. Note that the modelling on ρ can capture the heterogeneity induced by the block structure, in that the probability of an edge exhibiting a value may depend on its allocation.

Exact Integrated Completed Likelihood

In the previous section, we have shown that in a binary Markovian dynamic network the marginal likelihood p (X |Z, a, b) and the marginal prior p (Z|δ) have an exact analytical form. These terms can be recombined to obtain a particularly meaningful quantity, the so-called exact Integrated Completed Likelihood (ICL), defined as:

ICL ex = log [p (X |Z, a, b)] + log [p (Z|δ)] . (16)
Such a quantity corresponds to the exact value that the ICL criterion of [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] propose to maximise when choosing the number of groups in a finite mixture context. We stress that, although the definition in ( 16) relates to the specific case of binary dynamic networks, the same quantity can be evaluated analytically for all of the likelihood models listed in Table 1.

We propose to use ICL ex as a model-based optimality criterion for the clustering problem on the nodes of our dynamic network: in the space of all possible allocations, we seek a Ẑ maximising ICL ex . Note that, thanks to the discrete nature of the allocation variables, the optimal total number of groups K can be deduced automatically from Ẑ, along with the values K (t) t∈T .

We also emphasise that, from a Bayesian perspective, Ẑ corresponds to a Maximum A Posteriori (MAP) solution:

p (X |Z, a, b) p (Z|δ) = p (X , Z|a, b, δ) ∝ p (Z|X , a, b, δ) , (17)
where the proportionality is intended with respect to Z.

Greedy optimisation

In the Markovian dynamic stochastic block model described, we seek a clustering solution Ẑ maximising the ICL ex value defined in (16). A complete exploration of all possible allocations is not feasible, even for small datasets. However, similar combinatorial problems have been recently tackled successfully using heuristic greedy routines, resembling the well known Iterated Conditional Modes of [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF] et al. (2016) and for Gaussian finite mixtures in [START_REF] Bertoletti | Choosing the number of clusters in a finite mixture model using an exact integrated completed likelihood criterion[END_REF]. Here, we propose an extension of the same ideas to our dynamic network context.

The algorithm repeatedly sweeps over the network's nodes and attempts to reallocate them using a greedy behaviour. In each step, a random time frame t and a random node i are chosen, and the variation of ICL ex is evaluated for all the possible reallocations of the corresponding node. For every group g, the value (t,i)→g is evaluated, corresponding to the ICL ex value after moving the node identified by (t, i) to g. Eventually the node is reallocated to the group that yields the best increase in the objective function. This procedure continues until no reallocation can provide a further increase. The corresponding pseudocode is shown in Algorithm 1.

As input, the algorithm requires a starting partition, the hyperparameters' values, and a parameter K up denoting the largest admissible number of groups. During the optimisation, groups may be deleted (if they remain empty at all time frames) and created (if a node is reallocated to an empty group). In practice, the former case is very frequent whereas the latter is rather rare. Hence the best performance is achieved when the starting partition is composed of K up groups, with K up set to a fairly large value.

The algorithm is bound to find only a local optimum for the objective function, and its initialisation plays a crucial role. Due to the random update order, it is possible

Algorithm 1 GreedyIcl Initialise Z (t) i , ∀t ∈ T , ∀i ∈ V. Evaluate ICL ex and set = ICL ex and stop = ICL ex . Set stop = f alse. while !stop do Set U = {(t, i) : t ∈ T , i ∈ V}. Shuffle the elements of U. while U is not empty do (t, i) = pop (U). ĝ = arg max g=1,2,...,Kup (t,i)→g . = (t,i)→ĝ . Z (t) i = ĝ. end while if ≤ stop then stop = true else stop = . end if end while Return Z and .
that different starting partitions return different solutions, however we find that good initialisations often lead to a unique maximum.

As concerns the complexity of the algorithm, in the binary case one iteration can be

performed in O M + T N K 2
up , where M is the total number of edges appearing in the network. Once a node (t, i) is selected, the number of edges to and from group g are counted, for every g ∈ {1, . . . , K up }. Evidently this implies a cost of O (m), where m is the average value for the sum of in-degree and out-degree of a random node. Then, the quantities (t,i)→g are evaluated for every g. For each of these, the computational bottleneck is given by the calculation of the variation of the marginal likelihood value, which can be performed in O (K up ). This makes the average cost of updating an allocation

O m + K 2
up . Since all of the allocations are repeatedly updated, the overall complexity of one iteration is O T N m + T N K 2 up as previously claimed.

Final merge step

Once the GreedyIcl has converged, we additionally propose a hierarchical clustering on the optimal solution, following the approach of [START_REF] Côme | Model selection and clustering in stochastic block models based on the exact integrated complete data likelihood[END_REF]. We consider all possible pairs of groups and attempt to merge each pair into a single cluster. Such a move is accepted only if the ICL ex value increases, and everytime a merge move is performed the procedure is restarted. For each try, the computational bottleneck is given by the evaluation of the variation for the marginal prior, which can be performed in O K2 . Since the number of pairs is O K2 , and the number of accepted merge moves is capped at K, the overall complexity of this final step is O K5 . Here K denotes the number of groups for the optimal solution obtained through the GreedyIcl, which is normally much smaller than K up . We find that in practice this final merge step does not impact the overall computational time by much, and the additional computational burden may be neglected.

Simulation study

We propose an experiment to validate our methodology and compare it to the existing algorithm of [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF], using artificial data. The number of time frames, the number of nodes, and the true underlying number of groups are fixed throughout as follows: T = 4, N = 50 and K = 4, respectively. As concerns the transition probabilities matrices, the following general structure is assumed:

Π =     π ν ν ν ν π ν ν ν ν π ν ν ν ν π     (18)
where ν = (1-π)/(K -1), so that rows sum to one. Choosing a high π value, the clusters tend to be very stable and allocations do not change often over time. By contrast, a small π value represents the opposite scenario, corresponding to a highly instable system. In either case, the stationary distribution, which is used to generate the starting partition,

is α = 1 K , 1 K , 1 K , 1 K .
In our simulations we consider two different scenarios: π = 0.7 and π = 0.9, which are denoted low-stability and high-stability, respectively. Although these two cases may offer only a limited view of all the situations encompassed by the model, we believe these to be the most interesting examples and that most realised networks are in fact well captured by this representation.

As concerns the connection probability matrices, an affiliation structure is assumed:

    θ 0 0 0 0 0 θ 0 0 0 0 0 θ 0 0 0 0 0 θ 0     . (19)
A small perturbation is added independently to each of the entries of such matrix:

θ = θ 0 + 0.1u = 0 + 0.1u, ( 20
)
where u is drawn from a Uniform distribution in the interval [-1, 1] for each entry. While 0 is always fixed to 0.1, θ 0 is different for each matrix and determines the difficulty level in recovering the underlying latent clustering. We consider 9 different scenarios, each corresponding to a different choice of θ 0 , ranging from 0.1 to 0.9. The perturbation created by the random variable u is necessary to ensure identifiability of the model, as explained in [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF].

The experiment is executed as follows: for every choice of π and θ 0 , 100 random dynamic undirected networks are generated, each corresponding to a different random realisation of the connection probability matrix. In each network, the variational algorithm of [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF], implemented in the R package dynsbm (version 0.3), and our implementation of the greedy algorithm are executed. As concerns the variational approach, dynsbm is run 4 times for every value of K in the set {1, 2, . . . , 6}, and the best run is retained as optimal solution according to the approximate ICL criterion advocated in [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF]. Higher values of K make dynsbm particularly time consuming and hence are not considered, however, we note that in all of the runs the approximate ICL criterion never favours the model with 6 groups. As concerns our greedy algorithm (denoted GreedyIcl), we consider several types of initialisations, as follows:

• aggregated: an adjacency matrix of size N × N is obtained, by aggregating (summing) the T adjacency matrices of the generated network. Then, the kmeans algorithm is run on such a matrix using as number of centres a random draw from the discrete interval 0.5 * N , . . . , 0.75 * N . This associates a cluster membership to each of the N nodes. These allocations are assumed not to change over time, so that the initial partition is given by the output of the kmeans repeated over all the time frames.

• colbind: the only difference with the aggregated initialisation lies in the fact that kmeans is run on a matrix obtained by gathering the adjacency matrices one next to the other, obtaining a N × T N matrix. Note that this is the same initialisation used for the dynsbm algorithm.

• rowbind: in this case kmeans is instead run on a matrix obtained by stacking up the observed adjacency matrices. Since the size of this matrix is T N × N , the number of groups considered by kmeans is not capped at N as in the other cases:

hence we use a draw from the discrete interval 0.5 * T * N , . . . , 0.75 * T * N as number of centres. Also, in contrast with the previous cases, kmeans returns the allocations for every node at every time frame, hence allocations are not assumed to be unchanged over time.

• random: the GreedyIcl is initialised unsing a random partition, with K up chosen as in GreedyIcl rowbind.

Once initialised, the GreedyIcl algorithm is run once for each type of initialisation. The additional label

• all: is used to indicate the best solution obtained through GreedyIcl over all possible types of initialisations.

The optimal clusterings obtained using each method are compared to the true allocations using the Normalised Mutual Information (NMI) index [START_REF] Strehl | Cluster ensembles -a knowledge reuse framework for combining multiple partitions[END_REF].

This index takes values in the interval [0, 1] and describes how similar two partitions are.

These partitions must be vectors, hence the optimal clusterings obtained are vectorised by concatenating the partitions time-wise.

Figure 2 shows the results obtained regarding the clustering performance. It appears that the greedy optimisation coupled with kmeans initialisations and dynsbm achieve the best results, whereas the greedy optimisation with random intialisations perform poorly.

The average computing times are provided in Table 2.

As concerns model selection, our methodology outperforms the existing algorithm of

Average Normalised Mutual Information

Cluster stability = 0.7

Intra-cluster connection probability Average NMI q q q q q q q q q q q q q q q q q q q q q q q q dynsbm GreedyIcl all GreedyIcl aggregated GreedyIcl rowbind GreedyIcl colbind GreedyIcl random

Average Normalised Mutual Information

Cluster stability = 0.9

Intra-cluster connection probability Average NMI q q q q q q q q q q q q q q q q q q q q q q q q dynsbm GreedyIcl all GreedyIcl aggregated GreedyIcl rowbind GreedyIcl colbind GreedyIcl random

Figure 2: Simulation study. Average Normalised Mutual Information index between the true clustering and the optimal clustering solutions obtained through dynsbm (version 0.3) and several versions of GreedyIcl. Shaded regions represent the 90% quantile region associated to each scenario. These regions are plotted only for GreedyIcl all (black colour) and dynsbm (blue colour). It appears that the greedy algorithm with a good initialisation performs as well as the existing algorithm dynsbm.

Table 2: Simulation study. Average computing time required by each of the algorithms to be run on a generated dynamic network. Note that the variational expectation-maximisation of dynsbm and the greedy updates of our algorithm are both run exactly 4 times for each dataset. The initialisation may be different, but we find that this does not impact the computing time by much.

Algorithm

Average computing time (seconds)

dynsbm

Model-choice performance

Cluster stability = 0.7

Intra-cluster connection probability Proportion of correct estimations of K q q q q q q q q q q q q q q q q q q q q q q q q dynsbm GreedyIcl all GreedyIcl aggregated GreedyIcl rowbind

Model-choice performance

Cluster stability = 0.9

Intra-cluster connection probability Proportion of correct estimations of K q q q q q q q q q q q q q q q q q q q q q q q q dynsbm GreedyIcl all GreedyIcl aggregated GreedyIcl rowbind GreedyIcl colbind GreedyIcl random

Figure 3: Simulation study. For each combination of π and θ 0 , the proportion of networks where K is properly estimated is shown, for all of the algorithms considered.

We construct a binary directed dynamic network of emails transforming the data into an adjacency cube X , such that:

x (t) ij =
1, at least an email is sent from member i to member j at time frame t; 0, otherwise. We run our GreedyIcl all algorithm as in the simulation study, with each single run of the algorithm taking on average about 400 seconds. The overall best solution has 17

groups. An analysis of this clustering solution follows.

Activity levels and connectivity. As shown in the left panel of Figure 5, the clustering solution has very high stability over time, in that nodes do not change allocations frequently. The connection probability matrix, shown on the right panel, exhibits instead a much more complex situation. Evidently, a strong community structure is present, q q q q q q q q q q q q q q q q q q q q q q q q q q q 2000 2001 2002 On the right panel, the estimated connection probabilities are shown. Most groups have high within cluster connectivity, although a number of off-diagonal entries are also relevant. The groups are ordered according to their aggregated size (i.e. the sum of their sizes over all times), in descending order. Some entries cannot be properly estimated due to a small cluster size, hence are left blank.

since most of the diagonal elements of this matrix are fairly large.

We propose a brief characterisation of the role played by each group, using the information provided by this connection probability matrix, in conjunction with Tables 3,4 and 5.

• Group 1: this group contains inactive nodes only. These members do not send emails, but they may receive a few newsletters.

• Groups 2 to 8: these contain a massive portion of the nodes, and hence describe the most common connectivity profiles that may arise. The nodes in these groups are not particularly active (the connection probabilities are typically not large) and tend to receive more emails than they send. The sizes of these groups swell increasingly by acquiring nodes from group 1 until the bankruptcy, signalling increased activity in the months before the default.

• Groups 9, 11, 12, 13: these groups are mainly composed of managers, and correspond to active nodes, who send more than they receive. These groups represent the different profiles of the executives of the company. Also note that these groups have a high within cluster connectivity, hence it is reasonable to find the corresponding conversations particularly meaningful in terms of intelligence and company directives.

• Groups 10, 14, 15: three groups of very active nodes, although these are mainly composed of employees. Note that groups 10 and 15 only become relevant during the year 2001.

• Group 16: this group actually only ever contains one person (with a non-specified role). This person apparently has some special position and hence may be considered as an outlier.

• Group 17: the smallest of groups, containing 4 unique members, each for 1 timeframe only. Since this group has high connection probabilities towards many groups, it is reasonable to believe that nodes join this group whenever they send out newletters to all members. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

Avg. out-degree 0 3 3 1 3 4 3 5 5 3 7 10 13 12 11 11 57 Avg. in-degree 1 3 3 3 4 5 4 5 4 2 6 5 7 7 5 7 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q 2000 2001 2002

Proportions in group 1 by status

Time Proportion

Oct '01 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Employees Managers Figure 6: Enron dataset. On the left panel, the evolution of the number of non-empty groups is shown. This quantity appears to increase until the collapse, signalling an increase in heterogeneity in the network. On the right panel, the proportions of managers and employees allocated to the first group are shown. These proportions correspond to the nodes that are not particularly active.

Temporal dynamics. Our model allows for groups to become inactive, whenever their sizes decrease to zero. This means that the connection profile associated to a particular group is not exhibited by any of the nodes in the time frame considered. Hence, the number of non-empty groups can be used as a measure of heterogeneity in the network, and the temporal dynamics of this quantity can be used to assess how heterogeneity evolves over time. In the left panel of Figure 6, the evolution of the number of non-empty groups is shown. It appears that the network is particularly heterogeneous in the year before the collapse (K (t) = 15 or 16), whereas it eventually becomes rather homogeneous (K (t) = 10) after the default. A similar message is conveyed by the plot on the right panel of Figure 6. Here the number of inactive nodes (i.e. nodes allocated to group 1) is shown to decrease convincingly for both employees and managers. Eventually, after the collapse, the first group is repopulated as members quit their jobs and activities.

London bike sharing

Data obtained from bike sharing systems is particularly suited to perform network analysis. Statistical analysis of the flows of bikes may provide important information regarding the management of the system and could help increase the efficiency of the service. Cycle

Number of edges

Number of edges

Wed Thu Fri Sat Sun Mon Tue Wed q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 5000 10000 15000 20000 25000

Number of edges per day

Number of edges

Wed Thu Fri Sat Sun Mon Tue Wed q q q q q q q q Figure 7: The number of edges in each time frame for the London bikes dataset are shown on the left panel. The right panel shows instead the total number of edges in each day.

hires data can be easily rearranged and visualised as a dynamic network structure, where edges correspond to hires and nodes to stations. Similar approaches have been recently

proposed in a number of works, such as: [START_REF] Guigourès | Discovering patterns in time-varying graphs: a triclustering approach[END_REF], [START_REF] Randriamanamihaga | Clustering the Velib dynamic Origin/Destination flows using a family of Poisson mixture models[END_REF] and Matias et al. (2016). In a similar fashion, here we propose an application of our methodology to a dataset of cycle hires in London.

The data. The cycle hire data for the London bike sharing system is publicly available from Transport for London (2016). We use here the data from Wednesday, 5 June 2013 to Wednesday, 12 June 2013 (included). The discrete time frames correspond to blocks of three hours. For each time frame, we create a network adjacency matrix by transmuting bike trips into directed edges: an edge from station i to station j appears at time frame t if at least one bike is hired at station i during the corresponding three hours, and the same bike is then returned to station j. Bikes may be returned to the same station where they were hired: we consider this information not important for our analysis and hence discard all of the self-edges. The dynamic network so obtained is made of N = 566 nodes and T = 64 time frames. The data exhibits a very strong temporal heterogeneity, mainly due to the day-night cycle and the presence of the weekend days. In Figure 7 the observed number of edges is shown at every time frame. Weekdays typically exhibit two peaks in the activity level, corresponding to the start and the end of office hours.

This suggests that commuters are the main users during working days. By contrast, low

Number of non-empty groups

Number of non-empty groups

Wed Thu Fri Sat Sun Mon Tue Wed q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 8: This plot shows the number of non-empty groups at each time frame for the London bikes dataset. The total number of unique groups found is 43. The level of heterogeneity in the network has a noticeable temporal dynamic, which mimics the evolution of the number of links shown in Figure 7.

activity is observed throughout the weekend, with a single activity peak appearing in the late mornings.

Results. We run our algorithm GreedyIcl once for each of the initialisation methods, with each run taking about six hours long. The optimal clustering found has a total of K = 43 groups.

Heterogeneity. The temporal dynamic of the activity level is also exhibited in Figure 8 by the number of active groups at each time frame: as in the Enron dataset, we use this as a measure for the level of heterogeneity in the network. It appears that, as activity peaks in weekdays, the number of active groups doubles, signaling a heavy increase in heterogeneity within the network. Surprisingly, weekend days do not exhibit a markedly different heterogeneity level with respect to weekdays, even though the overall activity level is lower. This suggests that fewer nodes become active in the weekend, but their connection patterns are not particularly different than those seen in weekdays.

Characterisation of the groups. The groups can be roughly divided in two categories: the first 20 groups have a large aggregated size, mostly high stability, low connection probabilities and a good balance between expected out-degree and in-degree. The

Size of group 2

Group size Wed Thu Fri Sat Sun Mon Tue Wed q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 20 40 60 80 100

Size of group 6

Group size

Wed Thu Fri Sat Sun Mon Tue Wed q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 20 40 60 80 100

Size of group 12

Group size Wed Thu Fri Sat Sun Mon Tue Wed q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 20 40 60 80 100

Size of group 18

Group size Wed Thu Fri Sat Sun Mon Tue Wed q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 20 40 60 80 100

Size of group 21

Group size

Wed Thu Fri Sat Sun Mon Tue Wed q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 20 40 60 80 100

Size of group 33

Group size

Wed Thu Fri Sat Sun Mon Tue Wed q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 9: London bikes data. Temporal dynamics for the sizes of some of the groups found.

The three groups on the upper row belong to the first category (groups containing mostly inactive nodes): group 2 swells consistently every night; group 6 is stably present throughout; whereas group 12 is almost emptied in every weekday at the start of the office hours. The groups on the lower row contain instead very active stations: group 18 is emptied every night; group 21 is also emptied throughout the weekends; whereas group 33 activates every weekday at the start of office hours.

rest of the groups have smaller sizes, high instability, higher connection probabilities and in some cases very unbalanced out-degrees and in-degrees.

Figure 9 shows the temporal evolution of the size of a selection of groups. It appears that the trends in the activity level have a huge effect on the migrations of stations between groups. As nodes become more active, they leave the large groups (first category)

and move to smaller groups, which are better suited to capture the details of their new connection profile in the network. During the high-activity regime these nodes may end up visiting several groups based on the evolution of their connectivity patterns. Nonetheless a good portion of nodes is not affected by the variation of activity level and simply remains stably in the larger and inactive groups.

These arguments are very well supported by the plots in Figure 10. It appears that the stability of nodes is very heterogeneous in that there are many very stable nodes but also some very unstable ones (left panel). At the same time, the plot on the right hand

Stability of nodes

Number of swaps Frequencies q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1000 2000 3000 4000 0.0 0.2 0.4 0.6 0.8 1.0

Proportion of nodes changing allocation

Aggregated group size Proportion q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q

Figure 10: Stability level in the London bikes dataset. The left panel shows that although a good portion of nodes are very stable, there are some stations changing their allocations in more than 45 out of 64 time frames. The plot on the right hand side studies instead the stability of groups.

It appears that larger groups are very stable, and viceversa smaller groups are not. The size of the circles is proportional to the size of the group aggregated over time.

side reinforces the idea that larger groups are also very stable whereas smaller groups exhibits high instability.

As shown in Figure 11 the stability does not seem to be particularly related to the geographical position of stations. Furthermore, the number of groups containing stable stations (nodes that never change allocation) is relatively high.

One additional feature that distinguishes the two categories of groups found can be observed in Figure 12. Here, for every group, the expected out-degrees versus the expected in-degrees are plotted. While stable groups (first category) tend to exhibit a good balance between the two degrees, smaller and unstable groups can also have very unbalanced degrees. This means, evidently, that whenever a node joins one of this unbalanced groups, it may send many more edges than what it receives (or viceversa). In this bike sharing context, the corresponding stations may require special attention due to temporary excess of hiring demands (or excess of arrivals).

Conclusions

The statistical analysis of dynamic networks is particularly challenging, both from a modelling point of view (due to the temporal dynamics) and from an estimation point

Stations not changing allocation

Longitude Latitude q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -0. Most unstable stations (more than 35 swaps) Longitude Latitude q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 11: On the left panel, the spatial distribution of the stations not changing allocations is shown for the London bikes dataset. The colour of the circles correspond to the group they are allocated to at every time frame. On the right panel, instead, the stations performing more than 35 swaps are shown. In both plots the size corresponds to a combination of the out-degree and in-degree of the nodes.

Expected average degrees of groups

Average in-degree Average out-degree q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 12: Expected out-degree versus expected in-degree for the groups found in the London bikes dataset. The size of circles is proportional to the size of groups aggregated over time. All of the large groups of inactive nodes have balanced degrees, in that they tend to send and receive the same number of edges. Smaller groups, instead, can exhibit very unbalanced degrees, either in favour of the out-degree or in-degree.

of view, due to the inherent computational difficulties. In this paper we have focused on an extension of the stochastic block model to dynamic networks, where the temporal evolution of the nodal information is captured by a Markovian process. Our formulation allows one to integrate out all of the model parameters from both the likelihood and the prior, thereby obtaining an analytical formula for the marginal posterior of the allocation variables. In a model-based clustering context, such a marginal posterior is equivalent to the exact integrated completed likelihood, which is widely used as an optimality criterion for partitions.

Taking advantage of these results, we have proposed a greedy algorithm to estimate the optimal clustering of the nodes. Our algorithm resembles other tools proposed in the recent literature, and scales particularly well with the size of the data. However, only convergence to a local optimum is guaranteed, hence several restarts and a careful initialisation are particularly useful in order to find the global optimum.

One appealing feature of our approach is that model-choice is carried out automatically, since the number of groups can be deduced at each time frame from the optimal clustering solution.

Through a simulation study, we have validated both our optimality criterion and the greedy algorithm used to estimate the optimal partition. Also, we have compared our method to one based on a variational Expectation-Maximisation algorithm, showing that with careful initialisation our approach achieves better results when the underlying number of groups is unknown.

We have applied our methodology to the Enron email dataset, and found 17 underlying groups. Each of the groups found appears to have a specific role within the network and the original status of the members seem to be particularly related with the nodes' allocations. Also, we have analysed some relevant summaries from our optimal clustering underlining the dynamic evolution of the heterogeneity within the network and of the overall activity level.

We have also applied our methodology to a London Bikes dataset, to study the flows and connections between bike stations. Our method has returned an optimal clustering made of 43 groups. The analysis of such a complex structure has been particularly challenging but we have managed to extract some interesting information regarding the dynamics captured by the model.

We note that, while the practice of integrating out the likelihood parameters has been exploited in a number of recent papers, the collapsing of the transition probabilities in the Markov process is a rather new original technique. Not only this creates an ideal setup for our method, but it may also be generalised and exploited in arbitrary discrete hidden Markov models and their extensions.

Additionally, our work may be extended in a number of ways. We have used our algorithm only on binary networks, but this approach generalises to networks with other edge types. Incomplete weighted graphs can be easily handled using a Bernoulli presenceabsence indicator, and hence using the same framework based on the collapsing of the parameters.

The initialisation plays an important role in our method, yet finding a good starting partition is a great challenge. In fact, there is a lack of scalable methods that can handle complex dynamic objects such as networks. Furthermore, clustering cannot be attempted for each time frame independently since this result in label-switching problems. Here, we have proposed two new intialisation methods which may in some cases improve the results.

Extensions to the supervised classification case are also straightforward: if some of the allocations are known, these need not be updated during the greedy optimisation. From a Bayesian perspective, the optimal clustering solution obtained will then maximise the posterior predictive distribution, rather than the posterior. This strategy would be useful when dealing with nodes that join or leave the study dynamically, since, as suggested by [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF], a fixed group of inactive nodes may be created.

More generally, the optimisation of the exact Integrated Completed Likelihood is a particularly challenging task, due to the discrete search space and the pronounced multimodality of such objective function. The GreedyIcl algorithm introduced in this paper was shown to perform well, although future extensions may improve upon our work in terms of efficiency using alternative optimisation routines.

Figure 1 :

 1 Figure 1: Graphical model for the Markovian stochastic block model described.

  Matias and Miele 2016, as shown in Figure3.7 Enron datasetThe data. The Enron Corporation filed for bankruptcy in late 2001, leading to an unprecedented scandal and to dire consequences for the US stock market. The data we use in this paper consists of all the emails exchanged from January 2000 to March 2002 between the Enron members. This data was originally made public, and posted to the web, by the Federal Energy Regulatory Commission during its investigation.

  frame t corresponding to a one month period. The number of nodes is N = 148 and the number of time frames is T = 27. The number of edges observed at each time is shown in Figure4. Along with the email data, the status of each member within the company is known, and is one of the following: CEO, Director, Employee, In House Lawyer, Manager, Managing Director, N/A, President, Trader, Vice President. To simplify the exposition of the results, we gather every status other than Employee and N/A into a single class named Manager, so that only three status classes are considered.

Figure 4 :Figure 5 :

 45 Figure 4: Number of edges at each time frame for the Enron email dataset. The peak in October 2001 corresponds to the disclosure of bankruptcy.

Table 1 :

 1 A list of edge types that can be accounted for using the methodology proposed. For each case the corresponding Bayesian hierarchical structure is shown.

	Edge type	Likelihood	Conjugate prior
	Binary	Bernoulli	Beta
	Categorical	Multinomial	Dirichlet
	Counts (positive integers)	Poisson	Gamma
	Positive weights	Gamma	Gamma
	Truncated counts or proportions	Binomial	Beta
	Heavy tailed positive weights	Pareto	Gamma
	Real numbers	Normal	Normal-Gamma
	Real numbers with covariates	Multivariate Normal Normal-Gamma
	and N gh counts the total number of possible edges from group g to h:

Table 3 :

 3 Rounded average degrees for the Enron dataset, for each of the 17 clusters found.

Table 4 :

 4 Cluster counts separated by status for the Enron dataset. For this table, the same node at two different time frames is considered as two separate entities, hence the sum of all entries is T N .

		1	2	3	4	5 6 7 8 9 10 11 12 13 14 15 16 17
	Employee 576 90 36 45 92 50 61 22 7 46 1 18 0 17 18 0 1
	Manager 1027 231 191 102 34 83 53 16 89 40 87 45 41 6 5 0 2
	N/A	475 57 27 50 59 40 12 84 2 9 4 5 2 16 5 16 1

Table 5 :

 5 Group sizes at each time step for the Enron dataset.

Regarding the likelihood structure, a number of conjugate models can be considered, encompassing most edges' types. Table1provides a list of possible scenarios that may be of interest. In this paper, we focus on binary edges, and hence specify an independent Beta prior on the likelihood parameters Θ = {θ gh : g ∈ K, h ∈ K}. In such a case the compound distribution for X ("marginal likelihood") is:p (X |Z, a, b) = g∈K h∈K

Acknowledgements

The Insight Centre for Data Analytics is supported by Science Foundation Ireland under Grant Number SFI/12/RC/2289. Riccardo Rastelli and Nial Friel's research was supported by a Science Foundation Ireland grant: 12/IP/1424.


J Gain

H Long

G Cordonnier

M.-P Cani

EcoBrush: Interactive Control of Visually Consistent Large-Scale Ecosystems

Keywords: Descriptors (according to ACM CCS): I, 3, 7 [Computer Graphics]: Three-dimensional graphics and realism-

Figure 1: Terrain conditions (a), such as temperature, soil moisture, and sunlight exposure, are used to index a database of plant distributions (b) and synthesise an initial complete ecosystem (c), which can then be modified with semantic brushes, to adjust age, density and variability (d).

Introduction

Natural landscapes serve an important function in computer graphics and virtual environments. Indeed, in applications ranging from computer games and movies to landscape design and VR-based training, virtual landscapes are frequently a dominant visual element. One key aspect, particularly in larger scenes, is the depic-tion of plant ecosystems, from the modelling of individual plant specimens, to their collective placement, and final rendering. While models of individual plants [START_REF] Deussen | Digital design of nature: computer generated plants and organics[END_REF], and rendering of large-scale ecosystems [START_REF] Colditz | Realistic and interactive visualization of high-density plant ecosystems[END_REF] have received significant attention, the ecoplacement problem, which extends beyond plant position to include attributes such as species, age, height, vigour and canopy extent, is less well studied.

Broadly speaking, there are two approaches to eco-placement: bottom-up simulation of individual plant interactions [BSK * 15] or top-down randomised synthesis from statistical distributions [EVC * 15]. Simulation has the advantage that it explicitly encodes plant dynamics and supports known emergent ecosystem behaviour, such as succession, where slower growing but more suited plants gradually dominate, and self-thinning, where the overall number of plants is reduced in favour of fewer mature specimens. Unfortunately, it can take hours or even days to generate large-scale ecosystems, even with appropriate n-body parallelisation [NHP07]. Synthesis, on the other hand, uses random placement guided by statistical properties, and is capable of interactively populating large terrains with hundreds of thousands of plants. However, the difficulty lies in encoding ecosystem properties in the distributions so that the results are visually consistent with real landscapes.

Our key strategy is to complete the loop between simulation and synthesis, enabling many small-scale but high-fidelity simulations to be fed as inputs to distribution analysis and later synthesis. Figure 1 shows the typical workflow and key components of our system. The user supplies a base terrain and some non-biological (abiotic) input parameters, such as monthly rainfall, from which monthly terrain conditions for temperature, soil moisture, sunlight exposure and slope are derived. Off-line a database is pre-populated with histogram-based disk-distributions, which capture the statistical interactions between plant canopies and are derived from sandbox simulations of a particular biome. These samples capture the biological (biotic) competition for resources between individual plants. Based on the abiotic maps, distributions are drawn from the database and synthesised across the terrain to rapidly populate an initial ecosystem. A user can then overpaint the ecosystem to locally modify age, variability and plant density, while respecting the underlying terrain conditions.

Our technical contributions include:

1. Extensions to the state-of-the-art in both ecosystem simulation and synthesis. On the simulation side, we derive abiotic maps automatically rather than expecting artists to paint complex geomorphological attributes manually, and enhance plant simulation by separating between moisture (root) and sunlight (canopy) competition, as well as by incorporating gradiations of canopy density. On the synthesis side, we depart from the use of point processes, to include disk-based statistics. This is vital to reproducing the canopy-based spacing of plants and for placing shade-tolerant species under cover. 2. Two alternative approaches for using sandbox ecosystem simulations: First, uniform conditions can be sampled from the range of the input biome. Alternatively, terrain conditions can be clustered and cluster means used as simulation inputs, for a manageable set of seed distributions. This second solution enables us to generate very large ecosystems (10 × 10 km 2 ) without the need for a large distribution database. 3. Consideration of the balance between artistic freedom and realworld consistency. Since distributions are indexed by terrain conditions, we can offer freedom to the user -such as overpainting a region with another distribution or arbitrarily tuning plant density -while providing a suitability map that indicates consistency with terrain conditions. If desired, the user can then use a "healing brush" to progressively shift the new distribution towards the closest consistent state.

These key contributions enable the efficient creation and interactive editing of large-scale landscapes that accord with known emergent properties from plant ecology.

Related Work

There are many challenges to effectively portraying natural environments in computer graphics, including modelling the fine-scale geometry of individual plants and rendering the resulting scenes, which exhibit both large scale and high depth complexity. If we scope our problem as populating a terrain with a collection of plants, i.e. assigning attribute values to individual specimens, such as position, species, height, canopy extent, and vigour, then there are two broad approaches: bottom-up simulation and topdown statistical synthesis.

Simulation takes the Lagrangian approach of treating plants as particles. Using either simple radial interactions [DHL * 98, AD05], L-systems [LP * 02] or more complex agent-based models [Ch'13, BSK * 15], plants are seeded at a location, grow in size asymmetrically [AD05] under competition for resources, propagate on reaching maturity, and die due to lack of resources or old age. The fidelity of these simulations varies depending on which abiotic features, such as temperature, sunlight, moisture, wind, and slope, are considered and on the biotic complexity of individual plant models and their local interactions. For instance, ecosystem disturbance by humans, grazing animals, and fire is, with a few exceptions [BE03, Ch'13], otherwise neglected. In the simulation literature, Ch'ng's models [Ch'09, Ch'13] are the most complete, incorporating layered soils, various forms of seeding and seasonal differences in evergreen and deciduous species. Nevertheless, the botanical literature [FPR * 96, SIK07, SHG * 08] goes further to model in detail cycles such as evapotranspiration, photosynthesis and carbon dynamics. With some exceptions [START_REF] Sato | Seib-dgvm: A new dynamic global vegetation model using a spatially explicit individual-based approach[END_REF] this is in the service of Dynamic Global Vegetation Models, which are grid-based Eulerian simulations of climate impact on plant distributions. These typically provide dry-mass proportions for species in cells ranging from 30m to several kilometres on a side, and do not directly enable plant placement.

There is, however, a tension between simulation fidelity and performance. Complex models, even with approximate nearestneighbour acceleration, do not scale well to large, dense environ-ments and may take minutes or even hours to execute. A viable approach is to tailor models so that they consider abiotic and biotic factors only to the extent that emergent properties such as succession, self-thinning, clustering and shading are correctly exhibited. We adopt this strategy in precomputing ecosystems over small areas.

Statistical synthesis, on the other hand, generates plant positions randomly to fit a given distribution. Solutions range from straightforward half-toning of a density image [DHL * 98] and dart throwing [ACV * 14], which exhibit little more than basic plant separation, to Wang tiles, with an underlying Field-of-Neighbourhood plant distribution model [AD06], kernel-based deformation of probability density maps [LP * 02] and derived cross-species distribution histograms [EVC * 15], which can capture complex interspecies and environmental interactions. These approaches fall under the umbrella of the probability theory of statistical point processes [START_REF] Diggle | Statistical analysis of spatial and spatio-temporal point patterns[END_REF]. In fact, plant ecologists have long drawn on this theory [LIB * 09] to extract summary statistics for field-based sampling. It has also seen recent application in computer graphics [OG12, EVC * 15] for the analysis and synthesis of point distributions. While building on this vein of research, we extend it by explicitly incorporating the extent of plant canopies into distributions. This enables us to capture effects such as the coverage of shadeloving plants.

The issue of user interaction in ecosystem modelling has received less attention than is warranted. The ideal is direct, interactive and intuitive control over the final state of the ecosystem in a way that respects both the users' intentions and physical plausibility, perhaps through a sculpting, painting or sketching interface. Unfortunately, the interfaces of ecosystem simulators [DHL * 98, Ch'13] commonly take the form of painting input image maps (e.g., water, soil type, temperature). Users without a deep understanding of geomorphology and plant dynamics have difficulty designing final lengthy simulation outcomes from the poorly understood inputs of these systems. Bradbury et al. [BSK * 15] partially address this problem by allowing userdefined regions of a simulation to be paused, edited and then restarted to run at different simulation rates. WorldBrush [EVC * 15] provides more direct and immediate interaction through a palette of distributions that can be painted directly onto a landscape, with regions then modified through copy and paste, gradient interpolation, and seam-carved stretching. Unfortunately, WorldBrush neither enforces nor guides realism and instead relies on the user to maintain the link between abiotic conditions and selected distributions.

Overview

Our goal is to allow users to interactively design large-scale ecosystems that accord with observed landscapes and botanical principles. Regrettably, interactive simulation of ecosystems larger than a few hundred metres on a side is generally infeasible, because, all other things being equal, the number of plants is proportional to the simulation area. Instead, we use precomputed distributions derived from small sandbox simulations. Analyzing simulation results enables us to derive an ecosystem with the same statistical properties at least two orders of magnitude faster, enabling subsequent interactive editing while preserving consistency. The user provides a heightfield and other abiotic input parameters to a simulation process, which generates abiotic terrain maps for temperature, soil moisture, slope, and sunlight exposure. A clustering process then groups terrain regions according to similar abiotic conditions. As a separate preprocess, the range of possible abiotic conditions are sampled and small scale sandbox simulations are executed and then analysed to create a database of distributions. The abiotic maps, clusters, and distribution database are used to synthesise an ecoystem, which can then be modified by the user using semantic painting and resynthesised interactively.

More precisely, our system automatically populates the input terrain with an initial ecosystem, where plants match local abiotic conditions computed on the terrain. This ecosystem is generated from a meaningful set of sandbox simulations, analyzed and stored in a distribution database: the right plant distribution is synthesized on each region of the terrain, depending on the local abiotic conditions (see Figure 2). Users are then able to paint with a set of modifying brushes that can increase or decrease ecosystem age, variability, and density, even targeting individual species types. They can also select and paint specific distributions by treating the distribution database as a palette, while being guided by the system to maintain consistency.

To initially populate an ecosystem, the user supplies some simple inputs: a terrain heightfield, latitude, soil type, yearly temperature extremes, and average monthly rainfall. Rather than requiring users to paint abiotic maps, which is difficult to do without domain expertise, we use these inputs to simulate terrain conditions and automatically generate the maps. They include monthly abiotic condition maps for temperature (as affected by altitude), daily sunlight exposure (influenced by latitude and terrain self-shadowing), and soil moisture (a combination of rainfall, slope, soil absorption and runoff). An important insight of our approach is to cluster regions with similar conditions, enabling us to factorize the need for ecosystem simulations. This is done by collectively applying kmeans clustering to the abiotic maps. This first stage of the pipeline is detailed in Section 4.

The clusters of abiotic conditions can either be used to uniformly sample conditions for ecosystem simulations, such as temperature, sunlight exposure, and soil moisture, from the range of the corresponding biome; or alternatively the cluster means can be directly used to trigger a minimal set of such simulations. In both cases, we use small scale (100×100m 2 ) sandbox simulations (See Section 5). Also supplied to these simulations are viability parameters for a set of plant archetypes (or Plant Functional Types (PFT's) in plant ecosystem parlance), representing collections of similar species, such as "Tropical Broad-leaved Evergreens" or "Boreal Needleleaved Summergreens". This approach avoids having to exhaustively determine parameters for hundreds of species [FPR * 96].

Next, plant attributes (type, position and canopy extent) from simulation results are analysed to derive a new type of statistical distribution, where point processes are extended to distributions of overlapping disks representing plant canopies (See Section 6). Stored as pairwise interaction histograms relating expected distances between and within PFT categories, these distributions are placed in a database ready for subsequent synthesis to reproduce the simulation characteristics on a larger scale. Populating the database is kept manageable by restricting the simulation dimensions and executing them as a pre-process.

The ecosystem we generate from these distributions (where individual species are finally mapped from PFT's using a statistical post-process) can be used for rapidly and automatically populating large terrains without the need for any user intervention. We also believe this initial ecosystem with typical plant coverage is more helpful as a starting point for users than a blank canvas. We then provide a set of high-level painting tools enabling users to interactively and locally modify distributions and hence synthesised plant populations. The challenge is to design tools that provide freedom to the user while maintaining consistency in terms of abiotic constraints. This last interactive editing stage is described in Section 7.

Simulation and Clustering of Abiotic Terrain Conditions

In deriving abiotic conditions our aim is to require only a limited set of comprehensible user inputs, while still capturing significant terrain effects that impact ecosystems, including shadowing in defiles, increased runoff on slopes, water table access near rivers, temperature lapse with altitude, and the effects of latitude on growth cycles. Nevertheless, in the interests of efficiency our approach remains approximate and many secondary factors are not considered. In particular, our simulations are forward directed and do not consider feedback cycles, such as evaporation due to canopy interception and soil nutrients from decomposing litterfall.

The product of our simulation of abiotic conditions is a set of monthly abiotic maps over the terrain (see Figure 3) with averaged daily values for maximum temperature, sunlight exposure, and soil moisture, derived as follows:

Temperature -Our consideration of temperature is straightforward: the user provides maximum daily temperatures (averaged over a month) for the middle of winter and summer, at the lowest point on the terrain. These temperatures are linearly interpolated to the remaining months and reduced with altitude (at a userconfigurable lapse rate, which defaults to 6.5 • C per 1000m) to provide average daily maximum temperatures for each month. Sunlight Exposure -From a latitude and compass direction provided for the terrain, and factoring in the earth's axial tilt, it is possible to derive the sun's trajectory at the middle of a given month [START_REF] Sato | Seib-dgvm: A new dynamic global vegetation model using a spatially explicit individual-based approach[END_REF]. We cast rays from the sun, sampled at hourly intervals, to each heightfield position. If a ray is intersection free then the terrain position is not in shadow for that hour. This allows us to derive a map of the daily hours of direct sunlight on a monthly basis.

Soil Moisture -In reality the water content of soil is governed by a complex combination of precipitation, evaporation, runoff, seepage, soil type and water uptake by plants, among other factors. We approximate this with an equilibrium model that dynamically adjusts a single plant uptake value kp so that over multiple yearly cycles with the same precipitation the soil moisture achieves a rough balance. Also, instead of modelling soil type explicitly, the user configures a single maximal monthly soil absorption value, ka, and overall reservoir capacity, kc. Slope acts as a proxy for soil type by linearly weighting the absorption to represent the transition from loamy soils in valleys to rock on cliff faces. Any rainfall above this value contributes directly to runoff. Thus, we model the soil moisture M i,m for a given cell i and month m as:

M i,m = min(kc, M i,m-1 + min(Pm, A i ) -kp),
where Pm is the monthly precipitation provided by the user, and A i is the per cell absorption threshold. Any excess absorbed moisture, S i,m , above the soil capacity is spread as seepage over the next three months. Free surface water is then defined as a combination of runoff and seepage, as follows: From this runoff and seepage, a network of waterways can be established. For our purposes, plant growth is prevented in areas of open water but enhanced along riverbanks due to a more accessible water table. We employ the drainage model of Cordonnier et al. [CBC * 16], which provides water volume (F i,m ) taking into account excess initial moisture (F i,m ), inflow from neighbouring cells and outflow to the lowest neighbour. Any cell above a threshold (F i,m > ko) is marked as open water. These waterways are then expanded to neighbouring cells j within distance d i j < 0.5 W i using Tucker et al.'s [START_REF] Tucker | Modelling landscape evolution[END_REF] rough river-width formula: W i = 10 F 0.5 i (with F i in m 3 s -1 ). A second shell of cells within d i j < 0.5 km W i = Dm has its soil moisture (M j ) increased by a set amount ( fr) that is tailed off in proportion to relative slope and distance from the river, so as to represent water infiltration along riverbanks. For each nonriver cell j:

F i,m = max(0, Pm -A i ) + 3 ∑ j=1 1 3 S i,m-j
M j,m = M j,m + max(0, max i∈Nj Dm -d i j * (1 + kss i j ) Dm fr),
where s i j is the slope, N j is a set of cells in the neighbourhood of j, and km and ks are weighting constants for slope and distance, used to adjust the extent of riverbanks.

The final step in this pipeline is to identify clusters over the parameter space of abiotic conditions. These serve as seeds for an initial ecosystem and, in the absence of a distribution database, can supply a usefully constrained set of simulation parameters. Despite the large number of dimensions (12 each for temperature, moisture, illumination, 1 for slope, and 1 for age = 38), conditions are correlated and continuous to the extent that clusters are relatively contiguous and simple k-means clustering [Jai10] suffices. Typically, a choice of 6 -10 clusters results in cluster means that do not diverge too markedly from the underlying conditions (as evident in Figure 3). In order to supply plausible plant distributions to our analysis and synthesis procedure, a series of simulations is undertaken in small-scale 100 × 100m 2 regions with uniform input conditions. Conceptually, a sandbox simulation proceeds as follows: In each month, the vigour, V (p) ∈ [-1, 1], of every plant, p, is calculated. If positive, this vigour value weights the plant's growth rate; if negative, it provides a probability of mortality, as does a plant reaching the age of senescence. Vigour is a function of the Plant Functional Type (PFT) attributes of a specimen, the abiotic conditions in the plant's neighbourhood and biotic competition with other plants that impinge on a specimen's root or canopy extent. Note that unlike previous work we consider these separately.

Sandbox Simulation

If a plant has reached maturity and has sufficient vigour it will annually seed new plants randomly out to a radius determined by the PFT-specific propagation (cones, fruit, etc.). Each Plant Functional Type has a viability function that models response to temperature (B T ), sunlight (B L ), soil moisture (B M ), and slope (B S ), based in part on the corresponding abiotic maps. These are piecewise linear functions (see Figure 5) that range from -1 at the climatic extremes to a plateau of 1 over the optimal inner range. A plant's vigour is then the minimum over all viability function values: V (p) = min(B T (p), B L (p), B M (p), B S (p)). Taking the minimum in this way is motivated by the consequences of a limiting resource. It matters little, for example, that a plant has abundant water if it is starved of light.

Viability

Condition

Slope and temperature values are submitted directly to the corresponding viability functions (B S , B T ), but sunlight and moisture are affected by inter-plant canopy and root interaction (see Figure 4). From a plant's height and type, the circular extent of root and canopy can be derived (such ratios are known as allometries and widely used in plant biology simulation [SIK07, SHG * 08]). The intersection of these extents is integrated over a fine-scale grid (we use 20 × 20cm 2 cells). For sunlight, the tallest plant in a given cell receives the full sunlight and transmits a proportion, based on its type-specific alpha, to the undercanopy of shorter plants. For moisture, in cases where the minimum needs ∑ i E min of all incident plants i cannot be met, assignment is weighted according to relative root radius. Our motivation is that larger plants tend to dominate through deeper roots and better access to the water table. Finally, a plant's sunlight and moisture values are averaged over all cells within the canopy and root radii, respectively. This allows for partial shading and also seamlessly scales the greater resource requirements of larger plants.

Grass is a special case: a count of grass clumps over a reasonable expanse will likely number in the billions, and synthesising them on the same individual basis as plants is thus clearly infeasible. Rather, as a final pass after plant placement, we subsample the terrain grid and establish grass height per cell. The growth cycle is bypassed and vigour values, based on averaged yearly abiotic conditions, are translated directly into grass height. In terms of biotic factors, we ignore moisture competition on the basis that grass accesses only the topmost layer of soil, but consider sunlight competition by applying a concentric burn of grass heights (moderated by appropriate PFT alpha values) under plant canopies. Because abiotic conditions are sampled more coarsely than grass cells, we also found it necessary to introduce some random variation and low-pass averaging.

A B A A A B B B R A R B A -B > R A + R B R A < A -B ≤ R A + R B R A -R B < A -B ≤ R A R B < A -B ≤ R A + R B (a) (b) (c) (d)
Figure 6: The different pair-wise disk configurations for plant B with respect to A: (a) bands of distances between the disks, (b) less than half shaded, (c) more than half shaded, (d) fully shaded.

Disk-based Distribution Analysis and Synthesis

One of the main features of our method is to use the ecosystem simulations described above to derive statistical distributions, which can be used to synthesize plants over a much larger terrain at interactive rates. In contrast to conventional point processes, we analyze distributions of possibly overlapping circles of different radii, representing canopies instead of mere points. This extension is an original contribution to point-based analysis and synthesis methods.

Our point of departure is the WorldBrush system of Emilien et al. [EVC * 15]. WorldBrush is based on categories C i of scene elements (e.g., trees, plants, buildings, and roads) that are ranked to establish both dependencies and an order of evaluation. In our context, each plant functional type is separated by age (seedling, young, mature) into three categories, and all PFT categories are then prioritized according to their relative dominance over resources.

At its core, WorldBrush is built on pair-wise distance histograms between and within categories. Given categories Cx and Cy where x has higher priority than y, these encode the probability of a plant in Cx falling within a certain distance band of a plant from Cy. Such discretised pair correlation functions have proven capable of capturing irregularity, such as clustering and clumping, far better than noise functions [OG12]. Analysis involves counting the incidence of elements of Cx within concentric annular shells of points in Cy and normalizing over the total number of elements (nx, ny) and the area of the shells to generate a probability distribution. A population of points can then be synthesised by various forms of stochastic search (dart throwing, Metropolis-Hastings) depending on whether the correlation functions are used to score a given distribution or its iterative update.

As mentioned, a pure point process is unable to adequately capture the nuance of plant interaction. Simply drawing a plant's size from a distribution irrespective of position can lead to implausible configurations, such as canopies that intersect the trunks of neighbouring plants or shade-seeking plants that are exposed to direct sunlight. It is thus necessary to simultaneously consider both plant position and canopy radius. (This could, of course, be extended to roots, but we have chosen to focus on the canopy due to its obvious visibility.) This problem can be classified as a continuous marked point pro-cess [START_REF] Diggle | Statistical analysis of spatial and spatio-temporal point patterns[END_REF], in that additional attributes (marks) are associated with positions (points). An obvious approach would be to shrink or grow a plant's radius based on its local neighbourhood, but this necessitates expensive iterative adjustment spanning all categories that quickly devolves into a form of simulation. Another option would be to divide plant types into subcategories by canopy radius, significantly extending our existing age-based categorisation. In a similar vein, we could develop 2D distribution histograms with distance on one axis and radius on the other. Alas, these are costly solutions: the dependence between categories is O(N 2 ), while introducing an extra histogram axis will increase distributions by a multiple of the number of canopy bins.

Instead, our strategy is to measure distances between disks rather than points. By jointly considering the positions of two plants (B relative to A), the distance between them (d = ||A -B||), and their respective radii (R A and R B ), we define a new class of histogram bins (see Figure 6) that encodes the degree of canopy overlap: lessthan-half shaded, more-than-half shaded and fully shaded. For the non-overlapping case (Figure 6(a)) we populate the distribution histograms according to the closest approach between disks (replacing d with d -R A -R B in such calculations). Histogram normalisation also needs to be suitably adjusted to account for R B in both overlapping and non-overlapping cases. This has the virtues of conceptual simplicity and computational efficiency. The only caveat is that it introduces a form of non-uniform spatial distortion in that the areas represented by histogram bins are no longer consistent between distributions, but this seems to have negligible impact in practice.

We also need to adapt the synthesis process. We begin with a dart-throwing process to match expected plant density, while obeying some simple validity rules (d > R A , d > R B and f (X) = 0, where f (X) is the probability density function of the distribution X). This has the benefit of providing rapid initial feedback to the user. Next, eco-placement is refined by iterating between cycles of random position perturbation and radius adjustment of individual plants. A candidate perturbation is accepted with a probability that is proportional to the improvement in the probability density function. In spirit, this is closer to the synthesis of Öztireli and Gross [OG12] than the Metropolis-Hastings sampling of World-Brush [EVC * 15].

(2b) (1a) (1b) (2a)
Figure 7: Semantic brushes can be applied locally to an existing ecosystem to, for example, adjust age (increased from (1a) to (1b)), and density (increased (2a) to (2b)).

Semantic Brushes for Plant Distributions

By coupling simulations and distributions we are now able to automatically generate full-scale ecosystems over varied terrain conditions, but it is important to also provide meaningful interactive controls to artists and designers. The WorldBrush solution has users pick a particular distribution from a palette and paint it directly onto the terrain. In our case, this could spoil the plausibility of plant distributions. In contrast, we provide high-level semantic control through a range of context-sensitive brushes to modify the ecosystem age, variability, and plant density while maintaining consistency (see Figures 7 and9).

(1a) (1b) (2a) (2b)
Figure 8: The fit of an ecosystem (1a) is presented to users as a suitability map (1b), where closeness of fit ranges from green (good) to orange (poor). A semantic healing brush enables users to shift ecosystems closer to the underlying terrain conditions (2a), as reflected by improved suitability (2b).

First a palette of possible distributions needs to be provided to the user. Our sandbox simulation and disk-based distribution analysis pipeline can be interpreted as a mapping (E : T → D) from a space of terrain conditions (T ) to distributions (D). To support full ecosystem variety and semantic brush editing we approximate E by building a database of distributions that sample terrain conditions on a regular grid, g ∈ T . Since a naive sampling would be intractable due to the high-dimension of the parameter space, we make the following approximations: First, we decrease parameter dimensionality from 38 to 5 using a single seed value for moisture, sunlight and temperature, from which monthly patterns are derived prior to simulation. This approximation sacrifices some of the nuance of monthly variation and captures only broad seasonal effects but reduces the number of sandbox samples reasonably required to cover the search space and hence the cost of pre-processing. Second, we limit the sampling range to the climatic bounds of a particular biome. This makes the database biome specific, but is consistent with our selection of biome specialised Plant Functional Types. Third, we use a coarse sampling only, completed using interpolation between samples: given some non-sampled terrain conditions, we derive a distribution by fivefold linear interpolation of the sampling grid, with a set of interpolations for each dimension, using mass transport of distribution histograms [EVC * 15].

In the final ecosystem modelling phase, a design map of conditions over the terrain is drawn from T and used to index the distribution database. Interactive design starts from an ecosystem automatically populated by synthesising the obtained distributions according to the design map (see §4 and Figure 3). A set of semantic brushes of controllable radius can be used to modify the design map within a distance (d ≤ r 0 ) of the brush centre. We support feathering on the edges of the brush by gradiating the weight of any modification between an inner and outer brush radius (r 0 < d ≤ r 1 ).

The various brushes function as follows:

• Age -Since ecosystem age is one of the dimensions of T , an ecosystem can be made older or younger by shifting along this axis while holding other coordinates constant. • Variability -Moving conditions towards or away from the local mean of the design map, enables ecosystem variation to be smoothed or exaggerated. • Density -A common summary statistic of distributions is the ex- pected number of plants. By looking up distributions in the local T -space neighbourhood of a condition from the design map, a T -vector that increases or decreases plant density can be found. This can even be specialised to particular plant types or size categories.

The original abiotic conditions (the source map) can be regarded as a ground truth of sorts. Differences between the design and source maps are necessary for design freedom and can be motivated by ecosystem effects not catered for in our simulations, such as microclimates, local soil types, and fire disturbance. Nevertheless, in the interests of visual consistency it is useful to reflect this divergence back to the user in the form of a suitability map (see Figure 8), which uses colour gradiation based on the normalised L 2 distance between the source and design maps. We also provide a healing brush that shifts overpainted areas of the design map towards the source map.

Results and Discussion

In the interests of visual consistency, it is necessary to establish the validity of our underlying simulation engine. In line with previous work [DHL * 98, LP * 02, BL09, Ch'11] we have checked for expected emergent behaviour, such as self-thinning, where the number of plants is reduced over time as larger specimens crowd out smaller, and succession of a slower-growing but more climaticallysuited species over a fast-growing competitor. In figure 10 we show the introduction of a niche-suited invader species to a Mediterranean biome (patterned on Gritti et al. [START_REF] Gritti | Vulnerability of mediterranean basin ecosystems to climate change and invasion by exotic plant species[END_REF]), which results in dominance over other species.

The variety of achievable landscapes is shown in Figures 111213, all at 5 × 5 km 2 scale. The Plant Functional Type parameters for these biomes are summarised in Table 2 and the input abiotic maps and clusters appear in Figures 1 and3. In our interactive editor each PFT is provided a simplified iconic geometric model (see the key of Table 2) and rendered using OpenGL instancing with per instance transformations for size and placement. Export for final rendering is to Terragen. An important use-case is the generation of terrains in the absence of a database from simulations run on the means of a set of clustered conditions. Figure 12 demonstrates the results of this process, for an Alpine biome with 7 clusters, such as might be found in the Swiss Alps. Of note here is transition from temperate to boreal trees due to changes in mean temperature with increasing altitude.

One shortcoming of our simulations is that they do not cater for disturbance effects. This is particularly evident in Tropical Savannahs, which, without the impact of fires and grazing herds would be relatively dense because of the high rainfall. However, as shown in Figure 13, users can compensate for this with appropriate overpainting of the ecosystem.

The performance of our system was tested on a 3.3 GHz quadcore Intel Haswell, with 7 GB RAM, and an NVIDIA GTX 960 with 4 GB and the results are shown in Table 1. The system performs with reasonable interactive response (< 1s per update) for synthesis with up to 50, 000 plants, which roughly equates to a terrain area of 1 × 1 km 2 , depending on plant density. This does place a limit on brush radius in larger landscapes if interactivity is to be retained. However, it is possible to provide an initial placement estimate (the dart throwing phase) and then run refinements as a background thread (the perturbation phase).

A weakness of our system is the need for manual calibration of plant functional type parameters. Poor calibration tends to result either in homogenization through a lack of niche specialization, or a complete dearth of particular plants. Fortunately, the decision to avoid simulating individual species greatly reduces the search space. Nevertheless, the issue is complicated by the difficulty of mapping botanical data (such as the bioclimatic limits available for Global Vegetation Models) onto our model. This necessitates a cyclic process of adjusting parameters and checking on the emergent properties. Another key parameter is the density of sandbox sampling of terrain conditions. Undersampling leads to particular sub-niches of plant distributions being missed and disparities between the results of our system and a full terrain-spanning simulation, while oversampling unnecessarily increases pre-processing times. This could be improved with a mechanism for scattered sampling of the space of terrain conditions. In general, such parameter tuning is painstaking and could certainly benefit from automation, perhaps through some form of machine learning.

Number of plants

Conclusions

To summarize: we create a database of plant distributions by sampling the space of abiotic terrain conditions (moisture, sunlight, temperature, slope, and ecosystem age) for a given biome and supplying these samples as parameters to small-scale eco-position simulations. Statistical distributions are then derived from the resulting plant placements via disk-based distance analysis. By synthesising distributions drawn from the database a large ecosystem (5 × 5 km 2 in extent with up to a million plants) with varying terrain conditions can be populated in under a minute while retaining the visual consistency of the source simulations. Users can locally modify meaningful high-level attributes of the ecosystem, such as age, density, and variability, using an overpainting interface. Modifier brushes are interactive even on larger terrains, requiring 0.6s to update a 500 × 500m 2 subsection of an ecosystem.

There are several possibilities for future work. Most existing simulations, including our own, do no fully address either ecosystem disturbance (the impact of fires, storms, grazing, logging, disease, and the like) or two-way interaction with the terrain (such as litterfall enriching the soil, roots acting as barriers to erosion, and plant canopies intercepting rainfall). We compensate for this to an extent by allowing users to guide the ecosystems. For instance, fire damage can be approximated by applying an age brush, but this fails to account for fire-prompted seeding and germination or varying topkill and plant mortality. Certainly, there is an opportunity to improve simulation fidelity by borrowing more extensively from plant ecology models [SIK07, SHG * 08].

Another challenge lies in parameter selection and validation. It is very difficult to obtain accurate species information even for relatively simple plant models, such as ours. Our strategy has been to tune plant parameters based on expected emergent behaviour, such as ecological specialisation, succession, and self-thinning, but this suffers from being qualitative rather than quantitative. One alternative would be to validate against physical plot samples undertaken by plant ecologists. The difficulty here lies in finding a mapping from the specifics of this data to more generic parameter inputs and plant functional type outputs. Table 2: A summary of plant functional type parameters for our biome simulations, indicating their general characteristics (form, maximum height, longevity) and relative ranking of response to abiotic factors (on a ten-point scale). The key indicates the shape and colour used for ecosystem visualizations.

Figure 2 :

 2 Figure2: Method overview: The user provides a heightfield and other abiotic input parameters to a simulation process, which generates abiotic terrain maps for temperature, soil moisture, slope, and sunlight exposure. A clustering process then groups terrain regions according to similar abiotic conditions. As a separate preprocess, the range of possible abiotic conditions are sampled and small scale sandbox simulations are executed and then analysed to create a database of distributions. The abiotic maps, clusters, and distribution database are used to synthesise an ecoystem, which can then be modified by the user using semantic painting and resynthesised interactively.

  c 2017 The Author(s) Computer Graphics Forum c 2017 The Eurographics Association and John Wiley & Sons Ltd.

Figure 3 :

 3 Figure 3: Abiotic maps and final clusters: (a) temperature from altitude lapse, (b) sunlight exposure from latitude, affected by terrain shadowing, (c) soil moisture as a consequence of rainfall, absorption and runoff, and (d) clusters formed by grouping combined abiotic conditions. The simulation uses monthly abiotic maps, which are averaged here for illustration purposes.

c

  2017 The Author(s) Computer Graphics Forum c 2017 The Eurographics Association and John Wiley & Sons Ltd.

Figure 4 :

 4 Figure 4: Competition for sunlight and moisture. The circular canopy (green) and root extent (brown) of a plant are discretised on a grid and plants compete in individual cells of overlap (shaded) to derive averaged scores for moisture uptake and sunlight exposure.

Figure 5 :

 5 Figure 5: Piecewise linear functions define PFT-specific viability in terms of climatic extremes [E min , Emax] and an ideal range [I min , Imax] for each abiotic condition.

Figure 9 :

 9 Figure 9: Iterative editing of a Mediterranean canyon landscape. An initial cluster-based plant assignment (a) is modified by the user to introduce grass surrounds and young dense plants along the riverbank (b), followed by painting a varied spread of evergreen shrubs (in yellow and dark green) on the slopes (c), with final Terragen rendering (d).

Figure 11

 11 Figure 11 is a Mediterranean landscape characteristic of the

Figure 10 :

 10 Figure 10: Simulation validation: (a) a Mediterranean biome, (b) after the introduction of a fast-growing invader species (shown in purple), which overwhelms other plants.

Figure 11 :

 11 Figure 11: A Mediterranean biome typical of Southern France, exhibiting a variety of distributions. [scale: 5 × 5 km 2 , input ranges: temperature 4 -29 • C, rainfall 11 -127mm, mostly in winter, ~1.8 million plants]

Figure 12 :

 12 Figure 12: An Alpine biome based on the Swiss Alps. Plant types are elevation dependent: temperate species in the lower, warmer elevations (brown) and boreal species at higher, colder elevations (light green). [scale: 5 × 5 km 2 , input ranges: temperature 2 -23 • C at base, decreasing to -6 -15 • C at the higher elevations, rainfall 53 -161mm, ~1.5 million plants]

Figure 13 :

 13 Figure 13: A tropical Savannah biome inspired by Central Africa. User intervention is required to reduce the initial high plant density (top) caused by high rainfall and thereby reflect the impact of fire and elephant herds (bottom), which are not explicitly modelled in our simulations. [scale: 5 × 5 km 2 , input ranges: temperature 21 -26 • C, 53 -161mm rainfall, ~3.2 million plants reduced to ~1.4 million.]

Table 1 :

 1 Computation cost of synthesis with varying ecosystem size (in number of plants). The peturbation phase of synthesis allows up to 2 adjustments per plant.

		Time (s) Extent (m)
	44,377	0.82	1000
	232,921	5.79	2000
	543,335	18.49	3000
	991,981	49.29	4000
	1,568,848	103.60	5000
	2,245,103	180.15	6000

c 2017 The Author(s) Computer Graphics Forum c 2017 The Eurographics Association and John Wiley & Sons Ltd.

Acknowledgments

We gratefully acknowledge advice and assistance on histogram distributions from Ulysse Vimont, and on Botany from Adam West and Glenn Moncrieff.