Matthieu Rolland
email: matthieu.rolland@ifpen.fr

Cláudio P Fonte

Incertitude induced by testing a small number of catalytic pellets in fixed beds

Keywords: Fixed beds, downscaling, sampling, incertitude, criteria, catalyst testing Highlights

Heterogeneous catalyst pellets evaluation is performed in smaller and smaller reactors with the main objective of reducing catalyst development costs. Downscaling raises the question: is there a lower limit to the number of pellets that ought to be tested in a reactor so that the results do not depend on which pellets are sampled? Major sources of variability among a catalyst pellet population are dimensions and many parameters that can be grouped as "intrinsic kinetic activity" (porosity, tortuosity, active phase repartition and availability).

In this paper, we present a methodology to estimate the incertitude induced by variability on size and kinetics on the evaluation of the apparent kinetic constant in a fixed bed reactor. Analytical expressions are presented to predict this incertitude for sphere and cylinder shaped pellets, with a first order kinetic law in two limiting cases: absence of mass transfer limitations and presence of severe internal mass transfer limitations. The predicted incertitude scales as the square root of sample population: downsizing increases the incertitude. We propose criteria to evaluate the minimum number of pellets to sample, so that sampling induced variability is lower than an acceptable incertitude, expressed in °C. This acceptable incertitude could be for example a fraction of the experimental temperature incertitude.

• Methodology to evaluate the consequences of testing small catalyst samples in fixed beds

• Explicit equation to evaluate the incertitude for spherical and cylindrical pellets

• Criterion on minimum number of pellets required for a low variability catalytic test Heterogeneous catalyst pellets evaluation is performed in smaller and smaller reactors with numerous benefits such as reduction of the amount of feedstock, waste products and catalyst, reduction of investment cost per reactor, etc… Data acquired in the test is used for kinetic parameters estimation in order to select the best catalyst and to predict its performance in an industrial size reactor. Here we focus on the case of millimetric size pellets, used for example in refining processes, whose performance is evaluated in fixed beds. In some High-Throughput Experiment catalyst testing, the number of catalyst pellets can be as low as one hundred, while catalyst pellets are never exactly identical with variability in size, structure, impregnation quality… This raises a quite crucial question: is there a lower limit to the number of pellets that ought to be tested in the reactor so that the results do not depend on the laboratory test itself, i.e., mainly here on which pellets are sampled?

•
This question has received little attention in the literature. In 1988, Gierman (H. Gierman, 1988) has evoked sample inhomogeneity like impregnation quality as a potential limitation to scaling-down and proposed calculations for samples made of two populations, one normal and one abnormal. He estimated the maximum acceptable activity spread so that the average rate constant would be measured within 5% of the normal value depending of the reactor volume and abnormal population fraction. Regrettably, Gierman does not detail the methodology to estimate the average reactor kinetic. As expected, he concludes that larger reactors can level out excess activity better than small reactors. To our knowledge, the effect of catalyst sampling in small units has not been studied since.

In this paper we will examine major sources of sample inhomogeneity, present a methodology to estimate its effect on kinetic evaluation, and present analytical solutions and criterions for spherical and cylindrical catalyst pellets.

2 How can pellets from the same sample be different?

Before addressing the question of a lower limit on the number of grains, one must explore the potential sources of sample inhomogeneity that may affect chemical activity. We have identified variations in size, shape, structural properties and, for supported catalysts, impregnation quality.

Shape and size control depend only on the support manufacturing process. Samples of spherical beads produced by the oil drop method are very uniform as the shape and diameter is determined by surface tension and gravity. Another process to produce spherical beads is continuous coating of successive layers of support. It results in roughly spherical particles with large diameter distributions (Figure 1). Regarding extruded pellets, also called extrudates, extrusion guarantees a very precise shape and diameter (within a batch production). Due to extrusion die wear, some evolution can be observed in the long term, potentially inducing some batch to batch differences. Extruded pellets length results from random breakages occurring on the manufacturing lines like free fall from one equipment to the next: the extrudates length is therefore not homogeneous (Figure 1). Additional dimensional variability can be observed if the pellets are subjected to attrition, or agglomeration. Dimensional variations will lead to some effects on the pellets apparent kinetics, for cases where internal mass transfer limitations occur. (single column figure)

Little has been published about structural properties dependence on process fabrication: catalyst manufacturer would certainly not publish data hinting toward lack of homogeneity.

Similarly, little has been published about impregnation homogeneity except of course for eggshell and egg-yolk type catalyst. Both structural and impregnation quality variations result in variations in the intrinsic grain kinetic activity and may be accounted for in changes of the kinetic activity.

Other sources of variability in kinetic activity can result in variable exposure, in quantity and duration, to chemicals (oxygen or moisture), temperature ... The potential change in kinetic activity will then depend on the position of the pellet in the storage container: pellets near the lid or at the bottom of the container may not be affected the same way. Similarly, during container filling, segregation of pellets according to their size is quite likely. For this reason, mixing the source container before sampling is highly recommended.

In the rest of the paper, we assume that storage procedures were adequate and that before sampling the pellets, the source batch was well mixed so that the pellets present random dimensional and intrinsic kinetic variability around "average" values.

Evaluation of sampling induced incertitude for spherical pellets

In this part, we present a methodology to estimate the effect of sampling a low number of spherical pellets, less than 1000 to set a number (representing a few grams), on the catalyst activity evaluation. More precisely, we are interested in assessing the incertitude on the catalyst activity evaluation induced by some randomness on diameter and activity.

Reactor based model

The first part of the methodology deals with establishing a mathematical relationship between the apparent catalyst activity at the reactor scale (k reactor ) and the random characteristics of the M particles tested in the reactor.

The particle properties are randomly drawn assuming: 1) known probability distributions on diameter and activity; 2) sampling reservoir is infinite. Using Thiele modulus approach, for each particle an apparent kinetic constant is computed based on diameter and intrinsic activity. We further assume a first order kinetic law. For a sphere, the equations are:

apparent k k η = Eq. 1 ( )         - = φ φ φ η 3 1 3 tanh 1 1 Eq. 2 6 p eff d k D φ = Eq. 3
The reactor apparent kinetic constant, k reactor , is computed as follows. Let's split the fixed bed in N successive thin layers of thickness Δh. For these thin layers, we approximate the activity in each layer as the average activity of all particles in that layer1 .

∑

= = i M j j particle apparent i i layer k M k 1 _ , _ 1
Eq. 4

Assuming ideal reactor hydrodynamics, the concentration evolution through each layer can be computed as follow. Combining all the layers, we have then:

_ _ _ 1 d reactor outlet reactor inlet c N layer i i c c h k c u = ∆ = - ∑ ∫ Eq. 6
Assuming the apparent kinetics at the reactor scale is also of 1 st order, then Comparing equations 6 and 7, we deduce that, for all reaction orders, the kinetic constant measured at the reactor scale can be estimated as the average of all particles apparent activity:

∑ ∑ ∑ ∑ = = = = ≈ = = M j j particle apparent N i M j j particle apparent i N i i layer reactor k M k M N k N k i 1 _ , 1 1 _ , 1 _ 1 1 1 1
Eq. 8

In Equation 8, we further assumed that the particles are perfectly homogenously distributed in the fixed bed and that all layers possess the same number of particles.

Before getting on with the use of Eq. 8 to assess incertitude due to sampling, it is interesting to comment on the chosen assumptions. The model assumes ideal hydrodynamics (plug flow) and no external mass transfer limitations. Usually, reactors for kinetic parameters evaluation are chosen and designed to be operated close to those ideal conditions so that this assumption is not restrictive. Underlying basic assumptions of the work are 1) dividing the bed in thin layers with uniform concentration and number of particles (Eq. 5) and 2) evaluating the apparent activity at thin layer scale by averaging the activity of the catalyst in the layer (Eq. 4). These assumptions are equally at the basis of all mono-dimensional models of fixed bed reactors who have proved over decades to be quite good at representing fixed bed reactors. In chemical engineering, it is quite common to model fixed beds as a cascade of stirred tank reactors. In appendix, we show that this would lead to identical results.

The aforementioned derivation is based on the integration over each layer (Eq. 5) or reactor (Eq. 7) by separating variables. This is mathematically correct if and only if the right term of the equations, more particularly k layer , is not a function of concentration. This condition is met when the efficiency is not a function of the concentration, which is always the case for a 1 st order kinetic law (Eq. 1 to 3) and at low Thiele modulus for other reaction orders, in which case the efficiency is equal to 1. It is probably possible to find other cases where variable separation is possible but this further analysis is not in the scope of this paper.

In the rest of the paper, we will assume isothermal conditions but the methodology holds as well for non-isothermal reactors: if a temperature profile is known (by measuring or by simulation), it is possible to correct the activity of each particle using for example the Arrhenius law (in Eq. 1 and 3).

The approach is also valid for any solid catalyst stacked in a fixed bed reactor. A special mention has to be made about egg shell catalyst where the active phase is impregnated on only a limited depth, generally thin but normally thicker than the reactant penetration depth.

In those catalysts, concentration profiles are identical to those that would be calculated if the catalysts were fully impregnated so that proposed methodology is still valid without modifications.

Estimation of incertitude

In the previous part, we derived a set of equation to compute an apparent kinetic constant at reactor scale (k reactor ) based on the random characteristics of a sample made of M particles. Using this set of equations, it is possible to derive information about the statistics k reactor .

The incertitude on k reactor is denoted as I and is defined so that 95% of the measured values of k reactor are within I of the real kinetic constant. Assuming that the kinetic constant follows a Gaussian distribution, the incertitude I is equal to twice the standard deviation of the k reactor distribution. The determination of I requires therefore calculating the standard deviation of a rather complex function of two random variables: diameter, d, and intrinsic kinetic constant, k.

The incertitude can be computed using a Monte Carlo numerical approach and in some particular cases using analytical equations. In the Monte Carlo approach, numerous numerical experiments are repeated and the deviations to the homogeneous case are statistically analyzed. A single numerical experiment consists of the following steps:

-Selecting randomly M particles, each being defined by a random diameter and a random activity, -Estimate the overall kinetic constant at fixed bed scale, -Compare the activity to the homogeneous case.

The numerical experiment is repeated 500 times, which has been observed to be high enough for results to be statistically representative and independent of the number of repeats. The incertitude induced by the random sampling is set to twice the standard deviations of the activity distribution. The particle properties are randomly drawn assuming: 1) known probability distributions on diameter and activity; 2) sampling reservoir is infinite.

At low and high Thiele moduli, for an isothermal reactor, equations 1-3 can be simplified and it is possible to derive analytical expressions for incertitude. In the Appendix, the mathematical tools and derivations are summarized. At a low Thiele modulus, k apparent_particle is equal to the intrinsic activity and does not depend on particle diameter. Kinetic constants estimated on the reactor are equal to the average of the particle intrinsic kinetics. We then deduce that2 :

φ <<1,

M I I k =
Eq. 9

In this case, the incertitude does not depend on the number of pellets and scales as the inverse of the square root of the number of particles.

At a high Thiele modulus, a good approximation of Eq. 2 is:

φ >>1, φ η 1 =
Eq. 10

So that φ >>1,

j particle p j particle M j eff reactor d k D M k _ , _ 1 6 1 ∑ = =
Eq. 11

Using classical composite incertitude formulas and after a few mathematical operations (see the Appendix), the incertitude on the kinetic constant measured at reactor scale, if particle diameter and intrinsic kinetic are not correlated is given by:

φ >>1, 2 2 4 1       +         = k I d I M k I k p dp
Eq. 12 Similarly to low Thiele modules, I k scales as the inverse of the square root of the number of particles. The contribution of incertitude on diameter is larger than that of incertitude on intrinsic kinetic constant; this results from difference of exponent in Eq. 11.

Defining an acceptable incertitude on kinetic constant evaluation

The proposed methodology makes it possible to compute the incertitude on the apparent kinetic constant. In order to define a suitable particle number, a maximum acceptable incertitude ought to be defined beforehand by comparison with the other sources of incertitude in screening tests. These other sources of incertitude can result from measurement of process parameters (flow-rates, pressure, and temperature) or of feedstock properties (physical properties, chemical analysis).

An interesting way to present the incertitude, and common in the industrial practice, is to express it in the form of a temperature: what is the incertitude in the temperature of a reactor filled with identical pellets that would result in the same incertitude as the one induced by sampling a low number of non-ideal particles? The reasons for this choice are 1) that catalysts are very often compared on temperature difference basis; 2) the temperature difference can be straightforwardly compared to the incertitude on repeatability of temperature measurements in fixed bed pilot units. With this definition, an incertitude of 1°C means that in 95% of the cases, the kinetic at reactor scale is estimated within ± 1 °C of that of the reactor loaded with perfectly identical pellets.

With a first order development of the Arrhenius law, one can link absolute incertitude and incertitude expressed in temperature using:

2 ref T apparent a T I R E k I ∆ =
Eq. 13

Relative incertitude and incertitude expressed in temperature are directly proportional, with a slope related to activation energy and reference temperature. Transposition from one to the other is graphically straightforward.

It is worth emphasizing here that due to our definition of the incertitude expressed in temperature, the activation energy to be used in Eq. 13 is the apparent activation energy : it has to be estimated accounting for internal mass transfer limitations. A classical result is that for high Thiele modulus, the apparent activation energy is half of the intrinsic one.

Typically, an overall experiment incertitude should remain below 1 °C in order to ensure a correct ranking of catalysts. Thus an acceptable incertitude for sampling ought to be even lower, ideally in the 0.1-0.3°C range, depending on the quality of the control and measure of process parameters and catalyst improvement target. This acceptable range has to be expressed back in incertitude on the kinetic constant evaluation, which depends strongly on both activation energy and temperature of reference.

Figure 2 presents maximum acceptable relative values, I max , for various temperatures and activation energies (10 to 150 kJ/mol) assuming a maximum sampling induced incertitude of 0.2 °C. Large acceptable values are preferred to lower sensibility to sample inhomogeneity. A reaction presenting a low activation energy and/or operated at high temperature is more sensitive to eventual sample inhomogeneity. This conclusion may be counter-intuitive but results from our choice of expressing incertitude as a temperature. For a given experimental incertitude on the kinetic constant, a higher activation energy will results in lower incertitude expressed in °C. (single column figure)

Combining Eq. 13 with previous expressions of I (Eq. 9 and 12) for low and high Thiele moduli, and the remark on the apparent activation energy, a minimum number of pellets to be used can be estimated by:

φ <<1, 2 2 2 min max ref k apparent a T T I R M E I k ∆     =           Eq. 14 φ >>1, 2 2 2 2 min max 2 1 4 ref dp k apparent a T p T I I R M E I d k ∆           = +                    

Special cases

Constant activity, Random diameter: same mean value, different spread

For pellets whose intrinsic activity is not random (I k =0), but depends only on diameter spread, analytical solutions at low and high Thiele number (Eq. 9 and Eq. 12) reduce to:

φ <<1, I = 0 Eq. 15 φ >>1, p dp d I M k I 1 =
In these conditions, the relative incertitude (I/ k ) is proportional to the relative incertitude on diameter. The high Thiele modulus solution is an upper-bound estimation.

Random activity, fixed diameter

When diameter is well controlled, and intrinsic activity presents some randomness, Eq. 9 and 12 becomes:

φ <<1, M I I k = Eq. 16 φ >>1, M I I k 2 =
In this case, there is more incertitude at a low Thiele modulus: at high a Thiele modulus, internal mass transfer limitations dampen the effect of activity variations. The number 1/2 at the denominator comes from the exponent of k (1/2) in the Thiele modulus expression. For the same reason, the apparent activation energy is halved at high Thiele modulus. This leads to an interesting result: in absence of spread on the diameter, the expression for the minimum number of pellets are identical for low and high Thiele modulus (See Eq. 14 with I dp = 0). The same can be said about the incertitude expressed in °C. As expected, Monte Carlo simulations indicate that this observation still holds for intermediate Thiele moduli.

Results for spheres

Random diameter

In Figure 3, we present the incertitude induced by the distribution of diameter for 3 levels of activity corresponding to Thiele modulus ranging from 0.1 (negligible internal limitations) to 10 (strong internal limitations). In these simulations, the intrinsic activity is identical for all particles, and diameter distribution is a Gaussian distribution centered on 3 mm with a standard deviation of 0.35 mm. As expected, severe internal mass transfer limitation (high Thiele modulus) are required for diameter distribution to have an impact on the results. Figure 3 additionally shows that: 1) Monte-Carlo simulations based on 500 repeats are in very close agreement with the analytical solutions; 2) High Thiele modulus analytical solution can greatly overestimate incertitude at moderate Thiele Modulus (~1). In those conditions, Monte-Carlo simulations can help in assessing the risks. (single column figure)

Random activity

In Figure 4, we present the incertitude induced by the distribution of activity for Thiele modulus ranging from 0.1 (negligible internal limitations) to 10 (strong internal limitations). In these simulations, the particle diameter is constant at 3 mm and the activity follows a Gaussian distribution of standard deviation 10% of average activity. This confirms that the validity of the analytical solutions (Eq. 16) and again that High Thiele modulus analytical solution can greatly overestimate incertitude at moderate Thiele Modulus (~1).

Random activity, random diameter, independent variability

When introducing independent variability of diameter (standard deviation of 0.35 mm) and on the intrinsic kinetic constant (standard deviation of 10%), as shown in Figure 5, incertitudes do not add up. This case would correspond to a quite inhomogeneous batch on the active phase. At low and moderate Thiele modulus, the incertitude is equal to the incertitude on the intrinsic constant of the pellet. The effect of diameter only becomes visible with severe internal limitation. (single column figure)

Random activity, random diameter, coupled variability

During sample preparation, one may assume that larger particles can retain more (or less) active phase than smaller particles. Hence our idea to test the consequences of a coupling between diameter and intrinsic kinetic constant.

Figure 6 presents results obtained for 500 particles drawn from a Gaussian distribution (average 3mm, standard deviation 0.35 mm). The activity is a linear of the diameter: the slope is such that a change in diameter of 0.35 mm corresponds to a change of 10% of the catalytic activity (same as uncoupled case). When the coupling is negative (a larger particle retains less metals and its intrinsic activity is lower), a very large incertitude emerges for high Thiele modulus reactions: smaller particle are more active due to their lower internal limitations and higher kinetic constant. When the coupling is positive, i.e. a larger particle retains more metals and its intrinsic activity is higher, a quite interesting result emerges: on a limited Thiele modulus range, the incertitude is nearly 0. In that range, it is possible to sample a very low number of particles without increasing incertitude. This peculiar effect is maximal around a Thiele modulus of 1, at a location that depends on the random parameters chosen for diameter distribution and coupling. At a low Thiele modulus, the main source of incertitude is rare and very active large particles, and equally rare but very inactive small particles. At high Thiele modulus, the "large particle" source of incertitude is dampened by internal mass transfer limitations but not the "small particle" source, hence a lower but non-zero incertitude.

In the intermediate range of Thiele modulus, we propose that all particles present near equal activity as both diameter effects compensate (a smaller diameter results in both a smaller intrinsic activity and lower limitations): sensibility to diameter is very weak and incertitude quite small.

Extension of methodology by modified Thiele modulus

The methodology proposed for spheres can be extended to any particle shape provided it exists some way, either analytical or numerical, to compute the efficiency knowing the dimensions (random) of the particle. There exists some analytical expression for the Thiele modulus and efficiency for infinite cylinders:

eff cylinder D k d . 4 = φ
Eq. 17 where I 1 and I 0 are Bessel function of the 1 st and 0 th order respectively.

Aris (Aris, 1957) proposed a generalized Thiele modulus valid for any particle shape:

eff D k S V . = φ Eq. 19 With L d L d S V . 4 . 2 . + =
We have checked using finite element simulations that for cylinders Eq. 11 combined with infinite cylinder efficiency (Eq. 18) predicts the actual efficiency with an error of less than 6% and opted to use these analytical functions rather than tabulated numerical results.

It is very likely that the results still holds for other extruded pellets like lobed shapes.

Analytical solutions for cylinders

Exactly like the sphere case, the expression for efficiency (Eq. 18) can be simplified for low and high Thiele modulus.

φ <<1, 1 = η Eq. 20 φ >>1, φ η 1 =
As discussed in the introduction, we will assume now that there is no variability on the pellets diameter and some variability on activity and length.

Using the mathematical tools described in the Appendix, we can derive expressions for incertitude:

φ <<1, M I I k = Eq. 21 φ >>1, 2 2 2 2 4 1 4       +               = k I L I L S V M k I k L
The sampling incertitude can be expressed in a form quite similar to that of spherical beads with a shape factor correction. Combining with Eq. 13, we can derive criterion for cylinders:

φ <<1, 2 2 max 2 min .               = ∆ k I I T E R M k T ref apparent a Eq. 22 φ >>1,                     +                       = ∆ 2 2 2 2 2 max 2 min 4 1 4 . 2 k I L I L S V I T E R M k L T ref apparent a
As commented earlier, it is very likely that the results still holds for other extruded pellets like lobed shapes.

In this paper we present a novel methodology to estimate incertitude induced by sampling a limited number of particles in fixed bed kinetic studies. The methodology limits are given by its basic assumptions: first order kinetic law, ideal hydrodynamics, dimension and intrinsic activity as main source of variability. The predicted incertitude scales as the square roof of the size of the sample: downsizing increases the incertitude especially in presence of mass transfer limitations.

For spherical and cylindrical particles, analytical formulas have been derived for low and high Thiele modulus to 1) predict incertitude induced by sampling a limited number of pellets and 2) propose a criterion on the minimum number of pellets that ought to be tested to keep sampling incertitude below an acceptable level. Results may differ greatly depending on the applications. If dimensional variability is larger than intrinsic activity variability, the minimum number of catalyst pellets to test can be lowered by narrowing the catalyst size distribution (sieving or sorting for example).

