\documentclass[jidm,a4paper]{jidm}

\jidmVolume{7}
\jidmNumber{1}
\jidmYear{16}
\jidmMonth{February}
\setcounter{page}{3}

\usepackage{ifxetex}
\ifxetex
\usepackage{fontspec}
%\setmainfont{DejaVu Serif}
\setmainfont{TeX Gyre Termes}
\else
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\fi

\widowpenalty10000
\clubpenalty10000

\usepackage{graphicx,url}
\usepackage{subcaption}
\usepackage{graphicx,url,multirow,tikz,latexsym,amsmath}
\graphicspath{{fig/}}
\usepackage[subfolder]{gnuplottex}
\usepackage{clrscode3e}

\usetikzlibrary{spy,shapes.misc,positioning}

\usepackage{float}
\floatstyle{plain}
\newfloat{algorithm}{t}{lop}
\floatname{algorithm}{Algorithm}

\makeatletter
\def\gnewcommand{\g@star@or@long\new@command}
\def\g@star@or@long#1{% 
  \@ifstar{\let\l@ngrel@x\global#1}{\def\l@ngrel@x{\long\global}#1}}
\makeatother

% define newalgo and others, to provide exact same names along the article
\newcommand{\newalgo}[3]{ %1=Acronym, 2=Description, 3=Citation
	\expandafter\gnewcommand\csname ALGA#1\endcsname{#1}
	\expandafter\gnewcommand\csname ALGD#1\endcsname{#2}
	\expandafter\gnewcommand\csname ALGC#1\endcsname{#3}
}
\newcommand{\descalgo}[1] {%
		{\emph{\expandafter\csname ALGD#1\endcsname} (\csname ALGA#1\endcsname) \cite{\csname ALGC#1\endcsname}}%
}

\markboth{Thiago B. Oliveira, F{\'a}bio M. Costa, Vagner J. S. Rodrigues}
{Distributed Execution Plans for Multiway Spatial Join Queries using Multidimensional Histograms}

\title{Distributed Execution Plans for Multiway Spatial Join Queries using Multidimensional Histograms}

\author{Thiago Borges de Oliveira, F{\'a}bio Moreira Costa, Vagner Jos{\'e} do Sacramento Rodrigues}
\institute{Instituto de Inform{\'a}tica -- Universidade Federal de Goi{\'a}s (UFG), Brazil \\ \email{thborges@ufg.br, \{fmc, vsacramento\}@inf.ufg.br}}

\begin{abstract}
Multiway spatial join is a common and heavyweight type of query for spatial data processing on relational database management systems. This article presents a complete solution to process this type of query in distributed systems. We proposed a cost-based optimizer for multiway spatial join queries, based on a novel use of multidimensional histograms, which are used to represent two metrics that describe a dataset: cardinality and size of the spatial objects, together with one feature that represents the location of the dataset partitions in the cluster. A new method for histogram construction is proposed, together with a formula to estimate the cost of multiway spatial join queries and to select good executions plans. A greedy algorithm is proposed to schedule query execution while minimizing both the makespan and the communication cost. We adapted the Clone Join algorithm to improve its execution time and to process multiway spatial join queries. The evaluation demonstrated that the new method for histogram construction performs better than a well-know method, and provides a better estimate for the cost of execution plans in 75\% of the queries evaluated. The cost estimated by the proposed method is shown to be sufficiently close to real execution times, and our complete methodology was able to use the estimated costs to select the best execution plan for all queries used in the experiments. It also provided a nearly exact ordering of the query execution plans with respect to execution time.
\end{abstract}

\category{H.2.8}{Database Management}{Database Applications}[Spatial databases and GIS]
\category{H.2.4}{Database Management}{Systems}[Distributed databases \and Query processing]

\keywords{Distributed Processing, GIS, Multiway Spatial Join, Spatial Data}

\begin{document}

\begin{bottomstuff}
This work was partially supported by CNPq, process number 473.939/2012-6.
\end{bottomstuff}

\maketitle


\section{Introduction}

The amount of spatial data has greatly increased with the popularization of GPS-enabled devices. Spatial data, such as geotagged images, the Internet of Things (IoT) and ``Smart Cities'' sensors data, open data and census data, are continuously collected and organized in thematic datasets to support decision-making, in order to improve business intelligence and logistics efficiency in both companies and government. Spatial data is a complex, multidimensional data type that is handled by a relational database management system (RDBMS) through the use of queries with spatial predicates~\cite{carvalho2014jidm}.
 
An important type of query is the spatial join query, which finds correlated objects in two or more datasets by applying some spatial predicate like intersection or proximity~\cite{brinkhoff1996pps}. Spatial joins can be simple, when processing only two datasets, or complex or multiway, when processing more than two datasets in the same query~\cite{Mamoulis2001a}. Multiway spatial join queries are important in several application fields, including geography (e.g., to find animal species living in preservation areas that were damaged by fire), VLSI (e.g., to find circuits that formulate a specific topological configuration)~\cite{Mamoulis2001a} and digital pathology imaging~\cite{aji2012towards} (e.g., to analyze topological images of the brain in order to check whether a certain cancer is progressing).

The processing of a multiway spatial join is more complex than the processing of a simple spatial join, because a query can be executed in many different ways, called execution plans. An execution plan is a definition of the order in which datasets will be combined and the algorithms that will be used to find the final result of the query. The execution plan has a huge impact on the processing cost of a query, and thus, in industrial instances of practical size, a significant amount of effort has to be dedicated to determine the best plan for a query on relational databases~\cite{mishra1992join}. However, this effort is mainly dedicated to non-distributed databases or to distributed databases with scalar columns. Spatial data impose particular challenges to join processing and, in general, algorithms designed for scalar data in relational databases do not apply to spatial data due to the absence of total ordering in multidimensional data~\cite{Jacox2007}. 

Another particular challenge of spatial data is the complexity of computational geometry algorithms~--~which evaluate the predicates over spatial data. Even queries on small datasets have a high processing cost due to the complexity and the different types of spatial objects (e.g., points, lines, polygons). One way of reducing query execution time is to process spatial queries in distributed systems by partitioning the datasets using spatial columns~\cite{Jacox2007}. The partitioning of datasets using spatial columns imposes significant challenges to distributed systems due to the skewed nature of spatial datasets. The selection of execution plans in distributed systems, besides considering local CPU and I/O costs, must take into account the proper load balance of both data and queries. The balance has to respect the communication cost among machines and the even distribution of the workload among processing nodes, with regard to both the bandwidth limit on the network interfaces and the CPU load. 

A significant effort in the literature is dedicated to the distributed processing of spatial join queries~\cite{patel2000cja,Chung2005,Jacox2007,Savio2013}. However, the main focus has been on processing simple spatial joins. To the best of our knowledge, there is no work reported in the open literature that proposes the selection of efficient execution plans for multiway spatial join query execution in distributed systems. It appears that any such selection should be based on an evaluation of network cost, CPU cost, and load balance among the machines, considering the partitioning of datasets using spatial columns.

Choosing good execution plans and efficiently processing multiway spatial join queries in distributed systems is a step towards moving spatial data analysis to scalable platforms, as has already happened to relational and unstructured data. However, new methods and algorithms to estimate the cost of an execution plan need to be specified, taking into consideration the specifics of spatial data and the characteristics of distributed systems. Processing spatial data analysis in such environments can greatly improve the capabilities of spatial data processing, especially in today’s scenario of cloud computing platforms, taking advantage of scalability, elasticity, and pay-as-you-go offers.

This article builds on our previous work~\cite{deOliveira2015sbbd}, with improvements on the methodology for cost estimation of query execution plans, and a new algorithm for query scheduling and load balancing in distributed systems. The entire set of proposed methods constitutes an optimizer to process multiway spatial join queries in distributed systems. The main contributions of the present article are:

\begin{itemize}
\setlength{\itemsep}{0pt}
\item The identification of the characteristics of datasets and data distribution which are relevant to the efficient processing of multiway spatial join queries in distributed systems;
\item The definition of a multidimensional histogram data structure to organize these characteristics, observing data skewness in real spatial datasets;
\item A new method to construct a multidimensional histogram that improves the estimate of query processing cost;
\item An algorithm that estimates the cost of execution plans in distributed systems using multidimensional histograms; and
\item A greedy algorithm to schedule data copying between cluster nodes, which reduces bandwidth usage while maintaining load balance in the cluster.
\end{itemize}

The remainder of this article is organized as follows. In Section~\ref{sec:revisao} we present a survey of multiway spatial join processing, describing key concepts, existing estimation techniques for plan selection and algorithms to process multiway join queries, as well as a discussion of how to adapt these methods for distributed systems and related work. In Section~\ref{sec:histograma}, we propose a multidimensional histogram data structure, methods for constructing such histograms, an algorithm for estimating network I/O, a data copying scheduling algorithm and a basis for the selection of execution plans. Section~\ref{sec:avaliacao} presents an evaluation of the proposed methods and discusses the results. Section~\ref{sec:conclusao} presents conclusions and suggested future work.


\section{Concepts and Related Work} \label{sec:revisao}

In this section we cover some concepts and existing related work. We start with the fundamentals of multiway spatial join and, then, describe its processing steps: plan cost estimation using multidimensional histograms, plan enumeration and selection, and existing algorithms to process join queries. Finally, we present related work. 

\subsection{Multiway Spatial Join}

Multiway spatial join is a technique for combining successive join algorithms to process (complex) spatial join queries with more than two input datasets~\cite{Papadias1999a}. A multiway spatial join can be represented by graph $G(V,E)$, where each node represents a dataset and the edges represent the join predicates~\cite{Papadias1999a,Mamoulis2001a}. 

Different execution plans can be used to process a multiway spatial join query. The number of ways a query can be processed was investigated by \citeN{Mamoulis2001a} for serial (i.e., non-parallel, non-distributed) processing. This number can be determined by the query type, the number of input datasets, and the number of different join algorithms that can be used at each step. 

Figures~\ref{fig:multiwaybushy}, \ref{fig:multiwayabushy}, and \ref{fig:multiwayleft} illustrate the alternative plans to process the chain query graph presented in Figure~\ref{fig:multiwayqgraph}. In a chain query, the graph $G(V,E)$ is a tree, with any two vertices connected by exactly one path. In Figure~\ref{fig:multiwaybushy}, the datasets are pairwise joined in a first step, producing two intermediate results to be joined in a second step. In Figure \ref{fig:multiwayabushy}, three datasets are joined in a single step and the intermediate result is then joined with the remaining data set ``animals''. In Figure~\ref{fig:multiwayleft}, the datasets ``fire'' and ``riversides'' are joined and the intermediate results are processed in subsequent steps with the other two datasets, one at a time.

\begin{figure}[tb]
	\centering
	\tikzset{level distance=0.8cm,
		 dataset/.style={circle,fill=black!30,draw},
		 plan/.style={circle,fill=black!10,draw},
		 node distance=0.8cm}
	\scriptsize
	\begin{tabular}{|c|c|c|c|}
	\hline

	\begin{subfigure}[b]{.22\textwidth}
		\centering
		\begin{tikzpicture}
  			\node[dataset] (1) [label=above left:animals] {};
			\node[dataset] (2) [below of=1, label=below left:areas] {};
			\node[dataset] (3) [right of=2, label=below right:riversides] {};
			\node[dataset] (4) [right of=1, label=above right:fire] {};
			\path
				(1) edge node {} (2)
				(2) edge node {} (3)
				(3) edge node {} (4);
		\end{tikzpicture}
		\caption{Query Graph}
		\label{fig:multiwayqgraph}
	\end{subfigure}

	&

	\begin{subfigure}[b]{.22\textwidth}
		\centering
		\begin{tikzpicture}[level 1/.style={sibling distance=1.5cm},level 2/.style={sibling distance=.9cm}]
  			\node[] (1) [] {$\Join$}
				child{node[] (2) [] {$\Join$}
					child{node[plan] (4) [label=below:animals] {}}
					child{node[plan] (5) [label=below:areas] {}}
				}
				child{node[] (3) [] {$\Join$}
					child{node[plan] (6) [label=below:fire] {}}
					child{node[plan] (7) [label=below:riversides] {}}
				}
			;
			\path
				(1) edge node {} (2)
				(1) edge node {} (3)
				(2) edge node {} (4)
				(2) edge node {} (5)
				(3) edge node {} (6)
				(3) edge node {} (7);
		\end{tikzpicture}
		\caption{Bushy plan}
		\label{fig:multiwaybushy}
	\end{subfigure}

	&

	\begin{subfigure}[b]{.23\textwidth}
		\centering
		\begin{tikzpicture}
  			\node[] (1) [node distance=1.5cm] {$\Join$};
			\node[] (2) [below left of=1] {$\Join$};
			\node[plan] (3) [below right of=1, label=below:animals] {};
			\node[plan] (5) [below left of=2, label=below:areas] {};
			\node[plan] (6) [below of=2, label=below:riversides] {};
			\node[plan] (7) [below right of=2, label=below:fire] {};
			\path
				(1) edge node {} (2)
				(1) edge node {} (3)
				(2) edge node {} (5)
				(2) edge node {} (6)
				(2) edge node {} (7);
		\end{tikzpicture}
		\caption{Another bushy plan}
		\label{fig:multiwayabushy}
	\end{subfigure}

	&

	\begin{subfigure}[b]{.22\textwidth}
		\centering
		\begin{tikzpicture}
  			\node[] (1) [] {$\Join$};
			\node[] (2) [below left of=1] {$\Join$};
			\node[plan] (3) [below right of=1, label=below:animals] {};
			\node[] (4) [below left of=2] {$\Join$};
			\node[plan] (5) [below right of=2, label=below:areas] {};
			\node[plan] (6) [below left of=4, label=below:fire] {};
			\node[plan] (7) [below right of=4, label=below:riversides] {};
			\path
				(1) edge node {} (2)
				(1) edge node {} (3)
				(2) edge node {} (4)
				(2) edge node {} (5)
				(4) edge node {} (6)
				(4) edge node {} (7);
		\end{tikzpicture}
		\caption{Left-deep plan}
		\label{fig:multiwayleft}
	\end{subfigure}

	\\
	\hline
	\end{tabular}

\caption{Alternate execution plans to process a multiway spatial join query.}
\label{fig:multiwayplan}
\end{figure}

All the execution plans for a given query preserve the same query semantics, but may take different amounts of time to produce the final result. The graph that represents a query is used as an input for a plan selection algorithm, also called an {\em optimizer}. The optimizer considers some aspects of the datasets and the associated costs of the algorithms to select an execution plan that determines: ($i$) how the datasets will be combined, ($ii$) what the processing order is, and ($iii$) which algorithms will be used in each step. Despite the name, an optimizer, in general, is not an exact algorithm. In the next section, we discuss how the computational cost of an execution plan can be determined.


\subsection{Estimating the Cost of Spatial Join for Plan Selection} \label{sec:estimativacusto}

Despite the complexity of estimating the cost of execution plans, some studies have proposed formulae that estimate the cost of simple spatial join queries~\cite{roh2010hierarchically,Sivasubramaniam2001,Fornari2006a}, as well as methods to combine these formulae to estimate the cost of complex spatial join queries \cite{mamoulis2001selectivity}. These formulae assume that spatial objects fill the spatial extent uniformly and take the I/O and CPU costs of the join algorithm to estimate the cost of an execution plan. The difficulty of using these formulae on real datasets was latter studied by the same authors \cite{Mamoulis2001a}, who conclude that, when used on real spatial datasets, the formulae can induce bad execution plans, especially in the presence of dataset skewness.

\citeN{mamoulis2001selectivity} studied the applicability of the formulae in small regions of the dataset. The authors proposed the use of an uniform histogram that divides the spatial extent of the dataset in a map with fixed-size cells, each of them mapping the density of small regions. In this histogram, the cardinality value of each cell $x$, denoted by $C(x)$, is defined based on the number of spatial objects that have the center of their Minimum Bounding Rectangle (MBR) within the limits of~$x$.

The main advantages of uniform histograms are $i$) simplified construction, $ii$) time efficiency to estimate queries, and $iii$) incremental maintenance \cite{cormode2012}. In distributed systems, this type of histogram also provides a simple way to partition the datasets, a step needed for distributed processing. A deficiency of uniform histograms is the estimation error when data is skewed. The error can be reduced by increasing the number of cells, although this increases the quantity of memory needed to store the histogram structure \cite{cormode2012}. There are other types of histograms for spatial data \cite{roh2010hierarchically,ray2014skew} that can be used to improve the estimation of skewed datasets. However, the complexity of estimating query costs with these histograms is greater, and their recursive nature creates important drawbacks in incremental maintenance and data partitioning.

An important challenge in the construction of multidimensional histograms for spatial data is how to hash the objects into cells. Although the MBR center method, proposed by \citeN{Mamoulis2001a}, is accurate for point data, other types of spatial objects, such as lines and polygons, bring additional challenges to histogram construction due to their spatial extent. A single line or polygon object can overlap more than one cell in the histogram and hashing it through the MBR center can lead to large errors in estimates.
 
In distributed systems, besides the cardinality of the datasets, the size of spatial objects (given by the number of points) and the place (server) where the partitions are located are also relevant to estimate the communication cost. These variables need to be stored as part of the histogram statistics. 


\subsection{Selecting an Execution Plan for Multiway Spatial Join} \label{sec:selecplanoexec}

The size of the space of possible plans to select from is nonlinear with respect to the number of datasets. \citeN{Mamoulis2001a} studied the number of possible plans, considering three different algorithms to process spatial join queries in a non-distributed system. The recurrence $P(n) = 1+2P(n-1) + \sum_{2\le k < n-1}^{} P(k)P(n-k)$, with $P(2) = 1$ gives the number of plans, $P(n)$, for a chain query with $n$ datasets \cite{Mamoulis2001a}. In asymptotic terms, $P(n) = \Omega(2^n)$. Cycle queries and clique queries have an even larger number of possible plans. As a result, if a query involves a sufficiently large number of datasets, an optimizer can take more time planning (through having to enumerate all possible plans and to compute their costs) than actually executing the query.

In distributed systems, besides the combination order of datasets, there are different strategies to copy data partitions during join execution. Some options are: $i$) copy smaller data partitions to the servers where the largest partitions are located; $ii$) copy partitions in such a way that network contention does not make processors sleep; and $iii$) copy partitions to maintain cluster balance. These choices can also be considered for the selection of the execution plan, further increasing the number of possible plans.

To quickly identify good execution plans while searching only a small fraction of the space, \citeN{Mamoulis2001a} proposed an heuristic that randomly transforms an execution plan (seed) using a set of pre-defined rules, such as associativity and commutativity, together with methods such as iterative improvement and simulated annealing. In their experiments, the heuristic method was able to find plans only slightly more expensive than the optimal execution plans found by the exhaustive method. The proposed algorithm is a good strategy for queries with a large number of datasets. 



\subsection{Algorithms for the Distributed Processing of Multiway Spatial Join} \label{sec:particionamentodist}

\newalgo{PBSM}{Partition-based Spatial Merge Join}{Dewitt1996}
\newalgo{RJ}{R-Tree Join}{brinkhoff1993eps}
\newalgo{SISJ}{Slot Index Spatial Join}{Mamoulis1999}
\newalgo{HJ}{Spatial Hash Join}{lo1996spatial}

\newalgo{RSPR}{Replicated Semi-packed Parallel R-Tree}{Mutenda1999}
\newalgo{PASJ}{Proximity Area Spatial Join}{Savio2013}
\newalgo{CJ}{Clone Join}{patel2000cja}
\newalgo{SJ}{Shadow Join}{patel2000cja}
\newalgo{NPSJ}{Non-blocking Parallel Spatial Join}{Naughton2002a}
\newalgo{DST}{Distributed Synchronous Traversal}{Cunha2015}
\newalgo{ST}{Synchronous Traversal}{Papadias2001a}

Processing spatial join queries in distributed or parallel systems requires some data partitioning methods that divides the dataset objects into groups called data partitions or cells, which are assigned to servers during the filtering or refinement step. The methods for data partitioning can be classified into two main categories: 

\begin{enumerate}
\setlength{\itemsep}{0pt}
\item Disjoint space partitioning: uses a grid to split the space extent. Each cell of the grid groups the spatial objects according to their intersection with the cells. The cells do not intersect each other. Objects that intersect more than one cell are replicated.
\item Non-disjoint space partitioning: partitions can overlap each other to accommodate the extent of the objects that intersect them. Spatial objects are not replicated. An example of this type of partitioning is the set of MBR's on a given level of an R-Tree index.
\end{enumerate}

For non-disjoint space partitioning, the relevant distributed algorithms are \descalgo{RSPR}, \descalgo{DST}, and \descalgo{PASJ}. All of them use distributed R-Tree indexes to process the spatial join. For disjoint space partitioning, the relevant algorithms are \descalgo{CJ}, \descalgo{SJ}, and \descalgo{NPSJ}. CJ is the simplest of them. It uses a fixed grid to partition the datasets and replicate spatial objects that intersect more than one grid cell. The difference between CJ and SJ is that the latter uses a more sophisticated technique to reduce object replication. The NPSJ algorithm also uses a fixed grid and replicates objects as CJ. However, it is designed as a non-blocking algorithm that produces results early on the execution. Another difference between NPSJ and the others is that it creates local R-Trees on each data partition and produces the results using an \descalgo{RJ}, while CJ and SJ use \descalgo{PBSM}.

An issue in distributed spatial join processing is the distribution of the data partitions. The simplest and most used distribution method is round-robin. It distributes the partitions evenly between the servers, favoring load balance in the cluster. Other methods that favor the co-location of spatial data have also been proposed, such as the Proximity Area method \cite{Savio2013} for non-disjoint partitioning, that distributes the partitions based on their location, favoring reduction of network bandwidth usage during spatial join execution.

In a round-robin distribution, data partitions that are located in the same geographic region, but from different datasets, will often be distributed to different servers. The join algorithm thus needs to copy the intersecting partitions from different servers to execute the join, but the predicate checking will be more balanced among the nodes of the cluster. As discussed in \citeN{Savio2013}, distributing the partitions based on their location will reduce network bandwidth usage at the expense of compromising the load balancing of spatial join execution. 

\subsection{Clone Join and Reference Point Method} \label{sec:cj}

\descalgo{CJ} is a distributed spatial join algorithm for joining two non-indexed datasets. The data partitioning used by the algorithm is a grid of disjoint cells, each cell representing a small part of the spatial extent. The spatial objects are assigned to each cell based on the intersection between their MBRs and the cell boundaries. Objects that intersect more than one cell are replicated. Each cell, and the corresponding spatial objects, form a data partition and are assigned to one server using a round-robin distribution or a similar hash function applied to the cell number. The join operation is performed in parallel in each server, using a local partition-based spatial merge join (PBSM)~\cite{Dewitt1996}, which is a parallel spatial join algorithm for shared-memory parallel systems. 

\citeN{patel2000cja} assume, at the beginning of the execution of the join, that the two datasets are partitioned with the same grid. This can be a problem with a database system, as distribution will occur every time a join is executed using the dataset. The datasets can also be inserted at distinct times or have different geographical space extents.

Because of object replication, an additional step is necessary to eliminate duplicated results reported by the refinement. A duplicated result is reported whenever two spatial objects, that intersect each other, are replicated on two or more cells, and the cells are distributed to different servers. Each server will separately report the intersection between the objects. \citeN{patel2000cja} used a {\em distinct} verification at the end of the algorithm to eliminate the duplicated results. This is a deficiency of the algorithm because the distinct verification is a costly operation in distributed systems due to its network-bound behavior. Objects with large spatial extents (such as lines and polygons) naturally increase the number of replicated objects, and, consequently, also increase the usage of resources (notably, bandwidth and execution time).

\citeN{dittrich2000data} proposed a solution for the result duplication problem called the reference point method. The method consists of identifying the possible cells that will report duplicate results and allowing only one of them to do so. The method was originally proposed for the PBSM algorithm and later adapted to NPSJ~\cite{Naughton2002a}, a spatial join algorithm for distributed systems. However, it can be applied to CJ as well, as is detailed in Section~\ref{sec:cj_with_rpm}.



\subsection{Related Work} \label{sec:correlatos}

Despite the importance of multiway spatial join queries and the widespread use of spatial datasets, little work has been proposed to efficiently process multiway spatial join queries in distributed systems. In this section, we discuss the most relevant research efforts in the field. Overall, they propose isolated contributions that can be combined to form a complete solution for distributed multiway spatial join query processing.

To the best of our knowledge, \citeN{Mamoulis2001a} represent the most relevant study in the literature for the non-distributed processing of multiway spatial join. The authors make an in-depth evaluation of algorithms, data partitioning techniques, plan cost estimation, plan selection and query execution. In order to build a solution for distributed systems, several of the algorithms and methods they propose were adapted in our work. For simple spatial joins with only two datasets \citeN{Fornari2006a} proposed an optimizer for spatial join queries in non-distributed systems and also a rule-based optimizer~\cite{Fornari2007}.

A method for the distributed processing of multiway spatial join queries is presented in~\citeN{aji2012towards} using a MapReduce framework. The authors propose a data partitioning method specific for MapReduce, and a distributed algorithm that processes the multiway query in a single step, combining all the datasets at once. The processing of all datasets in one step is due to a limitation imposed by the framework, which requires that intermediate results are stored, incurring in high I/O costs. Several studies that use the MapReduce framework do not evaluate the alternative plans for executing a query~\cite{zhang2009sjmr,Zhong2012a,gupta2013processing}. \citeN{gupta2013processing} discusses the spatial data partitioning and replication for the MapReduce framework. In order to reduce communication between nodes in the cluster, \citeN{zhang2009sjmr} and \citeN{Zhong2012a} propose MapReduce algorithms and data partitioning techniques to process simple spatial join with two datasets. \citeN{Cunha2015} propose an algorithm to process multiway spatial join queries in distributed systems, also combining all datasets at once, but with the independent framework for spatial data processing proposed by \citeN{de2011dsi} and \citeN{Savio2013}. 

While the most popular strategy to process multiway spatial join in distributed systems is the combination of all datasets in a single step following the \descalgo{ST} algorithm (which was originally proposed for non-distributed processing), most of the studies were limited by the use of the MapReduce framework. To the best of our knowledge, there are no studies that compare this strategy with the selection of the execution plan by considering the entire space of possible plans. For non-distributed systems, \citeN{Mamoulis2001a} made this comparison and concluded that the best plan for real datasets is a hybrid plan, with some parts being processed with ST and others with a join algorithm for two datasets at a time. 

A number of research efforts~\cite{patel2000cja,Naughton2002a,Chung2005,Savio2013} propose distributed algorithms for spatial join using two datasets at a time. Despite the fact that they are not algorithms for multiway spatial join, these algorithms can be adapted to process multiway spatial join queries in successive steps, as it has been proposed by~\citeN{Mamoulis2001a} for non-distributed systems.

As these studies focus on distributed processing of multiway spatial join in only one step or are specific to two datasets, we have not found comparable methods for distributed plan cost estimation and distributed plan selection. In this article, we propose a new method for cost estimation (see Section~\ref{sec:histograma}) based on a new multidimensional histogram for distributed systems, together with a new method for plan selection. To the best of our knowledge, this is the first attempt to propose and evaluate such methods for processing multiway spatial join queries in distributed systems, evaluating the cost of alternate execution plans.

%\subsection{Design Choices}

%In our work, we are using only the exhaustive method and left for future work the implementation of this randomized algorithm.

%Due to spatial join being a CPU-bound algorithm, we implement the round-robin method and execute experiments to measure the impact of communication on the total execution time of the queries. Despite this choice, the methods we propose in Section~\ref{sec:histograma} were designed for use with any data distribution, as is demonstrated in Section~\ref{sec:estimandocomunicacao}.

%Following the results reported by \citeN{patel2000cja}, we use a disjoint space partitioning. Each grid cell groups and assigns the spatial objects to the servers. For join execution, we use the CJ algorithm with the duplication avoidance technique, based on the Reference Point method. The CJ algorithm is simpler than NPSJ and duplication avoidance also transforms it into a non-blocking algorithm. The next section covers the CJ algorithm and the Reference Point Method in more detail.

%We improved the 2D histogram to represent this two new values and developed an algorithm to estimate the cost of a spatial join, based on the three metrics: cardinality, points, and place. The details are discussed in Section~\ref{sec:histograma}.


%We carried out experiments with histograms with fixed cells and demonstrated promising results even with real and skewed datasets, as shown in the evaluation presented in Section~\ref{sec:avaliacao}.



\section{Multidimensional Histogram and Optimizer}\label{sec:histograma}

In this section we describe the proposed data structure, called a multidimensional histogram, along with its access method. This data structure combines parameters that describe the datasets and their allocation in the distributed system.
%In this section we describe how to combine the parameters of datasets and data distribution in the data structure and access method denominated multidimensional histogram.

A multidimensional histogram divides the spatial extent of a dataset using a grid with cells of fixed size. For each cell, it records: $i$) the number of spatial objects within the boundaries of the cell (the cardinality of the cell), $ii$) the combined size of the objects in the cell (in terms of their total number of points), and $iii$) the place (server in the cluster) where the cell is located. The first proposed metric is used to estimate the cardinality of the join output, in a similar way to the non-distributed execution of spatial join queries. The second proposed metric is used to estimate the intermediate histograms used to select an execution plan for a spatial join query. We use the three values of $i$), $ii$) and $iii$) together to estimate network bandwidth usage, as it is detailed in Section~\ref{sec:estimandocomunicacao}.

