Guillaume Aupy
email: guillaume.aupy@inria.fr

Olivier Beaumont
email: olivier.beaumont@inria.fr

Lionel Eyraud-Dubois
email: lionel.eyraud-dubois@inria.fr

What Size Should your Buffers to Disks be?

Burst-Buffers are high throughput, small size intermediate storage systems typically based on SSDs or NVRAM that are designed to be used as a potential buffer between the computing nodes of a supercomputer and its main storage system consisting of hard drives. Their purpose is to absorb the bursts of I/O that many HPC applications experience (for example for saving checkpoints or data from intermediate results). In this paper, we propose a probabilistic model for evaluating the performance of Burst-Buffers. From a model of application and a data management strategy, we build a Markov-chain-based model of the system, that allows us to quickly answer issues about dimensioning of the system: for a given set of applications, and for a given Burst-Buffer size and bandwidth, how often does the buffer overflow? We also provide extensive simulation results to validate our modeling approach.

I. INTRODUCTION

Solving the bottleneck of I/O is a major point in current HPC systems. This point is especially striking when observing their recent evolution. For instance, when Los Alamos National Laboratory moved from Cielo to Trinity, the peak performance moved from 1.4 Petaflops to 40 Petaflops (×28) while the I/O bandwidth moved to 160 GB/s to 1.45TB/s (only ×9) [START_REF]The trinity project[END_REF]. The same kind of results can be observed at Argonne National Laboratory when moving from Intrepid (0.6 PF, 88 GB/s) to Mira (10PF, 240 GB/s) and to Aurora (expected 180PF and 1TB/s) [START_REF]I/o at argonne national laboratory[END_REF]. Indeed, the main storage technology is still based on hard drive, that have shown a better capacity to scale up in terms of storage capacity than speed.

On the other hand, the usage of HPC systems makes I/O more and more important. First, in the framework on the convergence between HPC and BigData [START_REF] Reed | Exascale computing and big data[END_REF], HPC systems are more and more used to run BigData applications, that require much more I/O bandwidth than traditional HPC applications. The main characteristic of BigData workload is that they are dominated by read operations. Second, the MTBF (Mean Time Between Failures) of HPC systems is decreasing [START_REF] Cappello | Toward exascale resilience: 2014 update[END_REF], [START_REF] Bosilca | Unified model for assessing checkpointing protocols at extreme-scale[END_REF] and Checkpoint/Restart strategies are used to ensure the reliability of computations over a failure prone system. Contrarily to BigData applications, Checkpoint/Restart strategies consume a lot of I/O bandwidth for storing checkpoints, and are dominated by write operations. An important feature of checkpoints is that in general, each new checkpoint for an application can erase the previous one, so that not all checkpoints have to be ultimately saved to the slow disk. Third, HPC applications themselves consume a lot of I/O bandwidth (see Section II-C). They typically follow a quasiperiodic pattern, with an alternance of compute and I/O phases, that cannot be overlapped. As in the case of checkpointing, such HPC applications mostly write data to the disk, and due to their bursty nature, they are known to generate interference and idle times at the I/O level [START_REF] Yildiz | On the root causes of cross-application i/o interference in hpc storage systems[END_REF] when several applications attempt to store their results simultaneously.

In order to cope with the limited I/O bandwidth of HPC system, Burst-Buffers have emerged as promising solution [START_REF] Ddn Storage | BURST BUFFER & BEYOND, I/O & Application Acceleration Technology[END_REF], [START_REF] Henseler | Architecture and design of cray datawarp[END_REF], [START_REF] Daley | Performance characterization of scientific workflows for the optimal use of burst buffers[END_REF]. On the technological side, the use of NVRAM or SSDs makes it possible to achieve much higher bandwidth than hard disks. Due to their high cost and limited capacity, Burst-Buffers are not expected to replace hard-drives, but rather act like a potential intermediate storage layer between the computing nodes and the hard drive storage. Such a layer can be used both to increase data locality and to cope with I/O bursts, using their higher I/O bandwidth to avoid interferences and slowdowns.

There is still no clear consensus on Burst-Buffer architecture (see Section II-A) to know whether they should be centralized or distributed over the platform, and whether they should act as a cache between the computing nodes and the storage system or if they could be bypassed by the applications. In this paper, we consider the simplest model where the Burst-Buffer acts as a potential intermediate centralized layer, with a higher I/O bandwidth but a smaller capacity than the hard disk storage system.

Our goal is to propose a probabilistic model for applications, that is amenable to theoretical analysis based on Markov chains and that provides an estimation of the I/O contention as a function of the size of the Burst-Buffer at a very reasonable cost. The rest of this paper is organized as follows.

In Section II, we present the related work on Burst-Buffer architecture, bandwidth allocation and HPC applications models. Then, we present in Section III the probabilistic model that will be used throughout the paper. In Section IV, we propose two methods to evaluate I/O bandwidth overflows respectively based on the use of (weighted) Chernoff bounds and on an ad-hoc algorithm. Then, we prove in Section V that a simple Markov chain can be used to evaluate the idle time induced by Burst-Buffer overflow and we provide simulation results to study the important parameters of the Burst-Buffers. In Section VI, we validate the probabilistic model introduced in Section III by comparing the results of the Markov chain approach with those of a discrete event simulator based on a more sophisticated applications model. At last, we prove in Section VII that the Markov chain approach can be used to evaluate other Burst-Buffer management strategies, such as those recently introduced to limit transfers between the Burst-Buffer and the storage system. We give a complete summary of assumptions and results in Section VIII and we propose concluding remarks in Section IX. Note that all our results are validated through thorough evaluation. The code source and scripts for those evaluations are available at https://gitlab.inria.fr/ordo-bdx/io-peak.

II. RELATED WORK A. Burst-Buffer Architectures and models

There are many implementations of Burst-Buffers. The two most studied characteristics are the location of the buffers and whether the buffers are shared between multiple applications.

Typically, Burst-Buffers can be located between the compute nodes and the Parallel File System (PFS). This is the case of DDN IME [START_REF] Ddn Storage | BURST BUFFER & BEYOND, I/O & Application Acceleration Technology[END_REF], [START_REF] Schenck | Evaluation and performance modeling of a burst buffer solution[END_REF] and Cray DataWarp [START_REF] Henseler | Architecture and design of cray datawarp[END_REF], [START_REF] Mubarak | Quantifying i/o and communication traffic interference on dragonfly networks equipped with burst buffers[END_REF], [START_REF] Daley | Performance characterization of scientific workflows for the optimal use of burst buffers[END_REF]. In this pseudo-centralized architecture, the Burst-Buffers are often colocated with the I/O nodes. The management strategy can then differ. Mubarak et al. [START_REF] Mubarak | Quantifying i/o and communication traffic interference on dragonfly networks equipped with burst buffers[END_REF] study the case where the buffers are shared between the different applications on the platform and used to accelerate transfers and to prevent I/O congestion. On the contrary, in Schenck et al. [START_REF] Schenck | Evaluation and performance modeling of a burst buffer solution[END_REF] and Daley et al. [START_REF] Daley | Performance characterization of scientific workflows for the optimal use of burst buffers[END_REF], applications decide the size of the buffer that should be dedicated to them.

Another solution is a distributed version of Burst-Buffers where the buffers are allocated close to the compute nodes [START_REF] Kimpe | Integrated in-system storage architecture for high performance computing[END_REF], [START_REF] Bautista-Gomez | Fti: high performance fault tolerance interface for hybrid systems[END_REF], [START_REF] Tang | Toward managing hpc burst buffers effectively: Draining strategy to regulate bursty i/o behavior[END_REF]. A solution consists in allocating the distributed buffers to the applications using the compute nodes close to buffers [START_REF] Cao | To share or not to share: comparing burst buffer architectures[END_REF]. However, other strategies focus on how to share them between the different applications [START_REF] Mubarak | Quantifying i/o and communication traffic interference on dragonfly networks equipped with burst buffers[END_REF], [START_REF] Kimpe | Integrated in-system storage architecture for high performance computing[END_REF], [START_REF] Bautista-Gomez | Fti: high performance fault tolerance interface for hybrid systems[END_REF]. This is particularly true in the context of fault-tolerance, where using a buffer on a different node can allow the implementation of hierarchical checkpointing strategies that provide more resilience than in-node buffer strategies [START_REF] Bautista-Gomez | Fti: high performance fault tolerance interface for hybrid systems[END_REF]. Furthermore, in the case where, because of their costs, the number of buffers in the machine has to be limited, one must choose on which node they should be deployed and between which subset of applications they should be shared.

In this paper, we focus on the pseudo-centralized model, where Burst-Buffers are shared between applications. We use a model of architecture similar to the one shown by Schenck et al. [START_REF] Schenck | Evaluation and performance modeling of a burst buffer solution[END_REF]. We consider that applications are running on Compute Nodes. They use or generate data (also called I/O in this work) that has to be sent to the Parallel File System. This is done by sending their data to the I/O nodes where Burst-Buffers are located.

More precisely, the architecture we consider is depicted in Figure 1.

• The I/O nodes communicate with the PFS using a maximum bandwidth B. • The Burst-Buffers are located next to the I/O nodes, their total size is denoted by S (in GB) and they communicate with the I/O nodes with a bandwidth B BB .

B. Algorithms to deal with Burst-Buffers

When it comes to using Burst-Buffers, several solutions have been proposed. We present and discuss the most common ones. A natural idea is to use Burst-Buffers as a cache to improve the I/O-performance of applications [START_REF] Schenck | Evaluation and performance modeling of a burst buffer solution[END_REF]. For instance, DDN [START_REF] Ddn Storage | BURST BUFFER & BEYOND, I/O & Application Acceleration Technology[END_REF] announces bandwidth performance 10-fold that of PFS using their Burst-Buffers. The idea is to move the I/O to the Burst-Buffer as a temporary stage between compute nodes and the Parallel File System (whether the data is incoming or outgoing). Thanks to the higher bandwidth of the Burst-Buffers, this has the advantage of improving the I/O transfer time while pipelining the (slowest) phase of sending/receiving data from the PFS with the compute phase of the application. However, as was noted by Han et al. [START_REF] Han | Accelerating a burst buffer via user-level i/o isolation[END_REF], this idea is not viable, (i) Burst-Buffers are based on technologies that are extremely expensive with respect to hard drives, (ii) they are currently based on SSD technology, that is known to have a limited rewrite lifespan [START_REF] Han | Accelerating a burst buffer via user-level i/o isolation[END_REF]. Thus, the large number of I/O operations in HPC applications would decrease their lifespan too fast. This is why we do not consider this solution in this work.

