M Carlavan
email: carlavan@i3s.unice.fr

L Blanc-Féraud

M Antonini

C Thiebaut
email: carole.thiebaut@cnes.fr

C Latry

Y Bobichon
email: yves.bobichon@thalesaleniaspace.com

A SATELLITE IMAGING CHAIN BASED ON THE COMPRESSED SENSING TECHNIQUE

INTRODUCTION

Satellite imaging has been the focus of intense work in remote sensing for the last years. The ability of satellite optical systems to produce high-resolution images has indeed been of a great interest in applications such as change detection or image classification. It has however outcomed to be quite challenging for the design of satellite acquisition chains. Due to the increase of the imaging resolution and the conservation of the swath, images acquired from newly built satellites1 can actually reach a size of 30000 × 30000 pixels. Dealing with such a volume of data has important consequences on embedded resources, which require more memory, more computing capacity and therefore more powerful electrical sources.

In a classical satellite imaging system, the observed image is sampled at the Shannon frequency to give N pixels and then compressed by some coding algorithms, like the JPEG standards. The purpose of the coding step is to represent the image on a limited number K << N of coefficients to match the low capacity of the on-board mass storage. However, using an expensive scheme to sample the whole image for, finally, retaining only K coefficients may appear to be wasteful. Many on-board resources could then be saved up if the compressed coefficients were directly acquired out of the sensor.

Recently, a new theory of sampling has been emerged in the signal processing community. This theory, introduced as the Compressive Sampling or Compressed Sensing (CS) [START_REF] Candes | An introduction to compressive sampling[END_REF], suggests that one can reconstruct perfectly a signal, supposed to be sparse in some basis, from a limited (i.e. fewer than Shannon) number of incoherent measurements. Although the design of a sensor able to produce these measurements is difficult and beyond the scope of this paper, the CS technique clearly appears to be adapted to the satellite imaging chain. It could indeed drastically simplify the process of image acquisition by providing a reduced number of measurements, directly outcomed from the sensor, therefore saving an important quantity of resources. It is also valuable to point out that the CS framework provides an acquisition technique whose performances depend mainly on the reconstruction algorithm done on the ground. In comparison, the current acquisition imaging chain is bounded by the efficiency of the compression scheme embedded on-board. If one wants to increase the quality of the final image, one has to design a new satellite imaging chain. This "universal" coding feature [START_REF] Candès | Near-optimal signal recovery from random projections: Universal encoding strategies?[END_REF] of the CS is thus very attractive.

Though it is not a general result (see [START_REF] Goyal | Compressive sampling and lossy compression[END_REF] for example), previous works [START_REF] Schulz | On the empirical rate-distortion performance of compressive sensing[END_REF] have shown that the CS technique may be competitive regarding to a wavelet-based compression scheme on smoothed classical test images. But to the best of our knowledge, no works have been dedicated to this comparison for high-resolution satellite imaging, taking into account the degradations of the satellite imaging acquisition chain (blur, instrumental and quantizing noises). Our contribution is then to present a numerical exploration of the Compressed Sensing technique applied to satellite imaging. We present a satellite imaging chain based on this technique and we evaluate its efficiency with respect to the classical imaging chain composed of a Nyquist sampling followed by a wavelet-based compression scheme. Note that this work only focusses on the problem of CS image reconstruction in application to satellite imaging and does not deal with the problem of designing a sensor able to produce the measurements. The conception of such sensor is a very difficult task. Some results have been obtained on this aspect [START_REF] Duarte | Single-pixel imaging via compressive sampling[END_REF] but much work still needs to be addressed to provide efficient high-resolution CS sensors.

The rest of the paper is organized as follows. In section 2, we present the Compressed Sensing theory. We propose, in section 3, an algorithm adapted to the satellite imaging which reconstructs the image from the CS measurements. The originality of the proposed algorithm is that it takes into account the main degradations of the imaging chain such as blur, instrumental and quantizing noise. We compare, in section 4, the reconstruction results of this method to the ones obtained using the classical acquisition chain based on a wavelet transform. We show global rate-distortion and visual results on a real satellite data. We conclude in section 5 and present perspectives for future works.

COMPRESSIVE SAMPLING FRAMEWORK

Compressive Sampling (CS) is a new emergent sampling theory introduced independently by E.J Candès [START_REF] Candès | Decoding by linear programming[END_REF] and D.L. Donoho [START_REF] Donoho | For most large underdetermined systems of linear equations the minimal l 1 -norm solution is also the sparsest solution[END_REF]. The motivation behind the CS technique is to perform in the same time the acquisition and the compression of the signal. We give a quick overview of this technique in this section but more information can be found in the referred works.

In a classical satellite imaging system, the image is firstly sampled at the Nyquist frequency to give N pixels. However, it is well-known that natural images have a sparse representation in wavelet basis [START_REF] Antonini | Image coding using wavelet transform[END_REF] such that almost all the information of the image holds on a limited number K << N of coefficients. This property of sparsity is strongly used by image compression systems like the JPEG-2000 standard which projects the image on a wavelet basis such that only K coefficients become significant. But this technique does not seem to be adapted to low-resources applications as all the image is usually required to perform the compression. One needs then to design an acquisition chain which has sufficient resources to sample and store the N pixels of the image which, considering its size, is very expensive to achieve.

The main result of the Compressive Sampling theory states that the original signal can be recovered from a small number of measurements [START_REF] Candès | Decoding by linear programming[END_REF] directly outcomed from the sensor. The key of the CS theory relies on the supposed sparsity of the original signal, meaning that it can be perfectly represented in some basis Ψ with only K non-null coefficients. Based on this property of sparsity, the authors of [START_REF] Candès | Decoding by linear programming[END_REF] showed that only M (with M << N ) measurements are required to perfectly reconstruct the original signal with a high probability. These M observations are obtained by the projection of the true image x 0 ∈ R N on a measurement matrix Φ :

R N → R M y = Φx 0 . (1)
This formulation is actually very common and the particularity of the CS lies in the design of the measurement matrix. As mentioned by [START_REF] Candès | Decoding by linear programming[END_REF], the matrix Φ should spread the information of the image and most of random matrices perfectly meet this condition. More generally, a property named the Restricted Isometry Property (RIP) has been introduced in [START_REF] Candès | Near-optimal signal recovery from random projections: Universal encoding strategies?[END_REF] to ensure that the image may be recovered form the incoherent measurements. To simplify, the RIP requires the sub-matrices of the measurement matrix Φ to be orthogonal enough to be able to recover the K-sparse components of the image x 0 .

If the RIP condition is satisfied, the image x 0 can be recovered by minimizing the l 1 -norm of its coefficients under the constraint that its projection on Φ is equal to the observed vector y

[9] Find x = arg min Ψx 1 subject to x ∈ R N y = Φx . ( 2
)
The optimization problem (2) is a particular instance of the Basis Pursuit (BP) problem [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] which can be solved using classical algorithms from the linear programming literature. It can be interpreted as follows. The measurement matrix Φ is a matrix which spreads the image in the measurement vector y. Although it is not orthogonal, if Φ satisfies the RIP, then the inverse solution Φ † y contains all the information of the image x but in disorder. Also remind that the representation of the image x in the basis Ψ is sparse or, in other words, strongly compact. Minimizing the l 1 -norm of its coefficients will then put the non-null coefficients back at the correct position, recovering therefore the whole image.

The more sparse is x 0 , the easier it will be for the algorithm (2) to place these non-null coefficients. Clearly, recovering the image x 0 highly depends on the link between the compactness of the decomposition basis Ψ and the diffusion of the measurement matrix Φ. More generally, the algorithm (2) efficiently recovers the original image only if matrices Φ and Ψ are completely uncorrelated.

A mutual coherence (MC) µ has been introduced in [START_REF] Candès | Sparsity and incoherence in compressive sampling[END_REF] to measure this correlation and more precisely, to measure the correlation between each vector basis φ i and ψ i of Φ and Ψ. This coherence belongs to [1,

√ N ] [1]; a small value of µ meaning that the matrices Ψ and Φ are completely uncorrelated. For example, if Φ is the Fourier basis, then the minimal coherence is obtained with Ψ = I (the sampling operator) and is equal to 1. For our application, we will detail the choice of matrices Φ and Ψ in the numerical section of the paper. More generally, solving (2) recovers x 0 exactly if

M ≥ Cµ 2 (Φ, Ψ)K log(N ), C < 1 is a constant. ( 3
)
In classical imaging systems, acquired images are usually degraded by both blur and instrumental noise. As shown in [START_REF] Jianwei | Deblurring from highly incomplete measurements for remote sensing[END_REF], the CS technique can be extended to this scenario. The acquisition model becomes [START_REF] Jianwei | Deblurring from highly incomplete measurements for remote sensing[END_REF]

