Merge branch 'master' into lsh-multiprobe

Merges e6bc4b4

src/mlpack/core/tree/CMakeLists.txt

@@ -43,6 +43,7 @@ set(SOURCES
   rectangle_tree/dual_tree_traverser_impl.hpp
   rectangle_tree/r_tree_split.hpp
   rectangle_tree/r_tree_split_impl.hpp
+  rectangle_tree/no_auxiliary_information.hpp
   rectangle_tree/r_tree_descent_heuristic.hpp
   rectangle_tree/r_tree_descent_heuristic_impl.hpp
   rectangle_tree/r_star_tree_descent_heuristic.hpp
@@ -51,6 +52,15 @@ set(SOURCES
   rectangle_tree/r_star_tree_split_impl.hpp
   rectangle_tree/x_tree_split.hpp
   rectangle_tree/x_tree_split_impl.hpp
+  rectangle_tree/x_tree_auxiliary_information.hpp
+  rectangle_tree/hilbert_r_tree_descent_heuristic.hpp
+  rectangle_tree/hilbert_r_tree_descent_heuristic_impl.hpp
+  rectangle_tree/hilbert_r_tree_split.hpp
+  rectangle_tree/hilbert_r_tree_split_impl.hpp
+  rectangle_tree/hilbert_r_tree_auxiliary_information.hpp
+  rectangle_tree/hilbert_r_tree_auxiliary_information_impl.hpp
+  rectangle_tree/discrete_hilbert_value.hpp
+  rectangle_tree/discrete_hilbert_value_impl.hpp
   statistic.hpp
   traversal_info.hpp
   tree_traits.hpp

src/mlpack/core/tree/rectangle_tree.hpp

@@ -19,10 +19,16 @@
 #include "rectangle_tree/dual_tree_traverser_impl.hpp"
 #include "rectangle_tree/r_tree_split.hpp"
 #include "rectangle_tree/r_star_tree_split.hpp"
+#include "rectangle_tree/no_auxiliary_information.hpp"
 #include "rectangle_tree/r_tree_descent_heuristic.hpp"
 #include "rectangle_tree/r_star_tree_descent_heuristic.hpp"
 #include "rectangle_tree/traits.hpp"
 #include "rectangle_tree/x_tree_split.hpp"
+#include "rectangle_tree/x_tree_auxiliary_information.hpp"
+#include "rectangle_tree/hilbert_r_tree_descent_heuristic.hpp"
+#include "rectangle_tree/hilbert_r_tree_split.hpp"
+#include "rectangle_tree/hilbert_r_tree_auxiliary_information.hpp"
+#include "rectangle_tree/discrete_hilbert_value.hpp"
 #include "rectangle_tree/typedef.hpp"
 
 #endif

src/mlpack/core/tree/rectangle_tree/discrete_hilbert_value.hpp

@@ -0,0 +1,263 @@
+/**
+ * @file discrete_hilbert_value.hpp
+ * @author Mikhail Lozhnikov
+ *
+ * Defintion of the DiscreteHilbertValue class, a class that calculates
+ * the ordering of points using the Hilbert curve.
+ */
+#ifndef MLPACK_CORE_TREE_RECTANGLE_TREE_DISCRETE_HILBERT_VALUE_HPP
+#define MLPACK_CORE_TREE_RECTANGLE_TREE_DISCRETE_HILBERT_VALUE_HPP
+
+#include <mlpack/core.hpp>
+
+namespace mlpack {
+namespace tree /** Trees and tree-building procedures. */ {
+
+/**
+ * The DiscreteHilbertValue class stores Hilbert values for all of the points in
+ * a RectangleTree node, and calculates Hilbert values for new points.  This
+ * implementation calculates the full discrete Hilbert value; for a
+ * d-dimensional vector filled with elements of size E, each Hilbert value will
+ * take dE space.
+ */
+template<typename TreeElemType>
+class DiscreteHilbertValue
+{
+ public:
+  //! Depending on the precision of the tree element type, we may need to use
+  //! uint32_t or uint64_t.
+  typedef typename std::conditional<sizeof(TreeElemType) * CHAR_BIT <= 32,
+                                    uint32_t,
+                                    uint64_t>::type HilbertElemType;
+
+  //! Default constructor.
+  DiscreteHilbertValue();
+
+  /**
+   * Construct this for the node tree. If the node is the root this method
+   * computes the Hilbert value for each point in the tree's dataset.
+   *
+   * @param node The node that stores this Hilbert value.
+   */
+  template<typename TreeType>
+  DiscreteHilbertValue(const TreeType* tree);
+
+  /**
+   * Create a Hilbert value object by copying from another one.
+   *
+   * @param other The Hilbert value object from which the value will be copied.
+   */
+  DiscreteHilbertValue(const DiscreteHilbertValue& other);
+
+  //! Free memory
+  ~DiscreteHilbertValue();
+
+  /**
+   * Compare two points. It returns 1 if the first point is greater than the
+   * second one, -1 if the first point is less than the second one and 0 if the
+   * Hilbert values of the points are equal. In order to do it this method
+   * computes the Hilbert values of the points.
+   *
+   * @param pt1 The first point.
+   * @param pt2 The second point.
+   */
+  template<typename VecType1, typename VecType2>
+  static int ComparePoints(const VecType1& pt1, const VecType2& pt2,
+                           typename boost::enable_if<IsVector<VecType1>>* = 0,
+                           typename boost::enable_if<IsVector<VecType2>>* = 0);
+
+  /**
+   * Compare two Hilbert values. It returns 1 if the first value is greater than
+   * the second one, -1 if the first value is less than the second one and 0 if
+   * the values are equal. This method does not compute the Hilbert values.
+   *
+   * @param val1 The first point.
+   * @param val2 The second point.
+   */
+  static int CompareValues(const DiscreteHilbertValue& val1,
+                           const DiscreteHilbertValue& val2);
+
+  /**
+   * Compare the largest Hilbert value of the node with the val value.  It
+   * returns 1 if the value of the node is greater than val, -1 if the value of
+   * the node is less than val and 0 if the values are equal. This method does
+   * not compute the Hilbert values.
+   *
+   * @param val The Hilbert value to compare with.
+   */
+  int CompareWith(const DiscreteHilbertValue& val) const;
+
+  /**
+   * Compare the largest Hilbert value of the node with the Hilbert value of the
+   * point. It returns 1 if the value of the node is greater than the value of
+   * the point, -1 if the value of the node is less than the value of the point
+   * and 0 if the values are equal.  This method computes the Hilbert value of
+   * the point.
+   *
+   * @param pt The point to compare with.
+   */
+  template<typename VecType>
+  int CompareWith(const VecType& pt,
+                  typename boost::enable_if<IsVector<VecType>>* = 0) const;
+
+  /**
+   * Compare the Hilbert value of the cached point with the Hilbert value of the
+   * given point. It returns 1 if the value of the node is greater than the
+   * value of the point, -1 if the value of the node is less than the value of
+   * the point and 0 if the values are equal.  This method computes the Hilbert
+   * value of the point.
+   *
+   * @param pt The point to compare with.
+   */
+
+  template<typename VecType>
+  int CompareWithCachedPoint(
+      const VecType& pt,
+      typename boost::enable_if<IsVector<VecType>>* = 0) const;
+
+  /**
+   * Update the largest Hilbert value of the node and insert the point in the
+   * local dataset if the node is a leaf.
+   *
+   * @param node The node in which the point is being inserted.
+   * @param point The number of the point being inserted.
+   */
+  template<typename TreeType, typename VecType>
+  size_t InsertPoint(TreeType *node,
+                     const VecType& pt,
+                     typename boost::enable_if<IsVector<VecType>>* = 0);
+
+  /**
+   * Update the largest Hilbert value of the node.
+   *
+   * @param node The node being inserted.
+   */
+  template<typename TreeType>
+  void InsertNode(TreeType* node);
+
+  /**
+   * Update the largest Hilbert value of the node and delete the point from the
+   * local dataset.
+   *
+   * @param node The node from which the point is being deleted.
+   * @param localIndex The index of the point in the local dataset.
+   */
+  template<typename TreeType>
+  void DeletePoint(TreeType* node, const size_t localIndex);
+
+  /**
+   * Update the largest Hilbert value of the node.
+   *
+   * @param node The node from which another node is being deleted.
+   * @param nodeIndex The index of the node being deleted.
+   */
+  template<typename TreeType>
+  void RemoveNode(TreeType* node, const size_t nodeIndex);
+
+  /**
+   * Copy the local Hilbert value's pointer.
+   *
+   * @param val The DiscreteHilbertValue object from which the dataset
+   *    will be copied.
+   */
+  DiscreteHilbertValue& operator=(const DiscreteHilbertValue& val);
+
+  /**
+   * Nullify the localHilbertValues pointer in order to prevent an invalid free.
+   */
+  void NullifyData();
+
+  /**
+   * Update the largest Hilbert value and the local Hilbert values of an
+   * intermediate node.  The children of the node (or the points that the node
+   * contains) should be arranged according to their Hilbert values.
+   *
+   * @param node The node in which the information should be updated.
+   */
+  template<typename TreeType>
+  void UpdateLargestValue(TreeType* node);
+
+  /**
+   * This method updates the largest Hilbert value of a leaf node and
+   * redistributes the Hilbert values of points according to their new position
+   * after the split algorithm.
+   *
+   * @param parent The parent of the node that was split.
+   * @param firstSibling The first cooperating sibling.
+   * @param lastSibling The last cooperating sibling.
+   */
+  template<typename TreeType>
+  void RedistributeHilbertValues(TreeType* parent,
+                                 const size_t firstSibling,
+                                 const size_t lastSibling);
+
+  /**
+   * Calculate the Hilbert value of the point pt.
+   *
+   * @param pt The point for which the Hilbert value should be calculated.
+   */
+  template<typename VecType>
+  static arma::Col<HilbertElemType> CalculateValue(
+      const VecType& pt,
+      typename boost::enable_if<IsVector<VecType>>* = 0);
+
+  /**
+   * Compare two Hilbert values. It returns 1 if the first value is greater than
+   * the second one, -1 if the first value is less than the second one and 0 if
+   * the values are equal. This method does not compute the Hilbert values.
+   *
+   * @param value1 The first value.
+   * @param value2 The second value.
+   */
+  static int CompareValues(const arma::Col<HilbertElemType>& value1,
+                           const arma::Col<HilbertElemType>& value2);
+
+  //! Return the number of values.
+  size_t NumValues() const { return numValues; }
+  //! Modify the number of values.
+  size_t& NumValues() { return numValues; }
+
+  //! Return the Hilbert values.
+  const arma::Mat<HilbertElemType>* LocalHilbertValues() const
+  { return localHilbertValues; }
+  //! Modify the Hilbert values.
+  arma::Mat<HilbertElemType>*& LocalHilbertValues()
+  { return localHilbertValues; }
+
+  //! Return the cached point (valueToInsert).
+  const arma::Col<HilbertElemType>* ValueToInsert() const
+  { return valueToInsert; }
+  //! Modify the cached point (valueToInsert).
+  arma::Col<HilbertElemType>* ValueToInsert() { return valueToInsert; }
+
+ private:
+  //! The number of bits that we can store.
+  static constexpr size_t order = sizeof(HilbertElemType) * CHAR_BIT;
+  //! The local Hilbert values.
+  arma::Mat<HilbertElemType>* localHilbertValues;
+  //! Indicates that the node owns the localHilbertValues variable.
+  bool ownsLocalHilbertValues;
+  //! The number of values in the localHilbertValues dataset.
+  size_t numValues;
+  /** The Hilbert value of the point that is being inserted.
+   * The pointer is the same in all nodes. The value is updated in InsertPoint()
+   * if it is invoked at the root level. This variable helps to avoid
+   * multiple computation of the Hilbert value of a point in the insertion
+   * process.
+   */
+  arma::Col<HilbertElemType>* valueToInsert;
+  //! Indicates that the node owns the valueToInsert.
+  bool ownsValueToInsert;
+
+ public:
+  template<typename Archive>
+  void Serialize(Archive& ar, const unsigned int /* version */);
+};
+
+} // namespace tree
+} // namespace mlpack
+
+// Include implementation.
+#include "discrete_hilbert_value_impl.hpp"
+
+#endif // MLPACK_CORE_TREE_RECTANGLE_TREE_DISCRETE_HILBERT_VALUE_HPP

src/mlpack/core/tree/rectangle_tree/discrete_hilbert_value_impl.hpp

@@ -0,0 +1,457 @@
+/**
+ * @file discrete_hilbert_value.hpp
+ * @author Mikhail Lozhnikov
+ *
+ * Defintion of the DiscreteHilbertValue class, a class that calculates
+ * the ordering of points using the Hilbert curve.
+ */
+#ifndef MLPACK_CORE_TREE_RECTANGLE_TREE_DISCRETE_HILBERT_VALUE_IMPL_HPP
+#define MLPACK_CORE_TREE_RECTANGLE_TREE_DISCRETE_HILBERT_VALUE_IMPL_HPP
+
+#include "discrete_hilbert_value.hpp"
+
+namespace mlpack {
+namespace tree /** Trees and tree-building procedures. */ {
+
+template<typename TreeElemType>
+DiscreteHilbertValue<TreeElemType>::DiscreteHilbertValue() :
+    localHilbertValues(NULL),
+    ownsLocalHilbertValues(false),
+    numValues(0),
+    valueToInsert(NULL),
+    ownsValueToInsert(false)
+{ }
+
+template<typename TreeElemType>
+DiscreteHilbertValue<TreeElemType>::~DiscreteHilbertValue()
+{
+  if (ownsLocalHilbertValues)
+    delete localHilbertValues;
+  if (ownsValueToInsert)
+    delete valueToInsert;
+}
+
+template<typename TreeElemType>
+template<typename TreeType>
+DiscreteHilbertValue<TreeElemType>::DiscreteHilbertValue(const TreeType* tree) :
+    localHilbertValues(NULL),
+    ownsLocalHilbertValues(false),
+    numValues(0),
+    valueToInsert(tree->Parent() ?
+        tree->Parent()->AuxiliaryInfo().HilbertValue().ValueToInsert() :
+        new arma::Col<HilbertElemType>(tree->Dataset().n_rows)),
+    ownsValueToInsert(tree->Parent() ? false : true)
+{
+  // Calculate the Hilbert value for all points.
+  if (!tree->Parent()) // This is the root node.
+    ownsLocalHilbertValues = true;
+  else if (tree->Parent()->Children()[0]->IsLeaf())
+  {
+    // This is a leaf node.
+    assert(tree->Parent()->NumChildren() > 0);
+    ownsLocalHilbertValues = true;
+  }
+
+  if (ownsLocalHilbertValues)
+  {
+    localHilbertValues = new arma::Mat<HilbertElemType>(tree->Dataset().n_rows,
+        tree->MaxLeafSize() + 1);
+  }
+}
+
+template<typename TreeElemType>
+DiscreteHilbertValue<TreeElemType>::
+DiscreteHilbertValue(const DiscreteHilbertValue& other) :
+    localHilbertValues(
+        const_cast<arma::Mat<HilbertElemType>*>(other.LocalHilbertValues())),
+    ownsLocalHilbertValues(other.ownsLocalHilbertValues),
+    numValues(other.NumValues()),
+    valueToInsert(
+        const_cast<arma::Col<HilbertElemType>*>(other.ValueToInsert())),
+    ownsValueToInsert(false)
+{ }
+
+template<typename TreeElemType>
+template<typename VecType>
+arma::Col<typename DiscreteHilbertValue<TreeElemType>::HilbertElemType>
+DiscreteHilbertValue<TreeElemType>::
+CalculateValue(const VecType& pt,typename boost::enable_if<IsVector<VecType>>*)
+{
+  typedef typename VecType::elem_type VecElemType;
+  arma::Col<HilbertElemType> res(pt.n_rows);
+  // Calculate the number of bits for the exponent.
+  const int numExpBits = std::ceil(std::log2(
+      std::numeric_limits<VecElemType>::max_exponent -
+      std::numeric_limits<VecElemType>::min_exponent + 1.0));
+
+  // Calculate the number of bits for the mantissa.
+  const int numMantBits = order - numExpBits - 1;
+
+  for (size_t i = 0; i < pt.n_rows; i++)
+  {
+    int e;
+    VecElemType normalizedVal = std::frexp(pt(i),&e);
+    bool sgn = std::signbit(normalizedVal);
+
+    if (pt(i) == 0)
+      e = std::numeric_limits<VecElemType>::min_exponent;
+
+    if (sgn)
+      normalizedVal = -normalizedVal;
+
+    if (e < std::numeric_limits<VecElemType>::min_exponent)
+    {
+      HilbertElemType tmp = (HilbertElemType) 1 <<
+          (std::numeric_limits<VecElemType>::min_exponent - e);
+
+      e = std::numeric_limits<VecElemType>::min_exponent;
+      normalizedVal /= tmp;
+    }
+
+    // Extract the mantissa.
+    HilbertElemType tmp = (HilbertElemType) 1 << numMantBits;
+    res(i) = std::floor(normalizedVal * tmp);
+
+    // Add the exponent.
+    assert(res(i) < ((HilbertElemType) 1 << numMantBits));
+    res(i) |= ((HilbertElemType)
+        (e - std::numeric_limits<VecElemType>::min_exponent)) << numMantBits;
+
+    assert(res(i) < ((HilbertElemType) 1 << (order - 1)) - 1);
+
+    // Negative values should be inverted.
+    if (sgn)
+    {
+      res(i) = ((HilbertElemType) 1 << (order - 1)) - 1 - res(i);
+      assert((res(i) >> (order - 1)) == 0);
+    }
+    else
+    {
+      res(i) |= (HilbertElemType) 1 << (order - 1);
+      assert((res(i) >> (order - 1)) == 1);
+    }
+  }
+
+  HilbertElemType M = (HilbertElemType) 1 << (order - 1);
+
+  // Since the Hilbert curve is continuous we should permutate and intend
+  // coordinate axes depending on the position of the point.
+  for (HilbertElemType Q = M; Q > 1; Q >>= 1)
+  {
+    HilbertElemType P = Q - 1;
+
+    for (size_t i = 0; i < pt.n_rows; i++)
+    {
+      if (res(i) & Q) // Invert.
+        res(0) ^= P;
+      else // Permutate.
+      {
+        HilbertElemType t = (res(0) ^ res(i)) & P;
+        res(0) ^= t;
+        res(i) ^= t;
+      }
+    }
+  }
+
+  // Gray encode.
+  for (size_t i = 1; i < pt.n_rows; i++)
+    res(i) ^= res(i - 1);
+
+  HilbertElemType t = 0;
+
+  // Some coordinate axes should be inverted.
+  for (HilbertElemType Q = M; Q > 1; Q >>= 1)
+    if (res(pt.n_rows - 1) & Q)
+      t ^= Q - 1;
+
+  for (size_t i = 0; i < pt.n_rows; i++)
+    res(i) ^= t;
+
+  // We should rearrange bits in order to compare two Hilbert values faster.
+  arma::Col<HilbertElemType> rearrangedResult(pt.n_rows, arma::fill::zeros);
+
+  for (size_t i = 0; i < order; i++)
+    for (size_t j = 0; j < pt.n_rows; j++)
+    {
+      size_t bit = (i * pt.n_rows + j) % order;
+      size_t row = (i * pt.n_rows + j) / order;
+
+      rearrangedResult(row) |= (((res(j) >> (order - 1 - i)) & 1) <<
+          (order - 1 - bit));
+    }
+
+  return rearrangedResult;
+}
+
+template<typename TreeElemType>
+int DiscreteHilbertValue<TreeElemType>::
+CompareValues(const arma::Col<HilbertElemType>& value1,
+              const arma::Col<HilbertElemType>& value2)
+{
+  for (size_t i = 0; i < value1.n_rows; i++)
+  {
+    if (value1(i) > value2(i))
+      return 1;
+    else if (value1(i) < value2(i))
+      return -1;
+  }
+
+  return 0;
+}
+
+
+
+template<typename TreeElemType>
+template<typename VecType1, typename VecType2>
+int DiscreteHilbertValue<TreeElemType>::
+ComparePoints(const VecType1& pt1,
+              const VecType2& pt2,
+              typename boost::enable_if<IsVector<VecType1>>*,
+              typename boost::enable_if<IsVector<VecType2>>*)
+{
+  arma::Col<HilbertElemType> val1 = CalculateValue(pt1);
+  arma::Col<HilbertElemType> val2 = CalculateValue(pt2);
+
+  return CompareValues(val1, val2);
+}
+
+template<typename TreeElemType>
+int DiscreteHilbertValue<TreeElemType>::
+CompareValues(const DiscreteHilbertValue& val1,
+              const DiscreteHilbertValue& val2)
+{
+  if (val1.NumValues() > 0 && val2.NumValues() == 0)
+    return 1;
+  else if (val1.NumValues() == 0 && val2.NumValues() > 0)
+    return -1;
+  else if (val1.NumValues() == 0 && val2.NumValues() == 0)
+    return 0;
+
+  return CompareValues(val1.LocalHilbertValues()->col(val1.NumValues() - 1),
+                       val2.LocalHilbertValues()->col(val2.NumValues() - 1));
+}
+
+template<typename TreeElemType>
+int DiscreteHilbertValue<TreeElemType>::
+CompareWith(const DiscreteHilbertValue& val) const
+{
+  return CompareValues(*this, val);
+}
+
+template<typename TreeElemType>
+template<typename VecType>
+int DiscreteHilbertValue<TreeElemType>::
+CompareWith(const VecType& pt,
+            typename boost::enable_if<IsVector<VecType>>*) const
+{
+  arma::Col<HilbertElemType> val = CalculateValue(pt);
+
+  if (numValues == 0)
+    return -1;
+
+  return CompareValues(localHilbertValues->col(numValues - 1),val);
+}
+
+template<typename TreeElemType>
+template<typename VecType>
+int DiscreteHilbertValue<TreeElemType>::
+CompareWithCachedPoint(const VecType& ,
+            typename boost::enable_if<IsVector<VecType>>*) const
+{
+  if (numValues == 0)
+    return -1;
+
+  return CompareValues(localHilbertValues->col(numValues - 1), *valueToInsert);
+}
+
+template<typename TreeElemType>
+template<typename TreeType, typename VecType>
+size_t DiscreteHilbertValue<TreeElemType>::
+InsertPoint(TreeType *node,
+            const VecType& pt,
+            typename boost::enable_if<IsVector<VecType>>*)
+{
+  size_t i = 0;
+
+  // All points are inserted to the root node.
+  if (!node->Parent())
+    *valueToInsert = CalculateValue(pt);
+  if (node->IsLeaf())
+  {
+    // Find an appropriate place.
+    for (i = 0; i < numValues; i++)
+      if (CompareValues(localHilbertValues->col(i), *valueToInsert) > 0)
+        break;
+
+    for (size_t j = numValues; j > i; j--)
+      localHilbertValues->col(j) = localHilbertValues->col(j-1);
+
+    localHilbertValues->col(i) = *valueToInsert;
+    numValues++;
+    // Propagate changes of the largest Hilbert value downward.
+    TreeType* root = node->Parent();
+
+    while (root != NULL)
+    {
+      root->AuxiliaryInfo().HilbertValue().UpdateLargestValue(root);
+
+      root = root->Parent();
+    }
+  }
+
+  return i;
+}
+
+template<typename TreeElemType>
+template<typename TreeType>
+void DiscreteHilbertValue<TreeElemType>::InsertNode(TreeType* node)
+{
+  DiscreteHilbertValue &val = node->AuxiliaryInfo().HilbertValue();
+
+  if (CompareWith(node,val) < 0)
+  {
+    localHilbertValues = val.LocalHilbertValues();
+    numValues = val.NumValues();
+  }
+}
+
+template<typename TreeElemType>
+template<typename TreeType>
+void DiscreteHilbertValue<TreeElemType>::
+DeletePoint(TreeType* node, const size_t localIndex)
+{
+
+  // Delete the Hilbert value from the local dataset
+  for (size_t i = numValues - 1; i > localIndex; i--)
+    localHilbertValues->col(i - 1) = localHilbertValues->col(i);
+
+  numValues--;
+}
+
+template<typename TreeElemType>
+template<typename TreeType>
+void DiscreteHilbertValue<TreeElemType>::
+RemoveNode(TreeType* node, const size_t nodeIndex)
+{
+  if (node->NumChildren() <= 1)
+  {
+    localHilbertValues = NULL;
+    numValues = 0;
+    return;
+  }
+  if (nodeIndex + 1 == node->NumChildren())
+  {
+    // Update the largest Hilbert value if the value exists
+    TreeType* child = node->Children()[nodeIndex - 1];
+    if (child->AuxiliaryInfo.HilbertValue().NumValues() != 0)
+    {
+      numValues = child->AuxiliaryInfo.HilbertValue().NumValues();
+      localHilbertValues =
+          child->AuxiliaryInfo.HilbertValue().LocalHilbertValues();
+    }
+    else
+    {
+      localHilbertValues = NULL;
+      numValues = 0;
+    }
+  }
+}
+
+template<typename TreeElemType>
+DiscreteHilbertValue<TreeElemType>& DiscreteHilbertValue<TreeElemType>::
+operator=(const DiscreteHilbertValue& val)
+{
+  localHilbertValues = const_cast<arma::Mat<HilbertElemType>* >
+      (val.LocalHilbertValues());
+  ownsLocalHilbertValues = false;
+  numValues = val.NumValues();
+
+  return *this;
+}
+
+template<typename TreeElemType>
+void DiscreteHilbertValue<TreeElemType>::NullifyData()
+{
+  ownsLocalHilbertValues = false;
+}
+
+template<typename TreeElemType>
+template<typename TreeType>
+void DiscreteHilbertValue<TreeElemType>::UpdateLargestValue(TreeType* node)
+{
+  if (!node->IsLeaf())
+  {
+    // Update the largest Hilbert value
+    localHilbertValues = node->Children()[node->NumChildren() -
+        1]->AuxiliaryInfo().HilbertValue().LocalHilbertValues();
+    numValues = node->Children()[node->NumChildren() -
+        1]->AuxiliaryInfo().HilbertValue().NumValues();
+  }
+}
+
+template<typename TreeElemType>
+template<typename TreeType>
+void DiscreteHilbertValue<TreeElemType>::RedistributeHilbertValues(
+    TreeType* parent,
+    const size_t firstSibling,
+    const size_t lastSibling)
+{
+  // We need to update the local dataset if points were redistributed.
+  size_t numPoints = 0;
+  for (size_t i = firstSibling; i <= lastSibling; i++)
+    numPoints += parent->Children()[i]->NumPoints();
+
+  // Copy the local Hilbert values.
+  arma::Mat<HilbertElemType> tmp(localHilbertValues->n_rows, numPoints);
+
+  size_t iPoint = 0;
+  for (size_t i = firstSibling; i<= lastSibling; i++)
+  {
+    DiscreteHilbertValue<TreeElemType> &value =
+        parent->Children()[i]->AuxiliaryInfo().HilbertValue();
+
+    for (size_t j = 0; j < value.NumValues(); j++)
+    {
+      tmp.col(iPoint) = value.LocalHilbertValues()->col(j);
+      iPoint++;
+    }
+  }
+  assert(iPoint == numPoints);
+
+  iPoint = 0;
+
+  // Redistribute the Hilbert values.
+  for (size_t i = firstSibling; i <= lastSibling; i++)
+  {
+    DiscreteHilbertValue<TreeElemType> &value =
+        parent->Children()[i]->AuxiliaryInfo().HilbertValue();
+
+    for (size_t j = 0; j < parent->Children()[i]->NumPoints(); j++)
+    {
+      value.LocalHilbertValues()->col(j) = tmp.col(iPoint);
+      iPoint++;
+    }
+    value.NumValues() = parent->Children()[i]->NumPoints();
+  }
+
+  assert(iPoint == numPoints);
+}
+
+template<typename TreeElemType>
+template<typename Archive>
+void DiscreteHilbertValue<TreeElemType>::
+Serialize(Archive& ar, const unsigned int /* version */)
+{
+  using data::CreateNVP;
+
+  ar & CreateNVP(localHilbertValues, "localHilbertValues");
+  ar & CreateNVP(ownsLocalHilbertValues, "ownsLocalHilbertValues");
+  ar & CreateNVP(numValues, "numValues");
+  ar & CreateNVP(valueToInsert, "valueToInsert");
+  ar & CreateNVP(ownsValueToInsert, "ownsValueToInsert");
+}
+
+} // namespace tree
+} // namespace mlpack
+
+#endif  //  MLPACK_CORE_TREE_RECTANGLE_TREE_DISCRETE_HILBERT_VALUE_IMPL_HPP

src/mlpack/core/tree/rectangle_tree/dual_tree_traverser.hpp

@@ -19,11 +19,12 @@ namespace tree {
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 template<typename RuleType>
 class RectangleTree<MetricType, StatisticType, MatType, SplitType,
-                    DescentType>::DualTreeTraverser
+                    DescentType, AuxiliaryInformationType>::DualTreeTraverser
 {
  public:
   /**

src/mlpack/core/tree/rectangle_tree/dual_tree_traverser_impl.hpp

@@ -20,10 +20,12 @@ namespace tree {
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 template<typename RuleType>
-RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+              AuxiliaryInformationType>::
 DualTreeTraverser<RuleType>::DualTreeTraverser(RuleType& rule) :
     rule(rule),
     numPrunes(0),
@@ -35,10 +37,12 @@ DualTreeTraverser<RuleType>::DualTreeTraverser(RuleType& rule) :
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 template<typename RuleType>
-void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
 DualTreeTraverser<RuleType>::Traverse(RectangleTree& queryNode,
                                       RectangleTree& referenceNode)
 {
@@ -63,14 +67,14 @@ DualTreeTraverser<RuleType>::Traverse(RectangleTree& queryNode,
     {
       // Restore the traversal information.
       rule.TraversalInfo() = traversalInfo;
-      const double childScore = rule.Score(queryNode.Points()[query],
+      const double childScore = rule.Score(queryNode.Point(query),
           referenceNode);
 
       if (childScore == DBL_MAX)
         continue;  // We don't require a search in this reference node.
 
       for(size_t ref = 0; ref < referenceNode.Count(); ++ref)
-        rule.BaseCase(queryNode.Points()[query], referenceNode.Points()[ref]);
+        rule.BaseCase(queryNode.Point(query), referenceNode.Point(ref));
 
       numBaseCases += referenceNode.Count();
     }

src/mlpack/core/tree/rectangle_tree/hilbert_r_tree_auxiliary_information.hpp

@@ -0,0 +1,120 @@
+/**
+ * @file hilbert_r_tree_auxiliary_information.hpp
+ * @author Mikhail Lozhnikov
+ *
+ * Definition of the HilbertRTreeAuxiliaryInformation class,
+ * a class that provides some Hilbert r-tree specific information
+ * about the nodes.
+ */
+#ifndef MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_AUXILIARY_INFORMATION_HPP
+#define MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_AUXILIARY_INFORMATION_HPP
+
+namespace mlpack {
+namespace tree {
+
+template<typename TreeType,
+         template<typename> class HilbertValueType>
+class HilbertRTreeAuxiliaryInformation
+{
+ public:
+  //! The element type held by the tree.
+  typedef typename TreeType::ElemType ElemType;
+  //! Default constructor
+  HilbertRTreeAuxiliaryInformation();
+
+  /**
+   * Construct this as an auxiliary information for the given node.
+   *
+   * @param node The node that stores this auxiliary information.
+   */
+  HilbertRTreeAuxiliaryInformation(const TreeType* node);
+
+  /**
+   * Create an auxiliary information object by copying from the other node.
+   *
+   * @param other The node from which the information will be copied.
+   */
+  HilbertRTreeAuxiliaryInformation(
+      const HilbertRTreeAuxiliaryInformation& other);
+
+  /**
+   * The Hilbert R tree requires to insert points according to their Hilbert
+   * value. This method should take care of it.  It returns false if it does
+   * nothing and true if it handles the insertion process.
+   *
+   * @param node The node in which the point is being inserted.
+   * @param point The number of the point being inserted.
+   */
+  bool HandlePointInsertion(TreeType* node, const size_t point);
+
+  /**
+   * The Hilbert R tree requires to insert nodes according to their Hilbert
+   * value. This method should take care of it.  It returns false if it does
+   * nothing and true if it handles the insertion process.
+   *
+   * @param node The node in which the nodeToInsert is being inserted.
+   * @param nodeToInsert The node being inserted.
+   * @param insertionLevel The level of the tree at which the nodeToInsert
+   *        should be inserted.
+   */
+  bool HandleNodeInsertion(TreeType* node,
+                           TreeType* nodeToInsert,
+                           bool insertionLevel);
+
+  /**
+   * The Hilbert R tree requires all points to be arranged according to their
+   * Hilbert value. This method should take care of saving this property after
+   * the deletion process.  It returns false if it does nothing and true if it
+   * handles the deletion process.
+   *
+   * @param node The node from which the point is being deleted.
+   * @param localIndex The index of the point being deleted.
+   */
+  bool HandlePointDeletion(TreeType* node, const size_t localIndex);
+
+  /**
+   * The Hilbert R tree requires all nodes to be arranged according to their
+   * Hilbert value. This method should take care of saving this property after
+   * the deletion process.  It returns false if it does nothing and true if it
+   * handles the deletion process.
+   *
+   * @param node The node from which the node is being deleted.
+   * @param nodeIndex The index of the node being deleted.
+   */
+  bool HandleNodeRemoval(TreeType* node, const size_t nodeIndex);
+
+  /**
+   * Update the auxiliary information in the node. The method returns true if
+   * the update should be propogated downward.
+   *
+   * @param node The node in which the auxiliary information being update.
+   */
+  bool UpdateAuxiliaryInfo(TreeType* node);
+
+  //! Clear memory.
+  void NullifyData();
+
+ private:
+  //! The largest Hilbert value of a point enclosed by the node.
+  HilbertValueType<ElemType> hilbertValue;
+
+ public:
+  //! Return the largest Hilbert value of a point covered by the node.
+  const HilbertValueType<ElemType>& HilbertValue() const
+  { return hilbertValue; }
+  //! Modify the largest Hilbert value of a point covered by the node.
+  HilbertValueType<ElemType>& HilbertValue() { return hilbertValue; }
+
+  /**
+   * Serialize the information.
+   */
+  template<typename Archive>
+  void Serialize(Archive& ar, const unsigned int /* version */);
+};
+
+} // namespace tree
+} // namespace mlpack
+
+#include "hilbert_r_tree_auxiliary_information_impl.hpp"
+
+#endif//MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_AUXILIARY_INFORMATION_HPP

src/mlpack/core/tree/rectangle_tree/hilbert_r_tree_auxiliary_information_impl.hpp

@@ -0,0 +1,169 @@
+/**
+ * @file hilbert_r_tree_auxiliary_information.hpp
+ * @author Mikhail Lozhnikov
+ *
+ * Implementation of the HilbertRTreeAuxiliaryInformation class, a class that
+ * provides some Hilbert r-tree specific information about the nodes.
+ */
+#ifndef MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_AUXILIARY_INFORMATION_IMPL_HPP
+#define MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_AUXILIARY_INFORMATION_IMPL_HPP
+
+#include "hilbert_r_tree_auxiliary_information.hpp"
+
+namespace mlpack {
+namespace tree {
+
+template<typename TreeType,
+         template<typename> class HilbertValueType>
+HilbertRTreeAuxiliaryInformation<TreeType, HilbertValueType>::
+HilbertRTreeAuxiliaryInformation()
+{ }
+
+template<typename TreeType,
+         template<typename> class HilbertValueType>
+HilbertRTreeAuxiliaryInformation<TreeType, HilbertValueType>::
+HilbertRTreeAuxiliaryInformation(const TreeType* node) :
+    hilbertValue(node)
+{ }
+
+template<typename TreeType,
+         template<typename> class HilbertValueType>
+HilbertRTreeAuxiliaryInformation<TreeType, HilbertValueType>::
+HilbertRTreeAuxiliaryInformation(
+    const HilbertRTreeAuxiliaryInformation& other) :
+    hilbertValue(other.HilbertValue())
+{ }
+
+template<typename TreeType,
+         template<typename> class HilbertValueType>
+bool HilbertRTreeAuxiliaryInformation<TreeType, HilbertValueType>::
+HandlePointInsertion(TreeType* node, const size_t point)
+{
+  if (node->IsLeaf())
+  {
+    // Get the position at which the point should be inserted, and then update
+    // the largest Hilbert value of the node.
+    size_t pos = hilbertValue.InsertPoint(node, node->Dataset().col(point));
+
+    // Move points.
+    for (size_t i = node->NumPoints(); i > pos; i--)
+      node->Point(i) = node->Point(i - 1);
+
+    // Insert the point.
+    node->Point(pos) = point;
+    node->Count()++;
+  }
+  else
+  {
+    // Calculate the Hilbert value.
+    hilbertValue.InsertPoint(node, node->Dataset().col(point));
+  }
+
+  return true;
+}
+
+template<typename TreeType,
+         template<typename> class HilbertValueType>
+bool HilbertRTreeAuxiliaryInformation<TreeType, HilbertValueType>::
+HandleNodeInsertion(TreeType* node, TreeType* nodeToInsert, bool insertionLevel)
+{
+  if (insertionLevel)
+  {
+    size_t pos;
+
+    // Find the best position for the node being inserted.
+    // The node should be inserted according to its Hilbert value.
+    for (pos = 0; pos < node->NumChildren(); pos++)
+      if (HilbertValueType<ElemType>::CompareValues(
+                 node->Children()[pos]->AuxiliaryInfo().HilbertValue(),
+                 nodeToInsert->AuxiliaryInfo().HilbertValue()) < 0)
+          break;
+
+    // Move nodes.
+    for (size_t i = node->NumChildren(); i > pos; i--)
+      node->Children()[i] = node->Children()[i - 1];
+
+    // Insert the node.
+    node->Children()[pos] = nodeToInsert;
+    nodeToInsert->Parent() = node;
+
+    // Update the largest Hilbert value.
+    hilbertValue.InsertNode(nodeToInsert);
+  }
+  else
+    hilbertValue.InsertNode(nodeToInsert); // Update the largest Hilbert value.
+
+  return true;
+}
+
+template<typename TreeType,
+         template<typename> class HilbertValueType>
+bool HilbertRTreeAuxiliaryInformation<TreeType, HilbertValueType>::
+HandlePointDeletion(TreeType* node, const size_t localIndex)
+{
+  // Update the largest Hilbert value.
+  hilbertValue.DeletePoint(node,localIndex);
+
+  for (size_t i = localIndex + 1; localIndex < node->NumPoints(); i++)
+    node->Point(i - 1) = node->Point(i);
+
+  node->NumPoints()--;
+  return true;
+}
+
+template<typename TreeType,
+         template<typename> class HilbertValueType>
+bool HilbertRTreeAuxiliaryInformation<TreeType, HilbertValueType>::
+HandleNodeRemoval(TreeType* node, const size_t nodeIndex)
+{
+  // Update the largest Hilbert value.
+  hilbertValue.RemoveNode(node,nodeIndex);
+
+  for (size_t i = nodeIndex + 1; nodeIndex < node->NumChildren(); i++)
+    node->Children()[i - 1] = node->Children()[i];
+
+  node->NumChildren()--;
+  return true;
+}
+
+template<typename TreeType,
+         template<typename> class HilbertValueType>
+bool HilbertRTreeAuxiliaryInformation<TreeType, HilbertValueType>::
+UpdateAuxiliaryInfo(TreeType* node)
+{
+  if (node->IsLeaf())  //  Should already be updated
+    return true;
+
+  TreeType* child = node->Children()[node->NumChildren() - 1];
+  if (hilbertValue.CompareWith(child->AuxiliaryInfo().HilbertValue()) < 0)
+  {
+    hilbertValue = node->AuxiliaryInfo().HilbertValue();
+    return true;
+  }
+  return false;
+}
+
+template<typename TreeType,
+         template<typename> class HilbertValueType>
+void HilbertRTreeAuxiliaryInformation<TreeType, HilbertValueType>::
+NullifyData()
+{
+  hilbertValue.NullifyData();
+}
+
+template<typename TreeType,
+         template<typename> class HilbertValueType>
+template<typename Archive>
+void HilbertRTreeAuxiliaryInformation<TreeType ,HilbertValueType>::
+Serialize(Archive& ar, const unsigned int /* version */)
+{
+  using data::CreateNVP;
+
+  ar & CreateNVP(hilbertValue, "hilbertValue");
+}
+
+
+} // namespace tree
+} // namespace mlpack
+
+#endif//MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_AUXILIARY_INFORMATION_IMPL_HPP

src/mlpack/core/tree/rectangle_tree/hilbert_r_tree_descent_heuristic.hpp

@@ -0,0 +1,53 @@
+/**
+ * @file hilbert_r_tree_descent_heuristic.hpp
+ * @author Mikhail Lozhnikov
+ *
+ * Definition of HilbertRTreeDescentHeuristic, a class that chooses the best
+ * child of a node in an R tree when inserting a new point.
+ */
+#ifndef MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_DESCENT_HEURISTIC_HPP
+#define MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_DESCENT_HEURISTIC_HPP
+
+#include <mlpack/core.hpp>
+
+namespace mlpack {
+namespace tree {
+
+/**
+ * This class chooses the best child of a node in a Hilbert R tree when
+ * inserting a new point.  This is done, in this class, by using the Hilbert
+ * value of the point to be inserted.
+ */
+class HilbertRTreeDescentHeuristic
+{
+ public:
+  /**
+   * Evaluate the node using a heuristic. Returns the number of the node with
+   * minimum largest Hilbert value that is greater than the Hilbert value of the
+   * point being inserted.
+   *
+   * @param node The node that is being evaluated.
+   * @param point The number of the point that is being inserted.
+   */
+  template<typename TreeType>
+  static size_t ChooseDescentNode(const TreeType* node, const size_t point);
+
+  /**
+   * Evaluate the node using a heuristic. Returns the number of the node with
+   * minimum largest Hilbert value that is greater than the largest Hilbert
+   * value of the point being inserted.
+   *
+   * @param node The node that is being evaluated.
+   * @param insertedNode The node that is being inserted.
+   */
+  template<typename TreeType>
+  static size_t ChooseDescentNode(const TreeType* node,
+                                  const TreeType* insertedNode);
+};
+
+} //  namespace tree
+} //  namespace mlpack
+
+#include "hilbert_r_tree_descent_heuristic_impl.hpp"
+
+#endif // MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_DESCENT_HEURISTIC_HPP

src/mlpack/core/tree/rectangle_tree/hilbert_r_tree_descent_heuristic_impl.hpp

@@ -0,0 +1,49 @@
+/**
+ * @file hilbert_r_tree_descent_heuristic_impl.hpp
+ * @author Mikhail Lozhnikov
+ *
+ * Implementation of HilbertRTreeDescentHeuristic, a class that chooses the best
+ * child of a node in an R tree when inserting a new point.
+ */
+#ifndef MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_DESCENT_HEURISTIC_IMPL_HPP
+#define MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_DESCENT_HEURISTIC_IMPL_HPP
+
+#include "hilbert_r_tree_descent_heuristic.hpp"
+
+namespace mlpack {
+namespace tree {
+
+template<typename TreeType>
+size_t HilbertRTreeDescentHeuristic::ChooseDescentNode(
+    const TreeType* node,
+    const size_t point)
+{
+  size_t bestIndex = 0;
+
+  for (bestIndex = 0; bestIndex < node->NumChildren() - 1; bestIndex++)
+    if (node->Children()[bestIndex]->AuxiliaryInfo().HilbertValue().
+        CompareWithCachedPoint(node->Dataset().col(point)) > 0)
+      break;
+
+  return bestIndex;
+}
+
+template<typename TreeType>
+size_t HilbertRTreeDescentHeuristic::ChooseDescentNode(
+    const TreeType* node,
+    const TreeType* insertedNode)
+{
+  size_t bestIndex = 0;
+
+  for (bestIndex = 0; bestIndex < node->NumChildren() - 1; bestIndex++)
+    if (node->Children()[bestIndex]->AuxiliaryInfo().HilbertValue().
+        CompareWith(node, node->AuxiliaryInfo().HilbertValue()) > 0)
+      break;
+
+  return bestIndex;
+}
+
+} //  namespace tree
+} //  namespace mlpack
+
+#endif // MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_DESCENT_HEURISTIC_IMPL_HPP

src/mlpack/core/tree/rectangle_tree/hilbert_r_tree_split.hpp

@@ -1 +1,95 @@
+/**
+ * @file hilbert_r_tree_split.hpp
+ * @author Mikhail Lozhnikov
+ *
+ * Defintion of the HilbertRTreeSplit class, a class that splits the nodes of an R
+ * tree, starting at a leaf node and moving upwards if necessary.
+ */
+#ifndef MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_SPLIT_HPP
+#define MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_SPLIT_HPP
+
+#include <mlpack/core.hpp>
+
+namespace mlpack {
+namespace tree /** Trees and tree-building procedures. */ {
+
+/**
+ * The splitting procedure for the Hilbert R tree.  The template parameter
+ * splitOrder is the order of the splitting policy. The Hilbert R tree splits a
+ * node on overflow, turning splitOrder nodes into (splitOrder + 1) nodes.
+ *
+ * @tparam splitOrder Number of nodes to split.
+ */
+template<size_t splitOrder = 2>
+class HilbertRTreeSplit
+{
+ public:
+  /**
+   * Split a leaf node using the "default" algorithm.  If necessary, this split
+   * will propagate upwards through the tree.
+   *
+   * @param node The node that is being split.
+   * @param relevels Not used.
+   */
+  template<typename TreeType>
+  static void SplitLeafNode(TreeType* tree, std::vector<bool>& relevels);
+
+  /**
+   * Split a non-leaf node using the "default" algorithm.  If this is a root
+   * node, the tree increases in depth.
+   *
+   * @param node The node that is being split.
+   * @param relevels Not used.
+   */
+  template<typename TreeType>
+  static bool SplitNonLeafNode(TreeType* tree, std::vector<bool>& relevels);
+
+ private:
+  /**
+   * Try to find splitOrder cooperating siblings in order to redistribute their
+   * children evenly. Returns true on success.
+   *
+   * @param parent The parent of of the overflowing node.
+   * @param iTree The number of the overflowing node.
+   * @param firstSibling The first cooperating sibling.
+   * @param lastSibling The last cooperating sibling.
+   */
+  template<typename TreeType>
+  static bool FindCooperatingSiblings(TreeType* parent,
+                                      const size_t iTree,
+                                      size_t& firstSibling,
+                                      size_t& lastSibling);
+
+  /**
+   * Redistribute the children of the cooperating siblings evenly among them.
+   *
+   * @param parent The parent of of the overflowing node.
+   * @param firstSibling The first cooperating sibling.
+   * @param lastSibling The last cooperating sibling.
+   */
+  template<typename TreeType>
+  static void RedistributeNodesEvenly(const TreeType* parent,
+                                      const size_t firstSibling,
+                                      const size_t lastSibling);
+
+  /**
+   * Redistribute the points of the cooperating siblings evenly among them.
+   *
+   * @param parent The parent of of the overflowing node.
+   * @param firstSibling The first cooperating sibling.
+   * @param lastSibling The last cooperating sibling.
+   */
+  template<typename TreeType>
+  static void RedistributePointsEvenly(TreeType* parent,
+                                       const size_t firstSibling,
+                                       const size_t lastSibling);
+};
+
+} // namespace tree
+} // namespace mlpack
+
+// Include implementation.
+#include "hilbert_r_tree_split_impl.hpp"
+
+#endif
 

src/mlpack/core/tree/rectangle_tree/hilbert_r_tree_split_impl.hpp

@@ -1 +1,336 @@
+/**
+ * @file hilbert_r_tree_split_impl.hpp
+ * @author Mikhail Lozhnikov
+ *
+ * Implementation of class (HilbertRTreeSplit) to split a RectangleTree.
+ */
+#ifndef MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_SPLIT_IMPL_HPP
+#define MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_SPLIT_IMPL_HPP
 
+#include "hilbert_r_tree_split.hpp"
+#include "rectangle_tree.hpp"
+#include <mlpack/core/math/range.hpp>
+
+namespace mlpack {
+namespace tree {
+
+template<size_t splitOrder>
+template<typename TreeType>
+void HilbertRTreeSplit<splitOrder>::SplitLeafNode(TreeType* tree,
+                                                  std::vector<bool>& relevels)
+{
+  // If we are splitting the root node, we need will do things differently so
+  // that the constructor and other methods don't confuse the end user by giving
+  // an address of another node.
+  if (tree->Parent() == NULL)
+  {
+    // We actually want to copy this way.  Pointers and everything.
+    TreeType* copy = new TreeType(*tree, false);
+    copy->Parent() = tree;
+    tree->Count() = 0;
+    tree->NullifyData();
+    // Because this was a leaf node, numChildren must be 0.
+    tree->Children()[(tree->NumChildren())++] = copy;
+    SplitLeafNode(copy, relevels);
+    return;
+  }
+
+  TreeType* parent = tree->Parent();
+  size_t iTree = 0;
+  for (iTree = 0; parent->Children()[iTree] != tree; iTree++);
+
+  // Try to find splitOrder cooperating siblings in order to redistribute points
+  // among them and avoid split.
+  size_t firstSibling, lastSibling;
+  if (FindCooperatingSiblings(parent, iTree, firstSibling, lastSibling))
+  {
+    RedistributePointsEvenly(parent, firstSibling, lastSibling);
+    return;
+  }
+
+  // We can not find splitOrder cooperating siblings since they are all full.
+  // We introduce new one instead.
+  size_t iNewSibling = (iTree + splitOrder < parent->NumChildren() ?
+                        iTree + splitOrder : parent->NumChildren());
+
+  for (size_t i = parent->NumChildren(); i > iNewSibling ; i--)
+    parent->Children()[i] = parent->Children()[i - 1];
+
+  parent->NumChildren()++;
+
+  parent->Children()[iNewSibling] = new TreeType(parent);
+
+  lastSibling = (iTree + splitOrder < parent->NumChildren() ?
+                 iTree + splitOrder : parent->NumChildren() - 1);
+  firstSibling = (lastSibling > splitOrder ? lastSibling - splitOrder : 0);
+
+  assert(lastSibling - firstSibling <= splitOrder);
+  assert(firstSibling >= 0);
+  assert(lastSibling < parent->NumChildren());
+
+  // Redistribute the points among (splitOrder + 1) cooperating siblings evenly.
+  RedistributePointsEvenly(parent, firstSibling, lastSibling);
+
+  if (parent->NumChildren() == parent->MaxNumChildren() + 1)
+    SplitNonLeafNode(parent, relevels);
+}
+
+template<size_t splitOrder>
+template<typename TreeType>
+bool HilbertRTreeSplit<splitOrder>::
+SplitNonLeafNode(TreeType* tree, std::vector<bool>& relevels)
+{
+  // If we are splitting the root node, we need will do things differently so
+  // that the constructor and other methods don't confuse the end user by giving
+  // an address of another node.
+  if (tree->Parent() == NULL)
+  {
+    // We actually want to copy this way.  Pointers and everything.
+    TreeType* copy = new TreeType(*tree, false);
+
+    copy->Parent() = tree;
+    tree->NumChildren() = 0;
+    tree->NullifyData();
+    tree->Children()[(tree->NumChildren())++] = copy;
+
+    SplitNonLeafNode(copy, relevels);
+    return true;
+  }
+
+  TreeType* parent = tree->Parent();
+
+  size_t iTree = 0;
+  for (iTree = 0; parent->Children()[iTree] != tree; iTree++);
+
+  // Try to find splitOrder cooperating siblings in order to redistribute
+  // children among them and avoid split.
+  size_t firstSibling, lastSibling;
+  if (FindCooperatingSiblings(parent, iTree, firstSibling, lastSibling))
+  {
+    RedistributeNodesEvenly(parent, firstSibling, lastSibling);
+    return false;
+  }
+
+  // We can not find splitOrder cooperating siblings since they are all full.
+  // We introduce new one instead.
+  size_t iNewSibling = (iTree + splitOrder < parent->NumChildren() ?
+                        iTree + splitOrder : parent->NumChildren());
+
+  for (size_t i = parent->NumChildren(); i > iNewSibling ; i--)
+    parent->Children()[i] = parent->Children()[i - 1];
+
+  parent->NumChildren()++;
+
+  parent->Children()[iNewSibling] = new TreeType(parent);
+
+  lastSibling = (iTree + splitOrder < parent->NumChildren() ?
+                 iTree + splitOrder : parent->NumChildren() - 1);
+  firstSibling = (lastSibling > splitOrder ?
+                  lastSibling - splitOrder : 0);
+
+  assert(lastSibling - firstSibling <= splitOrder);
+  assert(firstSibling >= 0);
+  assert(lastSibling < parent->NumChildren());
+
+  // Redistribute children among (splitOrder + 1) cooperating siblings evenly.
+  RedistributeNodesEvenly(parent, firstSibling, lastSibling);
+
+  if (parent->NumChildren() == parent->MaxNumChildren() + 1)
+    SplitNonLeafNode(parent, relevels);
+  return false;
+}
+
+template<size_t splitOrder>
+template<typename TreeType>
+bool HilbertRTreeSplit<splitOrder>::FindCooperatingSiblings(
+    TreeType* parent,
+    const size_t iTree,
+    size_t& firstSibling,
+    size_t& lastSibling)
+{
+  size_t start = (iTree > splitOrder - 1 ? iTree - splitOrder + 1 : 0);
+  size_t end = (iTree + splitOrder <= parent->NumChildren() ?
+                iTree + splitOrder : parent->NumChildren());
+
+  size_t iUnderfullSibling;
+
+  // Try to find empty space among cooperating siblings.
+  if (parent->Children()[iTree]->NumChildren() != 0)
+  {
+    for (iUnderfullSibling = start; iUnderfullSibling < end;
+         iUnderfullSibling++)
+      if (parent->Children()[iUnderfullSibling]->NumChildren() <
+          parent->Children()[iUnderfullSibling]->MaxNumChildren() - 1)
+        break;
+  }
+  else
+  {
+    for (iUnderfullSibling = start; iUnderfullSibling < end;
+         iUnderfullSibling++)
+      if (parent->Children()[iUnderfullSibling]->NumPoints() <
+          parent->Children()[iUnderfullSibling]->MaxLeafSize() - 1)
+        break;
+  }
+
+  if (iUnderfullSibling == end) // All nodes are full.
+    return false;
+
+  if (iUnderfullSibling > iTree)
+  {
+    lastSibling = (iTree + splitOrder - 1 < parent->NumChildren() ?
+                   iTree + splitOrder - 1 : parent->NumChildren() - 1);
+    firstSibling = (lastSibling > splitOrder - 1 ?
+                    lastSibling - splitOrder + 1 : 0);
+  }
+  else
+  {
+    lastSibling = (iUnderfullSibling + splitOrder - 1 < parent->NumChildren() ?
+        iUnderfullSibling + splitOrder - 1 : parent->NumChildren() - 1);
+    firstSibling = (lastSibling > splitOrder - 1 ?
+        lastSibling - splitOrder + 1 : 0);
+  }
+
+  assert(lastSibling - firstSibling <= splitOrder - 1);
+  assert(firstSibling >= 0);
+  assert(lastSibling < parent->NumChildren());
+
+  return true;
+}
+
+template<size_t splitOrder>
+template<typename TreeType>
+void HilbertRTreeSplit<splitOrder>::
+RedistributeNodesEvenly(const TreeType *parent,
+                        size_t firstSibling, size_t lastSibling)
+{
+  size_t numChildren = 0;
+  size_t numChildrenPerNode, numRestChildren;
+
+  for (size_t i = firstSibling; i <= lastSibling; i++)
+    numChildren += parent->Children()[i]->NumChildren();
+
+  numChildrenPerNode = numChildren / (lastSibling - firstSibling + 1);
+  numRestChildren = numChildren % (lastSibling - firstSibling + 1);
+
+  std::vector<TreeType*> children(numChildren);
+
+  // Copy children's children in order to redistribute them.
+  size_t iChild = 0;
+  for (size_t i = firstSibling; i <= lastSibling; i++)
+  {
+    for (size_t j = 0; j < parent->Children()[i]->NumChildren(); j++)
+    {
+      children[iChild] = parent->Children()[i]->Children()[j];
+      iChild++;
+    }
+  }
+
+  iChild = 0;
+  for (size_t i = firstSibling; i <= lastSibling; i++)
+  {
+    // Since we redistribute children of a sibling we should recalculate the
+    // bound.
+    parent->Children()[i]->Bound().Clear();
+
+    for (size_t j = 0; j < numChildrenPerNode; j++)
+    {
+      parent->Children()[i]->Bound() |= children[iChild]->Bound();
+      parent->Children()[i]->Children()[j] = children[iChild];
+      children[iChild]->Parent() = parent->Children()[i];
+      iChild++;
+    }
+    if (numRestChildren > 0)
+    {
+      parent->Children()[i]->Bound() |= children[iChild]->Bound();
+      parent->Children()[i]->Children()[numChildrenPerNode] = children[iChild];
+      children[iChild]->Parent() = parent->Children()[i];
+      parent->Children()[i]->NumChildren() = numChildrenPerNode + 1;
+      numRestChildren--;
+      iChild++;
+    }
+    else
+    {
+      parent->Children()[i]->NumChildren() = numChildrenPerNode;
+    }
+    assert(parent->Children()[i]->NumChildren() <=
+           parent->Children()[i]->MaxNumChildren());
+
+    // Fix the largest Hilbert value of the sibling.
+    parent->Children()[i]->AuxiliaryInfo().HilbertValue().UpdateLargestValue(
+        parent->Children()[i]);
+  }
+}
+
+template<size_t splitOrder>
+template<typename TreeType>
+void HilbertRTreeSplit<splitOrder>::
+RedistributePointsEvenly(TreeType* parent,
+                         const size_t firstSibling,
+                         const size_t lastSibling)
+{
+  size_t numPoints = 0;
+  size_t numPointsPerNode, numRestPoints;
+
+  for (size_t i = firstSibling; i <= lastSibling; i++)
+    numPoints += parent->Children()[i]->NumPoints();
+
+  numPointsPerNode = numPoints / (lastSibling - firstSibling + 1);
+  numRestPoints = numPoints % (lastSibling - firstSibling + 1);
+
+  std::vector<size_t> points(numPoints);
+
+  // Copy children's points in order to redistribute them.
+  size_t iPoint = 0;
+  for (size_t i = firstSibling; i <= lastSibling; i++)
+  {
+    for (size_t j = 0; j < parent->Children()[i]->NumPoints(); j++)
+      points[iPoint++] = parent->Children()[i]->Point(j);
+  }
+
+  iPoint = 0;
+  for (size_t i = firstSibling; i <= lastSibling; i++)
+  {
+    // Since we redistribute points of a sibling we should recalculate the
+    // bound.
+    parent->Children()[i]->Bound().Clear();
+
+    size_t j;
+    for (j = 0; j < numPointsPerNode; j++)
+    {
+      parent->Children()[i]->Bound() |= parent->Dataset().col(points[iPoint]);
+      parent->Children()[i]->Point(j) = points[iPoint];
+      iPoint++;
+    }
+    if (numRestPoints > 0)
+    {
+      parent->Children()[i]->Bound() |= parent->Dataset().col(points[iPoint]);
+      parent->Children()[i]->Point(j) = points[iPoint];
+      parent->Children()[i]->Count() = numPointsPerNode + 1;
+      numRestPoints--;
+      iPoint++;
+    }
+    else
+    {
+      parent->Children()[i]->Count() = numPointsPerNode;
+    }
+    assert(parent->Children()[i]->NumPoints() <=
+           parent->Children()[i]->MaxLeafSize());
+  }
+
+  // Fix the largest Hilbert values of the siblings.
+  parent->AuxiliaryInfo().HilbertValue().RedistributeHilbertValues(parent,
+      firstSibling, lastSibling);
+
+  TreeType* root = parent;
+
+  while (root != NULL)
+  {
+    root->AuxiliaryInfo().HilbertValue().UpdateLargestValue(root);
+    root = root->Parent();
+  }
+}
+
+} // namespace tree
+} // namespace mlpack
+
+#endif  //  MLPACK_CORE_TREE_RECTANGLE_TREE_HILBERT_R_TREE_SPLIT_IMPL_HPP

src/mlpack/core/tree/rectangle_tree/no_auxiliary_information.hpp

@@ -0,0 +1,117 @@
+/**
+ * @file no_auxiliary_information.hpp
+ * @author Mikhail Lozhnikov
+ *
+ * Definition of the NoAuxiliaryInformation class, a class that provides
+ * no additional information about the nodes.
+ */
+#ifndef MLPACK_CORE_TREE_RECTANGLE_TREE_NO_AUXILIARY_INFORMATION_HPP
+#define MLPACK_CORE_TREE_RECTANGLE_TREE_NO_AUXILIARY_INFORMATION_HPP
+
+namespace mlpack {
+namespace tree {
+
+template<typename TreeType>
+class NoAuxiliaryInformation
+{
+ public:
+  //! Construct the auxiliary information object.
+  NoAuxiliaryInformation() { };
+  //! Construct the auxiliary information object.
+  NoAuxiliaryInformation(const TreeType* /* node */) { };
+  //! Construct the auxiliary information object.
+  NoAuxiliaryInformation(const TreeType& /* node */) { };
+
+  /**
+   * Some tree types require to save some properties at the insertion process.
+   * This method allows the auxiliary information the option of manipulating
+   * the tree in order to perform the insertion process. If the auxiliary
+   * information does that, then the method should return true; if the method
+   * returns false the RectangleTree performs its default behavior.
+   *
+   * @param node The node in which the point is being inserted.
+   * @param point The global number of the point being inserted.
+   */
+  bool HandlePointInsertion(TreeType* , const size_t)
+  {
+    return false;
+  }
+
+  /**
+   * Some tree types require to save some properties at the insertion process.
+   * This method allows the auxiliary information the option of manipulating
+   * the tree in order to perform the insertion process. If the auxiliary
+   * information does that, then the method should return true; if the method
+   * returns false the RectangleTree performs its default behavior.
+   *
+   * @param node The node in which the nodeToInsert is being inserted.
+   * @param nodeToInsert The node being inserted.
+   * @param insertionLevel The level of the tree at which the nodeToInsert
+   *        should be inserted.
+   */
+  bool HandleNodeInsertion(TreeType* /* node */,
+                           TreeType* /* nodeToInsert */,
+                           bool /* insertionLevel */)
+  {
+    return false;
+  }
+
+  /**
+   * Some tree types require to save some properties at the deletion process.
+   * This method allows the auxiliary information the option of manipulating
+   * the tree in order to perform the deletion process. If the auxiliary
+   * information does that, then the method should return true; if the method
+   * returns false the RectangleTree performs its default behavior.
+   *
+   * @param node The node from which the point is being deleted.
+   * @param localIndex The local index of the point being deleted.
+   */
+  bool HandlePointDeletion(TreeType* /* node */, const size_t /* localIndex */)
+  {
+    return false;
+  }
+
+  /**
+   * Some tree types require to save some properties at the deletion process.
+   * This method allows the auxiliary information the option of manipulating
+   * the tree in order to perform the deletion process. If the auxiliary
+   * information does that, then the method should return true; if the method
+   * returns false the RectangleTree performs its default behavior.
+   *
+   * @param node The node from which the node is being deleted.
+   * @param nodeIndex The local index of the node being deleted.
+   */
+  bool HandleNodeRemoval(TreeType* /* node */, const size_t /* nodeIndex */)
+  {
+    return false;
+  }
+
+  /**
+   * Some tree types require to propagate the information upward.
+   * This method should return false if this is not the case. If true is
+   * returned, the update will be propogated upward.
+   *
+   * @param node The node in which the auxiliary information being update.
+   */
+  bool UpdateAuxiliaryInfo(TreeType* /* node */)
+  {
+    return false;
+  }
+
+  /**
+   * Nullify the auxiliary information in order to prevent an invalid free.
+   */
+  void NullifyData()
+  { }
+
+  /**
+   * Serialize the information.
+   */
+  template<typename Archive>
+  void Serialize(Archive &, const unsigned int /* version */) { };
+};
+
+} // namespace tree
+} // namespace mlpack
+
+#endif  //  MLPACK_CORE_TREE_RECTANGLE_TREE_NO_AUXILIARY_INFORMATION_HPP

src/mlpack/core/tree/rectangle_tree/r_star_tree_descent_heuristic.hpp

@@ -14,25 +14,25 @@ namespace mlpack {
 namespace tree {
 
 /**
- * When descending a Rectangle tree to insert a point, we need to have a way to
+ * When descending a RectangleTree to insert a point, we need to have a way to
  * choose a child node when the point isn't enclosed by any of them.  This
- * heuristic is used to do so.
+ * heuristic is used to do so using the rules for the R* tree.
  */
 class RStarTreeDescentHeuristic
 {
  public:
   /**
-   * Evaluate the node using a hueristic.  The heuristic guarantees two things:
+   * Evaluate the node using a heuristic.  The heuristic guarantees two things:
    *
-   * 1.  If point is contained in (or on) bound, the value returned is zero.
-   * 2.  If the point is not contained in (or on) bound, the value returned is
-   *     greater than zero.
+   * 1. If point is contained in (or on) bound, the value returned is zero.
+   * 2. If the point is not contained in (or on) bound, the value returned is
+   *    greater than zero.
    *
    * @param bound The bound used for the node that is being evaluated.
-   * @param point The point that is being inserted.
+   * @param point The index of the point that is being inserted.
    */
   template<typename TreeType>
-  static size_t ChooseDescentNode(const TreeType* node, const arma::vec& point);
+  static size_t ChooseDescentNode(const TreeType* node, const size_t point);
 
   template<typename TreeType>
   static size_t ChooseDescentNode(const TreeType* node,

src/mlpack/core/tree/rectangle_tree/r_star_tree_descent_heuristic_impl.hpp

@@ -16,7 +16,7 @@ namespace tree {
 template<typename TreeType>
 inline size_t RStarTreeDescentHeuristic::ChooseDescentNode(
     const TreeType* node,
-    const arma::vec& point)
+    const size_t point)
 {
   // Convenience typedef.
   typedef typename TreeType::ElemType ElemType;
@@ -41,9 +41,9 @@ inline size_t RStarTreeDescentHeuristic::ChooseDescentNode(
           ElemType newOverlap = 1.0;
           for (size_t k = 0; k < node->Bound().Dim(); k++)
           {
-            ElemType newHigh = std::max(point[k],
+            ElemType newHigh = std::max(node->Dataset().col(point)[k],
                 node->Children()[i]->Bound()[k].Hi());
-            ElemType newLow = std::min(point[k],
+            ElemType newLow = std::min(node->Dataset().col(point)[k],
                 node->Children()[i]->Bound()[k].Lo());
             overlap *= node->Children()[i]->Bound()[k].Hi() < node->Children()[j]->Bound()[k].Lo() || node->Children()[i]->Bound()[k].Lo() > node->Children()[j]->Bound()[k].Hi() ? 0 : std::min(node->Children()[i]->Bound()[k].Hi(), node->Children()[j]->Bound()[k].Hi()) - std::max(node->Children()[i]->Bound()[k].Lo(), node->Children()[j]->Bound()[k].Lo());
             newOverlap *= newHigh < node->Children()[j]->Bound()[k].Lo() || newLow > node->Children()[j]->Bound()[k].Hi() ? 0 : std::min(newHigh, node->Children()[j]->Bound()[k].Hi()) - std::max(newLow, node->Children()[j]->Bound()[k].Lo());
@@ -91,8 +91,8 @@ inline size_t RStarTreeDescentHeuristic::ChooseDescentNode(
       for (size_t j = 0; j < node->Bound().Dim(); j++)
       {
         v1 *= node->Children()[i]->Bound()[j].Width();
-        v2 *= node->Children()[i]->Bound()[j].Contains(point[j]) ? node->Children()[i]->Bound()[j].Width() : (node->Children()[i]->Bound()[j].Hi() < point[j] ? (point[j] - node->Children()[i]->Bound()[j].Lo()) :
-            (node->Children()[i]->Bound()[j].Hi() - point[j]));
+        v2 *= node->Children()[i]->Bound()[j].Contains(node->Dataset().col(point)[j]) ? node->Children()[i]->Bound()[j].Width() : (node->Children()[i]->Bound()[j].Hi() < node->Dataset().col(point)[j] ? (node->Dataset().col(point)[j] - node->Children()[i]->Bound()[j].Lo()) :
+            (node->Children()[i]->Bound()[j].Hi() - node->Dataset().col(point)[j]));
       }
 
       assert(v2 - v1 >= 0);

src/mlpack/core/tree/rectangle_tree/r_star_tree_split.hpp

@@ -18,31 +18,23 @@ namespace tree /** Trees and tree-building procedures. */ {
  * nodes overflow, we split them, moving up the tree and splitting nodes
  * as necessary.
  */
-template <typename TreeType>
 class RStarTreeSplit
 {
  public:
-  //! Default constructor.
-  RStarTreeSplit() { }
-
-  //! Construct this with the specified node.
-  RStarTreeSplit(const TreeType* /* node */) { }
-
-  //! Create a copy of the other.split.
-  RStarTreeSplit(const TreeType& /* other */) { }
-
   /**
    * Split a leaf node using the algorithm described in "The R*-tree: An
    * Efficient and Robust Access method for Points and Rectangles."  If
    * necessary, this split will propagate upwards through the tree.
    */
-  void SplitLeafNode(TreeType* tree, std::vector<bool>& relevels);
+  template <typename TreeType>
+  static void SplitLeafNode(TreeType *tree,std::vector<bool>& relevels);
 
   /**
    * Split a non-leaf node using the "default" algorithm.  If this is a root
    * node, the tree increases in depth.
    */
-  bool SplitNonLeafNode(TreeType* tree, std::vector<bool>& relevels);
+  template <typename TreeType>
+  static bool SplitNonLeafNode(TreeType *tree,std::vector<bool>& relevels);
 
  private:
   /**
@@ -68,14 +60,8 @@ class RStarTreeSplit
   /**
    * Insert a node into another node.
    */
+  template <typename TreeType>
   static void InsertNodeIntoTree(TreeType* destTree, TreeType* srcNode);
-
- public:
-  /**
-   * Serialize the split.
-   */
-  template<typename Archive>
-  void Serialize(Archive &, const unsigned int /* version */) { };
 };
 
 } // namespace tree

src/mlpack/core/tree/rectangle_tree/r_star_tree_split_impl.hpp

@@ -22,7 +22,7 @@ namespace tree {
  * new nodes into the tree, spliting the parent if necessary.
  */
 template<typename TreeType>
-void RStarTreeSplit<TreeType>::SplitLeafNode(TreeType *tree,std::vector<bool>& relevels)
+void RStarTreeSplit::SplitLeafNode(TreeType *tree,std::vector<bool>& relevels)
 {
   // Convenience typedef.
   typedef typename TreeType::ElemType ElemType;
@@ -41,7 +41,7 @@ void RStarTreeSplit<TreeType>::SplitLeafNode(TreeType *tree,std::vector<bool>& r
     tree->Children()[(tree->NumChildren())++] = copy;
     assert(tree->NumChildren() == 1);
 
-    copy->Split().SplitLeafNode(copy, relevels);
+    RStarTreeSplit::SplitLeafNode(copy,relevels);
     return;
   }
 
@@ -53,12 +53,12 @@ void RStarTreeSplit<TreeType>::SplitLeafNode(TreeType *tree,std::vector<bool>& r
     // We sort the points by decreasing distance to the centroid of the bound.
     // We then remove the first p entries and reinsert them at the root.
     TreeType* root = tree;
-    while(root->Parent() != NULL)
+    while (root->Parent() != NULL)
       root = root->Parent();
     size_t p = tree->MaxLeafSize() * 0.3; // The paper says this works the best.
     if (p == 0)
     {
-      tree->Split().SplitLeafNode(tree, relevels);
+      RStarTreeSplit::SplitLeafNode(tree,relevels);
       return;
     }
 
@@ -68,17 +68,19 @@ void RStarTreeSplit<TreeType>::SplitLeafNode(TreeType *tree,std::vector<bool>& r
     for (size_t i = 0; i < sorted.size(); i++)
     {
       sorted[i].d = tree->Metric().Evaluate(center,
-          tree->LocalDataset().col(i));
+          tree->Dataset().col(tree->Point(i)));
       sorted[i].n = i;
     }
 
     std::sort(sorted.begin(), sorted.end(), StructComp<ElemType>);
-    std::vector<int> pointIndices(p);
+    std::vector<size_t> pointIndices(p);
+
     for (size_t i = 0; i < p; i++)
     {
       // We start from the end of sorted.
-      pointIndices[i] = tree->Points()[sorted[sorted.size() - 1 - i].n];
-      root->DeletePoint(tree->Points()[sorted[sorted.size() - 1 - i].n],
+      pointIndices[i] = tree->Point(sorted[sorted.size() - 1 - i].n);
+
+      root->DeletePoint(tree->Point(sorted[sorted.size() - 1 - i].n),
           relevels);
     }
 
@@ -104,7 +106,7 @@ void RStarTreeSplit<TreeType>::SplitLeafNode(TreeType *tree,std::vector<bool>& r
     std::vector<SortStruct<ElemType>> sorted(tree->Count());
     for (size_t i = 0; i < sorted.size(); i++)
     {
-      sorted[i].d = tree->LocalDataset().col(i)[j];
+      sorted[i].d = tree->Dataset().col(tree->Point(i))[j];
       sorted[i].n = i;
     }
 
@@ -140,25 +142,25 @@ void RStarTreeSplit<TreeType>::SplitLeafNode(TreeType *tree,std::vector<bool>& r
       std::vector<ElemType> minG2(maxG1.size());
       for (size_t k = 0; k < tree->Bound().Dim(); k++)
       {
-        minG1[k] = maxG1[k] = tree->LocalDataset().col(sorted[0].n)[k];
+        minG1[k] = maxG1[k] = tree->Dataset().col(tree->Point(sorted[0].n))[k];
         minG2[k] = maxG2[k] =
-            tree->LocalDataset().col(sorted[sorted.size() - 1].n)[k];
+            tree->Dataset().col(tree->Point(sorted[sorted.size() - 1].n))[k];
 
         for (size_t l = 1; l < tree->Count() - 1; l++)
         {
           if (l < cutOff)
           {
-            if (tree->LocalDataset().col(sorted[l].n)[k] < minG1[k])
-              minG1[k] = tree->LocalDataset().col(sorted[l].n)[k];
-            else if (tree->LocalDataset().col(sorted[l].n)[k] > maxG1[k])
-              maxG1[k] = tree->LocalDataset().col(sorted[l].n)[k];
+            if (tree->Dataset().col(tree->Point(sorted[l].n))[k] < minG1[k])
+              minG1[k] = tree->Dataset().col(tree->Point(sorted[l].n))[k];
+            else if (tree->Dataset().col(tree->Point(sorted[l].n))[k] > maxG1[k])
+              maxG1[k] = tree->Dataset().col(tree->Point(sorted[l].n))[k];
           }
           else
           {
-            if (tree->LocalDataset().col(sorted[l].n)[k] < minG2[k])
-              minG2[k] = tree->LocalDataset().col(sorted[l].n)[k];
-            else if (tree->LocalDataset().col(sorted[l].n)[k] > maxG2[k])
-              maxG2[k] = tree->LocalDataset().col(sorted[l].n)[k];
+            if (tree->Dataset().col(tree->Point(sorted[l].n))[k] < minG2[k])
+              minG2[k] = tree->Dataset().col(tree->Point(sorted[l].n))[k];
+            else if (tree->Dataset().col(tree->Point(sorted[l].n))[k] > maxG2[k])
+              maxG2[k] = tree->Dataset().col(tree->Point(sorted[l].n))[k];
           }
         }
       }
@@ -208,7 +210,7 @@ void RStarTreeSplit<TreeType>::SplitLeafNode(TreeType *tree,std::vector<bool>& r
   std::vector<SortStruct<ElemType>> sorted(tree->Count());
   for (size_t i = 0; i < sorted.size(); i++)
   {
-    sorted[i].d = tree->LocalDataset().col(i)[bestAxis];
+    sorted[i].d = tree->Dataset().col(tree->Point(i))[bestAxis];
     sorted[i].n = i;
   }
 
@@ -222,9 +224,9 @@ void RStarTreeSplit<TreeType>::SplitLeafNode(TreeType *tree,std::vector<bool>& r
     for (size_t i = 0; i < tree->Count(); i++)
     {
       if (i < bestAreaIndexOnBestAxis + tree->MinLeafSize())
-        treeOne->InsertPoint(tree->Points()[sorted[i].n]);
+        treeOne->InsertPoint(tree->Point(sorted[i].n));
       else
-        treeTwo->InsertPoint(tree->Points()[sorted[i].n]);
+        treeTwo->InsertPoint(tree->Point(sorted[i].n));
     }
   }
   else
@@ -232,9 +234,9 @@ void RStarTreeSplit<TreeType>::SplitLeafNode(TreeType *tree,std::vector<bool>& r
     for (size_t i = 0; i < tree->Count(); i++)
     {
       if (i < bestOverlapIndexOnBestAxis + tree->MinLeafSize())
-        treeOne->InsertPoint(tree->Points()[sorted[i].n]);
+        treeOne->InsertPoint(tree->Point(sorted[i].n));
       else
-        treeTwo->InsertPoint(tree->Points()[sorted[i].n]);
+        treeTwo->InsertPoint(tree->Point(sorted[i].n));
     }
   }
 
@@ -251,7 +253,7 @@ void RStarTreeSplit<TreeType>::SplitLeafNode(TreeType *tree,std::vector<bool>& r
   // just in case, we use an assert.
   assert(par->NumChildren() <= par->MaxNumChildren() + 1);
   if (par->NumChildren() == par->MaxNumChildren() + 1)
-    par->Split().SplitNonLeafNode(par, relevels);
+    RStarTreeSplit::SplitNonLeafNode(par,relevels);
 
   assert(treeOne->Parent()->NumChildren() <= treeOne->MaxNumChildren());
   assert(treeOne->Parent()->NumChildren() >= treeOne->MinNumChildren());
@@ -269,8 +271,7 @@ void RStarTreeSplit<TreeType>::SplitLeafNode(TreeType *tree,std::vector<bool>& r
  * higher up the tree because they were already updated if necessary.
  */
 template<typename TreeType>
-bool RStarTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
-                                                std::vector<bool>& relevels)
+bool RStarTreeSplit::SplitNonLeafNode(TreeType *tree,std::vector<bool>& relevels)
 {
   // Convenience typedef.
   typedef typename TreeType::ElemType ElemType;
@@ -288,7 +289,7 @@ bool RStarTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
     tree->NullifyData();
     tree->Children()[(tree->NumChildren())++] = copy;
 
-    copy->Split().SplitNonLeafNode(copy, relevels);
+    RStarTreeSplit::SplitNonLeafNode(copy,relevels);
     return true;
   }
 
@@ -643,9 +644,7 @@ bool RStarTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
   // just in case, we use an assert.
   assert(par->NumChildren() <= par->MaxNumChildren() + 1);
   if (par->NumChildren() == par->MaxNumChildren() + 1)
-  {
-    par->Split().SplitNonLeafNode(par, relevels);
-  }
+    RStarTreeSplit::SplitNonLeafNode(par,relevels);
 
   // We have to update the children of each of these new nodes so that they
   // record the correct parent.
@@ -673,8 +672,7 @@ bool RStarTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
  * numberOfChildren.
  */
 template<typename TreeType>
-void RStarTreeSplit<TreeType>::InsertNodeIntoTree(TreeType* destTree,
-                                                  TreeType* srcNode)
+void RStarTreeSplit::InsertNodeIntoTree(TreeType* destTree, TreeType* srcNode)
 {
   destTree->Bound() |= srcNode->Bound();
   destTree->Children()[destTree->NumChildren()++] = srcNode;

src/mlpack/core/tree/rectangle_tree/r_tree_descent_heuristic.hpp

@@ -29,10 +29,10 @@ class RTreeDescentHeuristic
    *    is greater than zero.
    *
    * @param node The node that is being evaluated.
-   * @param point The point that is being inserted.
+   * @param point The index of the point that is being inserted.
    */
   template<typename TreeType>
-  static size_t ChooseDescentNode(const TreeType* node, const arma::vec& point);
+  static size_t ChooseDescentNode(const TreeType* node, const size_t point);
 
   /**
    * Evaluate the node using a heuristic.  The heuristic guarantees two things:

src/mlpack/core/tree/rectangle_tree/r_tree_descent_heuristic_impl.hpp

@@ -15,7 +15,7 @@ namespace tree {
 
 template<typename TreeType>
 inline size_t RTreeDescentHeuristic::ChooseDescentNode(const TreeType* node,
-                                                       const arma::vec& point)
+                                                       const size_t point)
 {
   // Convenience typedef.
   typedef typename TreeType::ElemType ElemType;
@@ -31,11 +31,11 @@ inline size_t RTreeDescentHeuristic::ChooseDescentNode(const TreeType* node,
     for (size_t j = 0; j < node->Children()[i]->Bound().Dim(); j++)
     {
       v1 *= node->Children()[i]->Bound()[j].Width();
-      v2 *= node->Children()[i]->Bound()[j].Contains(point[j]) ?
+      v2 *= node->Children()[i]->Bound()[j].Contains(node->Dataset().col(point)[j]) ?
           node->Children()[i]->Bound()[j].Width() :
-          (node->Children()[i]->Bound()[j].Hi() < point[j] ?
-              (point[j] - node->Children()[i]->Bound()[j].Lo()) :
-              (node->Children()[i]->Bound()[j].Hi() - point[j]));
+          (node->Children()[i]->Bound()[j].Hi() < node->Dataset().col(point)[j] ?
+              (node->Dataset().col(point)[j] - node->Children()[i]->Bound()[j].Lo()) :
+              (node->Children()[i]->Bound()[j].Hi() - node->Dataset().col(point)[j]));
     }
 
     assert(v2 - v1 >= 0);

src/mlpack/core/tree/rectangle_tree/r_tree_split.hpp

@@ -18,45 +18,40 @@ namespace tree /** Trees and tree-building procedures. */ {
  * nodes overflow, we split them, moving up the tree and splitting nodes
  * as necessary.
  */
-template<typename TreeType>
 class RTreeSplit
 {
  public:
-  //! Default constructor.
-  RTreeSplit() { }
-
-  //! Construct this with the specified node.
-  RTreeSplit(const TreeType* /* node */) { }
-
-  //! Create a copy of the other split.
-  RTreeSplit(const TreeType& /* other */) { }
-
   /**
    * Split a leaf node using the "default" algorithm.  If necessary, this split
    * will propagate upwards through the tree.
    */
-  void SplitLeafNode(TreeType* tree, std::vector<bool>& relevels);
+  template<typename TreeType>
+  static void SplitLeafNode(TreeType *tree,std::vector<bool>& relevels);
 
   /**
    * Split a non-leaf node using the "default" algorithm.  If this is a root
    * node, the tree increases in depth.
    */
-  bool SplitNonLeafNode(TreeType* tree, std::vector<bool>& relevels);
+  template<typename TreeType>
+  static bool SplitNonLeafNode(TreeType *tree,std::vector<bool>& relevels);
 
  private:
   /**
    * Get the seeds for splitting a leaf node.
    */
-  static void GetPointSeeds(const TreeType* tree, int& i, int& j);
+  template<typename TreeType>
+  static void GetPointSeeds(const TreeType *tree,int& i, int& j);
 
   /**
    * Get the seeds for splitting a non-leaf node.
    */
-  static void GetBoundSeeds(const TreeType* tree, int& i, int& j);
+  template<typename TreeType>
+  static void GetBoundSeeds(const TreeType *tree,int& i, int& j);
 
   /**
    * Assign points to the two new nodes.
    */
+  template<typename TreeType>
   static void AssignPointDestNode(TreeType* oldTree,
                                   TreeType* treeOne,
                                   TreeType* treeTwo,
@@ -66,6 +61,7 @@ class RTreeSplit
   /**
    * Assign nodes to the two new nodes.
    */
+  template<typename TreeType>
   static void AssignNodeDestNode(TreeType* oldTree,
                                  TreeType* treeOne,
                                  TreeType* treeTwo,
@@ -75,15 +71,8 @@ class RTreeSplit
   /**
    * Insert a node into another node.
    */
+  template<typename TreeType>
   static void InsertNodeIntoTree(TreeType* destTree, TreeType* srcNode);
-
- public:
-  /**
-   * Serialize the split.
-   */
-  template<typename Archive>
-  void Serialize(Archive &, const unsigned int /* version */) { };
-
 };
 
 } // namespace tree

src/mlpack/core/tree/rectangle_tree/r_tree_split_impl.hpp

@@ -21,8 +21,7 @@ namespace tree {
  * new nodes into the tree, spliting the parent if necessary.
  */
 template<typename TreeType>
-void RTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
-                                         std::vector<bool>& relevels)
+void RTreeSplit::SplitLeafNode(TreeType *tree,std::vector<bool>& relevels)
 {
   // If we are splitting the root node, we need will do things differently so
   // that the constructor and other methods don't confuse the end user by giving
@@ -36,7 +35,7 @@ void RTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
     tree->NullifyData();
     // Because this was a leaf node, numChildren must be 0.
     tree->Children()[(tree->NumChildren())++] = copy;
-    copy->Split().SplitLeafNode(copy, relevels);
+    RTreeSplit::SplitLeafNode(copy,relevels);
     return;
   }
 
@@ -47,7 +46,7 @@ void RTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
   // rectangles, only points.  We assume that the tree uses Euclidean Distance.
   int i = 0;
   int j = 0;
-  RTreeSplit<TreeType>::GetPointSeeds(tree,i, j);
+  RTreeSplit::GetPointSeeds(tree,i, j);
 
   TreeType* treeOne = new TreeType(tree->Parent());
   TreeType* treeTwo = new TreeType(tree->Parent());
@@ -67,7 +66,7 @@ void RTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
   // just in case, we use an assert.
   assert(par->NumChildren() <= par->MaxNumChildren() + 1);
   if (par->NumChildren() == par->MaxNumChildren() + 1)
-    par->Split().SplitNonLeafNode(par, relevels);
+    RTreeSplit::SplitNonLeafNode(par,relevels);
 
   assert(treeOne->Parent()->NumChildren() <= treeOne->MaxNumChildren());
   assert(treeOne->Parent()->NumChildren() >= treeOne->MinNumChildren());
@@ -86,8 +85,7 @@ void RTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
  * higher up the tree because they were already updated if necessary.
  */
 template<typename TreeType>
-bool RTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
-                                            std::vector<bool>& relevels)
+bool RTreeSplit::SplitNonLeafNode(TreeType *tree,std::vector<bool>& relevels)
 {
   // If we are splitting the root node, we need will do things differently so
   // that the constructor and other methods don't confuse the end user by giving
@@ -100,13 +98,13 @@ bool RTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
     tree->NumChildren() = 0;
     tree->NullifyData();
     tree->Children()[(tree->NumChildren())++] = copy;
-    copy->Split().SplitNonLeafNode(copy, relevels);
+    RTreeSplit::SplitNonLeafNode(copy,relevels);
     return true;
   }
 
   int i = 0;
   int j = 0;
-  RTreeSplit<TreeType>::GetBoundSeeds(tree,i, j);
+  RTreeSplit::GetBoundSeeds(tree,i, j);
 
   assert(i != j);
 
@@ -133,7 +131,7 @@ bool RTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
   assert(par->NumChildren() <= par->MaxNumChildren() + 1);
 
   if (par->NumChildren() == par->MaxNumChildren() + 1)
-    par->Split().SplitNonLeafNode(par, relevels);
+    RTreeSplit::SplitNonLeafNode(par,relevels);
 
   // We have to update the children of each of these new nodes so that they
   // record the correct parent.
@@ -159,9 +157,7 @@ bool RTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
  * The indices of these points will be stored in iRet and jRet.
  */
 template<typename TreeType>
-void RTreeSplit<TreeType>::GetPointSeeds(const TreeType* tree,
-                                         int& iRet,
-                                         int& jRet)
+void RTreeSplit::GetPointSeeds(const TreeType *tree,int& iRet, int& jRet)
 {
   // Here we want to find the pair of points that it is worst to place in the
   // same node.  Because we are just using points, we will simply choose the two
@@ -172,7 +168,8 @@ void RTreeSplit<TreeType>::GetPointSeeds(const TreeType* tree,
     for (size_t j = i + 1; j < tree->Count(); j++)
     {
       const typename TreeType::ElemType score = arma::prod(arma::abs(
-          tree->LocalDataset().col(i) - tree->LocalDataset().col(j)));
+          tree->Dataset().col(tree->Point(i)) -
+          tree->Dataset().col(tree->Point(j))));
 
       if (score > worstPairScore)
       {
@@ -189,9 +186,7 @@ void RTreeSplit<TreeType>::GetPointSeeds(const TreeType* tree,
  * indices of the bounds will be stored in iRet and jRet.
  */
 template<typename TreeType>
-void RTreeSplit<TreeType>::GetBoundSeeds(const TreeType* tree,
-                                         int& iRet,
-                                         int& jRet)
+void RTreeSplit::GetBoundSeeds(const TreeType *tree,int& iRet, int& jRet)
 {
   // Convenience typedef.
   typedef typename TreeType::ElemType ElemType;
@@ -204,10 +199,10 @@ void RTreeSplit<TreeType>::GetBoundSeeds(const TreeType* tree,
       ElemType score = 1.0;
       for (size_t k = 0; k < tree->Bound().Dim(); k++)
       {
-        const ElemType hiMax = std::max(tree->Children()[i]->Bound()[k].Hi(),
-                                        tree->Children()[j]->Bound()[k].Hi());
-        const ElemType loMin = std::min(tree->Children()[i]->Bound()[k].Lo(),
-                                        tree->Children()[j]->Bound()[k].Lo());
+        const ElemType hiMax = std::max(tree->Child(i).Bound()[k].Hi(),
+                                        tree->Child(j).Bound()[k].Hi());
+        const ElemType loMin = std::min(tree->Child(i).Bound()[k].Lo(),
+                                        tree->Child(j).Bound()[k].Lo());
         score *= (hiMax - loMin);
       }
 
@@ -222,11 +217,11 @@ void RTreeSplit<TreeType>::GetBoundSeeds(const TreeType* tree,
 }
 
 template<typename TreeType>
-void RTreeSplit<TreeType>::AssignPointDestNode(TreeType* oldTree,
-                                               TreeType* treeOne,
-                                               TreeType* treeTwo,
-                                               const int intI,
-                                               const int intJ)
+void RTreeSplit::AssignPointDestNode(TreeType* oldTree,
+                                     TreeType* treeOne,
+                                     TreeType* treeTwo,
+                                     const int intI,
+                                     const int intJ)
 {
   // Convenience typedef.
   typedef typename TreeType::ElemType ElemType;
@@ -240,24 +235,20 @@ void RTreeSplit<TreeType>::AssignPointDestNode(TreeType* oldTree,
   treeOne->Count() = 0;
   treeTwo->Count() = 0;
 
-  treeOne->InsertPoint(oldTree->Points()[intI]);
-  treeTwo->InsertPoint(oldTree->Points()[intJ]);
+  treeOne->InsertPoint(oldTree->Point(intI));
+  treeTwo->InsertPoint(oldTree->Point(intJ));
 
   // If intJ is the last point in the tree, we need to switch the order so that
   // we remove the correct points.
   if (intI > intJ)
   {
-    oldTree->Points()[intI] = oldTree->Points()[--end]; // Decrement end.
-    oldTree->LocalDataset().col(intI) = oldTree->LocalDataset().col(end);
-    oldTree->Points()[intJ] = oldTree->Points()[--end]; // Decrement end.
-    oldTree->LocalDataset().col(intJ) = oldTree->LocalDataset().col(end);
+    oldTree->Point(intI) = oldTree->Point(--end); // Decrement end.
+    oldTree->Point(intJ) = oldTree->Point(--end); // Decrement end.
   }
   else
   {
-    oldTree->Points()[intJ] = oldTree->Points()[--end]; // Decrement end.
-    oldTree->LocalDataset().col(intJ) = oldTree->LocalDataset().col(end);
-    oldTree->Points()[intI] = oldTree->Points()[--end]; // Decrement end.
-    oldTree->LocalDataset().col(intI) = oldTree->LocalDataset().col(end);
+    oldTree->Point(intJ) = oldTree->Point(--end); // Decrement end.
+    oldTree->Point(intI) = oldTree->Point(--end); // Decrement end.
   }
 
   size_t numAssignedOne = 1;
@@ -299,7 +290,7 @@ void RTreeSplit<TreeType>::AssignPointDestNode(TreeType* oldTree,
       ElemType newVolTwo = 1.0;
       for (size_t i = 0; i < oldTree->Bound().Dim(); i++)
       {
-        ElemType c = oldTree->LocalDataset().col(index)[i];
+        ElemType c = oldTree->Dataset().col(oldTree->Point(index))[i];
         newVolOne *= treeOne->Bound()[i].Contains(c) ?
             treeOne->Bound()[i].Width() : (c < treeOne->Bound()[i].Lo() ?
             (treeOne->Bound()[i].Hi() - c) : (c - treeOne->Bound()[i].Lo()));
@@ -333,17 +324,16 @@ void RTreeSplit<TreeType>::AssignPointDestNode(TreeType* oldTree,
     // to the appropriate rectangle.
     if (bestRect == 1)
     {
-      treeOne->InsertPoint(oldTree->Points()[bestIndex]);
+      treeOne->InsertPoint(oldTree->Point(bestIndex));
       numAssignedOne++;
     }
     else
     {
-      treeTwo->InsertPoint(oldTree->Points()[bestIndex]);
+      treeTwo->InsertPoint(oldTree->Point(bestIndex));
       numAssignedTwo++;
     }
 
-    oldTree->Points()[bestIndex] = oldTree->Points()[--end]; // Decrement end.
-    oldTree->LocalDataset().col(bestIndex) = oldTree->LocalDataset().col(end);
+    oldTree->Point(bestIndex) = oldTree->Point(--end); // Decrement end.
   }
 
   // See if we need to satisfy the minimum fill.
@@ -352,22 +342,22 @@ void RTreeSplit<TreeType>::AssignPointDestNode(TreeType* oldTree,
     if (numAssignedOne < numAssignedTwo)
     {
       for (size_t i = 0; i < end; i++)
-        treeOne->InsertPoint(oldTree->Points()[i]);
+        treeOne->InsertPoint(oldTree->Point(i));
     }
     else
     {
       for (size_t i = 0; i < end; i++)
-        treeTwo->InsertPoint(oldTree->Points()[i]);
+        treeTwo->InsertPoint(oldTree->Point(i));
     }
   }
 }
 
 template<typename TreeType>
-void RTreeSplit<TreeType>::AssignNodeDestNode(TreeType* oldTree,
-                                              TreeType* treeOne,
-                                              TreeType* treeTwo,
-                                              const int intI,
-                                              const int intJ)
+void RTreeSplit::AssignNodeDestNode(TreeType* oldTree,
+                                    TreeType* treeOne,
+                                    TreeType* treeTwo,
+                                    const int intI,
+                                    const int intJ)
 {
   // Convenience typedef.
   typedef typename TreeType::ElemType ElemType;
@@ -442,7 +432,7 @@ void RTreeSplit<TreeType>::AssignNodeDestNode(TreeType* oldTree,
         // For each of the new rectangles, find the width in this dimension if
         // we add the rectangle at index to the new rectangle.
         const math::RangeType<ElemType>& range =
-            oldTree->Children()[index]->Bound()[i];
+            oldTree->Child(index).Bound()[i];
         newVolOne *= treeOne->Bound()[i].Contains(range) ?
             treeOne->Bound()[i].Width() : (range.Contains(treeOne->Bound()[i]) ?
             range.Width() : (range.Lo() < treeOne->Bound()[i].Lo() ?
@@ -528,8 +518,7 @@ void RTreeSplit<TreeType>::AssignNodeDestNode(TreeType* oldTree,
  * numberOfChildren.
  */
 template<typename TreeType>
-void RTreeSplit<TreeType>::InsertNodeIntoTree(TreeType* destTree,
-                                              TreeType* srcNode)
+void RTreeSplit::InsertNodeIntoTree(TreeType* destTree, TreeType* srcNode)
 {
   destTree->Bound() |= srcNode->Bound();
   destTree->Children()[destTree->NumChildren()++] = srcNode;

src/mlpack/core/tree/rectangle_tree/rectangle_tree.hpp

@@ -14,6 +14,7 @@
 #include "../statistic.hpp"
 #include "r_tree_split.hpp"
 #include "r_tree_descent_heuristic.hpp"
+#include "no_auxiliary_information.hpp"
 
 namespace mlpack {
 namespace tree /** Trees and tree-building procedures. */ {
@@ -34,13 +35,16 @@ namespace tree /** Trees and tree-building procedures. */ {
  * @tparam SplitType The type of split to use when inserting points.
  * @tparam DescentType The heuristic to use when descending the tree to insert
  *    points.
+ * @tparam AuxiliaryInformationType An auxiliary information contained
+ *    in the node. This information depends on the type of the RectangleTree.
  */
 
 template<typename MetricType = metric::EuclideanDistance,
          typename StatisticType = EmptyStatistic,
          typename MatType = arma::mat,
-         template<typename> class SplitType = RTreeSplit,
-         typename DescentType = RTreeDescentHeuristic>
+         typename SplitType = RTreeSplit,
+         typename DescentType = RTreeDescentHeuristic,
+         template<typename> class AuxiliaryInformationType = NoAuxiliaryInformation>
 class RectangleTree
 {
   // The metric *must* be the euclidean distance.
@@ -89,10 +93,8 @@ class RectangleTree
   bool ownsDataset;
   //! The mapping to the dataset
   std::vector<size_t> points;
-  //! The local dataset
-  MatType* localDataset;
-  //! The class that performs the split of the node.
-  SplitType<RectangleTree> split;
+  //! A tree-specific information
+  AuxiliaryInformationType<RectangleTree> auxiliaryInfo;
 
  public:
   //! A single traverser for rectangle type trees.  See
@@ -188,26 +190,23 @@ class RectangleTree
   void SoftDelete();
 
   /**
-   * Set dataset to null. Used for memory management.  Be careful.
+   * Nullify the auxiliary information. Used for memory management.
+   * Be cafeful.
    */
   void NullifyData();
 
   /**
-   * Inserts a point into the tree. The point will be copied to the data matrix
-   * of the leaf node where it is finally inserted, but we pass by reference
-   * since it may be passed many times before it actually reaches a leaf.
+   * Inserts a point into the tree.
    *
-   * @param point The point (arma::vec&) to be inserted.
+   * @param point The index of a point in the dataset.
    */
   void InsertPoint(const size_t point);
 
   /**
    * Inserts a point into the tree, tracking which levels have been inserted
-   * into.  The point will be copied to the data matrix of the leaf node where
-   * it is finally inserted, but we pass by reference since it may be passed
-   * many times before it actually reaches a leaf.
+   * into.
    *
-   * @param point The point (arma::vec&) to be inserted.
+   * @param point The index of a point in the dataset.
    * @param relevels The levels that have been reinserted to on this top level
    *      insertion.
    */
@@ -229,9 +228,8 @@ class RectangleTree
                   std::vector<bool>& relevels);
 
   /**
-   * Deletes a point in the tree.  The point will be removed from the data
-   * matrix of the leaf node where it is store and the bounding rectangles will
-   * be updated.  However, the point will be kept in the centeral dataset. (The
+   * Deletes a point from the treeand, updates the bounding rectangle.
+   * However, the point will be kept in the centeral dataset. (The
    * user may remove it from there if he wants, but he must not change the
    * indices of the other points.) Returns true if the point is successfully
    * removed and false if it is not.  (ie. the point is not in the tree)
@@ -239,10 +237,9 @@ class RectangleTree
   bool DeletePoint(const size_t point);
 
   /**
-   * Deletes a point in the tree, tracking levels.  The point will be removed
-   * from the data matrix of the leaf node where it is store and the bounding
-   * rectangles will be updated.  However, the point will be kept in the
-   * centeral dataset. (The user may remove it from there if he wants, but he
+   * Deletes a point from the tree, updates the bounding rectangle,
+   * tracking levels. However, the point will be kept in the centeral dataset.
+   * (The user may remove it from there if he wants, but he
    * must not change the indices of the other points.) Returns true if the point
    * is successfully removed and false if it is not.  (ie. the point is not in
    * the tree)
@@ -291,10 +288,12 @@ class RectangleTree
   //! Modify the statistic object for this node.
   StatisticType& Stat() { return stat; }
 
-  //! Return the split object of this node.
-  const SplitType<RectangleTree>& Split() const { return split; }
+  //! Return the auxiliary information object of this node.
+  const AuxiliaryInformationType<RectangleTree> &AuxiliaryInfo() const
+  { return auxiliaryInfo; }
   //! Modify the split object of this node.
-  SplitType<RectangleTree>& Split() { return split; }
+  AuxiliaryInformationType<RectangleTree>& AuxiliaryInfo()
+  { return auxiliaryInfo; }
 
   //! Return whether or not this node is a leaf (true if it has no children).
   bool IsLeaf() const;
@@ -329,16 +328,6 @@ class RectangleTree
   //! Modify the dataset which the tree is built on.  Be careful!
   MatType& Dataset() { return const_cast<MatType&>(*dataset); }
 
-  //! Get the points vector for this node.
-  const std::vector<size_t>& Points() const { return points; }
-  //! Modify the points vector for this node.  Be careful!
-  std::vector<size_t>& Points() { return points; }
-
-  //! Get the local dataset of this node.
-  const MatType& LocalDataset() const { return *localDataset; }
-  //! Modify the local dataset of this node.
-  MatType& LocalDataset() { return *localDataset; }
-
   //! Get the metric which the tree uses.
   MetricType Metric() const { return MetricType(); }
 
@@ -430,7 +419,11 @@ class RectangleTree
    *
    * @param index Index of point for which a dataset index is wanted.
    */
-  size_t Point(const size_t index) const;
+  size_t Point(const size_t index) const { return points[index]; }
+
+  //! Modify the index of a particular point in this node.  Be very careful when
+  //! you do this!  You may make the tree invalid.
+  size_t& Point(const size_t index) { return points[index]; }
 
   //! Return the minimum distance to another node.
   ElemType MinDistance(const RectangleTree* other) const
@@ -502,26 +495,6 @@ class RectangleTree
   static bool HasSelfChildren() { return false; }
 
  private:
-  /**
-   * Private copy constructor, available only to fill (pad) the tree to a
-   * specified level.  TO BE REMOVED
-   */
-  RectangleTree(const size_t begin,
-                const size_t count,
-                bound::HRectBound<MetricType> bound,
-                StatisticType stat,
-                const int maxLeafSize = 20) :
-      begin(begin),
-      count(count),
-      bound(bound),
-      stat(stat),
-      maxLeafSize(maxLeafSize) { }
-
-  RectangleTree* CopyMe()
-  {
-    return new RectangleTree(begin, count, bound, stat, maxLeafSize);
-  }
-
   /**
    * Splits the current node, recursing up the tree.
    *

src/mlpack/core/tree/rectangle_tree/rectangle_tree_impl.hpp

@@ -19,9 +19,11 @@ namespace tree {
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+              AuxiliaryInformationType>::
 RectangleTree(const MatType& data,
               const size_t maxLeafSize,
               const size_t minLeafSize,
@@ -42,13 +44,10 @@ RectangleTree(const MatType& data,
     dataset(new MatType(data)),
     ownsDataset(true),
     points(maxLeafSize + 1), // Add one to make splitting the node simpler.
-    localDataset(new MatType(arma::zeros<MatType>(data.n_rows,
-                                                  maxLeafSize + 1)))
+    auxiliaryInfo(this)
 {
   stat = StatisticType(*this);
 
-  split = SplitType<RectangleTree>(this);
-
   // For now, just insert the points in order.
   RectangleTree* root = this;
 
@@ -59,9 +58,11 @@ RectangleTree(const MatType& data,
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+              AuxiliaryInformationType>::
 RectangleTree(MatType&& data,
               const size_t maxLeafSize,
               const size_t minLeafSize,
@@ -82,13 +83,10 @@ RectangleTree(MatType&& data,
     dataset(new MatType(std::move(data))),
     ownsDataset(true),
     points(maxLeafSize + 1), // Add one to make splitting the node simpler.
-    localDataset(new MatType(arma::zeros<MatType>(dataset->n_rows,
-                                                  maxLeafSize + 1)))
+    auxiliaryInfo(this)
 {
   stat = StatisticType(*this);
 
-  split = SplitType<RectangleTree>(this);
-
   // For now, just insert the points in order.
   RectangleTree* root = this;
 
@@ -99,14 +97,17 @@ RectangleTree(MatType&& data,
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+              AuxiliaryInformationType>::
 RectangleTree(
-    RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>*
-        parentNode, const size_t numMaxChildren) :
-    maxNumChildren(numMaxChildren > 0 ? numMaxChildren :
-        parentNode->MaxNumChildren()),
+    RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                  AuxiliaryInformationType>*
+        parentNode,const size_t numMaxChildren) :
+    maxNumChildren(numMaxChildren > 0 ? numMaxChildren : 
+                                        parentNode->MaxNumChildren()),
     minNumChildren(parentNode->MinNumChildren()),
     numChildren(0),
     children(maxNumChildren + 1),
@@ -120,11 +121,9 @@ RectangleTree(
     dataset(&parentNode->Dataset()),
     ownsDataset(false),
     points(maxLeafSize + 1), // Add one to make splitting the node simpler.
-    localDataset(new MatType(arma::zeros<MatType>(parentNode->Bound().Dim(),
-                                                  maxLeafSize + 1)))
+    auxiliaryInfo(this)
 {
   stat = StatisticType(*this);
-  split = SplitType<RectangleTree>(this);
 }
 
 /**
@@ -134,9 +133,11 @@ RectangleTree(
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+              AuxiliaryInformationType>::
 RectangleTree(
     const RectangleTree& other,
     const bool deepCopy) :
@@ -153,30 +154,19 @@ RectangleTree(
     parentDistance(other.ParentDistance()),
     dataset(deepCopy ? new MatType(*other.dataset) : &other.Dataset()),
     ownsDataset(deepCopy),
-    points(other.Points()),
-    localDataset(NULL)
+    points(other.points),
+    auxiliaryInfo(other.auxiliaryInfo)
 {
-  split = SplitType<RectangleTree>(other);
   if (deepCopy)
   {
     if (numChildren > 0)
     {
       for (size_t i = 0; i < numChildren; i++)
-      {
         children[i] = new RectangleTree(*(other.Children()[i]));
-      }
-    }
-    else
-    {
-      localDataset = new MatType(other.LocalDataset());
     }
   }
   else
-  {
     children = other.Children();
-    arma::mat& otherData = const_cast<arma::mat&>(other.LocalDataset());
-    localDataset = &otherData;
-  }
 }
 
 /**
@@ -185,10 +175,12 @@ RectangleTree(
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 template<typename Archive>
-RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+              AuxiliaryInformationType>::
 RectangleTree(
     Archive& ar,
     const typename boost::enable_if<typename Archive::is_loading>::type*) :
@@ -206,9 +198,11 @@ RectangleTree(
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+              AuxiliaryInformationType>::
 ~RectangleTree()
 {
   for (size_t i = 0; i < numChildren; i++)
@@ -217,7 +211,6 @@ RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
   if (ownsDataset)
     delete dataset;
 
-  delete localDataset;
 }
 
 /**
@@ -227,9 +220,11 @@ RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
     SoftDelete()
 {
   parent = NULL;
@@ -242,17 +237,19 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
 }
 
 /**
- * Set the local dataset to null.
+ * Nullify the auxiliary information.
  */
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
     NullifyData()
 {
-  localDataset = NULL;
+  auxiliaryInfo.NullifyData();
 }
 
 /**
@@ -262,9 +259,11 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
     InsertPoint(const size_t point)
 {
   // Expand the bound regardless of whether it is a leaf node.
@@ -277,31 +276,31 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
   // If this is a leaf node, we stop here and add the point.
   if (numChildren == 0)
   {
-    localDataset->col(count) = dataset->col(point);
-    points[count++] = point;
+    if (!auxiliaryInfo.HandlePointInsertion(this, point))
+      points[count++] = point;
+
     SplitNode(lvls);
     return;
   }
 
   // If it is not a leaf node, we use the DescentHeuristic to choose a child
   // to which we recurse.
-  const size_t descentNode = DescentType::ChooseDescentNode(this,
-      dataset->col(point));
+  auxiliaryInfo.HandlePointInsertion(this, point);
+  const size_t descentNode = DescentType::ChooseDescentNode(this, point);
   children[descentNode]->InsertPoint(point, lvls);
 }
 
 /**
  * Inserts a point into the tree, tracking which levels have been inserted into.
- * The point will be copied to the data matrix of the leaf node where it is
- * finally inserted, but we pass by reference since it may be passed many times
- * before it actually reaches a leaf.
  */
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
     InsertPoint(const size_t point, std::vector<bool>& relevels)
 {
   // Expand the bound regardless of whether it is a leaf node.
@@ -310,16 +309,17 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
   // If this is a leaf node, we stop here and add the point.
   if (numChildren == 0)
   {
-    localDataset->col(count) = dataset->col(point);
-    points[count++] = point;
+    if (!auxiliaryInfo.HandlePointInsertion(this, point))
+      points[count++] = point;
+
     SplitNode(relevels);
     return;
   }
 
   // If it is not a leaf node, we use the DescentHeuristic to choose a child
   // to which we recurse.
-  const size_t descentNode = DescentType::ChooseDescentNode(this,
-      dataset->col(point));
+  auxiliaryInfo.HandlePointInsertion(this, point);
+  const size_t descentNode = DescentType::ChooseDescentNode(this,point);
   children[descentNode]->InsertPoint(point, relevels);
 }
 
@@ -334,9 +334,11 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
     InsertNode(RectangleTree* node,
                const size_t level,
                std::vector<bool>& relevels)
@@ -345,12 +347,16 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
   bound |= node->Bound();
   if (level == TreeDepth())
   {
-    children[numChildren++] = node;
-    node->Parent() = this;
+    if (!auxiliaryInfo.HandleNodeInsertion(this, node, true))
+    {
+      children[numChildren++] = node;
+      node->Parent() = this;
+    }
     SplitNode(relevels);
   }
   else
   {
+    auxiliaryInfo.HandleNodeInsertion(this, node, false);
     const size_t descentNode = DescentType::ChooseDescentNode(this, node);
     children[descentNode]->InsertNode(node, level, relevels);
   }
@@ -363,9 +369,11 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
     DeletePoint(const size_t point)
 {
   // It is possible that this will cause a reinsertion, so we need to handle the
@@ -384,8 +392,9 @@ bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
     {
       if (points[i] == point)
       {
-        localDataset->col(i) = localDataset->col(--count); // Decrement count.
-        points[i] = points[count];
+        if (!auxiliaryInfo.HandlePointDeletion(this, i))
+          points[i] = points[--count];
+
         // This function wil ensure that minFill is satisfied.
         CondenseTree(dataset->col(point), lvls, true);
         return true;
@@ -408,9 +417,11 @@ bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
     DeletePoint(const size_t point, std::vector<bool>& relevels)
 {
   if (numChildren == 0)
@@ -419,8 +430,9 @@ bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
     {
       if (points[i] == point)
       {
-        localDataset->col(i) = localDataset->col(--count);
-        points[i] = points[count];
+        if (!auxiliaryInfo.HandlePointDeletion(this, i))
+          points[i] = points[--count];
+
         // This function will ensure that minFill is satisfied.
         CondenseTree(dataset->col(point), relevels, true);
         return true;
@@ -436,6 +448,7 @@ bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
   return false;
 }
 
+
 /**
  * Recurse through the tree to remove the node.  Once we find the node, we
  * shrink the rectangles if necessary.
@@ -443,16 +456,21 @@ bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
     RemoveNode(const RectangleTree* node, std::vector<bool>& relevels)
 {
   for (size_t i = 0; i < numChildren; i++)
   {
     if (children[i] == node)
     {
-      children[i] = children[--numChildren]; // Decrement numChildren.
+      if (!auxiliaryInfo.HandleNodeRemoval(this, i))
+      {
+        children[i] = children[--numChildren]; // Decrement numChildren.
+      }
       CondenseTree(arma::vec(), relevels, false);
       return true;
     }
@@ -472,10 +490,11 @@ bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 size_t RectangleTree<MetricType, StatisticType, MatType, SplitType,
-                     DescentType>::TreeSize() const
+                     DescentType, AuxiliaryInformationType>::TreeSize() const
 {
   int n = 0;
   for (int i = 0; i < numChildren; i++)
@@ -487,10 +506,11 @@ size_t RectangleTree<MetricType, StatisticType, MatType, SplitType,
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 size_t RectangleTree<MetricType, StatisticType, MatType, SplitType,
-                     DescentType>::TreeDepth() const
+                     DescentType, AuxiliaryInformationType>::TreeDepth() const
 {
   int n = 1;
   RectangleTree* currentNode = const_cast<RectangleTree*> (this);
@@ -507,10 +527,11 @@ size_t RectangleTree<MetricType, StatisticType, MatType, SplitType,
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 inline bool RectangleTree<MetricType, StatisticType, MatType, SplitType,
-                          DescentType>::IsLeaf() const
+                          DescentType, AuxiliaryInformationType>::IsLeaf() const
 {
   return (numChildren == 0);
 }
@@ -522,13 +543,14 @@ inline bool RectangleTree<MetricType, StatisticType, MatType, SplitType,
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 inline
 typename RectangleTree<MetricType, StatisticType, MatType, SplitType,
-    DescentType>::ElemType
+    DescentType, AuxiliaryInformationType>::ElemType
 RectangleTree<MetricType, StatisticType, MatType, SplitType,
-    DescentType>::FurthestPointDistance() const
+    DescentType, AuxiliaryInformationType>::FurthestPointDistance() const
 {
   if (!IsLeaf())
     return 0.0;
@@ -547,13 +569,14 @@ RectangleTree<MetricType, StatisticType, MatType, SplitType,
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 inline
 typename RectangleTree<MetricType, StatisticType, MatType, SplitType,
-    DescentType>::ElemType
+    DescentType, AuxiliaryInformationType>::ElemType
 RectangleTree<MetricType, StatisticType, MatType, SplitType,
-    DescentType>::FurthestDescendantDistance() const
+    DescentType, AuxiliaryInformationType>::FurthestDescendantDistance() const
 {
   // Return the distance from the centroid to a corner of the bound.
   return 0.5 * bound.Diameter();
@@ -566,10 +589,11 @@ RectangleTree<MetricType, StatisticType, MatType, SplitType,
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 inline size_t RectangleTree<MetricType, StatisticType, MatType, SplitType,
-                            DescentType>::NumPoints() const
+                       DescentType, AuxiliaryInformationType>::NumPoints() const
 {
   if (numChildren != 0) // This is not a leaf node.
     return 0;
@@ -583,10 +607,11 @@ inline size_t RectangleTree<MetricType, StatisticType, MatType, SplitType,
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 inline size_t RectangleTree<MetricType, StatisticType, MatType, SplitType,
-                            DescentType>::NumDescendants() const
+                  DescentType, AuxiliaryInformationType>::NumDescendants() const
 {
   if (numChildren == 0)
   {
@@ -607,10 +632,11 @@ inline size_t RectangleTree<MetricType, StatisticType, MatType, SplitType,
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 inline size_t RectangleTree<MetricType, StatisticType, MatType, SplitType,
-                            DescentType>::Descendant(const size_t index) const
+    DescentType, AuxiliaryInformationType>::Descendant(const size_t index) const
 {
   // I think this may be inefficient...
   if (numChildren == 0)
@@ -633,20 +659,6 @@ inline size_t RectangleTree<MetricType, StatisticType, MatType, SplitType,
   }
 }
 
-/**
- * Return the index of a particular point contained in this node.
- */
-template<typename MetricType,
-         typename StatisticType,
-         typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-inline size_t RectangleTree<MetricType, StatisticType, MatType, SplitType,
-                            DescentType>::Point(const size_t index) const
-{
-  return points[index];
-}
-
 /**
  * Split the tree.  This calls the SplitType code to split a node.  This method
  * should only be called on a leaf node.
@@ -654,9 +666,11 @@ inline size_t RectangleTree<MetricType, StatisticType, MatType, SplitType,
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
     SplitNode(std::vector<bool>& relevels)
 {
   if (numChildren == 0)
@@ -667,7 +681,7 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
 
     // If we are full, then we need to split (or at least try).  The SplitType
     // takes care of this and of moving up the tree if necessary.
-    split.SplitLeafNode(this, relevels);
+    SplitType::SplitLeafNode(this,relevels);
   }
   else
   {
@@ -677,7 +691,7 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
 
     // If we are full, then we need to split (or at least try).  The SplitType
     // takes care of this and of moving up the tree if necessary.
-    split.SplitNonLeafNode(this, relevels);
+    SplitType::SplitNonLeafNode(this,relevels);
   }
 }
 
@@ -685,9 +699,11 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+              AuxiliaryInformationType>::
 RectangleTree() :
     maxNumChildren(0), // Try to give sensible defaults, but it shouldn't matter
     minNumChildren(0), // because this tree isn't valid anyway and is only used
@@ -699,8 +715,7 @@ RectangleTree() :
     minLeafSize(0),
     parentDistance(0.0),
     dataset(NULL),
-    ownsDataset(false),
-    localDataset(NULL)
+    ownsDataset(false)
 {
   // Nothing to do.
 }
@@ -712,9 +727,11 @@ RectangleTree() :
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
     CondenseTree(const arma::vec& point,
                  std::vector<bool>& relevels,
                  const bool usePoint)
@@ -729,7 +746,10 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
       if (parent->Children()[i] == this)
       {
         // Decrement numChildren.
-        parent->Children()[i] = parent->Children()[--parent->NumChildren()];
+        if (!auxiliaryInfo.HandleNodeRemoval(parent, i))
+        {
+          parent->Children()[i] = parent->Children()[--parent->NumChildren()];
+        }
 
         // We find the root and shrink bounds at the same time.
         bool stillShrinking = true;
@@ -743,7 +763,18 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
         if (stillShrinking)
           stillShrinking = root->ShrinkBoundForBound(bound);
 
-        // Reinsert the points at the root node.
+        stillShrinking = true;
+        root = parent;
+        while (root->Parent() != NULL)
+        {
+          if (stillShrinking)
+            stillShrinking = root->AuxiliaryInfo().UpdateAuxiliaryInfo(root);
+          root = root->Parent();
+        }
+        if (stillShrinking)
+          stillShrinking = root->AuxiliaryInfo().UpdateAuxiliaryInfo(root);
+
+       // Reinsert the points at the root node.
         for (size_t j = 0; j < count; j++)
           root->InsertPoint(points[j], relevels);
 
@@ -768,7 +799,10 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
         if (parent->Children()[j] == this)
         {
           // Decrement numChildren.
-          parent->Children()[j] = parent->Children()[--parent->NumChildren()];
+          if (!auxiliaryInfo.HandleNodeRemoval(parent,j))
+          {
+            parent->Children()[j] = parent->Children()[--parent->NumChildren()];
+          }
           size_t level = TreeDepth();
 
           // We find the root and shrink bounds at the same time.
@@ -783,6 +817,17 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
           if (stillShrinking)
             stillShrinking = root->ShrinkBoundForBound(bound);
 
+          stillShrinking = true;
+          root = parent;
+          while (root->Parent() != NULL)
+          {
+            if (stillShrinking)
+              stillShrinking = root->AuxiliaryInfo().UpdateAuxiliaryInfo(root);
+            root = root->Parent();
+          }
+          if (stillShrinking)
+            stillShrinking = root->AuxiliaryInfo().UpdateAuxiliaryInfo(root);
+
           // Reinsert the nodes at the root node.
           for (size_t i = 0; i < numChildren; i++)
             root->InsertNode(children[i], level, relevels);
@@ -802,7 +847,7 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
       RectangleTree* child = children[0];
 
       // Required for the X tree.
-      if(child->NumChildren() > maxNumChildren)
+      if (child->NumChildren() > maxNumChildren)
       {
         maxNumChildren = child->MaxNumChildren();
         children.resize(maxNumChildren+1);
@@ -818,10 +863,11 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
       for (size_t i = 0; i < child->Count(); i++)
       {
         // In case the tree has a height of two.
-        points[i] = child->Points()[i];
-        localDataset->col(i) = child->LocalDataset().col(i);
+        points[i] = child->Point(i);
       }
 
+      auxiliaryInfo = child->AuxiliaryInfo();
+
       count = child->Count();
       child->SoftDelete();
       return;
@@ -829,9 +875,13 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
   }
 
   // If we didn't delete it, shrink the bound if we need to.
-  if (usePoint && ShrinkBoundForPoint(point) && parent != NULL)
+  if (usePoint &&
+      (ShrinkBoundForPoint(point) || auxiliaryInfo.UpdateAuxiliaryInfo(this)) &&
+      parent != NULL)
     parent->CondenseTree(point, relevels, usePoint);
-  else if (!usePoint && ShrinkBoundForBound(bound) && parent != NULL)
+  else if (!usePoint &&
+           (ShrinkBoundForBound(bound) || auxiliaryInfo.UpdateAuxiliaryInfo(this)) &&
+           parent != NULL)
     parent->CondenseTree(point, relevels, usePoint);
 }
 
@@ -841,9 +891,11 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
     ShrinkBoundForPoint(const arma::vec& point)
 {
   bool shrunk = false;
@@ -856,8 +908,8 @@ bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
         ElemType min = std::numeric_limits<ElemType>::max();
         for (size_t j = 0; j < count; j++)
         {
-          if (localDataset->col(j)[i] < min)
-            min = localDataset->col(j)[i];
+          if (dataset->col(points[j])[i] < min)
+            min = dataset->col(points[j])[i];
         }
 
         if (bound[i].Lo() < min)
@@ -875,8 +927,8 @@ bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
         ElemType max = std::numeric_limits<ElemType>::lowest();
         for (size_t j = 0; j < count; j++)
         {
-          if (localDataset->col(j)[i] > max)
-            max = localDataset->col(j)[i];
+          if (dataset->col(points[j])[i] > max)
+            max = dataset->col(points[j])[i];
         }
 
         if (bound[i].Hi() > max)
@@ -937,9 +989,11 @@ bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
-bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
+bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
     ShrinkBoundForBound(const bound::HRectBound<MetricType>& /* b */)
 {
   // Using the sum is safe since none of the dimensions can increase.
@@ -971,10 +1025,12 @@ bool RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 template<typename Archive>
-void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
     Serialize(Archive& ar,
               const unsigned int /* version */)
 {
@@ -990,8 +1046,6 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
     if (ownsDataset && dataset)
       delete dataset;
 
-    if (localDataset)
-      delete localDataset;
   }
 
   ar & CreateNVP(maxNumChildren, "maxNumChildren");
@@ -1026,8 +1080,7 @@ void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
     ownsDataset = true;
 
   ar & CreateNVP(points, "points");
-  ar & CreateNVP(localDataset, "localDataset");
-  ar & CreateNVP(split, "split");
+  ar & CreateNVP(auxiliaryInfo, "auxiliaryInfo");
 
   // Because 'children' holds mlpack types (that have Serialize()), we can't use
   // the std::vector serialization.

src/mlpack/core/tree/rectangle_tree/single_tree_traverser.hpp

@@ -19,11 +19,12 @@ namespace tree {
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 template<typename RuleType>
 class RectangleTree<MetricType, StatisticType, MatType, SplitType,
-                    DescentType>::SingleTreeTraverser
+                    DescentType, AuxiliaryInformationType>::SingleTreeTraverser
 {
  public:
   /**

src/mlpack/core/tree/rectangle_tree/single_tree_traverser_impl.hpp

@@ -20,10 +20,12 @@ namespace tree {
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 template<typename RuleType>
-RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+              AuxiliaryInformationType>::
 SingleTreeTraverser<RuleType>::SingleTreeTraverser(RuleType& rule) :
     rule(rule),
     numPrunes(0)
@@ -32,10 +34,12 @@ SingleTreeTraverser<RuleType>::SingleTreeTraverser(RuleType& rule) :
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 template<typename RuleType>
-void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType>::
+void RectangleTree<MetricType, StatisticType, MatType, SplitType, DescentType,
+                   AuxiliaryInformationType>::
 SingleTreeTraverser<RuleType>::Traverse(
     const size_t queryIndex,
     const RectangleTree& referenceNode)
@@ -45,7 +49,7 @@ SingleTreeTraverser<RuleType>::Traverse(
   if (referenceNode.IsLeaf())
   {
     for (size_t i = 0; i < referenceNode.Count(); i++)
-      rule.BaseCase(queryIndex, referenceNode.Points()[i]);
+      rule.BaseCase(queryIndex, referenceNode.Point(i));
 
     return;
   }

src/mlpack/core/tree/rectangle_tree/traits.hpp

@@ -21,10 +21,11 @@ namespace tree {
 template<typename MetricType,
          typename StatisticType,
          typename MatType,
-         template<typename> class SplitType,
-         typename DescentType>
+         typename SplitType,
+         typename DescentType,
+         template<typename> class AuxiliaryInformationType>
 class TreeTraits<RectangleTree<MetricType, StatisticType, MatType, SplitType,
-                               DescentType>>
+                               DescentType, AuxiliaryInformationType>>
 {
  public:
   /**

src/mlpack/core/tree/rectangle_tree/typedef.hpp

@@ -38,7 +38,8 @@ using RTree = RectangleTree<MetricType,
                             StatisticType,
                             MatType,
                             RTreeSplit,
-                            RTreeDescentHeuristic>;
+                            RTreeDescentHeuristic,
+                            NoAuxiliaryInformation>;
 
 /**
  * The R*-tree, a more recent variant of the R tree.  This template typedef
@@ -65,7 +66,8 @@ using RStarTree = RectangleTree<MetricType,
                                 StatisticType,
                                 MatType,
                                 RStarTreeSplit,
-                                RStarTreeDescentHeuristic>;
+                                RStarTreeDescentHeuristic,
+                                NoAuxiliaryInformation>;
 
 /**
  * The X-tree, a variant of the R tree with supernodes.  This template typedef
@@ -90,7 +92,44 @@ using XTree = RectangleTree<MetricType,
                             StatisticType,
                             MatType,
                             XTreeSplit,
-                            RTreeDescentHeuristic>;
+                            RTreeDescentHeuristic,
+                            XTreeAuxiliaryInformation>;
+
+/**
+ * The Hilbert R-tree, a variant of the R tree with an ordering along
+ * the Hilbert curve. This template typedef satisfies the TreeType policy API.
+ *
+ * @code
+ * @inproceedings{kamel1994r,
+ *   author = {Kamel, Ibrahim and Faloutsos, Christos},
+ *   title = {Hilbert R-tree: An Improved R-tree Using Fractals},
+ *   booktitle = {Proceedings of the 20th International Conference on Very Large Data Bases},
+ *   series = {VLDB '94},
+ *   year = {1994},
+ *   isbn = {1-55860-153-8},
+ *   pages = {500--509},
+ *   numpages = {10},
+ *   url = {http://dl.acm.org/citation.cfm?id=645920.673001},
+ *   acmid = {673001},
+ *   publisher = {Morgan Kaufmann Publishers Inc.},
+ *   address = {San Francisco, CA, USA}
+ * }
+ * @endcode
+ *
+ * @see @ref trees, RTree, DiscreteHilbertRTree
+ */
+template<typename TreeType>
+using DiscreteHilbertRTreeAuxiliaryInformation =
+      HilbertRTreeAuxiliaryInformation<TreeType,DiscreteHilbertValue>;
+
+template<typename MetricType, typename StatisticType, typename MatType>
+using HilbertRTree = RectangleTree<MetricType,
+                            StatisticType,
+                            MatType,
+                            HilbertRTreeSplit<2>,
+                            HilbertRTreeDescentHeuristic,
+                            DiscreteHilbertRTreeAuxiliaryInformation>;
+
 
 } // namespace tree
 } // namespace mlpack

src/mlpack/core/tree/rectangle_tree/x_tree_auxiliary_information.hpp

@@ -0,0 +1,186 @@
+/**
+ * @file no_auxiliary_information.hpp
+ * @author Mikhail Lozhnikov
+ *
+ * Definition of the XTreeAuxiliaryInformation class, a class that provides
+ * some x-tree specific information about the nodes.
+ */
+#ifndef MLPACK_CORE_TREE_RECTANGLE_TREE_X_TREE_AUXILIARY_INFORMATION_HPP
+#define MLPACK_CORE_TREE_RECTANGLE_TREE_X_TREE_AUXILIARY_INFORMATION_HPP
+
+namespace mlpack {
+namespace tree {
+
+/**
+ * The XTreeAuxiliaryInformation class provides information specific to X trees
+ * for each node in a RectangleTree.
+ */
+template<typename TreeType>
+class XTreeAuxiliaryInformation
+{
+ public:
+  //! Default constructor
+  XTreeAuxiliaryInformation() :
+    normalNodeMaxNumChildren(0),
+    splitHistory(0)
+  { };
+
+  /**
+   * Construct this whith the specified node.
+   *
+   * @param node The node that stores this auxiliary information.
+   */
+  XTreeAuxiliaryInformation(const TreeType* node) :
+      normalNodeMaxNumChildren(node->Parent() ?
+          node->Parent()->AuxiliaryInfo().NormalNodeMaxNumChildren() :
+          node->MaxNumChildren()),
+      splitHistory(node->Bound().Dim())
+  { };
+
+  /**
+   * Create an auxiliary information object by copying from the other node.
+   *
+   * @param other The node from which the information will be copied.
+   */
+  XTreeAuxiliaryInformation(const TreeType& other) :
+      normalNodeMaxNumChildren(
+          other.AuxiliaryInfo().NormalNodeMaxNumChildren()),
+      splitHistory(other.AuxiliaryInfo().SplitHistory())
+  { };
+
+  /**
+   * Some tree types require to save some properties at the insertion process.
+   * This method allows the auxiliary information the option of manipulating the
+   * tree in order to perform the insertion process. If the auxiliary
+   * information does that, then the method should return true; if the method
+   * returns false the RectangleTree performs its default behavior.
+   *
+   * @param node The node in which the point is being inserted.
+   * @param point The global number of the point being inserted.
+   */
+  bool HandlePointInsertion(TreeType* /* node */, const size_t /* point */)
+  {
+    return false;
+  }
+
+  /**
+   * Some tree types require to save some properties at the insertion process.
+   * This method allows the auxiliary information the option of manipulating the
+   * tree in order to perform the insertion process. If the auxiliary
+   * information does that, then the method should return true; if the method
+   * returns false the RectangleTree performs its default behavior.
+   *
+   * @param node The node in which the nodeToInsert is being inserted.
+   * @param nodeToInsert The node being inserted.
+   * @param insertionLevel The level of the tree at which the nodeToInsert
+   *        should be inserted.
+   */
+  bool HandleNodeInsertion(TreeType* /* node */,
+                           TreeType* /* nodeToInsert */,
+                           bool /* insertionLevel */)
+  {
+    return false;
+  }
+
+  /**
+   * Some tree types require to save some properties at the deletion process.
+   * This method allows the auxiliary information the option of manipulating
+   * the tree in order to perform the deletion process. If the auxiliary
+   * information does that, then the method should return true; if the method
+   * returns false the RectangleTree performs its default behavior.
+   * @param node The node from which the point is being deleted.
+   * @param localIndex The local index of the point being deleted.
+   */
+  bool HandlePointDeletion(TreeType* , const size_t)
+  {
+    return false;
+  }
+
+  /**
+   * Some tree types require to save some properties at the deletion process.
+   * This method allows the auxiliary information the option of manipulating
+   * the tree in order to perform the deletion process. If the auxiliary
+   * information does that, then the method should return true; if the method
+   * returns false the RectangleTree performs its default behavior.
+   * @param node The node from which the node is being deleted.
+   * @param nodeIndex The local index of the node being deleted.
+   */
+  bool HandleNodeRemoval(TreeType* , const size_t)
+  {
+    return false;
+  }
+
+  /**
+   * Some tree types require to propagate the information upward.
+   * This method should return false if this is not the case. If true is
+   * returned, the update will be propogated upward.
+   * @param node The node in which the auxiliary information being update.
+   */
+  bool UpdateAuxiliaryInfo(TreeType* )
+  {
+    return false;
+  }
+
+  /**
+   * Nullify the auxiliary information in order to prevent an invalid free.
+   */
+  void NullifyData()
+  { }
+
+  /**
+   * The X tree requires that the tree records it's "split history".  To make
+   * this easy, we use the following structure.
+   */
+  typedef struct SplitHistoryStruct
+  {
+    int lastDimension;
+    std::vector<bool> history;
+
+    SplitHistoryStruct(int dim) : lastDimension(0), history(dim)
+    {
+      for (int i = 0; i < dim; i++)
+        history[i] = false;
+    }
+
+    template<typename Archive>
+    void Serialize(Archive& ar, const unsigned int /* version */)
+    {
+      ar & data::CreateNVP(lastDimension, "lastDimension");
+      ar & data::CreateNVP(history, "history");
+    }
+  } SplitHistoryStruct;
+
+ private:
+    //! The max number of child nodes a non-leaf normal node can have.
+  size_t normalNodeMaxNumChildren;
+  //! A struct to store the "split history" for X trees.
+  SplitHistoryStruct splitHistory;
+
+ public:
+  //! Return the maximum number of a normal node's children.
+  size_t NormalNodeMaxNumChildren() const { return normalNodeMaxNumChildren; }
+  //! Modify the maximum number of a normal node's children.
+  size_t& NormalNodeMaxNumChildren() { return normalNodeMaxNumChildren; }
+  //! Return the split history of the node assosiated with this object.
+  const SplitHistoryStruct& SplitHistory() const { return splitHistory; }
+  //! Modify the split history of the node assosiated with this object.
+  SplitHistoryStruct& SplitHistory() { return splitHistory; }
+
+  /**
+   * Serialize the information.
+   */
+  template<typename Archive>
+  void Serialize(Archive& ar, const unsigned int /* version */)
+  {
+    using data::CreateNVP;
+
+    ar & CreateNVP(normalNodeMaxNumChildren, "normalNodeMaxNumChildren");
+    ar & CreateNVP(splitHistory, "splitHistory");
+  }
+
+};
+
+} // namespace tree
+} // namespace mlpack
+
+#endif  //  MLPACK_CORE_TREE_RECTANGLE_TREE_X_TREE_AUXILIARY_INFORMATION_HPP

src/mlpack/core/tree/rectangle_tree/x_tree_split.hpp

@@ -28,61 +28,25 @@ const double MAX_OVERLAP = 0.2;
  * nodes overflow, we split them, moving up the tree and splitting nodes
  * as necessary.
  */
-template<typename TreeType>
 class XTreeSplit
 {
  public:
-  //! Default constructor
-  XTreeSplit();
-
-  //! Construct this with the specified node.
-  XTreeSplit(const TreeType* node);
-
-  //! Create a copy of the other.split.
-  XTreeSplit(const TreeType& other);
-
   /**
    * Split a leaf node using the algorithm described in "The R*-tree: An
    * Efficient and Robust Access method for Points and Rectangles."  If
    * necessary, this split will propagate upwards through the tree.
    */
-  void SplitLeafNode(TreeType* tree, std::vector<bool>& relevels);
+  template<typename TreeType>
+  static void SplitLeafNode(TreeType *tree,std::vector<bool>& relevels);
 
   /**
    * Split a non-leaf node using the "default" algorithm.  If this is a root
    * node, the tree increases in depth.
    */
-  bool SplitNonLeafNode(TreeType* tree, std::vector<bool>& relevels);
-
-  /**
-   * The X tree requires that the tree records it's "split history".  To make
-   * this easy, we use the following structure.
-   */
-  typedef struct SplitHistoryStruct
-  {
-    int lastDimension;
-    std::vector<bool> history;
-
-    SplitHistoryStruct(int dim) : lastDimension(0), history(dim)
-    {
-      for (int i = 0; i < dim; i++)
-        history[i] = false;
-    }
-
-    template<typename Archive>
-    void Serialize(Archive& ar, const unsigned int /* version */)
-    {
-      ar & data::CreateNVP(lastDimension, "lastDimension");
-      ar & data::CreateNVP(history, "history");
-    }
-  } SplitHistoryStruct;
+  template<typename TreeType>
+  static bool SplitNonLeafNode(TreeType *tree,std::vector<bool>& relevels);
 
  private:
-  //! The max number of child nodes a non-leaf normal node can have.
-  size_t normalNodeMaxNumChildren;
-  //! A struct to store the "split history" for X trees.
-  SplitHistoryStruct splitHistory;
-
   /**
    * Class to allow for faster sorting.
    */
@@ -107,23 +71,8 @@ class XTreeSplit
   /**
    * Insert a node into another node.
    */
+  template<typename TreeType>
   static void InsertNodeIntoTree(TreeType* destTree, TreeType* srcNode);
-
- public:
-  //! Return the maximum number of a normal node's children.
-  size_t NormalNodeMaxNumChildren() const { return normalNodeMaxNumChildren; }
-  //! Modify the maximum number of a normal node's children.
-  size_t& NormalNodeMaxNumChildren() { return normalNodeMaxNumChildren; }
-  //! Return the split history of the node assosiated with this object.
-  const SplitHistoryStruct& SplitHistory() const { return splitHistory; }
-  //! Modify the split history of the node assosiated with this object.
-  SplitHistoryStruct& SplitHistory() { return splitHistory; }
-
-  /**
-   * Serialize the split.
-   */
-  template<typename Archive>
-  void Serialize(Archive& ar, const unsigned int /* version */);
 };
 
 } // namespace tree

src/mlpack/core/tree/rectangle_tree/x_tree_split_impl.hpp

@@ -14,32 +14,6 @@
 namespace mlpack {
 namespace tree {
 
-template<typename TreeType>
-XTreeSplit<TreeType>::XTreeSplit() :
-    normalNodeMaxNumChildren(0),
-    splitHistory(0)
-{
-
-}
-
-template<typename TreeType>
-XTreeSplit<TreeType>::XTreeSplit(const TreeType*node) :
-    normalNodeMaxNumChildren(node->Parent() ?
-        node->Parent()->Split().NormalNodeMaxNumChildren() :
-        node->MaxNumChildren()),
-    splitHistory(node->Bound().Dim())
-{
-
-}
-
-template<typename TreeType>
-XTreeSplit<TreeType>::XTreeSplit(const TreeType &other) :
-    normalNodeMaxNumChildren(other.Split().NormalNodeMaxNumChildren()),
-    splitHistory(other.Split().SplitHistory())
-{
-
-}
-
 /**
  * We call GetPointSeeds to get the two points which will be the initial points
  * in the new nodes We then call AssignPointDestNode to assign the remaining
@@ -47,8 +21,7 @@ XTreeSplit<TreeType>::XTreeSplit(const TreeType &other) :
  * new nodes into the tree, spliting the parent if necessary.
  */
 template<typename TreeType>
-void XTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
-                                         std::vector<bool>& relevels)
+void XTreeSplit::SplitLeafNode(TreeType *tree,std::vector<bool>& relevels)
 {
   // Convenience typedef.
   typedef typename TreeType::ElemType ElemType;
@@ -66,7 +39,7 @@ void XTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
     // Because this was a leaf node, numChildren must be 0.
     tree->Children()[(tree->NumChildren())++] = copy;
     assert(tree->NumChildren() == 1);
-    copy->Split().SplitLeafNode(copy, relevels);
+    XTreeSplit::SplitLeafNode(copy,relevels);
     return;
   }
 
@@ -84,7 +57,7 @@ void XTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
     size_t p = tree->MaxLeafSize() * 0.3;
     if (p == 0)
     {
-      tree->Split().SplitLeafNode(tree, relevels);
+      XTreeSplit::SplitLeafNode(tree,relevels);
       return;
     }
 
@@ -94,17 +67,19 @@ void XTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
     for (size_t i = 0; i < sorted.size(); i++)
     {
       sorted[i].d = tree->Metric().Evaluate(center,
-          tree->LocalDataset().col(i));
+          tree->Dataset().col(tree->Point(i)));
        sorted[i].n = i;
     }
 
     std::sort(sorted.begin(), sorted.end(), structComp<ElemType>);
-    std::vector<int> pointIndices(p);
+    std::vector<size_t> pointIndices(p);
+
     for (size_t i = 0; i < p; i++)
     {
       // We start from the end of sorted.
-      pointIndices[i] = tree->Points()[sorted[sorted.size() - 1 - i].n];
-      root->DeletePoint(tree->Points()[sorted[sorted.size() - 1 - i].n],
+      pointIndices[i] = tree->Point(sorted[sorted.size() - 1 - i].n);
+
+      root->DeletePoint(tree->Point(sorted[sorted.size() - 1 - i].n),
           relevels);
     }
 
@@ -140,7 +115,7 @@ void XTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
     // Since we only have points in the leaf nodes, we only need to sort once.
     std::vector<sortStruct<ElemType>> sorted(tree->Count());
     for (size_t i = 0; i < sorted.size(); i++) {
-      sorted[i].d = tree->LocalDataset().col(i)[j];
+      sorted[i].d = tree->Dataset().col(tree->Point(i))[j];
       sorted[i].n = i;
     }
 
@@ -173,25 +148,25 @@ void XTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
       std::vector<ElemType> minG2(maxG1.size());
       for (size_t k = 0; k < tree->Bound().Dim(); k++)
       {
-        minG1[k] = maxG1[k] = tree->LocalDataset().col(sorted[0].n)[k];
-        minG2[k] = maxG2[k] = tree->LocalDataset().col(
-            sorted[sorted.size() - 1].n)[k];
+        minG1[k] = maxG1[k] = tree->Dataset().col(tree->Point(sorted[0].n))[k];
+        minG2[k] = maxG2[k] = tree->Dataset().col(
+            tree->Point(sorted[sorted.size() - 1].n))[k];
 
         for (size_t l = 1; l < tree->Count() - 1; l++)
         {
           if (l < cutOff)
           {
-            if (tree->LocalDataset().col(sorted[l].n)[k] < minG1[k])
-              minG1[k] = tree->LocalDataset().col(sorted[l].n)[k];
-            else if (tree->LocalDataset().col(sorted[l].n)[k] > maxG1[k])
-              maxG1[k] = tree->LocalDataset().col(sorted[l].n)[k];
+            if (tree->Dataset().col(tree->Point(sorted[l].n))[k] < minG1[k])
+              minG1[k] = tree->Dataset().col(tree->Point(sorted[l].n))[k];
+            else if (tree->Dataset().col(tree->Point(sorted[l].n))[k] > maxG1[k])
+              maxG1[k] = tree->Dataset().col(tree->Point(sorted[l].n))[k];
           }
           else
           {
-            if (tree->LocalDataset().col(sorted[l].n)[k] < minG2[k])
-              minG2[k] = tree->LocalDataset().col(sorted[l].n)[k];
-            else if (tree->LocalDataset().col(sorted[l].n)[k] > maxG2[k])
-              maxG2[k] = tree->LocalDataset().col(sorted[l].n)[k];
+            if (tree->Dataset().col(tree->Point(sorted[l].n))[k] < minG2[k])
+              minG2[k] = tree->Dataset().col(tree->Point(sorted[l].n))[k];
+            else if (tree->Dataset().col(tree->Point(sorted[l].n))[k] > maxG2[k])
+              maxG2[k] = tree->Dataset().col(tree->Point(sorted[l].n))[k];
           }
         }
       }
@@ -239,14 +214,16 @@ void XTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
   std::vector<sortStruct<ElemType>> sorted(tree->Count());
   for (size_t i = 0; i < sorted.size(); i++)
   {
-    sorted[i].d = tree->LocalDataset().col(i)[bestAxis];
+    sorted[i].d = tree->Dataset().col(tree->Point(i))[bestAxis];
     sorted[i].n = i;
   }
 
   std::sort(sorted.begin(), sorted.end(), structComp<ElemType>);
 
-  TreeType* treeOne = new TreeType(tree->Parent(), NormalNodeMaxNumChildren());
-  TreeType* treeTwo = new TreeType(tree->Parent(), NormalNodeMaxNumChildren());
+  TreeType* treeOne = new TreeType(tree->Parent(),
+                            tree->AuxiliaryInfo().NormalNodeMaxNumChildren());
+  TreeType* treeTwo = new TreeType(tree->Parent(),
+                            tree->AuxiliaryInfo().NormalNodeMaxNumChildren());
 
   // The leaf nodes should never have any overlap introduced by the above method
   // since a split axis is chosen and then points are assigned based on their
@@ -256,9 +233,9 @@ void XTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
     for (size_t i = 0; i < tree->Count(); i++)
     {
       if (i < bestAreaIndexOnBestAxis + tree->MinLeafSize())
-        treeOne->InsertPoint(tree->Points()[sorted[i].n]);
+        treeOne->InsertPoint(tree->Point(sorted[i].n));
       else
-        treeTwo->InsertPoint(tree->Points()[sorted[i].n]);
+        treeTwo->InsertPoint(tree->Point(sorted[i].n));
     }
   }
   else
@@ -266,9 +243,9 @@ void XTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
     for (size_t i = 0; i < tree->Count(); i++)
     {
       if (i < bestOverlapIndexOnBestAxis + tree->MinLeafSize())
-        treeOne->InsertPoint(tree->Points()[sorted[i].n]);
+        treeOne->InsertPoint(tree->Point(sorted[i].n));
       else
-        treeTwo->InsertPoint(tree->Points()[sorted[i].n]);
+        treeTwo->InsertPoint(tree->Point(sorted[i].n));
     }
   }
 
@@ -288,16 +265,16 @@ void XTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
   par->Children()[par->NumChildren()++] = treeTwo;
 
   // We now update the split history of each new node.
-  treeOne->Split().SplitHistory().history[bestAxis] = true;
-  treeOne->Split().SplitHistory().lastDimension = bestAxis;
-  treeTwo->Split().SplitHistory().history[bestAxis] = true;
-  treeTwo->Split().SplitHistory().lastDimension = bestAxis;
+  treeOne->AuxiliaryInfo().SplitHistory().history[bestAxis] = true;
+  treeOne->AuxiliaryInfo().SplitHistory().lastDimension = bestAxis;
+  treeTwo->AuxiliaryInfo().SplitHistory().history[bestAxis] = true;
+  treeTwo->AuxiliaryInfo().SplitHistory().lastDimension = bestAxis;
 
   // We only add one at a time, so we should only need to test for equality just
   // in case, we use an assert.
   assert(par->NumChildren() <= par->MaxNumChildren() + 1);
   if (par->NumChildren() == par->MaxNumChildren() + 1)
-    par->Split().SplitNonLeafNode(par, relevels);
+    XTreeSplit::SplitNonLeafNode(par,relevels);
 
   assert(treeOne->Parent()->NumChildren() <=
       treeOne->Parent()->MaxNumChildren());
@@ -319,8 +296,7 @@ void XTreeSplit<TreeType>::SplitLeafNode(TreeType* tree,
  * higher up the tree because they were already updated if necessary.
  */
 template<typename TreeType>
-bool XTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
-                                            std::vector<bool>& relevels)
+bool XTreeSplit::SplitNonLeafNode(TreeType *tree,std::vector<bool>& relevels)
 {
   // Convenience typedef.
   typedef typename TreeType::ElemType ElemType;
@@ -337,7 +313,7 @@ bool XTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
     tree->NumChildren() = 0;
     tree->NullifyData();
     tree->Children()[(tree->NumChildren())++] = copy;
-    copy->Split().SplitNonLeafNode(copy, relevels);
+    XTreeSplit::SplitNonLeafNode(copy,relevels);
     return true;
   }
 
@@ -352,7 +328,8 @@ bool XTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
   std::vector<bool> axes(tree->Bound().Dim());
   std::vector<int> dimensionsLastUsed(tree->NumChildren());
   for (size_t i = 0; i < tree->NumChildren(); i++)
-    dimensionsLastUsed[i] = tree->Child(i).Split().SplitHistory().lastDimension;
+    dimensionsLastUsed[i] = 
+                    tree->Child(i).AuxiliaryInfo().SplitHistory().lastDimension;
   std::sort(dimensionsLastUsed.begin(), dimensionsLastUsed.end());
 
   size_t lastDim = dimensionsLastUsed[dimensionsLastUsed.size()/2];
@@ -363,7 +340,8 @@ bool XTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
   {
     axes[i] = true;
     for (size_t j = 0; j < tree->NumChildren(); j++)
-      axes[i] = axes[i] & tree->Child(j).Split().SplitHistory().history[i];
+      axes[i] = axes[i] & 
+                tree->Child(j).AuxiliaryInfo().SplitHistory().history[i];
     if (axes[i] == true)
     {
       minOverlapSplitDimension = i;
@@ -376,7 +354,8 @@ bool XTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
     {
       axes[i] = true;
       for (size_t j = 0; j < tree->NumChildren(); j++)
-        axes[i] = axes[i] & tree->Child(j).Split().SplitHistory().history[i];
+        axes[i] = axes[i] & 
+                  tree->Child(j).AuxiliaryInfo().SplitHistory().history[i];
       if (axes[i] == true)
       {
         minOverlapSplitDimension = i;
@@ -771,8 +750,8 @@ bool XTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
           (tree->Parent()->NumChildren() == 1))
       {
         // We make the root a supernode instead.
-        tree->Parent()->MaxNumChildren() = tree->MaxNumChildren() +
-            NormalNodeMaxNumChildren();
+        tree->Parent()->MaxNumChildren() = tree->MaxNumChildren() + 
+                              tree->AuxiliaryInfo().NormalNodeMaxNumChildren();
         tree->Parent()->Children().resize(tree->Parent()->MaxNumChildren() + 1);
         tree->Parent()->NumChildren() = tree->NumChildren();
         for (size_t i = 0; i < tree->NumChildren(); i++)
@@ -789,7 +768,8 @@ bool XTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
       }
 
       // If we don't have to worry about the root, we just enlarge this node.
-      tree->MaxNumChildren() += NormalNodeMaxNumChildren();
+      tree->MaxNumChildren() += 
+                              tree->AuxiliaryInfo().NormalNodeMaxNumChildren();
       tree->Children().resize(tree->MaxNumChildren() + 1);
       for (size_t i = 0; i < tree->NumChildren(); i++)
         tree->Child(i).Parent() = tree;
@@ -802,10 +782,10 @@ bool XTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
   }
 
   // Update the split history of each child.
-  treeOne->Split().SplitHistory().history[bestAxis] = true;
-  treeOne->Split().SplitHistory().lastDimension = bestAxis;
-  treeTwo->Split().SplitHistory().history[bestAxis] = true;
-  treeTwo->Split().SplitHistory().lastDimension = bestAxis;
+  treeOne->AuxiliaryInfo().SplitHistory().history[bestAxis] = true;
+  treeOne->AuxiliaryInfo().SplitHistory().lastDimension = bestAxis;
+  treeTwo->AuxiliaryInfo().SplitHistory().history[bestAxis] = true;
+  treeTwo->AuxiliaryInfo().SplitHistory().lastDimension = bestAxis;
 
   // Remove this node and insert treeOne and treeTwo
   TreeType* par = tree->Parent();
@@ -831,9 +811,7 @@ bool XTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
   assert(par->NumChildren() <= par->MaxNumChildren() + 1);
 
   if (par->NumChildren() == par->MaxNumChildren() + 1)
-  {
-    par->Split().SplitNonLeafNode(par, relevels);
-  }
+    XTreeSplit::SplitNonLeafNode(par,relevels);
 
   // We have to update the children of each of these new nodes so that they
   // record the correct parent.
@@ -861,27 +839,13 @@ bool XTreeSplit<TreeType>::SplitNonLeafNode(TreeType* tree,
  * numberOfChildren.
  */
 template<typename TreeType>
-void XTreeSplit<TreeType>::InsertNodeIntoTree(TreeType* destTree, TreeType* srcNode)
+void XTreeSplit::InsertNodeIntoTree(TreeType* destTree, TreeType* srcNode)
 {
   destTree->Bound() |= srcNode->Bound();
   destTree->Children()[destTree->NumChildren()] = srcNode;
   destTree->NumChildren()++;
 }
 
-/**
- * Serialize the split.
- */
-template<typename TreeType>
-template<typename Archive>
-void XTreeSplit<TreeType>::Serialize(Archive& ar,const unsigned int /* version */)
-{
-  using data::CreateNVP;
-
-  ar & CreateNVP(normalNodeMaxNumChildren, "normalNodeMaxNumChildren");
-  ar & CreateNVP(splitHistory, "splitHistory");
-
-}
-
 } // namespace tree
 } // namespace mlpack
 

src/mlpack/methods/lsh/lsh_search_impl.hpp

@@ -167,7 +167,8 @@ void LSHSearch<SortPolicy>::Train(const arma::mat& referenceSet,
   }
 
   // We will store the second hash vectors in this matrix; the second hash
-  // vector for table i will be held in row i.
+  // vector for table i will be held in row i.  We have to use int and not
+  // size_t, otherwise negative numbers are cast to 0.
   arma::Mat<size_t> secondHashVectors(numTables, referenceSet.n_cols);
 
   for (size_t i = 0; i < numTables; i++)
@@ -190,15 +191,20 @@ void LSHSearch<SortPolicy>::Train(const arma::mat& referenceSet,
     hashMat /= hashWidth;
 
     // Step V: Putting the points in the 'secondHashTable' by hashing the key.
-    // Now we hash every key, point ID to its corresponding bucket.
-    secondHashVectors.row(i) = arma::conv_to<arma::Row<size_t>>::from(
-        secondHashWeights.t() * arma::floor(hashMat));
+    // Now we hash every key, point ID to its corresponding bucket.  We must
+    // also normalize the hashes to the range [0, secondHashSize).
+    arma::rowvec unmodVector = secondHashWeights.t() * arma::floor(hashMat);
+    for (size_t j = 0; j < secondHashVectors.n_cols; ++j)
+    {
+      double shs = (double) secondHashSize; // Convenience cast.
+      if (unmodVector[j] >= 0.0)
+        secondHashVectors[j] = size_t(fmod(unmodVector[j], shs));
+      else
+        secondHashVectors[j] = secondHashSize -
+            size_t(fmod(-unmodVector[j], shs));
+    }
   }
 
-  // Normalize hashes (take modulus with secondHashSize).
-  secondHashVectors.transform([secondHashSize](size_t val)
-      { return val % secondHashSize; });
-
   // Now, using the hash vectors for each table, count the number of rows we
   // have in the second hash table.
   arma::Row<size_t> secondHashBinCounts(secondHashSize, arma::fill::zeros);

src/mlpack/methods/neighbor_search/kfn_main.cpp

@@ -62,7 +62,7 @@ PARAM_INT("k", "Number of furthest neighbors to find.", "k", 0);
 // The user may specify the type of tree to use, and a few pararmeters for tree
 // building.
 PARAM_STRING("tree_type", "Type of tree to use: 'kd', 'cover', 'r', 'r-star', "
-    "'x', 'ball'.", "t", "kd");
+    "'x', 'ball', 'hilbert-r'.", "t", "kd");
 PARAM_INT("leaf_size", "Leaf size for tree building.", "l", 20);
 PARAM_FLAG("random_basis", "Before tree-building, project the data onto a "
     "random orthogonal basis.", "R");
@@ -72,6 +72,12 @@ PARAM_INT("seed", "Random seed (if 0, std::time(NULL) is used).", "s", 0);
 PARAM_FLAG("naive", "If true, O(n^2) naive mode is used for computation.", "N");
 PARAM_FLAG("single_mode", "If true, single-tree search is used (as opposed to "
     "dual-tree search).", "s");
+PARAM_DOUBLE("epsilon", "If specified, will do approximate furthest neighbor "
+    "search with given relative error. Must be in the range [0,1).", "e", 0);
+PARAM_DOUBLE("percentage", "If specified, will do approximate furthest neighbor"
+    " search. Must be in the range (0,1] (decimal form). Resultant neighbors "
+    "will be at least (p*100) % of the distance as the true furthest neighbor.",
+    "p", 1);
 
 // Convenience typedef.
 typedef NSModel<FurthestNeighborSort> KFNModel;
@@ -138,6 +144,24 @@ int main(int argc, char *argv[])
     Log::Fatal << "Invalid leaf size: " << lsInt << ".  Must be greater than 0."
         << endl;
 
+  // Sanity check on epsilon.
+  double epsilon = CLI::GetParam<double>("epsilon");
+  if (epsilon < 0 || epsilon >= 1)
+    Log::Fatal << "Invalid epsilon: " << epsilon << ".  Must be in the range "
+        << "[0,1)." << endl;
+
+  // Sanity check on percentage.
+  const double percentage = CLI::GetParam<double>("percentage");
+  if (percentage <= 0 || percentage > 1)
+    Log::Fatal << "Invalid percentage: " << percentage << ".  Must be in the "
+        << "range (0,1] (decimal form)." << endl;
+
+  if (CLI::HasParam("percentage") && CLI::HasParam("epsilon"))
+    Log::Fatal << "Cannot provide both epsilon and percentage." << endl;
+
+  if (CLI::HasParam("percentage"))
+    epsilon = 1 - percentage;
+
   // We either have to load the reference data, or we have to load the model.
   NSModel<FurthestNeighborSort> kfn;
   const bool naive = CLI::HasParam("naive");
@@ -162,9 +186,11 @@ int main(int argc, char *argv[])
       tree = KFNModel::BALL_TREE;
     else if (treeType == "x")
       tree = KFNModel::X_TREE;
+    else if (treeType == "hilbert-r")
+      tree = KFNModel::HILBERT_R_TREE;
     else
       Log::Fatal << "Unknown tree type '" << treeType << "'; valid choices are "
-          << "'kd', 'cover', 'r', 'r-star', 'x' and 'ball'." << endl;
+          << "'kd', 'cover', 'r', 'r-star', 'x', 'ball' and 'hilbert-r'." << endl;
 
     kfn.TreeType() = tree;
     kfn.RandomBasis() = randomBasis;
@@ -175,7 +201,8 @@ int main(int argc, char *argv[])
     Log::Info << "Loaded reference data from '" << referenceFile << "' ("
         << referenceSet.n_rows << "x" << referenceSet.n_cols << ")." << endl;
 
-    kfn.BuildModel(std::move(referenceSet), size_t(lsInt), naive, singleMode);
+    kfn.BuildModel(std::move(referenceSet), size_t(lsInt), naive, singleMode,
+        epsilon);
   }
   else
   {
@@ -191,6 +218,7 @@ int main(int argc, char *argv[])
     kfn.SingleMode() = CLI::HasParam("single_mode");
     kfn.Naive() = CLI::HasParam("naive");
     kfn.LeafSize() = size_t(lsInt);
+    kfn.Epsilon() = epsilon;
   }
 
   // Perform search, if desired.

src/mlpack/methods/neighbor_search/knn_main.cpp

@@ -63,7 +63,7 @@ PARAM_INT("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("tree_type", "Type of tree to use: 'kd', 'cover', 'r', 'r-star', "
-    "'x', 'ball'.", "t", "kd");
+    "'x', 'ball', 'hilbert-r'.", "t", "kd");
 PARAM_INT("leaf_size", "Leaf size for tree building (used for kd-trees, R "
     "trees, and R* trees).", "l", 20);
 PARAM_FLAG("random_basis", "Before tree-building, project the data onto a "
@@ -74,6 +74,8 @@ PARAM_INT("seed", "Random seed (if 0, std::time(NULL) is used).", "s", 0);
 PARAM_FLAG("naive", "If true, O(n^2) naive mode is used for computation.", "N");
 PARAM_FLAG("single_mode", "If true, single-tree search is used (as opposed to "
     "dual-tree search).", "S");
+PARAM_DOUBLE("epsilon", "If specified, will do approximate nearest neighbor "
+    "search with given relative error.", "e", 0);
 
 // Convenience typedef.
 typedef NSModel<NearestNeighborSort> KNNModel;
@@ -137,10 +139,14 @@ int main(int argc, char *argv[])
   // Sanity check on leaf size.
   const int lsInt = CLI::GetParam<int>("leaf_size");
   if (lsInt < 1)
-  {
     Log::Fatal << "Invalid leaf size: " << lsInt << ".  Must be greater "
         "than 0." << endl;
-  }
+
+  // Sanity check on epsilon.
+  const double epsilon = CLI::GetParam<double>("epsilon");
+  if (epsilon < 0)
+    Log::Fatal << "Invalid epsilon: " << epsilon << ".  Must be non-negative. "
+        << endl;
 
   // We either have to load the reference data, or we have to load the model.
   NSModel<NearestNeighborSort> knn;
@@ -166,9 +172,11 @@ int main(int argc, char *argv[])
       tree = KNNModel::BALL_TREE;
     else if (treeType == "x")
       tree = KNNModel::X_TREE;
+    else if (treeType == "hilbert-r")
+      tree = KNNModel::HILBERT_R_TREE;
     else
       Log::Fatal << "Unknown tree type '" << treeType << "'; valid choices are "
-          << "'kd', 'cover', 'r', 'r-star', 'x' and 'ball'." << endl;
+          << "'kd', 'cover', 'r', 'r-star', 'x', 'ball' and 'hilbert-r'." << endl;
 
     knn.TreeType() = tree;
     knn.RandomBasis() = randomBasis;
@@ -180,7 +188,8 @@ int main(int argc, char *argv[])
         << referenceSet.n_rows << " x " << referenceSet.n_cols << ")."
         << endl;
 
-    knn.BuildModel(std::move(referenceSet), size_t(lsInt), naive, singleMode);
+    knn.BuildModel(std::move(referenceSet), size_t(lsInt), naive, singleMode,
+        epsilon);
   }
   else
   {
@@ -196,6 +205,7 @@ int main(int argc, char *argv[])
     knn.SingleMode() = CLI::HasParam("single_mode");
     knn.Naive() = CLI::HasParam("naive");
     knn.LeafSize() = size_t(lsInt);
+    knn.Epsilon() = epsilon;
   }
 
   // Perform search, if desired.

src/mlpack/methods/neighbor_search/neighbor_search.hpp

@@ -84,11 +84,13 @@ class NeighborSearch
    *      dual-tree search).  This overrides singleMode (if it is set to true).
    * @param singleMode If true, single-tree search will be used (as opposed to
    *      dual-tree search).
+   * @param epsilon Relative approximate error (non-negative).
    * @param metric An optional instance of the MetricType class.
    */
   NeighborSearch(const MatType& referenceSet,
                  const bool naive = false,
                  const bool singleMode = false,
+                 const double epsilon = 0,
                  const MetricType metric = MetricType());
 
   /**
@@ -108,11 +110,13 @@ class NeighborSearch
    *      dual-tree search).  This overrides singleMode (if it is set to true).
    * @param singleMode If true, single-tree search will be used (as opposed to
    *      dual-tree search).
+   * @param epsilon Relative approximate error (non-negative).
    * @param metric An optional instance of the MetricType class.
    */
   NeighborSearch(MatType&& referenceSet,
                  const bool naive = false,
                  const bool singleMode = false,
+                 const double epsilon = 0,
                  const MetricType metric = MetricType());
 
   /**
@@ -138,10 +142,12 @@ class NeighborSearch
    * @param referenceSet Set of reference points corresponding to referenceTree.
    * @param singleMode Whether single-tree computation should be used (as
    *      opposed to dual-tree computation).
+   * @param epsilon Relative approximate error (non-negative).
    * @param metric Instantiated distance metric.
    */
   NeighborSearch(Tree* referenceTree,
                  const bool singleMode = false,
+                 const double epsilon = 0,
                  const MetricType metric = MetricType());
 
   /**
@@ -152,10 +158,12 @@ class NeighborSearch
    * @param naive Whether to use naive search.
    * @param singleMode Whether single-tree computation should be used (as
    *      opposed to dual-tree computation).
+   * @param epsilon Relative approximate error (non-negative).
    * @param metric Instantiated metric.
    */
   NeighborSearch(const bool naive = false,
                  const bool singleMode = false,
+                 const double epsilon = 0,
                  const MetricType metric = MetricType());
 
 
@@ -270,6 +278,11 @@ class NeighborSearch
   //! Modify whether or not search is done in single-tree mode.
   bool& SingleMode() { return singleMode; }
 
+  //! Access the relative error to be considered in approximate search.
+  double Epsilon() const { return epsilon; }
+  //! Modify the relative error to be considered in approximate search.
+  double& Epsilon() { return epsilon; }
+
   //! Access the reference dataset.
   const MatType& ReferenceSet() const { return *referenceSet; }
 
@@ -294,6 +307,8 @@ class NeighborSearch
   bool naive;
   //! Indicates if single-tree search is being used (as opposed to dual-tree).
   bool singleMode;
+  //! Indicates the relative error to be considered in approximate search.
+  double epsilon;
 
   //! Instantiation of metric.
   MetricType metric;

src/mlpack/methods/neighbor_search/neighbor_search_impl.hpp

@@ -75,6 +75,7 @@ NeighborSearch<SortPolicy, MetricType, MatType, TreeType, TraversalType>::
 NeighborSearch(const MatType& referenceSetIn,
                const bool naive,
                const bool singleMode,
+               const double epsilon,
                const MetricType metric) :
     referenceTree(naive ? NULL :
         BuildTree<MatType, Tree>(referenceSetIn, oldFromNewReferences)),
@@ -83,12 +84,14 @@ NeighborSearch(const MatType& referenceSetIn,
     setOwner(false),
     naive(naive),
     singleMode(!naive && singleMode), // No single mode if naive.
+    epsilon(epsilon),
     metric(metric),
     baseCases(0),
     scores(0),
     treeNeedsReset(false)
 {
-  // Nothing to do.
+  if (epsilon < 0)
+    throw std::invalid_argument("epsilon must be non-negative");
 }
 
 // Construct the object.
@@ -103,6 +106,7 @@ NeighborSearch<SortPolicy, MetricType, MatType, TreeType, TraversalType>::
 NeighborSearch(MatType&& referenceSetIn,
                const bool naive,
                const bool singleMode,
+               const double epsilon,
                const MetricType metric) :
     referenceTree(naive ? NULL :
         BuildTree<MatType, Tree>(std::move(referenceSetIn),
@@ -113,12 +117,14 @@ NeighborSearch(MatType&& referenceSetIn,
     setOwner(naive),
     naive(naive),
     singleMode(!naive && singleMode),
+    epsilon(epsilon),
     metric(metric),
     baseCases(0),
     scores(0),
     treeNeedsReset(false)
 {
-  // Nothing to do.
+  if (epsilon < 0)
+    throw std::invalid_argument("epsilon must be non-negative");
 }
 
 // Construct the object.
@@ -132,6 +138,7 @@ template<typename SortPolicy,
 NeighborSearch<SortPolicy, MetricType, MatType, TreeType, TraversalType>::
 NeighborSearch(Tree* referenceTree,
                const bool singleMode,
+               const double epsilon,
                const MetricType metric) :
     referenceTree(referenceTree),
     referenceSet(&referenceTree->Dataset()),
@@ -139,12 +146,14 @@ NeighborSearch(Tree* referenceTree,
     setOwner(false),
     naive(false),
     singleMode(singleMode),
+    epsilon(epsilon),
     metric(metric),
     baseCases(0),
     scores(0),
     treeNeedsReset(false)
 {
-  // Nothing else to initialize.
+  if (epsilon < 0)
+    throw std::invalid_argument("epsilon must be non-negative");
 }
 
 // Construct the object without a reference dataset.
@@ -158,6 +167,7 @@ template<typename SortPolicy,
 NeighborSearch<SortPolicy, MetricType, MatType, TreeType, TraversalType>::
     NeighborSearch(const bool naive,
                    const bool singleMode,
+                   const double epsilon,
                    const MetricType metric) :
     referenceTree(NULL),
     referenceSet(new MatType()), // Empty matrix.
@@ -165,11 +175,14 @@ NeighborSearch<SortPolicy, MetricType, MatType, TreeType, TraversalType>::
     setOwner(true),
     naive(naive),
     singleMode(singleMode),
+    epsilon(epsilon),
     metric(metric),
     baseCases(0),
     scores(0),
     treeNeedsReset(false)
 {
+  if (epsilon < 0)
+    throw std::invalid_argument("epsilon must be non-negative");
   // Build the tree on the empty dataset, if necessary.
   if (!naive)
   {
@@ -364,7 +377,8 @@ Search(const MatType& querySet,
   if (naive)
   {
     // Create the helper object for the tree traversal.
-    RuleType rules(*referenceSet, querySet, *neighborPtr, *distancePtr, metric);
+    RuleType rules(*referenceSet, querySet, *neighborPtr, *distancePtr, metric,
+        epsilon);
 
     // The naive brute-force traversal.
     for (size_t i = 0; i < querySet.n_cols; ++i)
@@ -376,7 +390,8 @@ Search(const MatType& querySet,
   else if (singleMode)
   {
     // Create the helper object for the tree traversal.
-    RuleType rules(*referenceSet, querySet, *neighborPtr, *distancePtr, metric);
+    RuleType rules(*referenceSet, querySet, *neighborPtr, *distancePtr, metric,
+        epsilon);
 
     // Create the traverser.
     typename Tree::template SingleTreeTraverser<RuleType> traverser(rules);
@@ -402,7 +417,7 @@ Search(const MatType& querySet,
 
     // Create the helper object for the tree traversal.
     RuleType rules(*referenceSet, queryTree->Dataset(), *neighborPtr,
-        *distancePtr, metric);
+        *distancePtr, metric, epsilon);
 
     // Create the traverser.
     TraversalType<RuleType> traverser(rules);
@@ -527,7 +542,8 @@ Search(Tree* queryTree,
 
   // Create the helper object for the traversal.
   typedef NeighborSearchRules<SortPolicy, MetricType, Tree> RuleType;
-  RuleType rules(*referenceSet, querySet, *neighborPtr, distances, metric);
+  RuleType rules(*referenceSet, querySet, *neighborPtr, distances, metric,
+      epsilon);
 
   // Create the traverser.
   TraversalType<RuleType> traverser(rules);
@@ -598,7 +614,7 @@ Search(const size_t k,
   // Create the helper object for the traversal.
   typedef NeighborSearchRules<SortPolicy, MetricType, Tree> RuleType;
   RuleType rules(*referenceSet, *referenceSet, *neighborPtr, *distancePtr,
-      metric, true /* don't return the same point as nearest neighbor */);
+      metric, epsilon, true /* don't return the same point as nearest neighbor */);
 
   if (naive)
   {

src/mlpack/methods/neighbor_search/neighbor_search_rules.hpp

@@ -22,6 +22,7 @@ class NeighborSearchRules
                       arma::Mat<size_t>& neighbors,
                       arma::mat& distances,
                       MetricType& metric,
+                      const double epsilon = 0,
                       const bool sameSet = false);
   /**
    * Get the distance from the query point to the reference point.
@@ -120,6 +121,9 @@ class NeighborSearchRules
   //! Denotes whether or not the reference and query sets are the same.
   bool sameSet;
 
+  //! Relative error to be considered in approximate search.
+  const double epsilon;
+
   //! The last query point BaseCase() was called with.
   size_t lastQueryIndex;
   //! The last reference point BaseCase() was called with.

src/mlpack/methods/neighbor_search/neighbor_search_rules_impl.hpp

@@ -1,8 +1,8 @@
 /**
- * @file nearest_neighbor_rules_impl.hpp
+ * @file neighbor_search_rules_impl.hpp
  * @author Ryan Curtin
  *
- * Implementation of NearestNeighborRules.
+ * Implementation of NeighborSearchRules.
  */
 #ifndef MLPACK_METHODS_NEIGHBOR_SEARCH_NEAREST_NEIGHBOR_RULES_IMPL_HPP
 #define MLPACK_METHODS_NEIGHBOR_SEARCH_NEAREST_NEIGHBOR_RULES_IMPL_HPP
@@ -20,6 +20,7 @@ NeighborSearchRules<SortPolicy, MetricType, TreeType>::NeighborSearchRules(
     arma::Mat<size_t>& neighbors,
     arma::mat& distances,
     MetricType& metric,
+    const double epsilon,
     const bool sameSet) :
     referenceSet(referenceSet),
     querySet(querySet),
@@ -27,6 +28,7 @@ NeighborSearchRules<SortPolicy, MetricType, TreeType>::NeighborSearchRules(
     distances(distances),
     metric(metric),
     sameSet(sameSet),
+    epsilon(epsilon),
     lastQueryIndex(querySet.n_cols),
     lastReferenceIndex(referenceSet.n_cols),
     baseCases(0),
@@ -112,7 +114,8 @@ inline double NeighborSearchRules<SortPolicy, MetricType, TreeType>::Score(
   }
 
   // Compare against the best k'th distance for this query point so far.
-  const double bestDistance = distances(distances.n_rows - 1, queryIndex);
+  double bestDistance = distances(distances.n_rows - 1, queryIndex);
+  bestDistance = SortPolicy::Relax(bestDistance, epsilon);
 
   return (SortPolicy::IsBetter(distance, bestDistance)) ? distance : DBL_MAX;
 }
@@ -128,7 +131,8 @@ inline double NeighborSearchRules<SortPolicy, MetricType, TreeType>::Rescore(
     return oldScore;
 
   // Just check the score again against the distances.
-  const double bestDistance = distances(distances.n_rows - 1, queryIndex);
+  double bestDistance = distances(distances.n_rows - 1, queryIndex);
+  bestDistance = SortPolicy::Relax(bestDistance, epsilon);
 
   return (SortPolicy::IsBetter(oldScore, bestDistance)) ? oldScore : DBL_MAX;
 }
@@ -419,6 +423,8 @@ inline double NeighborSearchRules<SortPolicy, MetricType, TreeType>::
   queryNode.Stat().SecondBound() = bestDistance;
   queryNode.Stat().AuxBound() = auxDistance;
 
+  worstDistance = SortPolicy::Relax(worstDistance, epsilon);
+
   if (SortPolicy::IsBetter(worstDistance, bestDistance))
     return worstDistance;
   else

src/mlpack/methods/neighbor_search/ns_model.hpp

@@ -69,7 +69,11 @@ class MonoSearchVisitor : public boost::static_visitor<void>
 
   MonoSearchVisitor(const size_t k,
                     arma::Mat<size_t>& neighbors,
-                    arma::mat& distances);
+                    arma::mat& distances) :
+      k(k),
+      neighbors(neighbors),
+      distances(distances)
+  {};
 };
 
 /**
@@ -177,6 +181,16 @@ class NaiveVisitor : public boost::static_visitor<bool&>
   bool& operator()(NSType *ns) const;
 };
 
+/**
+ * EpsilonVisitor exposes the Epsilon method of the given NSType.
+ */
+class EpsilonVisitor : public boost::static_visitor<double&>
+{
+ public:
+  template<typename NSType>
+  double& operator()(NSType *ns) const;
+};
+
 /**
  * ReferenceSetVisitor exposes the referenceSet of the given NSType.
  */
@@ -215,7 +229,8 @@ class NSModel
     R_TREE,
     R_STAR_TREE,
     BALL_TREE,
-    X_TREE
+    X_TREE,
+    HILBERT_R_TREE
   };
 
  private:
@@ -239,7 +254,8 @@ class NSModel
                  NSType<SortPolicy, tree::RTree>*,
                  NSType<SortPolicy, tree::RStarTree>*,
                  NSType<SortPolicy, tree::BallTree>*,
-                 NSType<SortPolicy, tree::XTree>*> nSearch;
+                 NSType<SortPolicy, tree::XTree>*,
+                 NSType<SortPolicy, tree::HilbertRTree>*> nSearch;
 
  public:
   /**
@@ -266,6 +282,10 @@ class NSModel
   bool Naive() const;
   bool& Naive();
 
+  //! Expose Epsilon.
+  double Epsilon() const;
+  double& Epsilon();
+
   //! Expose leafSize.
   size_t LeafSize() const { return leafSize; }
   size_t& LeafSize() { return leafSize; }
@@ -282,7 +302,8 @@ class NSModel
   void BuildModel(arma::mat&& referenceSet,
                   const size_t leafSize,
                   const bool naive,
-                  const bool singleMode);
+                  const bool singleMode,
+                  const double epsilon = 0);
 
   //! Perform neighbor search.  The query set will be reordered.
   void Search(arma::mat&& querySet,

src/mlpack/methods/neighbor_search/ns_model_impl.hpp

@@ -18,15 +18,6 @@
 namespace mlpack {
 namespace neighbor {
 
-//! Save parameters for monochromatic neighbor search.
-MonoSearchVisitor::MonoSearchVisitor(const size_t k,
-                                     arma::Mat<size_t>& neighbors,
-                                     arma::mat& distances) :
-    k(k),
-    neighbors(neighbors),
-    distances(distances)
-{}
-
 //! Monochromatic neighbor search on the given NSType instance.
 template<typename NSType>
 void MonoSearchVisitor::operator()(NSType *ns) const
@@ -185,6 +176,15 @@ bool& NaiveVisitor::operator()(NSType* ns) const
   throw std::runtime_error("no neighbor search model initialized");
 }
 
+//! Expose the Epsilon method of the given NSType.
+template<typename NSType>
+double& EpsilonVisitor::operator()(NSType* ns) const
+{
+  if (ns)
+    return ns->Epsilon();
+  throw std::runtime_error("no neighbor search model initialized");
+}
+
 //! Expose the referenceSet of the given NSType.
 template<typename NSType>
 const arma::mat& ReferenceSetVisitor::operator()(NSType* ns) const
@@ -293,12 +293,25 @@ bool& NSModel<SortPolicy>::Naive()
   return boost::apply_visitor(NaiveVisitor(), nSearch);
 }
 
+template<typename SortPolicy>
+double NSModel<SortPolicy>::Epsilon() const
+{
+  return boost::apply_visitor(EpsilonVisitor(), nSearch);
+}
+
+template<typename SortPolicy>
+double& NSModel<SortPolicy>::Epsilon()
+{
+  return boost::apply_visitor(EpsilonVisitor(), nSearch);
+}
+
 //! Build the reference tree.
 template<typename SortPolicy>
 void NSModel<SortPolicy>::BuildModel(arma::mat&& referenceSet,
                                      const size_t leafSize,
                                      const bool naive,
-                                     const bool singleMode)
+                                     const bool singleMode,
+                                     const double epsilon)
 {
   // Initialize random basis if necessary.
   if (randomBasis)
@@ -348,23 +361,30 @@ void NSModel<SortPolicy>::BuildModel(arma::mat&& referenceSet,
   switch (treeType)
   {
     case KD_TREE:
-      nSearch = new NSType<SortPolicy, tree::KDTree>(naive, singleMode);
+      nSearch = new NSType<SortPolicy, tree::KDTree>(naive, singleMode,
+          epsilon);
       break;
     case COVER_TREE:
       nSearch = new NSType<SortPolicy, tree::StandardCoverTree>(naive,
-          singleMode);
+          singleMode, epsilon);
       break;
     case R_TREE:
-      nSearch = new NSType<SortPolicy, tree::RTree>(naive, singleMode);
+      nSearch = new NSType<SortPolicy, tree::RTree>(naive, singleMode, epsilon);
       break;
     case R_STAR_TREE:
-      nSearch = new NSType<SortPolicy, tree::RStarTree>(naive, singleMode);
+      nSearch = new NSType<SortPolicy, tree::RStarTree>(naive, singleMode,
+          epsilon);
       break;
     case BALL_TREE:
-      nSearch = new NSType<SortPolicy, tree::BallTree>(naive, singleMode);
+      nSearch = new NSType<SortPolicy, tree::BallTree>(naive, singleMode,
+          epsilon);
       break;
     case X_TREE:
-      nSearch = new NSType<SortPolicy, tree::XTree>(naive, singleMode);
+      nSearch = new NSType<SortPolicy, tree::XTree>(naive, singleMode, epsilon);
+      break;
+    case HILBERT_R_TREE:
+      nSearch = new NSType<SortPolicy, tree::HilbertRTree>(naive, singleMode,
+          epsilon);
       break;
   }
 
@@ -389,13 +409,16 @@ void NSModel<SortPolicy>::Search(arma::mat&& querySet,
   if (randomBasis)
     querySet = q * querySet;
 
-  Log::Info << "Searching for " << k << " nearest neighbors with ";
+  Log::Info << "Searching for " << k << " neighbors with ";
   if (!Naive() && !SingleMode())
     Log::Info << "dual-tree " << TreeName() << " search..." << std::endl;
   else if (!Naive())
     Log::Info << "single-tree " << TreeName() << " search..." << std::endl;
   else
     Log::Info << "brute-force (naive) search..." << std::endl;
+  if (Epsilon() != 0 && !Naive())
+    Log::Info << "Maximum of " << Epsilon() * 100 << "% relative error."
+        << std::endl;
 
   BiSearchVisitor<SortPolicy> search(querySet, k, neighbors, distances,
       leafSize);
@@ -408,13 +431,16 @@ void NSModel<SortPolicy>::Search(const size_t k,
                                  arma::Mat<size_t>& neighbors,
                                  arma::mat& distances)
 {
-  Log::Info << "Searching for " << k << " nearest neighbors with ";
+  Log::Info << "Searching for " << k << " neighbors with ";
   if (!Naive() && !SingleMode())
     Log::Info << "dual-tree " << TreeName() << " search..." << std::endl;
   else if (!Naive())
     Log::Info << "single-tree " << TreeName() << " search..." << std::endl;
   else
     Log::Info << "brute-force (naive) search..." << std::endl;
+  if (Epsilon() != 0 && !Naive())
+    Log::Info << "Maximum of " << Epsilon() * 100 << "% relative error."
+        << std::endl;
 
   MonoSearchVisitor search(k, neighbors, distances);
   boost::apply_visitor(search, nSearch);
@@ -438,6 +464,8 @@ std::string NSModel<SortPolicy>::TreeName() const
       return "ball tree";
     case X_TREE:
       return "X tree";
+    case HILBERT_R_TREE:
+      return "Hilbert R tree";
     default:
       return "unknown tree";
   }

src/mlpack/methods/neighbor_search/sort_policies/furthest_neighbor_sort.cpp

@@ -1,5 +1,5 @@
 /***
- * @file nearest_neighbor_sort.cpp
+ * @file furthest_neighbor_sort.cpp
  * @author Ryan Curtin
  *
  * Implementation of the simple FurthestNeighborSort policy class.
@@ -12,7 +12,7 @@ size_t FurthestNeighborSort::SortDistance(const arma::vec& list,
                                           const arma::Col<size_t>& indices,
                                           double newDistance)
 {
-  // The first element in the list is the nearest neighbor.  We only want to
+  // The first element in the list is the furthest neighbor.  We only want to
   // insert if the new distance is greater than the last element in the list.
   if (newDistance < list[list.n_elem - 1])
     return (size_t() - 1); // Do not insert.

src/mlpack/methods/neighbor_search/sort_policies/furthest_neighbor_sort.hpp

@@ -145,6 +145,23 @@ class FurthestNeighborSort
    */
   static inline double CombineWorst(const double a, const double b)
   { return std::max(a - b, 0.0); }
+
+  /**
+   * Return the given value relaxed.
+   *
+   * @param value Value to relax.
+   * @param epsilon Relative error (non-negative).
+   *
+   * @return double Value relaxed.
+   */
+  static inline double Relax(const double value, const double epsilon)
+  {
+    if (value == 0)
+      return 0;
+    if (value == DBL_MAX || epsilon >= 1)
+      return DBL_MAX;
+    return (1 / (1 - epsilon)) * value;
+  }
 };
 
 } // namespace neighbor

src/mlpack/methods/neighbor_search/sort_policies/nearest_neighbor_sort.hpp

@@ -150,6 +150,21 @@ class NearestNeighborSort
       return DBL_MAX;
     return a + b;
   }
+
+  /**
+   * Return the given value relaxed.
+   *
+   * @param value Value to relax.
+   * @param epsilon Relative error (non-negative).
+   *
+   * @return double Value relaxed.
+   */
+  static inline double Relax(const double value, const double epsilon)
+  {
+    if (value == DBL_MAX)
+      return DBL_MAX;
+    return (1 / (1 + epsilon)) * value;
+  }
 };
 
 } // namespace neighbor

src/mlpack/methods/preprocess/preprocess_split_main.cpp

@@ -68,7 +68,7 @@ int main(int argc, char** argv)
   const double testRatio = CLI::GetParam<double>("test_ratio");
 
   // Check on label parameters.
-  if (CLI::HasParam("input_labels"))
+  if (CLI::HasParam("input_labels_file"))
   {
     if (!CLI::HasParam("training_labels_file"))
     {

src/mlpack/methods/range_search/range_search_main.cpp

@@ -70,7 +70,7 @@ PARAM_DOUBLE("min", "Lower bound in range.", "L", 0.0);
 // The user may specify the type of tree to use, and a few parameters for tree
 // building.
 PARAM_STRING("tree_type", "Type of tree to use: 'kd', 'cover', 'r', 'r-star', "
-    "'x', 'ball'.", "t", "kd");
+    "'x', 'ball', 'hilbert-r'.", "t", "kd");
 PARAM_INT("leaf_size", "Leaf size for tree building.", "l", 20);
 PARAM_FLAG("random_basis", "Before tree-building, project the data onto a "
     "random orthogonal basis.", "R");
@@ -173,9 +173,11 @@ int main(int argc, char *argv[])
       tree = RSModel::BALL_TREE;
     else if (treeType == "x")
       tree = RSModel::X_TREE;
+    else if (treeType == "hilbert-r")
+      tree = RSModel::HILBERT_R_TREE;
     else
       Log::Fatal << "Unknown tree type '" << treeType << "; valid choices are "
-          << "'kd', 'cover', 'r', 'r-star', 'x' and 'ball'." << endl;
+          << "'kd', 'cover', 'r', 'r-star', 'x', 'ball' and 'hilbert-r'." << endl;
 
     rs.TreeType() = tree;
     rs.RandomBasis() = randomBasis;

src/mlpack/methods/range_search/rs_model.cpp

@@ -22,7 +22,8 @@ RSModel::RSModel(TreeTypes treeType, bool randomBasis) :
     rTreeRS(NULL),
     rStarTreeRS(NULL),
     ballTreeRS(NULL),
-    xTreeRS(NULL)
+    xTreeRS(NULL),
+    hilbertRTreeRS(NULL)
 {
   // Nothing to do.
 }
@@ -122,6 +123,11 @@ void RSModel::BuildModel(arma::mat&& referenceSet,
       xTreeRS = new RSType<tree::XTree>(move(referenceSet), naive,
           singleMode);
       break;
+
+    case HILBERT_R_TREE:
+      hilbertRTreeRS = new RSType<tree::HilbertRTree>(move(referenceSet), naive,
+          singleMode);
+      break;
   }
 
   if (!naive)
@@ -231,6 +237,10 @@ void RSModel::Search(arma::mat&& querySet,
     case X_TREE:
       xTreeRS->Search(querySet, range, neighbors, distances);
       break;
+
+    case HILBERT_R_TREE:
+      hilbertRTreeRS->Search(querySet, range, neighbors, distances);
+      break;
   }
 }
 
@@ -273,6 +283,10 @@ void RSModel::Search(const math::Range& range,
     case X_TREE:
       xTreeRS->Search(range, neighbors, distances);
       break;
+
+    case HILBERT_R_TREE:
+      hilbertRTreeRS->Search(range, neighbors, distances);
+      break;
   }
 }
 
@@ -293,6 +307,8 @@ std::string RSModel::TreeName() const
       return "ball tree";
     case X_TREE:
       return "X tree";
+    case HILBERT_R_TREE:
+      return "Hilbert R tree";
     default:
       return "unknown tree";
   }
@@ -313,6 +329,8 @@ void RSModel::CleanMemory()
     delete ballTreeRS;
   if (xTreeRS)
     delete xTreeRS;
+  if (hilbertRTreeRS)
+    delete hilbertRTreeRS;
 
   kdTreeRS = NULL;
   coverTreeRS = NULL;
@@ -320,4 +338,5 @@ void RSModel::CleanMemory()
   rStarTreeRS = NULL;
   ballTreeRS = NULL;
   xTreeRS = NULL;
+  hilbertRTreeRS = NULL;
 }

src/mlpack/methods/range_search/rs_model.hpp

@@ -29,7 +29,8 @@ class RSModel
     R_TREE,
     R_STAR_TREE,
     BALL_TREE,
-    X_TREE
+    X_TREE,
+    HILBERT_R_TREE
   };
 
  private:
@@ -60,6 +61,8 @@ class RSModel
   RSType<tree::BallTree>* ballTreeRS;
   //! X tree based range search object (NULL if not in use).
   RSType<tree::XTree>* xTreeRS;
+  //! Hilbert R tree based range search object (NULL if not in use).
+  RSType<tree::HilbertRTree>* hilbertRTreeRS;
 
  public:
   /**

src/mlpack/methods/range_search/rs_model_impl.hpp

@@ -53,6 +53,10 @@ void RSModel::Serialize(Archive& ar, const unsigned int /* version */)
     case X_TREE:
       ar & CreateNVP(xTreeRS, "range_search_model");
       break;
+
+    case HILBERT_R_TREE:
+      ar & CreateNVP(hilbertRTreeRS, "range_search_model");
+      break;
   }
 }
 
@@ -70,6 +74,8 @@ inline const arma::mat& RSModel::Dataset() const
     return ballTreeRS->ReferenceSet();
   else if (xTreeRS)
     return xTreeRS->ReferenceSet();
+  else if (hilbertRTreeRS)
+    return hilbertRTreeRS->ReferenceSet();
 
   throw std::runtime_error("no range search model initialized");
 }
@@ -88,6 +94,8 @@ inline bool RSModel::SingleMode() const
     return ballTreeRS->SingleMode();
   else if (xTreeRS)
     return xTreeRS->SingleMode();
+  else if (hilbertRTreeRS)
+    return hilbertRTreeRS->SingleMode();
 
   throw std::runtime_error("no range search model initialized");
 }
@@ -106,6 +114,8 @@ inline bool& RSModel::SingleMode()
     return ballTreeRS->SingleMode();
   else if (xTreeRS)
     return xTreeRS->SingleMode();
+  else if (hilbertRTreeRS)
+    return hilbertRTreeRS->SingleMode();
 
   throw std::runtime_error("no range search model initialized");
 }
@@ -124,6 +134,8 @@ inline bool RSModel::Naive() const
     return ballTreeRS->Naive();
   else if (xTreeRS)
     return xTreeRS->Naive();
+  else if (hilbertRTreeRS)
+    return hilbertRTreeRS->Naive();
 
   throw std::runtime_error("no range search model initialized");
 }
@@ -142,6 +154,8 @@ inline bool& RSModel::Naive()
     return ballTreeRS->Naive();
   else if (xTreeRS)
     return xTreeRS->Naive();
+  else if (hilbertRTreeRS)
+    return hilbertRTreeRS->Naive();
 
   throw std::runtime_error("no range search model initialized");
 }

src/mlpack/methods/rann/krann_main.cpp

@@ -64,7 +64,7 @@ PARAM_INT("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("tree_type", "Type of tree to use: 'kd', 'cover', 'r', or "
-    "'x', 'r-star'.", "t", "kd");
+    "'x', 'r-star', 'hilbert-r'.", "t", "kd");
 PARAM_INT("leaf_size", "Leaf size for tree building (used for kd-trees, R "
     "trees, and R* trees).", "l", 20);
 PARAM_FLAG("random_basis", "Before tree-building, project the data onto a "
@@ -172,9 +172,11 @@ int main(int argc, char *argv[])
       tree = RANNModel::R_STAR_TREE;
     else if (treeType == "x")
       tree = RANNModel::X_TREE;
+    else if (treeType == "hilbert-r")
+      tree = RANNModel::HILBERT_R_TREE;
     else
       Log::Fatal << "Unknown tree type '" << treeType << "'; valid choices are "
-          << "'kd', 'cover', 'r', 'r-star' and 'x'." << endl;
+          << "'kd', 'cover', 'r', 'r-star', 'x' and 'hilbert-r'." << endl;
 
     rann.TreeType() = tree;
     rann.RandomBasis() = randomBasis;

src/mlpack/methods/rann/ra_model.hpp

@@ -40,7 +40,8 @@ class RAModel
     COVER_TREE,
     R_TREE,
     R_STAR_TREE,
-    X_TREE
+    X_TREE,
+    HILBERT_R_TREE
   };
 
  private:
@@ -73,6 +74,8 @@ class RAModel
   RAType<tree::RStarTree>* rStarTreeRA;
   //! Non-NULL if the X tree is used.
   RAType<tree::XTree>* xTreeRA;
+  //! Non-NULL if the Hilbert R tree is used.
+  RAType<tree::HilbertRTree>* hilbertRTreeRA;
 
  public:
   /**

src/mlpack/methods/rann/ra_model_impl.hpp

@@ -22,7 +22,8 @@ RAModel<SortPolicy>::RAModel(const TreeTypes treeType, const bool randomBasis) :
     coverTreeRA(NULL),
     rTreeRA(NULL),
     rStarTreeRA(NULL),
-    xTreeRA(NULL)
+    xTreeRA(NULL),
+    hilbertRTreeRA(NULL)
 {
   // Nothing to do.
 }
@@ -40,6 +41,8 @@ RAModel<SortPolicy>::~RAModel()
     delete rStarTreeRA;
   if (xTreeRA)
     delete xTreeRA;
+  if (hilbertRTreeRA)
+    delete hilbertRTreeRA;
 }
 
 template<typename SortPolicy>
@@ -64,6 +67,8 @@ void RAModel<SortPolicy>::Serialize(Archive& ar,
       delete rStarTreeRA;
     if (xTreeRA)
       delete xTreeRA;
+    if (hilbertRTreeRA)
+      delete hilbertRTreeRA;
 
     // Set all the pointers to NULL.
     kdTreeRA = NULL;
@@ -71,6 +76,7 @@ void RAModel<SortPolicy>::Serialize(Archive& ar,
     rTreeRA = NULL;
     rStarTreeRA = NULL;
     xTreeRA = NULL;
+    hilbertRTreeRA = NULL;
   }
 
   // We only need to serialize one of the kRANN objects.
@@ -91,6 +97,9 @@ void RAModel<SortPolicy>::Serialize(Archive& ar,
     case X_TREE:
       ar & data::CreateNVP(xTreeRA, "ra_model");
       break;
+    case HILBERT_R_TREE:
+      ar & data::CreateNVP(hilbertRTreeRA, "ra_model");
+      break;
   }
 }
 
@@ -107,6 +116,8 @@ const arma::mat& RAModel<SortPolicy>::Dataset() const
     return rStarTreeRA->ReferenceSet();
   else if (xTreeRA)
     return xTreeRA->ReferenceSet();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->ReferenceSet();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -125,6 +136,8 @@ bool RAModel<SortPolicy>::Naive() const
     return rStarTreeRA->Naive();
   else if (xTreeRA)
     return xTreeRA->Naive();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->Naive();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -143,6 +156,8 @@ bool& RAModel<SortPolicy>::Naive()
     return rStarTreeRA->Naive();
   else if (xTreeRA)
     return xTreeRA->Naive();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->Naive();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -161,6 +176,8 @@ bool RAModel<SortPolicy>::SingleMode() const
     return rStarTreeRA->SingleMode();
   else if (xTreeRA)
     return xTreeRA->SingleMode();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->SingleMode();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -179,6 +196,8 @@ bool& RAModel<SortPolicy>::SingleMode()
     return rStarTreeRA->SingleMode();
   else if (xTreeRA)
     return xTreeRA->SingleMode();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->SingleMode();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -197,6 +216,8 @@ double RAModel<SortPolicy>::Tau() const
     return rStarTreeRA->Tau();
   else if (xTreeRA)
     return xTreeRA->Tau();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->Tau();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -215,6 +236,8 @@ double& RAModel<SortPolicy>::Tau()
     return rStarTreeRA->Tau();
   else if (xTreeRA)
     return xTreeRA->Tau();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->Tau();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -233,6 +256,8 @@ double RAModel<SortPolicy>::Alpha() const
     return rStarTreeRA->Alpha();
   else if (xTreeRA)
     return xTreeRA->Alpha();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->Alpha();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -251,6 +276,8 @@ double& RAModel<SortPolicy>::Alpha()
     return rStarTreeRA->Alpha();
   else if (xTreeRA)
     return xTreeRA->Alpha();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->Alpha();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -269,6 +296,8 @@ bool RAModel<SortPolicy>::SampleAtLeaves() const
     return rStarTreeRA->SampleAtLeaves();
   else if (xTreeRA)
     return xTreeRA->SampleAtLeaves();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->SampleAtLeaves();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -287,6 +316,8 @@ bool& RAModel<SortPolicy>::SampleAtLeaves()
     return rStarTreeRA->SampleAtLeaves();
   else if (xTreeRA)
     return xTreeRA->SampleAtLeaves();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->SampleAtLeaves();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -305,6 +336,8 @@ bool RAModel<SortPolicy>::FirstLeafExact() const
     return rStarTreeRA->FirstLeafExact();
   else if (xTreeRA)
     return xTreeRA->FirstLeafExact();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->FirstLeafExact();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -323,6 +356,8 @@ bool& RAModel<SortPolicy>::FirstLeafExact()
     return rStarTreeRA->FirstLeafExact();
   else if (xTreeRA)
     return xTreeRA->FirstLeafExact();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->FirstLeafExact();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -341,6 +376,8 @@ size_t RAModel<SortPolicy>::SingleSampleLimit() const
     return rStarTreeRA->SingleSampleLimit();
   else if (xTreeRA)
     return xTreeRA->SingleSampleLimit();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->SingleSampleLimit();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -359,6 +396,8 @@ size_t& RAModel<SortPolicy>::SingleSampleLimit()
     return rStarTreeRA->SingleSampleLimit();
   else if (xTreeRA)
     return xTreeRA->SingleSampleLimit();
+  else if (hilbertRTreeRA)
+    return hilbertRTreeRA->SingleSampleLimit();
 
   throw std::runtime_error("no rank-approximate nearest neighbor search model "
       "initialized");
@@ -424,6 +463,8 @@ void RAModel<SortPolicy>::BuildModel(arma::mat&& referenceSet,
     delete rStarTreeRA;
   if (xTreeRA)
     delete xTreeRA;
+  if (hilbertRTreeRA)
+    delete hilbertRTreeRA;
 
   if (randomBasis)
     referenceSet = q * referenceSet;
@@ -472,6 +513,10 @@ void RAModel<SortPolicy>::BuildModel(arma::mat&& referenceSet,
       xTreeRA = new RAType<tree::XTree>(std::move(referenceSet), naive,
           singleMode);
       break;
+    case HILBERT_R_TREE:
+      hilbertRTreeRA = new RAType<tree::HilbertRTree>(std::move(referenceSet),
+          naive, singleMode);
+      break;
   }
 
   if (!naive)
@@ -549,6 +594,10 @@ void RAModel<SortPolicy>::Search(arma::mat&& querySet,
       // No mapping necessary.
       xTreeRA->Search(querySet, k, neighbors, distances);
       break;
+    case HILBERT_R_TREE:
+      // No mapping necessary.
+      hilbertRTreeRA->Search(querySet, k, neighbors, distances);
+      break;
   }
 }
 
@@ -583,6 +632,9 @@ void RAModel<SortPolicy>::Search(const size_t k,
     case X_TREE:
       xTreeRA->Search(k, neighbors, distances);
       break;
+    case HILBERT_R_TREE:
+      hilbertRTreeRA->Search(k, neighbors, distances);
+      break;
   }
 }
 
@@ -601,6 +653,8 @@ std::string RAModel<SortPolicy>::TreeName() const
       return "R* tree";
     case X_TREE:
       return "X tree";
+    case HILBERT_R_TREE:
+      return "Hilbert R tree";
     default:
       return "unknown tree";
   }

src/mlpack/tests/CMakeLists.txt

@@ -28,9 +28,11 @@ add_executable(mlpack_test
   kernel_pca_test.cpp
   kernel_traits_test.cpp
   kfn_test.cpp
+  akfn_test.cpp
   kmeans_test.cpp
   knn_test.cpp
   krann_search_test.cpp
+  aknn_test.cpp
   lars_test.cpp
   lbfgs_test.cpp
   lin_alg_test.cpp

src/mlpack/tests/activation_functions_test.cpp

@@ -25,7 +25,7 @@
 #include <mlpack/methods/ann/layer/hard_tanh_layer.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::ann;

src/mlpack/tests/ada_delta_test.cpp

@@ -12,7 +12,7 @@
 #include <mlpack/methods/logistic_regression/logistic_regression.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace arma;
 using namespace mlpack::optimization;

src/mlpack/tests/adaboost_test.cpp

@@ -8,7 +8,7 @@
 #include <mlpack/methods/adaboost/adaboost.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 #include "serialization.hpp"
 
 using namespace arma;

src/mlpack/tests/adam_test.cpp

@@ -11,7 +11,7 @@
 #include <mlpack/methods/logistic_regression/logistic_regression.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace arma;
 using namespace mlpack::optimization;

src/mlpack/tests/akfn_test.cpp

@@ -0,0 +1,240 @@
+/**
+ * @file akfn_test.cpp
+ *
+ * Tests for KFN (k-furthest-neighbors) with different values of epsilon.
+ */
+#include <mlpack/core.hpp>
+#include <mlpack/methods/neighbor_search/neighbor_search.hpp>
+#include <mlpack/core/tree/cover_tree.hpp>
+#include <boost/test/unit_test.hpp>
+#include "test_tools.hpp"
+
+using namespace mlpack;
+using namespace mlpack::neighbor;
+using namespace mlpack::tree;
+using namespace mlpack::metric;
+using namespace mlpack::bound;
+
+BOOST_AUTO_TEST_SUITE(AKFNTest);
+
+/**
+ * Test the dual-tree furthest-neighbors method with different values for
+ * epsilon. This uses both a query and reference dataset.
+ *
+ * Errors are produced if the results are not according to relative error.
+ */
+BOOST_AUTO_TEST_CASE(AproxVsExact1)
+{
+  arma::mat dataset;
+
+  if (!data::Load("test_data_3_1000.csv", dataset))
+    BOOST_FAIL("Cannot load test dataset test_data_3_1000.csv!");
+
+  KFN exact(dataset);
+  arma::Mat<size_t> neighborsExact;
+  arma::mat distancesExact;
+  exact.Search(dataset, 15, neighborsExact, distancesExact);
+
+  for (size_t c = 0; c < 4; c++)
+  {
+    KFN* akfn;
+    double epsilon;
+
+    switch (c)
+    {
+      case 0: // Use the dual-tree method with e=0.02.
+        epsilon = 0.02;
+        break;
+      case 1: // Use the dual-tree method with e=0.05.
+        epsilon = 0.05;
+        break;
+      case 2: // Use the dual-tree method with e=0.10.
+        epsilon = 0.10;
+        break;
+      case 3: // Use the dual-tree method with e=0.20.
+        epsilon = 0.20;
+        break;
+    }
+
+    // Now perform the actual calculation.
+    akfn = new KFN(dataset, false, false, epsilon);
+    arma::Mat<size_t> neighborsAprox;
+    arma::mat distancesAprox;
+    akfn->Search(dataset, 15, neighborsAprox, distancesAprox);
+
+    for (size_t i = 0; i < neighborsAprox.n_elem; i++)
+      REQUIRE_RELATIVE_ERR(distancesAprox(i), distancesExact(i), epsilon);
+
+    // Clean the memory.
+    delete akfn;
+  }
+}
+
+/**
+ * Test the dual-tree furthest-neighbors method with the exact method.  This
+ * uses only a reference dataset.
+ *
+ * Errors are produced if the results are not according to relative error.
+ */
+BOOST_AUTO_TEST_CASE(AproxVsExact2)
+{
+  arma::mat dataset;
+
+  if (!data::Load("test_data_3_1000.csv", dataset))
+    BOOST_FAIL("Cannot load test dataset test_data_3_1000.csv!");
+
+  KFN exact(dataset);
+  arma::Mat<size_t> neighborsExact;
+  arma::mat distancesExact;
+  exact.Search(15, neighborsExact, distancesExact);
+
+  KFN akfn(dataset, false, false, 0.05);
+  arma::Mat<size_t> neighborsAprox;
+  arma::mat distancesAprox;
+  akfn.Search(15, neighborsAprox, distancesAprox);
+
+  for (size_t i = 0; i < neighborsAprox.n_elem; i++)
+    REQUIRE_RELATIVE_ERR(distancesAprox[i], distancesExact[i], 0.05);
+}
+
+/**
+ * Test the single-tree furthest-neighbors method with the exact method.  This
+ * uses only a reference dataset.
+ *
+ * Errors are produced if the results are not according to relative error.
+ */
+BOOST_AUTO_TEST_CASE(SingleTreeVsExact)
+{
+  arma::mat dataset;
+
+  if (!data::Load("test_data_3_1000.csv", dataset))
+    BOOST_FAIL("Cannot load test dataset test_data_3_1000.csv!");
+
+  KFN exact(dataset);
+  arma::Mat<size_t> neighborsExact;
+  arma::mat distancesExact;
+  exact.Search(15, neighborsExact, distancesExact);
+
+  KFN akfn(dataset, false, true, 0.05);
+  arma::Mat<size_t> neighborsAprox;
+  arma::mat distancesAprox;
+  akfn.Search(15, neighborsAprox, distancesAprox);
+
+  for (size_t i = 0; i < neighborsAprox.n_elem; i++)
+    REQUIRE_RELATIVE_ERR(distancesAprox[i], distancesExact[i], 0.05);
+}
+
+/**
+ * Test the cover tree single-tree furthest-neighbors method against the exact
+ * method.  This uses only a random reference dataset.
+ *
+ * Errors are produced if the results are not according to relative error.
+ */
+BOOST_AUTO_TEST_CASE(SingleCoverTreeTest)
+{
+  arma::mat dataset;
+  dataset.randu(75, 1000); // 75 dimensional, 1000 points.
+
+  KFN exact(dataset);
+  arma::Mat<size_t> neighborsExact;
+  arma::mat distancesExact;
+  exact.Search(dataset, 15, neighborsExact, distancesExact);
+
+  StandardCoverTree<EuclideanDistance, NeighborSearchStat<FurthestNeighborSort>,
+      arma::mat> tree(dataset);
+
+  NeighborSearch<FurthestNeighborSort, LMetric<2>, arma::mat, StandardCoverTree>
+      coverTreeSearch(&tree, true, 0.05);
+
+  arma::Mat<size_t> neighborsCoverTree;
+  arma::mat distancesCoverTree;
+  coverTreeSearch.Search(dataset, 15, neighborsCoverTree, distancesCoverTree);
+
+  for (size_t i = 0; i < neighborsCoverTree.n_elem; ++i)
+    REQUIRE_RELATIVE_ERR(distancesCoverTree[i], distancesExact[i], 0.05);
+}
+
+/**
+ * Test the cover tree dual-tree furthest neighbors method against the exact
+ * method.
+ *
+ * Errors are produced if the results are not according to relative error.
+ */
+BOOST_AUTO_TEST_CASE(DualCoverTreeTest)
+{
+  arma::mat dataset;
+  data::Load("test_data_3_1000.csv", dataset);
+
+  KFN exact(dataset);
+  arma::Mat<size_t> neighborsExact;
+  arma::mat distancesExact;
+  exact.Search(dataset, 15, neighborsExact, distancesExact);
+
+  StandardCoverTree<EuclideanDistance, NeighborSearchStat<FurthestNeighborSort>,
+      arma::mat> referenceTree(dataset);
+
+  NeighborSearch<FurthestNeighborSort, LMetric<2>, arma::mat, StandardCoverTree>
+      coverTreeSearch(&referenceTree, false, 0.05);
+
+  arma::Mat<size_t> neighborsCoverTree;
+  arma::mat distancesCoverTree;
+  coverTreeSearch.Search(dataset, 15, neighborsCoverTree, distancesCoverTree);
+
+  for (size_t i = 0; i < neighborsCoverTree.n_elem; ++i)
+    REQUIRE_RELATIVE_ERR(distancesCoverTree[i], distancesExact[i], 0.05);
+}
+
+/**
+ * Test the ball tree single-tree furthest-neighbors method against the exact
+ * method.  This uses only a random reference dataset.
+ *
+ * Errors are produced if the results are not according to relative error.
+ */
+BOOST_AUTO_TEST_CASE(SingleBallTreeTest)
+{
+  arma::mat dataset;
+  dataset.randu(75, 1000); // 75 dimensional, 1000 points.
+
+  KFN exact(dataset);
+  arma::Mat<size_t> neighborsExact;
+  arma::mat distancesExact;
+  exact.Search(dataset, 15, neighborsExact, distancesExact);
+
+  NeighborSearch<FurthestNeighborSort, EuclideanDistance, arma::mat, BallTree>
+      ballTreeSearch(dataset, false, true, 0.05);
+
+  arma::Mat<size_t> neighborsBallTree;
+  arma::mat distancesBallTree;
+  ballTreeSearch.Search(dataset, 15, neighborsBallTree, distancesBallTree);
+
+  for (size_t i = 0; i < neighborsBallTree.n_elem; ++i)
+    REQUIRE_RELATIVE_ERR(distancesBallTree(i), distancesExact(i), 0.05);
+}
+
+/**
+ * Test the ball tree dual-tree furthest neighbors method against the exact
+ * method.
+ *
+ * Errors are produced if the results are not according to relative error.
+ */
+BOOST_AUTO_TEST_CASE(DualBallTreeTest)
+{
+  arma::mat dataset;
+  data::Load("test_data_3_1000.csv", dataset);
+
+  KFN exact(dataset);
+  arma::Mat<size_t> neighborsExact;
+  arma::mat distancesExact;
+  exact.Search(15, neighborsExact, distancesExact);
+
+  NeighborSearch<FurthestNeighborSort, EuclideanDistance, arma::mat, BallTree>
+      ballTreeSearch(dataset, false, false, 0.05);
+  arma::Mat<size_t> neighborsBallTree;
+  arma::mat distancesBallTree;
+  ballTreeSearch.Search(15, neighborsBallTree, distancesBallTree);
+
+  for (size_t i = 0; i < neighborsBallTree.n_elem; ++i)
+    REQUIRE_RELATIVE_ERR(distancesBallTree(i), distancesExact(i), 0.05);
+}
+
+BOOST_AUTO_TEST_SUITE_END();

src/mlpack/tests/aknn_test.cpp

@@ -0,0 +1,404 @@
+/**
+ * @file aknn_test.cpp
+ *
+ * Test file for KNN class with different values of epsilon.
+ */
+#include <mlpack/core.hpp>
+#include <mlpack/methods/neighbor_search/neighbor_search.hpp>
+#include <mlpack/methods/neighbor_search/unmap.hpp>
+#include <mlpack/methods/neighbor_search/ns_model.hpp>
+#include <mlpack/core/tree/cover_tree.hpp>
+#include <mlpack/core/tree/example_tree.hpp>
+#include <boost/test/unit_test.hpp>
+#include "test_tools.hpp"
+
+using namespace mlpack;
+using namespace mlpack::neighbor;
+using namespace mlpack::tree;
+using namespace mlpack::metric;
+using namespace mlpack::bound;
+
+BOOST_AUTO_TEST_SUITE(AKNNTest);
+
+/**
+ * Test the dual-tree nearest-neighbors method with different values for
+ * epsilon. This uses both a query and reference dataset.
+ *
+ * Errors are produced if the results are not according to relative error.
+ */
+BOOST_AUTO_TEST_CASE(AproxVsExact1)
+{
+  arma::mat dataset;
+
+  if (!data::Load("test_data_3_1000.csv", dataset))
+    BOOST_FAIL("Cannot load test dataset test_data_3_1000.csv!");
+
+  KNN exact(dataset);
+  arma::Mat<size_t> neighborsExact;
+  arma::mat distancesExact;
+  exact.Search(dataset, 15, neighborsExact, distancesExact);
+
+  for (size_t c = 0; c < 4; c++)
+  {
+    KNN* aknn;
+    double epsilon;
+
+    switch (c)
+    {
+      case 0: // Use the dual-tree method with e=0.02.
+        epsilon = 0.02;
+        break;
+      case 1: // Use the dual-tree method with e=0.05.
+        epsilon = 0.05;
+        break;
+      case 2: // Use the dual-tree method with e=0.10.
+        epsilon = 0.10;
+        break;
+      case 3: // Use the dual-tree method with e=0.20.
+        epsilon = 0.20;
+        break;
+    }
+
+    // Now perform the actual calculation.
+    aknn = new KNN(dataset, false, false, epsilon);
+    arma::Mat<size_t> neighborsAprox;
+    arma::mat distancesAprox;
+    aknn->Search(dataset, 15, neighborsAprox, distancesAprox);
+
+    for (size_t i = 0; i < neighborsAprox.n_elem; i++)
+      REQUIRE_RELATIVE_ERR(distancesAprox(i), distancesExact(i), epsilon);
+
+    // Clean the memory.
+    delete aknn;
+  }
+}
+
+/**
+ * Test the dual-tree nearest-neighbors method with the exact method.  This uses
+ * only a reference dataset.
+ *
+ * Errors are produced if the results are not according to relative error.
+ */
+BOOST_AUTO_TEST_CASE(AproxVsExact2)
+{
+  arma::mat dataset;
+
+  if (!data::Load("test_data_3_1000.csv", dataset))
+    BOOST_FAIL("Cannot load test dataset test_data_3_1000.csv!");
+
+  KNN exact(dataset);
+  arma::Mat<size_t> neighborsExact;
+  arma::mat distancesExact;
+  exact.Search(15, neighborsExact, distancesExact);
+
+  KNN aknn(dataset, false, false, 0.05);
+  arma::Mat<size_t> neighborsAprox;
+  arma::mat distancesAprox;
+  aknn.Search(15, neighborsAprox, distancesAprox);
+
+  for (size_t i = 0; i < neighborsAprox.n_elem; i++)
+    REQUIRE_RELATIVE_ERR(distancesAprox(i), distancesExact(i), 0.05);
+}
+
+/**
+ * Test the single-tree nearest-neighbors method with the exact method.  This
+ * uses only a reference dataset.
+ *
+ * Errors are produced if the results are not according to relative error.
+ */
+BOOST_AUTO_TEST_CASE(SingleTreeAproxVsExact)
+{
+  arma::mat dataset;
+
+  if (!data::Load("test_data_3_1000.csv", dataset))
+    BOOST_FAIL("Cannot load test dataset test_data_3_1000.csv!");
+
+  KNN exact(dataset);
+  arma::Mat<size_t> neighborsExact;
+  arma::mat distancesExact;
+  exact.Search(15, neighborsExact, distancesExact);
+
+  KNN aknn(dataset, false, true, 0.05);
+  arma::Mat<size_t> neighborsAprox;
+  arma::mat distancesAprox;
+  aknn.Search(15, neighborsAprox, distancesAprox);
+
+  for (size_t i = 0; i < neighborsAprox.n_elem; i++)
+    REQUIRE_RELATIVE_ERR(distancesAprox[i], distancesExact[i], 0.05);
+}
+
+/**
+ * Test the cover tree single-tree nearest-neighbors method against the exact
+ * method.  This uses only a random reference dataset.
+ *
+ * Errors are produced if the results are not according to relative error.
+ */
+BOOST_AUTO_TEST_CASE(SingleCoverTreeTest)
+{
+  arma::mat dataset;
+  dataset.randu(75, 1000); // 75 dimensional, 1000 points.
+
+  KNN exact(dataset);
+  arma::Mat<size_t> neighborsExact;
+  arma::mat distancesExact;
+  exact.Search(dataset, 15, neighborsExact, distancesExact);
+
+  StandardCoverTree<EuclideanDistance, NeighborSearchStat<NearestNeighborSort>,
+      arma::mat> tree(dataset);
+
+  NeighborSearch<NearestNeighborSort, LMetric<2>, arma::mat, StandardCoverTree>
+      coverTreeSearch(&tree, true, 0.05);
+
+  arma::Mat<size_t> neighborsCoverTree;
+  arma::mat distancesCoverTree;
+  coverTreeSearch.Search(dataset, 15, neighborsCoverTree, distancesCoverTree);
+
+  for (size_t i = 0; i < neighborsCoverTree.n_elem; ++i)
+    REQUIRE_RELATIVE_ERR(distancesCoverTree[i], distancesExact[i], 0.05);
+}
+
+/**
+ * Test the cover tree dual-tree nearest neighbors method against the exact
+ * method.
+ *
+ * Errors are produced if the results are not according to relative error.
+ */
+BOOST_AUTO_TEST_CASE(DualCoverTreeTest)
+{
+  arma::mat dataset;
+  data::Load("test_data_3_1000.csv", dataset);
+
+  KNN exact(dataset);
+  arma::Mat<size_t> neighborsExact;
+  arma::mat distancesExact;
+  exact.Search(dataset, 15, neighborsExact, distancesExact);
+
+  StandardCoverTree<EuclideanDistance, NeighborSearchStat<NearestNeighborSort>,
+      arma::mat> referenceTree(dataset);
+
+  NeighborSearch<NearestNeighborSort, EuclideanDistance, arma::mat,
+      StandardCoverTree> coverTreeSearch(&referenceTree, false, 0.05);
+
+  arma::Mat<size_t> neighborsCoverTree;
+  arma::mat distancesCoverTree;
+  coverTreeSearch.Search(&referenceTree, 15, neighborsCoverTree,
+      distancesCoverTree);
+
+  for (size_t i = 0; i < neighborsCoverTree.n_elem; ++i)
+    REQUIRE_RELATIVE_ERR(distancesCoverTree[i], distancesExact[i], 0.05);
+}
+
+/**
+ * Test the ball tree single-tree nearest-neighbors method against the exact
+ * method.  This uses only a random reference dataset.
+ *
+ * Errors are produced if the results are not according to relative error.
+ */
+BOOST_AUTO_TEST_CASE(SingleBallTreeTest)
+{
+  arma::mat dataset;
+  dataset.randu(50, 300); // 50 dimensional, 300 points.
+
+  KNN exact(dataset);
+  arma::Mat<size_t> neighborsExact;
+  arma::mat distancesExact;
+  exact.Search(dataset, 15, neighborsExact, distancesExact);
+
+  NeighborSearch<NearestNeighborSort, EuclideanDistance, arma::mat, BallTree>
+      ballTreeSearch(dataset, false, true, 0.05);
+
+  arma::Mat<size_t> neighborsBallTree;
+  arma::mat distancesBallTree;
+  ballTreeSearch.Search(dataset, 15, neighborsBallTree, distancesBallTree);
+
+  for (size_t i = 0; i < neighborsBallTree.n_elem; ++i)
+    REQUIRE_RELATIVE_ERR(distancesBallTree(i), distancesExact(i), 0.05);
+}
+
+/**
+ * Test the ball tree dual-tree nearest neighbors method against the exact
+ * method.
+ *
+ * Errors are produced if the results are not according to relative error.
+ */
+BOOST_AUTO_TEST_CASE(DualBallTreeTest)
+{
+  arma::mat dataset;
+  data::Load("test_data_3_1000.csv", dataset);
+
+  KNN exact(dataset);
+  arma::Mat<size_t> neighborsExact;
+  arma::mat distancesExact;
+  exact.Search(15, neighborsExact, distancesExact);
+
+  NeighborSearch<NearestNeighborSort, EuclideanDistance, arma::mat, BallTree>
+      ballTreeSearch(dataset, false, false, 0.05);
+  arma::Mat<size_t> neighborsBallTree;
+  arma::mat distancesBallTree;
+  ballTreeSearch.Search(15, neighborsBallTree, distancesBallTree);
+
+  for (size_t i = 0; i < neighborsBallTree.n_elem; ++i)
+    REQUIRE_RELATIVE_ERR(distancesBallTree(i), distancesExact(i), 0.05);
+}
+
+/**
+ * Make sure sparse nearest neighbors works with kd trees.
+ */
+BOOST_AUTO_TEST_CASE(SparseKNNKDTreeTest)
+{
+  // The dimensionality of these datasets must be high so that the probability
+  // of a completely empty point is very low.  In this case, with dimensionality
+  // 70, the probability of all 70 dimensions being zero is 0.8^70 = 1.65e-7 in
+  // the reference set and 0.9^70 = 6.27e-4 in the query set.
+  arma::sp_mat queryDataset;
+  queryDataset.sprandu(70, 200, 0.2);
+  arma::sp_mat referenceDataset;
+  referenceDataset.sprandu(70, 500, 0.1);
+  arma::mat denseQuery(queryDataset);
+  arma::mat denseReference(referenceDataset);
+
+  typedef NeighborSearch<NearestNeighborSort, EuclideanDistance, arma::sp_mat,
+      KDTree> SparseKNN;
+
+  SparseKNN aknn(referenceDataset, false, false, 0.05);
+  arma::mat distancesSparse;
+  arma::Mat<size_t> neighborsSparse;
+  aknn.Search(queryDataset, 10, neighborsSparse, distancesSparse);
+
+  KNN exact(denseReference);
+  arma::mat distancesExact;
+  arma::Mat<size_t> neighborsExact;
+  exact.Search(denseQuery, 10, neighborsExact, distancesExact);
+
+  for (size_t i = 0; i < neighborsExact.n_cols; ++i)
+    for (size_t j = 0; j < neighborsExact.n_rows; ++j)
+      REQUIRE_RELATIVE_ERR(distancesSparse(j, i), distancesExact(j, i), 0.05);
+}
+
+/**
+ * Ensure that we can build an NSModel<NearestNeighborSearch> and get correct
+ * results.
+ */
+BOOST_AUTO_TEST_CASE(KNNModelTest)
+{
+  typedef NSModel<NearestNeighborSort> KNNModel;
+
+  arma::mat queryData = arma::randu<arma::mat>(10, 50);
+  arma::mat referenceData = arma::randu<arma::mat>(10, 200);
+
+  // Build all the possible models.
+  KNNModel models[14];
+  models[0] = KNNModel(KNNModel::TreeTypes::KD_TREE, true);
+  models[1] = KNNModel(KNNModel::TreeTypes::KD_TREE, false);
+  models[2] = KNNModel(KNNModel::TreeTypes::COVER_TREE, true);
+  models[3] = KNNModel(KNNModel::TreeTypes::COVER_TREE, false);
+  models[4] = KNNModel(KNNModel::TreeTypes::R_TREE, true);
+  models[5] = KNNModel(KNNModel::TreeTypes::R_TREE, false);
+  models[6] = KNNModel(KNNModel::TreeTypes::R_STAR_TREE, true);
+  models[7] = KNNModel(KNNModel::TreeTypes::R_STAR_TREE, false);
+  models[8] = KNNModel(KNNModel::TreeTypes::X_TREE, true);
+  models[9] = KNNModel(KNNModel::TreeTypes::X_TREE, false);
+  models[10] = KNNModel(KNNModel::TreeTypes::BALL_TREE, true);
+  models[11] = KNNModel(KNNModel::TreeTypes::BALL_TREE, false);
+  models[12] = KNNModel(KNNModel::TreeTypes::HILBERT_R_TREE, true);
+  models[13] = KNNModel(KNNModel::TreeTypes::HILBERT_R_TREE, false);
+
+  for (size_t j = 0; j < 3; ++j)
+  {
+    // Get a baseline.
+    KNN aknn(referenceData);
+    arma::Mat<size_t> neighborsExact;
+    arma::mat distancesExact;
+    aknn.Search(queryData, 3, neighborsExact, distancesExact);
+
+    for (size_t i = 0; i < 14; ++i)
+    {
+      // We only have std::move() constructors so make a copy of our data.
+      arma::mat referenceCopy(referenceData);
+      arma::mat queryCopy(queryData);
+      if (j == 0)
+        models[i].BuildModel(std::move(referenceCopy), 20, false, false, 0.05);
+      if (j == 1)
+        models[i].BuildModel(std::move(referenceCopy), 20, false, true, 0.05);
+      if (j == 2)
+        models[i].BuildModel(std::move(referenceCopy), 20, true, false);
+
+      arma::Mat<size_t> neighborsAprox;
+      arma::mat distancesAprox;
+
+      models[i].Search(std::move(queryCopy), 3, neighborsAprox, distancesAprox);
+
+      BOOST_REQUIRE_EQUAL(neighborsAprox.n_rows, neighborsExact.n_rows);
+      BOOST_REQUIRE_EQUAL(neighborsAprox.n_cols, neighborsExact.n_cols);
+      BOOST_REQUIRE_EQUAL(neighborsAprox.n_elem, neighborsExact.n_elem);
+      BOOST_REQUIRE_EQUAL(distancesAprox.n_rows, distancesExact.n_rows);
+      BOOST_REQUIRE_EQUAL(distancesAprox.n_cols, distancesExact.n_cols);
+      BOOST_REQUIRE_EQUAL(distancesAprox.n_elem, distancesExact.n_elem);
+      for (size_t k = 0; k < distancesAprox.n_elem; ++k)
+        REQUIRE_RELATIVE_ERR(distancesAprox[k], distancesExact[k], 0.05);
+    }
+  }
+}
+
+/**
+ * Ensure that we can build an NSModel<NearestNeighborSearch> and get correct
+ * results, in the case where the reference set is the same as the query set.
+ */
+BOOST_AUTO_TEST_CASE(KNNModelMonochromaticTest)
+{
+  typedef NSModel<NearestNeighborSort> KNNModel;
+
+  arma::mat referenceData = arma::randu<arma::mat>(10, 200);
+
+  // Build all the possible models.
+  KNNModel models[14];
+  models[0] = KNNModel(KNNModel::TreeTypes::KD_TREE, true);
+  models[1] = KNNModel(KNNModel::TreeTypes::KD_TREE, false);
+  models[2] = KNNModel(KNNModel::TreeTypes::COVER_TREE, true);
+  models[3] = KNNModel(KNNModel::TreeTypes::COVER_TREE, false);
+  models[4] = KNNModel(KNNModel::TreeTypes::R_TREE, true);
+  models[5] = KNNModel(KNNModel::TreeTypes::R_TREE, false);
+  models[6] = KNNModel(KNNModel::TreeTypes::R_STAR_TREE, true);
+  models[7] = KNNModel(KNNModel::TreeTypes::R_STAR_TREE, false);
+  models[8] = KNNModel(KNNModel::TreeTypes::X_TREE, true);
+  models[9] = KNNModel(KNNModel::TreeTypes::X_TREE, false);
+  models[10] = KNNModel(KNNModel::TreeTypes::BALL_TREE, true);
+  models[11] = KNNModel(KNNModel::TreeTypes::BALL_TREE, false);
+  models[12] = KNNModel(KNNModel::TreeTypes::HILBERT_R_TREE, true);
+  models[13] = KNNModel(KNNModel::TreeTypes::HILBERT_R_TREE, false);
+
+  for (size_t j = 0; j < 2; ++j)
+  {
+    // Get a baseline.
+    KNN exact(referenceData);
+    arma::Mat<size_t> neighborsExact;
+    arma::mat distancesExact;
+    exact.Search(3, neighborsExact, distancesExact);
+
+    for (size_t i = 0; i < 14; ++i)
+    {
+      // We only have a std::move() constructor... so copy the data.
+      arma::mat referenceCopy(referenceData);
+      if (j == 0)
+        models[i].BuildModel(std::move(referenceCopy), 20, false, false, 0.05);
+      if (j == 1)
+        models[i].BuildModel(std::move(referenceCopy), 20, false, true, 0.05);
+
+      arma::Mat<size_t> neighborsAprox;
+      arma::mat distancesAprox;
+
+      models[i].Search(3, neighborsAprox, distancesAprox);
+
+      BOOST_REQUIRE_EQUAL(neighborsAprox.n_rows, neighborsExact.n_rows);
+      BOOST_REQUIRE_EQUAL(neighborsAprox.n_cols, neighborsExact.n_cols);
+      BOOST_REQUIRE_EQUAL(neighborsAprox.n_elem, neighborsExact.n_elem);
+      BOOST_REQUIRE_EQUAL(distancesAprox.n_rows, distancesExact.n_rows);
+      BOOST_REQUIRE_EQUAL(distancesAprox.n_cols, distancesExact.n_cols);
+      BOOST_REQUIRE_EQUAL(distancesAprox.n_elem, distancesExact.n_elem);
+      for (size_t k = 0; k < distancesAprox.n_elem; ++k)
+        REQUIRE_RELATIVE_ERR(distancesAprox[k], distancesExact[k], 0.05);
+    }
+  }
+}
+
+BOOST_AUTO_TEST_SUITE_END();

src/mlpack/tests/arma_extend_test.cpp

@@ -7,7 +7,7 @@
 
 #include <mlpack/core.hpp>
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace arma;

src/mlpack/tests/armadillo_svd_test.cpp

@@ -2,7 +2,7 @@
 #include <mlpack/methods/cf/svd_wrapper.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 BOOST_AUTO_TEST_SUITE(ArmadilloSVDTest);
 

src/mlpack/tests/aug_lagrangian_test.cpp

@@ -10,7 +10,7 @@
 #include <mlpack/core/optimizers/aug_lagrangian/aug_lagrangian.hpp>
 #include <mlpack/core/optimizers/aug_lagrangian/aug_lagrangian_test_functions.hpp>
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::optimization;

src/mlpack/tests/binarize_test.cpp

@@ -9,7 +9,7 @@
 #include <mlpack/core/math/random.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace arma;

src/mlpack/tests/cf_test.cpp

@@ -10,7 +10,7 @@
 #include <iostream>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 #include "serialization.hpp"
 
 BOOST_AUTO_TEST_SUITE(CFTest);

src/mlpack/tests/cli_test.cpp

@@ -22,7 +22,7 @@
 #define DEFAULT_INT 42
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 #define BASH_RED "\033[0;31m"
 #define BASH_GREEN "\033[0;32m"

src/mlpack/tests/convolution_test.cpp

@@ -13,7 +13,7 @@
 #include <mlpack/methods/ann/convolution_rules/svd_convolution.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::ann;

src/mlpack/tests/convolutional_network_test.cpp

@@ -23,7 +23,7 @@
 #include <mlpack/methods/ann/cnn.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::ann;

src/mlpack/tests/cosine_tree_test.cpp

@@ -9,7 +9,7 @@
 #include <mlpack/core/tree/cosine_tree/cosine_tree.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 BOOST_AUTO_TEST_SUITE(CosineTreeTest);
 

src/mlpack/tests/decision_stump_test.cpp

@@ -8,7 +8,7 @@
 #include <mlpack/methods/decision_stump/decision_stump.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::decision_stump;

src/mlpack/tests/det_test.cpp

@@ -7,7 +7,7 @@
  */
 #include <mlpack/core.hpp>
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 // This trick does not work on Windows.  We will have to comment out the tests
 // that depend on it.

src/mlpack/tests/distribution_test.cpp

@@ -7,7 +7,7 @@
 #include <mlpack/core.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::distribution;

src/mlpack/tests/emst_test.cpp

@@ -6,7 +6,7 @@
 #include <mlpack/core.hpp>
 #include <mlpack/methods/emst/dtb.hpp>
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 #include <mlpack/core/tree/cover_tree.hpp>
 

src/mlpack/tests/fastmks_test.cpp

@@ -9,7 +9,7 @@
 #include <mlpack/methods/fastmks/fastmks_model.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 #include "serialization.hpp"
 
 using namespace mlpack;

src/mlpack/tests/feedforward_network_test.cpp

@@ -24,7 +24,7 @@
 #include <mlpack/core/optimizers/rmsprop/rmsprop.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::ann;

src/mlpack/tests/gmm_test.cpp

@@ -15,7 +15,7 @@
 #include <mlpack/methods/gmm/eigenvalue_ratio_constraint.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::gmm;

src/mlpack/tests/hmm_test.cpp

@@ -8,7 +8,7 @@
 #include <mlpack/methods/gmm/gmm.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::hmm;

src/mlpack/tests/hoeffding_tree_test.cpp

@@ -12,7 +12,7 @@
 #include <mlpack/methods/hoeffding_trees/binary_numeric_split.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 #include "serialization.hpp"
 
 #include <stack>

src/mlpack/tests/ind2sub_test.cpp

@@ -6,7 +6,7 @@
  */
 #include <mlpack/core.hpp>
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 BOOST_AUTO_TEST_SUITE(ind2subTest);
 

src/mlpack/tests/init_rules_test.cpp

@@ -14,7 +14,7 @@
 #include <mlpack/methods/ann/init_rules/zero_init.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::ann;

src/mlpack/tests/kernel_pca_test.cpp

@@ -10,7 +10,7 @@
 #include <mlpack/methods/kernel_pca/kernel_pca.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 BOOST_AUTO_TEST_SUITE(KernelPCATest);
 

src/mlpack/tests/kernel_test.cpp

@@ -19,7 +19,7 @@
 #include <mlpack/core/metrics/mahalanobis_distance.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::kernel;

src/mlpack/tests/kernel_traits_test.cpp

@@ -9,7 +9,7 @@
 #include <mlpack/core.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::kernel;

src/mlpack/tests/kfn_test.cpp

@@ -7,7 +7,7 @@
 #include <mlpack/methods/neighbor_search/neighbor_search.hpp>
 #include <mlpack/core/tree/cover_tree.hpp>
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::neighbor;

src/mlpack/tests/kmeans_test.cpp

@@ -18,7 +18,7 @@
 #include <mlpack/methods/neighbor_search/neighbor_search.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::kmeans;

src/mlpack/tests/knn_test.cpp

@@ -10,7 +10,7 @@
 #include <mlpack/core/tree/cover_tree.hpp>
 #include <mlpack/core/tree/example_tree.hpp>
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::neighbor;
@@ -888,7 +888,9 @@ BOOST_AUTO_TEST_CASE(DualBallTreeTest)
   }
 }
 
-// Make sure sparse nearest neighbors works with kd trees.
+/**
+ * Make sure sparse nearest neighbors works with kd trees.
+ */
 BOOST_AUTO_TEST_CASE(SparseKNNKDTreeTest)
 {
   // The dimensionality of these datasets must be high so that the probability
@@ -975,7 +977,7 @@ BOOST_AUTO_TEST_CASE(KNNModelTest)
   arma::mat referenceData = arma::randu<arma::mat>(10, 200);
 
   // Build all the possible models.
-  KNNModel models[12];
+  KNNModel models[14];
   models[0] = KNNModel(KNNModel::TreeTypes::KD_TREE, true);
   models[1] = KNNModel(KNNModel::TreeTypes::KD_TREE, false);
   models[2] = KNNModel(KNNModel::TreeTypes::COVER_TREE, true);
@@ -988,6 +990,8 @@ BOOST_AUTO_TEST_CASE(KNNModelTest)
   models[9] = KNNModel(KNNModel::TreeTypes::X_TREE, false);
   models[10] = KNNModel(KNNModel::TreeTypes::BALL_TREE, true);
   models[11] = KNNModel(KNNModel::TreeTypes::BALL_TREE, false);
+  models[12] = KNNModel(KNNModel::TreeTypes::HILBERT_R_TREE, true);
+  models[13] = KNNModel(KNNModel::TreeTypes::HILBERT_R_TREE, false);
 
   for (size_t j = 0; j < 2; ++j)
   {
@@ -997,7 +1001,7 @@ BOOST_AUTO_TEST_CASE(KNNModelTest)
     arma::mat baselineDistances;
     knn.Search(queryData, 3, baselineNeighbors, baselineDistances);
 
-    for (size_t i = 0; i < 12; ++i)
+    for (size_t i = 0; i < 14; ++i)
     {
       // We only have std::move() constructors so make a copy of our data.
       arma::mat referenceCopy(referenceData);
@@ -1041,7 +1045,7 @@ BOOST_AUTO_TEST_CASE(KNNModelMonochromaticTest)
   arma::mat referenceData = arma::randu<arma::mat>(10, 200);
 
   // Build all the possible models.
-  KNNModel models[12];
+  KNNModel models[14];
   models[0] = KNNModel(KNNModel::TreeTypes::KD_TREE, true);
   models[1] = KNNModel(KNNModel::TreeTypes::KD_TREE, false);
   models[2] = KNNModel(KNNModel::TreeTypes::COVER_TREE, true);
@@ -1054,6 +1058,8 @@ BOOST_AUTO_TEST_CASE(KNNModelMonochromaticTest)
   models[9] = KNNModel(KNNModel::TreeTypes::X_TREE, false);
   models[10] = KNNModel(KNNModel::TreeTypes::BALL_TREE, true);
   models[11] = KNNModel(KNNModel::TreeTypes::BALL_TREE, false);
+  models[12] = KNNModel(KNNModel::TreeTypes::HILBERT_R_TREE, true);
+  models[13] = KNNModel(KNNModel::TreeTypes::HILBERT_R_TREE, false);
 
   for (size_t j = 0; j < 2; ++j)
   {
@@ -1063,7 +1069,7 @@ BOOST_AUTO_TEST_CASE(KNNModelMonochromaticTest)
     arma::mat baselineDistances;
     knn.Search(3, baselineNeighbors, baselineDistances);
 
-    for (size_t i = 0; i < 12; ++i)
+    for (size_t i = 0; i < 14; ++i)
     {
       // We only have a std::move() constructor... so copy the data.
       arma::mat referenceCopy(referenceData);

src/mlpack/tests/krann_search_test.cpp

@@ -10,7 +10,7 @@
 #include <mlpack/core/tree/cover_tree.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 #include <mlpack/methods/rann/ra_search.hpp>
 #include <mlpack/methods/rann/ra_model.hpp>
@@ -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[10];
+  KNNModel models[12];
   models[0] = KNNModel(KNNModel::TreeTypes::KD_TREE, false);
   models[1] = KNNModel(KNNModel::TreeTypes::KD_TREE, true);
   models[2] = KNNModel(KNNModel::TreeTypes::COVER_TREE, false);
@@ -636,13 +636,15 @@ BOOST_AUTO_TEST_CASE(RAModelTest)
   models[7] = KNNModel(KNNModel::TreeTypes::R_STAR_TREE, true);
   models[8] = KNNModel(KNNModel::TreeTypes::X_TREE, false);
   models[9] = KNNModel(KNNModel::TreeTypes::X_TREE, true);
+  models[10] = KNNModel(KNNModel::TreeTypes::HILBERT_R_TREE, false);
+  models[11] = KNNModel(KNNModel::TreeTypes::HILBERT_R_TREE, true);
 
   arma::Mat<size_t> 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 < 10; ++i)
+    for (size_t i = 0; i < 12; ++i)
     {
       // We only have std::move() constructors so make a copy of our data.
       arma::mat referenceCopy(referenceData);

src/mlpack/tests/lars_test.cpp

@@ -10,7 +10,7 @@
 #include <mlpack/methods/lars/lars.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::regression;

src/mlpack/tests/layer_traits_test.cpp

@@ -13,7 +13,7 @@
 #include <mlpack/methods/ann/layer/multiclass_classification_layer.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::ann;

src/mlpack/tests/lbfgs_test.cpp

@@ -10,7 +10,7 @@
 #include <mlpack/core/optimizers/lbfgs/test_functions.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack::optimization;
 using namespace mlpack::optimization::test;

src/mlpack/tests/lin_alg_test.cpp

@@ -10,7 +10,7 @@
 #include <mlpack/core/math/lin_alg.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace arma;
 using namespace mlpack;

src/mlpack/tests/linear_regression_test.cpp

@@ -7,7 +7,7 @@
 #include <mlpack/methods/linear_regression/linear_regression.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::regression;

src/mlpack/tests/load_save_test.cpp

@@ -9,7 +9,7 @@
 #include <mlpack/core.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::data;

src/mlpack/tests/local_coordinate_coding_test.cpp

@@ -10,7 +10,7 @@
 #include <mlpack/methods/local_coordinate_coding/lcc.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 #include "serialization.hpp"
 
 using namespace arma;

src/mlpack/tests/log_test.cpp

@@ -8,7 +8,7 @@
 #include <mlpack/core.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 

src/mlpack/tests/logistic_regression_test.cpp

@@ -9,7 +9,7 @@
 #include <mlpack/core/optimizers/sgd/sgd.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::regression;

src/mlpack/tests/lrsdp_test.cpp

@@ -8,7 +8,7 @@
 #include <mlpack/core/optimizers/sdp/lrsdp.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::optimization;

src/mlpack/tests/lsh_test.cpp

@@ -6,7 +6,7 @@
 #include <mlpack/core.hpp>
 #include <mlpack/core/metrics/lmetric.hpp>
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 #include <mlpack/methods/lsh/lsh_search.hpp>
 #include <mlpack/methods/neighbor_search/neighbor_search.hpp>

src/mlpack/tests/lstm_peephole_test.cpp

@@ -9,7 +9,7 @@
 #include <mlpack/methods/ann/layer/lstm_layer.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::ann;

src/mlpack/tests/math_test.cpp

@@ -8,7 +8,7 @@
 #include <mlpack/core/math/random.hpp>
 #include <mlpack/core/math/range.hpp>
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace math;

src/mlpack/tests/matrix_completion_test.cpp

@@ -8,7 +8,7 @@
 #include <mlpack/methods/matrix_completion/matrix_completion.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::matrix_completion;

src/mlpack/tests/maximal_inputs_test.cpp

@@ -9,7 +9,7 @@
 #include <mlpack/methods/sparse_autoencoder/maximal_inputs.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 

src/mlpack/tests/mean_shift_test.cpp

@@ -8,7 +8,7 @@
 #include <mlpack/methods/mean_shift/mean_shift.hpp>
 
 #include <boost/test/unit_test.hpp>
-#include "old_boost_test_definitions.hpp"
+#include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::meanshift;