Finding influential location regions based on reverse k-neighbor queries

Knowledge-Based Systems 47 (2013) 35–52

Contents lists available at SciVerse ScienceDirect

Knowledge-Based Systems journal homepage: www.elsevier.com/locate/knosys

Finding influential location regions based on reverse k-neighbor queries Marta Fort ⇑, J. Antoni Sellarès Universitat de Girona, Informàtica i Matemàtica Aplicada, Campus Montilivi, Edifici P4, E-17071 Girona, Spain

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 7 September 2012 Received in revised form 20 March 2013 Accepted 27 March 2013 Available online 8 April 2013

In this paper we introduce and solve several problems that arise in the single facility location field. A reverse k-influential location problem finds a region such that the location of a new facility, desirable or obnoxious, in the region guarantees a minimum k-influential value associated to the importance, attractiveness or repulsiveness, of the facility as a solution to a reverse k-nearest or farthest neighbor query. Solving reverse k-influential location problems help decision makers to progress towards suitable locations for a new facility. We present a parallel approach, to be ran on a graphics processing unit, for approximately solving reverse k-influential location problems, and also provide and discuss experimental results showing the efficiency and scalability of our approach. Ó 2013 Elsevier B.V. All rights reserved.

Keywords: Computer science Decision-making support systems Reverse k-nearest/farthest neighbor query Facility location Graphics Processing Unit (GPU)

1. Introduction Facility location analysis helps decision makers to identify optimal locations for facilities according to some predefined criteria. Facility locations selection is often a critical and difficult problem to solve for many organizations. The main objective of a single facility location problem is to place a new facility with respect to a given set of sites such that certain optimality criterion is satisfied [13,12,23]. Facilities can be attractive or desirable, such as supermarkets, schools or hospitals, that customer want them as close as possible, or obnoxious or undesirable, such as nuclear power plants, garbage dumps or chemical plants, that customers wish to push them as far as possible. The optimality criterion of most classical single facility location problems is based on distances. Many variations have been considered. For desirable facilities typically the maximum distance or the sum of distances from the facility to the sites is minimized. For obnoxious facilities, the minimum distance or the sum of distances from the facility to the sites is maximized. An inherent limitation of classical single facility location problems is that only one optimal location is returned as the answer, which in many cases is very expensive to compute and in general is also too restrictive, there possibly exist sub-optimal solutions that could be considered even more appropriated by the experts, than the optimal one.

⇑ Corresponding author. Tel.: +34 972 41 82 42; fax: +34 972 41 84 00. E-mail addresses: (J.A. Sellarès).

[email protected]

(M.

Fort),

[email protected]

0950-7051/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.knosys.2013.03.013

In practice the different importance of the sites, either attractiveness or repulsiveness, must be taken into account. A weight is assigned to each site to represent its importance. Experts take available information about different factors related to the site, for instance the daily mean number of customers of a cinema or the buying power of the people living in a residential zone, and then aggregate these factors in a unique weight [10,11]. Although the proximity or remoteness of the new facility to the existing sites is an indicator of its potential impact, during the process of deciding the location for the new facility it is also very important to take into account the ‘‘influence’’ of this location among its, nearest or farthest, neighbor sites. Reverse neighbor queries model, in a natural way, how the potential addition of a new facility influences its neighbors sites [18,19,31]. The solution sites of a k-reverse nearest query are the most highly influenced sites by a query facility, by considering their k-nearest neighbors. Meanwhile the solution sites of a reverse farthest query are the most poorly influenced ones by a query facility, by being their k-farthest neighbors.

1.1. Reverse k-influential location problems The k-influential value of a given facility is defined as the sum of the weights, associated to the importance, of the reverse k-neighbor sites of the facility. The objective of a reverse k-influential location problem is to find a region such that the location of a new facility in this region guarantees a minimum k-influential value. Such a minimum value should be chosen by experts according to the problem being solved. Observe that reverse k-influential location problems take

36

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

into account the influence capacity of the new facility among its weighted neighbor sites and return a whole region of suitable sub-optimal locations. Next we present several motivational examples. 1.1.1. Monochromatic case Example 1. If we are interested in locating a new center of cultural interest (museum, zoo, thematic park), to help it attracting visitors, we will consider that it is best located at points having a 3-nearest influential value, with respect to the centers of cultural interest of the region, of at least 10. Example 2. If we have to locate a new garbage plant, to best distribute this kind of undesirable facilities, it will be best located at points having a 4-farthest influential value, with respect the garbage plants of the region, that is at least 72. 1.1.2. Bichromatic case Example 3. If we want to locate a new restaurant, guaranteeing a minimal number of potential customers, the new restaurant will be best located at points having a 5-nearest influential value of at least 15, with respect to the set of cinemas of the area. Example 4. If we are interested in locating a new residential zone, to assure the comfort of the future owners, we will consider that the best location corresponds to points having a 5-farthest influential value, with respect the set of industrial zones of the region, that is at least 93. Reverse k-influential location problems can be combined with classical single facility location problems as it is shown in the following motivational examples. 1.1.3. Monochromatic case Example 5. In the case of Example 1, where we search for a suitable region to locate a new center of cultural interest according to a minimum 3-nearest influence value with respect to the set of cultural centers, we can additionally be interested in locating the new center in the point of this region so that the maximum distance or the sum of the distances to the set of its reverse 3-nearest neighbor centers is minimum. Example 6. For the Example 2, where we want to find a region where to locate a new garbage plant considering a given minimum 4-farthest influence value with respect to the set of garbage plants, we can locate the new plant in the point of the suitable region so that the minimum distance or the sum of distances to the set of its reverse 4-farthest neighbor garbage plants is maximum. 1.1.4. Bichromatic case Example 7. For the Example 3, where we want to find a region where to locate a a new restaurant according to a minimum 5-nearest influence value with respect to the set of cinemas, in addition we locate the restaurant in the point of this region so that the maximum distance or the sum of distances to the set of its reverse 5-nearest cinema neighbors is minimum. Example 8. In the case of Example 4, where we search for a suitable region where to locate a new residential zone considering a minimum 5-farthest influence value with respect to the set of industrial zones, we will also locate the new residential zone in a

point of the suitable region so that the minimum distance or the sum of distances to the set of its reverse 5-farthest industrial zone neighbors is maximum. 1.2. Graphics processing units The increasing programmability and high computational rates of Graphics Processing Units (GPUs) make them a compelling platform for computationally demanding tasks, where a big amount of data or operations have to be done. The parallel processing capability of the GPU allows it to divide complex computing tasks into thousands of smaller tasks that can be run concurrently. This ability is enabling researchers to address many challenging computational problems up to several orders of magnitude faster than usual CPUs. GPUs have quickly become an industry standard that power millions of desktops, notebooks, workstations and supercomputers around the world. Compute Unified Device Architecture (CUDA) is a parallel architecture for computing on NVIDIA’s GPUs, which is accessible to software developers through variants of industry standard programming languages, like CUDA C or OpenCL. These general purpose languages give access to CUDA architecture avoiding a graphics context. Since its introduction, CUDA has been widely deployed through a large number of applications and published research papers. General-Purpose computing on the GPU (GPGPU) is capturing, to drastically decrease execution times, the attention of researchers in many computational fields which range from numeric computing operations and physical simulations to bioinformatics, data mining or geometric processing [25,22,33,34,32]. 1.3. Our contributions In this paper we introduce and approximately solve several kinfluential location region problems, and also combined reverse k-influential location and classical distance based single facility location problems. By working towards practical solutions, we use Computational Geometry techniques together with the GPU parallel computing capabilities to solve the problems. Additionally, to improve the understanding of the proposed problems, we have designed an interface that visualizes the computation and visualization facilitates decision makers, with an iterative what-if-analysis process in which some parameters could obtain suitable locations for a new facility. The short response times for computing and visualizing the results facilitate the interaction with the user. Finally, we provide and discuss experimental results obtained with the implementation of our algorithms that show the efficiency and scalability of our approach. 1.4. Organization of the paper The rest of the paper is organized as follows. In Section 2 we summarize the necessary background to understand the paper. Section 3 formally introduces influential location problems together with related work. The approach for solving influential location problems on the GPU, for the monochromatic and bichromatic cases, is described in Sections 4 and 5. Sections 6 and 7 present the visualization of the results and provide the experimental results obtained with the presented algorithms. Finally, conclusions and farther comments are given in Section 8. 2. Background In this section we start by formally defining neighbor and reverse neighbor queries. Next several important aspects of CUDA

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