y = ΦHx 0 + z, (4)
where H : R N → R N is the blur matrix and z ∈ R M is an additive noise. In the classical case of an additive white Gaussian noise of variance σ 2 z , the reconstruction algorithm may write [START_REF] Jianwei | Deblurring from highly incomplete measurements for remote sensing[END_REF] Find x = arg min Ψx

1 subject to x ∈ R N y -ΦHx 2 2 ≤ M σ 2 z . ( 5
)
Similarly to (2), the optimization problem ( 5) is a particular instance of Basis Pursuit Denoising (BPDN) which can also be solved using linear programming techniques [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]. Based on the CS framework, we formulate in the next part the acquisition model in the case of the satellite chain and we propose an algorithm to reconstruct the image from the measurements vector.

APPLICATION TO SATELLITE IMAGING

As said previously, we assume that we have at our disposal a sensor able to produce incoherent measurements, in the sense of the CS framework. We are interested in evaluating the quality of the reconstructed image in comparison to the image obtained using the current acquisition chain based on wavelet compression [START_REF] Antonini | Image coding using wavelet transform[END_REF]. A classical satellite imaging chain is represented figure 1. It is mainly composed by three parts: The acquisition which captures the scene, the compression which reduces the weight of the measured image and the restoration which reconstructs the image from the measurements and reduces the degradations collected during the processing chain. We assume that the transmission does not introduce any degradations on the measured coefficients.

The acquisition process affects the quality of the true image by adding blur and instrumental noise. The blur H is mainly caused by the natural environment and the imperfection of the acquisition components. It is supposed to be known. The instrumental noise z is also the composition of several noise sources such as a photon noise and an electronic noise. It can be globally modeled as a zero-mean Gaussian distribution with a known variance σ 2 z .

In addition to blur and instrumental noise, the measurements are also degraded by quantizing noise. A satellite is not able to continuously transmit the acquired images as it sometimes flies over "shadow zones" where a transmission is not possible. It has to save the acquired images on its on-board mass storage to transmit them later. But the on-board storage capacity of a satellite is highly limited. The measurements are thus always coded to fit a target compression rate, which is also beneficial for fast transmission. This compression introduces an irreversible error known as the coding or quantizing noise. This quantization Q can be modeled as a scalar uniform quantization which quantizing step ∆ i depends on the coefficient y 0i , i ∈ {1, . . . , M } regarded

Q [y 0i ] = ∆ i y 0i ∆ i + 1 2 , (6)
where ⌊.⌋ is the floor function which returns the greatest integer less than or equal to its argument. The quantizing step ∆ i can be transmitted with the image as in the JPEG standard [START_REF] Wallace | The jpeg still picture compression standard[END_REF] or can be deduced during the decoding algorithm for more recent methods [START_REF] Taubman | High performance scalable image compression with ebcot[END_REF][START_REF] Said | A new, fast, and efficient image codec based on set partitioning in hierarchical trees[END_REF]. We assume in the following that the quantizing steps ∆ i are known. Let b = Q[y 0 ]y 0 be the quantizing error. From ( 6), we have for each coordinate b i of b

- ∆ i 2 ≤ b i < ∆ i 2 , ∀i ∈ {1, . . . , M } (7)
or equivalently

b ∈ B, with B := b ∈ R M , - ∆ i 2 ≤ b i < ∆ i 2 . ( 8
)
Using the previous definition of b, we propose to modelize the the observed measurements as

y = Q [ΦHx 0 + z] = ΦHx 0 + z + b, (9)
where y is the measurements vector. The extension of the reconstruction algorithm ( 5) to the acquistion model ( 9) is simple. First, simply remark that the problem (5) can also be written

Find x = arg min Ψx 1 subject to x ∈ R N , z ∈ R M z 2 2 ≤ M σ 2 z y = ΦHx + z . ( 10
)
The variable b needs to be add to the problem [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] to take into account the presence of the coding noise. Using ( 8) and ( 9), the reconstruction problem writes

Find x = arg min Ψx 1 subject to x ∈ R N , z ∈ R M , b ∈ R M z 2 2 ≤ M σ 2 z , b ∈ B, y = ΦHx + z + b . ( 11
)
The problem ( 11) can be further simplified by noting that the variable b can be replaced by y -(ΦHx + z). We finally propose to formulate the reconstruction problem as

Find x = arg min Ψx 1 subject to x ∈ R N , z ∈ R M z 2 2 ≤ M σ 2 z , y -(ΦHx + z) ∈ B . ( 12
)
The optimization problem ( 12) is a convex problem constrained on convex sets and thus admits a unique (convex) set of solutions [START_REF] Boyd | Convex Optimization[END_REF]. However, the presence of the linear operators Ψ, Φ and H make it difficult to solve. We propose here to use the Alternating Direction Method of Multipliers (ADMM) proposed in [START_REF] Afonso | An augmented lagrangian approach to the constrained optimization formulation of imaging inverse problems[END_REF]. The idea of this algorithm is to insert split variables such that the linear operators in [START_REF] Jianwei | Deblurring from highly incomplete measurements for remote sensing[END_REF] are moved into a linear constraint. Its presentation is however beyond the scope of this paper, and we refer the interested reader to [START_REF] Afonso | An augmented lagrangian approach to the constrained optimization formulation of imaging inverse problems[END_REF] (and references therein) for futher information on the numerical implementation of this algorithm.

NUMERICAL RESULTS