Figure~\ref{fig:hmunicipios} shows a visualization of a multidimensional histogram generated for a dataset that represents the political limits of Brazilian municipalities. The histogram captures important aspects of the dataset. At the top left corner of Figure~\ref{fig:hmunicipioscard}, the cardinality is low because this region corresponds to the Amazon Forest, where there are few large municipalities, resulting in lower cardinalities for the cells. In contrast, Figure~\ref{fig:hmunicipiospoints} shows the points metric for the same municipalities. The number of points is much higher than in other regions because the area(extent) of the objects is larger and more points are needed to represent their contour. Figure~\ref{fig:hmunicipiosplaces} is a map of the location of spatial data in the cluster. For this sample histogram, the legend indicates the location of each cell (in servers 1 to 8). A value of zero (black) represents a cell that contains no objects and does not need to be sent to a server.

\begin{figure}[tb]
	\centering
	\scriptsize
	\begin{subfigure}[b]{.33\textwidth}
		\centering
		\begin{gnuplot}[terminal=epslatex, terminaloptions={size 3.5cm,2.7cm monochrome}]
		unset key
		unset xtics
		unset ytics
		set ztics 10
		set xlabel "Longitude" offset 0,0.3
		set ylabel "Latitude" offset -0.9,-1
		set zlabel "Cardinality" rotate by 90 offset 4,0
		set lmargin at screen 0.2
		set bmargin at screen 0.25
		set rmargin at screen 0.9
		set tmargin at screen 0.99
		set view 40,20
		set hidden3d
		set ticslevel 0
		set dgrid3d 50,50 qnorm 2
		splot 'fig/hist-card-municipios.dat' with lines
		\end{gnuplot}
		\caption{}
		\label{fig:hmunicipioscard}
	\end{subfigure}
	\begin{subfigure}[b]{.33\textwidth}
		\centering
		\begin{gnuplot}[terminal=epslatex, terminaloptions={size 3.5cm,2.7cm monochrome}]
		unset key
		unset xtics
		unset ytics
		set xlabel "Longitude" offset 0,0.3
		set ylabel "Latitude" offset -0.9,-1
		set zlabel "Points" rotate by 90 offset 2.5,0
		set lmargin at screen 0.2
		set bmargin at screen 0.25
		set rmargin at screen 0.9
		set tmargin at screen 0.99
		set view 40,20
		set hidden3d
		set ticslevel 0
		set dgrid3d 50,50 qnorm 2
		splot 'fig/hist-points-municipios.dat' with lines
		\end{gnuplot}
		\caption{}
		\label{fig:hmunicipiospoints}
	\end{subfigure}
	\begin{subfigure}[b]{.32\textwidth}
		\begin{gnuplot}[terminal=epslatex, terminaloptions={size 3.5cm,2.7cm monochrome}]
		unset key
		#set xtics rotate
		unset xtics
		unset ytics
		set xlabel "Longitude" offset 0,2
		set ylabel "Latitude" rotate by 90 offset 1,0
		set lmargin at screen 0.2
		set bmargin at screen 0.2
		set rmargin at screen 0.8
		set tmargin at screen 0.8
		set pm3d map
		set dgrid3d 50,50 qnorm 2
		splot 'fig/hist-places-municipios-rr.dat'
		\end{gnuplot}
		\caption{}
		\label{fig:hmunicipiosplaces}
	\end{subfigure}

	\caption{Multidimensional histogram for a two-dimensional dataset that represents the political limits of Brazilian municipalities. The graphs represent: the cardinality (\subref{fig:hmunicipioscard}) and the points (\subref{fig:hmunicipiospoints}) metrics; and the location of partitions (\subref{fig:hmunicipiosplaces}).}
	\label{fig:hmunicipios}
