wrote emst_tests, fixed a few problems created by blindly replacing index_t with size_t

src/mlpack/methods/emst/CMakeLists.txt

@@ -33,3 +33,11 @@ target_link_libraries(union_find_test
   mlpack
   boost_unit_test_framework
 )
+
+add_executable(emst_test
+  emst_test.cc
+)
+target_link_libraries(emst_test
+  mlpack
+  boost_unit_test_framework
+)

src/mlpack/methods/emst/README

@@ -1,19 +0,0 @@
-This folder contains an implementation of the DualTreeBoruvka algorithm for finding a Euclidean Minimum Spanning Tree.  NOTE: this algorithm is awaiting publication and is not for distribution.  
-
-Compile with fl-build emst_main
-
-Run with ./emst_main --data="filename" for basic use.  This will create a file "output.txt" with the minimum spanning tree represented as an edge list.  
-
-Other options:
-
-bool --using_thor -> defaults to 0, thor is not yet supported, so setting it to 1 will cause an error message to be printed
-
-string --data -> the name of the file with the data points
-
-int --dtb/leaf_size -> defaults to 1, the number of points in the leaves of the tree.  For the DualTreeBoruvka algorithm, 1 gives the fastest performance.  I recommend only changing this parameter to conserve memory.
-
-bool --do_naive -> defaults to 0.  If it is 1, then the algorithm will compute the MST with both DualTreeBoruvka and a naive implementation of Boruvka's algorithm.  It will compare the two results and exit with an error if they are different.
-
-string --naive/output_filename -> The name of the file where the naive edge list should be printed.  Defaults to naive_output.txt
-
-string --dtb/output_filename -> The name of the file where the DTB edge list should be printed.  Defaults to output.txt

src/mlpack/methods/emst/dtb.hpp

@@ -4,7 +4,7 @@
  * @author Bill March (march@gatech.edu)
  *
  * Contains an implementation of the DualTreeBoruvka algorithm for finding a
- * Euclidean Minimum Spanning Tree.
+ * Euclidean Minimum Spanning Tree using the kd-tree data structure.
  *
  * Citation: March, W. B.; Ram, P.; and Gray, A. G.  Fast Euclidean Minimum Spanning
  * Tree: Algorithm, Analysis, Applications.  In KDD, 2010.
@@ -21,17 +21,9 @@
 #include <mlpack/core/tree/binary_space_tree.hpp>
 #include <mlpack/core/kernels/lmetric.hpp>
 
