From 0df2570d299fe0557eaff8a6b6116a77f800de96 Mon Sep 17 00:00:00 2001 From: Garry Boyer Date: Fri, 7 Sep 2007 23:49:16 +0000 Subject: [PATCH] breadth-first seems operational --- fastlib/thor/allnn.cc | 19 +++ fastlib/thor/build.py | 3 + fastlib/thor/dfs_impl.h | 4 +- fastlib/thor/gnp.h | 7 +- fastlib/thor/rbfs.h | 121 ++++++++++++++ fastlib/thor/rbfs_impl.h | 286 ++++++++++++++++++++++++++++++++++ fastlib/thor/thor.h | 1 + fastlib/u/garryb/nbr/fdkde.cc | 153 +++++++++++++----- fastlib/u/garryb/nbr/range.cc | 2 - fastlib/u/garryb/nbr/tkde.cc | 2 - 10 files changed, 546 insertions(+), 52 deletions(-) create mode 100644 fastlib/thor/rbfs.h create mode 100644 fastlib/thor/rbfs_impl.h diff --git a/fastlib/thor/allnn.cc b/fastlib/thor/allnn.cc index ca67c1e526..6e07656a65 100644 --- a/fastlib/thor/allnn.cc +++ b/fastlib/thor/allnn.cc @@ -147,6 +147,17 @@ class Allnn { return r_node.bound().MinDistanceSq(q_point.vec()) <= distance_sq; } + bool StartVisitingQueryPoint(const Param& param, + const QPoint& q_point, + index_t q_index, + const RNode& r_node, + QResult* q_result, + GlobalResult* global_result) { + distance_sq = q_result->distance_sq; + neighbor_i = q_result->neighbor_i; + return true; + } + void VisitPair(const Param& param, const QPoint& q_point, index_t q_index, const RPoint& r_point, index_t r_index) { @@ -168,6 +179,14 @@ class Allnn { q_result->distance_sq = distance_sq; q_result->neighbor_i = neighbor_i; } + + void FinishVisitingQueryPoint(const Param& param, + const QPoint& q_point, index_t q_index, + const RNode& r_node, + QResult* q_result, GlobalResult* global_result) { + q_result->distance_sq = distance_sq; + q_result->neighbor_i = neighbor_i; + } }; class Algorithm { diff --git a/fastlib/thor/build.py b/fastlib/thor/build.py index 4d853bb9c8..0e65a6fdad 100644 --- a/fastlib/thor/build.py +++ b/fastlib/thor/build.py @@ -6,6 +6,8 @@ cachearray_impl.h cache.h dfs.h dfs_impl.h +rbfs.h +rbfs_impl.h distribcache.h gnp.h kdtree.h @@ -30,3 +32,4 @@ librule(name = "thor", deplibs = ["fastlib:fastlib_int"]) binrule(name = "allnn", sources = ["allnn.cc"], deplibs = [":thor"]) +binrule(name = "allnnbfs", sources = ["allnnbfs.cc"], deplibs = [":thor"]) diff --git a/fastlib/thor/dfs_impl.h b/fastlib/thor/dfs_impl.h index 296de5a186..1ae9b9f43a 100644 --- a/fastlib/thor/dfs_impl.h +++ b/fastlib/thor/dfs_impl.h @@ -131,8 +131,6 @@ void DualTreeDepthFirst::Pair_( &q_node_mut->postponed)) { DEBUG_MSG(1.0, "Extrinsic prune"); } else { - global_result_.UndoDelta(param_, delta); - if (q_node->is_leaf() && r_node->is_leaf()) { DEBUG_MSG(1.0, "Base case"); BaseCase_(q_node, r_node, unvisited, q_node_mut); @@ -155,7 +153,7 @@ void DualTreeDepthFirst::Pair_( if (GNP::Algorithm::ConsiderPairIntrinsic( param_, *q_child, *r_node, &child_delta, &global_result_, &q_child_mut->postponed)) { - Pair_(q_child, r_node, delta, unvisited, q_child_mut); + Pair_(q_child, r_node, child_delta, unvisited, q_child_mut); } // We must VERY carefully apply both the horizontal and vertical join diff --git a/fastlib/thor/gnp.h b/fastlib/thor/gnp.h index 860d2234f3..a5f99b36bd 100644 --- a/fastlib/thor/gnp.h +++ b/fastlib/thor/gnp.h @@ -131,10 +131,13 @@ class BlankGlobalResult { void Init(const Param& param) {} template void Accumulate(const Param& param, const BlankGlobalResult& other) {} - template + // Delta-refinement is currently unimplemented, because we discovered it + // was buggy and neglected. if this is desired, then someone may want to + // reimplement it, but I'm not sure if it'd be that useful. + /*template void ApplyDelta(const Param& param, const Delta& delta) {} template - void UndoDelta(const Param& param, const Delta& delta) {} + void UndoDelta(const Param& param, const Delta& delta) {}*/ template void Postprocess(const Param& param) {} template diff --git a/fastlib/thor/rbfs.h b/fastlib/thor/rbfs.h new file mode 100644 index 0000000000..5a4d7b0835 --- /dev/null +++ b/fastlib/thor/rbfs.h @@ -0,0 +1,121 @@ +/** + * @file rbfs.h + * + * Depth-first dual-tree solver. + */ + +#ifndef THOR_RBFS_H +#define THOR_RBFS_H + +#include "gnp.h" +#include "cachearray.h" + +/** + * Depth-first dual-tree solver. + */ +template +class DualTreeRecursiveBreadth { + FORBID_COPY(DualTreeRecursiveBreadth); + + private: + struct QueueItem { + typename GNP::Delta delta; + index_t r_index; + }; + + struct Queue { + ArrayList q; + typename GNP::QSummaryResult summary_result; + typename GNP::QPostponed postponed; + + void Init(const typename GNP::Param& param); + + /** + * Consider a query-reference pair, possibly making an intrinsic prune. + */ + void Consider(const typename GNP::Param& param, + const typename GNP::QNode& q_node, const typename GNP::RNode& r_node, + index_t r_index, + typename GNP::GlobalResult *global_result); + + /** + * Add an existing item back to the queue, + * probably because it is a leaf. + */ + void Reconsider(const typename GNP::Param& param, + const QueueItem& item); + + void Done(const typename GNP::Param& param, + const typename GNP::QPostponed& parent_postponed, + const typename GNP::QNode& q_node); + }; + + private: + typename GNP::Param param_; + typename GNP::GlobalResult global_result_; + + CacheArray q_points_; + CacheArray q_nodes_; + CacheArray q_results_; + + CacheArray r_points_; + CacheArray r_nodes_; + const typename GNP::RNode *r_root_; + + bool do_naive_; + DualTreeRecursionStats stats_; + + public: + DualTreeRecursiveBreadth() {} + ~DualTreeRecursiveBreadth(); + + /** + * Solves the GNP. + * + * Results are stored in q_results and in this->global_result. + * The datanode contains possible parameters, and records some + * recursion statistics when debugging is enabled. + * All the other arguments are the GNP input, and are not modified. + */ + void Doit( + const typename GNP::Param& param_in, + index_t q_root_index, + index_t q_node_end_index, + DistributedCache *q_points, + DistributedCache *q_nodes, + DistributedCache *r_points, + DistributedCache *r_nodes, + DistributedCache *q_results); + + /** + * Gets the global result after computation. + */ + const typename GNP::GlobalResult& global_result() const { + return global_result_; + } + + const DualTreeRecursionStats& stats() const { + return stats_; + } + + private: + COMPILER_NOINLINE + void Begin_(index_t q_root_index); + COMPILER_NOINLINE + bool BeginExploringQueue_( + const typename GNP::QNode& q_node, Queue *parent_queue); + void Divide_(index_t q_node_i, Queue *parent_queue); + void DivideReferences_(index_t q_node_i, Queue *parent_queue); + void BaseCase_( + const typename GNP::QNode& q_node, + const typename GNP::RNode& r_node); + /** + * Postprocesses results and pushes down any postponed prunes. + */ + void PushDownPostprocess_(const typename GNP::QNode& q_node, + const typename GNP::QPostponed& postponed); +}; + +#include "rbfs_impl.h" + +#endif diff --git a/fastlib/thor/rbfs_impl.h b/fastlib/thor/rbfs_impl.h new file mode 100644 index 0000000000..d9f843c95d --- /dev/null +++ b/fastlib/thor/rbfs_impl.h @@ -0,0 +1,286 @@ +/** + * @file rbfs_impl.h + * + * Depth-first dual-tree solver template implementations. + */ + +template +DualTreeRecursiveBreadth::~DualTreeRecursiveBreadth() { + r_nodes_.StopRead(0); +} + +template +void DualTreeRecursiveBreadth::Doit( + const typename GNP::Param& param_in, + index_t q_root_index, + index_t q_end_index, + DistributedCache *q_points, + DistributedCache *q_nodes, + DistributedCache *r_points, + DistributedCache *r_nodes, + DistributedCache *q_results) { + param_.Copy(param_in); + + q_nodes_.Init(q_nodes, BlockDevice::M_READ); + r_points_.Init(r_points, BlockDevice::M_READ); + r_nodes_.Init(r_nodes, BlockDevice::M_READ); + + const typename GNP::QNode *q_root = q_nodes_.StartRead(q_root_index); + q_results_.Init(q_results, BlockDevice::M_OVERWRITE, + q_root->begin(), q_root->end()); + q_points_.Init(q_points, BlockDevice::M_READ, + q_root->begin(), q_root->end()); + q_nodes_.StopRead(q_root_index); + + global_result_.Init(param_); + + r_root_ = r_nodes_.StartRead(0); + + do_naive_ = false; + + Begin_(q_root_index); +} + +template +void DualTreeRecursiveBreadth::Begin_(index_t q_root_index) { + typename GNP::Delta delta; + CacheRead q_root(&q_nodes_, q_root_index); + + stats_.Init(); + stats_.tuples_analyzed = q_root->count() * r_root_->count(); + stats_.n_queries = q_root->count(); + + Queue queue; + + queue.Init(param_); + queue.Consider(param_, *q_root, *r_root_, 0, &global_result_); + + Divide_(q_root_index, &queue); +} + +template +void DualTreeRecursiveBreadth::PushDownPostprocess_( + const typename GNP::QNode& q_node, + const typename GNP::QPostponed& postponed) { + if (q_node.is_leaf()) { + index_t q_i = q_node.begin(); + CacheWriteIter q_result(&q_results_, q_i); + CacheReadIter q_point(&q_points_, q_i); + + for (; q_i < q_node.end(); q_i++, q_result.Next(), q_point.Next()) { + q_result->ApplyPostponed(param_, postponed, *q_point, q_i); + q_result->Postprocess(param_, *q_point, q_i, *r_root_); + global_result_.ApplyResult(param_, *q_point, q_i, *q_result); + } + } else { + for (int k = 0; k < GNP::QNode::CARDINALITY; k++) { + CacheRead q_child(&q_nodes_, q_node.child(k)); + + PushDownPostprocess_(*q_child, postponed); + } + } +} + +template +bool DualTreeRecursiveBreadth::BeginExploringQueue_( + const typename GNP::QNode& q_node, Queue *parent_queue) { + if (parent_queue->q.size() == 0 + || !GNP::Algorithm::ConsiderQueryTermination( + param_, q_node, parent_queue->summary_result, + global_result_, &parent_queue->postponed)) { + // Distribute mass results to the leaves + PushDownPostprocess_(q_node, parent_queue->postponed); + return false; + } else { + return true; + } +} + +template +void DualTreeRecursiveBreadth::Queue::Init( + const typename GNP::Param& param) { + q.Init(); + summary_result.Init(param); + postponed.Init(param); +} + +template +void DualTreeRecursiveBreadth::Queue::Consider( + const typename GNP::Param& param, + const typename GNP::QNode& q_node, const typename GNP::RNode& r_node, + index_t r_index, + typename GNP::GlobalResult *global_result) { + QueueItem *item = q.AddBack(); + item->r_index = r_index; + if (likely(GNP::Algorithm::ConsiderPairIntrinsic(param, q_node, r_node, + &item->delta, global_result, &postponed))) { + summary_result.ApplyDelta(param, item->delta); + } else { + q.PopBack(); + } +} + +template +void DualTreeRecursiveBreadth::Queue::Reconsider( + const typename GNP::Param& param, + const QueueItem& item) { + new(q.AddBack())QueueItem(item); + summary_result.ApplyDelta(param, item.delta); +} + +template +void DualTreeRecursiveBreadth::Queue::Done( + const typename GNP::Param& param, + const typename GNP::QPostponed& parent_postponed, + const typename GNP::QNode& q_node) { + postponed.ApplyPostponed(param, parent_postponed); + summary_result.ApplyPostponed(param, postponed, q_node); +} + +template +void DualTreeRecursiveBreadth::DivideReferences_( + index_t q_node_i, Queue* parent_queue) { + const typename GNP::QNode q_node(*q_nodes_.StartRead(q_node_i)); + q_nodes_.StopRead(q_node_i); + + if (!BeginExploringQueue_(q_node, parent_queue)) { + return; + } + + Queue child_queue; + child_queue.Init(param_); + + DEBUG_ONLY(stats_.node_node_considered += parent_queue->q.size()); + + for (index_t i = 0; i < parent_queue->q.size(); i++) { + const QueueItem *item = &parent_queue->q[i]; + CacheRead r_node(&r_nodes_, item->r_index); + + if (likely(GNP::Algorithm::ConsiderPairExtrinsic( + param_, q_node, *r_node, item->delta, parent_queue->summary_result, + global_result_, &parent_queue->postponed))) { + if (!r_node->is_leaf()) { + for (int k_r = 0; k_r < GNP::RNode::CARDINALITY; k_r++) { + index_t r_child_i = r_node->child(k_r); + CacheRead r_child(&r_nodes_, r_child_i); + + child_queue.Consider(param_, q_node, *r_child, r_child_i, + &global_result_); + } + } else { + BaseCase_(q_node, *r_node); + } + } + } + + child_queue.Done(param_, parent_queue->postponed, q_node); + DivideReferences_(q_node_i, &child_queue); +} + +template +void DualTreeRecursiveBreadth::Divide_( + index_t q_node_i, Queue* parent_queue) { + const typename GNP::QNode q_node(*q_nodes_.StartRead(q_node_i)); + q_nodes_.StopRead(q_node_i); + + if (q_node.is_leaf()) { + DivideReferences_(q_node_i, parent_queue); + return; + } + + if (!BeginExploringQueue_(q_node, parent_queue)) { + return; + } + + Queue child_queues[GNP::QNode::CARDINALITY]; + const typename GNP::QNode *q_children[GNP::QNode::CARDINALITY]; + + for (int k = 0; k < GNP::QNode::CARDINALITY; k++) { + q_children[k] = q_nodes_.StartRead(q_node.child(k)); + child_queues[k].Init(param_); + } + + DEBUG_ONLY(stats_.node_node_considered += parent_queue->q.size()); + + for (index_t i = 0; i < parent_queue->q.size(); i++) { + const QueueItem *item = &parent_queue->q[i]; + CacheRead r_node(&r_nodes_, item->r_index); + + if (likely(GNP::Algorithm::ConsiderPairExtrinsic( + param_, q_node, *r_node, item->delta, parent_queue->summary_result, + global_result_, &parent_queue->postponed))) { + if (!r_node->is_leaf()) { + for (int k_r = 0; k_r < GNP::RNode::CARDINALITY; k_r++) { + index_t r_child_i = r_node->child(k_r); + CacheRead r_child(&r_nodes_, r_child_i); + + for (int k_q = 0; k_q < GNP::QNode::CARDINALITY; k_q++) { + child_queues[k_q].Consider(param_, *q_children[k_q], *r_child, + r_child_i, &global_result_); + } + } + } else { + for (int k_q = 0; k_q < GNP::QNode::CARDINALITY; k_q++) { + child_queues[k_q].Reconsider(param_, *item); + } + } + } + } + + // Release the locks on the children to ease cache pressure in the FIFO + for (int k = 0; k < GNP::QNode::CARDINALITY; k++) { + child_queues[k].Done(param_, parent_queue->postponed, *q_children[k]); + q_nodes_.StopRead(q_node.child(k)); + } + + for (int k = 0; k < GNP::QNode::CARDINALITY; k++) { + Divide_(q_node.child(k), &child_queues[k]); + } +} + +template +void DualTreeRecursiveBreadth::BaseCase_( + const typename GNP::QNode& q_node, + const typename GNP::RNode& r_node) { + DEBUG_ONLY(stats_.node_point_considered += q_node.count()); + + typename GNP::PairVisitor visitor; + visitor.Init(param_); + + CacheRead first_q_point(&q_points_, q_node.begin()); + CacheWrite first_q_result(&q_results_, q_node.begin()); + CacheRead first_r_point(&r_points_, r_node.begin()); + size_t q_point_stride = q_points_.n_elem_bytes(); + size_t q_result_stride = q_results_.n_elem_bytes(); + size_t r_point_stride = r_points_.n_elem_bytes(); + index_t q_end = q_node.end(); + const typename GNP::QPoint *q_point = first_q_point; + typename GNP::QResult *q_result = first_q_result; + + for (index_t q_i = q_node.begin(); q_i < q_end; ++q_i) { + if (visitor.StartVisitingQueryPoint(param_, *q_point, q_i, r_node, + q_result, &global_result_)) { + const typename GNP::RPoint *r_point = first_r_point; + index_t r_i = r_node.begin(); + index_t r_left = r_node.count(); + + for (;;) { + visitor.VisitPair(param_, *q_point, q_i, *r_point, r_i); + if (unlikely(--r_left == 0)) { + break; + } + r_i++; + r_point = mem::PointerAdd(r_point, r_point_stride); + } + + visitor.FinishVisitingQueryPoint(param_, *q_point, q_i, r_node, + q_result, &global_result_); + + DEBUG_ONLY(stats_.point_point_considered += r_node.count()); + } + + q_point = mem::PointerAdd(q_point, q_point_stride); + q_result = mem::PointerAdd(q_result, q_result_stride); + } +} + diff --git a/fastlib/thor/thor.h b/fastlib/thor/thor.h index 4668fb1ef8..7688841b4d 100644 --- a/fastlib/thor/thor.h +++ b/fastlib/thor/thor.h @@ -9,6 +9,7 @@ #include "gnp.h" #include "dfs.h" +#include "rbfs.h" #include "thor_utils.h" #include "kdtree.h" diff --git a/fastlib/u/garryb/nbr/fdkde.cc b/fastlib/u/garryb/nbr/fdkde.cc index 2c8c76280b..2e4f473c36 100644 --- a/fastlib/u/garryb/nbr/fdkde.cc +++ b/fastlib/u/garryb/nbr/fdkde.cc @@ -1,6 +1,9 @@ #include "fastlib/fastlib_int.h" #include "thor/thor.h" +#define SOLVER_TYPE DualTreeRecursiveBreadth +//#define SOLVER_TYPE DualTreeDepthFirst + /** * Approximate kernel density estimation. * @@ -30,8 +33,6 @@ class FdKde { Kernel kernel; /** The amount of relative error prooportional to local lower bound. */ double rel_error_local; - /** The amount of relative error distributed uniformly. */ - double rel_error_global; /** The dimensionality of the data sets. */ index_t dim; @@ -43,13 +44,14 @@ class FdKde { double rel_error; /** Amount of error that is local error */ double p_local; + /** Amount of error that is global error */ + double p_global; /** The band width, h. */ double bandwidth; OT_DEF_BASIC(Param) { OT_MY_OBJECT(kernel); OT_MY_OBJECT(rel_error_local); - OT_MY_OBJECT(rel_error_global); OT_MY_OBJECT(dim); OT_MY_OBJECT(count); OT_MY_OBJECT(mul_constant); @@ -64,16 +66,16 @@ class FdKde { */ void Init(datanode *module) { bandwidth = fx_param_double_req(module, "h"); - p_local = fx_param_double(module, "p_local", 0.5); + p_local = fx_param_double(module, "p_local", 0); + p_global = 1 - p_local; rel_error = fx_param_double(module, "rel_error", 0.1); + rel_error_local = rel_error * p_local; } /** this is called after things are set. */ void SetDimensions() { kernel.Init(bandwidth, dim); mul_constant = 1.0 / (kernel.CalcNormConstant(dim) * (count - 1)); - rel_error_local = rel_error * p_local; - rel_error_global = rel_error * (1.0 - p_local) / count; } }; @@ -146,22 +148,27 @@ class FdKde { public: /** The density contribution postponed. */ DRange d_density; + index_t n_pruned; OT_DEF_BASIC(QPostponed) { OT_MY_OBJECT(d_density); + OT_MY_OBJECT(n_pruned); } public: void Init(const Param& param) { d_density.Init(0, 0); + n_pruned = 0; } void Reset(const Param& param) { d_density.Init(0, 0); + n_pruned = 0; } void ApplyPostponed(const Param& param, const QPostponed& other) { d_density += other.d_density; + n_pruned += other.n_pruned; } }; @@ -186,14 +193,17 @@ class FdKde { struct QResult { public: DRange density; + index_t n_pruned; OT_DEF_BASIC(QResult) { OT_MY_OBJECT(density); + OT_MY_OBJECT(n_pruned); } public: void Init(const Param& param) { density.Init(0, 0); + n_pruned = 0; } void Postprocess(const Param& param, @@ -208,6 +218,7 @@ class FdKde { const QPostponed& postponed, const QPoint& q, index_t q_index) { density += postponed.d_density; + n_pruned += postponed.n_pruned; } }; @@ -215,52 +226,69 @@ class FdKde { public: /** Bound on density from leaves. */ DRange density; + double used_error; + index_t n_pruned; OT_DEF_BASIC(QSummaryResult) { OT_MY_OBJECT(density); + OT_MY_OBJECT(used_error); + OT_MY_OBJECT(n_pruned); } public: void Init(const Param& param) { /* horizontal init */ density.Init(0, 0); - } - - void StartReaccumulate(const Param& param, const QNode& q_node) { - /* vertical init */ - density.InitEmptySet(); - } - - void Accumulate(const Param& param, const QResult& result) { - // TODO: applying to single result could be made part of QResult, - // but in some cases may require a copy/undo stage - density |= result.density; - } - - void Accumulate(const Param& param, - const QSummaryResult& result, index_t n_points) { - density |= result.density; - } - - void FinishReaccumulate(const Param& param, - const QNode& q_node) { - /* no post-processing steps necessary */ + used_error = 0; + n_pruned = 0; } /** horizontal join operator */ void ApplySummaryResult(const Param& param, const QSummaryResult& summary_result) { density += summary_result.density; + used_error += summary_result.used_error; + n_pruned += summary_result.n_pruned; } void ApplyDelta(const Param& param, const Delta& delta) { density += delta.d_density; + // delta's don't affect used error } void ApplyPostponed(const Param& param, const QPostponed& postponed, const QNode& q_node) { density += postponed.d_density; + used_error += postponed.d_density.width() / 2; + n_pruned += postponed.n_pruned; + } + + void StartReaccumulate(const Param& param, const QNode& q_node) { + /* vertical init */ + density.InitEmptySet(); + used_error = 0; + n_pruned = param.count; + } + + void Accumulate(const Param& param, const QResult& result) { + // TODO: applying to single result could be made part of QResult, + // but in some cases may require a copy/undo stage + density |= result.density; + used_error = max(used_error, result.density.width() / 2); + n_pruned = min(n_pruned, result.n_pruned); + } + + void Accumulate(const Param& param, + const QSummaryResult& result, index_t n_points) { + density |= result.density; + used_error = max(used_error, result.used_error); + n_pruned = min(n_pruned, result.n_pruned); + } + + void FinishReaccumulate(const Param& param, + const QNode& q_node) { + /* no post-processing steps necessary */ } }; @@ -286,13 +314,13 @@ class FdKde { sum_density += other.sum_density; foo += other.foo; } - void ApplyDelta(const Param& param, const Delta& delta) {} - void UndoDelta(const Param& param, const Delta& delta) {} void Postprocess(const Param& param) {} void Report(const Param& param, datanode *datanode) { fx_format_result(datanode, "avg_density_lo", "%g", sum_density.lo / param.count); fx_format_result(datanode, "avg_density_hi", "%g", sum_density.hi / param.count); fx_format_result(datanode, "avg_density", "%g", sum_density.mid() / param.count); + fx_format_result(datanode, "avg_rel_error", "%g", + sum_density.width() / sum_density.lo / 2); fx_format_result(datanode, "foo", "%g", foo / param.count); } void ApplyResult(const Param& param, @@ -323,17 +351,23 @@ class FdKde { const QSummaryResult& unapplied_summary_results, QResult* q_result, GlobalResult* global_result) { + q_result->n_pruned += r_node.count(); + DRange distance_sq_range = DRange( r_node.bound().MinDistanceSq(q.vec()), r_node.bound().MaxDistanceSq(q.vec())); DEBUG_ASSERT(distance_sq_range.width() >= 0); DRange d_density = param.kernel.RangeUnnormOnSq(distance_sq_range); - double summary_density_lo = unapplied_summary_results.density.lo - + d_density.lo + q_result->density.lo; + double density_lo = + (d_density.lo + unapplied_summary_results.density.lo + q_result->density.lo); + double allocated_error = - param.rel_error_local * d_density.lo - + param.rel_error_global * summary_density_lo * r_node.count(); + (param.rel_error * density_lo + - (q_result->density.width() / 2)) + / (param.count - q_result->n_pruned) * param.p_global; + allocated_error *= r_node.count(); + allocated_error += param.rel_error_local * d_density.lo; if (d_density.width() < allocated_error * 2) { q_result->density += d_density * r_node.count(); @@ -345,6 +379,21 @@ class FdKde { return true; } + /** + * This is the lame form of the function used by breadth-first. + * + * Since breadth-first tries to avoid getting to leaves anyways, it + * doesn't want to bother with giving you summary results, so it doesn't. + */ + bool StartVisitingQueryPoint(const Param& param, + const QPoint& q, index_t q_index, + const RNode& r_node, + QResult* q_result, + GlobalResult* global_result) { + density = 0; + return true; + } + void VisitPair(const Param& param, const QPoint& q, index_t q_index, const RPoint& r, index_t r_index) { @@ -353,11 +402,18 @@ class FdKde { } void FinishVisitingQueryPoint(const Param& param, - const QPoint& q, index_t q_index, - const RNode& r_node, + const QPoint& q, index_t q_index, const RNode& r_node, const QSummaryResult& unapplied_summary_results, - QResult* q_result, - GlobalResult* global_result) { + QResult* q_result, GlobalResult* global_result) { + q_result->density += density; + } + + /** + * Once again, the lame form for breadth-first. + */ + void FinishVisitingQueryPoint(const Param& param, + const QPoint& q, index_t q_index, const RNode& r_node, + QResult* q_result, GlobalResult* global_result) { q_result->density += density; } }; @@ -386,9 +442,16 @@ class FdKde { delta->d_density = param.kernel.RangeUnnormOnSq(distance_sq_range); delta->d_density *= r_node.count(); - DEBUG_ASSERT(delta->d_density.lo <= delta->d_density.hi); + DEBUG_ASSERT_MSG(delta->d_density.lo <= delta->d_density.hi * (1 + 1.0e-7), + "delta density lo %f > hi %f", + delta->d_density.lo, delta->d_density.hi); - return likely(delta->d_density.hi != 0); + if (likely(delta->d_density.hi != 0)) { + return true; + } else { + q_postponed->n_pruned += r_node.count(); + return false; + } } static bool ConsiderPairExtrinsic( @@ -400,11 +463,15 @@ class FdKde { const GlobalResult& global_result, QPostponed* q_postponed) { double allocated_error = - param.rel_error_local * delta.d_density.lo - + param.rel_error_global * q_summary_result.density.lo * r_node.count(); + (param.rel_error * q_summary_result.density.lo + - q_summary_result.used_error) + / (param.count - q_summary_result.n_pruned) * param.p_global; + allocated_error *= r_node.count(); + allocated_error += param.rel_error_local * delta.d_density.lo; if (delta.d_density.width() < allocated_error * 2) { q_postponed->d_density += delta.d_density; + q_postponed->n_pruned += r_node.count(); return false; } @@ -484,7 +551,7 @@ class FdKde { GlobalResult global_result_1; fx_timer_start(module, "kde_1"); - thor::RpcDualTree >( + thor::RpcDualTree >( fx_submodule(module, "gnp", "kde_1"), 200, *param, &tree, &tree, &results, &global_result_1); fx_timer_stop(module, "kde_1"); @@ -498,7 +565,7 @@ class FdKde { param->kernel.Init(sqrt(param->kernel.bandwidth_sq() * 2)); fx_timer_start(module, "kde_2"); - thor::RpcDualTree >( + thor::RpcDualTree >( fx_submodule(module, "gnp", "kde_2"), 200, *param, &tree, &tree, &results, &global_result_2); fx_timer_stop(module, "kde_2"); diff --git a/fastlib/u/garryb/nbr/range.cc b/fastlib/u/garryb/nbr/range.cc index 15a0058286..4e2bca4a76 100644 --- a/fastlib/u/garryb/nbr/range.cc +++ b/fastlib/u/garryb/nbr/range.cc @@ -125,8 +125,6 @@ class Range { void Accumulate(const Param& param, const GlobalResult& other) { count += other.count; } - void ApplyDelta(const Param& param, const Delta& delta) {} - void UndoDelta(const Param& param, const Delta& delta) {} void Postprocess(const Param& param) {} void Report(const Param& param, datanode *datanode) { fx_format_result(datanode, "per_point_avg", "%g", diff --git a/fastlib/u/garryb/nbr/tkde.cc b/fastlib/u/garryb/nbr/tkde.cc index 0e6fd637f5..c064d3254a 100644 --- a/fastlib/u/garryb/nbr/tkde.cc +++ b/fastlib/u/garryb/nbr/tkde.cc @@ -435,8 +435,6 @@ class Tkde { n_under_threshold += other.n_under_threshold; n_unknown += other.n_unknown; } - void ApplyDelta(const Param& param, const Delta& delta) {} - void UndoDelta(const Param& param, const Delta& delta) {} void Postprocess(const Param& param) {} void Report(const Param& param, datanode *datanode) { fx_format_result(datanode, "n_unknown", "%"LI"d",