\end{figure}


\subsection{Proportional Overlap Method} \label{sec:fractionary}

To build a histogram, each object in a dataset needs to be hashed into the histogram grid in order to find the cell or the set of cells that it intersects. The set of objects hashed to a cell forms a data partition. We use the terms {\em cell} and {\em data partition} interchangeably in the remaining text.

We propose a new method to build the multidimensional histogram, called the Proportional Overlap method. Instead of using the center of the MBR to hash a spatial object into the histogram, as proposed by \citeN{Mamoulis2001a}, we use the proportional covering area of each MBR. To calculate the cardinality metric, we increment the value stored in each cell based on the proportion of two areas ($a/b$), where $a$ is the area of intersection between the MBR of the spatial object and the histogram cell boundaries and $b$ is the area of the MBR. Formally, let $H_r$ be the histogram of the dataset $R$, and let $o \in R$ be an object of the dataset. The cardinality $C(x)$ of a cell $x \in H_r$ is given in~(\ref{eqn:cardinality}). Figure~\ref{fig:pom_and_hc} illustrates the hash process of a spatial object in a histogram. In the MBR center method (left), the center of the object's MBR identifies the histogram cell. The cardinality metric of the identified cell is thus increased by unity. In contrast, the proportional overlap method (right) increments the cardinality of each cell that intersects the object's MBR by adding the proportion obtained in~(\ref{eqn:cardinality}).%
%
\belowdisplayskip=-5pt%
\begin{align}\label{eqn:cardinality}
C(x) = \sum_{o \in R} area(MBR(o) \cap MBR(x)) / area(MBR(o))
\end{align}%
%
\begin{figure}[tb]
\centering
  \begin{tikzpicture}[scale=.3,every node/.style={minimum size=1cm},on grid]
		
    %slanting: production of a set of n 'laminae' to be piled up. N=number of grids.
    \begin{scope}[
    	yshift=0,every node/.append style={
    	    yslant=0.5,xslant=-1},yslant=0.5,xslant=-1
    	             ]
        \fill[white,fill opacity=.9] (0,0) rectangle (5,5);
        \draw[black,very thick] (0,0) rectangle (5,5);
        \draw[step=5mm, black] (0,0) grid (5,5);
        \fill[red,fill opacity=.1] (1.7,1.7) rectangle (3,2.5); %shadow
        \fill[red,fill opacity=.5] (2,2) rectangle (2.5,2.5); %shadow
        \fill[red] (2.35,2.05) circle (2pt);
    \end{scope}
    	
    \begin{scope}[
    	yshift=90,every node/.append style={
    	yslant=0.5,xslant=-1},yslant=0.5,xslant=-1
    	]
    	\fill[white,fill opacity=.9] (0,0) rectangle (5,5);
    	%\draw[step=10mm, black] (1,1) grid (4,4);
    	%\draw[black,very thick] (1,1) rectangle (4,4);
		\fill[fill opacity=.9,blue] (1.7,1.7) -- (1.7,1.9) -- (2,2) -- (2,2.5) -- (3,2.5) -- (3,2) -- (1.7,1.7);
        \draw[red, dashed] (1.7,1.7) rectangle (3,2.5);
        \fill[red] (2.35,2.05) circle (2pt);
    	\draw[dashed,black] (0,0) rectangle (5,5);
    \end{scope}

    %putting arrows and labels:
    \draw[-latex] (6.2,2) node[right]{Histogram}
         to[out=180,in=70] (4,2);

    \draw[-latex](5.9,5)node[right]{Spatial Extent}
        to[out=180,in=70] (3.6,5);

    \draw[-latex](3.5,7)node[right]{Spatial Object}
        to[out=180,in=60] (0.7,5.5);

	\draw[-latex](6,6)node[right]{MBR}
        to[out=180,in=45] (0.9,5.3);

	\draw[] (0.0,-0.5) node[right]{}; %just to align the figures
  \end{tikzpicture}%
  \begin{tikzpicture}[scale=.3,every node/.style={minimum size=1cm},on grid]
		
    %slanting: production of a set of n 'laminae' to be piled up. N=number of grids.
    \begin{scope}[
    	yshift=0,every node/.append style={
    	    yslant=0.5,xslant=-1},yslant=0.5,xslant=-1
    	             ]
        \fill[white,fill opacity=.9] (0,0) rectangle (5,5);
        \draw[black,very thick] (0,0) rectangle (5,5);
        \draw[step=5mm, black] (0,0) grid (5,5);
        \fill[red,fill opacity=.3] (1.7,1.7) rectangle (3,2.5); %shadow
        %\fill[fill opacity=.4,red] (1.7,1.7) -- (1.2,1.7) -- (2,2) -- (2,2.5) -- (3,2.5) -- (3,2) -- (1.7,1.7); % object
    \end{scope}
    	
    \begin{scope}[
    	yshift=90,every node/.append style={
    	yslant=0.5,xslant=-1},yslant=0.5,xslant=-1
    	]
    	\fill[white,fill opacity=.9] (0,0) rectangle (5,5);
    	%\draw[step=10mm, black] (1,1) grid (4,4);
    	%\draw[black,very thick] (1,1) rectangle (4,4);
		\fill[fill opacity=.9,blue] (1.7,1.7) -- (1.7,1.9) -- (2,2) -- (2,2.5) -- (3,2.5) -- (3,2) -- (1.7,1.7);
        \draw[red, dashed] (1.7,1.7) rectangle (3,2.5);
        %:\fill[red] (2.35,2.05) circle (3pt);
    	\draw[dashed,black] (0,0) rectangle (5,5);
    \end{scope}

    \draw[-latex] (2.2,-0.3) to[out=70,in=50] (1.1,2.3);
    \draw[] (-5.0,-0.5) node[right]{$area(MBR(o) \cap MBR(x)) / area(MBR(o))$};

  \end{tikzpicture}

  \caption{A spatial object being hashed using the MBR Center Method (left) and the Proportional Overlap Method (right). In the equation, $o$ is a spatial object and $x$ is a histogram cell.}
  \label{fig:pom_and_hc}