In this part, we evaluate the performances of the CS technique for satellite imaging in comparison to the chain based on a wavelet compression scheme. For the numerical experiments, we choose the measurement matrix Φ to be the noiselet transform [START_REF] Coifman | Noiselets[END_REF] and set Ψ to be the gradient operator such that ||Ψx | 1 is the Total Variation (TV) [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. We made the choice of the TV as it is almost equivalent to a Haar basis which, as required by the CS framework, shares a small mutual coherence with the noiselet transform [START_REF] Candes | An introduction to compressive sampling[END_REF]. Note that a TV prior is also well adapted to the regularity of high-resolution satellite images.

As mentioned previously, we compare the CS acquisition technique to the classical acquisition chain which consists in sampling the real image at the Shannon frequency followed by a compression scheme. The considered compression algorithm is close to the JPEG-2000 standard and uses the bi-orthogonal Cohen-Daubechies-Feauveau (CDF) 9/7 wavelet transform described in [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF] followed by the same quantization process as the one defined in [START_REF] Candès | Decoding by linear programming[END_REF]. As in the CS technique, we can design an algorithm to reconstruct the image from the noisy observed wavelet coefficients y

Find x = arg min Ψx 1 subject to x ∈ R N , z ∈ R N y -W (Hx + z) ∈ B z 2 2 ≤ N σ 2 z , ( 13
)
where W is the CDF 9/7 wavelet transform and Ψ is the gradient operator. Note that the formulation ( 13) is not expressed using any matrix Φ as, in this case, the measurement matrix is the sampling operator (Φ = I). We will compare the results of techniques ( 12) and ( 13) visually but also in a rate-distortion sense. As both techniques offer different ways to control the target coding rate, we now detail the choice of the chain in each case.

For the CS technique, we take benefit from the fact that the image can ideally be reconstructed from less measurements than Shannon. More precisely, for a low target rate, we will restrict the number of measurements M to be small and when the target rate is high, we will increase this number, the maximum number of measurements being equal to the number of pixels N . This particular choice comes from the fact that the distribution of the CS coefficients is quite large and that a high quantization has to be applied on these coefficients to reach low target rates [START_REF] Fletcher | On the rate-distortion performance of compressed sensing[END_REF]. It seems then more appropriate to tune the number of measurements M instead of tuning the quantizing steps, for a given coding rate. Consequently, we will always take ∆ i = 1, ∀i ∈ {1, . . . , M } for all coding rates. Note that these measurements will be taken randomly and that the position of the retained coefficients can be known at each side of the chain by transmitting the seed of the random generator.

The imaging chain based on a wavelet scheme does not however offer such feature. More precisely, all the coefficients have to be retained to be able to reconstruct the image. Consequently, for this technique, we will keep all the coefficients and we will tune the quantizing steps to reach the target coding rate. For more simplicity, we will take the same quantizing step for all coefficients ∆ i = ∆, ∀i ∈ {1, . . . , N }.

The reconstruction error will be evaluated using the Peak Signal-To-Noise Ratio (PSNR) defined, for 12 bits images, as

P SN R(x 0 , x) = 10 log 10 4095 2 1 N x 0 -x 2 2 . ( 14
)
where x 0 is the reference image and x is the reconstructed one. As mentioned previously, we will evaluate the results in a rate-distortion sense. The distortion will be evaluated using the reconstruction error [START_REF] Taubman | High performance scalable image compression with ebcot[END_REF]. For the evaluation of the coding rate, we assume that the quantized coefficients will be encoded using an entropy encoder. The coding rate can then be measured using the entropy (expressed in bits/symbol) of the coefficients y [START_REF] Shannon | A mathematical theory of communication[END_REF]. It is defined as

H(y) = - ∞ m=-∞ p y (m) log 2 (p y (m)), (15)
where p y (m) is the probability for a quantized coefficient to get the symbol m. Note that for the imaging chain based on the CS technique, we only retain M coefficients. The entropy of the quantized coefficients will be thus multiplied by the ratio between the number of measurements and the number of pixels for that case.

We simulate the two imaging chains on the reference image depicted figure 2. This image is a simulation of the PLEAIDES-HR imaging chain. The blur H used in this simulation is then the Point Spread Function (PSF) of the optics of the PLEAIDES-HR satellite which has been provided by the CNES. For the simulation, we set the standard deviation of the instrumental noise to be equal to 10, σ z = 10. Results are shown on figures 3 and 4. From the rate-distortion function displayed on figure 3, we see that the CS technique does not give competitive reconstruction results in comparison to the wavelet-based technique, and stands 5-6 dB below this technique, for all compression rates. Visually, the reconstructed images are not very good as well. We can see on figure 4 that the CS reconstruction algorithm overregularizes the solution and creates large patterns, therefore losing the details of the image. Although it seems clear that the CS is a good acquisition technique as it better spreads the information than a wavelet transform, it also appears that high-resolution satellites images are not sparse enough, in usual basis, to be used with this technique. Moreover, as said previously the CS coefficients have a large distribution (larger than wavelet coefficients) making their coding difficult to perform, even when one only retains a limited number of these coefficients. We have however strong thoughts that the CS could be an efficient acquisition strategy for satellite images as it has already shown interesting results in application where the image is naturally strongly sparse, such as in MRI application [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF]. Following this idea, an imaging chain based on the CS technique may be interesting for galaxy observation missions which naturally give sparse images, as in astronomy where the CS exhibits great performances [START_REF] Bobin | Compressed sensing in astronomy[END_REF]. But in the case of earth observation missions, the approximative sparsity of the images does not seem to be sufficient to make the CS technique competitive regarding to the classical wavelet approach.

CONCLUSION

In this paper, we have experimentaly studied the performances of the CS acquisition technique in application to satellite imaging. We have proposed an algorithm to take into account the degradations of the satellite imaging chain (blur, instrumental and quantizing noise) and showed reconstruction results, visually and in a rate-distortion sense, on a real satellite data. We performed a comparison of this method to the classical acquisition method based on a wavelet transform and we showed that the CS acquisition technique does not give competitive results for earth observation imaging but may be interesting for observation missions which give sparse images. There is also room for improvements by considering more properly the distribution of satellite images to enhance usual priors and quantizing strategies which need to better fit the characteristics of the CS coefficients distribution. This will be the focus of future works.

Fig. 1 .

 1 Fig. 1. Considered acquisition imaging chain.

Fig. 2 .

 2 Fig. 2. Reference image, Cannes harbour (1024 × 1024 pixels). The dynamic range of this image is 12 bits.

Fig. 3 .

 3 Fig. 3. Rate-distortion function for the two acquisition techniques. The dashed curve is the PSNR w.r.t. the compression rate for the CS acquisition technique while the solid curve is the PSNR w.r.t. the compression rate for the wavelet-based method.

such as the european PLEIADES-HR satellite.


Miguel Liroz-Gistau
email: miguel.lirozgistau@inria.fr

Reza Akbarinia
email: reza.akbarinia@inria.fr

Esther Pacitti
email: esther.pacitti@lirmm.fr

Fabio Porto
email: fporto@lncc.br

Patrick Valduriez
email: patrick.valduriez@inria.fr

Dynamic Workload-Based Partitioning for Large-Scale Databases

Applications with very large databases, where data items are continuously appended, are becoming more and more common. Thus, the development of efficient workload-based data partitioning is one of the main requirements to offer good performance to most of those applications that have complex access patterns, e.g. scientific applications. However, the existing workload-based approaches, which are executed in a static way, cannot be applied to very large databases. In this paper, we propose DynPart, a dynamic partitioning algorithm for continuously growing databases. DynPart efficiently adapts the data partitioning to the arrival of new data elements by taking into account the affinity of new data with queries and fragments. In contrast to existing static approaches, our approach offers a constant execution time, no matter the size of the database, while obtaining very good partitioning efficiency. We validated our solution through experimentation over real-world data; the results show its effectiveness.

Introduction

We are witnessing the proliferation of applications that have to deal with huge amounts of data. The major software companies, such as Google, Amazon, Microsoft or Facebook have adapted their architectures in order to support the enormous quantity of information that they have to manage. Scientific applications are also struggling with those kinds of scenarios and significant research efforts are directed to deal with it [START_REF] Ailamaki | Managing scientific data[END_REF]. An example of these applications is the management of astronomical catalogs; for instance those generated by the Dark Energy Survey (DES) [START_REF]The dark energy survey[END_REF] project on which we are collaborating. In this project, huge tables with billions of tuples and hundreds of attributes (corresponding to dimensions, mainly double precision real numbers) store the collected sky data. Data are appended to the catalog database as new observations are performed and the resulting database size is estimated to reach 100TB very soon. Scientists around the globe can query the database with queries that may contain a considerable number of attributes.

The volume of data that such applications hold poses important challenges for data management. In particular, efficient solutions are needed to partition and distribute the data in several servers. An efficient partitioning scheme would try to minimize the number of fragments accessed in the execution of a query, thus reducing the overhead associated to handle the distributed execution. Vertical partitioning solutions, such as column-oriented databases [START_REF] Stonebraker | C-store: a column-oriented DBMS[END_REF], may be useful for physical design on each node, but fail to provide an efficient distributed partitioning, in particular for applications with high dimensional queries, where joins would require data transfers between nodes. Traditional horizontal approaches, such as hashing or range-based partitioning [START_REF] Chang | Bigtable: a distributed storage system for structured data[END_REF][START_REF] Cooper | PNUTS: Yahoo!'s hosted data serving platform[END_REF], are unable to capture the complex access patterns of scientific applications, especially since they usually make use of complex relations over big sets of columns, and are hard to be predefined.

One solution is to use partitioning techniques based on the workload. Graphbased partitioning is an effective approach for that purpose [START_REF] Curino | Schism: a workload-driven approach to database replication and partitioning[END_REF]. A graph (or hypergraph) that represents the relations between queries and data elements is built and the problem is reduced to that of minimum k-way cut problem, for which several libraries are available. However, this method always requires to explore the entire graph in order to obtain the partitioning. This strategy works well for static applications, but scenarios where new data are inserted to the database continuously, which is the most common case for scientific computing, introduce an important problem. Each time a new set of data is appended, the partitioning should be redone from scratch, and as the size of the database grows, the execution time of such operation may become prohibitive.

In this paper, we are interested in dynamic partitioning of large databases that grow continuously. After modeling the problem of data partitioning in dynamic datasets, we propose a dynamic workload-based algorithm, called Dyn-Part, that efficiently adapts the partitioning to the arrival of new data elements. Our algorithm is designed based on a heuristic that we developed by taking into account the affinity of new data with queries and fragments. In contrast to the static workload-based algorithms, the execution time of our algorithm does not depend on the total size of the database, but only on that of the new data and this makes it appropriate for continuously growing databases. We validated our solution through experimentation over real-world data sets. The results show that it obtains high performance gains in terms of partitioning execution time compared to one of the most efficient static partitioning algorithms.

The remainder of this paper is organized as follows. In Section 2, we describe our assumptions and define formally the problem we address. In Section 3, we propose our solution for dynamic data partitioning. Section 4 reports on the results of our experimental validation and Section 5 concludes.

Problem Definition

In this section, we state the problem we are addressing and specify our assumptions. We start by defining the problem of static partitioning, and then extend it for a dynamic situation where the database can evolve over time.

Static Partitioning

The static partitioning is done over a set of data items and for a workload. Let D = {d 1 , ..., d n } be the set of data items. The workload consists of a set of queries W = {q 1 , ..., q m }. We use q(D) ⊆ D to denote the set of data items that a query q accesses when applied to the data set D. Given a data item d ∈ D, we say that it is compatible with a query q, denoted as comp(q, d), if d ∈ q(D). Queries are associated with a relative frequency f :

W → [0, 1], such that q∈W f (q) = 1.
Partitioning of a data set is defined as follows.

Definition 1. Partitioning of a data set D consists of dividing the data of D into a set of fragments, π(D) = {F 1 , ..., F p }, such that there is no intersection between the fragments and the union of all fragments is equal to D.

Let q(F ) denote the set of data items in fragment F that are compatible with q. Given a partitioning π(D), the set of relevant fragments of a query q, denoted as rel (q, π(D)), is the set of fragments that contain some data accessed by q, i.e. rel (q, π(D)) = {F ∈ π(D) : q(F ) = ∅}.

To avoid a high imbalance on the size of the fragments, we use an imbalance factor, denoted by ǫ. The size of the fragments at each time should satisfy the following condition:

|F | ≤ |D| |π(D)| (1 + ǫ).
In this paper, we are interested in minimizing the number of query accesses to fragments. Note that the minimum number of relevant fragments of a query q is minfr (q, π(

D)) = |q(D)| (|D| / |π(D)|)(1+ǫ)
. We define the efficiency of a partitioning for a workload based on its efficiency for queries. Let us first define the efficiency of a partitioning for a query as follows:

Definition 2. Given a query q, then the efficiency of a partitioning π(D) for q, denoted as eff (q, π(D)) is computed as:

eff (q, π(D)) = minfr (q, π(D)) |rel (q, π(D))| (1)
In the equation above, when the number of accessed fragments is equal to the minimum possible, i.e. minfr (q, π(D)), the efficiency is 1. Using eff (q, π(D)), we define the efficiency of a partitioning π(D) for a workload W as follows.

Definition 3. The efficiency of a partitioning π(D) for a workload W , denoted as eff

(W, π(D)), is eff (W, π(D)) = q∈W f (q) × eff (q, π(D)) (2)
Given a set of data items D and a workload W , the goal of static partitioning is to find a partitioning π(D) such that eff (W, π(D)) is maximized.

Dynamic Partitioning

Let us assume now that the data set D grows over time. For a given time t, we denote the set of data items of D at t as D(t).

During the application execution, there are some events, namely data insertions, by which new data items are inserted into D. Let T ev = (t 1 , . . . , t m ) be the sequence of time points corresponding to those events. In this paper, we assume that the workload is stable and neither the queries nor their frequencies change. However, the queries may access new data items as the data set grows.

Let us now define the problem of dynamic partitioning as follows. Let T ev = (t 1 , . . . , t m ) be the sequence of time points corresponding to data insertion events; D(t 1 ), . . . , D(t m ) be the set of data items at t 1 , . . . , t m respectively; and W be a given workload. Then, the goal is to find a set of partitionings π(D(t 1 )), ..., π(D(t m )) such that the sum of the efficiencies of the partitionings for W is maximized. In other words, our objective is as follows:

Objective:

Maximize t∈Tev q∈W (f (q) × eff (q, π(D(t))))
3 Affinity Based Dynamic Partitioning

In this section, we propose an algorithm, called DynPart, that deals with dynamic partitioning of data sets. It is based on a principle that we developed using the partitioning efficiency measure described in the previous section.

Principle

Let d be a new inserted data item. From the definition of partitioning efficiency, we infer that if we place d in a fragment F , then the total efficiency varies according to the following approximation 4 :

eff (W, π(D ∪ {d})) ≈ eff (W, π(D))- q:q(F )=∅∧comp(q,d) f (q) minfr (q, π(D)) |rel (q, π(D))| (|rel (q, π(D))| + 1) , (3)
Thus, the partitioning efficiency is reduced whenever there are queries that did not access F but after the insertion of d to F have to access it, thereby increasing the number of relevant fragments. The lower the number of those queries, the less the resulting loss of efficiency. Based on this idea, we define the affinity between the data d and fragment F , which we denote as

aff (d, F ) = q:q(F ) =∅∧comp(q,d) f (q) minfr (q, π(D)) |rel (q, π(D))| (|rel (q, π(D))| + 1) (4)
Using (4), we can develop a heuristic algorithm that places the new data items in the fragments based on the maximization of the affinity between the data items and the fragments.

Algorithm

Our DynPart algorithm takes a set of new data items D ′ as input and selects the best fragments to place them. For each new data item d ∈ D ′ , it proceeds as follows (see the pseudo-code in Algorithm 1). First, it finds the set of queries that are compatible with the data item. This can be done by executing the queries of W on D ′ or by comparing their predicates with every new data item. Then, for each compatible query q, DynPart finds the relevant fragments of q, and increases the fragments affinity by using the expression in [START_REF] Chang | Bigtable: a distributed storage system for structured data[END_REF]. Initially the affinity of fragments is set to zero.

Algorithm 1 DynPart algorithm

for each d ∈ D ′ do for each q : comp(q, d) do for each

F ∈ rel (q, π(D)) do if feasible(F ) then //aff (F ) is initialized to 0 aff (F ) ← aff (F ) + f (q) minfr (q,π(D)) |rel(q,π(D))|(|rel (q,π(D))|+1) if ∃F ∈ π(D) : aff (F ) > 0 then dests ← arg max F ∈π(D) aff (F ) else dests ← π(D) F dest ← select from arg min F ∈dests |F | move d to F dest and update metadata
After computing the affinity of the relevant fragments, DynPart has to choose the best fragment for d. Not all of the fragments satisfy the imbalance constraints, thus we must only consider those that do meet the restrictions. We define the function feasible(F ) to determine whether a fragment can hold more data items or not. Accordingly, DynPart selects from the set of feasible fragments the one with the highest affinity. If there are multiple fragments that have the highest affinity, then the smallest fragment is selected, in order to keep the partitioning as balanced as possible.

DynPart works over a set of new data items D ′ , instead of a single data item. This particularly reduces the amortized cost of finding the set of queries that are compatible with each of the inserted items. This is why, in practice, we wait until a given number of items have been inserted and then execute our algorithm for partitioning the new data.

Let comp avg be the average number of compatible queries per data item, and rel avg be the average number of relevant fragments per query. Then, the average execution time of the algorithm is O(comp avg × rel avg × |D ′ |), where |D ′ | is the number of new data items to be appended to the fragments. The complexity can be O(|W |×|π(D)|×|D ′ |) in the worst case. However, in practice, the averages are usually much smaller than the worst case values. The reason is that the queries usually access a small portion of the data, thus the average number of compatible queries is low. In any case, in order to reduce the number of queries, we may use a threshold on the frequency, so that only queries above that threshold are considered. In addition, the partitioning efficiency of our approach is good, so the average number of relevant fragments per query is low.

Experimental Evaluation

To validate our dynamic partitioning algorithm, we conducted a thorough experimental evaluation over real-world data. In Section 4.1, we describe our experimental setup. In Section 4.2, we report on the execution time of our algorithm and compare it with a well known static workload-based algorithm. In Section 4.3, we study the effect our heuristic on the partitioning efficiency.

Set-up

For our experimental evaluation we used the data from the Sloan Digital Sky Survey catalog, Data Realease 8 (DR8) [2], as it is being used in LIneA in Brazil 5 . It consists of a relational database with several observations for both stars and galaxies. We obtained a workload sample from the SDSS SkyServer SQL query log data, which stores the information about the real accesses performed by users. In total, the database comprises almost 350 million tuples, that take 1.2 TB of space. The workload consists of a total of 27000 queries.

All queries were executed on the database and the tuple ids accessed by each of them were recorded. Only tuples accessed by at least one query were considered. We simulated the insertions on the database by selecting a subset of the tuples as the initial state and appending the rest of the tuples in groups. We varied the following two parameters: 1) the number of tuples inserted to the database before each execution of our algorithm, |D ′ |; and 2) the number of fragments in which the database is partitioned, |π(D)|. On each of the experiments, the specific numbers are detailed. Throughout the experiments we chose an imbalance factor of 0.15. All experiments were executed in a 3.0 GHz Intel Core 2 Duo E8400, running Ubuntu 11.10 64 bits with 4GB of memory.