If the criterion on minimum number of pellets is not satisfied, we recommend first to try to increase the amount of catalyst loaded in the reactor paying attention to thermal issues. A next step would be to experimentally quantify the effect of sample inhomogeneity by comparing test response with reconstructed catalyst samples, for example mostly small vs. mostly large particles. If sample induced incertitude is not acceptable, then next step can be either to change the reactor or to include a sieving/sorting operation in the catalyst preparation workflow. The second option not only corresponds to capital expenditure but also higher complexity and operating costs.

Is it possible to go further in downscaling below 100 grains, let's say 50 grains? It is our belief, that on the sampling point of view at least, the answer is yes. With such limited sample population, it is easy to perform a thorough check of all the catalyst pellets before loading the reactor: optical measurement of the size of each grain, optical check of color disparity … Any outlier would be readily detected before the test. On the opposite, with such a limited number of particles, another emerging source of incertitude is the hydrodynamics that can affect drastically mass transfer, as it has been recently hinted using CFD (M. Rolland, 2013;Götz et al., 2015). We believe controlling the hydrodynamics (or modeling it) will be the major challenge toward to further downscaling.

Note about modeling with a cascade of continuous stirred tank reactors

A similar conclusion to that proposed in §3.1 can be obtained by splitting the reactor in a cascade of continuous stirred tank reactors (CSTR). In this case, assuming a 1 st order kinetic and that the activity in each CSTR reactors is the average of the apparent activity of the particles in that reactor, the concentration evolution (Eq. 5) can be re-written as : CSTR cascade equations are similar to those of layers cascade.

Q V k c c i CSTR i CSTR inlet i CSTR outlet

Composite incertitude formula

The incertitude of f(x, y), function of two independent random variables x and y, is given by 2 2 2 2 . .

y y x x f I y f I x f I         ∂ ∂ +         ∂ ∂ = .
6.3 Special cases:

1) f(x) = a.
x with a constant, then

x x x f I a I a I x f I . . . 2 2 2 2 = =       ∂ ∂ = 2) 2 1 2 1 ) , ( x x x x f + =
, with x 1 and x 2 random variables of identical incertitude I x :

( ) ( ) The apparent kinetic constant at reactor level is given by Eq. 8 which can be written in that case:

2 . 1 . 1 . . 2 2 2 2 2 2 2 2 2 1 x x x x x f I I I I x f I x f I = + =         ∂ ∂ +         ∂ ∂ = 3)

  Incertitude scales with the square root of sample population Nomenclature c x concentration [mol/m 3 ], (x = inter, outlet) D eff effective diffusion in the catalyst pellet [m 2 /s] d p or d particle diameter [m], random variable kinetic constant evaluation (at reactor scale) [s -1 ] k I absolute incertitude on intrinsic kinetic constant [s -1 ] max I maximum allowable absolute incertitude on intrinsic kinetic constant evaluation (at reactor scale) [s -1 ] T I ∆ same as I but expressed in temperature [°C]

Figure 1 :

 1 Figure 1 : Left: Sample of spherical beads with near spherical shape and various diameter (around 2 mm) -Right: Sample of trilobic extrudates with various length (around 3 mm). Glass container internal diameter: 16 mm.

Figure 2 :

 2 Figure 2: Maximum acceptable relative incertitude as a function of activation energy for 3 temperatures (80, 200 and 350°C) assuming a sampling induced incertitude of 0.2°C.

Figure 3 :

 3 Figure 3 : Incertitude for random diameters and various Thiele Modulus (0.1: negligible internal limitations, 1: medium internal limitations and 10: large internal limitations).

Figure 4 :

 4 Figure 4 : Incertitude for random activity and various Thiele Modulus (0.1: negligible internal limitations, 1: medium internal limitations and 10: large internal limitations).

Figure 5 :

 5 Figure 5 : Incertitude induced by independent variability on diameter (Gaussian distribution, average of 3 mm, standard deviation of 0.35 mm) and intrinsic kinetic constant (Gaussian distribution, standard deviation of 10%) for Thiele modulus of 0.1, 1 and 10.

Figure 6 :

 6 Figure 6 : Relative Incertitude as a function of Thiele modulus for a fixed bed of 500 particles for positive and negative coupling between diameter and activity.

  the limit of small reactors (V/Q << 1), a first order development gives

  Sum of N variables of identical incertitude I x . By extension of case 2, of N variables of identical incertitude I x . Combining case 1 and 3, we easily deduce that averaging N times a random variable reduces the incertitude by N . the particle apparent activity. We used the letter f instead of k to make the derivation easier.k reactor is the average of a composite random variable f.

  the diameter d is not a random variable.Let us define the apparent kinetic function g:

Strictly speaking, one should account for the solid fraction, but as it cancels out we prefer not to include it in the equations for clarity reasons.

Averaging N times a random variable reduces the incertitude by N . See appendix for derivation.

Acknowledgements

The authors acknowledge C. De Bellefon for discussions on the topic.

Aris, R., 1957. On shape factors for irregular particles-I: The steady state problem. Diffusion and reaction. Chemical Engineering Science 6 (6), 262-268. 10.1016/0009-2509(57)85028-3. Götz, D., Kuhn, M., Claus, P., 2015. Numerical modelling and performance studies of the original and advanced TEMKIN reactor in laboratory scale testing of industrial egg shell catalysts for the selective hydrogenation of acetylene. Chemical Engineering Research and Design 94, 594-604. 10.1016/j.cherd.2014.10.005. H. Gierman, 1988. Design of Laboratory Hydrotreating Reactors. Applied Catalysis 43, 277-286. M. Rolland, 2013. Significant Effect of catalyst Position and Orientation in Fixed-Beds Reactors: World Congress of Chemical Engineering 9. Seoul, Korea, Aug 18-23.


Julien Pilourdault

Vincent Leroy

Sihem Amer-Yahia

Distributed Evaluation of Top-k Temporal Joins

Keywords: Temporal data, Top-k, Join, Distributed processing CCS Concepts, Information systems → Join algorithms, Query operators

We study a particular kind of join, coined Ranked Temporal Join (RTJ), featuring predicates that compare time intervals and a scoring function associated with each predicate to quantify how well it is satisfied. RTJ queries are prevalent in a variety of applications such as network traffic monitoring, task scheduling, and tweet analysis. RTJ queries are often best interpreted as top-k queries where only the best matches are returned. We show how to exploit the nature of temporal predicates and the properties of their associated scoring semantics to design TKIJ , an efficient query evaluation approach on a distributed Map-Reduce architecture. TKIJ relies on an offline statistics computation that, given a time partitioning into granules, computes the distribution of intervals' endpoints in each granule, and an online computation that generates query-dependent score bounds. Those statistics are used for workload assignment to reducers. This aims at reducing data replication, to limit I/O cost. Additionally, high-scoring results are distributed evenly to enable each reducer to prune unnecessary results. Our extensive experiments on synthetic and real datasets show that TKIJ outperforms state-of-the-art competitors and provides very good performance for n-ary RTJ queries on temporal data.

INTRODUCTION

Temporal data is pervasive. Store receipts, tweets, traffic data, and temperature measures generated by weather sensors or by wearables, are just some examples. Such data is best represented as intervals with start and end timestamps. For instance, in a data center, large amounts of traf-To appear in SIGMOD '16 fic requests in the form of [IP address,start timestamp,end timestamp,. . .], are continuously generated by sensors. A system administrator who wishes to monitor traffic emanating from different countries, would formulate a query that returns pairs of requests (x, y) where x ends before y starts and x and y originate from different countries. In tweet analysis where intervals represent the lifespan of a hashtag, the query x meets y would find pairs of discussion topics where one started when the other ended.

Chronological relations between intervals are typically expressed using the Allen interval algebra [START_REF] Allen | Maintaining knowledge about temporal intervals[END_REF], which defines a set of Boolean predicates such as before, meets, starts and overlaps. In this work, we argue that the ability to evaluate interval predicates approximately and assign scores to resulting interval pairs, appears as a natural requirement to finding interesting results. Indeed, in monitoring traffic causality, an (x, y) pair where x ends just before y, may be preferred to those where x ends much earlier than the start of y. Similarly, in tweet analysis, a Boolean interpretation of x meets y is likely to return an empty result and will not detect discussion topics that started roughly after another ended.

In this paper, we formalize n-ary Ranked Temporal Join (RTJ) queries that feature any number of temporal predicates and return the k best results. To the best of our knowledge, this is the first work that formalizes n-ary RTJ queries and devises a general-purpose distributed evaluation approach. RTJ queries raise a number of new challenges. First, an appropriate ranked semantics needs to be devised for a variety of temporal predicates in order to capture the strength of relationship between intervals as a function of the desired semantics for each predicate. Most related studies focus on a common query that involves the intersects predicate [START_REF] Dignös | Overlap interval partition join[END_REF][START_REF] Enderle | Joining interval data in relational databases[END_REF][START_REF] Lu | On spatially partitioned temporal join[END_REF]. The idea is to retrieve tuples that share a common period of validity, in order to combine events whose time range has a non-empty intersection. Our semantics is richer since we are interested in any temporal predicate between intervals. For instance, we aim to capture predicates from the Allen algebra (see the three first columns of Figure 2). We illustrate that in the following example.

Network Traffic Monitoring. Consider the two interval collections C1 and C2 in Figure 1. Assume that each collection gathers requests from a different country. An analyst interested in monitoring traffic between two countries would seek (x, y) interval pairs, x ∈ C1, y ∈ C2, that form a sequence, i.e., where y starts as x ends. The best sequences are those satisfying meets(x, y) reflecting strict equality between the end of x and the beginning of y. In our example, only (x4, y4) qualifies for a Boolean interpretation. Yet, given the uncertainty on clocks in different network equipments, a small overlap or a short free period between two consecutive tasks would be more realistic. In order to express that with Boolean semantics, a query must evaluate a disjunction of before, overlaps and meets and would build many useless interval pairs. A ranked semantics on the other hand, would build (x, y) pairs where x almost meets y, allowing a tolerance on intervals' endpoints. For instance, tuples (x, y) where x ends at most 2 timestamps before or after y starts, could be considered high-scoring. Using such a semantics, we can return the top-3 {(x4, y4), (x1, y3), (x1, y1)}.

A key observation we rely on is that an approximate interpretation of a predicate amounts to approximating the strength of relationship between its intervals' endpoints. That is compatible with the flexible scoring approach proposed in [START_REF] Dubois | Fuzziness and uncertainty in temporal reasoning[END_REF] to approximate Allen predicates [START_REF] Allen | Maintaining knowledge about temporal intervals[END_REF]. It is based on associating a degree of satisfaction with equality and inequality of 2 intervals' endpoints. In this paper, we adapt this framework to allow scoring any temporal predicate. We also allow to use any monotone aggregation function to compute the score of a query result for n-ary RTJ queries. Consequently, a 3-way query involving the predicates starts(x, y) and meets(y, z) would return (x, y, z) tuples and associate to each tuple its degree of satisfaction of the query as an aggregation of individual predicate-dependent scores for each of starts(x, y) and meets(y, z).