\end{figure}

Calculation of the points metric, in turn, is performed by summing the number of points of each of the objects that intersects a cell. This reflects the behavior of the CJ algorithm, which replicates an object in each of its intersecting cells. The number of points $P(x)$ in a cell $x$ is thus given in~(\ref{eqn:points}).%
\belowdisplayskip=0pt%
\begin{align}\label{eqn:points}
P(x) = \sum_{o \in R\;|\;MBR(o)\;\cap\;x \neq \emptyset} Points(o)
\end{align}

After the hashing process, we determine the location of the data partitions -- their physical distribution in the cluster. The location feature reflects the initial distribution of the partitions in the cluster. As such, it is updated whenever the partitions are replicated as part of the join processing.

\subsection{Clone Join with the Reference Point Method} \label{sec:cj_with_rpm}

As described in Section~\ref{sec:cj}, \descalgo{CJ} is designed to process simple spatial join queries, which involve only two datasets. In this section we describe an adaptation of CJ to efficiently process multiway spatial join queries. 

The CJ algorithm has two drawbacks when dealing with multiway spatial join query execution: $i$) the {\em distinct} operator used to prevent duplicated results has a high communication cost, and $ii$) the {\em distinct} operator prevents pipeline execution of multiway joins, i.e., it acts as a barrier, delaying the start of the next step of a plan until the end of the previous step, thus causing slow performance.