The second natural idea proposed in the literature is to use Burst-Buffers to prevent I/O congestion [START_REF] Liu | On the role of burst buffers in leadership-class storage systems[END_REF], [START_REF] Herbein | Scalable i/oaware job scheduling for burst buffer enabled hpc clusters[END_REF] while maintaining their lifespan. To achieve this goal, the applications use the direct link to the PFS (see Figure 1) when its bandwidth B is not exceeded. When the bandwidth is exceeded by the set of transfers, then the higher bandwidth of the Burst-Buffer is used to complement the bandwidth of the PFS. This is one of the solutions advocated by DDN in [START_REF] Ddn Storage | BURST BUFFER & BEYOND, I/O & Application Acceleration Technology[END_REF]. The intuition behind this strategy is that the average use of PFS bandwidth is usually small enough, but that Burst-Buffers are crucial to deal with applications' (simultaneous) bursts. This corresponds to the model depicted in Figure 1 that will be used throughout this paper.

Finally, a large part of the literature on Burst-Buffers shows how to use them with a specific application workflow [START_REF] Daley | Performance characterization of scientific workflows for the optimal use of burst buffers[END_REF], [START_REF] Schenck | Evaluation and performance modeling of a burst buffer solution[END_REF]. Specifically, they consider systems where applications have dedicated and pre-allocated Burst-Buffers, and where the application can explicitly control its data transfers and the use of the Burst-Buffer. This must to be done for each application and is very platform dependent. In practice, only few applications have the human-power to implement this. By opposition, our work is only architecture dependent and does not require any additional work from applications developers. However, we believe our results can also be used by applications developers if they want to estimate the size of Burst-Buffers that they would need based on their application caracteristics.

C. Application Model

Many recent HPC studies have independently observed patterns in the I/O behavior of HPC applications. The quasi-periodicity of HPC applications has been well documented [START_REF] Carns | 24/7 characterization of petascale I/O workloads[END_REF], [START_REF] Aupy | Periodic i/o scheduling for super-computers[END_REF]. HPC applications alternate between computation and I/O transfer phases, this pattern being repeated almost identically over time.

Hu et al. [START_REF] Hu | Storage wall for exascale supercomputing[END_REF] summarized the four key characteristics of HPC applications observed in the literature.

1) Periodicity: Applications alternate between compute phases and I/O phases. Furthermore they do so in a quasi periodic fashion. 2) Burstiness: In addition to the observed periodicity, sometimes, short I/O bursts occur. 3) Synchronization: I/O transfers of an application are performed in a synchronized way between the different parallel processes. 4) Repeatability: The same jobs are often run many times with different inputs, hence their compute-I/O pattern of an application can be reasonably predicted before execution. When modeling applications, most Burst-Buffer related work use workload models based on these patterns [START_REF] Herbein | Scalable i/oaware job scheduling for burst buffer enabled hpc clusters[END_REF], [START_REF] Han | Accelerating a burst buffer via user-level i/o isolation[END_REF], [START_REF] Mubarak | Quantifying i/o and communication traffic interference on dragonfly networks equipped with burst buffers[END_REF]. In addition to this, Mubarak et al. [START_REF] Mubarak | Quantifying i/o and communication traffic interference on dragonfly networks equipped with burst buffers[END_REF] introduce a random background traffic representing HPC workloads such as graph computations and linear algebra solvers, based on the work of Yuan et al. [START_REF] Yuan | A new routing scheme for jellyfish and its performance with hpc workloads[END_REF].

III. MODEL

We consider a large computing platform, and we focus on modeling and analyzing the behavior of the storage system. We will thus ignore the computing nodes of the platform, and simply consider that we have a (fixed) set of n applications A i currently executing on the platform. We model the long-term storage system as a single file server with input bandwidth B (in a real system with several file servers, B would represent their aggregate bandwidth). In addition to this file server, the platform contains a Burst-Buffer, whose input bandwidth is significantly larger, and whose size is S. Figure 1 provides a schematic view of the model we consider.

Following the related work discussed in Section II-C, application A i alternates computation and I/O phases. When it performs computations it does not induce any load on the storage system. When it sends intermediate results to the disks, it does so with a fixed bandwidth denoted by b i . We do not assume that the applications follow a perfect periodic pattern. In the simulations proposed in Section VI, we add a certain amount of noise to model pseudo periodicity: for a phase of expected duration T i , we actually set its duration to a value chosen uniformly at random in [(1 -)T i , (1 + )T i ]. Finally, we denote the quasi-period of application A i by d i , and the proportion of time spent doing data transfers by p i , so that the expected duration of a data transfer for application A i is p i d i (and the phases where only computation happens have an expected duration of (1 -p i )d i ). Note that we do not specify whether an application overlaps the communication with some computation or not.

To obtain realistic values, we concentrate on the workflows described in the APEX report [START_REF] Lanl | Apex workflows[END_REF]. At LANL, most of the load comes from 4 applications, whose characteristics are given in Table I. To obtain the checkpointing period (not provided in [START_REF] Lanl | Apex workflows[END_REF]), we rely on theoretical works [START_REF] Aupy | Checkpointing algorithms and fault prediction[END_REF] to determine the optimal checkpointing period from the MTBF and the checkpointing time.

In order to obtain theoretical results, we model application data transfers with a random process. To achieve this, we omit the quasi-period d i in the model. More specifically, we consider discretized time units and we assume that during each of these time units, application A i sends data with probability p i (with bandwidth b i ). In order to have a time unit corresponding to the caracteristic size of the system, we set it as the average value of the data transfer times (p

i d i ) in what follows.
Therefore, in our model, all applications share a common time unit, and there is no correlation (no memory) between what happens at time step t and t + 1. This assumption is of course crucial to build a Markov chain model. However, if the length of a data transfer for A i is much longer than the time unit, the fact that A i is involved in a communication at time t strongly influences the probability that it is involved in the same communication at time t + 1. On the other hand, if the period of the pattern for A j is much shorter than the time unit, then the I/O bandwidth consumed during one time unit with our model is very imprecise, since it is either sending or not sending during the whole time unit, whereas such an A j actually performs several communication and computation phases. For this reason, we assume in this work that all applications share a similar characteristic time.

In order to validate the ability of the model to predict the load of the Burst-Buffer, we rely on simulations based on a discrete event simulator based on the actual characteristics of the applications, both on synthetic data and on the trace based on [START_REF] Lanl | Apex workflows[END_REF]. In this trace, the time period varies in a ratio of about 3. In Section VI, we consider time periods that vary in a ratio of 10 and show that this assumption is valid to predict Burst-Buffer needs, which validates the prediction capacity of the Markov-based model.

An important quantity in our model is the expected load EXPECTEDLOAD, defined by EXPECTEDLOAD = i p i b i . It represents the average bandwidth requirement over a long period of time. If many applications send data at the same time and exceed the bandwidth B to the PFS, the excess can be sent to the Burst-Buffer, and later on sent back to the disk when some bandwidth is available. If the Burst-Buffer is full, this will induce some sort of contention and idle time, which is what we try to avoid. Section V presents a model of the system behavior to analyze the waste implied by these overflow events.

IV. INSTANTANEOUS LOAD ESTIMATION

With the model of Section III, since the expected load EXPECTEDLOAD is given by EXPECTEDLOAD = i p i b i and the available bandwidth is given by B, we can immediately notice that if EXPECTEDLOAD > B, the I/O bandwidth will not be able to cope with all the demands and the Burst-Buffer, whatever its size, will eventually be saturated. To cope with this situation, we will introduce in Section V a delay mechanism where the applications freeze and stop sending data to the Burst-Buffer during one time unit so as to avoid Burst-Buffer overflow.

In this section, we are rather interested in estimating the probability distribution of the INSTANTLOAD of I/Os, where INSTANTLOAD = i∈A(t) b i denotes the bandwidth required by the set A(t) of applications involved in a transfer at instant t. We discretize the possible values of the bandwidth with respect to B. More precisely, let us set B = 100 and let us assume that all b i values are integer values. This corresponds to rounding the value b i to the closest integer and induces a loss of precision of at most 1%. Given the difficulty to precisely estimate the I/O bandwidth required by an application A i when transferring data to the disk, we consider that this approximation is smaller than the noise of the measurements.

Our goal is in particular to detect when it exceeds the capacity of I/O bandwidth to the Burst-Buffer B BB . We propose two different approaches to achieve this goal. We rely in Section IV-A on generalizations of Chernoff bounds to the case of a weighted sum of independent variables, and we propose in Section IV-B an ad-hoc strategy to estimate precisely the discretized probability distribution, and we compare both strategies in Section IV-C.

A. Chernoff based bounds

In order to analyze the probability distribution of INSTANT-LOAD, let us denote by X i the random variable indicating whether A i is sending I/Os. We assume that X i s are independent random variables, what is in general the case if workflows and coupled codes are considered as distinct applications. Using the model described in Section III, we have P r(X i = 1) = p i and P r(X i = 0) = (1 -p i ) and the random variable X which represents the instant load is given by X = i b i X i and its expected value is E(X) = i b i p i . We cannot apply directly Chernoff bounds since X is a weighted sum, but several extensions have been proposed in the literature to deal with weighted sums (see [START_REF] Chung | Connected components in random graphs with given expected degree sequences[END_REF], [START_REF]Concentration inequalities and martingale inequalities: a survey[END_REF], [START_REF] Habib | Probabilistic methods for algorithmic discrete mathematics[END_REF], [START_REF] Vöcking | Selfish load balancing[END_REF]). Theorem 1 provides a way to estimate the probability that the instant load differs by an additive value λ from its expected value.

Theorem 1 ( [START_REF]Concentration inequalities and martingale inequalities: a survey[END_REF]). Let X 1 , . . . , X n be independent random variables such that P r(X i = 1) = p i and P r(

X i = 0) = (1 -p i ) and let X = i b i X i . Let us set ν = i b 2 i p i and b = max{b 1 , . . . , b n }. Then, P r(X ≤ E(X) -λ) ≤ exp - λ 2 2ν and P r(X > E(X) + λ) ≤ exp - λ 2 2(ν + bλ/3) .

B. Ad-Hoc Estimation

Chernoff bounds are known to model well the asymptotic behaviors, but due to their simplicity and generality, they fail to provide useful estimations in practical cases as they overestimate the probability that X differs from its expected value. Therefore, we propose in this section an ad-hoc way to estimate the probability distribution of X.