architecture and parallel programming paradigms are discussed. Finally existent GPU-based algorithms providing fast solutions to some multiple proximity queries are summarized. 2.1. k-Neighbor and reverse k-neighbor queries k-Nearest neighbor queries and monochromatic and bichromatic reverse k-nearest neighbor queries, and their farthest variants, have attracted considerable attention in recent years [17–19,29,31]. We denote F a set of n facilities and C a set of m customers lying on an axis parallel rectangular region of interest R in the Euclidean plane, and d(p, q) represents the Euclidean distance between points p and q. 2.1.1. k-Neighbor queries  A k-nearest neighbor (k-NN) query finds the subset NNk(F, q) of k closest facilities to the query q 2 R:

NNk ðF; qÞ ¼ F 0  F;

jF 0 j ¼ k;

such that df 0 ðqÞ < df ðqÞ; 8f 0 2 F 0 ; f 2 F  F 0 .  A k-farthest neighbor (k-FN) query finds the subset FNk(F, q) of k farthest facilities to the query q 2 R:

FNk ðF; qÞ ¼ F 0  F;

jF 0 j ¼ k;

such that df 0 ðqÞ > df ðqÞ; 8f 0 2 F 0 ; f 2 F  F 0 . The points of NNk(F, q) are the k points that have the highest influence on q, and the points of FNk(F, q) are the k points that have least influence on q. To find the two nearest hotels to a specific location (see Fig. 1a), we solve a 2-NN query. Meanwhile if we are interested in finding the five farthest apartment blocks from a given school, we are solving a 5-FN query. 2.1.2. Monochromatic reverse k-neighbor queries  A monochromatic reverse k-nearest neighbor (R-k-NN) query retrieves the subset RNNk(F, q) of facilities of F that have the query point q 2 R in the set of their k-nearest neighbors in F [ {q}:

RNNk ðF; qÞ ¼ ff 2 Fjq 2 NN k ðF [ fqg; f Þg:  A monochromatic reverse k-farthest neighbor (R-k-FN) query determines the set RFNk(F, q) of facilities in F that have the query point q 2 R in the set of its k-farthest neighbors in F [ {q}:

RFN k ðF; qÞ ¼ ff 2 Fjq 2 FNk ðF [ fqg; f Þg: The points of RNNk(F, q) are the most highly influenced by q, and the points of RFNk(F, q) are the most poorly influenced by q. Solving a R-4-NN query can be useful to find the hotels that are most affected by the opening of a new hotel because it is between their 4-nearest neighbors (see Fig. 1b). Concerning the R-3-FN

Fig. 1. Given a set of hotels H = {h1, h2, h3}: (a) NN2(H, q) = {h1, h2} and (b) RNN2(H, q) = {h1}.

37

query we can be interested in finding the garbage plants that have the new one among the 3-farthest neighbors. The number of reverse k-nearest neighbors of a query point is at most 6k [3], and no known upper bound exists on the number of reverse k-farthest neighbors of a query point [19]. 2.1.3. Bichromatic reverse k-neighbor queries  A bichromatic reverse k-nearest neighbor (BR-k-NN) query for a point q 2 R seeks to determine the set BRNNk(F, C, q) of customers c 2 C for which q is a k-nearest neighbor in F [ {q}:

BRNN k ðF; C; qÞ ¼ fc 2 Cjq 2 NN k ðF [ fqg; cÞg:  A bichromatic reverse k-farthest neighbor (BR-k-FN) query for a point q 2 R seeks to determine the set BRFNk(FR, C, q) of customers c 2 C for which q is a k-farthest neighbor in F [ {q}:

BRFNk ðF; C; qÞ ¼ fc 2 Cjq 2 FN k ðF [ fqg; cÞg: The points of BRNNk(F, C, q) are the most highly influenced by q, and the points of BRFNk(F, C, q) are the most poorly influenced by q. Solving a BR-2-NN we can find the cinemas that have the new restaurant among their 2-nearest restaurants to determine the number of potential customers of the restaurant (see Fig. 2). On the other hand, to find the industrial zones that have the new residential zone among their 4-farthest neighbor, we have to solve a BR-4-FN. The number of bichromatic reverse k-nearest neighbors of a query point may be arbitrarily large [27]. Meanwhile, there is no known upper bound on the number of bichromatic reverse k-farthest neighbors of a query point [19]. From now on, we skip the word nearest or farthest if a definition or statement holds for both criteria. Thus we talk about monochromatic or bichromatic k-neighbors or monochromatic or bichromatic reverse k-neighbors without specifying the nearest or farthest criterion. In this sense, we refer to k-N, R-k-N and BRk-N queries and Nk(F, q), RNk(F, q), BRNk(F, C, q) subsets. When needed the used criterion will be specified. 2.2. CUDA architecture GPU architectures contain many processors which can execute instructions in parallel. A program executing on the GPU using CUDA architecture is decomposed in many parallel independent threads. Threads are organized into grids, of user defined sizes, and grids into blocks with a user specified number of threads. During the parallel execution each thread executes the instructions contained in the so called kernels. Each thread computation is independent from the others, however, there exist some readmodify-write operations called atomic functions. They read and return the value stored in a memory position, operate on it and store the result without allowing, during the whole process, any other access to that memory position. Thus when such functions

Fig. 2. The 2-nearest cinema neighbors of restaurant r are c3 and c4. The reverse 2nearest cinema neighbors of a new restaurant r are c1 and c2.

38

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

are used the operations done in a specific memory position are not done in parallel but serialized. There exist atomic functions to accumulate integer values and to maintain the maximum or minimum value. CUDA allow us to synchronize the execution of all the threads in a block, they wait until all the threads of a block have finished the part of the kernel code before the synchronization call, and next they keep on executing the following instructions. Several types of memories can be used. Data stored in global memory are accessible by every thread and are visible from the CPU. Global memory is where more data can be allocated, but is the slowest access time memory. Shared memory is only available for all threads in a block, is fast but not accessible from CPU. Finally, registers are the fastest memory on the GPU and store the local variables of each thread. Since the number of accesses to global or shared memory are reflected in the execution times of the algorithms, we represent by rgm ; wgm and t gm the time needed to read, write and transfer from CPU m values to global memory. When shared memory is used they are denoted rsm ; wsm , meanwhile tsm makes no sense because information cannot be transferred from CPU to shared memory. Finally we denote d the extra time needed by the atomic operations. As pointed out in Lieberman et al. [20], serial and parallel programming paradigms are focused on very different purposes making efficient serial algorithms unattractive candidates for parallel machines. When working with serial algorithms the total number of operations have to be minimized to guaranty time-efficiency, they have to be work-efficient. When working in parallel, workefficiency does not guaranty time-efficiency, elaborated algorithms tend to use complex data structures that cannot be built nor managed in parallel machines like GPUs, this makes that in some cases naive algorithms can be the best option. 2.3. Proximity queries In a proximity query we want to preprocess a point set into a data structure so as to efficiently answer questions based on distances. Given a set S of n points in Euclidean plane, basic multiple proximity queries are:

range searching are used and generalized to solve the farthest problem in the present work. Next they are briefly described. Multiple disk range searching queries are solved using the grid containing S, and considering a thread per disk. In a first step each thread locates the disk center in the grid and explores the grid cells intersecting the corresponding disk solving a range counting query. The number of points of S contained inside the disk for the nearest case, or outside it in the farthest one is determined to allocate enough space in the GPU to store the total output. The total output size is determined, and finally the range. To solve monochromatic and bichromatic all k-neighbor queries a thread per query is used. The thread locates its corresponding query point in the grid structure, and explores the grid cells in a spiral like search. The search works outwards starting at the cell containing q when the nearest case is handled, and looking inwards starting in the grid boundary, for the farthest case. 3. Influential location problems In this section we formalize the concept of the influential value of a point and formulate several influential location problems. 3.1. k-Influential value of a query point During the process of deciding the location for a new facility it is very important to take into account the importance of the sites facilities or customers, that we model assigning weights, and the influence of the new location among its neighbor sites, that we capture with the solutions of reverse k-neighbor queries. 3.1.1. Monochromatic case We assign a positive-integer weight w(f) to each facility f 2 F to model its importance.  The k-influential value of point q 2 R, Ik(F, q), is the sum of the weights of the facilities of RNk(F, q):

X

Ik ðF; qÞ ¼

wðf Þ:

f 2RNk ðF;qÞ

 Multiple disk range searching queries. Let D be a family of disks. A disk range searching query for a query disk D 2 D counts or reports the points in S \ D. In a multiple disk range searching query, instead of a unique range D 2 D, we have a subfamily D0 of ranges of D and we want to solve a disk range searching query for each range D 2 D0 .  All k-neighbor queries. Let Q be a set of m query points. A bichromatic all k-nearest/farthest neighbor query finds the subset of k points of S which are nearest/farthest to each point q of Q. The case where Q = S is known as monochromatic all k-neighbors problem.

 The k-influential map, MapIk(F), is a function that maps every point q 2 R to its k-influential value Ik(F, q). 3.1.2. Bichromatic case We assign a positive-integer weight w(c) to each customer c 2 C to model its importance.  The k-influential value of point q 2 R, Ik(F, C, q), is the sum of the weights of the customers of BRNk(F, C, q):

Ik ðF; C; qÞ ¼

X

wðcÞ:

c2BRNk ðF;C;qÞ

2.3.1. Solving multiple proximity queries on the GPU The uniform grid is one of the simplest data structures used to solve proximity queries between other geometric problems. It has been experimentally shown that for many applications uniform grids are efficient for both evenly and unevenly spaced data and moreover they are appropriate for parallel processing [1]. For a given set of points S, Leite et al. [21] build and store an uniform grid over a bounding box of S on the GPU by using CUDA, and then use it to solve k-nearest neighbor queries. Later, Fort and Sellarés [14] used the same grid structure to solve multiple disk range searching queries finding the points contained inside a disk. These approaches to find the k-neighbor and to solve the disk

 The k-influential map, MapIk(F, C), is a function that maps every point q 2 R to the k-influential value Ik(F, C, q). To indicate if we refer to k-near or k-far influential values we will substitute Ik by NearIk or FarIk, respectively. 3.2. k-Influential location problems We describe a problem, for the monochromatic and bichromatic cases, that tries to guarantee a minimum influential value of a new facility on the set of its neighbors.

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

Monochromatic case  RegIk(F, I0) problem computes the region RegIk(F, I0) where a new facility q should be placed, so that its k-influential value Ik(F, q) is at least I0. Bichromatic case  RegIk(F, C, I0) problem computes the region RegIk(F, C, I0) where a new facility q should be placed to have its influential value Ik (F, C, q) at least I0. 3.2.1. Combining influential location and single facility location problems In this section we describe a family of problems that are a combination of the location problems presented previously with classical single facility location problems. 3.2.1.1. Monochromatic case.  MinMaxNeark(F, I0) problem locates a new facility in RegNearIk (F, I0) so that the maximum distance to the set of its reverse k-nearest neighbors is minimum.  MinSumNeark(F, I0) problem locates a new facility in RegNearIk (F, I0) so that the sum of the distances to the set of its reverse k-nearest neighbors is minimum.  MaxMinFark(F, I0) problem locates a new facility in RegFarIk(F, I0) so that the minimum distance to the set of its reverse k-farthest neighbors is maximum.  MaxSumFark(F, I0) problem locates a new facility in RegFarIk(F, I0) so that the sum of distances to the set of its reverse k-farthest neighbors is maximum.

3.2.1.2. Bichromatic case.  MinMaxNeark(F, C, I0) problem locates a new facility in region RegNearIk(F, C, I0) so that the maximum distance to the set of its reverse k-nearest customer neighbors is minimum.  MinSumNeark(F, C, I0) problem locates a new facility in region RegNearIk(F, C, I0) so that the sum of distances to the set of its reverse k-nearest customer neighbors is minimum.  MaxMinFark(F, C, I0) problem locates a new facility in region RegFarIk(F, C, I0) so that the minimum distance to the set of its reverse k-farthest customer neighbors is maximum.  MaxSumFark(F, C, I0) problem locates a new facility in region RegFarIk(F, C, I0) so that the sum of distances to the set of its reverse k-farthest customer neighbors is maximum. 3.3. Related work There are some previous papers closely related with the study of reverse influential locations. Cabello et al. [4] show that, given a set of n facilities and a set of m P 2 customers, the problem of maximizing the number of bichromatic reverse 1-nearest customer neighbors of a new facility is 3SUM hard [16], and give an algorithm for finding the set of all possible optimal placements of the new facility in O(m2) time. They also solve the problem of maximizing the number of bichromatic reverse 1-farthest customer neighbors of a new facility, to be placed in a given bounded region, in O(m log m) time. No results of the algorithms implementation are reported. Wong et al. [28] provide an algorithm for solving the problem of maximizing the number of bichromatic reverse k-nearest customer neighbors of a new facility, whose time complexity in the worse case is super-quadratic in the number m of customers. The performance of the algorithm deteriorates when k is large and in some cases may take hours to return the optimal region. Zhou et al. [35] present an algorithm for solving the same problem that works by partitioning the space into quadrants and searches only in those quadrants that potentially contain an

39

optimal region. Experimental results show that their algorithm is much more efficient than the previous algorithms. Yan et al. [30] present a grid-based approximation algorithm for a related bichromatic problem, its running time takes from seconds to a couple of minutes. It returns a prefixed given number of sub-optimal small regions, named influential locations, that collectively attract as many clients as possible. An important point to emphasize is that the above-mentioned existing work focuses on computing the optimal location for a new facility, which in many cases is very expensive and in general too restrictive, instead of finding a whole region of suitable sub-optimal locations as in our proposal. There also exist several papers dealing with a related problem which tries to maximize the area of the Voronoi region of a new facility [8,7,6,15]. Instead of trying to locate a new facility so that the weight of its reverse k-neighbors is maximized, they maximize the area of the Euclidean plane region having the new facility as one of the k-neighbors, in fact all of them consider k = 1 and the nearest case, except for Fort and Sellarés [15] who also handle k > 1 and the farthest criterion. Combinations of the problem of maximizing the number of bichromatic reverse 1-nearest customer neighbors of a new facility with classical distance based facility location problems are presented by Cabello et al. [4]. Given a set of n facilities and a set of m customers, m P n P 2, algorithms are provided for solving in O(m2+) time the problem of minimizing the maximum distance or maximizing the minimum distance to the customers associated to the solution that maximizes the number of bichromatic reverse 1-nearest customer neighbors of a new facility. There are no empirical studies of the algorithms, no implementation of these algorithms are known. 3.4. Theoretical results In this section we present an exact algorithm for computing k-influential maps. We also give some properties of monochromatic and bichromatic RegFarIk regions. 3.4.1. k-Influential maps computation Next an algorithm for computing MapNeark(F), the k-influential map in the monochromatic case for the nearest criterion, is described. The problem is solved by using an arrangement of disks with an approach similar to the ones used in Korn et al. [18] and Cabello et al. [4]. For each facility f 2 F, consider the closed disk Df centered at f and covering its k-nearest neighbors in F. Observe that f is a reverse k-nearest neighbor of a point q if and only if q 2 Df. The set of disks Df, f 2 F, can be computed using a monochromatic all k-nearest neighbors query algorithm. The computational complexity of the naive algorithm is O(n2), but more efficient methods exist [5,9]. Denote D the arrangement of the n disks Df, f 2 F. If a point q belongs to a cell d of D determined by the intersection of t disks, it means that the t facilities of F, the ones corresponding to these disks, have q between its k-nearest neighbors. Consequently, the influential value Neark(F, q) of all the points q of cell d is the same. The arrangement D can be computed in O(n2) time [2]. Finally, the k-influential value of each cell can be obtained by traversing the arrangement D through adjacent cells. To determine the k-influential value, id, of a cell d the weights of the t facilities determining d have to be added. Then, the k-influential value of a cell d0 , id0 , adjacent to d can be obtained from id by adding/substracting to id the weight of the facility whose disk enters/exits when going from d to d0 . In this way MapNeark(F), the k-influential value of all the O(n2) cells of D can be obtained in O(n2) time. For the remaining three cases, monochromatic using the farthest criterion and bichromatic using both nearest and farthest criteria, similar algorithms can be designed.

40

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

Fig. 3. Example of RegFarIk(F, I0) that does not intersect the boundary of the region of interest R.

In all cases, the complexity of the algorithm for computing k-influential maps is of quadratic order with respect the number n of facilities or the number m of customers. This has motivated us to explore an efficient alternative GPU parallel approach for solving this problem. 3.4.2. Properties of RegFarIk regions It can be proven that any not empty monochromatic RegFarI1(F, I0) or bichromatic RegFarI1(F, C, I0) region, generically denoted by RegFarI1(I0), intersects the boundary of the region of interest R. In fact, by slightly modifying the proof given by Cabello et al. [4]   for the bichromatic case, we can prove that RegFarI1 IM intersects 0 the boundary, where IM 0 denotes the maximal influential value in the corresponding influential map. Then for those I0 > IM 0 , RegFarIk(I0) will be empty, and for any I0 < IM we have that 0     M RegFarI1 I0  RegFarI1 ðI0 Þ, and consequently since RegFarI1 IM 0 intersects the boundary, RegFarI1(I0) also intersects it. However, the previous property is not true in general for k > 1 as can be seen in Fig. 3, where RegFarI47(F, 47) is represented for a set F of 48 facilities each of them of weight 1 and it does not intersect the boundary of R. An equivalent example can be provided for the bichromatic case. 4. Solving influential location problems on the GPU: the monochromatic case We approximately solve monochromatic location problems discretizing the region of interest R by superimposing an axis-parallel rectangular grid of size G = H  W on R. We denote Q the set of grid points, corresponding to the geometric center of the grid cells. The points of Q are not explicitly computed, they are indexed between 0 and G  1, both included, from bottom to top by row and within each row from left to right, so that given an index between 0 and G  1 we can determine its corresponding grid point from the grid dimensions and the region of interest vertices. First we compute the monochromatic k-influential value of each grid point of Q, and then we find among the points of Q those solving the desired location problem. 4.1. k-Influential map computation We want to compute the k-influential value of each grid point q of Q. The k-influential value of a point q is the sum of the weights of the facilities fi having q as one of the k-reverse neighbors. We use

the fact that fi is a k-nearest reverse neighbor of a point q if and only if fi is closer to q than to the kth-nearest facility of Fn{fi}. Analogously, the point q 2 Q is one of the k-farthest reverse neighbors of fi if the distance from fi to q is bigger than the distance to the kth-farthest neighbor of Fn{fi}. Thus to compute Ik(F, q) we start finding the distance to the kth-neighbor, and then it is properly compared with the distance from each q to fi. The k-neighbor distance is computed by storing the facilities of F in a uniform grid and using the algorithms mentioned in Section 2.3. The distance of each fi to the rest of facilities in Fn{fi} is found in parallel. Once the distances to the kth neighbor are known, we compute the k-influential value of all the grid points of Q and store them in a 1D-array g of size G which stores in the ith position the kinfluential value of the grid points with index i. The k-influential values are found in parallel using a kernel and we provide to each thread the: grid dimensions, region of interest vertices coordinates, facilities coordinates, facilities weights and array g where the influential values will be stored. There exist two different ways of parallelizing the k-influential value computation. We can either consider a thread per grid point of Q or a thread per facility in F. The strategy to be chosen depends on the discretization size and the number n of facilities of F, as can be seen from the complexity analysis and the experimental results provided in Section 7. Next both algorithms are presented. 4.1.1. A thread per grid point All points of Q are processed in parallel by considering a thread per point qj. The thread with global index j processes the point qj and finds its coordinates from index j, the grid size and the region of interest vertices. Each thread explores all the facilities, and decides whether the weight of each facility fi has to be added to Ik(F, qj) by comparing the distance from fi to qj, d(fi, qj), with the k-neighbor distance of fi, dk(fi). A facility weight is added to the NearIk(q, F) when d(fi, qj) < dk(fi), and to FarIk(q, F) when d(fi, qj) > dk(fi). No facilities can be discarded a priori because the distance from qj to facility fi has to be compared with its corresponding dk(fi) which depends on fi and not on qj. In fact, the influential value is calculated accumulating the facilities weights in an auxiliary variable whose final value is stored in global memory in the corresponding j position of the influential values array. Notice that each thread has to explore all the facilities, for this reason the facilities are uploaded to shared memory in blocks of size B, the block size used when executing this kernel. In fact the thread occupying position i in the block reads from global memory and loads to local memory the facilities with indices i + kB for k = 0    n/B  1 in n/B different steps. Then it reads the n facilities from local memory. In this way each facility is read from global memory only once per block and from local memory B times per block. To obtain the desired results we have to add several synchronization points. The threads in a block have to wait until all the current B facilities have been loaded, and then until all the threads in the block have processed the currently loaded facilities before loading the next B ones. Obtaining the influential values using a thread per grid point of Q takes O(Neig(k, n) + Gn), where O(Neig(k, n)) represents the total work of the algorithm used to find the k-neighbor problem for the monochromatic case. The number of memory accesses after  obtaining the k-neighbor distances is O wgG þ rgGn=B þ wsGðn=Bþ1Þ þrsGn Þ, where w and r denote written and read values, the superindices indicate whether they are done to global or shared memory and the subindices refer to the number of accesses done. Notice that to obtain the total number of accesses to memory we should consider the ones done to build the uniform grid and to obtain the kth-neighbor for each facility.

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

4.1.2. A thread per facility In this case we have as many threads as facilities, and the array where the influential values are stored has to be initialized to zeros. To compute NearIk(qj, F), the weight wi has to be added to the influential values of all those points qj whose distance to fi is smaller than dk(fi). Similarly, to compute FarIk(qj, F) the weight wi has to be added to the influential values of all those points qj whose distance to fi is greater than dk(fi). We solve the nearest problem with a multiple circular in range query, and when the farthest case is considered a multiple circular out range query is handled. Each thread solves a circular range searching query with center on the facility fi, and radius equal to dk(fi). During the query we determine which points are contained in the range and add the weight wi to their k-influential value. To solve the multiple range searching query we use the algorithm of Fort and Sellarès presented in Section 2.3. However we avoid storing explicitly the points of set Q in a uniform grid, taking into account that they are obtained from an axis-parallel rectangular grid. Thus we work with this inherent uniform grid without building it. While computing k-influential values, handling all the facilities in parallel, we have to avoid to add simultaneously two weights to the same k-influential value. It is done by using atomic functions, which produces a slow down in the algorithm time performance, but guarantees correct results. The total work of the process is O(Neig(k, n) + R) where R is the total sum of the number of points contained in the ranges, which depends on n, G and the kth neighbor distances. In the worst case R = O(Gn). In this case the number of accesses to memory to be added are O rg2n þ wgdR , where d refers to the execution time required by each single atomic function. 4.2. Solving RegIk(F, I0) problem Once the k-influential values of the points of Q are computed and stored in the array g, to solve RegIk(F, I0) problem we have to find those points of Q where the influential value is at least I0. We count the number of grid points where the influential value is at least I0, allocate an integer array r of this size and finally store in it the indices of the points conforming the region RegIk(F, I0). Obtaining the optimal points takes a total work of O(G + P), where P denotes the number of points conforming the region. The total number of   accesses to memory of this part of the algorithm is O wgP þ r gn . As a particular case if we were interested in finding the points where the optimal influential value is obtained according to the considered discretization, we can compute the maximum influential value stored in g by using existing algorithms to find the maximum of a set of integer values that run under CUDA architecture [24]. After obtaining the maximal value, we proceed as before finding the points whose k-influential value is at least the obtained one. 4.3. Solving combined problems We have presented four different combined facility location problems, MinMaxNeark(F, I0) and MinSumNeark(F, I0) deal with nearest neighbors, and MaxMinFark(F, I0) and MaxSumFark(F, I0) with farthest neighbors. Notice that the optimal point is always found in the optimal region RegIk(F, I0). The process to solve these four facility location problems is almost the same, thus we only explain how the ones related to the nearest neighbors are solved. In order to solve the MinMaxNeark(F, I0) and the MinSumNeark(F, I0) problems for the monochromatic case, we have to store not only the k-influential value, but also the maximal distance or the distances sum to the facilities having each grid point as one of the reverse k-neighbors. We have to distinguish between two

41

different cases depending on wether the optimal region RegNearIk(F, I0) is known or not. 4.3.1. The region RegNearIk(F, I0) is not known In this case we can solve the MinMaxNeark(F, I0) or MinSumNeark(F, I0) problem by computing the maximum or the sum of the k-reverse neighbors during the process used to find the k-influential values. If a thread per spatial point is used, the maximum or sum of the distances of the facilities having qj as one of the k-reverse neighbors is computed by using an extra register. Thus, during this process we not only maintain the sum of all those weights wi of the facilities fi that have qj as a k-reverse neighbor but also the maximum or the sum of d(qj, fi). When all the facilities have been considered, we store in the j position of a float array the obtained maximum distance or distances sum. The algorithm time complexity remains the same  and the number of accesses to memory is incremented in O wgG , since we need to store the obtained distance value at the corresponding j position. In the case that a thread per facility fi is used, we proceed as in the case of finding the k-influential values. We solve a range searching problem for each facility fi and consider all those points qj contained in the range associated to fi. Now, instead of only accumulating the weights, we also sum the distances d(qj, fi) or maintain their maximum at the j position of the array where they are stored. As before, to obtain the correct results we have to use atomic functions to sum values or maintaining the maximum. The algorithm complexity does not change, because we only need an extra atomic operation (addition or maximum obtention) for each point contained in a range. In terms of memory reads and   writes it represents an extra O wgdR , where R is the total sum of the number of points contained in the ranges, which in the worst case becomes R = O(Gn). In order to find the solution to the problem, after computing the k-influential values together with the maximal distance or distances sum for all the points of Q, we have to find the point of the obtained region with minimal associated value related to the distances. Thus we do not need to find all the points defining the region but select only those points that, having k-influential value at least I0, also have associated the minimal maximum distance or distances sum. 4.3.2. The region RegNearIk(F, I0) is known In the case that we already know the region RegIk(F, I0), we do not need to re-start the algorithm but take advantage of this information and consider only the points, which generally are a small number, conforming the maximal region. To solve the problem we have to compute the maximal distance or the sum of the distances of all those facilities having one of the points in the maximal region as a k-reverse neighbor. Again the previous two strategies can be used. We could use a thread per grid point or a thread per facility. In the former case, when we use a thread per grid point, no changes need to be done. However we must take into account that we have in an integer array r the indices of the points contained in RegIk(F, I0). Thus thread with index i will not represent qi but qr(i), where r(i) stores the index of a point contained in the region RegNearIk(F, I0). The algorithm has a total work O(nP), and the number of   memory accesses is O rgPþPn=B þ rsPn þ wgP þ wsPn=B , where again P is the number of points in RegIk(F, I0). In the latter case, when a thread per facility is used, the set of points defining the maximal region is not placed in a grid, we have a set of indices but we do not know in advance any information of their spatial distribution. In this case instead of building a uniform grid containing them we take advantage of the fact that the

42

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

number of points conforming the optimal region is small, and all the points are checked for each facility. Thus the points are transferred to shared memory as it is done in the former strategy dealing with facilities. The total work of this strategy becomes again O(nP) and the total number of memory accesses is   O rgnP=Bþn þ r snP þ wgdR0 þ wgnP=B , where R0 denotes the number of points contained in one of the considered ranges, which in the worst case is R0 = O(nP).

5. Solving reverse influential location problems on the GPU: the bichromatic case To solve the bichromatic problems, we need to introduce small changes in all the presented algorithms dealing with monochromatic problems. Since the main algorithm is the one that computes the influential map, the influential value of all grid points in Q, we explain the changes it needs. The rest of the algorithms have to be similarly adapted. Now Ik(F, C, qj) is the sum of the weights of those customers ci 2 C being reverse k-neighbors of point qj 2 Q. Thus, we have to add wi whenever qj is one of the k-neighbors of ci taking into account the set of facilities F [ {qj}. Consequently we have to compare the distance from ci and its k-neighbor of F with the distance between ci and qj, and if necessary to sum or maintain their maximum. We start by finding for each customer ci the distance to its kth-facility, dk(ci), by considering all the customers of C in parallel. Then, to find the k-influential values, we consider the grid points Q, the set C of customers and the distance to the kth-facility. For each grid point qj, its distance to ci is compared with the distance dk(ci) according to the used proximity criterion. As before, the kth-influential value can also be computed by handling all grid points or all customers in parallel. Again the option to be used depends on the discretization size and the number m of customers. Computing the k-influential value for each customer has a total work of O(Neig(k, n, m) + Gm), where m is the number of customers and O(Neig(k, n, m)) represents the total work of the used algorithm to find the k-neighbor client of each customer. The number of associated memory accesses in the last part of the algorithm, where the k-influential values are computed from the kth neigh  bor distances, is O wgG þ wsGð2m=Bþ1Þ þ r s2Gm þ rgGm=B . When a thread per query is used, the work of this final part is O(Neig(k, n, m) + R), where R is the number of points contained in all the ranges which depends on m, G and the kth neighbor distances, and in the worst case it becomes R = O(Gm). In this case the number of accesses to   memory to be added are O rg2m þ wgdR .

6. Results visualization process We have designed an interface that allows us to visualize: (a) influential maps; (b) solution regions of reverse influential location problems; and (c) solution points of combined reverse influential location with single facility location problems. This visualization improves the understanding of the problems and facilitates decision makers, with an iterative what-if-analysis process in which some parameters could be changed (such as the minimum desired influential value, the visualized area or the discretization size), to obtain an approximate optimal location for a new facility. Each iteration is intended to provide additional detail to ensure a better solution. The short response times for computing and visualizing the results facilitate the interaction with the user.

6.1. Influential maps visualization To visualize the influential map, we use a color gradation according to the k-influential values, the darker a point is painted, the smaller is its k-influential value. In order to properly use all the gradation, the maximal k-influential value, IM 0 , is found and the lighter color is associated to this k-influential value. Fig. 4 presents the k-near influential values of a set of 100 red1 facilities for a grid of size 500  500 for: (a) k = 1, (b) k = 5 and (c) k = 25. The maximal influential value obtained in each case which are represented by the lighter color corresponds to a k-influential values of 27, 177 and 260, respectively. We can see that each circle contains exactly k facilities, and hence, the grid points it contains have the circle center as one of the k-reverse neighbors. Fig. 5 shows the k-farthest influential values of again a set of n = 100 facilities using a grid of size 500  500, for k = 1, 5 and 25. In this case we are interested in the grid points that are not contained in the circles. Notice that in the first image where k = 1 there exist only four different facilities tangent to the circles, thus all the points in F have as farthest facility one of these four facilities. For k = 5 the exterior of each circle contains five facilities, thus the circles are smaller than in the former case but bigger than the ones in Fig. 5c for k = 25. In this case the lighter color corresponds to a k-farthest influential value of 224, 280 and 410 for k = 1, k = 5 and k = 25, respectively. Fig. 6 presents a bichromatic influential map for n = 100 red facilities and m = 100 blue clients, k = 5 and the grid of size 500  500. Again a color gradation is used. Fig. 6a is associated to a nearest problem and Fig. 6b to a farthest one. In these cases the maximal influential values correspond to 143 and 389, respectively. In Fig. 7 we show the 1-influential values of 100 facilities by using a two different discretization sizes 100  100, and 500  500, in image (a) and (b), respectively. Comparing both images we can see how the boundaries of the obtained circles are much more pixelized for the smaller grid. In the case that we are interested in changing the discretization size all the computations have to be rebuild. Fig. 8 shows the zoom option on a 5-influential map of a set of 100 sites and 100 clients with a grid of size 1000  1000. The zoom is done in the left top corner of the domain in order to easily identify the area. This is one of the options of our application, we can zoom on any desired region until the maximal precision according to the used grid is obtained. Zooming on a desired area has an immediate response.

6.2. Influential location regions visualization To visualize the influential location region RegIk(F, I0) or RegIk(F, C, I0), we paint in a different color gradation the points defining the desired region. Accordingly, we paint in a green gradation the points whose influential value is bigger than the minimum desired one, I0, and in purple the ones with smaller k-influential values. Notice that, once the influential values have been computed obtaining the visualization of the influential location regions for a new minimal desired influential value is immediate and needs no recomputations. The green areas of Fig. 9 correspond to the influential locations associated to minimum 1-influential values 15 and 30 when a grid of size 500  500 is considered, in this case the maximal influential

1 For interpretation of color in Figs. 4–6 and 9–18, the reader is referred to the web version of this article.

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

Fig. 4. Nearest influential maps for with n = 100 and a grid of size 500  500 for: (a) k = 1; (b) k = 5; and (c) k = 25.

Fig. 5. Farthest influential maps, using a grid of size 500  500, n = 100 and: (a) k = 1; (b) k = 5; and (c) k = 25.

Fig. 6. 5-Influential maps with n = m = 100 and a grid of size 500  500 for the (a) nearest and (b) farthest case.

Fig. 7. 1-near influential maps with n = 100 and a grid of size: (a) 100  100 and (b) 500  500.

43

44

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

Fig. 8. A grid of size 1000  1000, n = m = 100 and 5-near influential map with the left top corner zoomed.

Fig. 9. (a) RegNear1(F, 15) and (b) RegNear1(F, 30), with n = 100 and a grid of size 500  500.

Fig. 10. Bichromatic case (a) RegNearI1(F, C, 15) and (b) RegFarI15(F, C, 200), in both cases n = m = 100, and the grid is of size 500  500.

Fig. 11. Grid of size 500  500 and n = 100, (a) MinMaxNear1(F, 15) and (b) MinMaxFar1(F, 300).

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

45

Fig. 12. Grid of size 500  500, n = m = 100. (a) MinSumNear15(F, C, 200) and (b) MinMaxNear15(F, C, 200).

Fig. 13. Running times for the monochromatic k-nearest influential map computation.

value obtained with this discretization is 45. Notice how RegNear1(F, 30)  RegNear1(F, 15). In Fig. 10 two different images of bichromatic influential location regions are provided. In both images n = m = 100 and a grid of size 500  500 is considered. The green parts are the ones conforming the influential location RegNearI1(F, C, 15) in image (a), where the maximal influential value is 55. In Fig. 10b, the green points define RegFarI15(F, C, 200), in this case the maximal influential value is of 449.

6.3. Solution of combined problems visualization To visualize the solution of a combined problem we keep on visualizing all the k-influential map, with the points not contained in the desired influential location with a dark purple gradation, those contained in the region with a green color gradation and the solution of the desired optimization problem in a light yellow color, in the images the solution is remarked with a square.

46

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

Fig. 14. Running times for the monochromatic k-farthest influential map computation.

Table 1 k-Influential values computational times, in (s), once the k-neighbor distance are already known for k = 10, using one thread per facility (Fac) and one thread per grid point (GP). n Fac

100 500 1000 5000 10,000 25,000 30,000

Near

Far

Fac

GP

Fac

GP

Fac

GP

Fac

GP

0.0050 0.0014 0.0008 0.0006 0.0006 0.0003 0.0003

0.0001 0.0002 0.0002 0.0010 0.0020 0.0048 0.0058

0.4420 0.1088 0.0604 0.0190 0.0167 0.0168 0.0167

0.0031 0.0101 0.0182 0.0836 0.1663 0.4144 0.4971

0.0099 0.0077 0.0070 0.0084 0.0127 0.0124 0.0126

0.0001 0.0002 0.0003 0.0010 0.0020 0.0048 0.0053

0.9309 0.6923 0.6285 0.7417 1.1396 1.4525 1.4296

0.0034 0.0104 0.0186 0.0836 0.1663 0.4151 0.4960

100  100

1000  1000

Given a set F of n = 100 facilities, Fig. 11a shows in green the RegNearI1(F, 15) and the yellow point contained in the square is the point solving the MinMaxNear1(F, 15) problem, in this case the optimal is achieved in a single point by considering a grid of size 500  500. Similarly, Fig. 11b presents the MinMaxFar1(F, 300) solution, which is defined by the yellow point on the left top corner of the image. Again the green region is RegFarI1(F, 300), in this case the maximal influential value is of 370. Finally in Fig. 12 we present a bichromatic example of the combined problems. Fig. 12a presents the solution to the MinSumNear15(F, C, 200) problem, and Fig. 12b solves the MinMax Near15(F, C, 200). In both cases we can see a single optimal point where the minimum is obtained. Notice that both optimal points

100  100

1000  1000

are contained in the same connected component of RegNearI15(F, C, 200), but in different positions.

7. Performance evaluation In this section we provide experimental results of the presented algorithms. We have considered several sets of facilities in the monochromatic case, and facilities and customers in the bichromatic case, as well as grids of different sizes in order to show the scalability and behavior of our algorithms in realistic settings. The number of facilities in the presented experimental results varies from n = 100 to 60,000 for both the monochromatic and

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

47

Fig. 15. Running times for the bichromatic k-nearest influential map computation for 10,000 facilities.

bichromatic cases, the number of customers in the bichromatic case also goes from m = 100 to 60,000. The k value varies from 1 to 100, and the discretization size goes from G = 100  100 to 1000  1000 and k value from 1 to 500. We have assigned random integer weights between 1 and 10 to the facilities, which are randomly placed in a region of interest of size [0, 0], [100, 100]. Accordingly, for instances, when a grid of size 500  500 is used each pixel represents an area of 0.04 units2, and for a grid of size 1000  1000 it becomes 0.01 units2. Our application allows, given a set of facilities to change the k value and the grid size, to recompute the obtained results. Notice that the considered values, in terms of set or discretization sizes and k values, are large enough to show the scalability of the presented algorithms when solving the considered problems in realistic environments. The provided times are obtained with the implementation of our algorithms using CUDA C, and the executions are done using a i5-200 3.30 GHz CPU with 8 GB of memory and a Nvidia GTX 580. The provided times are obtained as the median of 10 executions. When the GPU is used the time needed to transfer information from CPU to GPU is an important issue, because transferring information is expensive in terms of processing time. However, in our case we only transfer the facilities coordinates once in the monochromatic case and the facilities and customers coordinates in the bichromatic case. Notice that we do not transfer the grid points to the GPU, we only need to know the region of interest vertices, and the discretization grid dimensions. All the presented problems

can be solved by using the GPU without needing GPU–CPU communication. Thus, the transferring times involved in our algorithms are does needed to transfer sets of 100, 1000, 10,000, 30,000 or 60,000 points, 2d coordinates, from CPU to the GPU which are of 0.05 (ms), 0.05 (ms), 0.09 (ms), 0.15 (ms) or 0.22 (ms) respectively. Before starting presenting our running times we want to mention that there exist two related works solving bichromatic problems that also provide experimental results. Wong et al. [28], provide two different algorithms to find the region with maximal k-nearest influential value. They test their algorithms for small sets of clients and facilities, from 50 to 250 and k = 1 which take between 0.001 and 0.1 (s). Then they consider bigger sets, with 60,000 facilities and k = 1, the running time becomes 800 (s), meanwhile and for k = 5 and 20,000 facilities it takes 95,000 (s). On the other hand, Yan et al. [30] provide another method to find influential locations that uses adaptive grids, their running times for k = 1 and the nearest case are about 6 (s) for a set of 62,000 clients and 20,000 facilities, mean while, for k = 10 their running time increments to 800 (s). The times provided in the following subsections show how our algorithms have much better running times than the ones they provide. 7.1. Monochromatic case We provide a wide analysis of results obtained for the monochromatic case showing the differences between the proposed methods and strategies. The analysis shows how depending on

48

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

Fig. 16. Running times for the bichromatic k-farthest influential map computation for 10,000 facilities.

the set sizes and the grid dimensions it is more convenient to choose the strategy that uses a thread per grid point or per facility. 7.1.1. Influential maps computation The times needed to find the k-influential maps for several k, n and grid size values are presented in Figs. 13 and 14. Fig. 13 shows the times needed to find the k-nearest influential maps when a thread per facility and a thread per grid point are considered. Fig. 14 presents the times related to the k-farthest influential maps, again when a thread per facility and a thread per grid point are used. The times provided in these figures include the k-neighbor distance obtention. Accordingly, if a new strategy to find the k-neighbor appears, the presented running times would improve. Finding the k-influential values is really fast, as can be seen in Table 1 which shows the time needed for each method to find the kinfluential values once the k-neighbor distances are known, for k = 10. Analyzing Fig. 13a we can see that when using a thread per facility and fixed a grid resolution, the time needed to compute k-influential values decreases when the number of facilities increases from n = 10 to n = 5000, then the time starts increasing again. This happens because the more facilities we have the smaller the k-neighbor distance is and, consequently, the smaller is the number of grid points that have to be explored to solve the problem. Then, for n = 10.000 or greater, the time starts incrementing again due to the k-nearest see a more regular behavior in the one thread for grid point strategy. The time increases according to k, for small values of k the time slightly changes, but for k P 50, again

the k-nearest neighbor computation affects the running times. Concerning the grid discretization the bigger the grid, the bigger the running time. Exactly the same happens when the number of facilities increases. Comparing both methods, the first method provides better results in almost all the cases for sets of at least 5000 facilities. In Table 1 we provide the time needed to compute the kinfluential value to properly compare both strategies. In general, when the number of facilities is big, the strategy that uses one thread per grid point is not as efficient as the other because each grid point has to explore all the facilities. However, if the number of facilities is much bigger than the discretization size, and specially for big values of k, using a thread for facility has the problem that too many atomic operations are done in the same grid point, producing a big number of serialized operations. This problem appears for n = 30,000 using a grid of size 50  50 and k = 50 or 100; or a grid of a size 100  100 and k = 100. Related to this figure we also want to mention that the time needed for the set of n = 100 facilities varies from 0.3 (ms) to 2 (ms), this maximum is achieved for k = 50 and a grid of size 1000  1000. Next, in Fig. 14 we can see the running times for the farthest case for both strategies. Fig. 14a provides times obtained when one thread per facility is used, and Fig. 14b represents the one thread per grid point strategy. In this case the high of the columns for a grid of size 1000  1000 and sets of 25,000 and 30,000 facilities have been specified in numbers and they have been truncated to preserve the scale and allow comparing both strategies for the smaller grids, k-values and sets of sites (something similar has been done for the bichromatic case).

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

49

Fig. 17. Running times for the bichromatic k-nearest influential map computation for 10,000 customers.

For the farthest case when a thread per grid point is used we obtain faster results than for the one thread per grid point strategy. The draw back of the one thread per facility strategy are the atomic increments, and in general the k-farthest influential location values are bigger than the k-nearest ones, this means that the atomic increments are concentrated in the same gird point making that this strategy gets, in comparison, worst results with the farthest criterion. The strategy to be used in the farthest case should be one the thread per grid point one. Comparing the times provided in Figs. 13 and 14 we can see that, when a thread per facility is used the farthest case is always slower than in the nearest case due to the fact that the number of cells to explore when solving the farthest case tends to be bigger than when solving the nearest one and moreover the algorithm is not so simple. However, this is not true for the one thread per grid point, in both cases computing the farthest k-influential value takes the same time (see Table 1), and depending on the grid and k values we can see that finding the k-farthest neighbor distance is faster than finding the k-nearest neighbor distance. The presented algorithms can also handle bigger sets of facilities, for instances, for k = 500, n = 30.000 and the four considered grid sizes, the obtained times in the nearest case using thread per facility are 12.28, 12.39, 12.56 and 13.11 (s), using one thread per grid point they become 12.40, 12.41, 12.49 and 12.85 (s). For the farthest case they become, respectively, 4.92, 4.93, 6.06 and 9.52 (s) and 4.92, 4.93, 5.05 and 5.42 (s). When 60.000 facilities are used the corresponding times for a grid of size 1000  1000 and k values of 1, 100, 500 become, for

the nearest case, 1.52, 3.55 and 23.12 (s) using a thread per facility; and 1.55, 4.39 and 23.40 (s) using a thread per grid point. For the farthest case for a thread per facility they are 6.62, 8.31 and 22.52 (s), using a thread per grid point 1.29, 2.35 and 15.38 (s). 7.1.2. Influential location region computation After computing the influential maps we determine the maximum obtained influential value in order to use a proper color gradation to represent the popularity map and either to facilitate or guide the minimum required influential value used for the specific considered problem. In order to find this maximum influential value we have implemented a simple reduction type algorithm, that only depends on the grid size, which uses a single block. Each thread finds the maximum of several positions, and to end the process the first thread finds the global maximum once all the other threads have finished. Finding the maximum of a set of integer values with this algorithm takes 0.11 (ms), 0.15 (ms), 0.36 (ms), and 1.09 (ms), for a grid of size 50  50, 100  100, 500  500 and 1000  1000 respectively. However, a more elaborated reduction type algorithm can be used if better times are needed. Visualizing the influential location region is really fast, takes less than a millisecond. 7.1.3. Combined problems solution Concerning the computational times, we have provided two different strategies to solve a combined problem. One strategy assumes that the RegIk(F, I0) is known, then the distances associated to its points are computed, and the optimal point is found. The

50

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

Fig. 18. Running times for the bichromatic k-farthest influential map computation for 10,000 customers.

other strategy computes the distances during the k-influential value computation and then the combined problem is directly solved. When using the nearest criterion we can be interested in the MinMaxNeark(F, I0) or the MinSumNeark(F, I0) problem, with the farthest criterion the problems become MaxMinFark(F, I0) and MaxSumFark(F, I0). Next we provide times related to MinMaxNeark (F, I0) problem as examples of these running times. Solving a MinMaxNear1(F, 15) problem for a set of n = 1000 facilities using a grid of size 1000  1000 takes 0.03 (s) using either one thread per facility or one thread per query, meanwhile only finding the k-influential map takes 0.015 (s) using one thread per facility and 0.012 (s) when a thread per grid point is used. Solving MinMaxNear2(F, 15) takes, 0.032 (s) using one thread per grid point when finding the 2-influential map takes 0.019 (s), meanwhile when using a thread per grid point, the times become 0.033 (s) and 0.018 (s). Notice that the one thread per facility is clearly less affected when solving directly the combined problem and this is due to the fact that its drawback are the atomic functions serialization, and since both, the increment of the influential value and the maximum or sum maintenance are done in the same position, the serialization is produced only once. Meanwhile when using one thread per grid point, the experimental times correspond to the time needed to find the k-neighbor distance plus about twice the time needed to found the k-influential value, because the number of operations done per kernel is multiplied by two. On the other hand when the points conforming RegIk(F, I0) are known, the distances have to be find only for these points. For

M big values of I0  IM 0 , where I0 is the maximal influential value, RegIk(F, I0) tends to be a small part of R, and using the region is really useful. In general this strategy makes sense for I0 values determining small RegIk(F, I0). With the same example provided in the previous paragraph once RegIk(F, 15) has been determined, solving the MiMaxnNear1(F, 15) problem takes 0.0072 (s) using one thread per facility and 0.0054 (s) using one thread per grid point. Meanwhile solving MinMaxNear2(F, 15) takes, respectively 0.075 (s) and 0.011 (s). This important increment of time happens because for k = 1 the maximum influential value is 52, and for k = 2 it becomes 65, and in both cases we consider I0 = 15, it makes that for k = 2 the region Reg2(F, 15) contains 450.234 points meanwhile Reg1(F, 15) contains only 136.543 points.

7.2. Bichromatic case In this section we provide some times obtained for the bichromatic case when computing the k-influential maps, obtaining k-influential location regions only depends on the discretization size, and solving combined problems on the are of the region when this is known or on the k-influential map computation otherwise. Figs. 15 and 16 present the times needed to find the influential maps for several k, m and grid sizes for n = 10,000 facilities, when the algorithm uses one thread per grid point or per customer, respectively. These times again include the k-neighbor distance obtention. For each customer we find the k-neighbor in the set of 10,000 facilities, and the k-influential value is computed. The

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

analysis would be the same as the one provided for the monochromatic case, however in this case the time needed to find the k-nearest neighbor is that of finding for each of the m customers the k-nearest facility in the set of 10,000 facilities. Again, the one thread per facility strategy gives better results for big grid sizes considering the nearest case, but when the farthest criterion is considered the one thread per grid point always provides better running times. Next, in Figs. 17 and 18 a set of m = 10,000 customers is considered and we provide the times to solve several k-influential maps for several k values, grid sizes and sets from 100 to 30,000 facilities, for the nearest and farthest case. Notice, that in this case the time changes due to: (a) the time needed to find the kth distance, we always solve 10,000 queries but considering sets with a different number of facilities and (b) the kth distance which affects only the k-influential value computation when a thread per facility is used. These Figures also support the already provided conclusions. Again we can handle bigger sets of both customers and facilities, for a set of 60,000 customers and 10,000 facilities for a grid size of 1000  1000 and k = 1, 100, 500, the obtained times using the nearest criterion and with a thread per facility are: 0.049, 0.062 and 4.37 (s), when a thread per grid point is used they become 0.24, 0.76 and 4.43 (s). For the farthest case, they are, respectively 1.09, 1.44 and 4.23 (s) and 0.22, 0.43 and 2.99 (s). When 10,000 customers and 60,000 facilities are used the obtained times, in the order provided in the previous paragraph are: 0.08, 2.25 and 14.27 (s), 1.06, 2.49 and 11.63 (s), for the nearest case; and 6.35, 8.74 and 16.15 (s) and 1.04, 1.55 and 6.31 (s) for the farthest one. 8. Conclusions and farther comments In the paper we presented practical and scalable algorithms to compute the region RegIk, which contains the points such that the location of a new facility in it guarantees a minimum k-influential value, taking into account the already existent facilities in the monochromatic case or the existent facilities and the potential customers in the bichromatic case. We have proposed strategies to locate desirable and undesirable or obnoxious facilities by computing the influential value associated to the k-nearest or k-farthest facilities. The provided algorithms are based on a uniform discretization of R by using a grid. Two different parallel computing approaches, by using one thread per grid point or per facility are presented together with the advantages, disadvantages and requirements in terms of GPU programming of each strategy. The used discretization makes that we may obtain sub-optimal maximal regions if the discretization is not fine enough to have a grid point in the optimal one. However, as the experimental results show, we can consider big enough grids to obtain quite accurate results. We can guarantee that the optimal region is detected if it contains an axis parallel rectangular area equal to the one defining a grid cell. We have also provided algorithms to solve combined facility location problems such as MinMaxNeark, MinSumNeark, MaxMinFark and MaxSumFark problems. Images obtained with the presented way to visualize the obtained solutions, such as RegIk regions or the points solving combined facility location problems, improving the understanding of the problem and revealing complicated structures that would be hard to capture otherwise, are provided. Finally experimental results on realistic settings corroborating the complexity analysis presented during the algorithms description are provided. Acknowledgments We thank the reviewers for their suggestions and comments. Authors were partially supported by the Spanish MCI Grant TIN2010-20590-C02-02.

51

References [1] V. Akman, W.R. Franklin, M. Kankanhalli, C. Narayanaswami, Geometric computing and the uniform grid data technique, Computer Aided Design 21 (7) (1989) 410–420. [2] N.M. Amato, M.T. Goodrich, E.A. Ramos, Computing the arrangement of curve segments: divide-and-conquer algorithms via sampling, in: Proc. 11th ACMSIAM Sympos. Discrete Algorithms, 2000, pp. 705–706. [3] R. Benetis, S. Jensen, G. Kariauskas, S. Saltenis, Nearest and reverse nearest neighbor queries for moving objects, The VLDB Journal 15 (2006) 229– 249. [4] S. Cabello, J.M. Díaz-Báñez, S. Langerman, C. Seara, I. Ventura, Facility location problems in the plane based on reverse nearest neighbor queries, European Journal of Operational Research 202 (1) (2009) 99–106. [5] P.B. Callahan, S.R. Kosaraju, A decomposition of multidimensional point sets with applications to k-nearest-neighbors and n-body potential fields, Journal of the ACM 42 (1) (1995) 67–90. [6] O. Cheong, A. Efrat, S. Har-Peled, Finding a guard that sees most and a shop that sells most, Discrete and Computational Geometry 37 (4) (2007) 545– 563. [7] F. Dehne, R. Klein, R. Seidel, Maximizing a Voronoi region: the convex case, International Journal of Computational Geometry and Applications 15 (5) (2005) 463–476. [8] M. Denny, Solving geometric optimization problems using graphics hardware, Computer Graphics Forum 22 (3) (2003) 441–452. [9] M.T. Dickerson, D. Eppstein, Algorithms for proximity problems in higher dimensions, Computational Geometry: Theory and Applications 5 (5) (1996) 277–291. [10] T. Drezner, Optimal continuous location of a retail facility, facility attractiveness, and market share: an interactive model, Journal of Retailing 70 (1) (1994) 49–64. [11] T. Drezner, Z. Drezner, Validating the gravity-based competitive location model using inferred attractiveness, Annals OR 111 (1-4) (2002) 227– 237. [12] T. Drezner, Competitive facility location, Encyclopedia of Optimization (2009) 396–401. [13] Z. Drezner, Hamacher, Facility Location – Applications and Theory, Springer Verlag, 2002. [14] M. Fort, J.A. Sellarès, Multiple Range Searching on the GPU, Actas XIV Encuentros de Geometra Computacional (ISSN: 2014-2323), Alcal de Henares, 2011, pp. 75–78. [15] M. Fort J.A. Sellarès, GPU-Based Influence Regions Optimization. ICCSA, LNCS 7333, 2012, pp. 253–266. [16] A. Gajentaan, M. Overmars, On a class of O(n2) problems in computational geometry, Computational Geometry Theory and Application 5 (1995) 165– 185. [17] H.V. Jagadish, B.C. Ooi, K.L. Tan, C. Yu, R. Zhang, iDistance: an adaptive B+-tree based indexing method for nearest neighbor search, ACM Transactions on Database System 30 (2) (2005) 364–397. [18] F. Korn, S. Muthukrishnan, D. Srivastava, Reverse nearest neighbor aggregates over data streams, in: VLDB ’02: Proceedings of the 28th International Conference on Very Large Data Bases, VLDB Endowment, 2002, pp. 814– 825. [19] Y. Kumar, R. Janardan, P. Gupta, Efficient algorithms for reverse proximity query problems, in: Proc. 16th ACM-GIS, vol. 39, 2008, pp. 1–10. [20] M.D. Lieberman, J. Sankaranarayanan, H. Samet, A fast similarity join algorithm using graphics processing units, International Conference on Data Engineering (2008) 1111–1120. [21] P. Leite, J.M. Teixeira, T. Farias, V. Teichrieb, J. Kelner, Massively parallel nearest neighbor queries for dynamic point clouds on the GPU, in: 21st International Symposium on Computer Architecture and High Performance Computing, 2009, pp. 19–25. [22] L. Mussi, F. Daolio, S. Cagnoni, Evaluation of parallel particle swarm optimization algorithms within the CUDA architecture, Information Sciences 181 (20) (2011) 4642–4657. [23] S. Nickel, J. Puerto, Location Theory – A Unified Approach, Springer, 2009. [24] M. Harris, Optimizing Parallel Reduction in Cuda, White Paper, NVIDIA Developer Technology, 2007 [25] J.D. Owens, D. Luebke, N. Govindaraju, M. Harris, J. Kröger, Aaron E. Lefohn, Timothy J. Purcell, A survey of general-purpose computation on graphics hardware, Computer Graphics Forum 26 (1) (2007) 80–113. [27] M. Smid, Closest-point problems in computational geometry, in: J.-R. Sack, J. Urrutia (Eds.), Handbook of Computational Geometry, Elsevier, 2000, pp. 877– 935. [28] R.C. Wong, M.T. Özsu, P.S. Yu, A.W. Fu, L. Liu, Efficient method for maximizing bichromatic reverse nearest neighbor, Proceedings of the VLDB Endow 2 (2009) 1126–1137. [29] W. Wu, F. Yang, C.Y. Chan, K. Tan, FINCH: evaluating reverse k-NearestNeighbor queries on location data, Proceedings of the VLDB 1 (1) (2008) 1056– 1067. [30] D. Yan, R.C. Wong, W. Ng, Efficient methods for finding influential locations with adaptive grids, ACM-CIKM (2011) 1475–1484. [31] B. Yao, F. Li, P. Kumar, Reverse furthest neighbors in spatial databases, in: Proc. Int. Conference on Data Engineering, 2009, pp. 664–675.

52

M. Fort, J.A. Sellarès / Knowledge-Based Systems 47 (2013) 35–52

[32] A. Yaseen, Y. Li, Accelerating knowledge-based energy evaluation in protein structure modeling with graphics processing units, Journal of Parallel and Distributed Computing. 72 (2) (2012) 297–307. [33] Y. Zhang, F. Mueller, X. Cui, T. Potok, Data-intensive document clustering on graphics processing unit (GPU) clusters, Journal of Parallel and Distributed Computing 71 (2) (2011) 211–224.

[34] J. Zhao, Q. Liu, W. Wang, Z. Wei, P. Shi, A parallel immune algorithm for traveling salesman problem and its application on cold rolling scheduling, Information Sciences 181 (7) (2011) 1212–1223. [35] Z. Zhou, W. Wu, X. Li, M.L. Lee, W. Hsu, MaxFirst for MaxBRkNN, in: IEEE 27th International Conference on Data Engineering, 2011, pp. 828–839.