To efficiently process a multiway spatial join query using CJ, we can use the reference point method \cite{dittrich2000data} to identify duplicated results during the predicate checking, rather than as a separate step at the end. The predicate checking part of CJ consists of a simple loop that iterates over each pair of candidate spatial objects \cite{patel2000cja}. The reference point method consists of defining a reference point for each candidate pair of spatial objects and verifying whether or not it is within the boundaries of the cell being processed \cite{dittrich2000data}. As the cell boundaries are disjoint in CJ, only the cell containing the reference point will report the candidate pair. The required changes in the CJ algorithm are:

\begin{enumerate}
\setlength{\itemsep}{0pt}
\item Distribute and maintain cell boundaries as metadata on each server. Such metadata will be used to check if the reference point is inside a cell belonging to the server that is checking the predicate;
\item Change the CJ predicate checking algorithm:
	\begin{enumerate}
	\item Whenever a pair of objects has to be verified, compute the reference point using the MBRs of the two objects, $r$ and $s$: $rp = (max(r.xl, s.xl), min(r.yh, s.yh))$, where $xl$ is the lowest $x$ value, and $yh$ is the highest $y$ value for each of the two MBRs;
	\item Check if the reference point is within the boundaries of the cell that is being processed: if the reference point is not within the cell boundaries, ignore it (as it will be reported by another server). Otherwise, perform the predicate checking for the pair.
	\end{enumerate}
\item Remove the distinct operator at the end of the algorithm.
\end{enumerate}

The main advantage of these changes are that no communication is necessary to define the reference point. The only increase in communication is in data distribution, when the cell boundaries are transferred. This is a constant cost increase, however, as data partitions need to be transferred anyway. This cost is also compensated, as no distinct verification is needed at the end of the algorithm. The previous complexity of this verification is a function of the number of replicated spatial objects. Furthermore, the modified algorithm also has the non-blocking characteristic proposed in more recent algorithms, such as the \descalgo{NPSJ}, which further improves multiway spatial join query execution due to early reporting of results in subsequent steps.

%%THIS SECTION NEEDS SOME RE-WRITING, ESPECIALLY FROM THE THIRD PARAGRAPH FOR CONFORMANCE WITH THE CURRENT VERSION OF THE ALGORITHM (WAIT FOR THIAGO TO DO THE RE-WRITING, THEN FINISH REVIEWING THE SECTION (FMC)
\subsection{Query Scheduling and Plan Cost Estimation using Multidimensional Histogram} \label{sec:estimandocomunicacao}

To process multiway spatial join queries in a distributed system, we need to move the partitions of the involved datasets, which are already distributed in the cluster. This movement should be performed according to the semantics of the query predicate. In this section, we describe a method to estimate the cost of data movement and the overall cost of an execution plan based on the creation and use of intermediate histograms.

The data movement (communication) cost of a plan can be estimated using the metrics and features of the multidimensional histograms derived from the involved datasets. This is performed as follows. Let $H_a$ and $H_b$ be two histograms for datasets $a$ and $b$, respectively; let $C(c_x)$ be the cardinality of cell $c_x$; let $P(c_x)$ be the total number of points in cell $c_x$; and let $L(c_x)$ be the set of servers where cell $c_x$ is stored. An intermediate histogram, along with the communication cost, can be obtained using Algorithm~\ref{alg:communicationcost}, whose code is given later.

Algorithm~\ref{alg:communicationcost} performs three tasks: $i$) it determines a data copying strategy based on the pre-distribution of the datasets; $ii$) it estimates the communication cost of joining $a$ and $b$; and $iii$) estimate an intermediate histogram, $H_r$, for the result of the spatial join between datasets $a$ and $b$. This intermediate histogram will be used to estimate the next join steps of the multiway spatial join query. To construct the $H_r$ histogram, the algorithm evaluates each pair of cells from $H_a$ and $H_b$ that intersect each other. The algorithm estimates the metrics of this histogram using the proportional overlap method (lines \ref{alg_dens_start}-\ref{alg_dens_end}), and fills the resulting histogram with the estimated values (lines \ref{alg_comm_fill_cr} to \ref{alg_comm_next_points}).

The algorithm uses a greedy strategy to determine how data copying will be performed and to define the location feature of the $H_r$ histogram (lines \ref{alg_comm_start} to \ref{alg_comm_end}). The $\proc{ChooseMinCommServer}$ procedure (line \ref{alg_comm_start}) returns the server that requires the lowest communication cost to locate the cell pair $ca$ and $cb$ in the same place: If $L(ca) \cap L(cb) \neq \emptyset$, the two cells have a server in common and this server is returned; otherwise, a server of $L(ca)$ is returned if $P(ca) < P(cb)$, or vice-versa. The standard deviation of values in the $Pnts$ vector gives an estimate of the load balance for the execution plan. If the difference between the servers with the smallest and the largest number of points reaches a certain $\id{tradeoff}$ (line \ref{alg_comm_tradeoff}), instead of choosing the server with smallest amount of IO, the server with the smallest number of points is chosen (line \ref{alg_comm_tradeoff_choose}). This strategy attempts to maintain low communication cost, but if the balance in the cluster drops, communication is sacrificed at the expense of maintaining load balance. The number of points to be compared is then updated in the vector $Pnts$ (line \ref{alg_comm_pnts_amnt}), and the location of $ca$ and $cb$ cells are updated based on the chosen server (lines \ref{alg_comm_upd_ca} and \ref{alg_comm_upd_cb}). The communication cost is updated based on the number of points to be copied and the size of each point in the system (two double floating points)\footnote{We designed the system transfer protocol to send spatial object data in binary format. Its overhead is a constant proportional to the number of cells and objects to be transferred.}. This is performed for each server in the cluster, and stored in the $IO$ vector (lines \ref{alg_comm_copy_amnt} to \ref{alg_comm_end}). The return of the algorithm (line \ref{alg_comm_return}) is used for plan selection (in Section~\ref{sec:escolhaplanodst}).

