Reporting flock patterns on the GPU

XV Spanish Meeting on Computational Geometry, June 26-28, 2013 Reporting ock patterns on the GPU Marta Fort, J.Antoni Sellarès, and Nacho Valladares ∗ Departament Informàtica, Matemàtica Aplicada i Estadística. Universitat de Girona. Abstract In this paper we study the problem of nding ock patterns in a set of trajectories of moving entities. A ock refers to a large enough subset of entities that move close to each other for a given time interval. We present a parallel approach, to be run on a Graphics Processing Unit, for reporting maximal ocks. We also provide experimental results that show the e- ciency and scalability of our approach. 1 Introduction A trajectory is a sequence of sampled time-stamped point locations describing the path of a moving entity over a period of time. Assuming that the movement of an entity between two consecutive positions is done at constant speed and without changing direction, its trajectory is modelled by the polygonal line, that can self-intersect, whose vertices are the trajectory points. In this paper we study the ock pattern [1, 2, 3], that identies a group of entities moving close together during a given time interval (see Figure 1). e2 e3 e0 e1 x y tτ−1 D1 D2 t0 Figure 1: Example of two ocks. The increasing programmability and high computational rates of Graphics Processing Units (GPUs), together with CUDA (Compute Unied Device Architecture) and some programming languages which use this architecture, make them attractive to solve problems which can be treated in parallel as an alternative to CPUs. GPUs are used in dierent computational tasks where a big amount of data or operations have to be done, whenever they can be processed or done in parallel. Some recent works show that demanding algorithms in dierent elds can take advantage of the GPU parallel processing [4, 5, 6]. ∗ Email:(mfort,sellares,ivalladares)@imae.udg.edu. Authors research supported by TIN2010-20590-C02-02. In this paper, we present an ecient parallel GPU- based algorithm, designed under CUDA architecture, for reporting maximal ocks. The experimental results obtained with the implementation of our algorithm show the signicance of the presented approach according to its performance and scalability. 2 The ock pattern Let E={e0, . . . , en−1} be a set of n moving entities. The trajectory Ti of the entity ei is a sequence of τ points in the plane Ti:ei 0, . . . , ei τ−1 , where ei j denotes the position of ei at time tj with 0≤j < τ . We assume that the positions are sampled synchronously for all the entities, and that entity ei moves between two consecutive positions ei j , ei j+1 with constant speed and without changing direction. Consequently, the trajectory Ti is described by the polygonal line, which may self-intersect, whose vertices are the trajectory points. Flocks identify groups of entities whose trajectories are close together during a minimum period of time. Formally, according to [2], given a set E of entities, a minimum number of entities µ∈N , a number of time-steps δ∈N , a time interval Iδ j= [tj, tj+δ−1] , and a distance ϵ∈R : Denition 1 A ock f(δ, µ, ϵ) in a time interval Iδ j , consists of at least µ entities such that for every discrete time step tℓ∈Iδ j , there is a disk of radius ϵ that contains all the µ entities. The number of reported or determined ocks is a critical issue that can adversely aect the response time and the proper interpretation of the nal results. For a time-step tℓ∈Iδ j there may be several circles with radius ϵ that yield a ock with the same subset of entities of E , so we consider that two ocks are dierent if they involve two dierent subsets. Since there can be Θ(n2) combinatorially distinct ways to place a circle of radius ϵ among n points in the plane, at each timestep there can be Θ(n2) ock candidates involving different subsets of entities [9]. An algorithm having as output the position of the entities involved in a ock, would have an output size of Θ(n3) per time-step. Moreover, it is easy to observe that a set of entities could be present in many ocks, and even one single entity can be involved in several ocks. For example, a ock of µ+1 entities contains µ+1 ocks of µ entities. 19 Reporting ock patterns To overcome this problem, as in [3], instead of nding all the existent ocks, we are interested only in maximal ocks, so that the entities of one ock may not be a subset of the entities of another ock. Even in the case of nding maximal ocks, at every time-step there may still asymptotically be Θ(n2) ock candidates, although in many real situations the number of candidates decreases considerably. Denition 2 A maximal ock mf(δ, µ, ϵ) in the time interval Iδ j , is a maximal ock f(δ, µ, ϵ) in Iδ j with respect to the subset inclusion. We will denote Mδ j the family of maximal ocks in the time interval Iδ j , and Mδ the family of all maximal ocks, that is, the ocks of Mδ j , for j= 0, . . . , τ −1 . 2.1 Characterization Vieira et al. [3] prove that, for each pair of entity locations, there exists two such disks. Thus, 2n2 disks per time step must be tested. Since we are interested in those ocks that are maximal, we prove that keeping just one of the two disks is enough. Consequently, only n2 disks per time step must be tested, thus the number of tests is reduced by half. The following lemma, similar to one presented in [3], allows us to bound the number of disks to be tested in each time step. Lemma 3 Let P be a set of n points and D′ be a disk of radius r that covers a subset P′⊆P , |P′|=k≤n . There exist another disk D′′ of radius r such that: a) has at least two points pi, pj∈P′ on its boundary; b) covers a superset P′′ of P′ , P′⊆P′′ ⊆P ; c) it is located at the right side of points pi, pj . 3 Reporting ock patterns In this paper study the problem of reporting the family Mδ of all maximal ocks. To report Mδ , we must nd the family Mδ j in the time interval Iδ j , for j= 0, . . . , τ −1 . 3.1 Computing the family Mδ j Next, we describe the two main steps of the algorithm used to compute the family Mδ j in the time interval Iδ j . First step. At each time step ti in the time interval Iδ j , we compute the family Fi of potential ocks containing at least µ entities. Each family Fi is obtained by using the characterization given in Lemma 3. Next, we nd the maximal sets of each family Fi . Abusing of notation, we still denote Fi the obtained subfamily of potential maximal ocks. At the end of the step, we have the array of families of potential maximal ocks Fδ j= [Fj, . . . , Fj+δ−1] in the time interval Iδ j . Second step. We rst compute the family Iδ j obtained intersecting the δ families of potential maximal ocks in Fδ j . That is, Iδ j contains the sets which are the intersection of δ sets, one of each family Fi in Fδ j . Finally, we nd the subfamily of maximal sets of Iδ j that contain at least µ elements. The obtained subfamily is the desired Mδ j . 4 GPU implementation In this Section, we describe how to implement the dierent steps of the algorithm that computes Mδ j , when GPU parallel techniques are used. 4.1 The multi-grid structure G We denote by Ei={e0 i, . . . , en−1 i} the set of locations of the entities of E at time step ti , i∈ {0, . . . , τ − 1} . In [3], a single grid is used to perform quick disk range searching queries over each set Ei . Instead of using a grid for each time step, we build a multi-grid structure, denoted G , in parallel in the GPU, that allows us to search simultaneously in each one of the δ sets of Eδ j= [Ej, . . . , Ej+δ−1] . The multi-grid structure contains δ regular grids, where the edges of the grid cells have length ϵ . Each grid j contains the set of points corresponding to Ej . In order to maximize the parallel performance, the top-right corner is dened so that the grid is a square, that is, it has the same number of rows and columns (see Figure 2). The multi-grid structure G is composed of 7 arrays allocated in the GPU global memory. In Figure 2 we can see an example for δ= 3 . The bottom-left and the top-right corners of each grid are stored in the array bδ of size δ , where bδ[j] stores the two corner points of the grid containing Ej . The number of cells per row (or per column) of each grid is stored in ecδ and its prex sum is stored in epδ . Both arrays are integer arrays of size δ . To represent the δ grids, we use the array of integers cδ , whose size depends on the size of the dierent grids which directly depend on the distribution of the points. For a given grid cell, we store in the corresponding cδ position the number of points contained on that cell. The array pδ of the same size of cδ is obtained as the prex sum of cδ . Note that, for a better understanding, in Figure 2 the arrays cδ and pδ are showed as several 2D grids, but in the GPU memory they are stored one 1D array each, by storing one grid after the other starting from left to right and from top to bottom, thus, the rst grid corresponds to the time step ti and the last one of the time step ti+δ . With this structure we can easily access, in a parallel way, to any element of the δ grids, by using the fact that the set of points of a cell w starts at position pδ[w] of eδ , and is stored in cδ[w] consecutive positions of eδ . In each eδ position we store the two oats of the 20 XV Spanish Meeting on Computational Geometry, June 26-28, 2013 corresponding location. Finally, eeδ is an integer array of size nδ containing at position j the entity id of the eδ[j] location. cδ pδ t0t1t2 2 1 11 1 3 0 0 0 0 0 0 0 0 0 2 2 2 3 3 4 5 5 5 5 6 6 6 p0 0 eδp1 0p2 0 t0 ecδ epδ 231 025 t0t1t2 0 eeδ21 p0 1 p1 1p2 1p1 2p2 2p0 2 021 021 t1t2 bδb0b1b2 ǫ e0 e1 e2 ǫ Figure 2: Example of multi-grid G for δ= 3 . We construct the G structure in parallel as follows. We allocate CPU arrays ecδ′ , epδ′ , b′ δ of size δ , and p′ and eeδ′ of size nδ . While we load the input P into p′ the arrays ecδ′ , epδ′ and b′ δ are accordingly lled. Then ecδ , epδ , bδ and p are allocated in GPU memory and the ecδ′ , epδ′ , b′ δ and p′ values are transferred to GPU memory to ecδ , epδ , bδ and p respectively. Then, cδ and pδ are allocated in the GPU memory according to ecδ and epδ values. Additionally, the arrays eδ and eeδ are allocated with size nδ also in GPU memory. Because GPU has no dynamic memory we have to construct G in two steps. First, we ll cδ by counting the number of points contained in each cell. This is done in parallel by launching a kernel with nδ threads, that is, a thread per point and per time step. Each thread idx reads its p[idx] value and determines its position in the grid by using its position on the plane, the grid corresponding to the time step which the point belongs to, and the ecδ and epδ arrays. The corresponding cδ position is incremented in 1 . Because many threads may correspond to the same cδ position we use atomic operations to ensure no thread interferences. Once cδ is computed, pδ is obtained as the prex sum of cδ and we allocate an auxiliary array v with the same size of cδ , initialized to zeros. Now, we launch a parallel kernel with nδ threads where each thread idx reads its p[idx] value, determines its position w in the grid and stores p[idx] at eδ[pδ[w] + v[w]] . Every time a thread stores its associated point into eδ the corresponding v value is incremented in one by using an atomic operation. 4.2 Finding potential ocks The computation of the array of potential ocks Fδ j of a given interval Iδ j is done in two main steps. First, we need to nd the center of the disks which are obtained applying the characterization given in Lemma 3, and second, we need to report the entities contained in those disks. Additionally, because GPUs do not have dynamic memory, we split the process in several steps. First, we allocate an array of integers Dc of size nδ to count how many disks there exist. This is done in parallel by launching a kernel with nδ threads, that is, a thread per point within Iδ j . Using the structure G each thread idx goes through its neighbors and counts the number of points which are at distance ≤2ϵ . The results are stored in Dc[idx] . Because this is computed in parallel, each disk will be reported twice, once for each of the two points dening the disk. In order to avoid this, we force threads to check those points whose entity id is bigger than itself. Then, Dp of size nδ , is computed as the pre- x sum of Dc , and the array D is allocated with size Dp[nδ−1]+Dc[nδ−1] , where we will store the centers of disks. The same process is repeated, but, instead of counting, each thread reports the center of the disk. Let us recall that only the right disk per paired point has to be reported. In the second step, we actually report the potential ocks. To do this we use G to perform the range searching queries, using the values of the array D as the center of the queries with radius ϵ . Note that, because in order that the regular grid structure works for range searching queries, we must be able to locate any point within the grid. It may happen that some disk centers were placed outside the grid for a given time step. Thus, before we perform the range searching queries we have to reconstruct G so that the center of the disks were contained in G . We could modify G , but it is faster to rebuild it than to modify it. To count the number of points contained on each disk we rst allocate Pc , the array of integers of size equal to the total number of disks reported which is Dp[nδ−1]+Dc[nδ−1] , initialize it to zeros and launch a kernel with one thread per disk. Each thread idx reads the center of its disk D[idx] , locates D[idx] on the grid and checks the distance between the center of the disk and the points of its cell neighbors. Each thread counts how many points are at distance ≤ϵ to its disk center, and if there are at least µ the number is stored in Pc[idx] . Finally, we compute Pp as the prex sum of Pc . If there is at least one disk per time step with at least µ entities, we continue to report the potential ocks. The last position of Pp plus the last position of Pc gives us ω , the size of the array used to store the potential ocks. Thus, the array P is allocated as an integer array with size ω and is lled in parallel using a thread per disk. The process is the same as before but this time, instead of counting, each thread stores in P the entities contained inside each disk. Note that, Pc , Pp and P denote a structure contain- 21 Reporting ock patterns ing δ families of sets. This is the array of potential ocks Fδ j , that is, δ families of sets containing the potential ocks of δ consecutive time steps starting at time step tj . We denote Fj the family of potential ocks of the time step tj . The last step is to reduce each Fi∈ Fδ j so that it only contains the sets which are maximal, with respect to the partial order induced by the subset relation in each family. To this end, we use the parallel algorithm of Fort et al. [7]. When the process nishes, we have in Fδ j the maximal potential ocks with µ or more entities at each time step of the interval Iδ j . 4.3 Reporting Mδ ocks We want to report all the ocks Mδ j for all Iδ j, j = [0, . . . , τ −δ] . Thus, we could perform the same algorithm over each Iδ j , but many information computed for a given time step within a time interval can be reused for the following time intervals. The idea is to compute the potential ocks for a given interval and compute their intersections in a way that we can reuse the information for the next time interval. To this aim we proceed as follows. We start at the interval Iδ 0 and we compute Fδ 0 . Then, in order to obtain Mδ 0=Fδ−1∩Fδ−2∩···∩F0 we start at Fδ−1 and we intersect Fδ−2∩Fδ−1 . Next, we intersect Fδ−3∩Fδ−2∩Fδ−1 and so on. After each intersection, we compute the maximal sets for each resulting intersection discarding those sets which are non maximal or containing less than µ entities. When we are done, we add, to Mδ , the family Mδ 0 , which contains the ocks of Iδ 0 . For the following time steps, instead of recomputing all the potential ocks and all the intersections, we reuse the previous computations. If for a given time step tj , a family of potential ocks Fk, j ≤k < j +δ is not inside the array Fδ j , we compute Fδ j+δ−1 . That is, instead of computing the potential ocks of one single time step we compute them in δ intervals. Then, we perform the corresponding intersections. This technique does not decrease the number of intersection we have to do, but, in practice, the δ families we intersect in each step are smaller than the ones containing the potential ocks, leading to faster results. We repeat the process until j=τ−δ . The intersection process is also performed in parallel using the GPU algorithm presented by Fort et al in. [8]. 5 Experimental results The experimental results are obtained using an Intel Core2 CPU 6400 with a Nvidia GTX 480. We split the running time into 2 main steps, the potential ocks reporting process and the intersection process. The accumulated value, corresponding to the columns height, gives the total running time of the algorithm. The algorithm has been tested with a data set extracted from [10], that consists of 145 trajectories of 2 school buses around Athens metropolitan area in Greece with a total of 66,096 time steps. The results (Figure 3) shows that when we vary ϵ while we maintain µ= 5 and δ= 10 , the number of ocks reported increases as long as we increase ϵ . The most aected part is the intersection process. This is because the more ocks we report the more intersections we have to compute. When the parameter µ is varied, while ϵ= 1200 and δ= 10 , the condition to form ock is more restrictive. Thus, the number of ocks decreases and, consequently, the running times too. Figure 3: Running times varying ϵ and µ . References [1] J. Gudmundsson, M. J. van Kreveld, B. Speckmann, E- cient Detection of Patterns in 2D Trajectories of Moving Points, GeoInformatica 11 (2) (2007) 195215. [2] M. Benkert, J. Gudmundsson, F. Hübner, T. Wolle, Reporting ock patterns, Computational Geometry 41 (3) (2008) 111125. [3] M. R. Vieira, P. Bakalov, V. J. Tsotras, On-line discovery of ock patterns in spatio-temporal data, in: D. Agrawal, W. G. Aref, C.-T. Lu, M. F. Mokbel, P. Scheuermann, C. Shahabi, O. Wolfson (Eds.), GIS, ACM, 2009, pp. 286 295. [4] J. D. Owens, D. Luebke, N. Govindaraju, M. Harris, J. Krüger, A. E. Lefohn, T. J. Purcell, A survey of general-purpose computation on graphics hardware, Computer Graphics Forum 26 (1) (2007) 80113. [5] N. Coll, M. Fort, N. Madern, J. A. Sellarès, Multi-visibility maps of triangulated terrains, International Journal of Geographical Information Science 21 (10) (2007) 11151134. [6] M. Fort, J. A. Sellarès, N. Valladres, Computing popular places using graphics processors, in: Proc. SSTDM'10 in cooperation with IEEE ICDM'10, IEEE Computer Society, 2010, pp. 233241. [7] M. Fort, J.A. Sellaès, N. Valladares, Finding extremal sets on the GPU (Submitted). [8] M. Fort, J.A. Sellaès, N. Valladares, Intersecting two families of sets on the GPU (Submitted). [9] J. Gudmundsson, M. van Kreveld, B. Speckmann, Ecient Detection of Motion Patterns in Spatio-Temporal Data Sets, in: D. Pfoser, I. F. Cruz, M. Ronthaler (Eds.), GIS, ACM, 2004, pp. 250257. [10] Y. Theodoridis, R-Tree portal, http://www.rtreeportal.org (2011). 22