It is important to note that, unlike existing work on returning approximate answers [START_REF] Candan | Rankloud: Scalable multimedia data processing in server clusters[END_REF], we focus on computing exact answers of queries that make use of scored join predicates. The second challenge we tackle is to devise an efficient query evaluation strategy that guarantees to return the best answers for a variety of temporal predicates. The key difficulty here is to develop a general-purpose algorithm that works with a variety of predicates and ranked semantics and yet, that is able to exploit the nature of those predicates to devise an efficient evaluation.

The efficient processing of interval joins has been studied before [START_REF] Chawda | Processing interval joins on map-reduce[END_REF][START_REF] Dignös | Overlap interval partition join[END_REF][START_REF] Enderle | Joining interval data in relational databases[END_REF][START_REF] Lu | On spatially partitioned temporal join[END_REF]. The closest to our work is the recent one by Chawda et al. [START_REF] Chawda | Processing interval joins on map-reduce[END_REF] with a focus on processing queries on Map-Reduce [START_REF] Dean | Mapreduce: Simplified data processing on large clusters[END_REF]. However, their algorithms focus on a Boolean semantics and are not directly applicable in our case. In our work, scores constitute both a challenge and an opportunity. They are a challenge because, unlike Boolean semantics, every combination of intervals is potentially an answer. The opportunity lies in the ability to leverage statistics on input data in order to avoid computing low-scoring results. In this paper, we propose TKIJ , a 3-stage query evaluation approach for RTJ queries on Map-Reduce. A first offline process partitions time into granules and computes the distribution of intervals' endpoints in those granules. A second online process leverages those statistics to compute query-dependent score bounds. A third process assigns to each reducer enough high-scoring results to enable pruning unnecessary answers and hence balance the load between reducers. Each reducer executes the input RTJ query locally on the data it receives. While this query processing bears similarities to traditional top-k and rank-join processing [START_REF] Fagin | Combining fuzzy information from multiple systems[END_REF][START_REF] Fagin | Optimal aggregation algorithms for middleware[END_REF][START_REF] Finger | Robust and efficient algorithms for rank join evaluation[END_REF][START_REF] Ilyas | Supporting top-k join queries in relational databases[END_REF][START_REF] Ntarmos | Rank join queries in nosql databases[END_REF], it differs greatly due to the dependence of score computation on query predicates. As a result, partial predicate-dependent scores cannot be pre-computed and each reducer relies on the score bounds produced in the statistics collection phase to avoid unnecessary computation. A final Map-Reduce job merges results obtained from each reducer.

In summary, this paper makes the following contributions:

1. We formalize n-ary RTJ queries that combine interval collections with any number of arbitrary temporal predicates and return the k best results.

2. We design TKIJ , an efficient query evaluation approach on Map-Reduce that leverages scores to avoid computing unnecessary results and to balance the load between reducers.

3. We run extensive experiments on synthetic and real network traffic datasets. Our experiments show the efficiency of our pruning technique and the great benefit of using scores to distribute the load between reducers.

Because TKIJ executes only combinations that ensure to return the correct top-k answers, it scales to very large collections (up to 5M tuples per collection) and to high k values (up to 10 5 ).

Section 2 describes our data model and defines the problem of evaluating RTJ queries. TKIJ is presented in Section 3. Experiments are detailed in Section 4. We summarize the related work in Section 5 and conclude in Section 6.

DATA MODEL AND PROBLEM

We are given m collections of intervals C1, . . . , Cm. Each interval x has a unique identifier, a start time x and an end time x.

Boolean temporal predicates. The general form of a temporal predicate between two intervals x and y is denoted p(x, y) and is expressed as a Boolean conjunction of equalities and inequalities between their endpoints x, x, y, y. This allows to capture a wide range of predicates among which the seminal Allen algebra [START_REF] Allen | Maintaining knowledge about temporal intervals[END_REF]. The first 3 columns of Figure 2 summarize Allen temporal predicates and their semantics. For example, meets(x, y) imposes that y starts when x finishes while starts(x, y) requires that x and y start at the same time and that x ends before y. It is important to note that we aim to capture all Allen predicates but also any predicate comparing interval's endpoints. While we do not aim to provide a long list of new predicates, we discuss examples that we are using in our experiments. For example, in network traffic analysis, we introduce justBefore(x, y) that is satisfied iff y > x and y -x ≤ AVGz(z -z). The intuition is that the elapsed time between x and y is no greater than the average interval length. A special case is the predicate shiftMeets(x, y) that is satisfied iff (y -x) = AVGz(z -z). In tweet analysis, a possibly useful predicate would be sparks(x, y) which is satisfied iff (y -y) > 10 * (x -x) and y > x. As a result, the (x, y) pairs satisfying sparks(x, y) would identify hashtag pairs where the preceding hashtag lasted 10 times shorter that the following. That could be useful in determining causality of long-lasting events such as finding all shortlasting hashtags before the long-lasting #JeSuisCharlie.

Scored temporal predicates. Since we are interested in capturing the degree at which a temporal predicate is verified by a pair of intervals, we propose to associate a score to each predicate. Here again, we aim to be general and we adopt the flexible approach for scoring Allen predicates [START_REF] Dubois | Fuzziness and uncertainty in temporal reasoning[END_REF] and adapt it to our settings. This approach relies on two primitive approximation comparators on intervals' endpoints. Those comparators, equals(x, y) and greater(x, y), are used to express the degree of equality or inequality of intervals' endpoints x and y, where x ∈ {x, x}, y ∈ {y, y} as a graded value in [0, 1]. They rely on two parameters λ and ρ that provide flexibility in controlling the tolerance degree when comparing intervals' endpoints. Figure 3 shows how equals(x, y) and greater(x, y) are used with λ and ρ to express that tolerance. For instance, by defining that whenever |x -y| ≤ λ, equals(x, y) returns 1, λ sets a tolerance for exact endpoint equality. ρ, on the other hand, is used to define score values. A large ρ value defines a wide range of score values and a small ρ produces a more abrupt curve and fewer possible score values. Since temporal predicates are expressed as equalities and inequalities on intervals' endpoints, their approximation can be achieved using a conjunction of equals() and greater() with appropriate λ and ρ values. This allows us to associate a scored variant to each temporal predicate. We denote that variant s-p(x, y) and refer to it as scored temporal predicate, abusing the term "predicate" to mean "function". Indeed, while p(x, y) returns a Boolean value, s-p(x, y), returns a score in [0,1]. For example, we can define the scored version of starts(x, y) as s-starts(x, y) = min{equals(x, y), greater(y, x)}. We also propose to allow different values of λ and ρ for equals() and greater() for different predicates. That provides a finer control of the score values produced by each predicate. A Boolean interpretation of a predicate becomes a special case of our scored interpretation. For example, strict endpoint equality can be obtained by setting both λ equals and ρ equals to 0. Thus, we can define, s-justBefore(x, y) with λgreater and ρgreater set to 0, ρ equals to any value and λ equals to AVGz(z -z) (Figure 4).

Temporal join queries. We are interested in expressing n-ary join queries on interval collections C1, . . . , Cm. We express a query Q as a weakly connected oriented simple graph1 of the form (V, E). Each each vertex vi ∈ V is mapped to a collection Ci. Each edge (i, j) ∈ E between two vertices vi and vj is labeled with a scored temporal predicate s-p (i,j) () between the two collections Ci and Cj corresponding to vi and vj.

The evaluation of an n-ary join query Q returns a set of tuples of the form (x1, . . . , xn) where xi ∈ Ci. The score of each tuple in the query result is computed using a function S that aggregates the partial scores assigned by each predicate s-p (i,j) () associated with each query edge (i, j) ∈ E. S could be any monotone function such as the weighted sum as it is commonly the case in ranked aggregation [START_REF] Fagin | Optimal aggregation algorithms for middleware[END_REF][START_REF] Finger | Robust and efficient algorithms for rank join evaluation[END_REF][START_REF] Ilyas | Supporting top-k join queries in relational databases[END_REF][START_REF] Ntarmos | Rank join queries in nosql databases[END_REF].

For example, we can express a 3-way query that returns a tuple (x, y, z) where x ∈ C1, y ∈ C2 and z ∈ C3 and the score of (x, y, z) is computed as an aggregation of its partial scores for query predicates s-starts(x, y) and s-meets(y, z).

Although we use the term "join" to refer to our queries, their expressivity goes beyond traditional relational joins. Ranked Temporal Join (RTJ) problem. Given an n-ary temporal join query Q = (V, E) expressed over a set of collections C1, . . . , Cm corresponding to query vertices in V and temporal predicates s-p (i,j) () associated to each edge (i, j) ∈ E, our problem is to find a top-k set of tuples of the form (x1, . . . , xn), xi ∈ Ci, ranked by (descending) order of S (i,j)∈E (s-p (i,j) (xi, xj)).

TEMPORAL JOIN PROCESSING

We present TKIJ, our approach for evaluating Top-K Interval Joins, that efficiently finds the set of k best results for an RTJ query Q. We first provide an overview of TKIJ , then we give each step in detail.

Overview of TKIJ

Figure 5 summarizes TKIJ . Given a set of interval collections C1 . . . Cm, TKIJ executes a query-independent preprocessing phase to collect statistics on intervals' distribution. This phase partitions time into granules and computes buckets for each collection (a). A bucket associated to a collection Ci corresponds to a pair of granules, and contains the number of intervals of Ci starting at one granule and ending at another. Given a query Q, these statistics are used to evaluate bucket combinations that should be processed in order to obtain top-k results (b). TKIJ relies on a constraint programming solver to compute score bounds for each bucket combination and uses those bounds to prune combinations that do not contain top-k results. The third phase is the actual join processing which relies on two Map-Reduce jobs. The first job assigns a subset of buckets to each reducer rj (c) which then processes locally the RTJ query, returning local top-k results (d). This assignment aims at reducing data replication to limit I/O, and leverages score bounds to distribute high-scoring results evenly so that each reducer can quickly prune low-ranking results. The second Map-Reduce job merges all local results into a single query output (e).

Statistics collection

TKIJ pre-processes each dataset once in order to collect statistics which are then used to optimize the execution of any RTJ query on this dataset. These statistics maintain a matrix Bi representing the distribution of endpoints of intervals in each collection Ci. TKIJ partitions the time range of each Ci into a set of contiguous granules. We adopt a uniform partitioning which has been shown to be appropriate for temporal joins [START_REF] Chawda | Processing interval joins on map-reduce[END_REF][START_REF] Dignös | Overlap interval partition join[END_REF][START_REF] Gao | Join operations in temporal databases[END_REF].

As illustrated in Figure 5a, each matrix entry records the cardinality of a bucket, where a bucket b i,l,l ′ = (g i,l , g i,l ′ ) contains all intervals of Ci that start in g i,l and end in g i,l ′ :