\begin{algorithm}[tb]
\begin{codebox}
\Procname{$\proc{DMSJ\_StepCostsAndHr}(H_a, H_b)$}
\li $H_r \gets$ new empty histogram, based on the cell limits of $H_a$ and $H_b$
\li $IO[] \gets 0$ \label{alg_iovec}
\li $Pnts[] \gets 0$ \label{alg_pntsio}
\li \For each $ca \in H_a$, $cb \in H_b$, such that $ca \cap cb \neq \emptyset$
\li \Do $d_a \gets area(ca \cap cb) / area(ca) * C(ca)$ \label{alg_dens_start}
\li		$d_b \gets area(ca \cap cb) / area(cb) * C(cb)$
\li		$p_a \gets P(ca) / C(ca) * d_a$
\li 	$p_b \gets P(cb) / C(ca) * d_b$ \label{alg_dens_end}
\li		$cr \gets$ the cell of $H_r$ to be filled, based on the boundaries of $ca, cb$ \label{alg_comm_fill_cr}
\li		$C(cr) \gets d_a * d_b$
\li		$P(cr) \gets p_a$ or/and $p_b$, based on the predicate for the next step \label{alg_comm_next_points}
%\li		\Comment Determine the server where the cell pair will be processed
\li 	$s = \proc{ChooseMinCommServer}(IO, ca, cb)$  \label{alg_comm_start}
\li		$balance = (\proc{Max}(Pnts) - \proc{Min}(Pnts)) / \proc{Max}(Pnts)$
\li		\If $balance > \id{tradeoff}$ \label{alg_comm_tradeoff}
\li			\Then $s \gets \proc{MinIdx}(Pnts)$ \label{alg_comm_tradeoff_choose}
		\End
\li		$Pnts[s] \mathrel{{+}{=}} P(ca) + P(cb)$ \label{alg_comm_pnts_amnt}
\li		$L(ca) = L(ca) \cup s$ \label{alg_comm_upd_ca}
\li		$L(cb) = L(cb) \cup s$ \label{alg_comm_upd_cb}
\li 	\For $s \gets 1$ \To $servers$ \label{alg_comm_copy_amnt}
\li 	\Do \If $s \notin L(ca)$
\li			\Then $IO[s] \mathrel{{+}{=}} P(ca)$
			\End
\li			\If $s \notin L(cb)$
\li			\Then $IO[s] \mathrel{{+}{=}} P(cb)$ \label{alg_comm_end}
			\End
		\End
\li return $IO, Pnts, H_r$ \label{alg_comm_return}
\End
\end{codebox}
\caption{This algorithm determines a schedule to execute a spatial join query between datasets $a$ and $b$, returning both the CPU and network costs, and creates an intermediate histogram, $H_r$.}
\label{alg:communicationcost}
\end{algorithm}


\subsection{Selection of Distributed Execution Plans} \label{sec:escolhaplanodst}

Execution of the spatial predicate of a multiway spatial join query uses CPU intensive algorithms from computational geometry as a part of the computation of the query results. In a distributed system, this CPU cost should be evenly shared among the computing nodes in order to reduce the total execution time of the query (i.e., the makespan). Thus, we need to consider cluster balancing while selecting a plan to execute a query.

We use an exhaustive method to enumerate all the possible execution plans for a query. Let $P$ be the set of plans. For each execution plan $p \in P$, Algorithm~\ref{alg:plancost}, given later, is invoked to compute the cost of all steps in the plan. For the first step, the algorithm uses the histograms of two datasets, $p.steps[0].H_a$ and $p.steps[i].H_b$ (lines \ref{alg_cost_firststep_start} and \ref{alg_cost_firststep_end}), as input parameters to invoke Algorithm~\ref{alg:communicationcost} (line \ref{alg_call_algstepcost}). The returned $Pnts$ vector is aggregated in the $PlanPnts$ vector (lines \ref{alg_cost_agg_start} and \ref{alg_cost_agg_end}), and the intermediate histogram, $H_r$, is considered as the first histogram in the next step (line \ref{alg_cost_set_leftH}). Finally, in line \ref{alg_cost_set_allplans}, the procedure stores the $PlanPnts$ vector in $AllPlans[p]$, which is used to select the best plan for the query.

We suggest that the best distributed execution plan for a query is the plan with cost $O$ defined in~(\ref{eqn:optimal}). The expression $Max(AllPlans[p])$ in~(\ref{eqn:optimal}) selects the largest value of $PlanPnts$ for each execution plan $p$. This value refers to the server with the largest estimated load. This server will most likely be the bottleneck because the execution of the predicate checking algorithm will last longer in it. To conclude the selection, the expression $Min(...)$ in~(\ref{eqn:optimal}) selects from among the set of maximum values, the plan with the least number of points. In other words, we select the execution plan for which the bottleneck server has the lowest estimated load compared to the bottleneck servers from all the other plans. Thus, according to the estimate, the selected plan will be the one with the shortest makespan.%
%
\begin{align}\label{eqn:optimal}
		O = Min(\bigcup_{p \in P} Max(AllPlans[p]))
\end{align}%
%
The execution time need to select a plan using Algorithm~\ref{alg:plancost} and~($\ref{eqn:optimal}$) is just a small fraction of the total execution time of the sample queries used in the evaluation (see Section~\ref{sec:avaliacao}). However, this execution time increases substantially for queries that involve a large number of datasets, as the number of plans grows exponentially. For such queries, the use of pruning methods for relational queries, as proposed by~\citeN{ioannidis1990randomized} and latter adapted for spatial queries by~\citeN{Mamoulis2001a}, alleviates the problem. Such methods work by pruning the set of plans that need to be evaluated. Thus, our methodology can still be applied to the restricted set of plans that results from pruning. 

\begin{algorithm}[tb]
\begin{codebox}
\Procname{$\proc{DMSJ\_PlanCost}(p)$}
\li $H_{left} \gets p.steps[0].H_a$ \label{alg_cost_firststep_start}
\li $PlanPnts[] = 0$
\li \For $i \gets 0 \To \id{length}(p.steps)$
\li \Do  $H_{right} \gets p.steps[i].H_b$ \label{alg_cost_firststep_end}
\li      $\_, Pnts, H_r \gets \proc{DMSJ\_StepCostsAndHr}(H_{left}, H_{right})$ \label{alg_call_algstepcost}
\li      \For $s \gets 1 \To servers$ \label{alg_cost_agg_start}
\li      \Do $PlanPnts[s] \mathrel{{+}{=}} Pnts[s]$ \label{alg_cost_agg_end}
         \End
\li      $H_{left} \gets H_r$ \label{alg_cost_set_leftH}
    \End
\li $AllPlans[p] \gets PlanPnts$ \label{alg_cost_set_allplans}
\End
\end{codebox}
\caption{Algorithm to compute the cost of a query execution plan $p$.}
\label{alg:plancost}
\end{algorithm}


\section{Evaluation}\label{sec:avaliacao}