We describe in Algorithm 1 the code used to compute the probability distribution DIST of X = i b i X i , where DIST(k) is the probability that A 1 , . . . , A n contribute to exactly k. To compute DIST, we add applications A i one after the other and we update at step i the probability distribution using only applications A 1 , . . . , A i . The size of vector DIST is i b i , i.e. the maximal possible instant value which is upper bounded by nB = 100n, such that the overall complexity of Algorithm 1 is 100n 2 , where n denote the number of applications and is therefore expected to be very small (n = 20 is already large in practical settings), so that computing DIST is immediate.

(k) = (1 -pi)DIST(k) + piDIST(k -bi) if k ≥ bi ; 5 DIST(k) = (1 -pi)DIST(k) if k < bi ;

C. Estimating the distribution

We can use the Chernoff bounds proposed in Section IV-A to estimate the distribution of the instantaneous load. Indeed, let us set

• f u (x) = exp -( bipi-x) 2 2ν if x ≤ b i p i and f u (x) = 1 otherwise ; • f l (x) = exp - (x-bipi) 2 2(ν+b(x-bipi)/3) if x ≥ b i p i and f l (x) = 0 otherwise. Then, ∀x, f l (x) ≤ P r(X ≤ x) ≤ f u (x).
To estimate the precision of these inequalities, we have compared them to the result provided by Algorithm 1 on some randomly generated set of applications. We target an expected load of M 100(= B) and we generate applications in the following way. We add applications one by one by generating the values p i uniformly between 0 and p max and the values b i values uniformly between 0 and 1 until the expected load is reached. Then, we rescale the p i s so as to achieve exactly M = 100(= B). Figure 2 provides the comparison of the obtained distribution estimations for different values of p max . We can see that the Chernoff bounds are always correct, but very imprecise. Since computing the exact distribution with Algorithm 1 is in practice very fast, we will rely on this algorithm in the rest of the paper.

V. BURST-BUFFER LOAD ESTIMATION

In Section IV-B, we proposed Algorithm 1 to compute the probability DIST(k) of having INSTANTLOAD = k ∀0 ≤ k ≤ b i . As shown in Section IV-C, DIST(k) can be used directly to evaluate the probability of exceeding the I/O bandwidth to the disk. In Section V-A, we use DIST in order to evaluate the probability of exceeding the capacity of the Burst-Buffer, using a Markov Chain to model the load of Burst-Buffer.

A. Markov Chain Model

In this section, we propose to model the occupation of the Burst-Buffer when using Algorithm 2.

Following the model described in Section II-B, Algorithm 2 transfers as much data as possible to the disk. If INSTANTLOAD ≤ B, then the Burst-Buffer is emptied as fast as possible, provided that the I/O bandwidth to the disk B is not exceeded (see Section III and Figure 1). If INSTANTLOAD > B, then the occupation of the Burst-Buffer increases. If the capacity of the Burst-Buffer is exceeded, then all transfers are stopped for one time unit so as to empty the Burst-Buffer at rate B.

The occupation of the Burst-Buffer when using Algorithm 2 can be easily modeled by a Markov Chain [START_REF] Stewart | Introduction to the numerical solutions of Markov chains[END_REF]. Let us denote by Y j , 0 ≤ j ≤ S the "normal" state of the chain when the size of Burst-Buffer is exactly j and let us define i b i additional "overflow" states Y j , S +1 ≤ j ≤ S + i b i that correspond to the states where the capacity of the Burst-Buffer is exceeded.

Using DIST, we can compute the probability to move from State When the chain is in an overflow state Y i , i ≥ S, then applications are kept idle for one time unit and the Burst-Buffer is emptied at maximal rate, so that we jump to (normal or overflow) state Y i-B . Therefore, the probability P j,l of transition between state Y j and state Y l is defined as follows:

Y j to State Y l , ∀ 0 ≤ j, l ≤ S + i b i . When in normal state Y i , i ≤ S, if INSTANTLOAD = k, then • If i + k ≤ B,
• If 0 ≤ j ≤ S and 1 ≤ l ≤ S + i b i , then P j,l = DIST(l -j + B). • If 0 ≤ j ≤ B, then P j,0 = B-j k=0 DIST(k). • If S + 1 ≤ j ≤ S + i b i , then P j,j-B = 1.

B. Steady State Load Estimation

The Markov chain Y introduced in Section V-A can be used to estimate the idle time induced by Algorithm 2. Let us first prove with Theorem 2 that Y has a stationary probability π.

Theorem 2. There exists a unique stationary distribution π for the Markov chain Y .

Proof. Since p i < 1 for all i, we have DIST(0) > 0. This implies that each state Y i has a non-zero probability to jump to state Y max(i-B,0) , and thus that state 0 is accessible from any state Y i . The Markov chain thus has only one final class (the one containing state 0). This implies that Y has one and only one stationary distribution. Additionally, since P 0,0 > 0, this class is aperiodic and this distribution can be obtained as the limit of π l = P l π 0 , independently of π 0 .

The size of the transition matrix P is S + b i × S + b i and π is the solution of the linear system π = π * P . In all the experiments that follow, since we set B = 100 in Section IV-B, the size of P is therefore of a few hundreds. Moreover, in all our experiments, we use an iterative algorithm π l = π l-1 P that always converged in less than 10 iterations. Therefore, the computation of the stationary distribution is extremely fast. Theorem 3. The idle time induced by Algorithm 2 is given by k>S π(k). Proof. π(k) represents the probability of being in state Y k . For k ≤ S, the load of the Burst-Buffer is smaller than its size S, and no idle time is induced. On the other hand, states Y k , k > S correspond to an overflow of the buffer. In this case, Algorithm 2 stops all transfers from the applications with probability 1 in order to empty the Burst-Buffer for one timeunit. Therefore, the fraction of idle time corresponds to the probability of being in a state Y k , k > S and is therefore given by k>S π(k).

C. Idle Time Prediction

This section presents the results obtained with the Markov chain model presented above. We generate applications in the same way as in Section IV-C and we observe the resulting stationary distribution. We fix the bandwidth B of the disk to 100, and we vary the expected load M of the system, by considering different values of α = M B . As a reminder, once M is fixed, we add applications by picking uniformly at random p i ∈ [0, p max ] and b i ∈ [0, B] until i b i p i ≥ M and then we rescale the p i s so that i b i p i = M = αB. Small values of α result in an underloaded system in which the buffer does not fill very often. On the other hand, values of α larger than one mean that idle times are required even in presence of a very large buffer in order to allow enough time to write to the disks.

Figure 3 shows the resulting proportion of idle time (see Theorem 3) for different values of p max and α. We can observe that obtained results do not strongly depend on p max . and that, as expected, the idle time decreases with the buffer size, converging to 0 for values of α below 1, and to a positive value for α above 1.

VI. MODEL VALIDATION AND SIMULATIONS

A. Comments on the probabilistic model

In this section, we aim at evaluating the influence of the application model described in Section III on the quality of the prediction of the idle time induced by Algorithm 2. Indeed, we have proved in Section V that this model is theoretically and numerically tractable since a simple Markov chain enables to compute the idle time at low computational cost. Nevertheless, the model of Section III is not strictly equivalent to the models proposed in Section II-C. Indeed, it is generally assumed that HPC applications follow quasi periodic pattern alternating compute and communication phases and that the rate of emission is proportional to the number of I/O nodes and is therefore constant over time. Our model captures several crucial characteristics of the quasi periodic model, (i) alternating communication and computation phases, (ii) the quasi periodic pattern modeled by p i probability of I/O communication (iii) the constant bandwidth b i used by a given application A i when performing I/Os. Nevertheless, as underlined in Section III, the Markov-based model requires a common time unit for all applications and induces more randomness, which may induce modeling errors. It is therefore crucial to validate our probabilistic model against a more realistic (but theoretically intractable) model using a discrete event simulator that will be presented in Section VI-B.

B. Discrete Event Simulator

In order to validate the ability of the application model presented in Section III to predict the idle time, we build a discrete event simulator. The platform model is the same as the one introduced in [START_REF] Schenck | Evaluation and performance modeling of a burst buffer solution[END_REF] and described in Figure 1, and is characterized by the I/O bandwidth B between the processing nodes and the disks, that is shared by the transfers between the Burst-Buffer and the disk, so that the overall bandwidth to the disk is bounded by B.

On the application side, simulations rely on a more sophisticated model. Each application is characterized by a time period d i , a probability of sending I/Os p i and a used bandwidth b i , and we add another parameter u that will be used to add noise to the pattern of A i . More precisely, a pattern for A i is built as follows

• A i computes during (1 -p i )d i x where x is chosen uniformly at random in [1 -u, 1 + u] at each computing burst, • A i sends data at rate b i during p i d i x where x is chosen uniformly at random in [1 -u, 1 + u] at each burst of I/Os. Using above model, we generate a pattern that is almost periodic (exactly periodic if u = 0), and the quasi-period of A i is different for each application, so that the model is close to what has been observed in HPC systems (see Section II-C and references therein).

Our discrete event simulator is described in Algorithm 3. comp i is a boolean value that indicates whether the next burst for A i is a I/O (comp i = 0) or a compute (comp i = 1) burst.

Algorithm 3: Discrete Event Simulator

Input: The set of active applications A k 1 , . . . , A k l , with their characteristics pi, bi and di and their random starting times ti a large time bound Tmax, the regularity of bursts length u Output: The total idle time IDLETIME 1 ∀i, compi = 1 NEXTEVENT = min(t1, . . . , tn) ; 2 while NEXTEVENT ≤ Tmax do In order to compare the results obtained with the Markov chain model to those of the discrete event simulator, we once again generate applications for different values of p max and α as was done in Section V-C. In addition to p i and b i values, we also generate the values d i uniformly at random between some value d min and βd min , where β = 1, 2, 5, 10. β = 1 corresponds to the case where all applications share the same period, whereas β = 10 models applications with different dynamics. The value of d min is chosen in each case so that the expected value of d i is 10. The results are depicted on Figure 4. For each scenario, we run Algorithm 3 10 times up to T max = 1000, and we show on the plots the confidence intervals for the mean idle time with a filled ribbon. We can see that our model is able to predict the behavior of the Idle Time induced by limited size of the Burst-Buffer, even for large β values. The Markovian model tends to overestimate the Idle time, especially in the case when α = 1, i.e. when the load of the system actually corresponds, on average, to the capacity of the storage system. This validates the ability of the Markov-based model introduced in Section III to predict the load of the b.

b) On the APEX trace: We performed simulations using data presented in [START_REF] Lanl | Apex workflows[END_REF] for LANL and described in Section III. Using APEX trace, the expected load is small with respect to the bandwidth of the system (it corresponds to a value α = 9%). In order to validate the model into a more stressed context, we also scaled p i values so as to obtain α = 0.75, 1 and 1.25. The results are depicted on Figure 5. Again, the idle time predicted by the Markovian model has the correct behavior, but slightly overestimates the idle time, especially when α = 1.