Execution time

In this section, we study the execution time of the DynPart algorithm (DP in the figure) and compare it with a static graph partitioning algorithm (SP). For the later, we use PaToH6 , an hyper-graph partitioner. Figure 1(a) shows the comparison of the execution time for 16 fragments. In the case of the dynamic algorithm, we executed the algorithm with two values for |D ′ |: 0.5 and 1 million tuples. Similar results are obtained for different values of |π(D)|. As the difference between execution times of the static and the dynamic algorithms is significant, we use a logarithmic scale for the y-axis in order to show the results. The results are only depicted until a database size of 20 million tuples, as the memory requirements for the static partitioning are bigger than the memory of our servers. The dynamic algorithm does not cause any problem as the memory footprint depends on |D ′ |, which is constant throughout the experiment.

As we can seen, execution time increases for the static algorithm as the size of the database increases, provided that the size of the graph increases accordingly.

For the DynPart algorithm, on the other hand, the execution time stays at the same level, as it is always executed for the same number of data items. We executed our algorithm for different sizes of D ′ . Figure 1(b) shows the average execution time of the DynPart algorithm as |D ′ | increases for different number of fragments. As expected, the execution time is linearly related to the number of tuples before execution. Also, the higher number of fragments, the higher the execution time. This increase is not linear since the number of relevant fragments does not increase at the same pace. In fact, the number of relevant fragments does not exceed 8 for |π(D)| = 256 and 16 for |π(D)| = 1024.

