Modify interface and utils

4 files changed, 148 insertions, 110 deletions

diff --git a/src/Collapse/include/gudhi/Flag_complex_sparse_matrix.h b/src/Collapse/include/gudhi/Flag_complex_sparse_matrix.h
index d472bf15..e90d284d 100644
--- a/src/Collapse/include/gudhi/Flag_complex_sparse_matrix.h
+++ b/src/Collapse/include/gudhi/Flag_complex_sparse_matrix.h
@@ -16,6 +16,7 @@
 #include <gudhi/graph_simplicial_complex.h>
 
 #include <boost/functional/hash.hpp>
+#include <boost/graph/adjacency_list.hpp>
 
 #include <Eigen/Sparse>
 
@@ -32,6 +33,7 @@
 #include <list>
 #include <algorithm>
 
+
 namespace Gudhi {
 
 namespace collapse {
@@ -47,26 +49,28 @@ namespace collapse {
  * Algorithm.
  *
  * \tparam Vertex_handle type must be a signed integer type. It admits a total order <.
- * \tparam Filtration_value type for the value of the filtration function. Must be comparable with <.
+ * \tparam Filtration type for the value of the filtration function. Must be comparable with <.
  */
-template<typename Vertex_handle, typename Filtration_value>
+template<typename Vertex_handle, typename Filtration>
 class Flag_complex_sparse_matrix {
  private:
   using Edge = std::pair<Vertex_handle, Vertex_handle>;  // This is an ordered pair, An edge is stored with convention of the first
                                          // element being the smaller i.e {2,3} not {3,2}. However this is at the level
                                          // of row indices on actual vertex lables
-  using Filtered_edge = std::pair<Edge, Filtration_value>;
+  using Filtered_edge = std::pair<Edge, Filtration>;
 
   using Row_index = std::size_t;
 
   using Map_vertex_to_index = std::unordered_map<Vertex_handle, Row_index>;
 
-  using Sparse_row_matrix = Eigen::SparseMatrix<Filtration_value, Eigen::RowMajor>;
+  using Sparse_row_matrix = Eigen::SparseMatrix<Filtration, Eigen::RowMajor>;
 
   using Row_indices_vector = std::vector<Row_index>;
 
  public:
-  using Filtered_sorted_edge_list = std::vector<std::tuple<Filtration_value, Vertex_handle, Vertex_handle>>;
+  using Filtered_sorted_edge_list = std::vector<std::tuple<Filtration, Vertex_handle, Vertex_handle>>;
+  using Filtration_value = Filtration;
+  using Proximity_graph = Gudhi::Proximity_graph<Flag_complex_sparse_matrix>;
 
  private:
   std::unordered_map<Row_index, Vertex_handle> row_to_vertex_;
@@ -88,9 +92,9 @@ class Flag_complex_sparse_matrix {
   // Boolean vector to indicate if the index is critical or not.
   std::vector<bool> dominated_edge_indicator_;  // domination indicator
 
-  //! Stores the Map between vertices_<B>row_to_vertex_  and row indices <B>row_to_vertex_ -> row-index</B>.
+  //! Stores the Map between vertices<B>row_to_vertex_  and row indices <B>row_to_vertex_ -> row-index</B>.
   /*!
-    So, if the original simplex tree had vertices_ 0,1,4,5 <br>
+    So, if the original simplex tree had vertices 0,1,4,5 <br>
     <B>row_to_vertex_</B> would store : <br>
     \verbatim
     Values =  | 0 | 1 | 4 | 5 |
@@ -106,10 +110,10 @@ class Flag_complex_sparse_matrix {
   */
   std::unordered_map<Vertex_handle, Row_index> vertex_to_row_;
 
-  //! Stores the Sparse matrix of Filtration_value values representing the Original Simplicial Complex.
+  //! Stores the Sparse matrix of Filtration values representing the Original Simplicial Complex.
   /*!
     \code
-    Sparse_row_matrix   = Eigen::SparseMatrix<Filtration_value, Eigen::RowMajor> ;
+    Sparse_row_matrix   = Eigen::SparseMatrix<Filtration, Eigen::RowMajor> ;
     \endcode
     ;
       */
@@ -120,7 +124,7 @@ class Flag_complex_sparse_matrix {
   //! Stores <I>true</I> for dominated rows and  <I>false</I> for undominated rows.
   /*!
     Initialised to a vector of length equal to the value of the variable <B>rows</B> with all <I>false</I> values.
-    Subsequent removal of dominated vertices_ is reflected by concerned entries changing to <I>true</I> in this vector.
+    Subsequent removal of dominated vertices is reflected by concerned entries changing to <I>true</I> in this vector.
   */
   std::vector<bool> domination_indicator_;  //(domination indicator)
 
@@ -128,9 +132,9 @@ class Flag_complex_sparse_matrix {
   std::vector<Filtered_edge> f_edge_vector_;
 
 
-  //! Stores the number of vertices_ in the original Simplicial Complex.
+  //! Stores the number of vertices in the original Simplicial Complex.
   /*!
-    This stores the count of vertices_ (which is also the number of rows in the Matrix).
+    This stores the count of vertices (which is also the number of rows in the Matrix).
   */
   Row_index rows;
 
@@ -202,7 +206,7 @@ class Flag_complex_sparse_matrix {
   }
 
   template<typename FilteredEdgeInsertion>
-  void set_edge_critical(Row_index indx, Filtration_value filt, FilteredEdgeInsertion filtered_edge_insert) {
+  void set_edge_critical(Row_index indx, Filtration filt, FilteredEdgeInsertion filtered_edge_insert) {
 #ifdef DEBUG_TRACES
     std::cout << "The curent index  with filtration value " << indx << ", " << filt << " is primary critical" <<
     std::endl;
@@ -290,7 +294,7 @@ class Flag_complex_sparse_matrix {
     return common;
   }
 
-  void insert_vertex(Vertex_handle vertex, Filtration_value filt_val) {
+  void insert_vertex(Vertex_handle vertex, Filtration filt_val) {
     auto rw = vertex_to_row_.find(vertex);
     if (rw == vertex_to_row_.end()) {
       // Initializing the diagonal element of the adjency matrix corresponding to rw_b.
@@ -302,7 +306,7 @@ class Flag_complex_sparse_matrix {
     }
   }
 
-  void insert_new_edges(Vertex_handle u, Vertex_handle v, Filtration_value filt_val)
+  void insert_new_edges(Vertex_handle u, Vertex_handle v, Filtration filt_val)
   {
     // The edge must not be added before, it should be a new edge.
     insert_vertex(u, filt_val);
@@ -336,19 +340,27 @@ class Flag_complex_sparse_matrix {
     <B>domination_indicator_</B> are initialised by init() function which is
     called at the begining of this. <br>
   */
-  Flag_complex_sparse_matrix(const Filtered_sorted_edge_list& edge_t)
+  template<typename DistanceMatrix>
+  Flag_complex_sparse_matrix(const DistanceMatrix& distance_matrix)
   : rows(0) {
-    for (size_t bgn_idx = 0; bgn_idx < edge_t.size(); bgn_idx++) {
-      Vertex_handle u = std::get<1>(edge_t[bgn_idx]);
-      Vertex_handle v = std::get<2>(edge_t[bgn_idx]);
-      f_edge_vector_.push_back({{u, v}, std::get<0>(edge_t[bgn_idx])});
-      vertices_.emplace(u);
-      vertices_.emplace(v);
+    Vertex_handle num_vertices = std::distance(std::begin(distance_matrix), std::end(distance_matrix));
+
+    // This one is not part of the loop
+    vertices_.emplace(0);
+    // Only browse the lower part of the distance matrix
+    for (Vertex_handle line_index = 1; line_index < num_vertices; line_index++) {
+      for (Vertex_handle row_index = 0; row_index < line_index; row_index++) {
+#ifdef DEBUG_TRACES
+        std::cout << "Insert edge: fn[" << row_index << ", " << line_index << "] = "
+                  << distance_matrix[line_index][row_index] << std::endl;
+#endif  // DEBUG_TRACES
+        f_edge_vector_.push_back({{row_index, line_index}, distance_matrix[line_index][row_index]});
+      }
+      vertices_.emplace(line_index);
     }
   }
 
-  template<typename OneSkeletonGraph>
-  Flag_complex_sparse_matrix(const OneSkeletonGraph& one_skeleton_graph)
+  Flag_complex_sparse_matrix(const Proximity_graph& one_skeleton_graph)
   : rows(0) {
     // Insert all vertices_
     for (auto v_it = boost::vertices(one_skeleton_graph); v_it.first != v_it.second; ++v_it.first) {
@@ -362,6 +374,18 @@ class Flag_complex_sparse_matrix {
       Vertex_handle v = target(edge, one_skeleton_graph);
       f_edge_vector_.push_back({{u, v}, boost::get(Gudhi::edge_filtration_t(), one_skeleton_graph, edge)});
     }
+  }
+
+  // Performs edge collapse in a decreasing sequence of the filtration value.
+  template<typename FilteredEdgeInsertion>
+  void filtered_edge_collapse(FilteredEdgeInsertion filtered_edge_insert) {
+    Row_index endIdx = 0;
+
+    u_set_removed_redges_.clear();
+    u_set_dominated_redges_.clear();
+    critical_edge_indicator_.clear();
+
+    std::cout << "Sort it - " << f_edge_vector_.size() << std::endl;
     // Sort edges
     auto sort_by_filtration = [](const Filtered_edge& edge_a, const Filtered_edge& edge_b) -> bool
     {
@@ -373,17 +397,9 @@ class Flag_complex_sparse_matrix {
 #else
     std::stable_sort(f_edge_vector_.begin(), f_edge_vector_.end(), sort_by_filtration);
 #endif
-  }
-
-  // Performs edge collapse in a decreasing sequence of the filtration value.
-  template<typename FilteredEdgeInsertion>
-  void filtered_edge_collapse(FilteredEdgeInsertion filtered_edge_insert) {
-    Row_index endIdx = 0;
-
-    u_set_removed_redges_.clear();
-    u_set_dominated_redges_.clear();
-    critical_edge_indicator_.clear();
+    std::cout << "Sorted" << std::endl;
 
+    std::cout << vertices_.size() << std::endl;
     // Initializing sparse_row_adjacency_matrix_, This is a row-major sparse matrix.
     sparse_row_adjacency_matrix_ = Sparse_row_matrix(vertices_.size(), vertices_.size());
 
@@ -392,7 +408,7 @@ class Flag_complex_sparse_matrix {
       Edge edge = std::get<0>(fec);
       Vertex_handle u = std::get<0>(edge);
       Vertex_handle v = std::get<1>(edge);
-      Filtration_value filt = std::get<1>(fec);
+      Filtration filt = std::get<1>(fec);
 
       // Inserts the edge in the sparse matrix to update the graph (G_i)
       insert_new_edges(u, v, filt);
diff --git a/src/Collapse/test/collapse_unit_test.cpp b/src/Collapse/test/collapse_unit_test.cpp
index 38adfa8a..3a07e088 100644
--- a/src/Collapse/test/collapse_unit_test.cpp
+++ b/src/Collapse/test/collapse_unit_test.cpp
@@ -11,6 +11,7 @@
 #include <iostream>
 #include <tuple>
 #include <vector>
+#include <array>
 
 #define BOOST_TEST_DYN_LINK
 #define BOOST_TEST_MODULE "collapse"
@@ -32,7 +33,7 @@ bool find_edge_in_list(const Filtered_edge& edge, const Filtered_sorted_edge_lis
   }
   return false;
 }
-
+/*
 void trace_and_check_collapse(const Filtered_sorted_edge_list& edges, const Filtered_sorted_edge_list& removed_edges) {
   std::cout << "BEFORE COLLAPSE - Total number of edges: " << edges.size() << std::endl;
   BOOST_CHECK(edges.size() > 0);
@@ -68,70 +69,104 @@ void trace_and_check_collapse(const Filtered_sorted_edge_list& edges, const Filt
 }
 
 BOOST_AUTO_TEST_CASE(collapse) {
-  /*
-      1   2
-      o---o
-      |   |
-      |   |
-      |   |
-      o---o
-      0   3
-  */
+  //  1   2
+  //  o---o
+  //  |   |
+  //  |   |
+  //  |   |
+  //  o---o
+  //  0   3
   Filtered_sorted_edge_list edges {{1., 0, 1}, {1., 1, 2}, {1., 2, 3}, {1., 3, 0}};
   trace_and_check_collapse(edges, {});
   
-  /*
-      1   2
-      o---o
-      |\ /|
-      | x |
-      |/ \|
-      o---o
-      0   3
-  */
+  //  1   2
+  //  o---o
+  //  |\ /|
+  //  | x |
+  //  |/ \|
+  //  o---o
+  //  0   3
   edges.push_back({2., 0, 2});
   edges.push_back({2., 1, 3});
   trace_and_check_collapse(edges, {{2., 1, 3}});
 
-  /*
-      1   2   4
-      o---o---o
-      |\ /|   |
-      | x |   |
-      |/ \|   |
-      o---o---o
-      0   3   5
-  */
+  //  1   2   4
+  //  o---o---o
+  //  |\ /|   |
+  //  | x |   |
+  //  |/ \|   |
+  //  o---o---o
+  //  0   3   5
   edges.push_back({3., 2, 4});
   edges.push_back({3., 4, 5});
   edges.push_back({3., 5, 3});
   trace_and_check_collapse(edges, {{2., 1, 3}});
 
-  /*
-      1   2   4
-      o---o---o
-      |\ /|\ /|
-      | x | x |
-      |/ \|/ \|
-      o---o---o
-      0   3   5
-  */
+  //  1   2   4
+  //  o---o---o
+  //  |\ /|\ /|
+  //  | x | x |
+  //  |/ \|/ \|
+  //  o---o---o
+  //  0   3   5
   edges.push_back({4., 2, 5});
   edges.push_back({4., 4, 3});
   trace_and_check_collapse(edges, {{2., 1, 3}, {4., 4, 3}});
 
-  /*
-      1   2   4
-      o---o---o
-      |\ /|\ /|
-      | x | x |  + [0,4] and [1,5]
-      |/ \|/ \|
-      o---o---o
-      0   3   5
-  */
+  //  1   2   4
+  //  o---o---o
+  //  |\ /|\ /|
+  //  | x | x |  + [0,4] and [1,5]
+  //  |/ \|/ \|
+  //  o---o---o
+  //  0   3   5
   edges.push_back({5., 1, 5});
   edges.push_back({5., 0, 4});
   trace_and_check_collapse(edges, {{2., 1, 3}, {4., 4, 3}, {5., 0, 4}});
-}
+}*/
+
+
+BOOST_AUTO_TEST_CASE(collapse_from_distance_matrix) {
+  //  1   2
+  //  o---o
+  //  |\ /|
+  //  | x |
+  //  |/ \|
+  //  o---o
+  //  0   3
+  // Lower diagonal distance matrix
+  std::array<std::array<double, 4>, 4> distance_matrix = {{{0., 0., 0., 0.},
+                                                           {1., 0., 0., 0.},
+                                                           {2., 1., 0., 0.},
+                                                           {1., 2., 1., 0.} }};
+  
+  std::cout << "COLLAPSE - keep edges: " << std::endl;
+  Flag_complex_sparse_matrix flag_complex_sparse_matrix(distance_matrix);
+  Filtered_sorted_edge_list collapse_edges;
+  flag_complex_sparse_matrix.filtered_edge_collapse(
+    [&collapse_edges](std::pair<Vertex_handle, Vertex_handle> edge, Filtration_value filtration) {
+      std::cout << "f[" << std::get<0>(edge) << ", " << std::get<1>(edge) << "] = " << filtration << std::endl;
+        collapse_edges.push_back({filtration, std::get<0>(edge), std::get<1>(edge)});
+      });
+  std::cout << "AFTER COLLAPSE - Total number of edges: " << collapse_edges.size() << std::endl;
+  BOOST_CHECK(collapse_edges.size() == 5);
 
+  Filtered_sorted_edge_list edges {{1., 0, 1}, {1., 1, 2}, {1., 2, 3}, {1., 0, 3}, {2., 0, 2}, {2., 1, 3}};
+  
+  for (auto edge_from_collapse : collapse_edges) {
+    std::cout << "f[" << std::get<1>(edge_from_collapse) << ", " << std::get<2>(edge_from_collapse) << "] = "
+              << std::get<0>(edge_from_collapse) << std::endl;
+    // Check each edge from collapse is in the input
+    BOOST_CHECK(find_edge_in_list(edge_from_collapse, edges));
+  }
 
+  Filtered_sorted_edge_list removed_edges {{2., 1, 3}};
+  
+  std::cout << "CHECK COLLAPSE - Total number of removed edges: " << removed_edges.size() << std::endl;
+  for (auto removed_edge : removed_edges) {
+    std::cout << "f[" << std::get<1>(removed_edge) << ", " << std::get<2>(removed_edge) << "] = "
+              << std::get<0>(removed_edge) << std::endl;
+    // Check each removed edge from collapse is in the input
+    BOOST_CHECK(!find_edge_in_list(removed_edge, collapse_edges));
+  }
+}
\ No newline at end of file
diff --git a/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp b/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp
index f6926224..f4a460ab 100644
--- a/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp
+++ b/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp
@@ -11,10 +11,8 @@
 #include <gudhi/Flag_complex_sparse_matrix.h>
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Persistent_cohomology.h>
-#include <gudhi/Rips_edge_list.h>
-#include <gudhi/distance_functions.h>
 #include <gudhi/reader_utils.h>
-#include <gudhi/Points_off_io.h>
+#include <gudhi/graph_simplicial_complex.h>
 
 #include <boost/program_options.hpp>
 
@@ -23,8 +21,8 @@ using Filtration_value = Simplex_tree::Filtration_value;
 using Vertex_handle = Simplex_tree::Vertex_handle;
 
 using Flag_complex_sparse_matrix = Gudhi::collapse::Flag_complex_sparse_matrix<Vertex_handle, Filtration_value>;
+using Proximity_graph = Flag_complex_sparse_matrix::Proximity_graph;
 
-using Rips_edge_list = Gudhi::rips_edge_list::Rips_edge_list<Filtration_value>;
 using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
 using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>;
 using Distance_matrix = std::vector<std::vector<Filtration_value>>;
@@ -82,28 +80,22 @@ int main(int argc, char* argv[]) {
   program_options(argc, argv, csv_matrix_file, filediag, threshold, dim_max, p, min_persistence);
 
   Distance_matrix distances;
-  Distance_matrix sparse_distances;
 
   distances = Gudhi::read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_matrix_file);
   std::cout << "Read the distance matrix succesfully, of size: " << distances.size() << std::endl;
 
-  Flag_complex_sparse_matrix::Filtered_sorted_edge_list edge_t;
-
-  Rips_edge_list Rips_edge_list_from_file(distances, threshold);
-  Rips_edge_list_from_file.create_edges(edge_t);
-  
-  if (edge_t.size() <= 0) {
-    std::cerr << "Total number of egdes are zero." << std::endl;
-    exit(-1);
-  }
-
-  std::cout << "Total number of edges before collapse are: " << edge_t.size() << std::endl;
+  Proximity_graph proximity_graph = Gudhi::compute_proximity_graph<Simplex_tree>(boost::irange((size_t)0,
+                                                                                               distances.size()),
+                                                                                 threshold,
+                                                                                 [&distances](size_t i, size_t j) {
+                                                                                   return distances[j][i];
+                                                                                 });
 
   // Now we will perform filtered edge collapse to sparsify the edge list edge_t.
-  Flag_complex_sparse_matrix mat_filt_edge_coll(edge_t);
+  Flag_complex_sparse_matrix flag_complex(proximity_graph);
 
   Simplex_tree stree;
-  mat_filt_edge_coll.filtered_edge_collapse(
+  flag_complex.filtered_edge_collapse(
     [&stree](std::vector<Vertex_handle> edge, Filtration_value filtration) {
         // insert the 2 vertices with a 0. filtration value just like a Rips
         stree.insert_simplex({edge[0]}, 0.);
diff --git a/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp b/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp
index e322d3cd..b9130d4c 100644
--- a/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp
+++ b/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp
@@ -12,11 +12,9 @@
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Persistent_cohomology.h>
 #include <gudhi/distance_functions.h>
-#include <gudhi/reader_utils.h>
 #include <gudhi/Points_off_io.h>
 #include <gudhi/graph_simplicial_complex.h>
 
-#include <boost/graph/adjacency_list.hpp>
 #include <boost/program_options.hpp>
 
 #include<utility>  // for std::pair
@@ -31,15 +29,12 @@ using Point = std::vector<Filtration_value>;
 using Vector_of_points = std::vector<Point>;
 
 using Flag_complex_sparse_matrix = Gudhi::collapse::Flag_complex_sparse_matrix<Vertex_handle, Filtration_value>;
+using Proximity_graph = Flag_complex_sparse_matrix::Proximity_graph;
 
 using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
 using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>;
 using Distance_matrix = std::vector<std::vector<Filtration_value>>;
 
-using Adjacency_list = boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS,
-                                             boost::property<Gudhi::vertex_filtration_t, double>,
-                                             boost::property<Gudhi::edge_filtration_t, double>>;
-
 void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag,
                      Filtration_value& threshold, int& dim_max, int& p, Filtration_value& min_persistence);
 
@@ -76,9 +71,9 @@ int main(int argc, char* argv[]) {
   std::cout << "Successfully read " << point_vector.size() << " point_vector.\n";
   std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n";
 
-  Adjacency_list proximity_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(),
-                                                                                threshold,
-                                                                                Gudhi::Euclidean_distance());
+  Proximity_graph proximity_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(),
+                                                                                 threshold,
+                                                                                 Gudhi::Euclidean_distance());
 
   if (num_edges(proximity_graph) <= 0) {
     std::cerr << "Total number of egdes are zero." << std::endl;
@@ -89,7 +84,7 @@ int main(int argc, char* argv[]) {
 
   Simplex_tree stree;
   mat_filt_edge_coll.filtered_edge_collapse(
-    [&stree](std::vector<Vertex_handle> edge, Filtration_value filtration) {
+    [&stree](const std::vector<Vertex_handle>& edge, Filtration_value filtration) {
         // insert the 2 vertices with a 0. filtration value just like a Rips
         stree.insert_simplex({edge[0]}, 0.);
         stree.insert_simplex({edge[1]}, 0.);