diff --git a/src/mlpack/core/tree/CMakeLists.txt b/src/mlpack/core/tree/CMakeLists.txt index 2a849ae114..77e9026db4 100644 --- a/src/mlpack/core/tree/CMakeLists.txt +++ b/src/mlpack/core/tree/CMakeLists.txt @@ -50,6 +50,14 @@ set(SOURCES hollow_ball_bound_impl.hpp hrectbound.hpp hrectbound_impl.hpp + octree.hpp + octree/octree.hpp + octree/octree_impl.hpp + octree/single_tree_traverser.hpp + octree/single_tree_traverser_impl.hpp + octree/dual_tree_traverser.hpp + octree/dual_tree_traverser_impl.hpp + octree/traits.hpp rectangle_tree.hpp rectangle_tree/rectangle_tree.hpp rectangle_tree/rectangle_tree_impl.hpp diff --git a/src/mlpack/core/tree/binary_space_tree/traits.hpp b/src/mlpack/core/tree/binary_space_tree/traits.hpp index 8e6841a591..4ebdb04cdc 100644 --- a/src/mlpack/core/tree/binary_space_tree/traits.hpp +++ b/src/mlpack/core/tree/binary_space_tree/traits.hpp @@ -14,7 +14,7 @@ namespace mlpack { namespace tree { /** - * This is a specialization of the TreeType class to the BinarySpaceTree tree + * This is a specialization of the TreeTraits class to the BinarySpaceTree tree * type. It defines characteristics of the binary space tree, and is used to * help write tree-independent (but still optimized) tree-based algorithms. See * mlpack/core/tree/tree_traits.hpp for more information. diff --git a/src/mlpack/core/tree/octree.hpp b/src/mlpack/core/tree/octree.hpp new file mode 100644 index 0000000000..ed72ae5002 --- /dev/null +++ b/src/mlpack/core/tree/octree.hpp @@ -0,0 +1,17 @@ +/** + * @file octree.hpp + * @author Ryan Curtin + * + * Include all the necessary files to use the Octree class. + */ +#ifndef MLPACK_CORE_TREE_OCTREE_HPP +#define MLPACK_CORE_TREE_OCTREE_HPP + +#include +#include "bounds.hpp" +#include "octree/octree.hpp" +#include "octree/traits.hpp" +#include "octree/single_tree_traverser.hpp" +#include "octree/dual_tree_traverser.hpp" + +#endif diff --git a/src/mlpack/core/tree/octree/dual_tree_traverser.hpp b/src/mlpack/core/tree/octree/dual_tree_traverser.hpp new file mode 100644 index 0000000000..ec9774a30f --- /dev/null +++ b/src/mlpack/core/tree/octree/dual_tree_traverser.hpp @@ -0,0 +1,78 @@ +/** + * @file dual_tree_traverser.hpp + * @author Ryan Curtin + * + * Define the dual-tree traverser for the Octree. + */ +#ifndef MLPACK_CORE_TREE_OCTREE_DUAL_TREE_TRAVERSER_HPP +#define MLPACK_CORE_TREE_OCTREE_DUAL_TREE_TRAVERSER_HPP + +#include +#include "octree.hpp" + +namespace mlpack { +namespace tree { + +template +template +class Octree::DualTreeTraverser +{ + public: + /** + * Instantiate the given dual-tree traverser with the given rule set. + */ + DualTreeTraverser(RuleType& rule); + + /** + * Traverse the two trees. This does not reset the statistics of the + * traversals (it just adds to them). + */ + void Traverse(Octree& queryNode, Octree& referenceNode); + + //! Get the number of pruned nodes. + size_t NumPrunes() const { return numPrunes; } + //! Modify the number of pruned nodes (i.e. to reset it). + size_t& NumPrunes() { return numPrunes; } + + //! Get the number of visited node combinations. + size_t NumVisited() const { return numVisited; } + //! Modify the number of visited node combinations. + size_t& NumVistied() { return numVisited; } + + //! Get the number of times a node was scored. + size_t NumScores() const { return numScores; } + //! Modify the number of times a node was scored. + size_t& NumScores() { return numScores; } + + //! Get the number of times a base case was computed. + size_t NumBaseCases() const { return numBaseCases; } + //! Modify the number of times a base case was computed. + size_t& NumBaseCases() { return numBaseCases; } + + private: + //! The rule type to use. + RuleType& rule; + + //! The number of prunes. + size_t numPrunes; + //! The number of visited node combinations. + size_t numVisited; + //! The number of times a node was scored. + size_t numScores; + //! The number of times a base case was calculated. + size_t numBaseCases; + + //! Traversal information, held in the class so that it isn't continually + //! being reallocated. + typename RuleType::TraversalInfoType traversalInfo; +}; + +} // namespace tree +} // namespace mlpack + +// Include implementation. +#include "dual_tree_traverser_impl.hpp" + +#endif diff --git a/src/mlpack/core/tree/octree/dual_tree_traverser_impl.hpp b/src/mlpack/core/tree/octree/dual_tree_traverser_impl.hpp new file mode 100644 index 0000000000..7058adb002 --- /dev/null +++ b/src/mlpack/core/tree/octree/dual_tree_traverser_impl.hpp @@ -0,0 +1,147 @@ +/** + * @file dual_tree_traverser_impl.hpp + * @author Ryan Curtin + * + * Implementation of the dual-tree traverser for the octree. + */ +#ifndef MLPACK_CORE_TREE_OCTREE_DUAL_TREE_TRAVERSER_IMPL_HPP +#define MLPACK_CORE_TREE_OCTREE_DUAL_TREE_TRAVERSER_IMPL_HPP + +// In case it hasn't been included yet. +#include "dual_tree_traverser.hpp" + +namespace mlpack { +namespace tree { + +template +template +Octree::DualTreeTraverser:: + DualTreeTraverser(RuleType& rule) : + rule(rule), + numPrunes(0), + numVisited(0), + numScores(0), + numBaseCases(0) +{ + // Nothing to do. +} + +template +template +void Octree::DualTreeTraverser:: + Traverse(Octree& queryNode, Octree& referenceNode) +{ + // Increment the visit counter. + ++numVisited; + + // Store the current traversal info. + traversalInfo = rule.TraversalInfo(); + + if (queryNode.IsLeaf() && referenceNode.IsLeaf()) + { + const size_t begin = queryNode.Point(0); + const size_t end = begin + queryNode.NumPoints(); + for (size_t q = begin; q < end; ++q) + { + // First, see if we can prune the reference node for this query point. + rule.TraversalInfo() = traversalInfo; + const double score = rule.Score(q, referenceNode); + if (score == DBL_MAX) + { + ++numPrunes; + continue; + } + + const size_t rBegin = referenceNode.Point(0); + const size_t rEnd = rBegin + referenceNode.NumPoints(); + for (size_t r = rBegin; r < rEnd; ++r) + rule.BaseCase(q, r); + + numBaseCases += referenceNode.NumPoints(); + } + } + else if (!queryNode.IsLeaf() && referenceNode.IsLeaf()) + { + // We have to recurse down the query node. Order does not matter. + for (size_t i = 0; i < queryNode.NumChildren(); ++i) + { + rule.TraversalInfo() = traversalInfo; + const double score = rule.Score(queryNode.Child(i), referenceNode); + if (score == DBL_MAX) + { + ++numPrunes; + continue; + } + + Traverse(queryNode.Child(i), referenceNode); + } + } + else if (queryNode.IsLeaf() && !referenceNode.IsLeaf()) + { + // We have to recurse down the reference node, so we need to do it in an + // ordered manner. + arma::vec scores(referenceNode.NumChildren()); + std::vector tis; + for (size_t i = 0; i < referenceNode.NumChildren(); ++i) + { + rule.TraversalInfo() = traversalInfo; + scores[i] = rule.Score(queryNode, referenceNode.Child(i)); + tis.push_back(rule.TraversalInfo()); + } + + // Sort the scores. + arma::uvec scoreOrder = arma::sort_index(scores); + for (size_t i = 0; i < scoreOrder.n_elem; ++i) + { + if (scores[scoreOrder[i]] == DBL_MAX) + { + // We don't need to check any more---all children past here are pruned. + numPrunes += scoreOrder.n_elem - i; + break; + } + + rule.TraversalInfo() = tis[scoreOrder[i]]; + Traverse(queryNode, referenceNode.Child(scoreOrder[i])); + } + } + else + { + // We have to recurse down both the query and reference nodes. Query order + // does not matter, so we will do that in sequence. However we will + // allocate the arrays for recursion at this level. + arma::vec scores(referenceNode.NumChildren()); + std::vector + tis(referenceNode.NumChildren()); + for (size_t j = 0; j < queryNode.NumChildren(); ++j) + { + // Now we have to recurse down the reference node, which we will do in a + // prioritized manner. + for (size_t i = 0; i < referenceNode.NumChildren(); ++i) + { + rule.TraversalInfo() = traversalInfo; + scores[i] = rule.Score(queryNode.Child(j), referenceNode.Child(i)); + tis[i] = rule.TraversalInfo(); + } + + // Sort the scores. + arma::uvec scoreOrder = arma::sort_index(scores); + for (size_t i = 0; i < scoreOrder.n_elem; ++i) + { + if (scores[scoreOrder[i]] == DBL_MAX) + { + // We don't need to check any more---all children past here are pruned. + numPrunes += scoreOrder.n_elem - i; + break; + } + + rule.TraversalInfo() = tis[scoreOrder[i]]; + Traverse(queryNode.Child(j), referenceNode.Child(scoreOrder[i])); + } + } + } +} + +} // namespace tree +} // namespace mlpack + +#endif diff --git a/src/mlpack/core/tree/octree/octree.hpp b/src/mlpack/core/tree/octree/octree.hpp index 544fa754e7..0a8ffd2489 100644 --- a/src/mlpack/core/tree/octree/octree.hpp +++ b/src/mlpack/core/tree/octree/octree.hpp @@ -9,11 +9,12 @@ #include #include "../hrectbound.hpp" +#include "../statistic.hpp" namespace mlpack { namespace tree { -template class Octree @@ -24,6 +25,14 @@ class Octree //! The type of element held in MatType. typedef typename MatType::elem_type ElemType; + //! A single-tree traverser; see single_tree_traverser.hpp. + template + class SingleTreeTraverser; + + //! A dual-tree traverser; see dual_tree_traverser.hpp. + template + class DualTreeTraverser; + private: //! The children held by this node. std::vector children; @@ -36,11 +45,19 @@ class Octree size_t count; //! The minimum bounding rectangle of the points held in the node (and its //! children). - HRectBound bound; + bound::HRectBound bound; //! The dataset. MatType* dataset; //! The parent (NULL if this node is the root). Octree* parent; + //! The statistic. + StatisticType stat; + //! The distance from the center of this node to the center of the parent. + ElemType parentDistance; + //! The distance to the furthest descendant, cached to speed things up. + ElemType furthestDescendantDistance; + //! An instantiated metric. + MetricType metric; public: /** @@ -97,7 +114,7 @@ class Octree * @param data Dataset to create tree from. This will be copied! * @param maxLeafSize Maximum number of points in a leaf node. */ - Octree(const MatType& data, const size_t maxLeafSize = 20); + Octree(MatType&& data, const size_t maxLeafSize = 20); /** * Construct this as the root node of an octree on the given dataset. This @@ -183,6 +200,199 @@ class Octree const double width, const size_t maxLeafSize = 20); + /** + * Copy the given tree. Be careful! This may use a lot of memory. + * + * @param other Tree to copy from. + */ + Octree(const Octree& other); + + /** + * Move the given tree. The tree passed as a parameter will be emptied and + * will not be usable after this call. + * + * @param other Tree to move. + */ + Octree(Octree&& other); + + /** + * Initialize the tree from a boost::serialization archive. + * + * @param ar Archive to load tree from. Must be an iarchive, not an oarchive. + */ + template + Octree( + Archive& ar, + const typename boost::enable_if::type* = 0); + + /** + * Destroy the tree. + */ + ~Octree(); + + //! Return the dataset used by this node. + const MatType& Dataset() const { return *dataset; } + + //! Get the pointer to the parent. + Octree* Parent() const { return parent; } + //! Modify the pointer to the parent (be careful!). + Octree*& Parent() { return parent; } + + //! Return the bound object for this node. + const bound::HRectBound& Bound() const { return bound; } + //! Modify the bound object for this node. + bound::HRectBound& Bound() { return bound; } + + //! Return the statistic object for this node. + const StatisticType& Stat() const { return stat; } + //! Modify the statistic object for this node. + StatisticType& Stat() { return stat; } + + //! Return the number of children in this node. + size_t NumChildren() const; + + //! Return the metric that this tree uses. + MetricType Metric() const { return MetricType(); } + + /** + * Return the index of the nearest child node to the given query point. If + * this is a leaf node, it will return NumChildren() (invalid index). + */ + template + size_t GetNearestChild( + const VecType& point, + typename boost::enable_if>::type* = 0) const; + + /** + * Return the index of the furthest child node to the given query point. If + * this is a leaf node, it will return NumChildren() (invalid index). + */ + template + size_t GetFurthestChild( + const VecType& point, + typename boost::enable_if >::type* = 0) const; + + /** + * Return whether or not the node is a leaf. + */ + bool IsLeaf() const { return NumChildren() == 0; } + + /** + * Return the index of the nearest child node to the given query node. If it + * can't decide, it will return NumChildren() (invalid index). + */ + size_t GetNearestChild(const Octree& queryNode) const; + + /** + * Return the index of the furthest child node to the given query node. If it + * can't decide, it will return NumChildren() (invalid index). + */ + size_t GetFurthestChild(const Octree& queryNode) const; + + /** + * Return the furthest distance to a point held in this node. If this is not + * a leaf node, then the distance is 0 because the node holds no points. + */ + ElemType FurthestPointDistance() const; + + /** + * Return the furthest possible descendant distance. This returns the maximum + * distance from the centroid to the edge of the bound and not the empirical + * quantity which is the actual furthest descendant distance. So the actual + * furthest descendant distance may be less than what this method returns (but + * it will never be greater than this). + */ + ElemType FurthestDescendantDistance() const; + + //! Return the minimum distance from the center of the node to any bound edge. + ElemType MinimumBoundDistance() const; + + //! Return the distance from the center of this node to the center of the + //! parent node. + ElemType ParentDistance() const { return parentDistance; } + //! Modify the distance from the center of this node to the center of the + //! parent node. + ElemType& ParentDistance() { return parentDistance; } + + /** + * Return the specified child. If the index is out of bounds, unspecified + * behavior will occur. + */ + const Octree& Child(const size_t child) const { return *children[child]; } + + /** + * Return the specified child. If the index is out of bounds, unspecified + * behavior will occur. + */ + Octree& Child(const size_t child) { return *children[child]; } + + /** + * Return the pointer to the given child. This allows the child itself to be + * modified. + */ + Octree*& ChildPtr(const size_t child) { return children[child]; } + + //! Return the number of points in this node (0 if not a leaf). + size_t NumPoints() const; + + //! Return the number of descendants of this node. + size_t NumDescendants() const; + + /** + * Return the index (with reference to the dataset) of a particular + * descendant. + */ + size_t Descendant(const size_t index) const; + + /** + * Return the index (with reference to the dataset) of a particular point in + * this node. If the given index is invalid (i.e. if it is greater than + * NumPoints()), the indices returned will be invalid. + */ + size_t Point(const size_t index) const; + + //! Return the minimum distance to another node. + ElemType MinDistance(const Octree* other) const; + //! Return the maximum distance to another node. + ElemType MaxDistance(const Octree* other) const; + //! Return the minimum and maximum distance to another node. + math::RangeType RangeDistance(const Octree* other) const; + + //! Return the minimum distance to the given point. + template + ElemType MinDistance( + const VecType& point, + typename boost::enable_if>::type* = 0) const; + //! Return the maximum distance to the given point. + template + ElemType MaxDistance( + const VecType& point, + typename boost::enable_if>::type* = 0) const; + //! Return the minimum and maximum distance to another node. + template + math::RangeType RangeDistance( + const VecType& point, + typename boost::enable_if>::type* = 0) const; + + //! Store the center of the bounding region in the given vector. + void Center(arma::vec& center) const { bound.Center(center); } + + //! Serialize the tree. + template + void Serialize(Archive& ar, const unsigned int /* version */); + + protected: + /** + * A default constructor. This is meant to only be used with + * boost::serialization, which is allowed with the friend declaration below. + * This does not return a valid treee! The method must be protected, so that + * the serialization shim can work with the default constructor. + */ + Octree(); + + //! Friend access is given for the default constructor. + friend class boost::serialization::access; + private: /** * Split the node, using the given center and the given maximum width of this @@ -190,8 +400,11 @@ class Octree * * @param center Center of the node. * @param width Width of the current node. + * @param maxLeafSize Maximum number of points allowed in a leaf. */ - void SplitNode(const arma::vec& center, const double width); + void SplitNode(const arma::vec& center, + const double width, + const size_t maxLeafSize); /** * Split the node, using the given center and the given maximum width of this @@ -200,8 +413,18 @@ class Octree * @param center Center of the node. * @param width Width of the current node. * @param oldFromNew Mappings from old to new. + * @param maxLeafSize Maximum number of points allowed in a leaf. */ void SplitNode(const arma::vec& center, const double width, - std::vector& oldFromNew); + std::vector& oldFromNew, + const size_t maxLeafSize); }; + +} // namespace tree +} // namespace mlpack + +// Include implementation. +#include "octree_impl.hpp" + +#endif diff --git a/src/mlpack/core/tree/octree/octree_impl.hpp b/src/mlpack/core/tree/octree/octree_impl.hpp index cd8234ace3..987d139de5 100644 --- a/src/mlpack/core/tree/octree/octree_impl.hpp +++ b/src/mlpack/core/tree/octree/octree_impl.hpp @@ -9,19 +9,40 @@ #include "octree.hpp" +namespace mlpack { +namespace tree { + //! Construct the tree. template Octree::Octree(const MatType& dataset, - const double maxLeafSize) : + const size_t maxLeafSize) : + begin(0), + count(dataset.n_cols), + bound(dataset.n_rows), dataset(new MatType(dataset)), - + parent(NULL), + parentDistance(0.0) { - // Calculate empirical center of data. - bound |= *dataset; - arma::vec center = bound.Center(); - double maxWidth = bound.MaxWidth(); + if (count > 0) + { + // Calculate empirical center of data. + bound |= *this->dataset; + arma::vec center; + bound.Center(center); - SplitNode(center, maxWidth); + double maxWidth = 0.0; + for (size_t i = 0; i < bound.Dim(); ++i) + if (bound[i].Hi() - bound[i].Lo() > maxWidth) + maxWidth = bound[i].Hi() - bound[i].Lo(); + + SplitNode(center, maxWidth, maxLeafSize); + + furthestDescendantDistance = 0.5 * bound.Diameter(); + } + else + { + furthestDescendantDistance = 0.0; + } // Initialize the statistic. stat = StatisticType(*this); @@ -33,19 +54,37 @@ Octree::Octree( const MatType& dataset, std::vector& oldFromNew, const size_t maxLeafSize) : + begin(0), + count(dataset.n_cols), + bound(dataset.n_rows), dataset(new MatType(dataset)), - + parent(NULL), + parentDistance(0.0) { - // Calculate empirical center of data. - bound |= *dataset; - arma::vec center = bound.Center(); - double maxWidth = bound.MaxWidth(); - - oldFromNew.resize(data.n_cols); - for (size_t i = 0; i < data.n_cols; ++i) + oldFromNew.resize(this->dataset->n_cols); + for (size_t i = 0; i < this->dataset->n_cols; ++i) oldFromNew[i] = i; - SplitNode(center, maxWidth, oldFromNew); + if (count > 0) + { + // Calculate empirical center of data. + bound |= *this->dataset; + arma::vec center; + bound.Center(center); + + double maxWidth = 0.0; + for (size_t i = 0; i < bound.Dim(); ++i) + if (bound[i].Hi() - bound[i].Lo() > maxWidth) + maxWidth = bound[i].Hi() - bound[i].Lo(); + + SplitNode(center, maxWidth, oldFromNew, maxLeafSize); + + furthestDescendantDistance = 0.5 * bound.Diameter(); + } + else + { + furthestDescendantDistance = 0.0; + } // Initialize the statistic. stat = StatisticType(*this); @@ -56,32 +95,621 @@ template Octree::Octree( const MatType& dataset, std::vector& oldFromNew, + std::vector& newFromOld, const size_t maxLeafSize) : + begin(0), + count(dataset.n_cols), + bound(dataset.n_rows), dataset(new MatType(dataset)), - + parent(NULL), + parentDistance(0.0) { - // Calculate empirical center of data. - bound |= *dataset; - arma::vec center = bound.Center(); - double maxWidth = bound.MaxWidth(); - - oldFromNew.resize(data.n_cols); - for (size_t i = 0; i < data.n_cols; ++i) + oldFromNew.resize(this->dataset->n_cols); + for (size_t i = 0; i < this->dataset->n_cols; ++i) oldFromNew[i] = i; - SplitNode(center, maxWidth, oldFromNew); + if (count > 0) + { + // Calculate empirical center of data. + bound |= *this->dataset; + arma::vec center; + bound.Center(center); + + double maxWidth = 0.0; + for (size_t i = 0; i < bound.Dim(); ++i) + if (bound[i].Hi() - bound[i].Lo() > maxWidth) + maxWidth = bound[i].Hi() - bound[i].Lo(); + + SplitNode(center, maxWidth, oldFromNew, maxLeafSize); + + furthestDescendantDistance = 0.5 * bound.Diameter(); + } + else + { + furthestDescendantDistance = 0.0; + } + + // Initialize the statistic. + stat = StatisticType(*this); + + // Map the newFromOld indices correctly. + newFromOld.resize(this->dataset->n_cols); + for (size_t i = 0; i < this->dataset->n_cols; i++) + newFromOld[oldFromNew[i]] = i; +} + +//! Construct the tree. +template +Octree::Octree(MatType&& dataset, + const size_t maxLeafSize) : + begin(0), + count(dataset.n_cols), + bound(dataset.n_rows), + dataset(new MatType(std::move(dataset))), + parent(NULL), + parentDistance(0.0) +{ + if (count > 0) + { + // Calculate empirical center of data. + bound |= *this->dataset; + arma::vec center; + bound.Center(center); + + double maxWidth = 0.0; + for (size_t i = 0; i < bound.Dim(); ++i) + if (bound[i].Hi() - bound[i].Lo() > maxWidth) + maxWidth = bound[i].Hi() - bound[i].Lo(); + + SplitNode(center, maxWidth, maxLeafSize); + + furthestDescendantDistance = 0.5 * bound.Diameter(); + } + else + { + furthestDescendantDistance = 0.0; + } // Initialize the statistic. stat = StatisticType(*this); } +//! Construct the tree. +template +Octree::Octree( + MatType&& dataset, + std::vector& oldFromNew, + const size_t maxLeafSize) : + begin(0), + count(dataset.n_cols), + bound(dataset.n_rows), + dataset(new MatType(std::move(dataset))), + parent(NULL), + parentDistance(0.0) +{ + oldFromNew.resize(this->dataset->n_cols); + for (size_t i = 0; i < this->dataset->n_cols; ++i) + oldFromNew[i] = i; + + if (count > 0) + { + // Calculate empirical center of data. + bound |= *this->dataset; + arma::vec center; + bound.Center(center); + + double maxWidth = 0.0; + for (size_t i = 0; i < bound.Dim(); ++i) + if (bound[i].Hi() - bound[i].Lo() > maxWidth) + maxWidth = bound[i].Hi() - bound[i].Lo(); + + SplitNode(center, maxWidth, oldFromNew, maxLeafSize); + + furthestDescendantDistance = 0.5 * bound.Diameter(); + } + else + { + furthestDescendantDistance = 0.0; + } + + // Initialize the statistic. + stat = StatisticType(*this); +} + +//! Construct the tree. +template +Octree::Octree( + MatType&& dataset, + std::vector& oldFromNew, + std::vector& newFromOld, + const size_t maxLeafSize) : + begin(0), + count(dataset.n_cols), + bound(dataset.n_rows), + dataset(new MatType(std::move(dataset))), + parent(NULL), + parentDistance(0.0) +{ + oldFromNew.resize(this->dataset->n_cols); + for (size_t i = 0; i < this->dataset->n_cols; ++i) + oldFromNew[i] = i; + + if (count > 0) + { + // Calculate empirical center of data. + bound |= *this->dataset; + arma::vec center; + bound.Center(center); + + double maxWidth = 0.0; + for (size_t i = 0; i < bound.Dim(); ++i) + if (bound[i].Hi() - bound[i].Lo() > maxWidth) + maxWidth = bound[i].Hi() - bound[i].Lo(); + + SplitNode(center, maxWidth, oldFromNew, maxLeafSize); + + furthestDescendantDistance = 0.5 * bound.Diameter(); + } + else + { + furthestDescendantDistance = 0.0; + } + + // Initialize the statistic. + stat = StatisticType(*this); + + // Map the newFromOld indices correctly. + newFromOld.resize(this->dataset->n_cols); + for (size_t i = 0; i < this->dataset->n_cols; i++) + newFromOld[oldFromNew[i]] = i; +} + +//! Construct a child node. +template +Octree::Octree( + Octree* parent, + const size_t begin, + const size_t count, + const arma::vec& center, + const double width, + const size_t maxLeafSize) : + begin(begin), + count(count), + bound(parent->dataset->n_rows), + dataset(parent->dataset), + parent(parent) +{ + // Calculate empirical center of data. + bound |= dataset->cols(begin, begin + count - 1); + + // Now split the node. + SplitNode(center, width, maxLeafSize); + + // Calculate the distance from the empirical center of this node to the + // empirical center of the parent. + arma::vec trueCenter, parentCenter; + bound.Center(trueCenter); + parent->Bound().Center(parentCenter); + parentDistance = metric.Evaluate(trueCenter, parentCenter); + + furthestDescendantDistance = 0.5 * bound.Diameter(); + + // Initialize the statistic. + stat = StatisticType(*this); +} + +//! Construct a child node. +template +Octree::Octree( + Octree* parent, + const size_t begin, + const size_t count, + std::vector& oldFromNew, + const arma::vec& center, + const double width, + const size_t maxLeafSize) : + begin(begin), + count(count), + bound(parent->dataset->n_rows), + dataset(parent->dataset), + parent(parent) +{ + // Calculate empirical center of data. + bound |= dataset->cols(begin, begin + count - 1); + + // Now split the node. + SplitNode(center, width, oldFromNew, maxLeafSize); + + // Calculate the distance from the empirical center of this node to the + // empirical center of the parent. + arma::vec trueCenter, parentCenter; + bound.Center(trueCenter); + parent->Bound().Center(parentCenter); + parentDistance = metric.Evaluate(trueCenter, parentCenter); + + furthestDescendantDistance = 0.5 * bound.Diameter(); + + // Initialize the statistic. + stat = StatisticType(*this); +} + +//! Copy the given tree. +template +Octree::Octree(const Octree& other) : + begin(other.begin), + count(other.count), + bound(other.bound), + dataset((other.parent == NULL) ? new MatType(*other.dataset) : NULL), + parent(NULL), + stat(other.stat), + parentDistance(other.parentDistance), + furthestDescendantDistance(other.furthestDescendantDistance), + metric(other.metric) +{ + // If we have any children, we need to create them, and then ensure that their + // parent links are set right. + for (size_t i = 0; i < other.NumChildren(); ++i) + { + children.push_back(new Octree(other.Child(i))); + children[i]->parent = this; + children[i]->dataset = this->dataset; + } +} + +//! Move the given tree. +template +Octree::Octree(Octree&& other) : + children(std::move(other.children)), + begin(other.begin), + count(other.count), + bound(std::move(other.bound)), + dataset(other.dataset), + parent(other.parent), + stat(std::move(other.stat)), + parentDistance(other.parentDistance), + furthestDescendantDistance(other.furthestDescendantDistance), + metric(std::move(other.metric)) +{ + // Update the parent pointers of the direct children. + for (size_t i = 0; i < children.size(); ++i) + children[i]->parent = this; + + other.begin = 0; + other.count = 0; + other.dataset = new MatType(); + other.parentDistance = 0.0; + other.furthestDescendantDistance = 0.0; + other.parent = NULL; +} + +template +Octree::Octree() : + begin(0), + count(0), + bound(0), + dataset(new MatType()), + parent(NULL), + parentDistance(0.0), + furthestDescendantDistance(0.0) +{ + // Nothing to do. +} + +template +template +Octree::Octree( + Archive& ar, + const typename boost::enable_if::type*) : + Octree() // Create an empty tree. +{ + // De-serialize the tree into this object. + ar >> data::CreateNVP(*this, "tree"); +} + +template +Octree::~Octree() +{ + // Delete the dataset if we aren't the parent. + if (!parent) + delete dataset; + + // Now delete each of the children. + for (size_t i = 0; i < children.size(); ++i) + delete children[i]; + children.clear(); +} + +template +size_t Octree::NumChildren() const +{ + return children.size(); +} + +template +template +size_t Octree::GetNearestChild( + const VecType& point, + typename boost::enable_if>::type*) const +{ + // It's possible that this could be improved by caching which children we have + // and which we don't, but for now this is just a brute force search. + ElemType bestDistance = DBL_MAX; + size_t bestIndex = NumChildren(); + for (size_t i = 0; i < NumChildren(); ++i) + { + const double dist = children[i]->MinDistance(point); + if (dist < bestDistance) + { + bestDistance = dist; + bestIndex = i; + } + } + + return bestIndex; +} + +template +template +size_t Octree::GetFurthestChild( + const VecType& point, + typename boost::enable_if>::type*) const +{ + // It's possible that this could be improved by caching which children we have + // and which we don't, but for now this is just a brute force search. + ElemType bestDistance = -1.0; // Initialize to invalid distance. + size_t bestIndex = NumChildren(); + for (size_t i = 0; i < NumChildren(); ++i) + { + const double dist = children[i]->MaxDistance(point); + if (dist > bestDistance) + { + bestDistance = dist; + bestIndex = i; + } + } + + return bestIndex; +} + +template +size_t Octree::GetNearestChild( + const Octree& queryNode) const +{ + // It's possible that this could be improved by caching which children we have + // and which we don't, but for now this is just a brute force search. + ElemType bestDistance = DBL_MAX; + size_t bestIndex = NumChildren(); + for (size_t i = 0; i < NumChildren(); ++i) + { + const double dist = children[i]->MaxDistance(queryNode); + if (dist < bestDistance) + { + bestDistance = dist; + bestIndex = i; + } + } + + return bestIndex; +} + +template +size_t Octree::GetFurthestChild( + const Octree& queryNode) const +{ + // It's possible that this could be improved by caching which children we have + // and which we don't, but for now this is just a brute force search. + ElemType bestDistance = -1.0; // Initialize to invalid distance. + size_t bestIndex = NumChildren(); + for (size_t i = 0; i < NumChildren(); ++i) + { + const double dist = children[i]->MaxDistance(queryNode); + if (dist > bestDistance) + { + bestDistance = dist; + bestIndex = i; + } + } + + return bestIndex; +} + +template +typename Octree::ElemType +Octree::FurthestPointDistance() + const +{ + // If we are not a leaf, then this distance is 0. Otherwise, return the + // furthest descendant distance. + return (children.size() > 0) ? 0.0 : furthestDescendantDistance; +} + +template +typename Octree::ElemType +Octree::FurthestDescendantDistance() const +{ + return furthestDescendantDistance; +} + +template +typename Octree::ElemType +Octree::MinimumBoundDistance() const +{ + return bound.MinWidth() / 2.0; +} + +template +size_t Octree::NumPoints() const +{ + // We have no points unless we are a leaf; + return (children.size() > 0) ? 0 : count; +} + +template +size_t Octree::NumDescendants() const +{ + return count; +} + +template +size_t Octree::Descendant( + const size_t index) const +{ + return begin + index; +} + +template +size_t Octree::Point(const size_t index) + const +{ + return begin + index; +} + +template +typename Octree::ElemType +Octree::MinDistance(const Octree* other) + const +{ + return bound.MinDistance(other->Bound()); +} + +template +typename Octree::ElemType +Octree::MaxDistance(const Octree* other) + const +{ + return bound.MaxDistance(other->Bound()); +} + +template +math::RangeType::ElemType> +Octree::RangeDistance(const Octree* other) + const +{ + return bound.RangeDistance(other->Bound()); +} + +template +template +typename Octree::ElemType +Octree::MinDistance( + const VecType& point, + typename boost::enable_if>::type*) const +{ + return bound.MinDistance(point); +} + +template +template +typename Octree::ElemType +Octree::MaxDistance( + const VecType& point, + typename boost::enable_if>::type*) const +{ + return bound.MaxDistance(point); +} + + +template +template +math::RangeType::ElemType> +Octree::RangeDistance( + const VecType& point, + typename boost::enable_if>::type*) const +{ + return bound.RangeDistance(point); +} + +//! Serialize the tree. +template +template +void Octree::Serialize( + Archive& ar, + const unsigned int /* version */) +{ + using data::CreateNVP; + + // + + // If we're loading and we have children, they need to be deleted. + if (Archive::is_loading::value) + { + for (size_t i = 0; i < children.size(); ++i) + delete children[i]; + children.clear(); + + if (!parent) + delete dataset; + } + + ar & CreateNVP(begin, "begin"); + ar & CreateNVP(count, "count"); + ar & CreateNVP(bound, "bound"); + ar & CreateNVP(stat, "stat"); + ar & CreateNVP(parentDistance, "parentDistance"); + ar & CreateNVP(furthestDescendantDistance, "furthestDescendantDistance"); + ar & CreateNVP(metric, "metric"); + + // Due to quirks of boost::serialization, depending on how the user + // serializes the tree, it's possible that the root of the tree will + // accidentally be serialized twice. So if we are a first-level child, we + // avoid serializing the parent. The true (non-duplicated) parent will fix + // the parent link. + bool hasFakeParent = false; + if (Archive::is_saving::value && parent != NULL && parent->parent == NULL) + { + Octree* fakeParent = NULL; + hasFakeParent = true; + ar & CreateNVP(fakeParent, "parent"); + ar & CreateNVP(hasFakeParent, "hasFakeParent"); + } + else + { + ar & CreateNVP(parent, "parent"); + ar & CreateNVP(hasFakeParent, "hasFakeParent"); + } + + // Only serialize the dataset if we don't have a fake parent. Otherwise, the + // real parent will come and set it later. + if (!hasFakeParent) + ar & CreateNVP(dataset, "dataset"); + + size_t numChildren = 0; + if (Archive::is_saving::value) + numChildren = children.size(); + ar & CreateNVP(numChildren, "numChildren"); + if (Archive::is_loading::value) + children.resize(numChildren); + + for (size_t i = 0; i < numChildren; ++i) + { + std::ostringstream oss; + oss << "child" << i; + ar & CreateNVP(children[i], oss.str()); + } + + // Fix the child pointers, if they were set to a fake parent. + if (Archive::is_loading::value && parent == NULL) + { + for (size_t i = 0; i < children.size(); ++i) + { + children[i]->dataset = this->dataset; + children[i]->parent = this; + } + } +} //! Split the node. template void Octree::SplitNode( const arma::vec& center, - const double width) + const double width, + const size_t maxLeafSize) { + // No need to split if we have fewer than the maximum number of points in this + // node. + if (count <= maxLeafSize) + return; + // We must split the dataset by sequentially creating each of the children. // We do this in two steps: first we make a pass to count the number of points // that will fall into each child; then in the second pass we rearrange the @@ -101,7 +729,7 @@ void Octree::SplitNode( // the points fall on. The last dimension represents the most significant // bit in the assignment; the bit is '1' if it falls to the right of the // center. - if (dataset(d, begin + i) > center(d)) + if ((*dataset)(d, begin + i) > center(d)) assignments(i) |= (1 << d); } @@ -130,7 +758,7 @@ void Octree::SplitNode( for (size_t d = 0; d < center.n_elem; ++d) { // Is the dimension "right" (1) or "left" (0)? - if ((i >> d) & 1 == 0) + if (((i >> d) & 1) == 0) childCenter[d] = center[d] - childWidth; else childCenter[d] = center[d] + childWidth; @@ -148,8 +776,14 @@ template void Octree::SplitNode( const arma::vec& center, const double width, - std::vector& oldFromNew) + std::vector& oldFromNew, + const size_t maxLeafSize) { + // No need to split if we have fewer than the maximum number of points in this + // node. + if (count <= maxLeafSize) + return; + // We must split the dataset by sequentially creating each of the children. // We do this in two steps: first we make a pass to count the number of points // that will fall into each child; then in the second pass we rearrange the @@ -169,7 +803,7 @@ void Octree::SplitNode( // the points fall on. The last dimension represents the most significant // bit in the assignment; the bit is '1' if it falls to the right of the // center. - if (dataset(d, begin + i) > center(d)) + if ((*dataset)(d, begin + i) > center(d)) assignments(i) |= (1 << d); } @@ -183,8 +817,9 @@ void Octree::SplitNode( // really a problem. We use non-contiguous submatrix views to extract the // columns in the correct order. dataset->cols(begin, begin + count - 1) = dataset->cols(begin + ordering); + std::vector oldFromNewCopy(oldFromNew); // We need the old indices. for (size_t i = 0; i < count; ++i) - oldFromNew[ordering[i] + begin] = i + begin; + oldFromNew[i + begin] = oldFromNewCopy[ordering[i] + begin]; // Now that the dataset is reordered, we can create the children. size_t childBegin = begin; @@ -200,7 +835,7 @@ void Octree::SplitNode( for (size_t d = 0; d < center.n_elem; ++d) { // Is the dimension "right" (1) or "left" (0)? - if ((i >> d) & 1 == 0) + if (((i >> d) & 1) == 0) childCenter[d] = center[d] - childWidth; else childCenter[d] = center[d] + childWidth; @@ -212,3 +847,8 @@ void Octree::SplitNode( childBegin += childCounts[i]; } } + +} // namespace tree +} // namespace mlpack + +#endif diff --git a/src/mlpack/core/tree/octree/single_tree_traverser.hpp b/src/mlpack/core/tree/octree/single_tree_traverser.hpp new file mode 100644 index 0000000000..a82149f3fb --- /dev/null +++ b/src/mlpack/core/tree/octree/single_tree_traverser.hpp @@ -0,0 +1,53 @@ +/** + * @file single_tree_traverser.hpp + * @author Ryan Curtin + * + * Definition of the single tree traverser for the octree. + */ +#ifndef MLPACK_CORE_TREE_OCTREE_SINGLE_TREE_TRAVERSER_HPP +#define MLPACK_CORE_TREE_OCTREE_SINGLE_TREE_TRAVERSER_HPP + +#include +#include "octree.hpp" + +namespace mlpack { +namespace tree { + +template +template +class Octree::SingleTreeTraverser +{ + public: + /** + * Instantiate the traverser with the given rule set. + */ + SingleTreeTraverser(RuleType& rule); + + /** + * Traverse the reference tree with the given query point. This does not + * reset the number of pruned nodes. + * + * @param queryIndex Index of query point. + * @param referenceNode Node in reference tree. + */ + void Traverse(const size_t queryIndex, Octree& referenceNode); + + //! Get the number of pruned nodes. + size_t NumPrunes() const { return numPrunes; } + //! Modify the number of pruned nodes. + size_t& NumPrunes() { return numPrunes; } + + private: + //! The instantiated rule. + RuleType& rule; + //! The number of reference nodes that have been pruned. + size_t numPrunes; +}; + +} // namespace tree +} // namespace mlpack + +// Include implementation. +#include "single_tree_traverser_impl.hpp" + +#endif diff --git a/src/mlpack/core/tree/octree/single_tree_traverser_impl.hpp b/src/mlpack/core/tree/octree/single_tree_traverser_impl.hpp new file mode 100644 index 0000000000..49ba9c2d73 --- /dev/null +++ b/src/mlpack/core/tree/octree/single_tree_traverser_impl.hpp @@ -0,0 +1,67 @@ +/** + * @file single_tree_traverser_impl.hpp + * @author Ryan Curtin + * + * Implementation of the single tree traverser for octrees. + */ +#ifndef MLPACK_CORE_TREE_OCTREE_SINGLE_TREE_TRAVERSER_IMPL_HPP +#define MLPACK_CORE_TREE_OCTREE_SINGLE_TREE_TRAVERSER_IMPL_HPP + +// In case it hasn't been included yet. +#include "single_tree_traverser.hpp" + +namespace mlpack { +namespace tree { + +template +template +Octree::SingleTreeTraverser:: + SingleTreeTraverser(RuleType& rule) : + rule(rule) +{ + // Nothing to do. +} + +template +template +void Octree::SingleTreeTraverser:: + Traverse(const size_t queryIndex, Octree& referenceNode) +{ + // If we are a leaf, run the base cases. + if (referenceNode.NumChildren() == 0) + { + const size_t refBegin = referenceNode.Point(0); + const size_t refEnd = refBegin + referenceNode.NumPoints(); + for (size_t r = refBegin; r < refEnd; ++r) + rule.BaseCase(queryIndex, r); + } + else + { + // Do a prioritized recursion, by scoring all candidates and then sorting + // them. + arma::vec scores(referenceNode.NumChildren()); + for (size_t i = 0; i < scores.n_elem; ++i) + scores[i] = rule.Score(queryIndex, referenceNode.Child(i)); + + // Sort the scores. + arma::uvec sortedIndices = arma::sort_index(scores); + + for (size_t i = 0; i < sortedIndices.n_elem; ++i) + { + // If the node is pruned, all subsequent nodes in sorted order will also + // be pruned. + if (scores[sortedIndices[i]] == DBL_MAX) + { + numPrunes += (sortedIndices.n_elem - i); + break; + } + + Traverse(queryIndex, referenceNode.Child(sortedIndices[i])); + } + } +} + +} // namespace tree +} // namespace mlpack + +#endif diff --git a/src/mlpack/core/tree/octree/traits.hpp b/src/mlpack/core/tree/octree/traits.hpp new file mode 100644 index 0000000000..4a0f9564c9 --- /dev/null +++ b/src/mlpack/core/tree/octree/traits.hpp @@ -0,0 +1,66 @@ +/** + * @file traits.hpp + * @author Ryan Curtin + * + * Specialization of the TreeTraits class for the Octree class. + */ +#ifndef MLPACK_CORE_TREE_OCTREE_TRAITS_HPP +#define MLPACK_CORE_TREE_OCTREE_TRAITS_HPP + +#include + +namespace mlpack { +namespace tree { + +/** + * This is a specialization of the TreeTraits class to the Octree tree type. It + * defines characteristics of the octree, and is used to help write + * tree-independent (but still optimized) tree-based algorithms. See + * mlpack/core/tree/tree_traits.hpp for more information. + */ +template +class TreeTraits> +{ + public: + /** + * No octree nodes will overlap. + */ + static const bool HasOverlappingChildren = false; + + /** + * Points are not shared across nodes in the octree. + */ + static const bool HasDuplicatedPoints = false; + + /** + * There is no guarantee that the first point in a node is its centroid. + */ + static const bool FirstPointIsCentroid = false; + + /** + * Points are not contained at multiple levels of the octree. + */ + static const bool HasSelfChildren = false; + + /** + * Points are rearranged during building of the tree. + */ + static const bool RearrangesDataset = true; + + /** + * This is not necessarily a binary tree. + */ + static const bool BinaryTree = false; + + /** + * NumDescendants() represents the number of unique descendant points. + */ + static const bool UniqueNumDescendants = true; +}; + +} // namespace tree +} // namespace mlpack + +#endif diff --git a/src/mlpack/methods/neighbor_search/kfn_main.cpp b/src/mlpack/methods/neighbor_search/kfn_main.cpp index 2eb5bd6ee9..96454721e4 100644 --- a/src/mlpack/methods/neighbor_search/kfn_main.cpp +++ b/src/mlpack/methods/neighbor_search/kfn_main.cpp @@ -68,10 +68,10 @@ PARAM_INT_IN("k", "Number of furthest neighbors to find.", "k", 0); // building. PARAM_STRING_IN("tree_type", "Type of tree to use: 'kd', 'vp', 'rp', 'max-rp', " "'ub', 'cover', 'r', 'r-star', 'x', 'ball', 'hilbert-r', 'r-plus', " - "'r-plus-plus'.", "t", "kd"); + "'r-plus-plus', 'octree'.", "t", "kd"); PARAM_INT_IN("leaf_size", "Leaf size for tree building (used for kd-trees, " "vp trees, random projection trees, UB trees, R trees, R* trees, X trees, " - "Hilbert R trees, R+ trees and R++ trees).", "l", 20); + "Hilbert R trees, R+ trees, R++ trees, and octrees).", "l", 20); PARAM_FLAG("random_basis", "Before tree-building, project the data onto a " "random orthogonal basis.", "R"); PARAM_INT_IN("seed", "Random seed (if 0, std::time(NULL) is used).", "s", 0); @@ -262,10 +262,13 @@ int main(int argc, char *argv[]) tree = KFNModel::MAX_RP_TREE; else if (treeType == "ub") tree = KFNModel::UB_TREE; + else if (treeType == "octree") + tree = KFNModel::OCTREE; else Log::Fatal << "Unknown tree type '" << treeType << "'; valid choices are " << "'kd', 'vp', 'rp', 'max-rp', 'ub', 'cover', 'r', 'r-star', 'x', " - << "'ball', 'hilbert-r', 'r-plus' and 'r-plus-plus'." << endl; + << "'ball', 'hilbert-r', 'r-plus', 'r-plus-plus', and 'octree'." + << endl; kfn.TreeType() = tree; kfn.RandomBasis() = randomBasis; diff --git a/src/mlpack/methods/neighbor_search/knn_main.cpp b/src/mlpack/methods/neighbor_search/knn_main.cpp index 0e29ff838c..2127b79583 100644 --- a/src/mlpack/methods/neighbor_search/knn_main.cpp +++ b/src/mlpack/methods/neighbor_search/knn_main.cpp @@ -73,7 +73,8 @@ PARAM_STRING_IN("tree_type", "Type of tree to use: 'kd', 'vp', 'rp', 'max-rp', " "'r-plus-plus', 'spill'.", "t", "kd"); PARAM_INT_IN("leaf_size", "Leaf size for tree building (used for kd-trees, vp " "trees, random projection trees, UB trees, R trees, R* trees, X trees, " - "Hilbert R trees, R+ trees, R++ trees and spill trees).", "l", 20); + "Hilbert R trees, R+ trees, R++ trees, spill trees, and octrees).", "l", + 20); PARAM_DOUBLE_IN("tau", "Overlapping size (only valid for spill trees).", "u", 0); PARAM_DOUBLE_IN("rho", "Balance threshold (only valid for spill trees).", "b", @@ -276,11 +277,13 @@ int main(int argc, char *argv[]) tree = KNNModel::MAX_RP_TREE; else if (treeType == "ub") tree = KNNModel::UB_TREE; + else if (treeType == "octree") + tree = KNNModel::OCTREE; else Log::Fatal << "Unknown tree type '" << treeType << "'; valid choices are " << "'kd', 'vp', 'rp', 'max-rp', 'ub', 'cover', 'r', 'r-star', 'x', " - << "'ball', 'hilbert-r', 'r-plus', 'r-plus-plus' and 'spill'." - << endl; + << "'ball', 'hilbert-r', 'r-plus', 'r-plus-plus', 'spill', and " + << "'octree'." << endl; knn.TreeType() = tree; knn.RandomBasis() = randomBasis; diff --git a/src/mlpack/methods/neighbor_search/ns_model.hpp b/src/mlpack/methods/neighbor_search/ns_model.hpp index 4f4c47d6e1..b0ae690da5 100644 --- a/src/mlpack/methods/neighbor_search/ns_model.hpp +++ b/src/mlpack/methods/neighbor_search/ns_model.hpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include "neighbor_search.hpp" @@ -133,6 +134,9 @@ class BiSearchVisitor : public boost::static_visitor //! Bichromatic neighbor search specialized for SPTrees. void operator()(SpillKNN* ns) const; + //! Bichromatic neighbor search specialized for octrees. + void operator()(NSTypeT* ns) const; + //! Construct the BiSearchVisitor. BiSearchVisitor(const arma::mat& querySet, const size_t k, @@ -188,6 +192,9 @@ class TrainVisitor : public boost::static_visitor //! Train specialized for SPTrees. void operator()(SpillKNN* ns) const; + //! Train specialized for octrees. + void operator()(NSTypeT* ns) const; + //! Construct the TrainVisitor object with the given reference set, leafSize //! for BinarySpaceTrees, and tau and rho for spill trees. TrainVisitor(arma::mat&& referenceSet, @@ -287,7 +294,8 @@ class NSModel RP_TREE, MAX_RP_TREE, SPILL_TREE, - UB_TREE + UB_TREE, + OCTREE }; private: @@ -325,7 +333,8 @@ class NSModel NSType*, NSType*, SpillKNN*, - NSType*> nSearch; + NSType*, + NSType*> nSearch; public: /** diff --git a/src/mlpack/methods/neighbor_search/ns_model_impl.hpp b/src/mlpack/methods/neighbor_search/ns_model_impl.hpp index 062a6ec4d7..c2d6c6b021 100644 --- a/src/mlpack/methods/neighbor_search/ns_model_impl.hpp +++ b/src/mlpack/methods/neighbor_search/ns_model_impl.hpp @@ -96,6 +96,15 @@ void BiSearchVisitor::operator()(SpillKNN* ns) const throw std::runtime_error("no neighbor search model initialized"); } +//! Bichromatic neighbor search specialized for octrees. +template +void BiSearchVisitor::operator()(NSTypeT* ns) const +{ + if (ns) + return SearchLeaf(ns); + throw std::runtime_error("no neighbor search model initialized"); +} + //! Bichromatic neighbor search on the given NSType considering the leafSize. template template @@ -150,7 +159,7 @@ void TrainVisitor::operator()(NSTypeT* ns) const //! Train on the given NSType specialized for KDTrees. template -void TrainVisitor::operator ()(NSTypeT* ns) const +void TrainVisitor::operator()(NSTypeT* ns) const { if (ns) return TrainLeaf(ns); @@ -159,7 +168,7 @@ void TrainVisitor::operator ()(NSTypeT* ns) const //! Train on the given NSType specialized for BallTrees. template -void TrainVisitor::operator ()(NSTypeT* ns) const +void TrainVisitor::operator()(NSTypeT* ns) const { if (ns) return TrainLeaf(ns); @@ -168,7 +177,7 @@ void TrainVisitor::operator ()(NSTypeT* ns) const //! Train specialized for SPTrees. template -void TrainVisitor::operator ()(SpillKNN* ns) const +void TrainVisitor::operator()(SpillKNN* ns) const { if (ns) { @@ -184,6 +193,15 @@ void TrainVisitor::operator ()(SpillKNN* ns) const throw std::runtime_error("no neighbor search model initialized"); } +//! Train specialized for Octrees. +template +void TrainVisitor::operator()(NSTypeT* ns) const +{ + if (ns) + return TrainLeaf(ns); + throw std::runtime_error("no neighbor search model initialized"); +} + //! Train on the given NSType considering the leafSize. template template @@ -485,6 +503,9 @@ void NSModel::BuildModel(arma::mat&& referenceSet, case UB_TREE: nSearch = new NSType(searchMode, epsilon); break; + case OCTREE: + nSearch = new NSType(searchMode, epsilon); + break; } TrainVisitor tn(std::move(referenceSet), leafSize, tau, rho); diff --git a/src/mlpack/methods/range_search/range_search_main.cpp b/src/mlpack/methods/range_search/range_search_main.cpp index 9487832b19..990e3c10a0 100644 --- a/src/mlpack/methods/range_search/range_search_main.cpp +++ b/src/mlpack/methods/range_search/range_search_main.cpp @@ -72,10 +72,10 @@ PARAM_DOUBLE_IN("min", "Lower bound in range.", "L", 0.0); // building. PARAM_STRING_IN("tree_type", "Type of tree to use: 'kd', 'vp', 'rp', 'max-rp', " "'ub', 'cover', 'r', 'r-star', 'x', 'ball', 'hilbert-r', 'r-plus', " - "'r-plus-plus'.", "t", "kd"); + "'r-plus-plus', 'octree'.", "t", "kd"); PARAM_INT_IN("leaf_size", "Leaf size for tree building (used for kd-trees, " "vp trees, random projection trees, UB trees, R trees, R* trees, X trees, " - "Hilbert R trees, R+ trees and R++ trees).", "l", 20); + "Hilbert R trees, R+ trees, R++ trees, and octrees).", "l", 20); PARAM_FLAG("random_basis", "Before tree-building, project the data onto a " "random orthogonal basis.", "R"); PARAM_INT_IN("seed", "Random seed (if 0, std::time(NULL) is used).", "s", 0); @@ -191,10 +191,12 @@ int main(int argc, char *argv[]) tree = RSModel::MAX_RP_TREE; else if (treeType == "ub") tree = RSModel::UB_TREE; + else if (treeType == "octree") + tree = RSModel::OCTREE; else Log::Fatal << "Unknown tree type '" << treeType << "; valid choices are " << "'kd', 'vp', 'rp', 'max-rp', 'ub', 'cover', 'r', 'r-star', 'x', " - << "'ball', 'hilbert-r', 'r-plus' and 'r-plus-plus'." << endl; + << "'ball', 'hilbert-r', 'r-plus', 'r-plus-plus', and 'octree'." << endl; rs.TreeType() = tree; rs.RandomBasis() = randomBasis; diff --git a/src/mlpack/methods/range_search/rs_model.cpp b/src/mlpack/methods/range_search/rs_model.cpp index 6025471734..140b53cf90 100644 --- a/src/mlpack/methods/range_search/rs_model.cpp +++ b/src/mlpack/methods/range_search/rs_model.cpp @@ -29,7 +29,8 @@ RSModel::RSModel(TreeTypes treeType, bool randomBasis) : vpTreeRS(NULL), rpTreeRS(NULL), maxRPTreeRS(NULL), - ubTreeRS(NULL) + ubTreeRS(NULL), + octreeRS(NULL) { // Nothing to do. } @@ -164,6 +165,28 @@ void RSModel::BuildModel(arma::mat&& referenceSet, ubTreeRS = new RSType(move(referenceSet), naive, singleMode); break; + + case OCTREE: + // If necessary, build the octree. + if (naive) + { + octreeRS = new RSType(move(referenceSet), naive, + singleMode); + } + else + { + vector oldFromNewReferences; + RSType::Tree* octree = + new RSType::Tree(move(referenceSet), + oldFromNewReferences, leafSize); + octreeRS = new RSType(octree, singleMode); + + // Give the model ownership of the tree and the mappings. + octreeRS->treeOwner = true; + octreeRS->oldFromNewReferences = move(oldFromNewReferences); + } + + break; } if (!naive) @@ -301,6 +324,38 @@ void RSModel::Search(arma::mat&& querySet, case UB_TREE: ubTreeRS->Search(querySet, range, neighbors, distances); break; + + case OCTREE: + if (!octreeRS->Naive() && !octreeRS->SingleMode()) + { + // Build a query tree and search. + Timer::Start("tree_building"); + Log::Info << "Building query tree..." << endl; + vector oldFromNewQueries; + RSType::Tree queryTree(move(querySet), oldFromNewQueries, + leafSize); + Log::Info << "Tree built." << endl; + Timer::Stop("tree_building"); + + vector> neighborsOut; + vector> distancesOut; + octreeRS->Search(&queryTree, range, neighborsOut, distancesOut); + + // Remap the query points. + neighbors.resize(queryTree.Dataset().n_cols); + distances.resize(queryTree.Dataset().n_cols); + for (size_t i = 0; i < queryTree.Dataset().n_cols; ++i) + { + neighbors[oldFromNewQueries[i]] = neighborsOut[i]; + distances[oldFromNewQueries[i]] = distancesOut[i]; + } + } + else + { + // Search without building a second tree. + octreeRS->Search(querySet, range, neighbors, distances); + } + break; } } @@ -371,6 +426,10 @@ void RSModel::Search(const math::Range& range, case UB_TREE: ubTreeRS->Search(range, neighbors, distances); break; + + case OCTREE: + octreeRS->Search(range, neighbors, distances); + break; } } @@ -405,6 +464,8 @@ std::string RSModel::TreeName() const return "random projection tree (max split)"; case UB_TREE: return "UB tree"; + case OCTREE: + return "octree"; default: return "unknown tree"; } @@ -413,32 +474,20 @@ std::string RSModel::TreeName() const // Clean memory. void RSModel::CleanMemory() { - if (kdTreeRS) - delete kdTreeRS; - if (coverTreeRS) - delete coverTreeRS; - if (rTreeRS) - delete rTreeRS; - if (rStarTreeRS) - delete rStarTreeRS; - if (ballTreeRS) - delete ballTreeRS; - if (xTreeRS) - delete xTreeRS; - if (hilbertRTreeRS) - delete hilbertRTreeRS; - if (rPlusTreeRS) - delete rPlusTreeRS; - if (rPlusPlusTreeRS) - delete rPlusPlusTreeRS; - if (vpTreeRS) - delete vpTreeRS; - if (rpTreeRS) - delete rpTreeRS; - if (maxRPTreeRS) - delete maxRPTreeRS; - if (ubTreeRS) - delete ubTreeRS; + delete kdTreeRS; + delete coverTreeRS; + delete rTreeRS; + delete rStarTreeRS; + delete ballTreeRS; + delete xTreeRS; + delete hilbertRTreeRS; + delete rPlusTreeRS; + delete rPlusPlusTreeRS; + delete vpTreeRS; + delete rpTreeRS; + delete maxRPTreeRS; + delete ubTreeRS; + delete octreeRS; kdTreeRS = NULL; coverTreeRS = NULL; @@ -453,4 +502,5 @@ void RSModel::CleanMemory() rpTreeRS = NULL; maxRPTreeRS = NULL; ubTreeRS = NULL; + octreeRS = NULL; } diff --git a/src/mlpack/methods/range_search/rs_model.hpp b/src/mlpack/methods/range_search/rs_model.hpp index 659319d4da..bea33d5cb3 100644 --- a/src/mlpack/methods/range_search/rs_model.hpp +++ b/src/mlpack/methods/range_search/rs_model.hpp @@ -13,6 +13,7 @@ #include #include #include +#include #include "range_search.hpp" @@ -36,7 +37,8 @@ class RSModel VP_TREE, RP_TREE, MAX_RP_TREE, - UB_TREE + UB_TREE, + OCTREE }; private: @@ -84,6 +86,8 @@ class RSModel //! Universal B tree based range search object //! (NULL if not in use). RSType* ubTreeRS; + //! Octree-based range search object (NULL if not in use). + RSType* octreeRS; public: /** diff --git a/src/mlpack/methods/range_search/rs_model_impl.hpp b/src/mlpack/methods/range_search/rs_model_impl.hpp index 2a8502fda5..69b7b7150d 100644 --- a/src/mlpack/methods/range_search/rs_model_impl.hpp +++ b/src/mlpack/methods/range_search/rs_model_impl.hpp @@ -81,6 +81,10 @@ void RSModel::Serialize(Archive& ar, const unsigned int /* version */) case UB_TREE: ar & CreateNVP(ubTreeRS, "range_search_model"); break; + + case OCTREE: + ar & CreateNVP(octreeRS, "range_search_model"); + break; } } @@ -112,6 +116,8 @@ inline const arma::mat& RSModel::Dataset() const return maxRPTreeRS->ReferenceSet(); else if (ubTreeRS) return ubTreeRS->ReferenceSet(); + else if (octreeRS) + return octreeRS->ReferenceSet(); throw std::runtime_error("no range search model initialized"); } @@ -144,6 +150,8 @@ inline bool RSModel::SingleMode() const return maxRPTreeRS->SingleMode(); else if (ubTreeRS) return ubTreeRS->SingleMode(); + else if (octreeRS) + return octreeRS->SingleMode(); throw std::runtime_error("no range search model initialized"); } @@ -176,6 +184,8 @@ inline bool& RSModel::SingleMode() return maxRPTreeRS->SingleMode(); else if (ubTreeRS) return ubTreeRS->SingleMode(); + else if (octreeRS) + return octreeRS->SingleMode(); throw std::runtime_error("no range search model initialized"); } @@ -208,6 +218,8 @@ inline bool RSModel::Naive() const return maxRPTreeRS->Naive(); else if (ubTreeRS) return ubTreeRS->Naive(); + else if (octreeRS) + return octreeRS->Naive(); throw std::runtime_error("no range search model initialized"); } @@ -240,6 +252,8 @@ inline bool& RSModel::Naive() return maxRPTreeRS->Naive(); else if (ubTreeRS) return ubTreeRS->Naive(); + else if (octreeRS) + return octreeRS->Naive(); throw std::runtime_error("no range search model initialized"); } diff --git a/src/mlpack/methods/rann/krann_main.cpp b/src/mlpack/methods/rann/krann_main.cpp index 591d741f95..7830f89299 100644 --- a/src/mlpack/methods/rann/krann_main.cpp +++ b/src/mlpack/methods/rann/krann_main.cpp @@ -65,10 +65,11 @@ PARAM_INT_IN("k", "Number of nearest neighbors to find.", "k", 0); // The user may specify the type of tree to use, and a few parameters for tree // building. PARAM_STRING_IN("tree_type", "Type of tree to use: 'kd', 'ub', 'cover', 'r', " - "'x', 'r-star', 'hilbert-r', 'r-plus', 'r-plus-plus'.", "t", "kd"); + "'x', 'r-star', 'hilbert-r', 'r-plus', 'r-plus-plus', 'octree'.", "t", + "kd"); PARAM_INT_IN("leaf_size", "Leaf size for tree building (used for kd-trees, " - "UB trees, R trees, R* trees, X trees, Hilbert R trees, R+ trees and " - "R++ trees).", "l", 20); + "UB trees, R trees, R* trees, X trees, Hilbert R trees, R+ trees, " + "R++ trees, and octrees).", "l", 20); PARAM_FLAG("random_basis", "Before tree-building, project the data onto a " "random orthogonal basis.", "R"); PARAM_INT_IN("seed", "Random seed (if 0, std::time(NULL) is used).", "s", 0); @@ -182,10 +183,12 @@ int main(int argc, char *argv[]) tree = RANNModel::R_PLUS_PLUS_TREE; else if (treeType == "ub") tree = RANNModel::UB_TREE; + else if (treeType == "octree") + tree = RANNModel::OCTREE; else Log::Fatal << "Unknown tree type '" << treeType << "'; valid choices are " << "'kd', 'ub', 'cover', 'r', 'r-star', 'x', 'hilbert-r', " - << "'r-plus' and 'r-plus-plus'." << endl; + << "'r-plus', 'r-plus-plus', 'octree'." << endl; rann.TreeType() = tree; rann.RandomBasis() = randomBasis; diff --git a/src/mlpack/methods/rann/ra_model.hpp b/src/mlpack/methods/rann/ra_model.hpp index 48be9e915d..676a555de1 100644 --- a/src/mlpack/methods/rann/ra_model.hpp +++ b/src/mlpack/methods/rann/ra_model.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #include "ra_search.hpp" @@ -44,7 +45,8 @@ class RAModel HILBERT_R_TREE, R_PLUS_TREE, R_PLUS_PLUS_TREE, - UB_TREE + UB_TREE, + OCTREE }; private: @@ -85,6 +87,8 @@ class RAModel RAType* rPlusPlusTreeRA; //! Non-NULL if the UB tree is used. RAType* ubTreeRA; + //! Non-NULL if the octree is used. + RAType* octreeRA; public: /** diff --git a/src/mlpack/methods/rann/ra_model_impl.hpp b/src/mlpack/methods/rann/ra_model_impl.hpp index 30452ef7ef..8184c355bd 100644 --- a/src/mlpack/methods/rann/ra_model_impl.hpp +++ b/src/mlpack/methods/rann/ra_model_impl.hpp @@ -26,7 +26,8 @@ RAModel::RAModel(const TreeTypes treeType, const bool randomBasis) : hilbertRTreeRA(NULL), rPlusTreeRA(NULL), rPlusPlusTreeRA(NULL), - ubTreeRA(NULL) + ubTreeRA(NULL), + octreeRA(NULL) { // Nothing to do. } @@ -34,24 +35,16 @@ RAModel::RAModel(const TreeTypes treeType, const bool randomBasis) : template RAModel::~RAModel() { - if (kdTreeRA) - delete kdTreeRA; - if (coverTreeRA) - delete coverTreeRA; - if (rTreeRA) - delete rTreeRA; - if (rStarTreeRA) - delete rStarTreeRA; - if (xTreeRA) - delete xTreeRA; - if (hilbertRTreeRA) - delete hilbertRTreeRA; - if (rPlusTreeRA) - delete rPlusTreeRA; - if (rPlusPlusTreeRA) - delete rPlusPlusTreeRA; - if (ubTreeRA) - delete ubTreeRA; + delete kdTreeRA; + delete coverTreeRA; + delete rTreeRA; + delete rStarTreeRA; + delete xTreeRA; + delete hilbertRTreeRA; + delete rPlusTreeRA; + delete rPlusPlusTreeRA; + delete ubTreeRA; + delete octreeRA; } template @@ -66,24 +59,16 @@ void RAModel::Serialize(Archive& ar, // This should never happen, but just in case, be clean with memory. if (Archive::is_loading::value) { - if (kdTreeRA) - delete kdTreeRA; - if (coverTreeRA) - delete coverTreeRA; - if (rTreeRA) - delete rTreeRA; - if (rStarTreeRA) - delete rStarTreeRA; - if (xTreeRA) - delete xTreeRA; - if (hilbertRTreeRA) - delete hilbertRTreeRA; - if (rPlusTreeRA) - delete rPlusTreeRA; - if (rPlusPlusTreeRA) - delete rPlusPlusTreeRA; - if (ubTreeRA) - delete ubTreeRA; + delete kdTreeRA; + delete coverTreeRA; + delete rTreeRA; + delete rStarTreeRA; + delete xTreeRA; + delete hilbertRTreeRA; + delete rPlusTreeRA; + delete rPlusPlusTreeRA; + delete ubTreeRA; + delete octreeRA; // Set all the pointers to NULL. kdTreeRA = NULL; @@ -127,6 +112,9 @@ void RAModel::Serialize(Archive& ar, case UB_TREE: ar & data::CreateNVP(ubTreeRA, "ra_model"); break; + case OCTREE: + ar & data::CreateNVP(octreeRA, "ra_model"); + break; } } @@ -151,6 +139,8 @@ const arma::mat& RAModel::Dataset() const return rPlusPlusTreeRA->ReferenceSet(); else if (ubTreeRA) return ubTreeRA->ReferenceSet(); + else if (octreeRA) + return octreeRA->ReferenceSet(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -177,6 +167,8 @@ bool RAModel::Naive() const return rPlusPlusTreeRA->Naive(); else if (ubTreeRA) return ubTreeRA->Naive(); + else if (octreeRA) + return octreeRA->Naive(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -203,6 +195,8 @@ bool& RAModel::Naive() return rPlusPlusTreeRA->Naive(); else if (ubTreeRA) return ubTreeRA->Naive(); + else if (octreeRA) + return octreeRA->Naive(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -229,6 +223,8 @@ bool RAModel::SingleMode() const return rPlusPlusTreeRA->SingleMode(); else if (ubTreeRA) return ubTreeRA->SingleMode(); + else if (octreeRA) + return octreeRA->SingleMode(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -255,6 +251,8 @@ bool& RAModel::SingleMode() return rPlusPlusTreeRA->SingleMode(); else if (ubTreeRA) return ubTreeRA->SingleMode(); + else if (octreeRA) + return octreeRA->SingleMode(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -281,6 +279,8 @@ double RAModel::Tau() const return rPlusPlusTreeRA->Tau(); else if (ubTreeRA) return ubTreeRA->Tau(); + else if (octreeRA) + return octreeRA->Tau(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -307,6 +307,8 @@ double& RAModel::Tau() return rPlusPlusTreeRA->Tau(); else if (ubTreeRA) return ubTreeRA->Tau(); + else if (octreeRA) + return octreeRA->Tau(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -333,6 +335,8 @@ double RAModel::Alpha() const return rPlusPlusTreeRA->Alpha(); else if (ubTreeRA) return ubTreeRA->Alpha(); + else if (octreeRA) + return octreeRA->Alpha(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -359,6 +363,8 @@ double& RAModel::Alpha() return rPlusPlusTreeRA->Alpha(); else if (ubTreeRA) return ubTreeRA->Alpha(); + else if (octreeRA) + return octreeRA->Alpha(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -385,6 +391,8 @@ bool RAModel::SampleAtLeaves() const return rPlusPlusTreeRA->SampleAtLeaves(); else if (ubTreeRA) return ubTreeRA->SampleAtLeaves(); + else if (octreeRA) + return octreeRA->SampleAtLeaves(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -411,6 +419,8 @@ bool& RAModel::SampleAtLeaves() return rPlusPlusTreeRA->SampleAtLeaves(); else if (ubTreeRA) return ubTreeRA->SampleAtLeaves(); + else if (octreeRA) + return octreeRA->SampleAtLeaves(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -437,6 +447,8 @@ bool RAModel::FirstLeafExact() const return rPlusPlusTreeRA->FirstLeafExact(); else if (ubTreeRA) return ubTreeRA->FirstLeafExact(); + else if (octreeRA) + return octreeRA->FirstLeafExact(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -463,6 +475,8 @@ bool& RAModel::FirstLeafExact() return rPlusPlusTreeRA->FirstLeafExact(); else if (ubTreeRA) return ubTreeRA->FirstLeafExact(); + else if (octreeRA) + return octreeRA->FirstLeafExact(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -489,6 +503,8 @@ size_t RAModel::SingleSampleLimit() const return rPlusPlusTreeRA->SingleSampleLimit(); else if (ubTreeRA) return ubTreeRA->SingleSampleLimit(); + else if (octreeRA) + return octreeRA->SingleSampleLimit(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -515,6 +531,8 @@ size_t& RAModel::SingleSampleLimit() return rPlusPlusTreeRA->SingleSampleLimit(); else if (ubTreeRA) return ubTreeRA->SingleSampleLimit(); + else if (octreeRA) + return octreeRA->SingleSampleLimit(); throw std::runtime_error("no rank-approximate nearest neighbor search model " "initialized"); @@ -570,24 +588,16 @@ void RAModel::BuildModel(arma::mat&& referenceSet, } // Clean memory, if necessary. - if (kdTreeRA) - delete kdTreeRA; - if (coverTreeRA) - delete coverTreeRA; - if (rTreeRA) - delete rTreeRA; - if (rStarTreeRA) - delete rStarTreeRA; - if (xTreeRA) - delete xTreeRA; - if (hilbertRTreeRA) - delete hilbertRTreeRA; - if (rPlusTreeRA) - delete rPlusTreeRA; - if (rPlusPlusTreeRA) - delete rPlusPlusTreeRA; - if (ubTreeRA) - delete ubTreeRA; + delete kdTreeRA; + delete coverTreeRA; + delete rTreeRA; + delete rStarTreeRA; + delete xTreeRA; + delete hilbertRTreeRA; + delete rPlusTreeRA; + delete rPlusPlusTreeRA; + delete ubTreeRA; + delete octreeRA; if (randomBasis) referenceSet = q * referenceSet; @@ -652,6 +662,26 @@ void RAModel::BuildModel(arma::mat&& referenceSet, ubTreeRA = new RAType(std::move(referenceSet), naive, singleMode); break; + case OCTREE: + // Build tree, if necessary. + if (naive) + { + octreeRA = new RAType(std::move(referenceSet), naive, + singleMode); + } + else + { + std::vector oldFromNewReferences; + typename RAType::Tree* octree = + new typename RAType::Tree(std::move(referenceSet), + oldFromNewReferences, leafSize); + octreeRA = new RAType(octree, singleMode); + + // Give the model ownership of the tree. + octreeRA->treeOwner = true; + octreeRA->oldFromNewReferences = oldFromNewReferences; + } + break; } if (!naive) @@ -745,6 +775,37 @@ void RAModel::Search(arma::mat&& querySet, // No mapping necessary. ubTreeRA->Search(querySet, k, neighbors, distances); break; + case OCTREE: + if (!octreeRA->Naive() && !octreeRA->SingleMode()) + { + // Build a second tree and search. + Timer::Start("tree_building"); + Log::Info << "Building query tree..." << std::endl; + std::vector oldFromNewQueries; + typename RAType::Tree queryTree(std::move(querySet), + oldFromNewQueries, leafSize); + Log::Info << "Tree built." << std::endl; + Timer::Stop("tree_building"); + + arma::Mat neighborsOut; + arma::mat distancesOut; + octreeRA->Search(&queryTree, k, neighborsOut, distancesOut); + + // Unmap the query points. + distances.set_size(distancesOut.n_rows, distancesOut.n_cols); + neighbors.set_size(neighborsOut.n_rows, neighborsOut.n_cols); + for (size_t i = 0; i < neighborsOut.n_cols; ++i) + { + neighbors.col(oldFromNewQueries[i]) = neighborsOut.col(i); + distances.col(oldFromNewQueries[i]) = distancesOut.col(i); + } + } + else + { + // Search without building a second tree. + octreeRA->Search(querySet, k, neighbors, distances); + } + break; } } @@ -791,6 +852,9 @@ void RAModel::Search(const size_t k, case UB_TREE: ubTreeRA->Search(k, neighbors, distances); break; + case OCTREE: + octreeRA->Search(k, neighbors, distances); + break; } } @@ -817,6 +881,8 @@ std::string RAModel::TreeName() const return "R++ tree"; case UB_TREE: return "UB tree"; + case OCTREE: + return "octree"; default: return "unknown tree"; } diff --git a/src/mlpack/tests/CMakeLists.txt b/src/mlpack/tests/CMakeLists.txt index 9ad40927c8..099aa2989e 100644 --- a/src/mlpack/tests/CMakeLists.txt +++ b/src/mlpack/tests/CMakeLists.txt @@ -57,6 +57,7 @@ add_executable(mlpack_test network_util_test.cpp nmf_test.cpp nystroem_method_test.cpp + octree_test.cpp pca_test.cpp perceptron_test.cpp quic_svd_test.cpp diff --git a/src/mlpack/tests/knn_test.cpp b/src/mlpack/tests/knn_test.cpp index 8045b71681..6f956dfe42 100644 --- a/src/mlpack/tests/knn_test.cpp +++ b/src/mlpack/tests/knn_test.cpp @@ -1066,7 +1066,7 @@ BOOST_AUTO_TEST_CASE(KNNModelTest) arma::mat referenceData = arma::randu(10, 200); // Build all the possible models. - KNNModel models[26]; + KNNModel models[28]; models[0] = KNNModel(KNNModel::TreeTypes::KD_TREE, true); models[1] = KNNModel(KNNModel::TreeTypes::KD_TREE, false); models[2] = KNNModel(KNNModel::TreeTypes::COVER_TREE, true); @@ -1093,6 +1093,8 @@ BOOST_AUTO_TEST_CASE(KNNModelTest) models[23] = KNNModel(KNNModel::TreeTypes::MAX_RP_TREE, false); models[24] = KNNModel(KNNModel::TreeTypes::UB_TREE, true); models[25] = KNNModel(KNNModel::TreeTypes::UB_TREE, false); + models[26] = KNNModel(KNNModel::TreeTypes::OCTREE, true); + models[27] = KNNModel(KNNModel::TreeTypes::OCTREE, false); for (size_t j = 0; j < 2; ++j) { @@ -1102,7 +1104,7 @@ BOOST_AUTO_TEST_CASE(KNNModelTest) arma::mat baselineDistances; knn.Search(queryData, 3, baselineNeighbors, baselineDistances); - for (size_t i = 0; i < 26; ++i) + for (size_t i = 0; i < 28; ++i) { // We only have std::move() constructors so make a copy of our data. arma::mat referenceCopy(referenceData); @@ -1147,7 +1149,7 @@ BOOST_AUTO_TEST_CASE(KNNModelMonochromaticTest) arma::mat referenceData = arma::randu(10, 200); // Build all the possible models. - KNNModel models[26]; + KNNModel models[28]; models[0] = KNNModel(KNNModel::TreeTypes::KD_TREE, true); models[1] = KNNModel(KNNModel::TreeTypes::KD_TREE, false); models[2] = KNNModel(KNNModel::TreeTypes::COVER_TREE, true); @@ -1174,6 +1176,8 @@ BOOST_AUTO_TEST_CASE(KNNModelMonochromaticTest) models[23] = KNNModel(KNNModel::TreeTypes::MAX_RP_TREE, false); models[24] = KNNModel(KNNModel::TreeTypes::UB_TREE, true); models[25] = KNNModel(KNNModel::TreeTypes::UB_TREE, false); + models[26] = KNNModel(KNNModel::TreeTypes::OCTREE, true); + models[27] = KNNModel(KNNModel::TreeTypes::OCTREE, false); for (size_t j = 0; j < 2; ++j) { @@ -1183,7 +1187,7 @@ BOOST_AUTO_TEST_CASE(KNNModelMonochromaticTest) arma::mat baselineDistances; knn.Search(3, baselineNeighbors, baselineDistances); - for (size_t i = 0; i < 26; ++i) + for (size_t i = 0; i < 28; ++i) { // We only have a std::move() constructor... so copy the data. arma::mat referenceCopy(referenceData); diff --git a/src/mlpack/tests/krann_search_test.cpp b/src/mlpack/tests/krann_search_test.cpp index 3d3f9184a1..c0bcadb45c 100644 --- a/src/mlpack/tests/krann_search_test.cpp +++ b/src/mlpack/tests/krann_search_test.cpp @@ -625,7 +625,7 @@ BOOST_AUTO_TEST_CASE(RAModelTest) data::Load("rann_test_q_3_100.csv", queryData, true); // Build all the possible models. - KNNModel models[18]; + KNNModel models[20]; models[0] = KNNModel(KNNModel::TreeTypes::KD_TREE, false); models[1] = KNNModel(KNNModel::TreeTypes::KD_TREE, true); models[2] = KNNModel(KNNModel::TreeTypes::COVER_TREE, false); @@ -644,13 +644,15 @@ BOOST_AUTO_TEST_CASE(RAModelTest) models[15] = KNNModel(KNNModel::TreeTypes::R_PLUS_PLUS_TREE, true); models[16] = KNNModel(KNNModel::TreeTypes::UB_TREE, false); models[17] = KNNModel(KNNModel::TreeTypes::UB_TREE, true); + models[18] = KNNModel(KNNModel::TreeTypes::OCTREE, false); + models[19] = KNNModel(KNNModel::TreeTypes::OCTREE, true); arma::Mat qrRanks; data::Load("rann_test_qr_ranks.csv", qrRanks, true, false); // No transpose. for (size_t j = 0; j < 3; ++j) { - for (size_t i = 0; i < 18; ++i) + for (size_t i = 0; i < 20; ++i) { // We only have std::move() constructors so make a copy of our data. arma::mat referenceCopy(referenceData); diff --git a/src/mlpack/tests/octree_test.cpp b/src/mlpack/tests/octree_test.cpp new file mode 100644 index 0000000000..f23701602c --- /dev/null +++ b/src/mlpack/tests/octree_test.cpp @@ -0,0 +1,342 @@ +/** + * @file octree_test.cpp + * @author Ryan Curtin + * + * Test various properties of the Octree. + */ +#include +#include + +#include +#include "test_tools.hpp" +#include "serialization.hpp" + +using namespace mlpack; +using namespace mlpack::math; +using namespace mlpack::tree; +using namespace mlpack::metric; +using namespace mlpack::bound; + +BOOST_AUTO_TEST_SUITE(OctreeTest); + +/** + * Build a quad-tree (2-d octree) on 4 points, and guarantee four points are + * created. + */ +BOOST_AUTO_TEST_CASE(SimpleQuadtreeTest) +{ + // Four corners of the unit square. + arma::mat dataset("0 0 1 1; 0 1 0 1"); + + Octree<> t(dataset, 1); + + BOOST_REQUIRE_EQUAL(t.NumChildren(), 4); + BOOST_REQUIRE_EQUAL(t.Dataset().n_cols, 4); + BOOST_REQUIRE_EQUAL(t.Dataset().n_rows, 2); + BOOST_REQUIRE_EQUAL(t.NumDescendants(), 4); + BOOST_REQUIRE_EQUAL(t.NumPoints(), 0); + for (size_t i = 0; i < 4; ++i) + { + BOOST_REQUIRE_EQUAL(t.Child(i).NumDescendants(), 1); + BOOST_REQUIRE_EQUAL(t.Child(i).NumPoints(), 1); + } +} + +/** + * Build an octree on 3 points and make sure that only three children are + * created. + */ +BOOST_AUTO_TEST_CASE(OctreeMissingChildTest) +{ + // Only three corners of the unit square. + arma::mat dataset("0 0 1; 0 1 1"); + + Octree<> t(dataset, 1); + + BOOST_REQUIRE_EQUAL(t.NumChildren(), 3); + BOOST_REQUIRE_EQUAL(t.Dataset().n_cols, 3); + BOOST_REQUIRE_EQUAL(t.Dataset().n_rows, 2); + BOOST_REQUIRE_EQUAL(t.NumDescendants(), 3); + BOOST_REQUIRE_EQUAL(t.NumPoints(), 0); + for (size_t i = 0; i < 3; ++i) + { + BOOST_REQUIRE_EQUAL(t.Child(i).NumDescendants(), 1); + BOOST_REQUIRE_EQUAL(t.Child(i).NumPoints(), 1); + } +} + +/** + * Ensure that building an empty octree does not fail. + */ +BOOST_AUTO_TEST_CASE(EmptyOctreeTest) +{ + arma::mat dataset; + Octree<> t(dataset); + + BOOST_REQUIRE_EQUAL(t.NumChildren(), 0); + BOOST_REQUIRE_EQUAL(t.Dataset().n_cols, 0); + BOOST_REQUIRE_EQUAL(t.Dataset().n_rows, 0); + BOOST_REQUIRE_EQUAL(t.NumDescendants(), 0); + BOOST_REQUIRE_EQUAL(t.NumPoints(), 0); +} + +/** + * Ensure that maxLeafSize is respected. + */ +BOOST_AUTO_TEST_CASE(MaxLeafSizeTest) +{ + arma::mat dataset(5, 15, arma::fill::randu); + Octree<> t1(dataset, 20); + Octree<> t2(std::move(dataset), 20); + + BOOST_REQUIRE_EQUAL(t1.NumChildren(), 0); + BOOST_REQUIRE_EQUAL(t1.NumDescendants(), 15); + BOOST_REQUIRE_EQUAL(t1.NumPoints(), 15); + + BOOST_REQUIRE_EQUAL(t2.NumChildren(), 0); + BOOST_REQUIRE_EQUAL(t2.NumDescendants(), 15); + BOOST_REQUIRE_EQUAL(t2.NumPoints(), 15); +} + +/** + * Check that the mappings given are correct. + */ +BOOST_AUTO_TEST_CASE(MappingsTest) +{ + // Test with both constructors. + arma::mat dataset(3, 5, arma::fill::randu); + arma::mat datacopy(dataset); + std::vector oldFromNewCopy, oldFromNewMove; + + Octree<> t1(dataset, oldFromNewCopy, 1); + Octree<> t2(std::move(dataset), oldFromNewMove, 1); + + for (size_t i = 0; i < oldFromNewCopy.size(); ++i) + { + BOOST_REQUIRE_SMALL(arma::norm(datacopy.col(oldFromNewCopy[i]) - + t1.Dataset().col(i)), 1e-3); + BOOST_REQUIRE_SMALL(arma::norm(datacopy.col(oldFromNewMove[i]) - + t2.Dataset().col(i)), 1e-3); + } +} + +/** + * Check that the reverse mappings are correct too. + */ +BOOST_AUTO_TEST_CASE(ReverseMappingsTest) +{ + // Test with both constructors. + arma::mat dataset(3, 300, arma::fill::randu); + arma::mat datacopy(dataset); + std::vector oldFromNewCopy, oldFromNewMove, newFromOldCopy, + newFromOldMove; + + Octree<> t1(dataset, oldFromNewCopy, newFromOldCopy); + Octree<> t2(std::move(dataset), oldFromNewMove, newFromOldMove); + + for (size_t i = 0; i < oldFromNewCopy.size(); ++i) + { + BOOST_REQUIRE_SMALL(arma::norm(datacopy.col(oldFromNewCopy[i]) - + t1.Dataset().col(i)), 1e-3); + BOOST_REQUIRE_SMALL(arma::norm(datacopy.col(oldFromNewMove[i]) - + t2.Dataset().col(i)), 1e-3); + + BOOST_REQUIRE_EQUAL(newFromOldCopy[oldFromNewCopy[i]], i); + BOOST_REQUIRE_EQUAL(newFromOldMove[oldFromNewMove[i]], i); + } +} + +/** + * Make sure no children at the same level are overlapping. + */ +template +void CheckOverlap(TreeType& node) +{ + // Check each combination of children. + for (size_t i = 0; i < node.NumChildren(); ++i) + for (size_t j = i + 1; j < node.NumChildren(); ++j) + BOOST_REQUIRE_EQUAL(node.Child(i).Bound().Overlap(node.Child(j).Bound()), + 0.0); // We need exact equality here. + + for (size_t i = 0; i < node.NumChildren(); ++i) + CheckOverlap(node.Child(i)); +} + +BOOST_AUTO_TEST_CASE(OverlapTest) +{ + // Test with both constructors. + arma::mat dataset(3, 300, arma::fill::randu); + + Octree<> t1(dataset); + Octree<> t2(std::move(dataset)); + + CheckOverlap(t1); + CheckOverlap(t2); +} + +/** + * Make sure no points are further than the furthest point distance, and that no + * descendants are further than the furthest descendant distance. + */ +template +void CheckFurthestDistances(TreeType& node) +{ + arma::vec center; + node.Center(center); + + // Compare points held in the node. + for (size_t i = 0; i < node.NumPoints(); ++i) + { + // Handle floating-point inaccuracies. + BOOST_REQUIRE_LE(metric::EuclideanDistance::Evaluate(node.Dataset().col(node.Point(i)), + center), node.FurthestPointDistance() * (1 + 1e-5)); + } + + // Compare descendants held in the node. + for (size_t i = 0; i < node.NumDescendants(); ++i) + { + // Handle floating-point inaccuracies. + BOOST_REQUIRE_LE(metric::EuclideanDistance::Evaluate(node.Dataset().col(node.Descendant(i)), + center), node.FurthestDescendantDistance() * (1 + 1e-5)); + } + + for (size_t i = 0; i < node.NumChildren(); ++i) + CheckFurthestDistances(node.Child(i)); +} + +BOOST_AUTO_TEST_CASE(FurthestDistanceTest) +{ + // Test with both constructors. + arma::mat dataset(3, 500, arma::fill::randu); + + Octree<> t1(dataset); + Octree<> t2(std::move(dataset)); + + CheckFurthestDistances(t1); + CheckFurthestDistances(t2); +} + +/** + * The maximum number of children a node can have is limited by the + * dimensionality. So we test to make sure there are no cases where we have too + * many children. + */ +template +void CheckNumChildren(TreeType& node) +{ + BOOST_REQUIRE_LE(node.NumChildren(), std::pow(2, node.Dataset().n_rows)); + for (size_t i = 0; i < node.NumChildren(); ++i) + CheckNumChildren(node.Child(i)); +} + +BOOST_AUTO_TEST_CASE(MaxNumChildrenTest) +{ + for (size_t d = 1; d < 10; ++d) + { + arma::mat dataset(d, 1000 * d, arma::fill::randu); + Octree<> t(std::move(dataset)); + + CheckNumChildren(t); + } +} + +/** + * Test the copy constructor. + */ +template +void CheckSameNode(TreeType& node1, TreeType& node2) +{ + BOOST_REQUIRE_EQUAL(node1.NumChildren(), node2.NumChildren()); + BOOST_REQUIRE_NE(&node1.Dataset(), &node2.Dataset()); + + // Make sure the children actually got copied. + for (size_t i = 0; i < node1.NumChildren(); ++i) + BOOST_REQUIRE_NE(&node1.Child(i), &node2.Child(i)); + + // Check that all the points are the same. + BOOST_REQUIRE_EQUAL(node1.NumPoints(), node2.NumPoints()); + BOOST_REQUIRE_EQUAL(node1.NumDescendants(), node2.NumDescendants()); + for (size_t i = 0; i < node1.NumPoints(); ++i) + BOOST_REQUIRE_EQUAL(node1.Point(i), node2.Point(i)); + for (size_t i = 0; i < node1.NumDescendants(); ++i) + BOOST_REQUIRE_EQUAL(node1.Descendant(i), node2.Descendant(i)); + + // Check that the bound is the same. + BOOST_REQUIRE_EQUAL(node1.Bound().Dim(), node2.Bound().Dim()); + for (size_t d = 0; d < node1.Bound().Dim(); ++d) + { + BOOST_REQUIRE_CLOSE(node1.Bound()[d].Lo(), node2.Bound()[d].Lo(), 1e-5); + BOOST_REQUIRE_CLOSE(node1.Bound()[d].Hi(), node2.Bound()[d].Hi(), 1e-5); + } + + // Check that the furthest point and descendant distance are the same. + BOOST_REQUIRE_CLOSE(node1.FurthestPointDistance(), + node2.FurthestPointDistance(), 1e-5); + BOOST_REQUIRE_CLOSE(node1.FurthestDescendantDistance(), + node2.FurthestDescendantDistance(), 1e-5); +} + +BOOST_AUTO_TEST_CASE(CopyConstructorTest) +{ + // Use a small random dataset. + arma::mat dataset(3, 100, arma::fill::randu); + + Octree<> t(dataset); + Octree<> t2(t); + + CheckSameNode(t, t2); +} + +/** + * Test the move constructor. + */ +BOOST_AUTO_TEST_CASE(MoveConstructorTest) +{ + // Use a small random dataset. + arma::mat dataset(3, 100, arma::fill::randu); + + Octree<> t(std::move(dataset)); + Octree<> tcopy(t); + + // Move the tree. + Octree<> t2(std::move(t)); + + // Make sure the original tree has no data. + BOOST_REQUIRE_EQUAL(t.Dataset().n_rows, 0); + BOOST_REQUIRE_EQUAL(t.Dataset().n_cols, 0); + BOOST_REQUIRE_EQUAL(t.NumChildren(), 0); + BOOST_REQUIRE_EQUAL(t.NumPoints(), 0); + BOOST_REQUIRE_EQUAL(t.NumDescendants(), 0); + BOOST_REQUIRE_SMALL(t.FurthestPointDistance(), 1e-5); + BOOST_REQUIRE_SMALL(t.FurthestDescendantDistance(), 1e-5); + BOOST_REQUIRE_EQUAL(t.Bound().Dim(), 0); + + // Check that the new tree is the same as our copy. + CheckSameNode(tcopy, t2); +} + +/** + * Test serialization. + */ +BOOST_AUTO_TEST_CASE(SerializationTest) +{ + // Use a small random dataset. + arma::mat dataset(3, 500, arma::fill::randu); + Octree<> t(std::move(dataset)); + + Octree<>* xmlTree; + Octree<>* binaryTree; + Octree<>* textTree; + + SerializePointerObjectAll(&t, xmlTree, binaryTree, textTree); + + CheckSameNode(t, *xmlTree); + CheckSameNode(t, *binaryTree); + CheckSameNode(t, *textTree); + + delete xmlTree; + delete binaryTree; + delete textTree; +} + +BOOST_AUTO_TEST_SUITE_END(); diff --git a/src/mlpack/tests/range_search_test.cpp b/src/mlpack/tests/range_search_test.cpp index 0968b07594..361969157d 100644 --- a/src/mlpack/tests/range_search_test.cpp +++ b/src/mlpack/tests/range_search_test.cpp @@ -1249,7 +1249,7 @@ BOOST_AUTO_TEST_CASE(RSModelTest) arma::mat referenceData = arma::randu(10, 200); // Build all the possible models. - RSModel models[26]; + RSModel models[28]; models[0] = RSModel(RSModel::TreeTypes::KD_TREE, true); models[1] = RSModel(RSModel::TreeTypes::KD_TREE, false); models[2] = RSModel(RSModel::TreeTypes::COVER_TREE, true); @@ -1276,6 +1276,8 @@ BOOST_AUTO_TEST_CASE(RSModelTest) models[23] = RSModel(RSModel::TreeTypes::MAX_RP_TREE, false); models[24] = RSModel(RSModel::TreeTypes::UB_TREE, true); models[25] = RSModel(RSModel::TreeTypes::UB_TREE, false); + models[26] = RSModel(RSModel::TreeTypes::OCTREE, true); + models[27] = RSModel(RSModel::TreeTypes::OCTREE, false); for (size_t j = 0; j < 2; ++j) { @@ -1289,7 +1291,7 @@ BOOST_AUTO_TEST_CASE(RSModelTest) vector>> baselineSorted; SortResults(baselineNeighbors, baselineDistances, baselineSorted); - for (size_t i = 0; i < 26; ++i) + for (size_t i = 0; i < 28; ++i) { // We only have std::move() constructors, so make a copy of our data. arma::mat referenceCopy(referenceData); @@ -1333,7 +1335,7 @@ BOOST_AUTO_TEST_CASE(RSModelMonochromaticTest) arma::mat referenceData = arma::randu(10, 200); // Build all the possible models. - RSModel models[26]; + RSModel models[28]; models[0] = RSModel(RSModel::TreeTypes::KD_TREE, true); models[1] = RSModel(RSModel::TreeTypes::KD_TREE, false); models[2] = RSModel(RSModel::TreeTypes::COVER_TREE, true); @@ -1360,6 +1362,8 @@ BOOST_AUTO_TEST_CASE(RSModelMonochromaticTest) models[23] = RSModel(RSModel::TreeTypes::MAX_RP_TREE, false); models[24] = RSModel(RSModel::TreeTypes::MAX_RP_TREE, true); models[25] = RSModel(RSModel::TreeTypes::MAX_RP_TREE, false); + models[26] = RSModel(RSModel::TreeTypes::OCTREE, true); + models[27] = RSModel(RSModel::TreeTypes::OCTREE, false); for (size_t j = 0; j < 2; ++j) { @@ -1372,7 +1376,7 @@ BOOST_AUTO_TEST_CASE(RSModelMonochromaticTest) vector>> baselineSorted; SortResults(baselineNeighbors, baselineDistances, baselineSorted); - for (size_t i = 0; i < 26; ++i) + for (size_t i = 0; i < 28; ++i) { // We only have std::move() cosntructors, so make a copy of our data. arma::mat referenceCopy(referenceData);