VII. LAZY EMPTYING

A. Algorithm

To prevent from the phenomena of Write Amplification [START_REF] Hu | Write amplification analysis in flash-based solid state drives[END_REF], Han et al. [START_REF] Han | Accelerating a burst buffer via user-level i/o isolation[END_REF] consider lazy strategies for emptying the Burst-Buffer. In particular, they propose to empty the Burst-Buffer only when its load reaches a certain level. The rationale behind this is that performing a transfer between the Burst-Buffer and the disk has a fixed cost that does not only depend on the transferred data. Therefore, there is a clear interest in keeping the Burst-Buffer to disk link quiet as often as possible. Algorithm 4 is based on this idea. If the load of the Burst-Buffer is smaller than a given THRESHOLD < S, the Burst-Buffer is not emptied, even if some bandwidth remains available once the transfers between the computing nodes and the disk have been allocated.

In this context, we have two conflicting objectives. Indeed, we are interested in both maximizing the quiet time of the Burst-Buffer (denoted by QUIET) and to minimize the idle time induced by Burst-Buffer overflow (denoted by IDLETIME). Clearly, a high value of THRESHOLD is expected to make QUIET larger (since we keep the lazy strategy longer), but also to make IDLETIME larger (since we take more risks by not emptying the buffer). The details of the algorithm are presented in Algorithm 4. In order to enforce the properties of the associated Markov chain (see Section VII-B and Theorem 4), we keep a (very small probability) of emptying the Burst-Buffer, even if THRESHOLD < S.

B. Markov Chain associated to Algorithm 4

Similarly to what has been proposed in V-A, it is easy to associate a Markov chain to Algorithm 4, where state Y LAZY j corresponds to a load j of the Burst-Buffer and whose transition matrix P LAZY is defined as follows

• If THRESHOLD ≤ j ≤ S and 1 ≤ l ≤ S + i b i , then P LAZY j,l = DIST(l -j + B). • If 0 ≤ j < THRESHOLD:

if 0 ≤ l < j then P LAZY j,l = 0.01 × DIST(l -j + B);

Figure 1 :

 1 Figure 1: Modeling of the pseudo-centralized platform.

Algorithm 1 :

 1 FINDDISTRIBUTION(A 1 , . . . , A n ) Input: The characteristics of A1, . . . , An Output: DIST, such that DIST(k) = P r( biXi = k) 1 DIST is a initialized with zeros and its size is bi ;

Figure 2 :

 2 Figure 2: Comparison of distribution estimations for different values of p max , EXPECTEDLOAD = 100.

Algorithm 2 : 2 l 1 b 3 Transfer 6 Transfer 7 if LOAD + l 1 b k i -B ≤ S then 8 l 1 b 11 Transfer

 2213678111 ALLOCATETRANSFERSInput: The set of active applications A k 1 , . . . , A k l , the load of the Burst-Buffer LOAD Output: The set of transfers to the disk and to the Burst-Buffer1 if l 1 b k i ≤ B thenTransfer k i from the computing nodes to the disk ; min(LOAD, (B -l 1 b k i )) from the Burst-Buffer to the disk ; B from the computing nodes to the disk ; Transfer k i -B from the computing nodes to the Burst-Buffer; S -LOAD from the computing nodes to the Burst-Buffer;12 Stop all applications until the transfer is over ;13During each time unit, transfer B from the computing nodes to the disk ;14Complete with Burst-Buffer data at the end ;

  then the Burst-Buffer becomes empty and we jump to state Y 0 . • If i + k -B ≤ S, then we jump to normal state Y i+k-B . • If i+k -B > S, then we jump to overflow state Y i+k-B .

Figure

  Figure Proportion of idle time as a function of buffer size, for different values of p max and α.

3 i 4 5 if compi == 1 then 6 tioccupation at time ti 14 while

 345614 = argmin(t1, . . . , tn) ; Choose x uniformly at random in [1 -u, 1 + u] ; = NEXTEVENT + (1 -pi)dix; Burst-Buffer occupation at time ti ≥ S do 15 Stop all transfers for one time unit; synthetic data:

Figure 4 :

 4 Figure4: Comparison of predicted idle times between Markov chain (dashed lines) and discrete event simulator -Synthetic Data (full lines). ratiomax is the ratio between the bounds of the interval in which the different d i are chosen in each simulated set.

Figure 5 : 1 l=0

 51 Figure 5: Comparison of predicted idle times between Markov chain (dashed lines) and discrete event simulator -APEX data (full lines)

Theorem 4 . 2 l 1 b k i from the applications to the disk ; 3 if LOAD ≥ THRESHOLD then 4 Transfer 7 11 Transfer 16 Transfer 17 Stop all applications until the transfer is over ; 18 During 5 .

 42347111617185 There exists a unique stationary distribution π LAZY for the Markov chain Y LAZY .Proof. The difference between Y and Y LAZY j lies in the weights only, and not in the structure of the transitions. Therefore, the proof of Theorem 2 directly applies to Theorem 4 if we define π LAZY as the solution of the linear system π LAZY = π LAZY × P LAZY . Algorithm 4: Lazy ALLOCATETRANSFERS Input: The set of active applications A k 1 , . . . , A k l , the load of the Burst-Buffer LOAD 1 if l 1 b k i ≤ B then Transfer min(LOAD, (B -l 1 b k i )) from the Burst-Buffer to the disk ; With probability 0.01, transfer max(LOAD -(B -l 1 b k i ), 0) from the Burst-Buffer to the disk ; B from the computing nodes to the disk ; 12 if LOAD + l 1 b k i -B ≤ S then 13 Transfer l 1 b k i -B from the computing nodes to the Burst-Buffer; S -LOAD from the computing nodes to the Burst-Buffer; each time unit, transfer B from the computing nodes to the disk ; 19 Complete with Burst-Buffer data at the end ; The idle time induced by Algorithm 4 is given by IDLETIME = k>S π LAZY (k) and the quiet time induced by Algorithm 4 is given by QUIET =

Table I :

 I Characteristics of the APEX applications data set.

	Workflow	EAP	LAP	Silverton	VPIC
	Number of Instances	13	4	2	1
	b i (GB/s)	160	80	160	160
	d i Period (s)	5671 12682	15005	4483
	Checkpoint time (s)	20	25	280	23,4
	p i (×10 -3 )	3.51	1.97	18.7	5.11

ACKNOWLEDGEMENTS

This work was supported in part by the ANR Dash project (ANR-17-CE25-0004). The authors would like to thank Jean-Thomas Acquaviva for insightful discussions and the reviewers for useful comments.

S k=0 π LAZY (k) × l≥k P LAZY (k, l) .

Proof. The states Y LAZY k with k > S correspond to overflow states and systematically induce one time unit of idle time.

To estimate QUIET, we sum, for all possible "normal" states Y LAZY k , k ≤ S, the probability to not decrease the load of the Burst-Buffer, i.e. to jump to a state Y LAZY l , l ≥ k.

C. IDLETIME and QUIET times induced by Algorithm 4

To study the effect of the threshold value on the performance of the Burst-Buffer, we analyze the results obtained using the Markov chain model. We generate applications like previously, and we consider varying values for α, p max , and for the threshold ratio (i.e., the ratio between the THRESHOLD value and the buffer size S). A threshold ratio of 0 is equivalent to the basic case, and a threshold value of 1 means that the burst buffer is only emptied when it is full. Results for IDLETIME time are shown in Figure 6. In addition, we also measure QUIET time during which the burst buffer is not emptied (i.e., its occupation either increases or stays the same). This is shown on Figure 7. We can see that in most cases, even a small positive value of THRESHOLD allows to improve the QUIET, except when the load is high and p max is also high. In those cases, the Burst-Buffer is almost continuously kept full due to the constant arrival of application data, and both strategies behave similarly. In all other cases, the buffer has a non-zero probability of having a small enough load for the threshold to be useful. In many cases however, a too high value of THRESHOLD (above 20 or 40% of the Burst-Buffer size) induces significant cost on IDLETIME time. Since even a small THRESHOLD brings a benefit on the QUIET times (and for low load, increasing THRESHOLD further does not degrade IDLETIME time), we can see that this lazy strategy indeed has the potential to increase the performance of the Burst-Buffer. It is important to correctly dimension the THRESHOLD value, however a 20% value appears to be a good safe guess.

VIII. SUMMARY OF RESULTS AND RECOMMENDATIONS

In this section, we provide (i) a summary of the assumptions underlying the Markovian model and (ii) a summary of the results and findings that we obtained using it.

In this paper, we focus on using Burst-Buffers solely as a temporary storage unit for data destined to go on the PFS. The Burst-Buffer is not used as a cache for the PFS, but as a buffer when the bandwidth to PFS is saturated. In particular, we leave for future work the ability to store intermediate data on the Burst-Buffer, that may or may not be committed to the PFS. We assume that this is a system-wide process, and that there is no application-specific strategy to use the Burst-Buffer.

As far as applications are concerned, we assume that their I/O behavior is (quasi-)periodic (most of the I/O within an application happens in chunks) and rather predictable, in terms of quantity of data sent and of freqency of sending periods. In the Markovian model, however we assume that these sending periods happen at random times. We also assume that they all share a similar characteristic time, which is often almost the case when the sending pattern comes from checkpointing.

To validate the Markovian model, we compared its results for Idle Time prediction against an ad-hoc discrete event simulator. We performed the comparison both on synthetic data covering a wide range of parameters and on an actual workload trace from LANL.

We prove in this paper that the Markovian model is tractable and can be used to predict with a very low computing cost the instantaneous load distribution and the idle time. It can also be augmented to analyze more sophisticated management strategies for the Burst-Buffer, like lazy emptying strategies.

In the crucial context of Idle Time prediction, the comparison with simulations show that the Markovian model, despite its bursty nature, is in fact able to predict the behavior of the load of the Burst-Buffer. It therefore provides a low-cost algorithm to dimension the Burst-Buffer characteristics (size, strategy, ...) on HPC systems.

IX. CONCLUSION

In this paper, we propose a probabilistic model for the use of Burst-Buffers in the context of organizing the transfers of a set of HPC quasi-periodic applications. We prove that this model enables to estimate at low computational cost the load of the Burst-Buffer under several classical management strategies. Present work opens several perspectives. First, it would be of great interest to include in the probabilistic model the specific characteristics of BigData workload (i.e. the ability to prefetch data) and checkpointing strategies (i.e. the ability of removing old checkpoints without writing them eventually to the disk), that can improve the performance of the Burst-Buffer. Then, this work makes it possible to evaluate complex Burst-Buffer management strategies where the Burst-Buffer is used differently by the different applications. At last, extending this model to more distributed Burst-Buffer architectures is also an interesting perspective.