Partitioning Efficiency

One of the important issues to consider for the DynPart algorithm is how our heuristic algorithm affects the partitioning efficiency. We executed the DynPart algorithm as the database is fed with new data after an initial partitioning using the static approach. With |D ′ | = 1 M, Fig. 2 shows how the partitioning efficiency evolves as the database grows for different number of fragments. Similar results were obtained for other configurations of |D ′ |. The efficiency decreases as the database grows, as expected, but this reduction is very small. For example, in the worst case, |π(D)| = 256, the partitioning efficiency decreases 2.23 × 10 -3 in average for each 10 million new tuples.

Conclusions

In this paper, we proposed DynPart, a dynamic algorithm for partitioning continuously growing large databases. We modeled the partitioning problem for dynamic datasets and proposed a new heuristic to efficiently distribute new arriving data, based on its affinity with the different fragments in the application.

We validated our approach through implementation, and compared its execution time with that of a static graph-based partitioning approach. The results show that as the size of the database grows, the execution time of the static algorithm increases significantly, but that of our algorithm remains stable. They also show that although the DynPart algorithm is designed based on a heuristic approach, it does not degrade partitioning efficiency considerably.

The results of our experiments show that our dynamic partitioning strategy is able to efficiently deal with the data of our application. But, we believe that its use is not limited to this application, and it can be used for data partitioning in many other applications in which the data items are appended continuously.

Fig. 1 .

 1 Fig. 1. a) Execution times of the static (SP) and dynamic (DP) approaches as DB size increases (|π(D)| = 16), b) DynPart execution time vs. |D ′ |

Fig. 2 .

 2 Fig. 2. Comparison of partitioning efficiency as the size of the DB grows (|D ′ | = 100k)

Note that this approximation is an equality in all cases but when the increment in |q(D)| makes minfr (q, π(D)) to be increased by 1, which happens very rarely.

Data from the DES project is still unavailable, so we have used data from SDSS, which is a similar, previous project

http://bmi.osu.edu/ ~umit/software.html


Parikshit Dutta

Elise Arnaud

Emmanuel Prados

Mathieu Saujot

Benoit Lefévre

Sensitivity Analysis of a Large Scale Integrated Land Use and Transportation System

In this work we perform sensitivity analysis of an integrated land use and transportation model. The modeling system used is TRANUS and the region modeled is the Grenoble urban agglomeration. The total effect of the input land use parameters on the QoI is investigated, which is the mean squared error in land use assignments. It is observed that 3 amongst 100 uncertain input parameters have maximum contribution towards QoI uncertainty.

Introduction

In a context of climate change and dwindling fossil resources, urban sustainability has become a key policy issue. Given the complexity of modern urban areas, designing sustainable policies calls for more than sheer expert knowledge. This is especially true for transport or land use policies, because of the strong interplay between the land use and the transportation systems. In such a context, land use and transport integrated modeling offers invaluable analysis tools for planners working on transportation and urban projects. Unfortunately the lack of confidence in their results, is somewhat nowadays an obstacle to a large diffusion to these tools. In particular, it is well known that uncertainty is pervasive in land use and transportation models. Uncertainty can be of epistemic kind, arising from mis-measurement, data mis-specification, and poor sampling. Also it can be of the aleatoric kind, which exhibits itself during forward propagation of the model. Quantification of these uncertainties is important for accurate forecasting of land use and transportation demands. It is a difficult problem to quantify uncertainty because of the large number of variables involved, and uncertainty within them is poorly characterized. Hence, it is important to identify major sources of uncertainty and their contribution towards output space variability. These reasons assert the need of performing sensitivity analysis while modeling a large scale land use and transportation system.

Uncertainty evolution has been studied by researchers for various ILUTMs, applied to model different regions [START_REF] Pradhan | Uncertainty propagation in an integrated land use-transportation modeling framework: output variation via urbansim[END_REF][START_REF] Waddell | Modeling urban development for land use, transportation, and environmental planning[END_REF]. The role of uncertainty in decision making and use parameter estimation in this case has also been investigated, by the land use community [START_REF] Sevcíková | Uncertain benefits: Application of bayesian melding to the alaskan way viaduct in seattle[END_REF]. As far as sensitivity analysis is concerned, scientists have studied models involving few parameters in a small scale [START_REF] Krishnamurthy | Propagation of uncertainty in transportation land use models: Investigation of dram-empal and utpp predictions in austin, texas[END_REF]. But an intensive analysis, focusing on sensitivity analysis of large scale ILUTM, with several parameters as uncertain, has largely been ignored.

In this work we have performed sensitivity analysis during calibration of the land use part of TRANUS [START_REF] De La Barra | Tranus-j: putting large models into small computers[END_REF], which was used to model the city of Grenoble, France. TRANUS land use algorithm consists in the resolution of a system of around 20 nonlinear equations and inequalities. The solution of such a system represents an economic equilibrium between supply and demand. The Grenoble land use model consists of 22 sectors which are spread amongst 225 zones. The total number of uncertain parameters for the land use model is 100. Each uncertain parameter has been assessed upon their total effect [START_REF] Saltelli | Sensitivity analysis[END_REF] on the quantity of interest or QoI , which is the L 2 norm of the difference in calculated (X) and observed production (X obs ).

Q = X obs -X 2 (1)

Methodology

The 22 economic sectors in the Grenoble model consist of population classified into 7 categories, employment classified into 8 categories, and real estate classified into 7 categories which interact with each other; see Table 1 or [START_REF] Dutta | Uncertainty propagation and sensitivity analysis during calibration of an integrated land use and transport model[END_REF]. The model describes all the interactions via a set of parameters which are then increased by the number of interactions. In the current work, the parameters assumed to be uncertain are the maximum and minimum consumption amounts, as well as the elasticity parameters in the induced demand function; the penalization parameters and substitution elasticities in the substitute demand equation, and the utility level in the locational utility function. Referring to TRANUS mathematical description [START_REF] De La Barra | The mathematical and algorithmic structure of tranus[END_REF], the uncertain parameter vector is given by Ψ = [min mn , max mn , δ mn , ωmn , δ, θ n ] ⊤ , where, Ψ ∈ R 100 ; the index m and n covering all economics sectors. Due to interactions amongst themselves, these parameters are constrained. Sampling such a space, is difficult, especially when the dimensionality is high [START_REF] Hastie | High-dimensional problems: p n[END_REF]. Hence the sensitivity analysis problem was computationally hard. We performed a complexity analysis of the algorithm for land use part of TRANUS. Given N zones and M sectors, the total production has an order of complexity of O(N 4 M 2 ). If the convergence is achieved after Z iterations, then the overall complexity of the total production computation is O Z × N 4 M 2 . For a sensitivity analysis using P s samples for each parameter, each having sample size R s , the overall complexity for production is O P s R s × Z × N 4 M 2 for each parameter. The order then increases polynomially as the number of zones and sectors increase.

In this work, we have employed high performance computing, and use of parallel machines to perform sensitivity analysis. We have used symmetric multiprocessing [START_REF] Romano | Towards scalable parallelism in monte carlo particle transport codes using remote memory access[END_REF] to distribute TRANUS implementation in P s R s parallel machines. However, the scalability of the land use algorithm itself was not investigated in the current work. In future, our plan is to further parallelize the execution of TRANUS land use model.

Results

To carry out the sensitivity analysis, Monte Carlo method is used [START_REF] Doubilet | Probabilistic sensitivity analysis using monte carlo simulation. a practical approach[END_REF]. As little or no prior information about the input parameters are available, the joint density is assumed to be the product of independent Gaussian marginal densities. The Monte Carlo samples for sensitivity analysis are obtained using rejection sampling on the joint density [START_REF] Gilks | Markov chain Monte Carlo in practice[END_REF], taking into account the constraints amongst parameters. At first, total effect of the input parameters on QoI in eqn. ( 1) is studied. It is observed that, the minimum amount of land used for economic activities required by industry (min 1-13 ) has 81% total effect on the QoI, followed by elasticity parameter of below average income households due to cost of individual housing δ 8-11 at 14%. It is described by a pie chart in fig. 1. Rest of the parameters have negligible or lesser effect on QoI.

We observe similar behaviors as previously observed, if we consider the price and the production convergence indicators (see [START_REF] De La Barra | The mathematical and algorithmic structure of tranus[END_REF]) as QoI s instead of the L2 production error (1). For both cases only one parameter, i.e., δ 8-11 , is responsible for majority of the unexplained output variation, see fig. 2(a) and fig. 2(b).