To evaluate our proposed methodology, we chose a set of real spatial datasets obtained from the Brazilian Institute of Geography and Statistics\footnote{\url{www.ibge.gov.br}} (IBGE), from the LAPIG Laboratory\footnote{Image Processing and Geoprocessing Laboratory -- \url{www.lapig.iesa.ufg.br/lapig/}} of the Institute of Social and Environmental Studies (IESA) at UFG, and from the current version of the Digital Chart of the World \footnote{\url{http://gis-lab.info/qa/vmap0-eng.html}}(DCW). The selected datasets and their characteristics are detailed in Table~\ref{tab:datasets}. Datasets with small cardinality were selected to conduct experiments with selective execution plans, in which one of the datasets strongly restricts the query output. On the other hand, datasets comprising complex spatial objects, such as rivers and roads, were selected due to the difficulty and inaccuracy of their representation in histograms.


The queries used in the experiments are listed in Table~\ref{tab:consultas}. This set of queries involves all of the datasets and represents multiway spatial joins that are common in practical applications. All queries are of the chain type and the spatial predicate used in each individual join operation corresponds to the intersection between two datasets. The queries with three and four datasets have two and six execution plans, respectively. The execution plans for a query are presented in the form $Q{1..8}\;P_{1..6}$, for example, $Q1\;P_1$, $Q1\;P_2$. We argue that this is a representative set of queries as it uses real datasets and involves all types of complex spatial objects and their combinations. However, queries involving more than four datasets should be evaluated in the future, as well as other less common types of queries, such as clique and cycle queries.

\begin{table}[tb]
		\centering
		\caption{Datasets used in experiments.}
		\begin{tabular}{p{4cm}llrr}
				\hline
				\bf Name & \bf Abrev. & \bf Type & \bf Cardinality & \bf SHP File Size (MB)  \\
				%\hline
				%\multicolumn{5}{c}{\it Brazilian datasets (IBGE and LAPIG)} \\
				\hline
				Seaport & $P$ & Points & 120 & <1,0 \\
				Vegetation & $V$ & Polygons & 2.140 & 4,7 \\
				Counties & $M$ & Polygons & 5.564 & 38,8 \\
				Fire alerts & $A$ & Polygons & 32.578 & 11,2 \\
				Roads & $R$ & Lines & 51.646 & 15,2 \\
				Rivers & $H$ & Lines & 226.963 & 64,5 \\
				%\multicolumn{5}{c}{\it Word wide datasets (DCW)} \\
				%\hline
				Crops & $CU$ & Polygons & 123.746 & 69,3 \\
				Rails & $FM$ & Lines & 194.261 & 28,7 \\
				Inland Water & $RA$ & Polygons & 338.860 & 136,7 \\
				Elevation Contour & $CR$ & Lines & 703.574 & 572,5 \\
				Rivers & $HM$ & Lines & 943.638 & 243,2 \\
				\hline
		\end{tabular}
		\label{tab:datasets}
\end{table}



\begin{table}[tb]
	\centering
	\caption{Multiway spatial queries used in experiments.}
	\begin{tabular}{llp{5.8cm}r}
		\hline
		\bf Abrev. & \bf Query & \bf Main characteristic & \bf Join Cardin. \\
		\hline
		Q1 & $FM \Join HM \Join CU$ & Join with polygons and lines & 2313 \\
		Q2 & $A \Join HM \Join FM \Join CU$ & Small set output & 168 \\
		Q3 & $HM \Join FM \Join CR$ & Lines bigger datasets & 2659 \\
		Q4 & $HM \Join CR \Join RA$ & Bigger datasets & 771 \\
		Q5 & $V \Join A \Join M$ & Small datasets, colocated & 36.469 \\
		Q6 & $A \Join P \Join M \Join V$ & Null output set & 0 \\
		Q7 & $HM \Join CR \Join FM \Join R $ & All lines datasets & 3.139 \\
		Q8 & $RA \Join CU \Join A \Join M $ & All polygons datasets & 26.128\\
		\hline
	\end{tabular}
	\label{tab:consultas}
\end{table}

The experiments presented in the next sections were performed on the Azure Platform, using eight A2 virtual machines, with 2 vCPUs and 3.5 GB of RAM each. The machines were allocated at the same data center and were interconnected by a virtual network. To the best of our knowledge, there is no public specification provided by the provider that defines the bandwidth of such a network. However, we experienced a network bandwidth limit of approximately 200Mbps.



\subsection{Trade-off Between Load Balance and Communication}

This section presents an experiment designed to choose a value for the parameter $\id{tradeoff}$ in Algorithm~\ref{alg:communicationcost}. As discussed in Section~\ref{sec:histograma}, the algorithm performs a greedy search to find a schedule that is well-balanced and has a small communication cost. We executed the algorithm for each query plan, varying the $\id{tradeoff}$ from $0.02$ to $0.9$, and measured the total estimated communication cost and the standard deviation of the CPU cost -- the sum of the $IO$ vector and the standard deviation of $Pnts$ vector, respectively.

The results are presented in Figure~\ref{fig:tradeoff}. The $y$ axis indicates the CPU and the Network costs, which have been normalized. Due to space constraints, we present the results of only the four queries: $Q1$, $Q3$, $Q5$, and $Q8$. The other queries produced very similar results. For all execution plans, the point where communication cost matched load balance was approximately at the $\id{tradeoff} = 0.2$. Within this point ($\id{tradeoff} \le 0.2$), load balance continued improving but communication cost grew quickly. Beyond this point, the inverse occurred.

\begin{figure}[b]
\scriptsize
	\centering
	\begin{gnuplot}[terminal=epslatex, terminaloptions={size 15.5cm,3.5cm color colortext}]
		load 'defaultgnuplot.cfg'
		unset key
		set tmargin 0.8
		set bmargin 3.6
		set rmargin 1.0
		unset grid
		set style data lines
		set xtics rotate scale 0,0
		set ylabel "CPU or Network cost"
		set multiplot layout 1,4
		set xlabel "$\\id{tradeoff}$" offset 0,1
		set title '$Q1$' offset 0,-0.5
		plot 'fig/tradeoffq1.txt' u 1:2 lw 3 lt 2, '' u 1:3 lw 3 lt 3, '' u 1:4 lw 3 lt 4, '' u 1:5 lw 3 lt 5 
		set title '$Q3$'
		plot 'fig/tradeoffq3.txt' u 1:2 lw 3 lt 2, '' u 1:3 lw 3 lt 3, '' u 1:4 lw 3 lt 4, '' u 1:5 lw 3 lt 5 
		set title '$Q5$'
		plot 'fig/tradeoffq5.txt' u 1:2 lw 3 lt 2, '' u 1:3 lw 3 lt 3, '' u 1:4 lw 3 lt 4, '' u 1:5 lw 3 lt 5 
		set key box linestyle 12 at screen 0.9,0.1 vertical maxrows 1
		set key samplen 2
		set title '$Q8$'
		plot 'fig/tradeoffq8.txt' u 1:2 lw 3 lt 2 ti 'Balance $P_1$', '' u 1:3 lw 3 lt 3 ti 'Balance $P_2$', '' u 1:4 lw 3 lt 4 ti 'Comm. $P_1$', '' u 1:5 lw 3 lt 5 ti 'Comm. $P_2$'
	\end{gnuplot}
	\caption{Trade-off between communication cost and load balancing in the cluster for queries $Q1$, $Q3$, $Q5$, and $Q8$. The $y$ axis is normalized and corresponds to the CPU or Network cost, conforming to the legend.}
	\label{fig:tradeoff}
\end{figure}

Since $\id{tradeoff}$ values smaller than $0.2$ resulted in an improvement in load balance with a disproportional increase in communication cost, we used $\id{tradeoff} = 0.2$. However, as the graphs show, in systems with better network bandwidth, this value could be significantly reduced to improve query balance.

%This occurs because the metric for balance is an approximation, and using it to force a better estimated balance not always resulted in a better balance in the execution of the query.

 

\subsection{Evaluation of the Proportional Overlap Method} \label{sec:avaliafracionario}

This section presents an experiment designed to evaluate the Proportional Overlap Method, used to compute the cardinality and points metrics of multidimensional histograms, as described in Section~\ref{sec:fractionary}. The experiment consisted of creating a multidimensional histogram for all the datasets involved in the queries of Table~\ref{tab:consultas}, using both the Proportional Overlap Method and the MBR center method \cite{Mamoulis2001a}. Algorithm~\ref{alg:communicationcost} was executed for each execution plan to obtain an estimate of the total number of data (in bytes) to be transferred using each method. The estimates were compared with the real communication cost of each plan, obtained from real executions of the plans in the cluster. For this comparison, we used equation $error = abs(real-estimated) / real * 100$ to compute the distance between the estimated and the real communication costs.

The result of the experiment is presented in Figure~\ref{fig:comparamethist}. The chart shows the error for different execution plans using each of the two methods. The last two groups of bars, in turn, show the average and the maximum error considering all the execution plans. The error incurred with the MBR center method is $\approx 50.7\%$ on average for all plans, and, in the worst case, ($Q6\;P_5$) it was $90.19\%$. The Proportional Overlap method had an average error of $\approx 28.7\%$, and a maximum error of $78.80\%$ (also for $Q6\;P_3$). The error of proportional overlap method is lower in 24 of the 32 queries evaluated (75\%). These values show the effectiveness of the Proportional Overlap method to build better histograms to estimate the communication cost of execution plans. In the following experiments, we used only the Proportional Overlap method, since it provides more precision for the datasets used.


\begin{figure}[tb]
	\centering
	\begin{gnuplot}[terminal=epslatex, terminaloptions={size 15cm,4cm color colortext}]
		load 'defaultgnuplot.cfg'
		set key outside center top horizontal
		set key width -4
		set tmargin 1.0
		set bmargin 2.2
		set rmargin 1.0
		set grid noxtics
		set style fill solid 0.5 border -1
		set style data histogram 
		set style histogram clustered gap 1 title textcolor lt -1
		set xtics rotate scale 0,0
		set ylabel "error (\\%)" offset 2,0
		plot 'fig/histogramerror.txt' using 2:xtic(1) ti 'Proportional O.M.', '' u 3:xtic(1) ti 'MBR Center' 
	\end{gnuplot}
	\caption{Error in the estimation of network communication cost using the two histogram construction methods.}
	\label{fig:comparamethist}
\end{figure}


\subsection{Evaluation of Communication Cost Estimation} \label{sec:communication_evaluation}

In this section we describe an experiment designed to evaluate Algorithm~\ref{alg:communicationcost}, by generating intermediate histograms and computing the overall number of spatial object points and network I/O per server. Algorithm~\ref{alg:communicationcost} is an important part of the plan selection method. It estimates the overall number of points and the communication cost, as well as the load balance in the cluster. The experiment consisted of obtaining an estimate of the communication cost for each plan and server pair, as well as the server with maximum and minimum estimated communication cost. The method then compared the estimated values with real values measured in the current execution of the plans in the cluster. As the estimated communication cost was based on the number of points transferred, and the number of points was obtained from the histograms, a good estimate indicates that the methodology can effectively be used to determine the cost of execution plans for plan selection.

The results of the experiment are shown in Figure~\ref{fig:comparacomun}. To facilitate the comparison of values, the charts are organized in two sections, each with a separated $y$ axis. The section on the right-hand side shows the four plans ($Q3\;P_2$, $Q7\;P_2$, $Q7\;P_3$ and $Q8\;P_2$) with larger $y$ values. Each bar in the chart represents the average of estimated and real network communication cost, among all servers. The total number of communication for each plan can be calculated by multiplying the value of the corresponding bar by 8 (the number of servers used in the experiment). It ranges from $20.13 MB$ ($Q6\;P_6$) to $1.47 GB$ ($Q7\;P_3$). Each bar also has a superior and inferior limit, which represent, respectively, the server with the maximum and minimum number of network communication. A large variation of these limits indicates that the network load is unbalanced in the cluster.

As can be seen in Figure~\ref{fig:comparacomun}, the estimated values are different from the real ones, although they are very close. This behavior also holds for the error bars, which compare the estimated maximum and minimum errors with the maximum and minimum real values. This confirms that the algorithm provides a good estimate of the total network communication cost, as well as the network communication cost for each server. There is a set of plans that shows a distinct behavior: $Q5\;P_{1..2}$, $Q6\;P_{2..4}$, and $Q8\;P_{4..6}$. While, for all other plans, the estimated value is below the real value, these plans show the opposite. This behavior is due to the propagation of the error in the estimation of the size of the objects (see line \ref{alg_comm_next_points} of Algorithm~\ref{alg:communicationcost}). Because an average of the object sizes is propagated to the intermediate histogram, big objects that could be filtered in a real execution of the join remain in the estimation, causing this error.

The results reflect the fact that histograms are a simplified representation of the datasets. The method cannot produce precise estimates due to this simplification. However, the data show reasonable proximity between the estimates and the real values for all plans used in the experiment. We concluded that Algorithm~\ref{alg:communicationcost} and the associated methodology can provide a good estimate of the communication cost for each execution plan based on the distribution of the data and on the multidimensional histograms. 

\begin{figure}[tb]
	\centering
	\begin{gnuplot}[terminal=epslatex, terminaloptions={size 15cm,4cm color colortext}]
		load 'defaultgnuplot.cfg'
		set key top center
		set key width -4
		set tmargin 0.1
		set lmargin 3.5
		set bmargin 2.3
		set grid noxtics
		set style fill solid 0.5 border -1
		set style data histogram 
		set style histogram errorbar gap 1 title textcolor lt -1
		set boxwidth 0.7 relative
		set xtics rotate scale 0,0
		set ylabel "Mbytes" offset 2.2,0

		set multiplot layout 1,2
		
		set rmargin at screen 0.80
		plot 'fig/comparacomunica2.txt' using ($2/1024):($4/1024):($3/1024):xtic(1) ti 'Estimated', '' u ($5/1024):($6/1024):($7/1024):xtic(1) dt 1 ti 'Real' 

		unset key
		unset ytics
		set y2tics
        set grid y2tics
		set rmargin 4
		set lmargin at screen 0.81
		set ylabel "Mbytes" offset 13.6,0
		set border 9 back ls 11

		plot 'fig/comparacomunica.txt' using ($2/1024):($4/1024):($3/1024):xtic(1) ti '', '' u ($5/1024):($6/1024):($7/1024):xtic(1) dt 1

	\end{gnuplot}
	\caption{Average server communication: estimated and real for each query. The error bars represent the minimum and maximum network communication between all servers.}
	\label{fig:comparacomun}
\end{figure}


\subsection{Evaluation of the Selection of Distributed Execution Plans} \label{sec:plan_chosen}

In this section we describe an evaluation of the plan selection methodology described in Section~\ref{sec:escolhaplanodst}. The experiment consisted of selecting an execution plan using the estimates provided by Algorithm~\ref{alg:communicationcost} and~(\ref{eqn:optimal}). It then measured the actual execution time of all plans in the cluster to check whether or not the selected plan (based on the estimates) was indeed the plan with the shortest execution time.

Figure~\ref{fig:escolheplano} shows the estimated cost of all query plans, based on the number of points of the spatial objects. An arrow indicates the plan that was selected by the method -- the lowest bar in each group. Each bar in the chart indicates the estimated number of points to be compared (in millions) by the server with the highest number of points ($Max(AllPlans[p])$ according to~(\ref{eqn:optimal})).

The execution time for each of the $32$ plans is presented in Figure~\ref{fig:temporeal}. The chart uses a logarithmic scale due to the difference between the execution times of the best plan (in seconds) and the worst plan (in hundreds of seconds). The chart shows that the selected plan was indeed the best one for all the queries in the experiment. Furthermore, the method was able to establish a consistent relative ordering of all plans for almost all queries. The only two differences in the ordering were $Q8\;P_5$ and $Q8\;P_6$, which appeared in inverse order between estimate and real execution ($0.8$ seconds of difference between the plans), and $Q6\;P_{2/5, 4/3}$, which also differed from the estimates, but only by $0.2$ seconds on average. The error in the ordering could lead to the wrong choice of an execution. However, as can be observed in the chart, the error always occurred with very similar plans with respect to their execution times. A wrong choice in this scenario will still lead to a good plan, with a slightly worse execution time. This confirms that the overall number of points of the spatial objects is a good metric to determine the cost of execution plans. It also confirms that~(\ref{eqn:optimal}) selects a relatively good plan, based on the estimated cost of plans.

\begin{figure}[tb]
	\centering
	\begin{gnuplot}[terminal=epslatex, terminaloptions={size 15cm,3cm color colortext}]
		load 'defaultgnuplot.cfg'
		set key outside spacing 0.8
		set yrange [0:270]
		set rmargin 6
		set bmargin 0.5
		set tmargin 0
		set grid
		set style fill solid 0.5 border -1
		set style data histograms
		set style histogram clustered gap 2 title textcolor lt -1
		set boxwidth 0.8 relative
		set ylabel "$Max(AllPlans[p])/10^6$" offset 1.7,0
		unset xlabel 

		set arrow from screen .131,.5 to screen .131,.2 ls 1 lw 3
		set arrow from screen .225,.5 to screen .225,.2 ls 1 lw 3
		set arrow from screen .319,.5 to screen .319,.2 ls 1 lw 3
		set arrow from screen .425,.5 to screen .425,.38 ls 1 lw 3
		set arrow from screen .555,.5 to screen .555,.135 ls 1 lw 3
		set arrow from screen .662,.5 to screen .662,.12 ls 1 lw 3
		set arrow from screen .743,.5 to screen .743,.145 ls 1 lw 3
		set arrow from screen .825,.5 to screen .825,.18 ls 1 lw 3

		plot 'fig/estimams.txt' using ($2/1000000):xtic(1) ti col, '' u ($3/1000000) ti col, '' u ($4/1000000) ti col, '' u ($5/1000000) ti col, '' u ($6/1000000) ti col, '' u ($7/1000000) ti col
	\end{gnuplot}
	\caption{Estimated maximum number of spatial points in a server for each plan, after query scheduling by Algorithm~\ref{alg:communicationcost}. Red arrows indicate the best plan selected for each query by~(\ref{eqn:optimal}).}
	\label{fig:escolheplano}
\end{figure}

\begin{figure}[tb]
	\centering
	\begin{gnuplot}[terminal=epslatex, terminaloptions={size 15cm,2.5cm color colortext}]
		load 'defaultgnuplot.cfg'
		set key outside spacing 0.8
		set logscale y
		set rmargin 6
		set bmargin 0.5
		set tmargin 0
		set grid
		set style fill solid 0.5 border -1
		set style data histograms
		set style histogram clustered gap 2 title textcolor lt -1
		set boxwidth 0.8 relative
		set ylabel "$sec$" offset 4,0
		unset xlabel 

		plot 'fig/realms.txt' using ($2/1000):xtic(1) ti col, '' u ($3/1000) ti col, '' u ($4/1000) ti col, '' u ($5/1000) ti col, '' u ($6/1000) ti col, '' u ($7/1000) ti col
	\end{gnuplot}
	\caption{Execution time of all plans in the cluster.}
	\label{fig:temporeal}
\end{figure}



\section{Conclusion}\label{sec:conclusao}

In this article we proposed a comprehensive solution methodology for multiway spatial join query processing in distributed systems. We identified metadata related to the datasets and the data distribution, which are relevant for the distributed processing of multiway spatial join queries. We proposed a novel use of multidimensional histogram as a data structure to organize such metadata, and which is also used as the basis for a distributed data access method in computing clusters. A new method for histogram construction was proposed based on the proportional overlap of spatial objects and cell boundaries. The proposed method performs better than a well-know method, and provides a better estimate for the cost of execution plans in 75\% of the queries evaluated. The complete methodology for cost estimation was used to implement an optimizer that was able to choose the best execution plan for all the queries considered in the experiments (which use real spatial datasets). Furthermore, the methodology also provided a good ordering of query plans (only two errors, but producing very similar plans), which indicates that it can be used to select good plans for more complex queries.

We also proposed an improved version of Clone Join. The new algorithm was adapted with the reference point method to prevent duplicated results from being reported by the join operation. The modified version of CJ is also a non-blocking algorithm, meaning that it can be used to efficiently process multiway spatial join queries in distributed systems.

We plan to investigate new metrics to be included in the histogram, such as the current workload that results from the execution of concurrent multiway spatial join queries. Besides having the potential to produce more precise estimates, this may lead to a more realistic model for production-grade spatial databases and to address interesting challenges for the scheduling of distributed queries in a cluster. We also plan to replace the greedy algorithm by more elaborate combinatorial methods in order to achieve better query scheduling in the cluster.

\begin{received} 
\end{received}

\bibliographystyle{jidm}
\bibliography{jidmb}

\end{document}