Guillaume Damiand
email: guillaume.damiand@liris.cnrs.fr

Aldo Gonzalez-Lorenzo
email: aldo.gonzalez-lorenzo@liris.cnrs.fr

Jarek Rossignac

Florent Dupont
email: florent.dupont@liris.cnrs.fr

Hierarchical Representation for Rasterized Planar Face Complexes

Keywords: Planar polygonal meshes, Irregular representation, Hierarchical representation, Combinatorial maps, Compact representation, Topology preserving rasterization

A useful example of a Planar Face Complex (PFC) is a connected network of streets, each modeled by a zero-thickness curve. The union of these decomposes the plane into faces that may be topologically complex. The previously proposed rasterized representation of the PFC (abbreviated rPFC) (1) uses a fixed resolution pixel grid, (2) quantizes the geometry of the vertices and edges to pixel-resolution, (3) assumes that no street is contained in a single pixel, and (4) encodes the graph connectivity using a small and fixed number of bits per pixel by decomposing the exact topology into per-pixel descriptors. The hierarchical (irregular) version of the rPFC (abbreviated hPFC) proposed here improved on rPFC in several ways: (1) It uses an adaptively constructed tree to eliminate the "no street in a pixel" constraint of the rPFC, hence making it possible to represent exactly any PFC topology and (2) for PFCs of the models tested, and more generally for models with relatively large empty regions, it reduces the storage cost significantly.

Introduction

Consider a planar graph, G, that is embedded in the plane and comprises a connected network of finite and possibly curved edges and their bounding vertices. For example, each edge may represent a street and each vertex may represent a street junction. Their union decomposes the plane into faces that may be topologically complex. For example, G may have multi-edges (more than one edge joining any given pair of vertices). Furthermore, a face may contain, in its boundary, cracks (edges that bound a single face), dead-ends (vertices that bound a single edge-a crack), and loops (edges that start and end at the same vertex). The unbounded face is called exterior. An example is shown in Fig. 1. We use the term Planar Face Complex (PFC) for such an arrangement.

Many applications need to represent and to traverse a PFC. Examples include street networks in Geographic Information System (GIS) [START_REF] Schmidt | Forest point processes for the automatic extraction of networks in raster data[END_REF], geological models [START_REF] Castanié | VolumeExplorer: Roaming large volumes to couple visualization and data processing for oil and gas exploration[END_REF], overlapping SVG elements [START_REF] Dalstein | Vector graphics complexes[END_REF], and multi-material structures [START_REF] Kwok | Geometric analysis and computation using layered depth-normal images for three-dimensional microfabrication[END_REF].

Different solutions have been proposed to represent PFCs. Some approaches describe the connectivity of the graph [START_REF] Baumgart | Winged-edge polyhedron representation[END_REF][START_REF] Mäntylä | An Introduction to Solid Modeling[END_REF][START_REF] Alumbaugh | Compact array-based mesh data structures[END_REF] explicitly. This may yield a high storage cost for complex graphs. Other approaches use an image format (regular grid of pixels) to describe a rasterized approximation of the PFC [START_REF] Samet | Connected component labeling using quadtrees[END_REF][START_REF] Benouamer | Bridging the gap between CSG and Brep via a triple ray representation[END_REF][START_REF] Szymczak | Piecewise regular meshes: Construction and compression[END_REF][START_REF] Ju | Dual contouring of Hermite data[END_REF], which assigns each pixel to a different face, without attempting to capture the topology inside the shared pixels that contains one or more edges.

The recently proposed rasterized Planar Face Complex (rPFC) [START_REF] Damiand | Rasterized planar face complex[END_REF] unifies these approaches by defining a compact representation of the topology of the PFC that decomposes it into per-pixel descriptors, each using a short string of bits to encode the topology of the intersection of the PFC with a pixel.

The rPFC model has many advantages: (1) It represents graph connectivity exactly and hence supports exact topological graph traversal; (2) It provides spatial indexing to both quantized geometry and exact topology; (3) It can represent nontrivial topology in a pixel (such as dangling edges, multiple vertices and multiple connected components); (4) It requires only a few bits per pixel.

However, the rPFC has a drawback: It cannot represent a graph that has an edge that fits entirely inside a pixel. Hence, to represent a graph with some relatively small edges, we either must use a high-resolution grid (see Fig. 1), which increases storage cost, or must simplify the graph by collapsing small edges in a preprocessing step, which implies the loss of the original topology. Furthermore, when the rPFC encoding stores a topology descriptor for each private pixel (a pixel that lies entirely in a face), the rPFC storage of large clusters of private pixels is wasteful.

Motivation

Our overarching motivation is to reduce the storage size of this graph, while preserving the benefits provided by the previously proposed rPFC representation, namely (1) random access and traversal (RAT) at constant amortized time (CAT) cost and (2) constant cost localization of the edges and vertices that intersect any given pixel. We also wish to provide efficient support for distributed processing, window-stream processing, and progressive refinements.

We believe that the above characteristics are important for navigation, query, and maintenance applications of huge databases of planar graphs, which may represent the geometry and connectivity of streets, rivers, or utility networks.

Our second main motivation is to use the 2D representation as the main tool to define 3D compact representation. This paper is the first step, necessary for the definition of a compact representation of 3D meshes.

Contribution

The high-level, novel contribution reported here is the combination of a hierarchical representation with the rPFC (per-pixel) encoding of geometry and connectivity. In this paper, we define the hierarchical rasterized Planar Face Complex representation (hPFC) which addresses the drawback of the previous rPFC. The proposed hPFC is essentially a tree. Hence, our solution includes a quadtree as a special case. At the coarsest level, it is an rPFC, A. But some of the pixels of A, instead of containing the bit-string that encodes the local topology of the PFC, contain an index to a refined rPFC of the portion of the PFC inside that pixel. Such a more detailed rPFC, B, may, in turn contain pixels which, each, refer to even finer rPFC, C, and so on recursively.

This irregular representation allows us to remove the "no small edge" constraint imposed by the rPFC: When an edge fits entirely in a pixel, the pixel is subdivided.

Moreover, using an irregular (hierarchical) grid allows to reduce the storage cost of large clusters of private pixels.

For example, the rPFC shown in Fig. 1 uses a grid of 192 pixels and involves 79 crossings (points where an edge of G crosses a pixel border). A coarser grid would produce an invalid rasterization in which at least one edge is contained in a single pixel. As shown in Fig. 2, using the irregular grid of an hPFC solves the problem: The same PFC may be encoded as an hPFC that uses only a total of 42 pixels and involves only 30 crossings.

Organization

The paper is organized as follows. In Sect. 2, we review the rPFC model and discuss other relevant prior art. In Sect. 3, we present the hPFC model and the details of the operators needed to navigate through the PFC using its hPFC representation. In Sect. 4 we give a compact encoding of hPFC that provides a good time-complexity for traversal operators, while allowing to navigate through the graph without needing to decode the whole data-structure, but only the current pixel. In Sect. 5, we present experimental results, comparing our new solution with the previous regular version.

Prior Art

Data-Structures for Polygonal Meshes

A variety of edge-based data structures have been proposed in order to represent polygonal meshes, such as Combinatorial Maps, Corner Table, Doubly Connected Edge List, Half-Edge, Surface Mesh. . . [START_REF] Baumgart | Winged-edge polyhedron representation[END_REF][START_REF] Mäntylä | An Introduction to Solid Modeling[END_REF][START_REF] Lienhardt | N-Dimensional generalized combinatorial maps and cellular quasi-manifolds[END_REF][START_REF] Alumbaugh | Compact array-based mesh data structures[END_REF][START_REF] Dobkin | Primitives for the manipulation of threedimensional subdivisions[END_REF][START_REF] Aleardi | Succinct representations of planar maps[END_REF][START_REF] Kallmann | Star-vertices: a compact representation for planar meshes with adjacency information[END_REF][START_REF] Kettner | Using generic programming for designing a data structure for polyhedral surfaces[END_REF][START_REF] Snoeyink | Tripod: a minimalist data structure for embedded triangulations[END_REF][START_REF] Sieger | Design, implementation, and evaluation of the surface mesh data structure[END_REF]. They differ in their storage cost, in the type of operators that they support, and in the topological restrictions that they impose on the mesh. Many are reviewed in [START_REF] Botsch | Polygon Mesh Processing[END_REF][START_REF] Damiand | Combinatorial Maps: Efficient Data Structures for Computer Graphics and Image Processing[END_REF].

These data structures provide Random Access and Traversal (RAT) of the meshes, often in constant time, or sometimes in Constant Amortized Time (CAT). Their main drawback is to use a large number of bits per element (edge, vertex), which limits their applicability and performance for complex meshes. Some solutions used rasterized images, where each pixel stores the color of the region that contains its center. The image can be compressed, for example by using RLE (Run Length Encoding). But the digitization does not represent street networks, removes all cracks and dead ends, and can disconnect regions.

Rasterized images were used in [START_REF] Nehab | Random-access rendering of general vector graphics[END_REF] to accelerate the rendering of antialiased vector graphics. That approach uses a coarse lattice in which each cell contains a variable-length encoding of the graphics primitives that overlap it. The proposed hPFC extends this previous work by capturing the connectivity (incidence and ordering) of the graph in a constant-length per-cell format and hence providing support for RAT operators.

In [START_REF] Benouamer | Bridging the gap between CSG and Brep via a triple ray representation[END_REF][START_REF] Shade | Layered depth images[END_REF], a solution stores the crossing vertices between the mesh and an inter-pixel grid, and recompute (explicitly or implicitly) a simplified topology of the mesh. Such a representation can be used to accelerate Boolean operations [START_REF] Kwok | Geometric analysis and computation using layered depth-normal images for three-dimensional microfabrication[END_REF]. But it only represents an approximated topological description of the mesh.

Compact Representations of Polygonal Meshes

Several compression schemes propose to encode local mesh connectivity by using a few bits per element for polygon graphs [START_REF] Blelloch | Succinct representations of separable graphs[END_REF][START_REF] Blandford | Compact representations of simplicial meshes in two and three dimensions[END_REF] and for triangle meshes [START_REF] Campagna | Direct edges -a scalable representation for triangle meshes[END_REF][START_REF] Aleardi | Catalog-based representation of 2D triangulations[END_REF][START_REF] Aleardi | Explicit array-based compact data structures for triangulations[END_REF][START_REF] Gurung | SQuad: Compact representation for triangle meshes[END_REF][START_REF] Gurung | LR: Compact connectivity representation for triangle meshes[END_REF][START_REF] Yamanaka | A compact encoding of plane triangulations with efficient query supports[END_REF][START_REF] Alumbaugh | Compact array-based mesh data structures[END_REF][START_REF] Luffel | Grouper: A compact, streamable triangle mesh data structure[END_REF].