Bi[l][l ′ ] = |b i,l,l ′ | = |{x ∈ Ci, x ∈ g i,l ∧ x ∈ g i,l ′ }|.
As an example, given g1,1 = [10, 20] and g1,2 = [START_REF] Ilyas | Supporting top-k join queries in relational databases[END_REF][START_REF] Zhang | Efficient multi-way theta-join processing using mapreduce[END_REF], the matrix entry for b1,1,2 = (g1,1, g1,2) indicates 6 intervals starting in [START_REF] Doulkeridis | Processing of rank joins in highly distributed systems[END_REF][START_REF] Ilyas | Supporting top-k join queries in relational databases[END_REF] and ending in [START_REF] Ilyas | Supporting top-k join queries in relational databases[END_REF][START_REF] Zhang | Efficient multi-way theta-join processing using mapreduce[END_REF].

Range partitioning is a common approach in temporal join processing [START_REF] Chawda | Processing interval joins on map-reduce[END_REF][START_REF] Dignös | Overlap interval partition join[END_REF][START_REF] Gao | Join operations in temporal databases[END_REF][START_REF] Lu | On spatially partitioned temporal join[END_REF]. The rationale is that intervals having similar endpoints are likely to satisfy similar join predicates. For example, most previous studies, that focus on intersection joins, leverage partitions to avoid pairs of intervals that are guaranteed not to intersect. Similarly, TKIJ relies on these statistics to obtain information on the distribution of intervals within buckets and prune the search space of any RTJ query.

Statistics are computed in a single Map-Reduce phase. Each mapper reads a fraction of the data and maintains a local matrix per collection. Matrices are then aggregated in the reduce phase, and the reducer responsible for collection Ci outputs a final matrix Bi. While we focus in this paper on the case of statistics computed from scratch for a new dataset, we can easily handle updates by applying the same process on the inserted/deleted data.

Selection of bucket combinations

We now describe how TKIJ uses pre-computed statistics to estimate score bounds on candidate results. Then, we present how score bounds are used to avoid computing unnecessary results while we guarantee to return the exact topk results. We further develop several processing strategies that aim to tackle computational costs raised by this pruning step.

Estimating score bounds

Processing an RTJ query Q requires to return the top-k tuples (x1, . . . , xn), xi ∈ Ci according to a scoring function S. Since any tuple (x1, . . . , xn) is a potential answer, we investigate how to reduce the amount of data processed using scores. We use ω

= (b 1,l 1 ,l ′ 1 , . . . , b n,ln,l ′ n ) to denote a bucket combination, ω.nbRes = ∏ n i=1 |b i,l i ,l ′ i
| the total number of results that can be obtained from a bucket combination ω, and Ω the set of all combinations. We define score upper and lower bounds in each ω, denoted ω.UB and ω.LB .

Definition 1. The score upper-bound (resp. lower-bound)

