From 95c99641c07d79e030c1e2ba87b3059d64ac38db Mon Sep 17 00:00:00 2001 From: rriegel Date: Thu, 17 Jan 2008 08:39:56 +0000 Subject: [PATCH] platonic allnn complete --- fastlib/u/march/platonic_allnn/allnn.h | 802 +++++++++++-------- fastlib/u/march/platonic_allnn/allnn_main.cc | 45 +- fastlib/u/march/platonic_allnn/build.py | 54 +- 3 files changed, 535 insertions(+), 366 deletions(-) diff --git a/fastlib/u/march/platonic_allnn/allnn.h b/fastlib/u/march/platonic_allnn/allnn.h index f17bddf0b5..c9cf22059f 100644 --- a/fastlib/u/march/platonic_allnn/allnn.h +++ b/fastlib/u/march/platonic_allnn/allnn.h @@ -1,128 +1,196 @@ /** * @file allnn.h * - * Defines AllNN class to perform all-nearest-neighbors on two specified data sets. + * This file contains a "platonic" example of FASTlib code for a + * linkable library component. It implements a rudimentary dual-tree + * algorithm. For more details, see accompanying file allnn_main.cc. + * + * @see allnn_main.cc */ -// inclusion guards, please add them to your .h files -#ifndef ALLNN_H -#define ALLNN_H +// Header files should always have inclusion guards. It's a good idea +// to "sign" these guards with the containing folder or project name, +// in the off chance that someone else has a file with the same name. +#ifndef PLATONIC_ALLNN_H +#define PLATONIC_ALLNN_H -// We need to include fastlib. If you want to use fastlib, you need to have this line in addition to -// the deplibs section of your build.py +// You can include all core FASTlib components at once as follows. +// Your "deplibs" entry in build.py should mirror your includes. #include /** -* Performs all-nearest-neighbors. This class will build the trees and perform the recursive -* computation. -*/ + * A computation class for dual-tree and naive all-nearest-neighbors. + * + * This class builds trees for (assumed distinct) input query and + * reference sets on Init. The all-nearest-neighbors computation is + * then performed by calling ComputeNeighbors or ComputeNaive. + * + * This class is only intended to compute once per instantiation. + * + * Example use: + * + * @code + * AllNN allnn; + * struct datanode* allnn_module; + * ArrayList results; + * + * allnn_module = fx_submodule(NULL, "allnn", "allnn"); + * allnn.Init(query_set, reference_set, allnn_module); + * allnn.ComputeNeighbors(&results); + * @endcode + */ class AllNN { - - //////////////////////////// Nested Classes /////////////////////////////////////////////// + + ////////// Nested Classes ////////////////////////////////////////// + + private: /** - * Extra data for each node in the tree. For all nearest neighbors, each node only - * needs its upper bound on its nearest neighbor distances. + * Additional data stored at each node of a BinarySpaceTree, used by + * our QueryTree type. + * + * Dual-tree all-nearest-neighbors maintains an upper bound on + * nearest neighbor distance for each query node. */ class QueryStat { - - // Defines many useful things for a class, including a pretty printer and copy constructor - OT_DEF_BASIC(QueryStat) { - // Include this line for all non-pointer members - // There are other versions for arrays and pointers, see base/otrav.h - OT_MY_OBJECT(max_distance_so_far_); - } // OT_DEF_BASIC - + private: - /** - * The upper bound on the node's nearest neighbor distances. + * An upper bound on nearest neighbor distance for points within + * the node, to be modified after tree formation. */ double max_distance_so_far_; - + + // The object traversal macros establish a FASTlib-complient + // storage class, providing many tools including pretty printing + // and copy construction. See base/otrav.h for more details. + OT_DEF_BASIC(QueryStat) { + // Declare a non-pointer/array member variable to be traversed. + // See base/otrav.h for other kinds of declarations. + OT_MY_OBJECT(max_distance_so_far_); + } + public: - double max_distance_so_far() { - return max_distance_so_far_; - } - - + return max_distance_so_far_; + } + void set_max_distance_so_far(double new_dist) { - max_distance_so_far_ = new_dist; - } - - // In addition to any member variables for the statistic, all stat classes need two Init - // functions, one for leaves and one for non-leaves. - + max_distance_so_far_ = new_dist; + } + /** - * Initialization function used in tree-building when initializing a leaf node. For allnn, - * needs no additional information at the time of tree building. + * An Init function required by BinarySpaceTree to build + * statistics for a leaf node. + * + * All-nearest-neighbors fills statistics during computation, but + * we must set it to an appropriate starting value here. */ void Init(const Matrix& matrix, index_t start, index_t count) { // The bound starts at infinity max_distance_so_far_ = DBL_MAX; - } - - /** - * Initialization function used in tree-building when initializing a non-leaf node. For other algorithms, - * node statistics can be built using information from the children. - */ - void Init(const Matrix& matrix, index_t start, index_t count, const QueryStat& left, const QueryStat& right) { - // For allnn, non-leaves can be initialized in the same way as leaves - Init(matrix, start, count); - } - - }; //class AllNNStat - - // QueryTrees are BinarySpaceTrees where the data are bounded by Euclidean bounding boxes, - // the data are stored in a Matrix, and each node has a QueryStat for its bound. - typedef BinarySpaceTree, Matrix, QueryStat> QueryTree; - - // ReferenceTrees are the same as QueryTrees, but don't need node statistics for this algorithm. - typedef BinarySpaceTree, Matrix> ReferenceTree; - - /////////////////////////////// Members ////////////////////////////////////////////////// - private: - // These will store our data sets. - Matrix queries_; - Matrix references_; - // Pointers to the roots of the two trees. - QueryTree* query_tree_; - ReferenceTree* reference_tree_; - // The total number of prunes. - index_t number_of_prunes_; - // The module containing the parameters for this computation. - struct datanode* module_; - // A permutation of the indices for tree building. - ArrayList old_from_new_queries_; - ArrayList old_from_new_references_; - // The number of points in a leaf - index_t leaf_size_; - // The distance to the candidate nearest neighbor for each query - Vector neighbor_distances_; - // The indices of the candidate nearest neighbor for each query - ArrayList neighbor_indices_; - - - /////////////////////////////// Constructors ///////////////////////////////////////////// - - // Add this at the beginning of a class to prevent accidentally calling the copy constructor - FORBID_ACCIDENTAL_COPIES(AllNN); - + } - public: + /** + * An Init function required by BinarySpaceTree to build + * statistics from two child nodes. + * + * For all-nearest-neighbors, we reuse initialization for leaves. + */ + void Init(const Matrix& matrix, index_t start, index_t count, + const QueryStat& left, const QueryStat& right) { + Init(matrix, start, count); + } + + }; /* class AllNNStat */ + + // The tree directory defines several tools for the creation of + // custom tree types, especially for kd-trees. The DHrectBound<2> + // gives us the normal kind of kd-tree bounding boxes, using the + // 2-norm, and Matrix specifies the storage type of our data. + + /** kd-tree (binary with hrect bounds) with stats for queries. */ + typedef BinarySpaceTree, Matrix, QueryStat> QueryTree; + /** kd-tree without stats for references. */ + typedef BinarySpaceTree, Matrix> ReferenceTree; + + ////////// Members Variables /////////////////////////////////////// + + private: + /** Module used to pass parameters into the AllNN object. */ + struct datanode* module_; + + /** Copy of the query matrix given in Init. */ + Matrix queries_; + /** Copy of the reference matrix given in Init. */ + Matrix references_; + + /** Root of a tree formed on queries_. */ + QueryTree* query_tree_; + /** Root of a tree formed on references_. */ + ReferenceTree* reference_tree_; + + /** Maximum number of points in either tree's leaves. */ + index_t leaf_size_; + + /** Permutation mapping indices of queries_ to original order. */ + ArrayList old_from_new_queries_; + /** Permutation mapping indices of references_ to original order. */ + ArrayList old_from_new_references_; /** - * Constructors are generally very simple in FASTlib; most of the work is done by Init(). This is only - * responsible for ensuring that the object is ready to be destroyed safely. - */ + * Candidate nearest neighbor distances, modified during + * compuatation. Later, true nearest neighbor distances. + */ + Vector neighbor_distances_; + /** + * Candidate nearest neighbor indicies, modified during + * compuatation. Later, true nearest neighbor indices. + */ + ArrayList neighbor_indices_; + + /** Number of node-pairs pruned by the dual-tree algorithm. */ + index_t number_of_prunes_; + + + ////////// Constructors //////////////////////////////////////////// + + // It is easy to accidently call copy constructors in C++. The most + // common mistake is to define functions with object arguments + // passed by value: + // + // void foo(HugeObject x) {...} + // + // This recursively copies each member variable of the object, which + // would be disasterous, for instance, if stored query and reference + // matrices are huge. Core FASTlib components usually mitigate this + // by passing objects by const reference or by pointer: + // + // void bar(const HugeObject& x, HugeObject* y) {...} + // + // Non-const pointers are used when the outside object is modified. + // + // The following disables copy construction and assignment for + // objects of this class, which prevents functions like foo from + // compiling, saving you from poor performance and strange bugs. + FORBID_ACCIDENTAL_COPIES(AllNN); + + public: + // Default constructors should be kept very simple and should never + // allocate memory. Their two responsibilities are to ensure that + // it's safe to destroy the object without having otherwise used it + // (e.g. to set pointers to NULL) and to poison memory when in debug + // mode with BIG_BAD_NUMBER = 2146666666 = NaN as a double and + // BIG_BAD_POINTER = 0xdeadbeef. AllNN() { query_tree_ = NULL; reference_tree_ = NULL; - } - - /** - * The tree is the only member we are responsible for deleting. The others will take care of themselves. - */ + + DEBUG_POISON_PTR(module_); + DEBUG_ONLY(leaf_size_ = BIG_BAD_NUMBER); + DEBUG_ONLY(number_of_prunes_ = BIG_BAD_NUMBER); + } + + // Note that we don't delete the fx module; it's managed externally. ~AllNN() { if (query_tree_ != NULL) { delete query_tree_; @@ -130,297 +198,389 @@ class AllNN { if (reference_tree_ != NULL) { delete reference_tree_; } - } - - - /////////////////////////////// Helper Functions /////////////////////////////////////////////////// - + } + + ////////// Helper Functions //////////////////////////////////////// + /** - * Computes the minimum squared distance between the bounding boxes of two nodes + * Computes the minimum squared distance between the bounding boxes + * of two nodes. */ - double MinNodeDistSq_ (QueryTree* query_node, ReferenceTree* reference_node) { - // node->bound() gives us the DHrectBound class for the node - // It has a function MinDistanceSq which takes another DHrectBound + double MinNodeDistSq_(QueryTree* query_node, ReferenceTree* reference_node) { return query_node->bound().MinDistanceSq(reference_node->bound()); - } - - + } + + /** - * Performs exhaustive computation between two leaves. + * Performs exhaustive computation between two nodes. + * + * Note that naive also makes use of this function. */ - void ComputeBaseCase_(QueryTree* query_node, ReferenceTree* reference_node) { - - // DEBUG statements should be used frequently, since they incur no overhead - // when compiled in fast mode - - // Check that the pointers are not NULL + void GNPBaseCase_(QueryTree* query_node, ReferenceTree* reference_node) { + + // Debug checks should be used frequently. They incur no overhead + // when compiled in --mode=fast and very little otherwise. + + /* Make sure we didn't try to split children */ DEBUG_ASSERT(query_node != NULL); DEBUG_ASSERT(reference_node != NULL); - // Check that we really should be in the base case + + /* Make sure we should be in the base case */ DEBUG_WARN_IF(!query_node->is_leaf()); DEBUG_WARN_IF(!reference_node->is_leaf()); - - // Used to find the query node's new upper bound - double query_max_neighbor_distance = -1.0; - - // node->begin() is the index of the first point in the node, node->end is one past the last index - for (index_t query_index = query_node->begin(); query_index < query_node->end(); query_index++) { - - // Get the query point from the matrix + + /* Used to find the query node's new upper bound */ + double max_nearest_neighbor_distance = -1.0; + + /* Loop over all query-reference pairs */ + + // Trees don't store their points, but instead give index ranges. + // To make this feasible, they have to rearrange their input + // matrices, which is why we were sure to make copies. + for (index_t query_index = query_node->begin(); + query_index < query_node->end(); query_index++) { + + // MakeColumnVector aliases (i.e. points to but does not copy) a + // column from the matrix. + // + // A brief aside: BLAS/LAPACK is coded in Fortran and thus + // expects matrices to be column major. We side with their + // format for compatiblity, and accordingly, it is more cache + // friendly to store data points along columns, as is common in + // statistics, than along rows, as is more conventional. Vector query_point; queries_.MakeColumnVector(query_index, &query_point); - - // We'll do the same for the references - for (index_t reference_index = reference_node->begin(); reference_index < reference_node->end(); reference_index++) { - + + // It's not terrible form to leave TODO statements in code you + // intend to maintain, especially when coding under a deadline. + // These are easy to search for, though for some reason, Garry + // was more partial to "where's WALDO". More memorable, maybe? + + /* TODO: try pruning query points vs reference node */ + + for (index_t reference_index = reference_node->begin(); + reference_index < reference_node->end(); reference_index++) { + Vector reference_point; references_.MakeColumnVector(reference_index, &reference_point); - - // We'll use lapack to find the distance between the two vectors - double distance = la::DistanceSqEuclidean(query_point, reference_point); - - // If the reference point is closer than the current candidate, we'll update the candidate + + // BLAS can perform many vectors ops more quickly than C/C++. + double distance = + la::DistanceSqEuclidean(query_point, reference_point); + + /* Record points found to be closer than the best so far */ if (distance < neighbor_distances_[query_index]) { neighbor_distances_[query_index] = distance; neighbor_indices_[query_index] = reference_index; } - - - } // for reference_index - - // We need to find the upper bound distance for this query node + + } /* for reference_index */ + + /* Find the upper bound nn distance for this node */ if (neighbor_distances_[query_index] > query_max_neighbor_distance) { - query_max_neighbor_distance = neighbor_distances_[query_index]; + query_max_neighbor_distance = neighbor_distances_[query_index]; } - - } // for query_index - - // Update the upper bound for the query_node + + } /* for query_index */ + + /* Update the upper bound nn distance for the node */ query_node->stat().set_max_distance_so_far(query_max_neighbor_distance); - - - } // ComputeBaseCase_ - - + + } /* GNPBaseCase_ */ + + /** - * The recursive function + * Performs one node-node comparison in the GNP algorithm and + * recurses upon all child combinations if no prune. */ - void ComputeNeighborsRecursion_ (QueryTree* query_node, ReferenceTree* reference_node, double lower_bound_distance) { - - // DEBUG statements should be used frequently, either with or without messages - - // A DEBUG statement with no predefined message + void GNPRecursion_(QueryTree* query_node, ReferenceTree* reference_node, + double lower_bound_distance) { + + /* Make sure we didn't try to split children */ DEBUG_ASSERT(query_node != NULL); - // A DEBUG statement with a predefined message - DEBUG_ASSERT_MSG(reference_node != NULL, "reference node is null"); - // Make sure the bounding information is correct - DEBUG_ASSERT(lower_bound_distance == MinNodeDistSq_(query_node, reference_node)); - + DEBUG_ASSERT(reference_node != NULL); + + // The following asserts equality of two doubles and prints their + // values if it fails. Note that this *isn't* a particularly fast + // debug check, though; it negates the benefit of passing ahead a + // precomputed distance entirely. That's why we have --mode=fast. + + /* Make sure the precomputed bounding information is correct */ + DEBUG_SAME_DBL(lower_bound_distance, + MinNodeDistSq_(query_node, reference_node)); + if (lower_bound_distance > query_node->stat().max_distance_so_far()) { - // Pruned by distance + + /* + * A reference node with lower-bound distance greater than this + * query node's upper-bound nearest neighbor distance cannot + * contribute a reference closer than any of the queries' + * current neighbors, hence prune + */ number_of_prunes_++; - } - // node->is_leaf() works as one would expect - else if (query_node->is_leaf() && reference_node->is_leaf()) { - // Base Case - ComputeBaseCase_(query_node, reference_node); - } - else if (query_node->is_leaf()) { - // Only query is a leaf - - // We'll order the computation by distance - double left_distance = MinNodeDistSq_(query_node, reference_node->left()); - double right_distance = MinNodeDistSq_(query_node, reference_node->right()); - + + } else if (query_node->is_leaf() && reference_node->is_leaf()) { + + /* Cannot further split leaves, so process exhaustively */ + GNPBaseCase_(query_node, reference_node); + + } else if (query_node->is_leaf()) { + + /* Query node's a leaf, but we can split references */ + double left_distance = + MinNodeDistSq_(query_node, reference_node->left()); + double right_distance = + MinNodeDistSq_(query_node, reference_node->right()); + + /* + * Nearer reference node more likely to contribute neighbors + * (and thus tighten bounds), so visit it first + */ if (left_distance < right_distance) { - ComputeNeighborsRecursion_(query_node, reference_node->left(), left_distance); - ComputeNeighborsRecursion_(query_node, reference_node->right(), right_distance); + ComputeNeighborsRecursion_(query_node, + reference_node->left(), left_distance); + ComputeNeighborsRecursion_(query_node, + reference_node->right(), right_distance); + } else { + ComputeNeighborsRecursion_(query_node, + reference_node->right(), right_distance); + ComputeNeighborsRecursion_(query_node, + reference_node->left(), left_distance); } - else { - ComputeNeighborsRecursion_(query_node, reference_node->right(), right_distance); - ComputeNeighborsRecursion_(query_node, reference_node->left(), left_distance); - } - - } - - else if (reference_node->is_leaf()) { - // Only reference is a leaf - - double left_distance = MinNodeDistSq_(query_node->left(), reference_node); - double right_distance = MinNodeDistSq_(query_node->right(), reference_node); - - ComputeNeighborsRecursion_(query_node->left(), reference_node, left_distance); - ComputeNeighborsRecursion_(query_node->right(), reference_node, right_distance); - - // We need to update the upper bound based on the new upper bounds of the children + + } else if (reference_node->is_leaf()) { + + /* Reference node's a leaf, but we can split queries */ + double left_distance = + MinNodeDistSq_(query_node->left(), reference_node); + double right_distance = + MinNodeDistSq_(query_node->right(), reference_node); + + /* Order of recursion does not matter */ + ComputeNeighborsRecursion_(query_node->left(), + reference_node, left_distance); + ComputeNeighborsRecursion_(query_node->right(), + reference_node, right_distance); + + /* Update upper bound nn distance base new child bounds */ query_node->stat().set_max_distance_so_far( max(query_node->left()->stat().max_distance_so_far(), query_node->right()->stat().max_distance_so_far())); - - } - - else { - // Recurse on both as above - - double left_distance = MinNodeDistSq_(query_node->left(), reference_node->left()); - double right_distance = MinNodeDistSq_(query_node->left(), reference_node->right()); - + + } else { + + /* + * Neither node is a leaf, so split both + * + * The order we process the query node's children doesn't + * matter, but for each we should visit their nearer reference + * node first. + */ + + double left_distance = + MinNodeDistSq_(query_node->left(), reference_node->left()); + double right_distance = + MinNodeDistSq_(query_node->left(), reference_node->right()); + if (left_distance < right_distance) { - ComputeNeighborsRecursion_(query_node->left(), reference_node->left(), left_distance); - ComputeNeighborsRecursion_(query_node->left(), reference_node->right(), right_distance); + ComputeNeighborsRecursion_(query_node->left(), + reference_node->left(), left_distance); + ComputeNeighborsRecursion_(query_node->left(), + reference_node->right(), right_distance); + } else { + ComputeNeighborsRecursion_(query_node->left(), + reference_node->right(), right_distance); + ComputeNeighborsRecursion_(query_node->left(), + reference_node->left(), left_distance); } - else { - ComputeNeighborsRecursion_(query_node->left(), reference_node->right(), right_distance); - ComputeNeighborsRecursion_(query_node->left(), reference_node->left(), left_distance); - } - - left_distance = MinNodeDistSq_(query_node->right(), reference_node->left()); - right_distance = MinNodeDistSq_(query_node->right(), reference_node->right()); - + + left_distance = + MinNodeDistSq_(query_node->right(), reference_node->left()); + right_distance = + MinNodeDistSq_(query_node->right(), reference_node->right()); + if (left_distance < right_distance) { - ComputeNeighborsRecursion_(query_node->right(), reference_node->left(), left_distance); - ComputeNeighborsRecursion_(query_node->right(), reference_node->right(), right_distance); + ComputeNeighborsRecursion_(query_node->right(), + reference_node->left(), left_distance); + ComputeNeighborsRecursion_(query_node->right(), + reference_node->right(), right_distance); + } else { + ComputeNeighborsRecursion_(query_node->right(), + reference_node->right(), right_distance); + ComputeNeighborsRecursion_(query_node->right(), + reference_node->left(), left_distance); } - else { - ComputeNeighborsRecursion_(query_node->right(), reference_node->right(), right_distance); - ComputeNeighborsRecursion_(query_node->right(), reference_node->left(), left_distance); - } - - // Update the upper bound as above + + /* Update upper bound nn distance base new child bounds */ query_node->stat().set_max_distance_so_far( max(query_node->left()->stat().max_distance_so_far(), query_node->right()->stat().max_distance_so_far())); - + } - - } // ComputeNeighborsRecursion_ - - - - - ////////////////////////////////// Public Functions //////////////////////////////////////////////// - + + } /* GNPRecursion_ */ + + ////////// Public Functions //////////////////////////////////////// + + // Note that we initialize with const references below to keep from + // copying data until we want to. By the way, which side you put + // the &'s and *'s on is on the level of deep-seated religious + // belief: some people get real angry if you defy them, but you're + // really no worse a person either way. The compiler is agnostic. + /** - * Setup the class and build the trees. Note: we are initializing with const references to prevent - * local copies of the data. - */ - void Init(const Matrix& queries_in, const Matrix& references_in, struct datanode* module_in) { - - // set the module + * Read parameters, copy data into the class, and build the trees. + */ + void Init(const Matrix& queries_in, const Matrix& references_in, + struct datanode* module_in) { + module_ = module_in; - - // track the number of prunes - number_of_prunes_ = 0; - - // Get the leaf size from the module - leaf_size_ = fx_param_int(module_, "leaf_size", 20); - // Make sure the leaf size is valid - DEBUG_ASSERT(leaf_size_ > 0); - - // Copy the matrices to the class members since they will be rearranged. + + /* The data sets need to have the same number of points */ + DEBUG_SAME_SIZE(queries_.n_cols(), references_.n_cols()); + /* The data sets need to have the same dimensionality */ + DEBUG_SAME_SIZE(queries_.n_rows(), references_.n_rows()); + + /* Copy input matrices as they will be rearranged */ queries_.Copy(queries_in); references_.Copy(references_in); - - // The data sets need to have the same number of points - DEBUG_SAME_SIZE(queries_.n_rows(), references_.n_rows()); - - // Initialize the list of nearest neighbor candidates - neighbor_indices_.Init(queries_.n_cols()); - - // Initialize the vector of upper bounds for each point. - neighbor_distances_.Init(queries_.n_cols()); - neighbor_distances_.SetAll(DBL_MAX); - // We'll time tree building + leaf_size_ = fx_param_int(module_, "leaf_size", 20); + DEBUG_ASSERT(leaf_size_ > 0); + + // Timers are another handy tool provided by FASTexec. These are + // emitted automatically once you call fx_done. fx_timer_start(module_, "tree_building"); - // This call makes each tree from a matrix, leaf size, and two arrays that record the permutation of the data points - // Instead of NULL, it is possible to specify an array new_from_old_ - query_tree_ = tree::MakeKdTreeMidpoint(queries_, leaf_size_, &old_from_new_queries_, NULL); - reference_tree_ = tree::MakeKdTreeMidpoint(references_, leaf_size_, &old_from_new_references_, NULL); - - // Stop the timer we started above + // Input matrices are rearranged to an in-order traversal of + // either tree. To help in iterpretting results, the third + // argument is Init'd to a mapping from rearranged indices to the + // original order. The fourth argument, if provided, would + // initialize the reverse of said. + + /* Build the trees */ + query_tree_ = tree::MakeKdTreeMidpoint( + queries_, leaf_size_, &old_from_new_queries_, NULL); + reference_tree_ = tree::MakeKdTreeMidpoint( + references_, leaf_size_, &old_from_new_references_, NULL); + + // While we don't make use of this here, it is possible to start + // timers after stopping them. They continue where they left off. fx_timer_stop(module_, "tree_building"); - } // Init - - /** - * Initializes the AllNN structure for naive computation. This means that we simply ignore the tree building. - */ - void InitNaive(const Matrix& queries_in, const Matrix& references_in, struct datanode* module_in){ - - module_ = module_in; - - queries_.Copy(queries_in); - references_.Copy(references_in); - - DEBUG_SAME_SIZE(queries_.n_rows(), references_.n_rows()); - + /* Ready the list of nearest neighbor candidates to be filled. */ neighbor_indices_.Init(queries_.n_cols()); + + /* Ready the vector of upper bound nn distances for use. */ neighbor_distances_.Init(queries_.n_cols()); neighbor_distances_.SetAll(DBL_MAX); - - // The only difference is that we set leaf_size_ to be large enough that each tree has only one node - leaf_size_ = max(queries_.n_cols(), references_.n_cols()); - - query_tree_ = tree::MakeKdTreeMidpoint(queries_, leaf_size_, &old_from_new_queries_, NULL); - reference_tree_ = tree::MakeKdTreeMidpoint(references_, leaf_size_, &old_from_new_references_, NULL); - - } // InitNaive - - + + number_of_prunes_ = 0; + + } /* Init */ + + /** - * Computes the nearest neighbors and stores them in *results + * Initializes the AllNN structure for naive computation. + * + * We have no need to build trees for naive. + */ + void InitNaive(const Matrix& queries_in, const Matrix& references_in, + struct datanode* module_in){ + + module_ = module_in; + + /* The data sets need to have the same number of points */ + DEBUG_SAME_SIZE(queries_.n_cols(), references_.n_cols()); + /* The data sets need to have the same dimensionality */ + DEBUG_SAME_SIZE(queries_.n_rows(), references_.n_rows()); + + /* Copy input matrices */ + queries_.Copy(queries_in); + references_.Copy(references_in); + + /* + * A bit of a trick so we can still use BaseCase_: we'll expand + * the leaf size so that our trees only have one node. + */ + leaf_size_ = max(queries_.n_cols(), references_.n_cols()); + + /* Build the (single node) trees */ + query_tree_ = tree::MakeKdTreeMidpoint( + queries_, leaf_size_, &old_from_new_queries_, NULL); + reference_tree_ = tree::MakeKdTreeMidpoint( + references_, leaf_size_, &old_from_new_references_, NULL); + + /* Ready the list of nearest neighbor candidates to be filled. */ + neighbor_indices_.Init(queries_.n_cols()); + + /* Ready the vector of upper bound nn distances for use. */ + neighbor_distances_.Init(queries_.n_cols()); + neighbor_distances_.SetAll(DBL_MAX); + + number_of_prunes_ = 0; + + } /* InitNaive */ + + + /** + * Computes the nearest neighbors and stores them in results if + * provided. */ void ComputeNeighbors(ArrayList* results) { - - // Starting another timer + fx_timer_start(module_, "dual_tree_computation"); - - // Start on the root of each tree - ComputeNeighborsRecursion_(query_tree_, reference_tree_, MinNodeDistSq_(query_tree_, reference_tree_)); - - // Stopping the timer above + + /* Start recursion on the roots of either tree */ + Recursion_(query_tree_, reference_tree_, + MinNodeDistSq_(query_tree_, reference_tree_)); + fx_timer_stop(module_, "dual_tree_computation"); - - // We need to initialize the results list before filling it - results->Init(neighbor_indices_.size()); - // We need to map the indices back from how they have been permuted - for (index_t i = 0; i < neighbor_indices_.size(); i++) { - (*results)[old_from_new_queries_[i]] = old_from_new_references_[neighbor_indices_[i]]; - } - - // Save the total number of prunes to the fx module; this will appear on the command line + + // Save the total number of prunes to the FASTexec module; this + // will printed after calling fx_done or can be read back later. fx_format_result(module_, "number_of_prunes", "%d", number_of_prunes_); - - } // ComputeNeighbors - - + + if (results) { + EmitResults(results); + } + + } /* ComputeNeighbors */ + + /** - * Does the entire computation naively + * Computes the nearest neighbors naively. */ void ComputeNaive(ArrayList* results) { - - // timing for the naive computation + fx_timer_start(module_, "naive_time"); - ComputeBaseCase_(query_tree_, reference_tree_); - - // stopping the timer + /* BaseCase_ on the roots is equivalent to naive */ + BaseCase_(query_tree_, reference_tree_); + fx_timer_stop(module_, "naive_time"); - - // The same code as above - results->Init(neighbor_indices_.size()); - for (index_t i = 0; i < neighbor_indices_.size(); i++) { - (*results)[old_from_new_queries_[i]] = old_from_new_references_[neighbor_indices_[i]]; + + if (results) { + EmitResults(results); } - - } // ComputeNaive - -}; //class AllNN + } /* ComputeNaive */ + /** + * Initialize and fill an ArrayList of results. + */ + void EmitResults(ArrayList* results) { + results->Init(neighbor_indices_.size()); -#endif -// end inclusion guards + /* Map the indices back from how they have been permuted. */ + for (index_t i = 0; i < neighbor_indices_.size(); i++) { + (*results)[old_from_new_queries_[i]] = + old_from_new_references_[neighbor_indices_[i]]; + } + + } /* EmitResults */ + +}; /* class AllNN */ + +#endif /* PLATONIC_ALLNN_H */ diff --git a/fastlib/u/march/platonic_allnn/allnn_main.cc b/fastlib/u/march/platonic_allnn/allnn_main.cc index e84842cd6b..669964a9b5 100644 --- a/fastlib/u/march/platonic_allnn/allnn_main.cc +++ b/fastlib/u/march/platonic_allnn/allnn_main.cc @@ -1,20 +1,24 @@ /** - * @file allnn_main.h + * @file allnn_main.cc * - * This file contains a "platonic" example of FASTlib code. It is a - * rudimentary dual-tree algorithm for all-nearest-neighbors, but more - * importantly, it demonstrates several useful functions for common - * tasks as well as proper coding style. + * This file contains a "platonic" example of FASTlib code for a + * stand-alone executable. It makes use of an accompanying library + * implimenting a rudimentary dual-tree all-nearest-neighbors + * algorithm, but more importantly, it demonstrates useful functions + * for common tasks as well as proper coding style. * * Note however that the degree of documentation in this file well - * exceeds expectations. You should always provide Doxygen-formatted + * exceeds expectations. You should always provide Doxygen-parsed * comments (those starting with slash-star-star) for classes, their * members, and functions, but snippets of code only deserve - * documentation if it is not immediately clear what they do. (Use - * slash-slash or slash-star for line-by-line comments that Doxygen - * should ignore.) Here, we assume you are a total beginner with - * FASTlib and a novice with C++, so many additional explanations have - * been provided. + * documentation if it is not immediately clear what they do. Here, + * we assume you are a total beginner with FASTlib and a novice with + * C++, so many additional explanations have been provided. We will + * denote explanatory comments not needed in normal coding with + * slash-slash and the appropriate degree of code documentation with + * slash-star and slash-star-star. + * + * @see allnn.h */ // To begin, note the "@file" at the top of the previous comment @@ -51,7 +55,7 @@ int main(int argc, char* argv[]) { ////////// DUAL-TREE ALLNN ///////////////////////////////////////// AllNN allnn; - + // FASTexec organizes parameters and results into submodules. Think // of this as creating a new folder named "allnn_module" under the // rood directory (NULL) for the AllNN object to work inside. Here, @@ -76,7 +80,7 @@ int main(int argc, char* argv[]) { ////////// NAIVE ALLNN ///////////////////////////////////////////// - // Compare results with naive if run with "--do_naive=true". + /* Compare results with naive if run with "--do_naive=true" */ if (fx_param_bool(NULL, "do_naive", 0)) { // Our design of the AllNN class renders it usable only once; @@ -93,10 +97,11 @@ int main(int argc, char* argv[]) { ArrayList naive_results; naive_allnn.ComputeNaive(&naive_results); - // A quick sanity check, now that we have naive results. We don't - // want to run the for-loop unless debugging, hence the #ifdef. - // Most debug-only commands can instead by handled by DEBUG_ONLY - // or other debugging macros. See base/debug.h for more details. + /* Perform a quick sanity check now that we have naive results */ + + // We don't want to run the for-loop unless debugging, hence the + // #ifdef. For debug-only one-liners, use DEBUG_ONLY(expr) or the + // other debugging macros. See base/debug.h for more details. #ifdef DEBUG for (index_t i = 0; i < results.size(); ++i) { // Prints a message if results are different. Note the peculiar @@ -106,9 +111,9 @@ int main(int argc, char* argv[]) { "i = %"LI"d, results[i] = %"LI"d, naive_results[i] = %"LI"d", i, results[i], naive_results[i]); } -#endif // DEBUG +#endif /* DEBUG */ - } + } /* if do_naive */ ////////// OUTPUT RESULTS ////////////////////////////////////////// @@ -132,4 +137,4 @@ int main(int argc, char* argv[]) { // main should return 0 if the program terminates normally. return 0; -} +} /* main */ diff --git a/fastlib/u/march/platonic_allnn/build.py b/fastlib/u/march/platonic_allnn/build.py index fd9061cccb..4dc8d33f10 100644 --- a/fastlib/u/march/platonic_allnn/build.py +++ b/fastlib/u/march/platonic_allnn/build.py @@ -1,42 +1,46 @@ -# a librule creates a library, there must not be a main -# It is possible to have many librules in a single build.py +# A librule creates a library, or code that lacks a main function. +# You can define many librules in a single build.py. librule( - # What do you want the library to be called? - # Use this name to include the library in binrules elsewhere - # If the name is omitted, the library will be called the name of the directory + # What do you want the library to be called? You'll use this + # name to include the library in binrules elsewhere. If the + # omitted, the librule uses the name of build.py's directory. name = "allnn", - # Any .c or .cc files where library functions are defined - # This line can be omitted if there are no .cc or .c files + # Any .c or .cc files where library functions are defined. + # This line can be omitted if there are no .cc or .c files. #sources = ["allnn.cc"], - - # Any .h files where library functions are defined + + # Any .h files where library functions are defined. headers = ["allnn.h"], - - # libraries this library depends on - # fastlib:fastlib means the fastlib library is in the fastlib directory - # fastlib includes all the library functionality + + # Other libraries that this library depends upon. It's often + # easiest just to indicate "fastlib:fastlib", interpreted + # "directory:librule", to link with all of FASTlib's core + # components. deplibs = ["fastlib:fastlib"], - - # A file containing a main with test functions - # fl-build allnn_tests creates an executable called allnn_tests + + # You can specify a unit test file, which should contain a + # main function that runs a batch of tests. In the future, + # this will be compiled and run automatically, but for now, + # you can compile this explicitly with "fl-build allnn_tests". #tests = ["allnn_tests.cc"] ) -# a binrule creates an executable, there must be a main in one of the sources -# It is possible to have many binrules in a single build.py file +# A binrule creates an executable, or a stand-alone program that has a +# main function. It's possible to have many binrules in one build.py. binrule( - # the name of the executable + # The name of the executable. name = "allnn_main", - + + # The .c or .cc file containing main and any others you need. sources = ["allnn_main.cc"], - - # This line can be omitted if there are no headers + + # This line can be omitted if there are no headers. #headers = ["allnn_main.h"], - - # :allnn means allnn is in the same directory as this build.py file - deplibs = [":allnn", "fastlib:fastlib"] + + # The leading colon means to check this build.py for allnn. + deplibs = [":allnn"] )