Often, the connectivity information is broken into a chunk per face, per edge, or per vertex. For example, the 2D version of Tetstreamer [START_REF] Bischoff | TetStreamer: compressed back-to-front transmission of Delaunay tetrahedra meshes[END_REF] organizes triangle faces into topological rings and divides connectivity information into one bit per edge (for some edges) and one or several bits per vertex. But extensions of this approach to more general (PFC) graphs would be challenging and the representation more expensive. More importantly, such schemes assume that the bits of the mesh encoding is received in a specific order. This compressed format must be decompressed first and converted into a more expensive format that is suitable for RAT in CAT.

More recent representations offer a much low storage costs while still supporting RAT in CAT for the most common access and traversal operations. For example, the Zipper format is restricted to triangle meshes, but uses on average only 6 bits per triangle and can be constructed in linear space and time [START_REF] Gurung | Zipper: A compact connectivity data structure for triangle meshes[END_REF]. Such representations rely on a specific ordering of vertices. The streamable version, Grouper [START_REF] Luffel | Grouper: A compact, streamable triangle mesh data structure[END_REF], of this approach stores about two vertex-references per triangle.

Hierarchical Representations of Polygonal Meshes

Many hierarchical solutions have been defined in order to reduce the memory space used in order to represent a mesh such as for example quadtrees [START_REF] Samet | Connected component labeling using quadtrees[END_REF][START_REF] Ayala | Object representation by means of nonminimal division quadtrees and octrees[END_REF]. [START_REF] Hoppe | Progressive meshes[END_REF] proposes a progressive mesh representation, a new scheme which allows to store and to transmit arbitrary triangle meshes. Several other hierarchical and pyramidal models were defined and used for example to represent multiresolution terrain models [START_REF] Floriani | Multiresolution models for topographic surface description[END_REF]. In [START_REF] Gandoin | Progressive lossless compression of arbitrary simplicial complexes[END_REF], a compressed encoding of 3D triangular meshes is defined, based on a hierarchy and an encoding of split operators, which allows to encode both manifold meshes but also "triangle soups". In [START_REF] Lefebvre | Compressed random-access trees for spatially coherent data[END_REF], a compressed random-access tree is used in order to represent spatially coherent data. But these representations are either for grid of pixels, or for triangle meshes.

Quadtrees were also used to represent a set of points [START_REF] Finkel | Quad trees a data structure for retrieval on composite keys[END_REF][START_REF] Orenstein | Multidimensional tries used for associative searching[END_REF] or of line-segments [START_REF] Shneier | Two hierarchical linear feature representations: Edge pyramids and edge quadtrees[END_REF][START_REF] Tamminen | The excell method for efficient geometric access to data[END_REF][START_REF] Webber | On encoding boundaries with quadtrees[END_REF]. These representations do not capture connectivity. The MX quadtree [START_REF] Hunter | Operations on images using quad trees[END_REF] does capture the connectivity of simple polygons, but does not support junction vertices with more than two incident edges. Hence, these previously proposed hierarchical representations cannot be used to describe the connectivity of PFCs.

PM-quadtrees [START_REF] Samet | Storing a collection of polygons using quadtrees[END_REF] allow us to represent PFCs. In that approach, the model is split recursively into chunks that are each sufficiently simple. The three variants proposed in [START_REF] Samet | Storing a collection of polygons using quadtrees[END_REF] use different criteria to define sufficiently simple, but, similarly to [START_REF] Ayala | Object representation by means of nonminimal division quadtrees and octrees[END_REF], require to have no more than one vertex in any quadtree leaf. Our solution proposed here uses a slightly less constraining criterion: we allow more than one connected component per leaf, each having at most one vertex in the leaf. This added flexibility may reduce the need to subdivide a pixel, and hence the overall storage cost, in situations, such as thin regions or constrictions in a face, in which a vertex lies close to another vertex that is not connected to it by edges in the leaf. Another difference is the generality of our solution: we are not limited to quadtrees, but can instead use any image size as input and any image size for refined images. Finally, graph information in hPFC in not only stored in leaves, like in PM-quadtrees, but also in inner nodes. Indeed, each node of our hPFC is a rPFC that can mix word pixels describing local part of the PFC and refined pixels that link to refined images.

rPFC

In this section, we summarize the aspects of the rPFC definition that are important to understand our contribution. We encourage the reader to consult the original work [START_REF] Damiand | Rasterized planar face complex[END_REF] for further details.

Consider a regular grid of pixels (which we define to be open faces, which do not include their borders), separated by roads (relatively open edges, which do not contain their vertices), which meet at crossroads (vertices, each having four incident cross-roads and four incident pixels). Each pixel P is incident upon four roads: west, north, east and south (labeled 'W', 'N', 'E', and 'S'). A pixel is stabbed (by G) when its intersection with G is not empty. It is private (to a face) otherwise. A crossing is an intersection point between G and a road of the grid (Fig. 1). To simplify the rPFC construction, we bias the geometric tests to ensure that no vertex of G lies on a road or cross-road and that no edge of G contains a cross-road.

In our explanation, we make references to the crossings and think of them as vertices inserted into the graph at places where an edge crosses a road. We wish to stress that we do not actually insert these crossings. We just imagine them there so that we can reference them more easily in our explanations. We use the term junction when referring to the original vertices of G, so as to distinguish them from the crossings. In the graph produced by these imaginary insertions of crossings, each connected component of the intersection of G with a given pixel comprises exactly zero or one junction. A component with no junction is a stabbing: a segment of an edge of G between two consecutive crossings. A component with a junction comprises a set of edge segments that connect the junction to some of the crossings on the boundary of that pixel.

The rPFC model orders crossings and crossroads around the pixel, associates the symbol '+' with each crossroad, and associates a symbol with each crossing. It distinguishes the following crossing types and uses a different symbol for each type:

• crossing '=' leads to a dead end vertex;

• crossing '[' starts a component that has no junction;

• crossing ']' ends a component that has no junction;

• crossing '(' starts a component connected to a junction;

• crossing 'Y' adds a branch to a junction;

• crossing ')' ends a component connected to a junction. The topology of a pixel, P, is encoded by the word formed by the concatenation of the clockwise sequence of the above symbols, starting from the lower left crossing (cf. example in Fig. 3).

The first three '+' symbols are used to separate the crossings into four sub-lists, one per crossroad. These will be referenced using labels: 'W', 'N', 'E', and 'S'. The occurrence of the 4 th '+' symbol indicates the end of the word.

Assuming a valid resolution (no edge of G is contained in a single pixel), the rPFC of G is represented by the sequence of words listed in scanline order of the corresponding pixels. Private (i.e., empty) pixels are identified by the word "++++".

Random Access to the PFC uses the word associated with a chosen pixel P to verify whether P is private, and if not, to calculate the number of connected components and hence of the number of junction vertices contained in P. Different mechanisms may be used to provide direct or indirect access to the face ID of private pixels. Vertices in P are assigned local vertex IDs.

Traversal of the graph G of the PFC uses local identifiers (IDs) of darts (sidewalks in the "edges are streets" metaphor that are oriented to have the street on their left). These IDs are assigned in clockwise order around P and grouped by road. Simple algorithms exist for mapping (i.e., for finding the local identifier of) a dart to the opposite dart (the sidewalk on the other side of the street) and for mapping (the local ID with respect to P of) a dart that exits a pixel P to the (the local ID with respect to Q of) the next dart that (continues the sidewalk past the crossing and) enters the adjacent pixel Q. One may traverse the bounding loop (which is always unique since faces of G have no hole) of sidewalks around a face by following the sequence of next maps.

Other simple algorithms have been provided to compute the next dart of a dart that enters P and the previous dart for both entering and exiting darts. Finally, given the local ID of a vertex in a pixel, one may easily compute the ID of one of the outgoing or incoming darts, and inversely, given an outgoing or incoming dart, one may compute the corresponding vertex ID.

The inserted crossings vertices are hidden from the application, which operates on the darts and vertices of the original PFC.

All these computations have cost that is proportional to the length (symbol count) of the word associated with the pixel.

When needed, mechanisms for identifying the global ID of the face that is on the right of a sidewalk may be provided, for example by using a hash table and may requiring the traversal of the bounding loop of a face.

hPFC

In this section, we propose a formal definition of the hPFC, an algorithm for constructing it, and the two basic tools allowing us to implement the RAT operators which are the navigation through pixels and the crossings identification.

Definition of hPFC and Terminology

We use the term image to refer to a regular grid of n × m pixels. Each pixel of an image stores either a word (in which case we say that it is a word pixel) or the reference to another image (in which case we say that it is a refined pixel and that the referenced image is a refined image). Note that we use a different terminology than the one proposed in [START_REF] Samet | Storing a collection of polygons using quadtrees[END_REF][START_REF] Samet | Foundations of Multidimensional and Metric Data Structures[END_REF] because our approach is not restricted to quadtrees.

An hPFC is a tree of images, having a root image I 0 that has n 0 × m 0 pixels.

For compactness, we assume that all refined images of an hPFC have the same n×m pixel size. Note that the tree becomes a quadtree [START_REF] Samet | Connected component labeling using quadtrees[END_REF] for n 0 = m 0 = 1 and for n = m = 2.

The hPFC of the PFC of Fig. 2 is shown in Fig. 4, where

n 0 × m 0 is set to 4 × 3 and n × m is set to 2 × 2.

hPFC Construction

We briefly explain how we build a hPFC. The pixels are computed recursively. We set an initial grid of m 0 × n 0 pixels. For each pixel (identified by the coordinates of its diagonal), we check that there is no edge of the PFC contained in it. If so, we compute its word. Otherwise, we subdivide it into m r × n r equal pixels and repeat this operation. In practice, we filter the edges of the PFC for each pixel to avoid unnecessary computations of intersections.

Given a pixel P, we compute its word by first finding all its crossings. For this, we intersect the edges of the PFC with the edges of the pixel. For each crossing, we keep the intersected edge of the pixel, the coordinates of the intersection, and a reference to its endpoints contained inside the pixel. Then, we use the intersected edge and the coordinates for sorting the crossings, and we use the references to assign the symbols.

Navigation Through Pixels

To support constant cost moves of a random walk over the plane and to support the implementation the graph traversal operators (which are discussed in Sect. 2.4), we need to be able to identify the adjacent pixel(s) of any given pixel.