ω.UB (resp. ω.LB) of a bucket combination ω = (b 1,l 1 ,l ′ 1 , . . . , b n,ln,l ′ n ) is the upper-bound (resp. lower-bound) of S (i,j)∈E (s-p (i,j) (xi, xj)) where xi ∈ g i,l i , xi ∈ g i,l ′ i , ∀i ∈ 1 . . . n.
As an example, suppose that query Q features a predicate s-meets (1,2) (x, y) where x ∈ C1 and y ∈ C2, using scoring parameters (λ equals , ρ equals ) = (4, 8). Collected statistics show 6 intervals in bucket b1,1,2 = ( [START_REF] Doulkeridis | Processing of rank joins in highly distributed systems[END_REF][START_REF] Ilyas | Supporting top-k join queries in relational databases[END_REF], [START_REF] Ilyas | Supporting top-k join queries in relational databases[END_REF][START_REF] Zhang | Efficient multi-way theta-join processing using mapreduce[END_REF]) for C1 and 7 intervals in bucket b2,2,3 = ( [START_REF] Ilyas | Supporting top-k join queries in relational databases[END_REF][START_REF] Zhang | Efficient multi-way theta-join processing using mapreduce[END_REF], [START_REF] Zhang | Efficient multi-way theta-join processing using mapreduce[END_REF]40]) for C2. We build the bucket combination ω = (b1,1,2, b2,2,3). Then, we can derive bounds on the score S(x, y) = s-meets(x, y) of a result (x, y) ∈ ω. The maximum possible score is 1 (e.g. with (x, y) = ( [START_REF] Dubois | Fuzziness and uncertainty in temporal reasoning[END_REF][START_REF] Prud'homme | Choco3 Documentation[END_REF], [START_REF] Prud'homme | Choco3 Documentation[END_REF]35])), and the minimum score is 0.25 (with (x, y) = ( [START_REF] Fagin | Optimal aggregation algorithms for middleware[END_REF][START_REF] Ilyas | Supporting top-k join queries in relational databases[END_REF], [START_REF] Zhang | Efficient multi-way theta-join processing using mapreduce[END_REF]35])). Hence, ω.UB = 1, ω.LB = 0.25. Thus, 42 results in ω have a score in [0. [START_REF] Prud'homme | Choco3 Documentation[END_REF][START_REF] Afrati | Optimizing joins in a map-reduce environment[END_REF].

TKIJ relies on a constraint programming solver as a generic approach to compute score bounds for any combination of predicates. Computing score bounds for a bucket combination requires to solve the following problem:

Bounds Problem. Let ω = (b 1,l 1 ,l ′ 1 , . . . , b n,ln,l ′ n ). Find (x1, . . . ,
s.t. x i ∈ g i,l i ∀i ∈ 1 . . . n (1) xi ∈ g i,l ′ i ∀i ∈ 1 . . . n (2)
Each xi ∈ Ci is mapped to a decision variable xi. For each partial score s-p (i,j) (xi, xj), we create an intermediate variable s-pij. For all (i, j) ∈ E, we impose that the variables xi, xj and s-pij satisfy the constraints {Cij : s-pij = s-p (i,j) (xi, xj)}. Then, we create a variable score and impose Cs : score = S (i,j)∈E (s-pij). The solver then maximizes (and minimizes in the case of a lower-bound) score such that constraint Cs, all constraints Ci and all constraints on decision variables (1)(2) are satisfied. While we virtually allow any temporal predicates, in practice, predicate implementation depends on the range of constraints supported by the server used in the implementation.

Pruning bucket combinations

TKIJ leverages computed score bounds to reduce computation cost by eliminating results that are guaranteed not to be in the top-k. To do so, it computes Ω k,S ⊆ Ω, a subset of the search space that is sufficient to guarantee correctness. We define Ω k,S as follows: Definition 2. The set of Top Buckets Ω k,S is a subset of Ω satisfying the following conditions: This definition ensures that whenever a bucket combination ω is pruned, there are at least k results from Ω k,S with a score higher than results generated from ω. Note that Ω k,S is not unique: given a valid set of bucket combinations, any super-set is also valid. To compute Ω k,S , we design the getTopBuckets algorithm (Algorithm 1). Algorithm getTop-Buckets uses as input a set of bucket combinations whose score bounds are pre-computed. It first computes a lower bound kthResLB on the score of the k th result (Lines 1-6). Then, it keeps only bucket combinations whose score upperbound is greater than kthResLB (Lines 7-13). The process safely stops in Line 11 since no bucket combination outside the collected set has results with score above kthResLB.

• ∀ω ∈ Ω \ Ω k,S ∃Ψ ⊆ Ω k,S : -∀ω ′ ∈ Ψ ω ′ .LB ≥ ω.UB - ∑ ω ′ ∈Ψ ω ′ .nbRes ≥ k
Pruning unnecessary results is a two-step process, coined TopBuckets. A first step computes score bounds for bucket combinations using a solver. Then, a second step executes getTopBuckets that uses those bounds to eliminate unnecessary results. In our setting, all (x1, . . . , xn) combinations are potential answers, and we cannot employ traditional top-k techniques to prune the search space. TopBuckets addresses this challenge using pre-computed statistics to locate high-scoring answers. Still, a new challenge is to limit the overhead of TopBuckets due to computing score bounds for bucket combinations. In the following section, we discuss this challenge and possible solutions.

TopBuckets Strategies

A first straightforward approach to find Ω k,S is to build all possible bucket combinations Ω, compute their score bounds using a solver and then use getTopBuckets to prune useless ones. In this strategy, coined brute-force, a large number of n-tuples of buckets are assigned a score by the solver (with g granules per collection, |Ω| is O(g 2n )). Moreover, for each combination, 2n decision variables need to be assigned by the solver. Thus, as n or g increase, brute-force becomes inefficient.

To tackle that, we propose the loose strategy (Algorithm 2 with the flag onePhase set true.) loose first builds all bucket pairs (b i,l i ,l ′ i , b j,l j ,l ′ j ) for each scored predicate (i, j) ∈ E (Line 1). Score bounds are then computed by the solver (Line 3) for each pair. Then, loose builds n-tuples of buckets (Lines 4-5). For each ω = (b 1,l 1 ,l ′ 1 , . . . , b n,ln,l ′ n ), we obtain score bounds using bounds computed for each pair (b i,l i ,l ′ i , b j,l j ,l ′ j ). Compute ω.UB, ω.LB using solver return TopBuckets(Lm)

To calculate correct bounds, we rely on the monotonicity of S. Without loss of generality, suppose that S is monotonically increasing. In the expression of S, we replace each partial score s-p (i,j) (xi, xj) with the upper bound of the corresponding pair of bucket. Therefore S (i,j)∈E ((b i,l i ,l ′ i , b j,l j ,l ′ j ).UB) is a correct upper-bound for ω. Then, loose runs getTopBuckets on the generated bucket combinations (Line 6) and returns the selected combinations. The rationale behind loose is that processing time can decrease significantly because (i) fewer bucket combinations need to be assigned a score bound (their number is O(|E|•g 4 )) and (ii) the solver needs to assign only 4 variables when computing score bounds. A drawback of loose is that the aggregation of bounds using S may result in loose bounds. We illustrate that in the following example.

Example. Figure 6 depicts a dataset with three buckets b1, b2, b3 from C1, C2, C3. We have b1 = (g1, g2), b2 = (g2, g3), b3 = (g3, g3), where g1 = [START_REF] Doulkeridis | Processing of rank joins in highly distributed systems[END_REF][START_REF] Ilyas | Supporting top-k join queries in relational databases[END_REF], g2 = [START_REF] Ilyas | Supporting top-k join queries in relational databases[END_REF][START_REF] Zhang | Efficient multi-way theta-join processing using mapreduce[END_REF], g3 = [START_REF] Zhang | Efficient multi-way theta-join processing using mapreduce[END_REF]40]. Our query features scored predicates s-starts (1,2) (x, y) and s-starts (2,3) (y, z) and our aggregation function is the normalized sum. We use the scoring parameters {(λ equals , ρ equals ), (λgreater, ρgreater) = {(1, 3), (0, 4)}. loose first computes bounds for ω1 = (b1, b2). We have ω1.UB = 1 (because s-starts (1,2) (x1, y1) = 1), and ω1.LB = 0 (s-starts (1,2) (x1, y2) = 0). Then, loose computes bounds for ω1 = (b2, b3). We have ω2.UB = 1 (s-starts (2,3) (y2, z1) = 1), and ω2.LB = 0 (s-starts (2,3) (y2, z2) = 0). Then, loose merges combinations ω1, ω2 in ω3 = (b1, b2, b3), and computes ω3.UB = S(1, 1) = 1, ω3.LB = S(0, 0) = 0. Yet, brute-force, computes ω3.UB = 0.5 because there is no (x, y, z) such that s-starts (1,2) (x, y) = s-starts (2,3) (y, z) = 1 given the buckets' bounds: it is impossible to have both equals(x, y) = 1 and equals(y, z) = 1 with x ∈ g1, y ∈ g2 and z ∈ g3. Thus, in this example, loose returns an exact lower-bound and a loose upper-bound, while brute-force returns tight bounds. These observations lead to propose a third strategy twophase, that combines brute-force and loose. two-phase is executed by Algorithm 2 when the flag onePhase is set to false. First, two-phase computes loose bounds to eliminate some bucket combinations (Lines 1-7), identically to loose. Then, two-phase refines the bounds of the remaining combinations (Lines 8-9) to obtain exact bounds. The rationale behind two-phase is that its first phase may help reduce the number of bucket combinations that need to be assigned a score in the second phase, thus improving the solver's running time. Unlike loose, two-phase returns tight bounds When selecting bucket combinations, TKIJ runs TopBuckets using one of the three strategies presented. Each strategy (i) employs the solver and (ii) executes getTopBuckets once or twice using loose or tight score bounds on bucket combinations.

Distributed Top-k Join Processing

The TopBuckets process generates Ω k,S , a set of bucket combinations that are sufficient to accurately compute the top-k results. We now describe how TKIJ computes the top-k results (Steps (c)-(d)-(e) in Figure 5). We implement TKIJ on Map-Reduce [START_REF] Dean | Mapreduce: Simplified data processing on large clusters[END_REF]. Given a set of r reducers, TKIJ assigns each bucket combination ω ∈ Ω k,S to a single reducer rj, j ∈ 1 . . . r, that processes results in ω. The main challenge in distributed join processing is to devise an efficient workload assignment function. When performing large-scale joins, I/O often constitutes a major bottleneck. We first review existing assignment algorithms, then we consider the specifics of distributed top-k computation and show that it is essential to take scores into account when dividing the workload. We present DistributeTopBuckets, a novel function that focuses on assigning high-scoring results to each reducer, while minimizing I/O cost as a secondary objective. Finally, we present how an RTJ query is processed using appropriate Map-Reduce algorithms.

Existing I/O optimizations

When different reducers require the same chunk of data, this data is replicated in the shuffle phase of Map-Reduce, which increases input cost. Several distributed join algorithms, such as RCCIS [START_REF] Chawda | Processing interval joins on map-reduce[END_REF] and the work of Afrati et al. [START_REF] Afrati | Optimizing joins in a map-reduce environment[END_REF] specifically aim at reducing that cost. In TKIJ , this corresponds to different reducers being assigned bucket combinations involving the same bucket. Other approaches focus on assigning a balanced load to each reducer [START_REF] Chawda | Processing interval joins on map-reduce[END_REF][START_REF] Okcan | Processing theta-joins using mapreduce[END_REF]. This ensures that the number of results generated by each reducer is comparable, so that no reducer will have a larger workload in output dominated tasks. Finally, some algorithms optimize both input and output costs simultaneously [START_REF] Zhang | Efficient multi-way theta-join processing using mapreduce[END_REF]. All these approaches are not directly applicable to our settings. They achieve optimizations for specific queries (equi-join [START_REF] Afrati | Optimizing joins in a map-reduce environment[END_REF], 2-way θ-join [START_REF] Okcan | Processing theta-joins using mapreduce[END_REF], m-way θ-join [START_REF] Zhang | Efficient multi-way theta-join processing using mapreduce[END_REF]). One close related work to ours [START_REF] Chawda | Processing interval joins on map-reduce[END_REF] reduces I/O cost by leveraging the Boolean interpretation of Allen predicates. That is not directly applicable to scored predicates.

Top-k optimizations

TKIJ significantly differs from standard Map-Reduce-based join processes due to its ranked semantics. In TKIJ , each reducer processes a full RTJ query locally using the bucket combinations it receives (Figure 5d). Hence, it is important for all bucket b ∈ ω ▷ Assign buckets in ω to rj 6:

M ← M ∪ (b, rj) return M to ensure that each reducer quickly identifies high-scoring results as it is usually the case in top-k processing [START_REF] Fagin | Combining fuzzy information from multiple systems[END_REF][START_REF] Fagin | Optimal aggregation algorithms for middleware[END_REF][START_REF] Finger | Robust and efficient algorithms for rank join evaluation[END_REF][START_REF] Schnaitter | Evaluating rank joins with optimal cost[END_REF]. Therefore, the assignment of bucket combinations to reducers favors an even distribution of high-scoring results.

DTB algorithm

TKIJ relies on the DistributeTopBuckets algorithm (Algorithm 3) to assign bucket combinations from Ω k,S to reducers. Following the principles described above, DTB increases the probability that each reducer receives a fair share of high-scoring results. This step relies on the knowledge, for each bucket combination, of the number of results generated, as well as their score bounds. DTB first sorts Ω k,S by descending order of score upper-bound (Line 1) to access them according to their likelihood of generating highscoring results. It then assigns each bucket combination to a reducer using the getReducer function (Algorithm 4), which returns a reducer among the ones that were assigned the fewest bucket combinations so far. Furthermore, DTB opportunistically optimizes I/O cost. First, in the worst case, a reducer evaluates all the results it is assigned. If Ωr j is the set of bucket combinations assigned to a reducer rj, then ∑ ω∈|Ωr j | ω.nbRes would be computed. DTB first computes the average number of results assigned to reducers (Algorithm 3, Line 2). Then getReducer ensures that reducers that are already assigned more than twice the average number of results are discarded (Algorithm 4, Lines 3, 7). This heuristic limits imbalance in a worst-case scenario. In the case, where several reducers have received the same of bucket combinations, getReducer selects the reducer that was already assigned the largest fraction of current ω from previous overlapping bucket combinations (Algorithm 4, Lines 8-10). For a given w, getReducer evaluates for each reducer rj the input cost of assigning w to rj using inCost(rj, ω). We define inCost(rj, ω) =

∑ b i,l,l ′ ∈ω |b i,l,l ′ | • Φ(rj, b i,l,l ′ ) where Φ(rj, b i,l,l ′ ) returns 1 if b i,l,l ′
was already assigned to rj, else 0. This process favors assignments that reduce replication cost. Having selected the reducer rj that has to be assigned the current bucket combination ω, DTB stores all the assignments (b i,l,l ′ , rj) where b i,l,l ′ ∈ ω (Algorithm 3, Lines 5-6). These assignments determine to which reducers buckets are communicated, ensuring both output correctness and join processing efficiency.

Join Processing

The final phase of TKIJ first runs DTB using Ω k,S to determine data distribution among reducers. Then, a Map-Reduce phase processes the RTJ query locally. For each Algorithm 4 getReducer Input: ω bucket combination to assign, avgRes average number of results per reducer Output: reducer to which ω has to be assigned 1: min ω assigned = +∞ 2: for j = 1 to r ▷ Retrieve the minimum amount of bucket combinations assigned 3:

if rj.nbRes < 2 × avgRes 4:

min ω assigned = min{|Ωr j |, min ω assigned} 5: minCost = +∞ 6: for j = 1 to r ▷ Find the best reducer wrt inCost(ω, rj) 7: if rj.nbRes < 2 × avgRes ∧ |Ωr j | = min ω assigned 8:

if inCost(rj, ω) < minCost 9: bestReducer=rj 10: minCost = inCost(rj, ω) return bestReducer input interval x, a mapper computes the bucket b ix,lx,l ′

x in which x falls. Then, x is communicated to all reducers rj that received b ix,lx,l ′

x . That way, each reducer rj receives its share of input data, and a list of bucket combinations Ωr j ⊆ Ω k,S whose results are potential top-k candidates. Once each reducer has processed locally the RTJ query, we run an additional Map-Reduce phase (Step (e) in Figure 5), that merges local results and returns the final top-k answers.

EXPERIMENTS

Platform. We conduct experiments on an 8-node industrial cluster with 6 workers. Each worker has 1 Intel Xeon E5-2650L (8 cores), 32GB RAM, 5TB disk. Machines run Centos 6.6 with Cloudera 5.2.5 and Hadoop 2.5.0. All presented results are averages of 5 consecutive runs.

Statistics collection. We use the same number of granules g for each collection. We observed that only the number of interval |Ci| per collection had a significant impact on statistics collection time. Statistics collection lasted between 28s for |Ci| = 2×10 5 and 36s for |Ci| = 5×10 6 . Since this task is only executed once as a pre-processing for a given dataset, we do not include it in query evaluation time.

Selection of bucket combinations. We implement a distributed and multi-threaded version of TopBuckets. We split the set of buckets B1 into 6 equal-sized groups B1,j, j ∈ 1 . . . 6. Then each worker j ∈ 1 . . . 6 runs a local version of TopBuckets using buckets in B1,j and all buckets in Bi, i ∈ 2 . . . n. Thus, each worker has as input a disjoint set of possible bucket combinations. We use Choco [START_REF] Prud'homme | Choco3 Documentation[END_REF], a software for constraint programming to compute score bounds for each bucket combination. Each execution of the solver is handled in a separate thread. A second phase of TopBuckets merges local TopBuckets results on a selected worker and returns the set of best bucket combinations Ω k,S .

Distributed join processing. We run TKIJ using 24 reducers. Local query execution accesses the received bucket combinations by descending order of score upper-bounds. It then uses R-Trees to access intervals in memory. For an interval xi and a score value v, it queries the R-Tree and returns only intervals xj s.t. s-p (i,j) (xi, xj) ≥ v. Queries. Tables 1 and2 summarize queries and score parameters. We use S =

Id Scored Temporal Predicates In Q. xi ∈ Ci ∀i ∈ 1 . . . n. Q b,b s-before(x1, x2), s-before(x2, x3) Q f,f s-finishedBy(x1,
∑ (i,j)∈E s-p (i,j) (x i ,x j ) |E|
to compute the score of a query result (x1, . . . , xn). Unless otherwise specified, k = 100.

Summary of Results

We show that TKIJ processes various RTJ queries efficiently on both synthetic data and real network traffic logs. TKIJ scales to collections of up to 5 million intervals (|Ci| ∈ [1M, 5M ]) and efficiently returns the top-k results for k ∈ [10, 10 5 ]. We show that our workload distribution approach, DistributeTopBuckets, that fairly distributes highscoring results, outperforms a more naive approach based on the LPT algorithm. That is particularly useful for queries that return few high-scoring results, such as those featuring equality-based predicates (e.g starts). We observe that the efficiency of DistributeTopBuckets decreases with coarser statistics. We also show that as the number of collections or the number of predicates in a query vary, TopBuckets efficiently prunes buckets combinations. Experiments on network traffic data show that on queries featuring before or overlaps, TopBuckets can select few bucket combinations guaranteeing TKIJ to return high-scoring results. We observe that a higher number of granules g (finer statistics) helps prune unnecessary results, which improves overall join processing time. However, it also makes pruning computation with TopBuckets slower. Hence, we run experiments to find a sweet-spot value for g.

Synthetic Data

To generate synthetic data, we use the same parameters as in previous work [START_REF] Chawda | Processing interval joins on map-reduce[END_REF]. We use a pseudo-random uniform generator to get intervals' startpoints and lengths in specified ranges (respectively s = [0, 10 5 ] and w = [START_REF] Afrati | Optimizing joins in a map-reduce environment[END_REF]100]). Intervals' endpoints are integers. We vary the number |Ci| of intervals per collection, and the number n of collections. A collection of 5M intervals measures ≈ 113MB (text format).

Score Distribution

We conduct preliminary experiments to get an insight of the score distribution of a set of results using scored temporal predicates. Especially, we want to measure the share of results that have a high score. We want to verify that the fewer the number of high-scoring results, the less likely a worker will receive high scoring results from DTB . We compute all the combinations (x1, x2) ∈ C1 × C2 and we evaluate each result with scored predicates s-before(x1, x2), s-meets(x1, x2), s-overlaps(x1, x2) and s-starts(x1, x2).

On Figure 7, we plot scores for the top-50000 results. s-before is the predicate with the largest number of highscoring results (scores of top-50000 results equal 1.0), since a single inequality on endpoints is required. More highscoring results can be found with s-overlaps (≈ 18,000 results), that evaluates only inequality on endpoints, than with s-meets (≈ 9000), where an equality is required. Fewer results are assigned a high score when s-starts is used since it requires both equality and inequality on endpoints. Thus, we can expect a faster join processing on queries using only inequality-based predicates where more high-scoring results can be found, since high-scoring results favor early termination of local top-k processing.

Workload Distribution

We conduct experiments to validate our workload distribution approach. We analyze executions of TKIJ using DTB , our workload distribution algorithm, and using a more straightforward algorithm.

LPT. In the context of task scheduling, the LPT (Longest Processing Time) heuristic aims to minimize scheduling time on parallel machines [START_REF] Drozdowski | Scheduling for Parallel Processing[END_REF]. LPT executes tasks in descending order of processing time. In our context, a naive approach would minimize the maximum number of candidate join results that a reducer has to process, so as to reduce the duration of the longest task. With LPT , bucket combinations are analogous to tasks that we want to assign to a set of reducers. We sort the set of bucket combinations by descending order of number of results (ω.nbRes) and assign each one to the least loaded reducer.

Results. Figure 8a presents the running time of the join phase on all queries, where |Ci| varies from 1 × 10 6 to 1.6 × 10 6 . On Q b,b , running time is identical for LPT and DTB , since a single bucket combination is selected and a large number of results with maximum score can be quickly found during the join phase. On other queries, DTB outperforms LPT for two reasons. Firstly, LPT incurs a higher shuffle cost (on average 43% higher). When assigning a bucket combination to a reducer, DTB favors assignments that lessen shuffle cost. LPT favors the assignment of bucket combinations with a large number of results to the least loaded reducers. Hence, buckets have a higher probability to be sent to several reducers with LPT than with DTB . Secondly, LPT does not necessarily give a fair share of high-scoring results to each reducer. Figure 8b shows the running time of the longest reducer task (we omit Q b,b where LPT and DTB perform equally for the reason exposed above). DTB always outperforms LPT because it increases the probability that all reduce tasks terminate early since they can all find high-scoring results. This difference is exacerbated on query Q s,f,m with |Ci| = 1M . Here, the few results that satisfy best all 3 predicates featured in Q s,f,m are better distributed using DTB . On Figure 8c, we represent the minimum score of the k th result among the results returned by reducers. These results support our observation: the score of returned results is higher when distribution is defined using DTB , while unnecessary results with lower scores are returned with LPT .

b,b Q o,o Q f,f Q s,s Q s,f,m |C i |=1.6M |C i |=1.4M |C i |=1.2M |C i |=1M (
Q o ,o Q f, f Q s ,s Q s ,f ,m Q o ,o Q f, f Q s ,s Q s ,f ,m Q o ,o Q f, f Q s ,s Q s ,f ,m Q o ,o Q f, f Q s ,s Q s ,f ,m Max. Time Reducer (s) LPT DTB |C i |=1.6M |C i |=1.4M |C i |=1.2M |C i |=1M (b) Max Running Time of Reducers 0.85 0.9 0.95 1 Q o ,o Q f, f Q s ,s Q s ,f ,m Q o ,o Q f, f Q s ,s Q s ,f ,m Q o ,o Q f, f Q s ,s Q s ,f ,m Q o ,o Q f, f Q s ,s Q s ,f ,m Min. Score of k-th Result LPT DTB |C i |=1.6M |C i |=1.4M |C i |=1.2M |C i |=1M (

TopBuckets Strategies

We conduct experiments on the TopBuckets strategies exposed in Section 3.3. We vary the number of collections n using queries

Q b* , Q o* and Q m* .
Figure 9 summarizes the results. We do not report results where running time exceeds 1 hour. Experiments show the inefficiency of brute-force and two-phase. On these strategies, n-tuples of buckets, where n ∈ 3 . . . 5, are assigned a score bound by the solver. With brute-force, the running time of TopBuckets quickly increases (solid black box on Figure 9) with n, because the solver needs to compute score bounds for a large number of bucket combinations, each one requiring to assign 2n variables. The two-phase strategy only beats brute-force on Q b* (Figure 9a) where its first phase prunes a large share (more than 99% for any n) of possible bucket combinations, thus limiting the running time of the second phase, that computes exact bounds using a smaller set of bucket combinations. On others queries, two-phase does not improve running time: the first phase does not prune enough combinations (e.g. 52% on Q m* for n = 4) to lessen the cost of the second phase where remaining combinations need to be assigned tight score bounds. The loose strategy is the most efficient: (i) loose bounds do not impact significantly join processing time as a large share of potential results (e.g. 81% on Q o* for n = 4) remained pruned and (ii) TopBuckets scales with the number of collections n. In the remainder of our experiments, we use loose as the TopBuckets strategy.

Number of Granules

Since the second phase of TKIJ relies on collected statistics to prune the input space and distribute the workload, we expect TKIJ to depend on the granularity of statistics (e.g. coarse or fine-grained). We conduct experiments to validate this intuition, varying the number of granules g.

We present respectively on Figures 10a, 10b

(g=5, Q o,m ) (g=5, Q o,o ) (g=5, g>140, Q s,f,m ) Q b,b Q f,b Q o,o Q o,m Q s,f,m
(a) Running Time total running time of a range of queries, the load imbalance of the join phase computed using Max Time Reducer Average Time Reducer and the detailed running time for query Qo,m. We do not report results for executions where the total running time exceeds 1 hour. On queries that return the fewest high-scoring results (Qo,m, Q s,f,m ), we observe on Figure 10a that with a lower g (coarse statistics), running time degrades. TKIJ suffers here from poor workload distribution. This is expected as with fewer granules, we have fewer bucket combinations and especially fewer high-scoring ones. As TKIJ relies on a round-robin distribution of high-scoring buckets combinations, there is a lower probability to provide each reducer high-scoring results. Hence, we can observe on Figure 10b the imbalance that is more variable when g decreases. As illustrated on Figure 10c, when g increases, the local join processing is faster, thanks to a better workload distribution and to a larger pruning of unnecessary results. 81% of potential results are pruned for g = 20, while it is 96% for g = 100 (grey filled curve on Figure 10c). Yet, the pruning process TopBuckets is slower when g increases and thus worsen the overall response time. Note that because it features more predicates, query Q s,f,m requires to evaluate more bucket pairs during the TopBuckets process. Thus, Top-Buckets running time increases faster with g than on others queries, hence the impact on the overall response time. For queries Q b,b and Qo,o, coarse statistics have nearly no effect since a large number of high-scoring results can be found during the join phase (except when g = 5 where workload distribution fails for query Qo,o: a reducer does not quickly find high-scoring results). Finally, we observe that a number of granules g ≈ 40 provides the best trade-off for the various queries that we experimented. Future investigations include the design of benchmark approaches or the adaptation of optimization techniques [START_REF] Dignös | Overlap interval partition join[END_REF] to compute the optimal number of granules that minimizes execution time of TKIJ .

(g=5, Q o,m ) (g=5, Q o,o ) (g=5, g>140, Q s,f,m ) Q b,b Q f,b Q o,o Q o,m Q s,f,m Q s,f

Scalability

We vary the size of collections |Ci| to evaluate the scalability of TKIJ . We compare TKIJ to state-of-the-art competitors. As a baseline, we borrow algorithms from related work on processing interval joins on Map-Reduce [START_REF] Chawda | Processing interval joins on map-reduce[END_REF]. In this work, algorithms RCCIS and All-Matrix are designed to process interval joins using Allen predicates. In our settings, we use these algorithms to return only results that satisfy all the Boolean predicates of a RTJ query (i.e. a subset of top-k results). We also impose reducers to stop join processing if k results are found. Then, we merge and sort local results using a final Map-Reduce phase, identically to TKIJ . RCCIS handles only colocation predicates where intervals intersect (e.g. overlaps, meets). All-Matrix handles only sequence predicates (before, after ). RCCIS and All-Matrix also partition the temporal range using contiguous granules. For RCCIS , we set the number of granules to 24 which implies that 24 reducers are used. For All-Matrix , the number of reducers depends on the number of granules and of collections. We used 4 granules with n = 3, yielding 20 reducers. For TKIJ , we conduct a first set of experiments with 7) which incurs major pruning. That justifies the absence of a significant difference between the scored and Boolean cases.

Effect of k

We conducted experiments where k varies in [10, 10 5 ] on a range of queries (Q b,b , Qo,o, Q s,f,m , Q f,b , Qo,m) with |Ci| = 2 × 10 6 . We observed that TKIJ is almost constant on all

|C i | (x10 6 ) Q b,b Q f,b Q o,o Q o,m Q s,f,m Q jB,jB Q sM,sM
Parameters: g = 40, k = 100, P = P3, TopBuckets: loose

Network Traffic Data

Data

We use network traffic data collected on firewall logs of a data hosting company. Each log contains packets exchanged between servers and clients (≈ 5GB, 100M packets per day). Each packet has a timestamp (seconds). We selected one log and built a list of connections by grouping packets exchanged between a pair (server, client). Only consecutive packets whose timestamps are within a time interval [0, 60] are grouped. A connection [client, server, start, end ] represents the activity of client on server, where the first packet was sent or received at timestamp start and the last one at timestamp end. The dataset obtained includes 3,636,814 intervals (≈ 83MB), whose minimum, maximum, and average length are respectively 1, 86,459 and 54 seconds. Figure 12 shows the distribution of start points and lengths. Then, we copy each list of connections 3 times and process 3-way queries. We are interested in real-life scenarios occurring in network traffic analysis, hence we process queries Q jB,jB and Q sM,sM (Table 1). Query Q jB,jB returns the sequences of connections that closely follow each other, while Q sM,sM returns sequences where a delay was observed between two connections.

Scalability

We verify that TKIJ scales with various dataset sizes. When generating connections, we use various randomly selected samples on the log file used. We pick from 5% to 35% of input data. We obtain collections of connections whose number of intervals varies from 0.58 ×10 6 to 2.31 ×10 6 .

We present total running times on Figure 13. We note that running time increases faster than what was observed on synthetic data. Here, when we use larger samples on input data, we have more buckets containing at least an interval. When |Ci| = 0.58×10 6 , there are 151 buckets containing at least an interval, while there are 296 for |Ci| = 2.31×10 6 . Thus, TopBuckets has to process more bucket combinations with higher |Ci|. On query Q s,f,m , that features more predicates, the time taken by TopBuckets is dominant (e.g. 82% of overall response time for |Ci| = 1.38×10 6 ). Hence, overall response time increases faster on query Q s,f,m than on all ). Thus, we can find bucket combinations (e.g. ω = (b1, b2)) whose results (x, y) are guaranteed to have a high score s-overlaps(x, y). Then, TopBuckets returns less bucket combinations, reducing the search space while guaranteeing correctness.

b,b Q f,b Q o,o Q o,m Q s,f,m Q jB,jB Q sM,sM

Effect of k

We verify TKIJ on various values of k. We present running times on Figure 14. For all queries except Qo,o, we observe that TKIJ remains nearly constant when k ≤ 5000. Then, the running time increases slowly when k > 5000: as more results need to be returned, more intermediate results are built before termination especially with queries having fewer high-scoring results. On Qo,o we observe that TKIJ increases slightly between k = 1000 and k = 5000. That is explained by an increase (from 643 to 41,272) in the number of bucket combinations |Ω k,S | necessary to return the correct top-k results, which in turn increases the number of intermediate results.

RELATED WORK

Three research areas relate to our work, however none of them addresses the RTJ problem.

Interval Joins. The closest work to ours [START_REF] Chawda | Processing interval joins on map-reduce[END_REF] addresses the processing of multi-way joins on Map-Reduce for Allen Boolean predicates [START_REF] Allen | Maintaining knowledge about temporal intervals[END_REF]. The first algorithm, RCCIS, solves colocation queries, where all predicates require intervals to have a non-empty intersection. RCCIS reduces the amount of data shuffled in Map-Reduce by sending to the same reducer intervals that are most likely to be colocated. The second algorithm, All-Matrix handles sequence queries. Because such queries imply unavoidable replication, All-Matrix focuses on load balancing. None of those algorithms is applicable to solving the RTJ problem as they do not handle scored results. Although flexible interpretations of Allen predicates were proposed in approximate reasoning [START_REF] Dubois | Fuzziness and uncertainty in temporal reasoning[END_REF][START_REF] Ohlbach | Relations between fuzzy time intervals[END_REF][START_REF] Schockaert | Fuzzifying allen's temporal interval relations[END_REF], none of them designed join algorithms. A series of work on spatiotemporal data focused on efficiently retrieving overlapping objects with Boolean semantics. Objects are stored in parti-tions [START_REF] Dignös | Overlap interval partition join[END_REF][START_REF] Lu | On spatially partitioned temporal join[END_REF], or tree structures [START_REF] Enderle | Joining interval data in relational databases[END_REF], such as the R-Tree [START_REF] Beckmann | The r*-tree: An efficient and robust access method for points and rectangles[END_REF][START_REF] Guttman | R-trees: A dynamic index structure for spatial searching[END_REF] or the quadtree [START_REF] Finkel | Quad trees: A data structure for retrieval on composite keys[END_REF]. A more recent investigation proposed a compound index structure using segment trees to find intervals that intersect a query-interval in a key-value cloudstore [START_REF] Sfakianakis | Interval indexing and querying on key-value cloud stores[END_REF]. In summary, these studies focus on reducing I/O and support only Boolean overlaps and intersections.

Top-k Processing. Instance-optimal algorithms [START_REF] Fagin | Combining fuzzy information from multiple systems[END_REF][START_REF] Fagin | Optimal aggregation algorithms for middleware[END_REF] were proposed for rank-joins ranging from centralized implementations for HRJN [START_REF] Ilyas | Supporting top-k join queries in relational databases[END_REF] and Pull-Bound Rank Join [START_REF] Schnaitter | Evaluating rank joins with optimal cost[END_REF] in [START_REF] Finger | Robust and efficient algorithms for rank join evaluation[END_REF], to distributed ones in [START_REF] Doulkeridis | Processing of rank joins in highly distributed systems[END_REF][START_REF] Ntarmos | Rank join queries in nosql databases[END_REF]. Specifically, in NoSQL databases, the BFHM algorithm [START_REF] Ntarmos | Rank join queries in nosql databases[END_REF] places each tuple in a bucket that depends on its score. Tuples are compressed using Bloomfilters allowing to select the best buckets first then retrieve tuples to be joined from the database. Because a Bloomfilter yields false positives, BFHM may require several iterations. Our work differs in two aspects. First, our scores are predicate-dependent and are not known a priori. Second, reiterating join processing in our case would incur too high an overhead.

Load Balancing. Load balancing is a common concern in parallel top-k processing. RanKloud [START_REF] Candan | Rankloud: Scalable multimedia data processing in server clusters[END_REF] computes data statistics to retrieve an estimation of the kth join score. Then, only the part of input data whose score is above the estimated one is uniformly distributed and processed in parallel on a set of workers. Similarly, we rely on computing statistics to avoid processing useless data. However, while RanKloud outputs approximate results, our algorithm guarantees to return exact top-k results. Another work [START_REF] Doulkeridis | On saying "enough already!" in mapreduce[END_REF] proposes a partitioning scheme on Map-Reduce based on partitioning for parallel skyline query processing [START_REF] Vlachou | Angle-based space partitioning for efficient parallel skyline computation[END_REF]. The angle-based partitioning distributes evenly the volume that contains points near the best possible point, increasing the probability that skyline points will fall evenly in each partition. This idea is reintroduced in the context of top-k joins [START_REF] Doulkeridis | On saying "enough already!" in mapreduce[END_REF] and is shown to be superior to a cardinality-based partitioning. Although we share the same intuition, we cannot directly apply this technique since it assumes that partial scores are known a priori while in our case they are predicate-dependent.

CONCLUSION

In this paper, we introduce RTJ queries using scoring functions reflecting the degree of satisfaction of a join predicate. We design TKIJ , a parallel multi-way top-k interval join algorithm on Map-Reduce. TKIJ relies on an offline collection of statistics to prune the search space and a workload distribution scheme appropriate to top-k processing. We conduct experiments on synthetic data that validate our approach and show the efficiency of TKIJ on various queries. We observe the same effectiveness of TKIJ on real network traffic logs. Future investigations include the integration of interval attributes (e.g. IP address for a connection) in the join conditions, to build hybrid queries, similarly to previous work [START_REF] Chawda | Processing interval joins on map-reduce[END_REF]. We also plan to pursue investigations on various distributed platforms. Although robust and scalable, Map-Reduce suffer from limitations on complex queries [START_REF] Doulkeridis | A survey of large-scale analytical query processing in mapreduce[END_REF]. For instance, a field of investigation is the design of algorithms that take advantage of communication between workers.

Figure 1 :

 1 Figure 1: Motivating Example

Figure 2 :

 2 Figure 2: The Allen Algebra with Boolean and Scored Temporal Predicates.

1 Figure 3 :Figure 4 :

 134 Figure 3: Approximating equals and greater . Here, x ∈ {x, x}, y ∈ {y, y}.

Figure 6 :

 6 Figure 6: Example of Bucket Combinations

Algorithm 3 3 :

 33 DistributeTopBuckets (DTB) Input: Ω k,S Output: M : assignments (bucket, reducer) 1: Sort Ω k,S by descending order of ω.UB for all ω ∈ Ω k,S 4: rj = getReducer (avgRes, ω) 5:

Figure 7 :

 7 Figure 7: Synthetic Data -Score Distribution

  c) Minimum Score of k-th Result Parameters: g = 20, k = 1000, P = P2, TopBuckets : loose

Figure 8 :

 8 Figure 8: Synthetic Data -Workload Distribution

Figure 9 :

 9 Figure 9: Synthetic Data -Detailed execution time, all TopBuckets strategies

Figure 10 :

 10 Figure 10: Synthetic Data -Effect of number of granules g

Figure 11 :

 11 Figure 11: Synthetic Data -Scalability

Figure 12 :

 12 Figure 12: Network Traffic Data Distribution

Figure 13 :

 13 Figure 13: Network Traffic Data -Scalability

Param.: |Ci| = 1 .

 1 03 × 10 6 , g = 40, P = P3, TopBuckets: loose

Figure 14 :

 14 Figure 14: Network Traffic Data -Effect of k

  Algorithm 2 Strategies loose, two-phase Input: Boolean onePhase, Buckets {b i,l i ,l ′ i : i ∈ 1 . . . n} Output: Ω k,S 1: L2 ← all bucket pairs (bi, bj) s.t. (i, j) ∈ E

	2: for all ω in L2
	3:	Compute ω.UB, ω.LB using solver
	4: for ω in Ω
	5:	Compute ω.UB, ω.LB using score bounds from
		(bi, bj) ∈ L2
	6: Lm ←TopBuckets(Ω)
	7: if onePhase then return Lm
	8: for ω ∈ Lm
	9:

Table 2 :

 2 Scored Predicates Parameters

Table 2 )

 2 . We are hence using a Boolean interpretation of predicates. Because TKIJ must return k results, if only k ′ < k results satisfy the Boolean predicates (with S(t) = 1.0), k-k ′ other results that do not satisfy at least one predicate will be returned (with S(t) < 1.0). Thus, TKIJ may need to return more results than All-Matrix or RCCIS , that only return results fully satisfying all Boolean predicates. Hence, results of each algorithm are not directly comparable. For TKIJ , we conduct a second set of experiments using scored predicates with the score parameters P1.

	scoring re-

Results. We present total running times on Figure

11

. For query Q b,b (Figure

11a

), TKIJ remains nearly constant since TopBuckets returns only 1 bucket combination. Thus, TKIJ processes only a small share of input data. All-Matrix shuffles intervals belonging to all possible results, hence running time increases with |Ci|. For query Qo,o (Figure

11b

), Top-Buckets selects more combinations to process. As the number of intervals per bucket increases with |Ci|, more data is shuffled and processed during the local join processing, hence the running time increases linearly with |Ci|. Figure

11b also

shows that TKIJ outperforms RCCIS on |Ci| > 3.5×10

6

. In RCCIS, a first Map-Reduce phase builds intermediate results to determine which intervals need to be replicated to ensure output correctness in the join phase. Thus, its running time increases with |Ci|. Meanwhile, TKIJ decides which tuples should be combined on the basis of TopBuckets, which does not depend on |Ci| and is on average 93% faster than the first phase of RCCIS. On Qs,m (Figure

11c

), we do not observe this phenomenon anymore: RCCIS first phase is faster (there are fewer intermediate results), while TKIJ 's join phase is longer (there are fewer high-

no self loops and (i, j) ∈ E =⇒ (j, i) ̸ ∈ E

ACKNOWLEDGMENTS

This work was partially funded by the Datalyse PIA project.