An analysis, of the total effect when only one kind of parameter is uncertain; is also performed. It is observed that if only elasticity parameters or minimum or maximum consumption amounts are varied, δ 8-11 , min 1-13 and max 1-13 have the maximum total effect on the QoI , respectively. On the contrary, if the parameters of the demand substitution function, i.e. substitution elasticity (δ) or penalty (ω mn ) parameters and utility level parameter (θ n ), are assumed to be uncertain, then more than one parameter contribute significantly, towards QoI variance. Figure 3 show the plots.

As far as the computational time is concerned, one iteration of TRANUS calibration step required approximately 530 seconds. Considering 100 samples of sample size 100 each the total sensitivity analysis time, if performed serially would take 1472.22 hours for a single parameter. We used 100 parallel machines which decreased the processing time 100 times.

Conclusion

In this work, we have performed sensitivity analysis during the calibration of the TRANUS model applied to the city of Grenoble in France (only on the land use part). The uncertain parameters taken were parameters from demand and substitution functions. It was observed that the unexplained variance in the quantity of interest, mainly due to parameters min 1-13 , δ 8-11 and max 1-13 , in that order. In future, we plan to use this information to solve a inverse problem and perform parameter estimation during calibration of TRANUS model of Grenoble.

Fig. 2 .

 2 Fig. 2. Description of total effect when a) production and b) price convergence indicator are QoI s.

Fig. 3 .

 3 Fig. 3. Total effect when only certain type of variable is uncertain. Top row represents parameters of the elasticity function. Notice the dominance of few parameters towards QoI uncertainty. Bottom row represents parameters of the substitution function and utility level. The effects are distributed evenly amongst parameters.

Table 1 .

 1 Description of the sectors.

	Sector number	Description	Sector Type
	1	Industrial employment	Employment
	2	Public employment	Employment
	3	Offices, research& development employment Employment
	4	Retail employment	Employment
	5	School employment	Employment
	6	Household income (rich)	Income
	7	Household income (above average)	Income
	8	Household income (below average)	Income
	9	Household income (poor)	Income
	10	Students	Income
	11	Individual housing	Real estate
	12	Apartment Housing	Real estate
	13	Land for economic activities	Real estate
	14	Land for shops in urban area	Real estate
	15	Social housing	Real estate
	16	University	Employment
	17	Retail employment (less frequent)	Employment
	18	Supermarket employment	Employment
	19	Retired income (rich)	Income
	20	Retired income (poor)	Income
	21	Other commercial ventures	Employment
	22	Other commercial land use	Real estate

  Fig.1. Total effect of parameters on the QoI . Sector 8 and 11 refers to 'below average income households' and 'individual housing' and 1 and 13 refers to 'industrial employment' and 'land for economic activities'

				Rest
		4%	< 1%	Min consump. sectors 19 & 22 Demand elast. sectors 8 & 11
				Min. consump. sectors 1 & 13
		14%
				81%
	2%	Rest Demand elast. sectors 8 & 11		2%	< 1%	Rest Min. consump. sectors 19 & 22 Demand elast. sectors 8 & 11
	98%				97%
	(a) The production convergence indicator is the quantity	(b) The price convergence indicator is the quantity of
	of interest		interest


J.-Y Courtonne

J Alapetite

P.-Y Longaretti

D Dupré

E Arnaud

E Prados

Study of cereals flows at local scales: Examples in the Rhône-Alpes région, the Isère département and the SCOT de Grenoble

The purpose of this article is to put forward the role applied mathematics and computer science could play in the field of ecological accounting and particularly in that of material flow analysis. It is done based on a detailed study on modeling cereal flows at sub-national scales.

Material Flow Analysis at the scale of urban areas, applied to primary materials, half-products and finished products, as well as the understanding of supply chains (production, intermediate transformation, final consumption) are necessary first steps to tackle the analysis of environment impacts. In this perspective, the lack of data at local scale is a serious issue and it is often necessary to find ways of using data from larger geographical levels (département, région or country). Our study is based on the supply and use balances produced by Agreste 1 at the national level. We display the results using Sankey diagrams because they make it possible to have a compact and complete view of the flows (production, transformation, trade of primary or derived products, and final consumption), while also underlining incoherences.

The national diagram tells us the data we have to look for at sub-national levels. In France, the main use of cereal grains (exception made of exports to foreign countries which represent roughly half of the production) are:

• livestock feed (about 2/3 of interior uses),

• consumption of the food-processing industry: mills, production of breakfast cereals, pasta and couscous, malt, and about half of the starch industry,

• consumption of other industries: distilleries (bioethanol) and the other half of the starch industry.

The per capita final consumption calculated with Agreste and FranceAgriMer data is very close to the figures taken from Insee2 surveys, which is a way of validating the methodology. In this work, we did not look into the meat industry nor into the trade of meat products. This should be done in the future since a considerable amount of cereals are consumed indirectly through meat.

Our goal was to produce Sankey diagrams for cereals at the scales of the Rhône-Alpes région, the Isère département and the SCOT3 de Grenoble.

The informations we are looking for fall into 5 categories: primary production (1), intermediate consumption (2), stocks (3), trade (4), final consumption (5). For each of these categories, we found data sources (Agreste, FranceAgriMer, ProdCom and SitraM) that are however only available down to a given scale. To go below these scales we distributed the last available data among sub-territories using proxys. For example, we show the correlation at the région and département level between the CORINE Land Cover "non-irrigated arable land" category and the production of cereals (1). Thus we use this correlation at the level of the SCOT where production data does not exist.

(2a) The ProdCom database provides the country's production (usually in kg) of a large number of half products and finished products. Using conversion factors, it is possible to deduce the amount of primary material used by each industry. We distributed the national intermediate consumptions of these industries based on the number of employees in each region. We finally looked for the precise location of the factories to distribute the regional intermediate consumptions among sub-territories.

(2b)

Regional consumption of cereals for livestock feed are not known, mainly because the cereals are distributed in various forms: self-consumption at the farm, grains bought by the farmer, or consumption of industrial compound products. We estimated the total livestock consumption of cereals for each region based on the data on livestock and slaughter and on the nutritional needs of the animals. Below the level of the département, we used the data from the Agreste Recensement 2000, although it is incomplete due to some confidential figures.

(3) Also because of the confidentiality of the data, information about the stock is partial. Averaging the figures throughout several years solves this issue since the variations of stock tend to compensate one year after another. In this work we therefore studied the period 2001-2009.

(4) For the trade data, we used a revised version of the SitraM database4 in order to avoid double counting between international and national trade.

(5) The per capita consumption was considered to be equal all over France. For instance, given that bread is mostly locally produced, we can deduce the amount of flour needed to sustain the final consumption of the territory. This is not true however for some other products, like glucose which is present in hundreds of derived products.

At the end, the more detailed the scale the more we need to use estimation technics. Moreover, we fill the data starting from the edges of the supply chain (production and final consumption), going to the center (intermediate transformation).

Our results for Rhône-Alpes and Isère seem globally consistent. There are however several incoherences which seem to point in the direction of imperfect trade data. At the scale of the SCOT de Grenoble, the smallest one studied here, it is interesting to visualize the separation between production and final consumption: the territory produces a lot of cereals but has no mill to convert them into flour. We therefore deduced it is imported. At any scale, self-sufficiancy is an interesting indicator to look at, even if this is not necessarily an objective for the territory. This concept can have a weak meaning (the territory produces more cereals than its inhabitants consume) or a strong meaning (the territory has access to the necessary production and transformation capacities to provide its inhabitants with every derived product).

We then developed a methodology to assess empirically the order of magnitude of the TRM survey5 . This statistical survey about transport of goods by road is the main part of SitraM national data concerning cereals. In order to gain in reliability, we based the study on every French région and we aggregated the different kinds of cereals. We then calculated the regional final consumptions of flour, bread, biscuits, breakfast cereals, pasta and couscous (all together and expressed in their equivalent weight of grain) in two different ways. For the first calculation, we used production, transformation and trade data whereas we based the second calculation on the national per-capita average. The case of cereals is favorable because most of the trade concerns primary products: the inevitable uncertainties surrounding transformation factors do not play a big role. Assuming both results should be equal (same cereal diet all over France) and that the less robust data is the inter-regional trade by road, we made an estimation of the margin of error of the TRM data on cereals. Some régions are indeed a lot over or a lot under the expected consumption (with even a few negative consumptions). Depending on the région, an error of 20% and up to 210% in the TRM data would explain the results. This is however still to be confirmed by a theoretical study of the TRM survey. Whatever the conclusion, the SitraM database could still be utilized, as a minimum, to assess if a resource is rather local or non-local, or even how many kilometers it travels before arriving to its destination.