More precisely, given a pixel P, we want to find the first pixel Q adjacent to P in a given direction D (west, north, east or south). If P has several adjacent pixels in direction D (case of pixel a in Fig. 2 for west direction), Q will be the first one along the road between P and Q, considering the reverse orientation than D in pixel P (for pixel a, the first adjacent pixel considering south direction, i.e. pixel b).

A pixel P is denoted by a triplet (r, i, j), with r an image, and (i, j) the position of P in r.

Algorithm 1 allows to find Q, the first pixel adjacent to a given pixel P in the west direction (algorithms for other directions are all similar). This is the classical algorithm for quadtree navigation (see for example [START_REF] Gargantini | An effective way to represent quadtrees[END_REF][START_REF] Samet | An overview of quadtrees, octrees, and related hierarchical data structures[END_REF]), slightly modified to deal with subdivision of any size n × m.

Algorithm 1: Move a pixel to west.

Input: P: a pixel in an hPFC, that does not belong to the first column of I 0 . Output: Q: the first pixel to the west of P. This algorithm has two main parts. First, we do a up-the-tree ascension while the current pixel Q belongs to the west of its image (lines 2-3). Since P does not belong to the first column of I 0 , we are sure that this loop finishes. Then we can move the pixel Q from one position in west direction (line 4).

The last main part (lines 5-10) is a classical down-the-tree descent, while the current pixel Q is a refined pixel. During this descent, we use either the same y coordinate than the pixel at the same depth than Q in the path going from I 0 to P (line 8), or we go to the first pixel in the east road of Q (line 9) if such a pixel does not exist (case when Q.depth > P.depth).

In order to be able to find the parent of pixel Q in constant time, without storing explicitly the parent link in the tree, we represent each pixel by a stack of triplets, describing the full path from I 0 to Q. Thanks to this stack, going up in the tree is done directly by a pop operation, and going down by a push of the new triplet.

Let us consider the hPFC given in Fig. 4 (that represents the irregular grid shown in Fig. 2) and show some examples of the move to west algorithm.

• starting from pixel a = ((I 0 , 1, 0)), we first move to west to ((I 0 , 0, 0)), then we go down twice to reach b = ((I 0 , 0, 0),(I 1 ,1,0),(I 5 ,1,0)) which is the first adjacent pixel to a;

• starting from pixel c = ((I 0 , 1, 2),(I 4 ,0,1),(I 9 ,0,1)), we first go up twice to reach ((I 0 , 1, 2)), then we move to west to d = ((I 0 , 0, 2)) which is the first (and unique) adjacent pixel to c;

• starting from pixel e = ((I 0 , 1, 2),(I 4 ,1,1),(I 10 ,0,1)), we first go up once to reach ((I 0 , 1, 2),(I 4 ,1,1)), then we move to west to ((I 0 , 1, 2),(I 4 ,0,1)), then we go down to f = ((I 0 , 1, 2),(I 4 ,1,1),(I 9 ,1,1)) which is the first (and unique) adjacent pixel to e. Note in this case that we go directly to pixel (1,1) in image I 9 because we use the y coordinate of the pixel at depth 2 in the path going from I 0 to P (which is P itself).

Crossings Identification on an hPFC

The second tool required to support the graph traversal operators is the identification of a crossing of a pixel in its adjacent pixel. More precisely, this is the crossing problem: Given the local ID of crossing X of a pixel P and the label of the road D that contains it, we want: (1) to identify the unique word pixel Q that also contains X and (2) to identify the local ID of X in Q.

This problem is discussed in [START_REF] Damiand | Rasterized planar face complex[END_REF], but only for cases where P and Q are word pixels of the same image.

There are three different cases:

1. P.depth = Q.depth: P and Q have the same number of crossings, and correspondence is straightforward (cf. [START_REF] Damiand | Rasterized planar face complex[END_REF], this is for example the case of pixels e and f in Fig. 2);

2. If P.depth > Q.depth: We count the number of crossings C in all the pixels adjacent to Q, after P (considering the order given by the orientation of road D in P), plus the number of crossings after X in P. The local ID of X in Q is the crossing number C in the road shared with P (counting the number of crossings is done by navigating through the pixel grids thanks to the operations introduced in the previous subsection). An example is given in Fig. 5. Let us suppose X is the unique crossing in the west road of pixel a (in blue in the figure). d is the pixel adjacent to the west of a. The number of crossings in the east road of d before a is 4 (crossings of pixels b and c). The local ID, in d, of X is 4 in the east road of d (thus the 5 th crossing of this road because IDs started from 0); 3. If P.depth < Q.depth: We compute N, the position of crossing X in its road, in reverse orientation. We move Q through the pixels adjacent to P, using the orientation of road D in Q. We stop this loop when the number of crossings, C, in road D is smaller or equal than N. In this case, we have found pixel Q which shares crossing X with P, and the local ID of X in Q is the crossing number N in the road shared with P. During the loop, when C>N, N is decreased by C. An example is given in Fig. 6. Here, X is the 5 th crossing in the west road of pixel a (in blue in the figure), and N = 6, the position of X in the west road starting from the north. b is the first pixel adjacent to the west of a. We iterate through pixels c, d, while decreasing N from 6 to 4, 2 then 1. Entering in pixel e, we have N = 1 < 3 the number of crossings in the east road of e, thus we stop the loop. The local ID, in e, of X is 1 in the east road of e (thus the 2 nd crossing of this road).

Computational Space and Time Complexity

Time Complexity

The time complexity of Algo. 1 that allows us to navigate between adjacent pixels in the west direction is linear in the depth of the tree (the same for algorithms for other roads). Indeed, this algorithm mainly traverses the tree, possibly a first time bottom-up and a second time top-down.

The complexity of the operator that computes the local crossing number of a given crossing in its adjacent pixel (given in the previous subsection), is linear in the depth of the tree plus the maximal number of crossings in the pixel (indeed, the time complexity of the navigation between crossings inside a same pixel is linear in the number or crossings in the pixel).

With these two operators, it is possible to define the previous, next and opposite operators by using the technique explained in [START_REF] Damiand | Rasterized planar face complex[END_REF] and summarized in Sect. 2.4. The idea is to follow the stabbings of an edge to a junction. The complexity of these operators is thus linear in the maximum number of stabbings in an edge times depth of the tree. @1 @2 @3 @4 @5 @6 @7 @8 @9 @10 0110 0001 0000 1000 0000 0000 0000 0000 0000 0000

1 0 1 0 0 0 1 0 0 1 0 0 ++=++ (+Y++)+ Figure 7
. Compact encoding of the hPFC shown in Fig. 4. In level 0, the code of each pixel starts with a bit to differentiate refined and word pixels (in green), and each pixel is encoded by the same number of bits (both refined and word pixels). In other levels, these mask-bits are grouped at the beginning of each image. Different record sizes are used to encode word and refined pixels.

We assume that the depth of the tree, the maximal number of crossings per pixel and the number of stabbings per edge are small comparing to the number of edges, and thus can be considered as constants. This assumption is verified in our experiments and makes sense in general, even though pathological models could be designed to invalidate it. With these assumptions, the complexity of the three operators previous, next and opposite is constant.

Computational Space Complexity

The space complexity of hPFC is given by the number of pixels times the (constant) storage size of one pixel. The size of a refined pixel is the size used to store the address of the refined image, while the size of a word pixel is given by the encoding of its word. This size is set to be large enough to contain the largest encoding of a pixel. It may be large if the graph contains a vertex with a huge valence (count of incident edges).

Compact Encoding of hPFC

In this section, we propose a compact encoding of hPFC. The key point of our solution is to use a small amount of memory space, while allowing to navigate through an hPFC by decoding locally the traversed pixel (we do not need to decode the full model, contrary to several compressed encoding).

Main Principle and Notation

rPFCs are encoded level by level of the tree. The first level, numbered 0, contains only I 0 , the root of the tree, while other levels, between 1 and d (for depth) can contain several rPFCs. For example, the hPFC of Fig. 4 has 3 levels.

We use our own made memory manager that can allocate an arbitrary number of contiguous bits in memory, segment it into contiguous k-bit records, and support writing and reading k-bit integers in and out of these records using the integer index of a record.

All images at depth l are stored in a contiguous pool. Let s l denote the size (measured in bits) of the memory pool allocated to store all of the images at depth l.

Let n l r denote the number of refined pixels at depth l, and let n l w denote the number of word pixels at depth l.

Let k l w denote the number of bits required to store words at level l, i.e. the number of bits to encode the longest word in this level, and let k l @ denote the number of bits required to store the rPFC addresses for refined pixels at level l.

Data Structure for the Root Image I 0

The information associated with each pixel of the root image is stored using a fixed representation that allocates the same number of bits to each pixel. For each pixel, P, the first bit indicates whether P is a word pixel or a refined pixel. The remaining max(k 0 @ , k 0 w ) bits encode either the address of the rPFC in level 1 (for refined pixels) or the word (for word pixels) (cf. example in Fig. 7).

Hence, we have a direct access to any pixel (i, j) and can compute its memory address in the pool from k 0 @ , k 0 w and the dimensions n 0 and m 0 of I 0 .

The main advantage of this encoding is the direct access to any pixel. Its main drawback is that some more memory space is wasted due to the use of the same number of bits for word and refined pixels. But this loss is negligible because the number of pixels in I 0 is usually very small comparing to the number of pixels in refined images (cf. Sect. 5).

Data Structure for Refined Images

The representation of the pool of a refined (non zero) level d stores the words of all its images in a contiguous (compacted) array.

The encoding of each refined image I starts by a bit mask (one bit per pixel of I) that identifies the type of each pixel (word or refined). Following this mask, each pixel P of I is encoded: the word of P for word pixel and the address of the refined image for refined pixels. Contrary to I 0 , we use here an irregular encoding: word pixels are encoded by k d w bits while refined pixels use k d @ (cf. example in Fig. 7). The key advantage of this encoding is to use less memory than the representation used for the root image (no memory space is lost to align the encoding of the two types of pixels). Its main drawback is the cost of computing the location of one pixel: we need to decode first the bit mask for all pixels before P in order to compute the address of P in the current image. But the overhead is small because the size of refined images is usually very small (often 2 × 2 in our experiments, cf. Sect. 5).

Computing Sizes and In-Pool-Locations of Arrays

The numbers n l r , of refined pixels per level, n l w , of word pixels per level, and k l w , of bits to store the maximal word per level can all be computed for each level of the tree, during the building or the loading of an hPFC.

The number of bits needed to store rPFC addresses per level, k l @ , and the size of each level, s l , can be computed recursively (top-down) by using n l r , n l w and k l w . First, there is no refined pixel at the last level of the tree, and thus n d r = 0 and k d @ = 0. We can thus directly compute

