From 81755ea6c08b0cc72a8548b886ae5a009fdb647a Mon Sep 17 00:00:00 2001 From: Garry Boyer Date: Thu, 6 Sep 2007 17:15:26 +0000 Subject: [PATCH] hi --- fastlib/base/cc.cc | 7 +- fastlib/base/cc.h | 6 + fastlib/base/otrav.h | 5 +- fastlib/base/otrav_impl.h | 1 + fastlib/col/string.h | 7 + fastlib/math/kernel.h | 9 +- fastlib/thor/dfs.h | 8 +- fastlib/thor/dfs_impl.h | 32 +-- fastlib/thor/gnp.h | 67 ++++++ fastlib/thor/thor_utils.h | 8 + fastlib/thor/thor_utils_impl.h | 10 +- fastlib/thor/thortree.h | 12 +- fastlib/thor/thortree_impl.h | 5 +- fastlib/tree/bounds.h | 76 +++++-- fastlib/u/garryb/nbr/build.py | 4 + fastlib/u/garryb/nbr/paper/paper.tex | 1 + fastlib/u/garryb/thor.tex | 299 +++++++++++++++++++++++---- fastlib/u/rriegel/nbc/nbc.cc | 11 +- 18 files changed, 462 insertions(+), 106 deletions(-) diff --git a/fastlib/base/cc.cc b/fastlib/base/cc.cc index b470ace495..86a6931997 100644 --- a/fastlib/base/cc.cc +++ b/fastlib/base/cc.cc @@ -6,9 +6,14 @@ #include "cc.h" -static const double DBL_ZERO = 0.0; +namespace { + const double DBL_ZERO = 0.0; +}; + const double DBL_NAN = DBL_ZERO / DBL_ZERO; const double FLT_NAN = DBL_NAN; +const double DBL_INF = 1.0 / DBL_ZERO; +const double FLT_INF = DBL_INF; #if defined(DEBUG) || defined(PROFILE) namespace cc__private { diff --git a/fastlib/base/cc.h b/fastlib/base/cc.h index 7c7e9214b6..805e0b37ba 100644 --- a/fastlib/base/cc.h +++ b/fastlib/base/cc.h @@ -22,6 +22,12 @@ extern const double DBL_NAN; /** NaN value for floats. Use isnanf to check for this. */ extern const double FLT_NAN; +/** Infinity value for doubles. */ +extern const double DBL_INF; + +/** Infinity value for floats. */ +extern const double FLT_INF; + /* TODO: Decide if templated min and max would prefer const refs. * * I'm purposely not making these const refs for the case of integers diff --git a/fastlib/base/otrav.h b/fastlib/base/otrav.h index 626fb2c733..35122ef93e 100644 --- a/fastlib/base/otrav.h +++ b/fastlib/base/otrav.h @@ -73,7 +73,10 @@ * Note that OT_DEF calls this automatically. */ #define OT_GEN_COPY_CONSTRUCTOR(AClass) \ - public: AClass(const AClass& other) { ot__private::ZOTDeepCopier::Doit(other, this); } private: + public: AClass(const AClass& other) { \ + ot__private::ZOTDeepCopier::Doit(other, this); \ + } \ + private: /** * Automatically create a dstructor based on object traversal. diff --git a/fastlib/base/otrav_impl.h b/fastlib/base/otrav_impl.h index cf8130c1ea..2505f62a34 100644 --- a/fastlib/base/otrav_impl.h +++ b/fastlib/base/otrav_impl.h @@ -507,6 +507,7 @@ namespace ot__private { static void Doit(const T& src, T *dest) { ot__private::ZOTDeepCopier d; mem::Copy(dest, &src, 1); + TraverseObjectPostprocess(dest); d.MyObject(*dest); } diff --git a/fastlib/col/string.h b/fastlib/col/string.h index 012dd7a8e1..d807b64ba2 100644 --- a/fastlib/col/string.h +++ b/fastlib/col/string.h @@ -389,6 +389,13 @@ class String { return CompareNoCase(s) == 0; } + /** + * Checks if this string begins with another string. + */ + bool StartsWith(const char *s) const { + return strncmp(array_.begin(), s, strlen(s)) == 0; + } + /** * Compares two strings. * diff --git a/fastlib/math/kernel.h b/fastlib/math/kernel.h index c807b034bd..c75a264049 100644 --- a/fastlib/math/kernel.h +++ b/fastlib/math/kernel.h @@ -24,22 +24,29 @@ struct GaussianKernel { private: double inv_bandwidth_2sq_; + double bandwidth_sq_; OT_DEF_BASIC(GaussianKernel) { OT_MY_OBJECT(inv_bandwidth_2sq_); + OT_MY_OBJECT(bandwidth_sq_); } public: static const bool HAS_CUTOFF = false; public: + double bandwidth_sq() const { + return bandwidth_sq_; + } + /** * Initializes to a specific bandwidth. * * @param bandwidth_in the standard deviation sigma */ void Init(double bandwidth_in) { - inv_bandwidth_2sq_ = 1.0 / (2.0 * bandwidth_in * bandwidth_in); + bandwidth_sq_ = bandwidth_in * bandwidth_in; + inv_bandwidth_2sq_ = 1.0 / (2.0 * bandwidth_sq_); } /** diff --git a/fastlib/thor/dfs.h b/fastlib/thor/dfs.h index ae0da96352..7a84d8310b 100644 --- a/fastlib/thor/dfs.h +++ b/fastlib/thor/dfs.h @@ -42,9 +42,7 @@ class DualTreeDepthFirst { const typename GNP::RNode *r_root_; bool do_naive_; - uint64 n_naive_; - uint64 n_pre_naive_; - uint64 n_recurse_; + DualTreeRecursionStats stats_; public: DualTreeDepthFirst() {} @@ -74,6 +72,10 @@ class DualTreeDepthFirst { const typename GNP::GlobalResult& global_result() const { return global_result_; } + + const DualTreeRecursionStats& stats() const { + return stats_; + } private: COMPILER_NOINLINE diff --git a/fastlib/thor/dfs_impl.h b/fastlib/thor/dfs_impl.h index 09f1cd5544..296de5a186 100644 --- a/fastlib/thor/dfs_impl.h +++ b/fastlib/thor/dfs_impl.h @@ -52,14 +52,14 @@ void DualTreeDepthFirst::Begin_(index_t q_root_index) { CacheRead q_root(&q_nodes_, q_root_index); QMutables *q_root_mut = &q_mutables_[q_root_index]; - DEBUG_ONLY(n_naive_ = 0); - DEBUG_ONLY(n_pre_naive_ = 0); - DEBUG_ONLY(n_recurse_ = 0); - bool need_explore = GNP::Algorithm::ConsiderPairIntrinsic( param_, *q_root, *r_root_, &delta, &global_result_, &q_root_mut->postponed); + stats_.Init(); + stats_.tuples_analyzed = q_root->count() * r_root_->count(); + stats_.n_queries = q_root->count(); + if (need_explore) { typename GNP::QSummaryResult empty_summary_result; @@ -73,24 +73,6 @@ void DualTreeDepthFirst::Begin_(index_t q_root_index) { } PushDownPostprocess_(q_root_index, q_root_mut); - - //fx_timer_stop(datanode_, "execute"); - /*DEBUG_ONLY(fx_format_result(datanode_, "naive_ratio", "%f", - 1.0 * n_naive_ / q_root->count() / r_root_->count())); - DEBUG_ONLY(fx_format_result(datanode_, "naive_per_query", "%f", - 1.0 * n_naive_ / q_root->count())); - DEBUG_ONLY(fx_format_result(datanode_, "pre_naive_ratio", "%f", - 1.0 * n_pre_naive_ / q_root->count() / r_root_->count())); - DEBUG_ONLY(fx_format_result(datanode_, "pre_naive_per_query", "%f", - 1.0 * n_pre_naive_ / q_root->count())); - DEBUG_ONLY(fx_format_result(datanode_, "recurse_ratio", "%f", - 1.0 * n_recurse_ / q_root->count() / r_root_->count())); - DEBUG_ONLY(fx_format_result(datanode_, "recurse_per_query", "%f", - 1.0 * n_recurse_ / q_root->count()));*/ - -/* if (fx_param_bool(datanode_, "print", 0)) { - ot::Print(q_results_); - }*/ } template @@ -131,7 +113,7 @@ void DualTreeDepthFirst::Pair_( DEBUG_MSG(1.0, "Checking (%d,%d) x (%d,%d)", q_node->begin(), q_node->end(), r_node->begin(), r_node->end()); - DEBUG_ONLY(n_recurse_++); + DEBUG_ONLY(stats_.node_node_considered++); /* begin prune checks */ typename GNP::QSummaryResult mu(q_node_mut->summary_result); @@ -241,7 +223,7 @@ void DualTreeDepthFirst::BaseCase_( const typename GNP::QSummaryResult& unvisited, QMutables *q_node_mut) { - DEBUG_ONLY(n_pre_naive_ += q_node->count() * r_node->count()); + DEBUG_ONLY(stats_.node_point_considered += q_node->count()); q_node_mut->summary_result.StartReaccumulate(param_, *q_node); @@ -279,7 +261,7 @@ void DualTreeDepthFirst::BaseCase_( visitor.FinishVisitingQueryPoint(param_, *q_point, q_i, *r_node, unvisited, q_result, &global_result_); - DEBUG_ONLY(n_naive_ += r_node->count()); + DEBUG_ONLY(stats_.point_point_considered += r_node->count()); } q_node_mut->summary_result.Accumulate(param_, *q_result); diff --git a/fastlib/thor/gnp.h b/fastlib/thor/gnp.h index 317fd9433c..860d2234f3 100644 --- a/fastlib/thor/gnp.h +++ b/fastlib/thor/gnp.h @@ -217,4 +217,71 @@ class BlankAlgorithm { } }; +struct DualTreeRecursionStats { + public: + /** + * An rpc-style reductor suitable for recursion stats. + */ + struct Reductor { + /** + * Reduces two elements. + * + * @param right a new element to merge + * @param left the element to merge into + */ + void Reduce(const DualTreeRecursionStats& right, + DualTreeRecursionStats* left) const { + left->Add(right); + } + }; + + public: + double tuples_analyzed; + index_t n_queries; + int64 node_node_considered; + int64 node_point_considered; + int64 point_point_considered; + + OT_DEF(DualTreeRecursionStats) { + OT_MY_OBJECT(tuples_analyzed); + OT_MY_OBJECT(n_queries); + OT_MY_OBJECT(node_node_considered); + OT_MY_OBJECT(node_point_considered); + OT_MY_OBJECT(point_point_considered); + } + + public: + void Init() { + tuples_analyzed = 0; + n_queries = 0; + node_node_considered = 0; + node_point_considered = 0; + point_point_considered = 0; + } + + void Add(const DualTreeRecursionStats& other) { + tuples_analyzed += other.tuples_analyzed; + n_queries += other.n_queries; + node_node_considered += other.node_node_considered; + node_point_considered += other.node_point_considered; + point_point_considered += other.point_point_considered; + } + + void Report(datanode *module) { + fx_format_result(module, "p_node_node", + "%g", node_node_considered / tuples_analyzed); + fx_format_result(module, "p_node_point", + "%g", node_point_considered / tuples_analyzed); + fx_format_result(module, "p_point_point", + "%g", point_point_considered / tuples_analyzed); + + fx_format_result(module, "r_node_node", + "%g", 1.0 * node_node_considered / n_queries); + fx_format_result(module, "r_node_point", + "%g", 1.0 * node_point_considered / n_queries); + fx_format_result(module, "r_point_point", + "%g", 1.0 * point_point_considered / n_queries); + } +}; + #endif diff --git a/fastlib/thor/thor_utils.h b/fastlib/thor/thor_utils.h index 1abd7b6ce9..63202f44b8 100644 --- a/fastlib/thor/thor_utils.h +++ b/fastlib/thor/thor_utils.h @@ -50,6 +50,7 @@ class ThreadedDualTreeSolver { DistributedCache *r_nodes_cache_; DistributedCache *q_results_cache_; typename GNP::GlobalResult global_result_; + DualTreeRecursionStats stats_; Mutex mutex_; public: @@ -80,6 +81,13 @@ class ThreadedDualTreeSolver { return global_result_; } + /** + * Recursion statistics. + */ + const DualTreeRecursionStats& stats() const { + return stats_; + } + private: void ThreadBody_(); }; diff --git a/fastlib/thor/thor_utils_impl.h b/fastlib/thor/thor_utils_impl.h index 94cffdf3e3..8b6eb014a7 100644 --- a/fastlib/thor/thor_utils_impl.h +++ b/fastlib/thor/thor_utils_impl.h @@ -18,6 +18,7 @@ void thor::ThreadedDualTreeSolver::Doit( q_results_cache_ = q_results_cache_in; global_result_.Init(*param_); + stats_.Init(); if (n_threads > 1) { ArrayList threads; @@ -61,6 +62,7 @@ void thor::ThreadedDualTreeSolver::ThreadBody_() { q_results_cache_); mutex_.Lock(); + stats_.Add(solver.stats()); global_result_.Accumulate(*param_, solver.global_result()); mutex_.Unlock(); } @@ -202,10 +204,16 @@ void thor::RpcDualTree(datanode *module, int base_channel, q_results->WaitSync(fx_submodule(io_module, NULL, "q_results")); fx_timer_stop(module, "write_results"); +#ifdef DEBUG + DualTreeRecursionStats stats = solver.stats(); + rpc::Reduce(base_channel + 5, DualTreeRecursionStats::Reductor(), &stats); + stats.Report(fx_submodule(module, NULL, "recursion")); +#endif + GlobalResultReductor global_result_reductor; typename GNP::GlobalResult my_global_result = solver.global_result(); global_result_reductor.Init(¶m); - rpc::Reduce(base_channel + 5, global_result_reductor, &my_global_result); + rpc::Reduce(base_channel + 6, global_result_reductor, &my_global_result); work_queue->Report(work_module); my_global_result.Report(param, diff --git a/fastlib/thor/thortree.h b/fastlib/thor/thortree.h index 1180dac5ed..1af4054efb 100644 --- a/fastlib/thor/thortree.h +++ b/fastlib/thor/thortree.h @@ -24,7 +24,15 @@ struct TreeGrain { /** One past the last point. */ index_t point_end_index; - void InitBlank() { + template + void Init(const TSkeletonNode& node) { + node_index = node.index(); + node_end_index = node.end_index(); + point_begin_index = node.node().begin(); + point_end_index = node.node().end(); + } + + void InitInvalid() { node_index = -1; node_end_index = -1; point_begin_index = -1; @@ -116,7 +124,7 @@ class ThorTreeDecomposition { DEBUG_ASSERT(root_in != NULL); grain_by_owner_.Init(rpc::n_peers()); for (int i = 0; i < grain_by_owner_.size(); i++) { - grain_by_owner_[i].InitBlank(); + grain_by_owner_[i].InitInvalid(); } FillLinearization_(root_); } diff --git a/fastlib/thor/thortree_impl.h b/fastlib/thor/thortree_impl.h index facb63c7f5..1a5473ed2f 100644 --- a/fastlib/thor/thortree_impl.h +++ b/fastlib/thor/thortree_impl.h @@ -71,10 +71,7 @@ void ThorTreeDecomposition::FillLinearization_(DecompNode *node) { if (node->info().is_singleton() || !node->is_complete()) { grain = &grain_by_owner_[node->info().begin_rank]; - grain->node_index = node->index(); - grain->node_end_index = node->end_index(); - grain->point_begin_index = node->node().begin(); - grain->point_end_index = node->node().end(); + grain->Init(*node); } else { for (int k = 0; k < Node::CARDINALITY; k++) { FillLinearization_(node->child(k)); diff --git a/fastlib/tree/bounds.h b/fastlib/tree/bounds.h index 3f305d1aa7..6dd718e7e7 100644 --- a/fastlib/tree/bounds.h +++ b/fastlib/tree/bounds.h @@ -368,7 +368,7 @@ class DHrectBound { * to the specified power. */ double MinDistanceSq(const double *mpoint) const { - double sumsq = 0; + double sum = 0; const DRange *mbound = bounds_; index_t d = dim_; @@ -383,10 +383,10 @@ class DHrectBound { mbound++; mpoint++; - sumsq += math::Pow(v); + sum += math::Pow(v); } while (--d); - return math::Pow<2, t_pow>(sumsq) / 4; + return math::Pow<2, t_pow>(sum) / 4; } /** @@ -410,7 +410,7 @@ class DHrectBound { * */ double MinToMidSq(const DHrectBound& other) const { - double sumsq = 0; + double sum = 0; const DRange *a = this->bounds_; const DRange *b = other.bounds_; @@ -426,10 +426,10 @@ class DHrectBound { a++; b++; - sumsq += math::Pow(v); + sum += math::Pow(v); } - return math::Pow<2, t_pow>(sumsq) / 4; + return math::Pow<2, t_pow>(sum) / 4; } /** @@ -437,26 +437,26 @@ class DHrectBound { * to the specified power. */ double MaxDistanceSq(const Vector& point) const { - double sumsq = 0; + double sum = 0; DEBUG_ASSERT(point.length() == dim_); for (index_t d = 0; d < dim_; d++) { - sumsq += math::Pow( + sum += math::Pow( max(point[d] - bounds_[d].lo, bounds_[d].hi - point[d])); } - return math::Pow<2, t_pow>(sumsq); + return math::Pow<2, t_pow>(sum); } /** - * Calculates minimum bound-to-point squared distance, + * Calculates minimum bound-to-bound squared distance, * to the specified power. * * Example: bound1.MinDistanceSq(other) for minimum squared distance. */ double MinDistanceSq(const DHrectBound& other) const { - double sumsq = 0; + double sum = 0; const DRange *a = this->bounds_; const DRange *b = other.bounds_; index_t mdim = dim_; @@ -471,10 +471,42 @@ class DHrectBound { // (x * 2)^2 / 4 = x^2 double v = (v1 + fabs(v1)) + (v2 + fabs(v2)); - sumsq += math::Pow(v); + sum += math::Pow(v); } - return math::Pow<2, t_pow>(sumsq) / 4; + return math::Pow<2, t_pow>(sum) / 4; + } + + /** + * Calculates minimum and maximum bound-to-bound squared distance, + * to the specified power. + * + * Example: bound1.MinDistanceSq(other) for minimum squared distance. + */ + DRange RangeDistanceSq(const DHrectBound& other) const { + double sum_lo = 0; + double sum_hi = 0; + const DRange *a = this->bounds_; + const DRange *b = other.bounds_; + index_t mdim = dim_; + + DEBUG_SAME_INT(dim_, other.dim_); + + for (index_t d = 0; d < mdim; d++) { + double v1 = b[d].lo - a[d].hi; + double v2 = a[d].lo - b[d].hi; + // We invoke the following: + // x + fabs(x) = max(x * 2, 0) + // (x * 2)^2 / 4 = x^2 + double v_lo = (v1 + fabs(v1)) + (v2 + fabs(v2)); + double v_hi = min(v1, v2); + + sum_lo += math::Pow(v_lo); + sum_hi += math::Pow(v_hi); + } + + return DRange(math::Pow<2, t_pow>(sum_lo) / 4, + math::Pow<2, t_pow>(sum_hi)); } /** @@ -482,7 +514,7 @@ class DHrectBound { * to the specified power. */ double MinimaxDistanceSq(const DHrectBound& other) const { - double sumsq = 0; + double sum = 0; const DRange *a = this->bounds_; const DRange *b = other.bounds_; index_t mdim = dim_; @@ -494,10 +526,10 @@ class DHrectBound { double v2 = a[d].lo - b[d].lo; double v = max(v1, v2); v = (v + fabs(v)); /* truncate negatives to zero */ - sumsq += math::Pow(v); + sum += math::Pow(v); } - return math::Pow<2, t_pow>(sumsq) / 4; + return math::Pow<2, t_pow>(sum) / 4; } /** @@ -505,18 +537,18 @@ class DHrectBound { * to the specified power. */ double MaxDistanceSq(const DHrectBound& other) const { - double sumsq = 0; + double sum = 0; const DRange *a = this->bounds_; const DRange *b = other.bounds_; DEBUG_ASSERT(dim_ == other.dim_); for (index_t d = 0; d < dim_; d++) { - sumsq += math::Pow( + sum += math::Pow( max(b[d].hi - a[d].lo, a[d].hi - b[d].lo)); } - return math::Pow<2, t_pow>(sumsq); + return math::Pow<2, t_pow>(sum); } /** @@ -524,17 +556,17 @@ class DHrectBound { * to the specified power. */ double MidDistanceSq(const DHrectBound& other) const { - double sumsq = 0; + double sum = 0; const DRange *a = this->bounds_; const DRange *b = other.bounds_; DEBUG_ASSERT(dim_ == other.dim_); for (index_t d = 0; d < dim_; d++) { - sumsq += math::PowAbs(a[d].hi + a[d].lo - b[d].hi - b[d].lo); + sum += math::PowAbs(a[d].hi + a[d].lo - b[d].hi - b[d].lo); } - return math::Pow<2, t_pow>(sumsq) / 4; + return math::Pow<2, t_pow>(sum) / 4; } /** diff --git a/fastlib/u/garryb/nbr/build.py b/fastlib/u/garryb/nbr/build.py index 1d503c4667..5dc618b974 100644 --- a/fastlib/u/garryb/nbr/build.py +++ b/fastlib/u/garryb/nbr/build.py @@ -3,6 +3,10 @@ binrule(name = "tkde", sources = ["tkde.cc"], deplibs = ["thor:thor"]) +binrule(name = "fdkde", + sources = ["fdkde.cc"], + deplibs = ["thor:thor"]) + binrule(name = "allnn_rpc", headers = ["allnn.cc"], sources = ["allnn_rpc.cc"], diff --git a/fastlib/u/garryb/nbr/paper/paper.tex b/fastlib/u/garryb/nbr/paper/paper.tex index 1762e833a4..0e6750e582 100644 --- a/fastlib/u/garryb/nbr/paper/paper.tex +++ b/fastlib/u/garryb/nbr/paper/paper.tex @@ -233,6 +233,7 @@ Parallel programming is well-known for its significant cost in terms of developm However, we show that all GNP's can be phrased as an instance of a generalized $N$-body algorithm, leading to an intriguing solution: Invest in a parallel solver for the general problem, requiring algorithm developers to provide only problem-specific portions. In practice, this approach leads to efficient parallel and serial performance, and furthermore, the parallel algorithm in some cases takes less time to develop than a from-scratch serial algorithm. We describe in this paper both the formal problem generalization, and our first-version implementation \THOR that handles dual-tree algorithms with a locally-ordered depth-first expansion pattern. +\authorsnote{Make sure it's clear at this point what is novel in my paper. What should be clear is: (a) the mathematical model, and (b) its implementation in parallel.} \section{Fast Dual-Tree Algorithms} diff --git a/fastlib/u/garryb/thor.tex b/fastlib/u/garryb/thor.tex index 595a4e38a9..08f8359918 100644 --- a/fastlib/u/garryb/thor.tex +++ b/fastlib/u/garryb/thor.tex @@ -11,12 +11,23 @@ \newcommand{\itemt}[1]{\item {\bf #1} -} + +\newcommand{\authornote}[1]{\footnote{Note to self: #1}} +\newcommand{\authorsnote}[1]{\authornote{#1}} +\newcommand{\com}[1]{{\small \textit{((#1))}}} + \newcommand{\union}{\cup} \newcommand{\intersect}{\cap} \newcommand{\Union}{\bigcup} \newcommand{\Intersect}{\bigcap} \newcommand{\bigvec}[1]{\mathop{\overrightarrow{#1}}} +\newcommand{\otimeshat}{\widehat{\otimes}} +\newcommand{\odothat}{\widehat{\odot}} + +\newcommand{\prefsplit}[2]{#1 \succ #2} +\newcommand{\summary}{\hat{\sigma}} + \DeclareMathOperator*{\map}{map} \DeclareMathOperator*{\worst}{worst} \DeclareMathOperator*{\argmin}{argmin} @@ -274,18 +285,25 @@ \end{itemize} \end{slide} +\begin{slide}{GNP Basics} + Unlike many software frameworks, ours is inherently mathematical. + + Next, we dive into the formalisms behind our framework. +\end{slide} + \begin{slide}{GNP Basics - Notation} \begin{itemize} \itemt{Map Operator} $\map$ builds a list of results. Mathematically, $\map$ outputs a set of key-value pairs. \itemt{Generic Operators} $\bigodot, \bigotimes$ mean {\em any} commutative, associative operator. - \itemt{Trees} A tree node is {\em equivalent} to a set of points. - $\kdroot{X}$ is the entire data set, $X$ is a subset, - $X = \kdleft{X} \union \kdright{X}$ partitions a node into two children. - \itemt{Sets} $X, Y$ refer to two sets. $Q, R$ are used for queries and - references, when the $\map$ operator is in use. + \itemt{Trees} A tree node {\em summarizes} a set of points. + $\kdroot{X}$ is the entire data set. + $X$ is a subset with parent $\kdparent{X}$, + partitioned into children $X = \kdleft{X} \union \kdright{X}$. \itemt{Generic Statistics} $\outstat(Q)$ is a statistic on a set of points, such as bounding box, mean, or variance. + \itemt{Upper/Lower Bounds} $\distup(\outstat(X), \outstat(Y))$ is an upper-bound distance + between bounding boxes. $\distlo$ is the lower-bound. \end{itemize} \end{slide} @@ -305,7 +323,6 @@ \end{array}\] \end{slide} - \begin{slide}{GNP Basics - Expansion and Pruning} \begin{itemize} \itemt{Expansion} @@ -330,65 +347,273 @@ \end{itemize} \end{slide} -\begin{slide}{GNP Basics - Q/R Problems} +\begin{slide}{Q/R Problems} \begin{itemize} - \item A \defterm{query-reference} problem has $\bigodot = \map$, i.e. one result per query. - (Consider $\map$'s output ``under the hood'' to be a set of key-value pairs.) - \item Any second-order reduce problem can be transformed into a query-reference problem: - Replace outer operator with $\map$, and apply the original operator as a post-processing step. - \item We'll only talk about query-reference problems for the rest of the slides. + \itemt{Definition} A \defterm{query-reference} problem has $\bigodot = \map$, i.e. one result per query. + \itemt{Universality} {\it Any} second-order reduce problem can be transformed into a query-reference problem. + \begin{itemize} \item + Replace outer operator with $\map$, and treat original operator as a postprocess step. + \end{itemize} + \end{itemize} + Next, we'll focus on solving these problems formally. +\end{slide} + +\begin{slide}{Q/R Problems - Formalization} + Assume commutative, associative $\opqr$. + For each query $q$, solve $\outqr(q, \kdroot{R})$, given: + \begin{eqnarray*} + \outqr(q, R) &\equiv& \gqr(q, \inqr(q, R)), +x \\ + \inqr(q, R) &\equiv& \Opqr_{r \in R} \fqr(q, r). + \end{eqnarray*} + Sometimes, the same \defterm{mass result} $\inqrv(Q, R)$ + applies to many queries for a particular reference node. + This is a prune, formally expressed, + \[ + \forall q \in Q,~~ \inqr(q, R) \gets \inqrv(Q, R). + \] +\end{slide} + +\begin{slide}{Q/R Problems - Solving} + We may compute $\inqrv$ using the three rules: + \[\begin{array}{l} + \inqrv(Q, R) = \inqrv(Q, \kdleft{R}) \opqr \inqrv(Q, \kdright{R}), + \\ + \text{if prune occurs for } \kdparent{Q} \supset Q \text{, then } \inqrv(Q, R) = \inqrv(\kdparent{Q}, R), + \\ + \text{if } \canpruneqrv(\outstat(Q), \outstat(R)) \text{, then } \inqrv(Q, R) = \deltaqrv(\outstat(Q), \outstat(R)). + \end{array}\] + We have: + \begin{itemize} + \itemt{$\canpruneqrv$} An intrinsic prune check. + \itemt{$\deltaqrv$} The pruned value, computed with summary statistics. \end{itemize} \end{slide} -\begin{slide}{Math - QR Intrinsic} +\begin{slide}{Q/R Problems - Example} + Range count is a query-reference analog of two-point correlation, + \[\map_{q \in Q} \sum_{r \in R} I(\dist{q}{r} < h).\] + A THOR user simply implements the following in C++. + \begin{eqnarray*} + \gqr(q, \inqr) &\equiv& \inqr + \\ + \opqr &\equiv& + + \\ + \fqr(q,r) &\equiv& I(\dist{q}{r} < h) + \\ + \canpruneqrv(\sigma(Q), \sigma(R)) + &\equiv& + \begin{array}{l}\distup(\outstat(Q),\outstat(R)) < h \\ \vee \distlo(\outstat(Q),\outstat(R)) \geq h\end{array} + \\ + \deltaqrv(\outstat(Q),\outstat(R)) &\equiv& \left\{ \begin{array}{l} 0 \text{ if } \distup(\outstat(Q),\outstat(R)) < h \\ |R| \text{ if } \distlo(\outstat(Q),\outstat(R)) \geq h \end{array}\right. + \end{eqnarray*} +\end{slide} + +\begin{slide}{Q/R Problems - Extrinsic Pruning} + \begin{itemize} + \item We compute \defterm{query summary results} $\inmu(Q, \kdroot{R})$ + to achieve extrinsic pruning for query node $Q$. + The \defterm{ideal} value bounds the exact results of each query, + \begin{equation*} + \inmu(Q, R) \gets \Outopmu_{q \in Q} \fmu(q, \inqr(q, R)). + \end{equation*} + During computation, we make pessimistic bounds using only + summary statistics and completed computations. + \item Extrinsic prunes use a separate prune check $\canprunemu(\outstat(Q), \outstat(R), \inmu(Q, \kdroot{R}))$. + \item The subjective \defterm{quality} of $\lettermu$ reflects how + closely it approaches the ideal value. + Higher-quality values aid pruning. + \item Computational details omitted for brevity. + \end{itemize} +\end{slide} + +\begin{slide}{Q/R Problems - Conclusion} + \begin{itemize} + \itemt{Generality} + We maintain that our query-reference model can solve all second-order reduce problems. + Whether a problem can be solved efficiently is up to pruning rules and the execution pattern. + \itemt{Execution} + The aforementioned ``rules'' have few constraints on their execution. + Depth-first, breadth-first, and importantly {\em parallel} expansions are fair game! + \end{itemize} + We'll now focus on parallel executions. \end{slide} -\begin{slide}{Math - QR Extrinsic} - % What is mu - % How is it solved - % Quality of mu +\begin{slide}{Work Decomposition} + A dual-tree problem can be decomposed both query and reference trees. + We decompose by dividing the query tree only: + \begin{itemize} + \itemt{Private Data} + No writes to shared data. + \itemt{Quality of $\lettermu$} + For high-quality query summary results, one must consider all references. + Considering referenes independently makes this problematic. + \itemt{Scheduling} + Distant query-reference pairs complete orders of magnitude faster than + nearby ones. + Unpredictable runtimes unduly complicates scheduling. + \end{itemize} \end{slide} -\begin{slide}{Math - Summary} - % For those who didn't understand it, or - % For those who didn't get the point of it +\begin{slide}{Work Decomposition - An Item} + \begin{itemize} + \item A \defterm{work item} is a query node that is solved with a standard + serial expansion pattern. + \item We'll address the topic of decomposition two more times. + \end{itemize} \end{slide} -\begin{slide}{Task Decomposition} - % Query-based partitioning - % why it is non-obvious -\end{slide} - -\begin{slide}{Communication - Essential Trees} +\begin{slide}{Communication} + \begin{itemize} + \itemt{Fundamental Issue} Communication is a fundamental issue of parallel programming. + \itemt{Limiting Factor} Decades of experience show communication is the prime + limiting factor in parallel $N$-body methods. + \itemt{Essential Trees} Although writes are independent, many reference + points and tree nodes may be needed to perform a computation. These are + the \defterm{essential} points and nodes for a particular query node. + \end{itemize} % Introduce Domain Decomposition here \end{slide} -\begin{slide}{Communication - Nature of Communication} +\begin{slide}{Comm. - Decompositions} + \begin{itemize} + \itemt{Domain Decompositions} + A domain decomposition divides queries among processors. + \itemt{Goal} + One attempts to minimize locally essential trees, by maximizing + overlap between queries owned by a processor. (WALDO Figure goes here.) + \itemt{Past Techniques} + Neglecting scheduling details, we consider previous work: + \begin{itemize} + \itemt{ORB Tree} A $kd$-tree balanced by estimated work divides + points hierarchically among processors. + \itemt{Cost-Zones} Distribute points based on contiguous regions of + the Barnes-Hut tree traversal. + \end{itemize} + \itemt{THOR's Technique} + We combine the two by building the underlying $kd$-tree with + the domain decomposition in mind. + \end{itemize} \end{slide} -\begin{slide}{Communication - Distributed Cache System} +\begin{slide}{Comm. - Approach} + \begin{itemize} + \item We summarize previous approaches for $N$-body communication: + \begin{itemize} + \itemt{Sender Initiated} + Compute locally essential trees beforehand and send in one large pass. + {\bf Con}: Computing these trees is problem-specific. + \itemt{Use Shared Memory} + J. Singh argues you should just use large NUMA system. + {\bf Con}: NUMA systems are very expensive. + \itemt{Use Blocks} + Divide data and nodes into blocks, and request blocks on demand. + Similar to underlying NUMA systems, and allows disk-based algorithms. + \end{itemize} + \item We use blocks. Blocked methods are adaptive, have attractive + performance, and naturally allow disk-based trees. + Salmon and Warren successfully tried this a decade ago! + \end{itemize} \end{slide} -\begin{slide}{Communication - Data Layout} +\begin{slide}{Comm. - Distributed Cache} + \begin{itemize} + \itemt{Cached Array} + A cached array is distributed among machines block-by-block. + One cache one for points, one for nodes. + \itemt{Replication} + Blocks are stored on all machines that need it, but maintained by only one. + All machines have the directory. + \itemt{Consistency} + Global consistency is maintained only at infrequent sync points, once per GNP. + Between sync points, THOR writes to explicitly declared disjoint contiguous regions. + \itemt{Software Cache} + Blocks are accessed by explicitly locking and unlocking pages. + An unlocked page is placed in a queue for potential eviction. + Out-of-core blocks are requested from disk or other machines. + \end{itemize} \end{slide} -\begin{slide}{Scheduling - Problem} - % establish what the scheduling problme is - % per-task variance - % locality problem: multi-threaded and cluster +\begin{slide}{Comm. - Data Layout} + \begin{itemize} + \itemt{Point Order} + Points are stored in perfect left-to-right tree order. + For $kd$-trees, this is analagous to a Morton ordering. + \itemt{Point Blocking} + Each block of points corresponds to one tree node. + The tree-builder makes block-aligned splits until a single block + is reached, then performs midpoint splits. + \itemt{Node Order} + Previous work shows both pre-order and high-fanout orderings work well + for dual-tree algorithms. We use pre-order because it allows contiguous + groupings of nodes. + \end{itemize} \end{slide} -\begin{slide}{Scheduling - Traditional} +\begin{slide}{Scheduling} + \begin{itemize} + \itemt{Objective} + Assign tasks to processors to minimize overall run time, keeping in mind + the ``weakest link'' principle. + \itemt{Past Techniques} + The aforementioned ORB and cost-zones techniques require work estimates + per query from previous iterations. + Unfortunately, THOR cannot generally assume multiple iterations, though + future work can take this into account. + \itemt{Dynamic Scheduling} + THOR must use a dynamic scheduler to account for unforseen load imbalances. + In practice, very simple dynamic schedulers achieve good performance. + \end{itemize} \end{slide} -\begin{slide}{Scheduling - Dynamic} - % my work: - % previous work: fractiling +\begin{slide}{Scheduling in THOR} + An overview of THOR's scheduler: + \begin{itemize} + \itemt{Initial Assignment} + Divide query tree into nodes smaller than a processor is expected to handle. + Create an initial task decomposition based on the domain decomposition. + \itemt{Assignment Process} + A thread requests a work item from the master whenever it becomes idle. + Threads from the same machine are treated identically. + \itemt{Overflow} + If a machine completes its work, the spatially closest unfinished work + item is assigned. + \end{itemize} +\end{slide} + +\begin{slide}{Scheduling - Further Work} + THOR's primitive scheduler works well, but it could be better: + \begin{itemize} + \itemt{Restart} For pathological imbalances, a work item taking unusually + long should be restarted and distributed among processors. + \itemt{Distributed Assignment} Centralized scheduling becomes a bottleneck + as the number of processors grows larger. + \itemt{Hierarchical Redistribution} Rather than a greedy overflow policy + that doesn't optimize for holistic locality, overflow can be + handled by aggressively reassigning large chunks of work to large chunks + of processors. + \end{itemize} \end{slide} % QUESTION: Separate Setup/Results/Discussion sections? +\begin{slide}{Experiments} + Our experiments must be {\em realistic}: + \begin{itemize} + \itemt{Application} Our experiments must apply to problems someone + might want to solve. Finding nearest neighbors on uniform data is + silly -- just use a grid! + \itemt{Necessity} + A problem that takes a few minutes on one processor probably + doesn't need to be parallel. + If it's an hour on one processor, 128 processors is too many - consider + the bureaucratic red tape invovled. + \end{itemize} + We'll look at problems that are too large to be solved practically on + a single processor. +\end{slide} + \begin{slide}{Experiments - Setup} % algorithms and data-sets \end{slide} diff --git a/fastlib/u/rriegel/nbc/nbc.cc b/fastlib/u/rriegel/nbc/nbc.cc index 231aa96ebd..315de1669c 100644 --- a/fastlib/u/rriegel/nbc/nbc.cc +++ b/fastlib/u/rriegel/nbc/nbc.cc @@ -620,27 +620,20 @@ class Nbc { density_neg += delta.d_density_neg; } - bool ApplyPostponed(const Param& param, + void ApplyPostponed(const Param& param, const QPostponed& postponed, const QNode& q_node) { - bool change_made; - if (unlikely(postponed.label != LAB_EITHER)) { label &= postponed.label; DEBUG_ASSERT(label != LAB_NEITHER); - change_made = true; } if (unlikely(!postponed.moment_info_pos.is_empty())) { density_pos += postponed.moment_info_pos.ComputeKernelSumRange( param.kernel_pos, q_node.bound()); - change_made = true; } if (unlikely(!postponed.moment_info_neg.is_empty())) { density_neg += postponed.moment_info_neg.ComputeKernelSumRange( param.kernel_neg, q_node.bound()); - change_made = true; } - - return change_made; } }; @@ -682,7 +675,7 @@ class Nbc { void ApplyResult(const Param& param, const QPoint& q_point, index_t q_i, const QResult& result) { - fflush(stderr); + fflush(stderr); // ??? if (result.label == LAB_POS) { ++count_pos; } else if (result.label == LAB_EITHER) {