apaam/mlpack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

platonic allnn complete

Browse Source
This commit is contained in:
rriegel
2008-01-17 08:39:56 +00:00
parent 78ac134db3
commit 95c99641c0
3 changed files with 535 additions and 366 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib/u/march/platonic_allnn/allnn.h
+481 -321
View File
@@ -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 <fastlib/fastlib.h>
/**
* 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<index_t> 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<DHrectBound<2>, Matrix, QueryStat> QueryTree;
// ReferenceTrees are the same as QueryTrees, but don't need node statistics for this algorithm.
typedef BinarySpaceTree<DHrectBound<2>, 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<index_t> old_from_new_queries_;
ArrayList<index_t> 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<index_t> 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<DHrectBound<2>, Matrix, QueryStat> QueryTree;
/** kd-tree without stats for references. */
typedef BinarySpaceTree<DHrectBound<2>, 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<index_t> old_from_new_queries_;
/** Permutation mapping indices of references_ to original order. */
ArrayList<index_t> 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<index_t> 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<QueryTree>(queries_, leaf_size_, &old_from_new_queries_, NULL);
reference_tree_ = tree::MakeKdTreeMidpoint<ReferenceTree>(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<QueryTree>(
queries_, leaf_size_, &old_from_new_queries_, NULL);
reference_tree_ = tree::MakeKdTreeMidpoint<ReferenceTree>(
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<QueryTree>(queries_, leaf_size_, &old_from_new_queries_, NULL);
reference_tree_ = tree::MakeKdTreeMidpoint<ReferenceTree>(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<QueryTree>(
queries_, leaf_size_, &old_from_new_queries_, NULL);
reference_tree_ = tree::MakeKdTreeMidpoint<ReferenceTree>(
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<index_t>* 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<index_t>* 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<index_t>* 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 */
fastlib/u/march/platonic_allnn/allnn_main.cc
+25 -20
View File
@@ -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<index_t> 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 */
fastlib/u/march/platonic_allnn/build.py
+29 -25
View File
@@ -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"]
)