s d = n d w × k d w .
Then, each number of bits k l @ can be computed directly thanks to s l+1 . For each rPFC address at level l, we need to be able to address any position in the memory pool at level l + 1. Thus k l @ = log 2 (s l+1 ) . We can compute s l = n l w ×k l w +n l @ ×k l @ , and thus compute k l @ and s l for each level.

Illustration

Figure 7 shows an illustration of our compact encoding for the hPFC previously given in Fig. 4.

The first level has 12 pixels since the size of I 0 is 4 × 3; 4 of them being refined pixels having their encoding starting by 1.

The size of all refined images is 2 × 2, the encoding of each image starts by 4 bits that describe the type of all of its pixels.

The length of the maximal word in Level 2 is 7 and thus k 2 w = 16 (1 bit per '+', 4 bits per other symbol). The number of bits used to store Level 2 is thus s 2 = 6 × (4 × 16 + 4) = 408 bits (there are 6 images, each one with 4 pixels, each pixel uses 16 bits plus the 4 bits mask).

The length of the maximal word in Level 1 is 8 and thus k 1 w = 20. k 1 @ = 9 in order to address any position in Level 2. The number of bits used to store Level 1 is thus s 9 + 10 × 20 = 270 bits (there are 4 images, 6 refined pixels and 20 word pixels).

Similarly for Level 0, we have k 0 w = 20, k 0 @ = 9, s 0 = 12 × (20 + 1) = 252 bits (there are 12 pixels, each one encoded with 20 bits whatever their type, plus one bit for the type of the pixel).

Experiments

In [START_REF] Damiand | Rasterized planar face complex[END_REF], four different architectures were introduced to represent the words in a rPFC. The first version (V1) explicitly encodes the words using a fixed prefix code and puts them in a matrix with constant size entries. The other three versions concatenate the encoding of words in each row and store them in an array using different optimizations for the private pixels.

In this section we compare the hPFC against the first version of the rPFC using synthetic and real data sets. For the sake of simplicity, we set m 0 = n 0 and m r = n r .

Delaunay Triangulations

Our first data set consists of series of five random Delaunay triangulations with 500, 1,000, 2,000 and 4,000 vertices, from which we randomly removed 20% of the edges to increase irregularity. One of these triangulations is depicted in Fig. 8. More information about their number of cells is described in Table 1.

In our first experiment, we report results for different values for m 0 and m r . We have considered only the five Delaunay triangulations with 4,000 vertices and we have computed hPFCs with m 0 in {1,2,4,8,16,32,64,125,256} and m r in {2,4,8,16}. m 0 . Indeed, setting a large m 0 or m r produces unnecessary small pixels, many of which are empty. For these inputs, we conclude that m r must be 2 and m 0 must not be greater than 64.

Note that the oscillations in this graphic are produced by the combined effect of m 0 and m r on the size of the largest pixels. To illustrate this, consider the curve for m r = 16. The rPFCs with m 0 = 1 and m 1 = 16 are the same and thus occupy the same memory space. The memory space grows because when we refine the initial pixel grid, we multiply the number of pixels. On the other side, by refining further we end up avoiding later subdivision and thus the memory space decreases.

In our second experiment we process the twenty Delaunay triangulations and compare the memory space used by the rPFC and the hPFC. The grid size of the rPFCs is computed for each triangulation so that no edge fits inside a pixel. For the hPFCs, we fix (m 0 ,m r )=(64,2). The results are given in Fig. 10. For the 20 Delaunay triangulations, on average, hPFC is 1352, 1078, 4910 and 3620 times more compact than the rPFC (for 0.5K, 1K, 2K and 4K vertices, resp.). The reason for this is that a single short edge in a mesh forces the rPFC to use a large number of pixels. The rPFCs of the triangulations with 4,000 vertices contain in average more than 39.9 million pixels, 99.4% of them being private, while the hPFCs have only 10,292 pixels, and only 2.7% of them are private. Thus, even with the overhead for storing the hierarchical structure of the hPFC, it is 3,000 times more compact than the rPFC. In this experiment, the depth of the tree is between 5 and 8, the maximum length of words is between 24 and 30, and the average number of stabbings per edge is between 1.2 and 1.3.

To complement the previous experiment, we now compare the time complexity of the hPFC and the rPFC. For each triangulation, we traverse all its half-edges with a depth-breadth search. Figure 11 shows the results. It is clear that the hPFC is much faster than the rPFC. For the triangulations with 4,000 vertices, the rPFC traverses in average 151.03 pixels per edge of the PFC, while the hPFC traverses only 8.33 (word and refined) pixels.

We conclude that the hPFC is more efficient both in memory space and time than the rPFC for these triangulations. Note, however, that we were able to represent all these meshes with the rPFC, which is not the case for the following experiments.

GIS models

Our real data set consists of spatial databases of eight countries (Australia, Brazil, Canada, China, Germany, France, United Kingdom and Russia) available at http://www.gadm.org/. They are in shapefile format, a popular geospatial vector data format for geographic information system (GIS) softwares. The PFC of Australia is shown in Fig. 12.

Note that, unlike in the previous data set, these PFCs contain very short edges and thus the rPFC needs a huge number of pixels. Therefore, we have created five simplified versions of these meshes for comparing the hPFC with the rPFC. The simplified meshes are made so that they can fit in rPFCs with m 0 in {1,024, 2,048, 4,096, 8,192, 16,384}. Table 2 describes the average number of cells in each batch of simplified meshes. The original meshes, which contain on average 1,407,761 vertices, could only be represented with the hPFC. In this experiment, for the hPFCs that represent the original meshes, the depth of the tree is between 10 and 16, the maximum length of words is between 12 and 24, and the average number of stabbings per edge is between 0.3 and 0.5.

Figure 13 shows the average memory space and time for the simplified and the exact meshes. Regarding the memory space cost, the hPFC still outperforms the rPFC. The hPFC is 9, 17, 34, 75 and 188 times more compact than the rPFC (for the simplified meshes). For the exact (non simplified) GIS meshes, the rPFC would require 1,079 millions of pixels in average, while the hPFC only has 1.5 million (41.4% of private).

Regarding the time complexity, the roles are inverted. For the least simplified meshes, the rPFC traverses only 3.79 pixels per edge and the hPFC, 5.34. The GIS models have much shorter edges than the Delaunay triangulations, and thus, the overhead for navigating through the levels of the hPFC ends up making it 3.5 times slower.

Note that we could not compute the rPFCs for the original meshes because of the excessive number of pixels necessary.

We conclude that the hPFC is more compact than the rPFC, and it is necessary for representing real data sets with exact topology. Also, we conjecture that its time complexity depends on the ratio between the minimum and the mean length of the edges of the PFC.

Conclusion and Future Works

The hPFC hierarchical representation of Planar Face Complexes that is proposed here uses a fixed-resolution pixel-grid, quantizes the geometry of the vertices and edges to pixelresolution, and encodes graph connectivity by small number of bits per pixel. It improves the previously proposed rPFC, which requires that no edge be contained in a single pixel, by constructing a tree of rPFCs. This hierarchical, quadtree-like model, makes it possible to represent exactly any PFC topology and reduces storage cost significantly.

In future work, we plan to propose different compact representations for hPFC. We will explore different techniques that combine the principle of hPFC or of rPFC with more compact encodings or compression of the words, of the masks, and of the addresses. Indeed, there are numerous options for such improvements, and results of comparisons depend heavily on the nature of the PFC.

We would like to study if rPFC and hPFC can be used to represent non planar graphs, for example to describe a street network having bridges and tunnels. In some configurations, it maybe possible to represent a non-planar graph by a planar graph by inserting vertices, but support of traversal operators on such extensions is more challenging.

Lastly, we are working on the extension of rPFC and hPFC in 3D.

Figure 1 .

 1 Figure 1. A Planar Face Complex (PFC), that has 21 vertices, 24 edges and 5 faces. (The portion of the external face shown here is white.) The chosen pixel resolution is valid: It ensures that no edge lies entirely in a single pixel. The domain, which includes all but the external face, comprises 192 pixels, of which 126 are empty. The 79 crossings (where an edge enters a pixel) are shown as small yellow disks.

Figure 2 .

 2 Figure 2. The hPFC of the PFC shown in Fig. 1 uses only 42 pixels (6 being labeled). Only 17 of these are empty. It has only 30 crossings.

Figure 3 .

 3 Figure 3. A pixel having 11 crossings (yellow disks), and their corresponding symbols. The word describing the local topology inside this pixel is "=(Y+Y[(+)[+]])+".

Figure 4 .

 4 Figure 4. Tree representation of the hPFC of Fig. 2. The root image has 4 refined pixels. The tree has 3 depth layers and contains 42 word pixels (shown in blue) and 10 refined pixels. A "-" symbol is used to mark private pixels (which are contained in a finite face or in the external one).

1 Q 6 r 8 j 2 ←

 1682 ← P; 2 while Q is adjacent to the west border do 3 Q←parent(Q); 4 Q.i ←Q.i -1; 5 while Q is a refined pixel do ← address in pixel Q; 7 if Q.depth ≤ P.depth then y coordinate of pixel at depth Q.depth in the path going from I 0 to P; 9 else j 2 ← 0; 10 Q← (r 2 , n -1, j 2 ); 11 return Q;

Figure 5 .Figure 6 .

 56 Figure 5. Example of adjacent pixels with different depths: a.depth = 2 > d.depth = 0.

Figure 8 .

 8 Figure 8. One of the Delaunay triangulations used with 500 vertices, and a detail of it.

Figure 9 .

 9 Figure 9. Memory comparison between hPFCs with different parameters.

Figure 10 .

 10 Figure 10. Memory comparison between the rPFC and the hPFC for the synthetic data set.

Figure 11 .

 11 Figure 11. Time comparison between the rPFC and the hPFC for the synthetic data set.

Figure 12 .

 12 Figure 12. One of the GIS mesh used (Australia), and a detail of it.

Figure 13 .

 13 Figure 13. Memory and time comparison between the rPFC and the hPFC for the GIS data set.

Table 1 .

 1 Average number of edges and faces in each batch of Delaunay triangulations.

	# vertices	0.5K	1K	2K	4K
	# edges	1,205 2,405 4,805 9,605
	# faces	703 1,403 2,804 5,605

1 = 4 × 4 + 6 ×

Table 2 .

 2 Average number of vertices, edges and faces in each batch of simplified GIS models and the exact GIS models.

Acknowledgements

This project received funding from the European Unions Horizon 2020 Research and Innovation program under the Marie Sklodowska-Curie (grant 659526).