The methodology discussed in this article, although still perfectible, aims to be applicable to other sectors of the economy and to other territories without having to undertake time-consuming and expensive survey-based local studies. Applied mathematics and computer science could help improving the results in two ways. The first improvement could be provided by data mining, automatic search for correlations and inconsistencies in the data sources. The second help would be in the field of equation solving: reconstructing supply chains is like having to solve many sub-systems of equations (for each product: supply = use) that are linked together with some equalities (the transformation factors between primary and derived products). Applied mathematics could help to formulate and solve this problem in a more efficient and robust way than it is done currently.

Finally, once the flows in terms of quantities have been established, they can be translated into multiple environmental pressures such as the use of arable land, nitrogen, water or energy. We applied this with the example of the Water Footprint.

Institut National de la Statistique et des Etudes Economiques.

Schéma de cohérence territoriale : the SCoT is an urbanplanning document dedicated to a group of towns or urban areas. By extension, we use the term SCOT to refer to the territory.

Système d'Information sur le TRAnsport de Marchandises.

Transport Routier de Marchandises.


Miguel Liroz-Gistau
email: miguel.liroz_gistau@inria.fr

Reza Akbarinia
email: reza.akbarinia@inria.fr

Esther Pacitti
email: esther.pacitti@lirmm.fr

Fabio Porto
email: fporto@lncc.br

Patrick Valduriez
email: patrick.valduriez@inria.fr

DynPart: Dynamic Partitioning for Large-Scale Databases

Keywords: Bases de données distribuées, partitionement dynamique, partitionnement basé sur workload 1

Il y a de plus en plus des applications scientifiques avec de très grandes bases de données distribuées où des nouvelles données sont ajoutées à la base de données en permanence. Pour offrir une bonne performance à la plupart de ces applications qui ont normalement des schémas d'accès complexes, on a besoin de développer des méthodes efficaces pour le partitionnement de données basé sur le workload. Toutefois, les approches existant basées sur le workload, qui sont exécutées d'une manière statique, ne peuvent pas être appliquée aux bases de données très volumineuses et dynamiques. Dans cet article, nous proposons DynPart, un algorithme dynamique de partitionnement pour les bases de données en croissance permanente. DynPart s'adapte efficacement à l'arrivée de nouvelles données en prenant en compte l'affinité de ces données avec les requêtes et les fragments. Contrairement aux approches antérieures, notre approche offre un temps d'exécution constant, peu importe la taille de la base de données, tout en obtenant une très bonne efficacité de partitionnement. Nous avons validé notre solution par l'expérimentation sur des données réelles ; les résultats montrent sa bonne performance.

Introduction

We are witnessing the proliferation of applications that have to deal with huge amounts of data. The major software companies, such as Google, Amazon, Microsoft or Facebook have adapted their architectures in order to support the enormous quantity of information that they have to manage. Scientific applications are also struggling with those kinds of scenarios and significant research efforts are directed to deal with it [START_REF] Ailamaki | Managing scientific data[END_REF]. An example of these applications is the management of astronomical catalogs; for instance those generated by the Dark Energy Survey (DES) [START_REF]The Dark Energy Survey[END_REF] project on which we are collaborating. In this project, huge tables with billions of tuples and hundreds of attributes (corresponding to dimensions, mainly double precision real numbers) store the collected sky data. Data are appended to the catalog database as new observations are performed and the resulting database size is estimated to reach 100TB very soon. Scientists around the globe can query the database with queries that may contain a considerable number of attributes.

The volume of data that such applications hold poses important challenges for data management. In particular, efficient solutions are needed to partition and distribute the data in several servers. An efficient partitioning scheme would try to minimize the number of fragments accessed in the execution of a query, thus reducing the overhead associated to handle the distributed execution. Vertical partitioning solutions, such as column-oriented databases [START_REF] Stonebraker | C-store: a column-oriented DBMS[END_REF], may be useful for physical design on each node, but fail to provide an efficient distributed partitioning, in particular for applications with high dimensional queries, where joins would require data transfers between nodes. Traditional horizontal approaches, such as hashing or range-based partitioning [START_REF] Chang | Bigtable: a distributed storage system for structured data[END_REF][START_REF] Cooper | PNUTS: Yahoo!'s hosted data serving platform[END_REF], are unable to capture the complex access patterns of scientific applications, especially since they usually make use of complex relations over big sets of columns, and are hard to be predefined.

One solution is to use partitioning techniques based on the workload. Graphbased partitioning is an effective approach for that purpose [START_REF] Curino | Schism: a workload-driven approach to database replication and partitioning[END_REF]. A graph (or hypergraph) that represents the relations between queries and data elements is built and the problem is reduced to that of minimum k-way cut problem, for which several libraries are available. However, this method always requires to explore the entire graph in order to obtain the partitioning. This strategy works well for static applications, but scenarios where new data are inserted to the database continuously, which is the most common case for scientific computing, introduce an important problem. Each time a new set of data is appended, the partitioning should be redone from scratch, and as the size of the database grows, the execution time of such operation may become prohibitive.

In this paper, we are interested in dynamic partitioning of large databases that grow continuously. After modeling the problem of data partitioning in dynamic datasets, we propose a dynamic workload-based algorithm, called DynPart, that efficiently adapts the partitioning to the arrival of new data elements. Our algorithm is designed based on a heuristic that we developed by taking into account the affinity of new data with queries and fragments. In contrast to the static workloadbased algorithms, the execution time of our algorithm does not depend on the total size of the database, but only on that of the new data and this makes it appropriate for continuously growing databases. We validated our solution through experimentation over real-world data sets. The results show that it obtains high performance gains in terms of partitioning execution time compared to one of the most efficient static partitioning algorithms.

The remainder of this paper is organized as follows. In Section 2, we describe our assumptions and define formally the problem we address. In Section 3, we propose our solution for dynamic data partitioning. Section 4 reports on the results of our experimental validation and Section 5 concludes.

Problem Definition

In this section, we state the problem we are addressing and specify our assumptions. We start by defining the problem of static partitioning, and then extend it for a dynamic situation where the database can evolve over time.

Static Partitioning

The static partitioning is done over a set of data items and for a workload. Let D = {d 1 , ..., d n } be the set of data items. The workload consists of a set of queries W = {q 1 , ..., q m }. We use q(D) ⊆ D to denote the set of data items that a query q accesses when applied to the data set D. Given a data item d ∈ D, we say that it is compatible with a query q, denoted as comp(q, d), if d ∈ q(D). Queries are associated with a relative frequency f :

W → [0, 1], such that q∈W f (q) = 1.
Partitioning of a data set is defined as follows.

Definition 1. Partitioning of a data set D consists of dividing the data of D into a set of fragments, π(D) = {F 1 , ..., F p }, such that there is no intersection between the fragments and the union of all fragments is equal to D.

Let q(F ) denote the set of data items in fragment F that are compatible with q. Given a partitioning π(D), the set of relevant fragments of a query q, denoted as rel (q, π(D)), is the set of fragments that contain some data accessed by q, i.e. rel (q, π(D)) = {F ∈ π(D) : q(F ) = ∅}.

To avoid a high imbalance on the fragments' load, we use an imbalance factor, denoted by ǫ. We define the load of a dataset D as L(D) = q∈W f (q)|q(D)|. Then, the load of each fragment must satisfy:

L(F ) ≤ L(D) |π(D)| (1 + ǫ).
In this paper, we are interested in minimizing the number of query accesses to fragments. Note that the minimum number of relevant fragments of a query q is minfr(q, π(D)) = L(q(D))

(L(D) / |π(D)|)(1+ǫ)
. We define the efficiency of a partitioning for a workload based on its efficiency for queries. Let us first define the efficiency of a partitioning for a query as follows:

Definition 2. Given a query q, then the efficiency of a partitioning π(D) for q, denoted as eff (q, π(D)) is computed as:

eff (q, π(D)) = minfr(q, π(D)) |rel (q, π(D))| (1)
In the equation above, when the number of accessed fragments is equal to the minimum possible, i.e. minfr(q, π(D)), the efficiency is 1. Using eff (q, π(D)), we define the efficiency of a partitioning π(D) for a workload W as follows.

Definition 3. The efficiency of a partitioning π(D) for a workload W , denoted as eff (W, π(D)), is

eff (W, π(D)) = q∈W f (q) × eff (q, π(D)) (2)
Given a set of data items D and a workload W , the goal of static partitioning is to find a partitioning π(D) such that eff (W, π(D)) is maximized.

Dynamic Partitioning

Let us assume now that the data set D grows over time. For a given time we denote the set of data items of D at t as D(t).

During the application execution, there are some events, namely data insertions, by which new data items are inserted into D. Let T ev = (t 1 , . . . , t m ) be the sequence of time points corresponding to those events. In this paper, we assume that the workload is stable and neither the queries nor their frequencies change. However, the queries may access new data items as the data set grows.