-PARAM(size_t, "leaf_size", "Size of the leaves.", "naive", 1, false);
-
 namespace mlpack {
 namespace emst {
 
-/*
-const fx_submodule_doc dtb_submodules[] = {
-FX_SUBMODULE_DOC_DONE
-};
- */
-
 /**
 * A Stat class for use with fastlib's trees.  This one only stores two values.
  *
@@ -43,10 +35,11 @@ FX_SUBMODULE_DOC_DONE
  * points in this node.  If points in this node are in different components,
  * this value will be negative.
  */
+  
 class DTBStat {
  private:
   double max_neighbor_distance_;
-  size_t component_membership_;
+  int component_membership_;
 
  public:
   void set_max_neighbor_distance(double distance) {
@@ -57,11 +50,11 @@ class DTBStat {
     return max_neighbor_distance_;
   }
 
-  void set_component_membership(size_t membership) {
+  void set_component_membership(int membership) {
     component_membership_ = membership;
   }
 
-  size_t component_membership() {
+  int component_membership() {
     return component_membership_;
   }
 
@@ -108,8 +101,6 @@ class DTBStat {
  */
 class DualTreeBoruvka {
 
-//  FORBID_ACCIDENTAL_COPIES(DualTreeBoruvka);
-
  public:
   // For now, everything is in Euclidean space
   static const size_t metric = 2;
@@ -144,7 +135,7 @@ class DualTreeBoruvka {
   size_t number_r_recursions_;
   size_t number_both_recursions_;
 
-  int do_naive_;
+  bool do_naive_;
 
   DTBTree* tree_;
 
@@ -286,6 +277,9 @@ class DualTreeBoruvka {
       //pruned by component membership
 
       mlpack::Log::Assert(reference_node->stat().component_membership() >= 0);
+               
+      mlpack::Log::Info << query_node->stat().component_membership() << "q mem\n";
+      mlpack::Log::Info << reference_node->stat().component_membership() << "r mem\n";
 
       number_component_prunes_++;
     }
@@ -411,7 +405,7 @@ class DualTreeBoruvka {
 
   struct SortEdgesHelper_ {
     bool operator() (const EdgePair& pairA, const EdgePair& pairB) {
-      return (pairA.distance() > pairB.distance());
+      return (pairA.distance() < pairB.distance());
     }
   } SortFun;
 
@@ -432,29 +426,33 @@ class DualTreeBoruvka {
     SortEdges_();
 
     mlpack::Log::Assert(number_of_edges_ == number_of_points_ - 1);
-    results.set_size(3, number_of_edges_);
+    results.set_size(number_of_edges_, 3);
 
+    // need to unpermute the point labels
     if (!do_naive_) {
       for (size_t i = 0; i < (number_of_points_ - 1); i++) {
 
-        edges_[i].set_lesser_index(old_from_new_permutation_[edges_[i]
-          .lesser_index()]);
+        // Make sure the edge list stores the smaller index first to 
+        // make checking correctness easier
+        size_t ind1, ind2;
+        ind1 = old_from_new_permutation_[edges_[i].lesser_index()];
+        ind2 = old_from_new_permutation_[edges_[i].greater_index()];
+        
+        edges_[i].set_lesser_index(std::min(ind1, ind2));
+        edges_[i].set_greater_index(std::max(ind1, ind2));
 
-        edges_[i].set_greater_index(old_from_new_permutation_[edges_[i]
-          .greater_index()]);
-
-        results(0, i) = edges_[i].lesser_index();
-        results(1, i) = edges_[i].greater_index();
-        results(2, i) = sqrt(edges_[i].distance());
+        results(i, 0) = edges_[i].lesser_index();
+        results(i, 1) = edges_[i].greater_index();
+        results(i, 2) = sqrt(edges_[i].distance());
 
       }
     }
     else {
 
       for (size_t i = 0; i < number_of_edges_; i++) {
-        results(0, i) = edges_[i].lesser_index();
-        results(1, i) = edges_[i].greater_index();
-        results(2, i) = sqrt(edges_[i].distance());
+        results(i, 0) = edges_[i].lesser_index();
+        results(i, 1) = edges_[i].greater_index();
+        results(i, 2) = sqrt(edges_[i].distance());
       }
 
     }
@@ -532,9 +530,10 @@ class DualTreeBoruvka {
     fx_result_int(module_, "number_r_recursions", number_r_recursions_);
     fx_result_int(module_, "number_both_recursions", number_both_recursions_);*/
     // TODO, not sure how I missed this last time.
-    mlpack::Log::Info << "total_squared_length" << total_dist_ << std::endl;
-    mlpack::Log::Info << "number_of_points" << number_of_points_ << std::endl;
-    mlpack::Log::Info << "dimension" << data_points_.n_rows << std::endl;
+    mlpack::Log::Info << "Total squared length: " << total_dist_ << std::endl;
+    mlpack::Log::Info << "Number of points: " << number_of_points_ << std::endl;
+    mlpack::Log::Info << "Dimension: " << data_points_.n_rows << std::endl;
+    /*
     mlpack::Log::Info << "number_of_loops" << std::endl;
     mlpack::Log::Info << "number_distance_prunes" << std::endl;
     mlpack::Log::Info << "number_component_prunes" << std::endl;
@@ -542,6 +541,10 @@ class DualTreeBoruvka {
     mlpack::Log::Info << "number_q_recursions" << std::endl;
     mlpack::Log::Info << "number_r_recursions" << std::endl;
     mlpack::Log::Info << "number_both_recursions" << std::endl;
+     */
+    
+    mlpack::CLI::GetParam<double>("dtb/total_squared_length") = total_dist_;
+    
   } // OutputResults_
 
   /////////// Public Functions ///////////////////
@@ -554,11 +557,9 @@ class DualTreeBoruvka {
 
 
   /**
-   * Takes in a reference to the data set and a module.  Copies the data,
+   * Takes in a reference to the data set.  Copies the data,
    * builds the tree, and initializes all of the member variables.
    *
-   * This module will be checked for the optional parameters "leaf_size" and
-   * "do_naive".
    */
   void Init(const arma::mat& data) {
 
@@ -572,17 +573,20 @@ class DualTreeBoruvka {
       // This gives best pruning empirically
       // Use leaf_size=1 unless space is a big concern
       CLI::GetParam<int>("tree/leaf_size") =
-          CLI::GetParam<size_t>("naive/leaf_size");
+          CLI::GetParam<int>("emst/leaf_size");
 
-      Timers::StartTimer("naive/tree_building");
+      Timers::StartTimer("emst/tree_building");
 
       tree_ = new DTBTree(data_points_, old_from_new_permutation_);
 
-      Timers::StopTimer("naive/tree_building");
+      Timers::StopTimer("emst/tree_building");
+      
     }
     else {
+      
       tree_ = NULL;
       old_from_new_permutation_.resize(0);
+      
     }
 
     number_of_points_ = data_points_.n_cols;
@@ -614,23 +618,30 @@ class DualTreeBoruvka {
     Timers::StartTimer("emst/MST_computation");
 
     while (number_of_edges_ < (number_of_points_ - 1)) {
+      
       ComputeNeighbors_();
 
       AddAllEdges_();
 
       Cleanup_();
 
-      Log::Info << "number_of_loops = " << number_of_loops_ << std::endl;
+      
+      Log::Info << "Finished loop number: " << number_of_loops_ << std::endl;
+      Log::Info << number_of_edges_ << " edges found so far.\n\n";
+      /*
+      Log::Info << number_leaf_computations_ << " base cases.\n";
+      Log::Info << number_distance_prunes_ << " distance prunes.\n";
+      Log::Info << number_component_prunes_ << " component prunes.\n";
+      Log::Info << number_r_recursions_ << " reference recursions.\n";
+      Log::Info << number_q_recursions_ << " query recursions.\n";
+      Log::Info << number_both_recursions_ << " dual recursions.\n\n";
+      */
+      
     }
 
     Timers::StopTimer("emst/MST_computation");
 
-//    if (results != NULL) {
-
-      EmitResults_(results);
-
-//    }
-
+    EmitResults_(results);
 
     OutputResults_();
 

src/mlpack/methods/emst/emst_main.cpp

@@ -3,149 +3,78 @@
  *
  * Calls the DualTreeBoruvka algorithm from dtb.h
  * Can optionally call Naive Boruvka's method
- * See README for command line options.
+ * 
+ * For algorithm details, see: 
+ * March, W.B., Ram, P., and Gray, A.G.
+ * Fast Euclidean Minimum Spanning Tree: Algorithm, Analysis, Applications. 
+ * In KDD, 2010.
  *
  * @author Bill March (march@gatech.edu)
 */
 
-
-#include <mlpack/core.h>
 #include "dtb.hpp"
 
-PARAM_FLAG("using_thor", "For when an implementation of thor is around",
-    "emst");
-PARAM_STRING_REQ("input_file", "Data input file.", "emst");
-PARAM_STRING("output_file", "Data output file.", "emst", "emst_output.csv");
+#include <mlpack/core.h>
 
-PARAM_FLAG("do_naive", "Check against naive.", "naive");
+PARAM_STRING_REQ("input_file", "Data input file.", "emst");
+PARAM_STRING("output_file", "Data output file.  Stored as an edge list.", "emst", "emst_output.csv");
+PARAM_FLAG("do_naive", "Compute the MST using .", "naive");
 PARAM_STRING("output_file", "Naive data output file.", "naive",
     "naive_output.csv");
-
-PARAM(double, "total_squared_length", "Calculation result.", "dtb", 0.0, false);
+PARAM_INT("leaf_size", "Leaf size in the kd-tree.  Singleton leaves give the empirically best performance at the cost of greater memory requirements.", "emst", 1);
+PARAM_DOUBLE("total_squared_length", "Squared length of the computed tree.", "dtb", 0.0);
 
 using namespace mlpack;
 using namespace mlpack::emst;
 
 int main(int argc, char* argv[]) {
+
   CLI::ParseCommandLine(argc, argv);
 
-  // For when I implement a thor version
-  bool using_thor = CLI::GetParam<bool>("emst/using_thor");
 
+  ///////////////// READ IN DATA //////////////////////////////////
+  std::string data_file_name = CLI::GetParam<std::string>("emst/input_file");
 
-  if (using_thor) {
-    Log::Warn << "thor is not yet supported" << std::endl;
+  Log::Info << "Reading in data.\n";
+  
+  arma::mat data_points;
+  data_points.load(data_file_name.c_str());
+  
+  // Do naive
+  if (CLI::GetParam<bool>("naive/do_naive")) {
+    
+    Log::Info << "Running naive algorithm.\n";
+    
+    DualTreeBoruvka naive;
+    //CLI::GetParam<bool>("naive/do_naive") = true;
+    
+    naive.Init(data_points);
+    
+    arma::mat naive_results;
+    naive.ComputeMST(naive_results);
+    
+    std::string naive_output_filename =
+    CLI::GetParam<std::string>("naive/output_file");
+    
+    naive_results.save(naive_output_filename.c_str(), arma::csv_ascii, false,
+                       true);
+    
   }
   else {
-
-    ///////////////// READ IN DATA //////////////////////////////////
-
-    std::string data_file_name = CLI::GetParam<std::string>("emst/input_file");
-
-    arma::mat data_points;
-    data_points.load(data_file_name.c_str());
+  
+    Log::Info << "Data read, building tree.\n";
 
     /////////////// Initialize DTB //////////////////////
     DualTreeBoruvka dtb;
-
+    dtb.Init(data_points);
+    
+    Log::Info << "Tree built, running algorithm.\n\n";
+    
     ////////////// Run DTB /////////////////////
     arma::mat results;
 
     dtb.ComputeMST(results);
 
-    //////////////// Check against naive //////////////////////////
-    if (CLI::GetParam<bool>("naive/do_naive")) {
-
-      DualTreeBoruvka naive;
-      CLI::GetParam<bool>("naive/do_naive") = true;
-
-      naive.Init(data_points);
-
-      arma::mat naive_results;
-      naive.ComputeMST(naive_results);
-
-      /* Compare the naive output to the DTB output */
-
-      Timers::StartTimer("naive/comparison");
-
-
-      // Check if the edge lists are the same
-      // Loop over the naive edge list
-      int is_correct = 1;
-      /*
-      for (size_t naive_index = 0; naive_index < results.size();
-           naive_index++) {
-
-        int this_loop_correct = 0;
-        size_t naive_lesser_index = results[naive_index].lesser_index();
-        size_t naive_greater_index = results[naive_index].greater_index();
-        double naive_distance = results[naive_index].distance();
-
-        // Loop over the DTB edge list and compare against naive
-        // Break when an edge is found that matches the current naive edge
-        for (size_t dual_index = 0; dual_index < naive_results.size();
-             dual_index++) {
-
-          size_t dual_lesser_index = results[dual_index].lesser_index();
-          size_t dual_greater_index = results[dual_index].greater_index();
-          double dual_distance = results[dual_index].distance();
-
-          if (naive_lesser_index == dual_lesser_index) {
-            if (naive_greater_index == dual_greater_index) {
-              DEBUG_ASSERT(naive_distance == dual_distance);
-              this_loop_correct = 1;
-              break;
-            }
-          }
-
-        }
-
-        if (this_loop_correct == 0) {
-          is_correct = 0;
-          break;
-        }
-
-      }
-     */
-      if (is_correct == 0) {
-
-        Log::Warn << "Naive check failed!" << std::endl <<
-        "Edge lists are different." << std::endl << std::endl;
-
-        // Check if the outputs have the same length
-        if (CLI::GetParam<double>("naive/total_squared_length") !=
-           CLI::GetParam<double>("naive/total_squared_length")) {
-
-          Log::Fatal << "Total lengths are different! "
-             << " One algorithm has failed." << std::endl;
-
-          return 1;
-        }
-        else {
-          // NOTE: if the edge lists are different, but the total lengths are
-          // the same, the algorithm may still be correct.  The MST is not
-          // uniquely defined for some point sets.  For example, an equilateral
-          // triangle has three minimum spanning trees.  It is possible for
-          // naive and DTB to find different spanning trees in this case.
-          Log::Info << "Total lengths are the same.";
-          Log::Info << "It is possible the point set";
-          Log::Info << "has more than one minimum spanning tree." << std::endl;
-        }
-
-      }
-      else {
-        Log::Info << "Naive and DualTreeBoruvka produced the same MST." <<
-          std::endl << std::endl;
-      }
-
-      Timers::StopTimer("naive/comparison");
-
-      std::string naive_output_filename =
-          CLI::GetParam<std::string>("naive/output_file");
-
-      naive_results.save(naive_output_filename.c_str(), arma::csv_ascii, false,
-          true);
-    }
 
     //////////////// Output the Results ////////////////
 
@@ -153,8 +82,8 @@ int main(int argc, char* argv[]) {
         CLI::GetParam<std::string>("emst/output_file");
 
     results.save(output_filename.c_str(), arma::csv_ascii, false, true);
-
-  }// end else (if using_thor)
+    
+  }
 
   return 0;
 

src/mlpack/methods/emst/emst_test.cc

@@ -0,0 +1,138 @@
+/**
+ * @file emst_test.cc
+ *
+ * Test file for EMST methods
+ */
+#include <mlpack/core.h>
+#include "dtb.hpp"
+
+#define BOOST_TEST_MODULE EMST Test
+#include <boost/test/unit_test.hpp>
+
+using namespace mlpack;
+using namespace mlpack::emst;
+
+/***
+ * Simple emst test with small, synthetic dataset.  This is an
+ * exhaustive test, which checks that each method for performing the calculation
+ * (dual-tree, single-tree, naive) produces the correct results.  The dataset is
+ * in one dimension for simplicity -- the correct functionality of distance
+ * functions is not tested here.
+ */
+BOOST_AUTO_TEST_CASE(exhaustive_synthetic_test) {
+  // Set up our data.
+  arma::mat data(1, 11);
+  data[0] = 0.05; // Row addressing is unnecessary (they are all 0).
+  data[1] = 0.37;
+  data[2] = 0.15;
+  data[3] = 1.25;
+  data[4] = 5.05;
+  data[5] = -0.22;
+  data[6] = -2.00;
+  data[7] = -1.30;
+  data[8] = 0.45;
+  data[9] = 0.91;
+  data[10] = 1.00;
+
+
+  // Now perform the actual calculation.
+  arma::mat results;
+  
+  
+  DualTreeBoruvka dtb;
+  dtb.Init(data);
+  dtb.ComputeMST(results);
+
+  // Now the exhaustive check for correctness.  
+
+  BOOST_REQUIRE(results(0, 0) == 1);
+  BOOST_REQUIRE(results(0, 1) == 8);
+  BOOST_REQUIRE_CLOSE(results(0, 2), 0.08, 1e-5);
+  
+  BOOST_REQUIRE(results(1, 0) == 9);
+  BOOST_REQUIRE(results(1, 1) == 10);
+  BOOST_REQUIRE_CLOSE(results(1, 2), 0.09, 1e-5);
+  
+  BOOST_REQUIRE(results(2, 0) == 0);
+  BOOST_REQUIRE(results(2, 1) == 2);
+  BOOST_REQUIRE_CLOSE(results(2, 2), 0.1, 1e-5);
+  
+  BOOST_REQUIRE(results(3, 0) == 1);
+  BOOST_REQUIRE(results(3, 1) == 2);
+  BOOST_REQUIRE_CLOSE(results(3, 2), 0.22, 1e-5);
+  
+  BOOST_REQUIRE(results(4, 0) == 3);
+  BOOST_REQUIRE(results(4, 1) == 10);
+  BOOST_REQUIRE_CLOSE(results(4, 2), 0.25, 1e-5);
+  
+  BOOST_REQUIRE(results(5, 0) == 0);
+  BOOST_REQUIRE(results(5, 1) == 5);
+  BOOST_REQUIRE_CLOSE(results(5, 2), 0.27, 1e-5);
+  
+  BOOST_REQUIRE(results(6, 0) == 8);
+  BOOST_REQUIRE(results(6, 1) == 9);
+  BOOST_REQUIRE_CLOSE(results(6, 2), 0.46, 1e-5);
+  
+  BOOST_REQUIRE(results(7, 0) == 6);
+  BOOST_REQUIRE(results(7, 1) == 7);
+  BOOST_REQUIRE_CLOSE(results(7, 2), 0.7, 1e-5);
+  
+  BOOST_REQUIRE(results(8, 0) == 5);
+  BOOST_REQUIRE(results(8, 1) == 7);
+  BOOST_REQUIRE_CLOSE(results(8, 2), 1.08, 1e-5);
+  
+  BOOST_REQUIRE(results(9, 0) == 3);
+  BOOST_REQUIRE(results(9, 1) == 4);
+  BOOST_REQUIRE_CLOSE(results(9, 2), 3.8, 1e-5);
+  
+  
+}
+
+/**
+ * Test the dual tree method against the naive computation.
+ *
+ * Errors are produced if the results are not identical.
+ */
+BOOST_AUTO_TEST_CASE(dual_tree_vs_naive) {
+  arma::mat input_data;
+
+  // Hard-coded filename: bad!
+  // Code duplication: also bad!
+  if (!input_data.load("test_data_3_1000.csv", arma::auto_detect, false,
+      true))
+    BOOST_FAIL("Cannot load test dataset test_data_3_1000.csv!");
+
+  // Set up matrices to work with (may not be necessary with no ALIAS_MATRIX?).
+  arma::mat dual_data = arma::trans(input_data);
+  arma::mat naive_data = arma::trans(input_data);
+
+  // Reset parameters from last test.
+  DualTreeBoruvka dtb;
+  dtb.Init(dual_data);
+  
+  arma::mat dual_results;
+  dtb.ComputeMST(dual_results);
+  
+  // Set naive mode.
+  CLI::GetParam<bool>("naive/do_naive") = true;
+
+  DualTreeBoruvka dtb_naive;
+  dtb_naive.Init(naive_data);
+  
+  arma::mat naive_results;
+  dtb_naive.ComputeMST(naive_results);
+  
+  BOOST_REQUIRE(dual_results.n_cols == naive_results.n_cols);
+  BOOST_REQUIRE(dual_results.n_rows == naive_results.n_rows);
+  
+  for (size_t i = 0; i < dual_results.n_rows; i++) {
+
+    BOOST_REQUIRE(dual_results(i,0) == naive_results(i,0));
+    BOOST_REQUIRE(dual_results(i,1) == naive_results(i,1));
+    BOOST_REQUIRE_CLOSE(dual_results(i,2), naive_results(i,2), 1e-5);
+
+  }
+  
+  
+}
+

src/mlpack/methods/emst/test_data.csv

@@ -1,20 +0,0 @@
-0.3964647737602753 0.8404853694114252
-0.3533360972452435 0.4465834347965441
-0.3186927723118806 0.8864284332230312
-0.01558284940832877 0.5840902203172718
-0.1593686265318048 0.3837158748071943
-0.6910043733821958 0.05885891359273643
-0.899854306161604 0.1635459506303647
-0.1590715025818064 0.5330647140218545
-0.6041441897112385 0.5826990212072189
-0.2699711179070157 0.3904781954634089
-0.2934005701189513 0.7423774060339809
-0.298525606318119 0.07553807853778238
-0.4049826335833338 0.8573779427081831
-0.9419683232918992 0.6628306597899964
-0.8464757799300067 0.002755081426883521
-0.4623792450254847 0.5325960244382983
-0.78787662089292 0.2656122349713712
-0.9827522631010304 0.3067851306141804
-0.6008551364891055 0.608715653358658
-0.2124387982011875 0.8858951305876062