Let us now define the problem of dynamic partitioning as follows. Let T ev = (t 1 , . . . , t m ) be the sequence of time points corresponding to data insertion events; D(t 1 ), . . . , D(t m ) be the set of data items at t 1 , . . . , t m respectively; and W be a given workload. Then, the goal is to find a set of partitionings π(D(t 1 )), ..., π(D(t m )) such that the sum of the efficiencies of the partitionings for W is maximized. In other words, our objective is as follows:

Objective: Maximize t∈Tev q∈W (f (q) × eff (q, π(D(t))))

Affinity Based Dynamic Partitioning

In this section, we propose an algorithm, called DynPart, that deals with dynamic partitioning of data sets. It is based on a principle that we developed using the partitioning efficiency measure described in the previous section.

Principle

Let d be a new inserted data item. From the definition of partitioning efficiency, we infer that if we place d in a fragment F , then the total efficiency varies according to the following approximation 1 :

eff (W, π(D ∪ {d})) ≈ eff (W, π(D))- q:q(F )=∅∧comp(q,d) f (q) minfr(q, π(D)) |rel (q, π(D))| (|rel (q, π(D))| + 1) , (3)
Thus, the partitioning efficiency is reduced whenever there are queries that did not access F but after the insertion of d to F have to access it, thereby increasing the number of relevant fragments. The lower the number of those queries, the less the resulting loss of efficiency. Based on this idea, we define the affinity between the data d and fragment F , which we denote as

aff (d, F ) = q:q(F ) =∅∧comp(q,d) f (q) minfr(q, π(D)) |rel (q, π(D))| (|rel (q, π(D))| + 1) (4)
Using Equation 4, we develop a heuristic algorithm that assigns each new data item to the fragment which maximizes the affinity measure.

Algorithm

Our DynPart algorithm takes a set of new data items D ′ as input and selects the best fragments to place them. For each new data item d ∈ D ′ , it proceeds as follows (see the pseudo-code in Algorithm 1). First, it finds the set of queries that are compatible with the data item. This can be done by executing the queries of W on D ′ or by comparing their predicates with every new data item. Then, for each compatible query q, DynPart finds the relevant fragments of q, and increases the fragments affinity by using the expression in Equation 4. Initially the affinity of fragments is set to zero.

Algorithm 1 DynPart algorithm

for each d ∈ D ′ do for each q : comp(q, d) do for each

F ∈ rel (q, π(D)) do if feasible(F ) then //aff (F ) is initialized to 0 aff (F ) ← aff (F ) + f (q) minfr (q,π(D)) |rel (q,π(D))|(|rel (q,π(D))|+1) if ∃F ∈ π(D) : aff (F ) > 0 then dests ← arg max F ∈π(D) aff (F ) else dests ← {F ∈ π(D) : feasible(F )} F dest ← select from arg min F ∈dests L(F ) move d to F dest and

update metadata

After computing the affinity of the relevant fragments, DynPart has to choose the best fragment for d. Not all of the fragments satisfy the imbalance constraints, thus we must only consider those that do meet the restrictions. We define the function feasible(F ) to determine whether a fragment can hold more data items or not. Accordingly, DynPart selects from the set of feasible fragments the one with the highest affinity. If there are multiple fragments that have the highest affinity, then the least loaded fragment is selected, in order to keep the partitioning as balanced as possible.

DynPart works over a set of new data items ′ , instead of a single data item. This particularly reduces the amortized cost of finding the set of queries that are compatible with each of the inserted items. This is why, in practice, we wait until a given number of items have been inserted and then execute our algorithm for partitioning the new data.

Let comp avg be the average number of compatible queries per data item, and rel avg be the average number of relevant fragments per query. Then, the average execution time of the algorithm is O(comp avg × rel avg × |D ′ |), where |D ′ | is the number of new data items to be appended to the fragments. The complexity can be O(|W | × |π(D)| × |D ′ |) in the worst case. However, in practice, the averages are usually much smaller than the worst case values. The reason is that the queries usually access a small portion of the data, thus the average number of compatible queries is low. In any case, in order to reduce the number of queries, we may use a threshold on the frequency, so that only queries above that threshold are considered. In addition, the partitioning efficiency of our approach is good, so the average number of relevant fragments per query is low.

Experimental Evaluation

To validate our dynamic partitioning algorithm, we conducted a thorough experimental evaluation over real-world data. In Section 4.1, we describe our experimental setup. In Section 4.2, we report on the execution time of our algorithm and compare it with a well known static workload-based algorithm. In Section 4.3, we study the effect our heuristic on the partitioning efficiency.

Set-up

For our experimental evaluation we used the data from the Sloan Digital Sky Survey catalog, Data Realease 8 (DR8) [START_REF]Sloan Digital Sky Survey[END_REF], as it is being used in LIneA in Brazil 2 . It consists of a relational database with several observations for both stars and galaxies. We obtained a workload sample from the SDSS SkyServer SQL query log data, which stores the information about the real accesses performed by users. In total, the database comprises almost 350 million tuples, that take 1.2 TB of space. The workload consists of a total of 27000 queries.

All queries were executed on the database and the tuple ids accessed by each of them were recorded. Only tuples accessed by at least one query were considered. We simulated the insertions on the database by selecting a subset of the tuples as the initial state and appending the rest of the tuples in groups. We varied the

Execution time

In this section, we study the execution time of the DynPart algorithm (DP in the figure) and compare it with a static graph partitioning algorithm (SP). For the later, we use PaToH3 , an hyper-graph partitioner. Figure 1(a) shows the comparison of the execution time for 16 fragments. In the case of the dynamic algorithm, we executed the algorithm with two values for |D ′ |: 0.5 and 1 million tuples. Similar results are obtained for different values of |π(D)|. As the difference between execution times of the static and the dynamic algorithms is significant, we use a logarithmic scale for the y-axis in order to show the results. The results are only depicted until a database size of 20 million tuples, as the memory requirements for the static partitioning are bigger than the memory of our servers. The dynamic algorithm does not cause any problem as the memory footprint depends on |D ′ |, which is constant throughout the experiment.

As we can seen, execution time increases for the static algorithm as the size of the database increases, provided that the size of the graph increases accordingly. For the DynPart algorithm, on the other hand, the execution time stays at the same level, as it is always executed for the same number of data items.

We executed our algorithm for different sizes of D ′ . Figure 1(b) shows the average execution time of the DynPart algorithm as |D ′ | increases for different number of fragments. As expected, the execution time is linearly related to the number of tuples before execution. Also, the higher number of fragments, the higher the execution time. This increase is not linear since the number of relevant fragments does not increase at the same pace. In fact, the number of relevant fragments does not exceed 8 for |π(D)| = 256 and 16 for |π(D)| = 1024.

Partitioning Efficiency

One of the important issues to consider for the DynPart algorithm is how our heuristic algorithm affects the partitioning efficiency. We executed the DynPart algorithm as the database is fed with new data after an initial partitioning using the static approach. With |D ′ | = 1 M, Fig. 2 shows how the partitioning efficiency evolves as the database grows for different number of fragments. Similar results were obtained for other configurations of |D ′ |. The efficiency decreases as the database grows, as expected, but this reduction is very small. For example, in the worst case, |π(D)| = 256, the partitioning efficiency decreases 2.19 × 10 -3 in average for each 10 million new tuples.

Conclusions

In this paper, we proposed DynPart, a dynamic algorithm for partitioning continuously growing large databases. We modeled the partitioning problem for dynamic datasets and proposed a new heuristic to efficiently distribute new arriving data, based on its affinity with the different fragments in the application.

We validated our approach through implementation, and compared its execution time with that of a static graph-based partitioning approach. The results show that as the size of the database grows, the execution time of the static algorithm increases significantly, but that of our algorithm remains stable. They also demonstrate that although the DynPart algorithm is designed based on a heuristic approach, it does not degrade partitioning efficiency considerably.

Based on the results we can state that our dynamic partitioning strategy is able to efficiently deal with the data of our application. But, we believe that its use is not limited to this application, and it can be used for data partitioning in many other applications in which the data items are appended continuously.

Figure 1 :

 1 Figure 1: a) Execution times of the static (SP) and dynamic (DP) approaches as DB size increases (|π(D)| = 16), b) DynPart execution time vs. |D ′ |

  |π(D)| = 4 DP, |π(D)| = 16 DP, |π(D)| = 32 DP, |π(D)| = 256

Figure 2 :

 2 Figure 2: Comparison of partitioning efficiency as the size of the DB grows (|D ′ | = 1M)

Note that this approximation is an equality in all cases but when the increment in |q(D)| makes minfr (q, π(D)) to be increased by 1, which happens very rarely.

Data from the DES project is still unavailable, so we have used data from SDSS, which is a similar, previous project

http://bmi.osu.edu/~umit/software.html