Try to add inf values

2 files changed, 171 insertions, 0 deletions

diff --git a/src/Collapse/example/CMakeLists.txt b/src/Collapse/example/CMakeLists.txt
index 96354489..95e1696f 100644
--- a/src/Collapse/example/CMakeLists.txt
+++ b/src/Collapse/example/CMakeLists.txt
@@ -21,3 +21,10 @@ add_test(NAME Edge_collapse_conserve_persistence_1 COMMAND $<TARGET_FILE:Edge_co
 
          add_test(NAME Edge_collapse_conserve_persistence_2 COMMAND $<TARGET_FILE:Edge_collapse_conserve_persistence>
          "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2.")
+
+# Tests
+add_executable ( inf edge_collapse_conserve_persistence_with_inf.cpp )
+
+if (TBB_FOUND)
+  target_link_libraries(inf ${TBB_LIBRARIES})
+endif()
diff --git a/src/Collapse/example/edge_collapse_conserve_persistence_with_inf.cpp b/src/Collapse/example/edge_collapse_conserve_persistence_with_inf.cpp
new file mode 100644
index 00000000..0555725f
--- /dev/null
+++ b/src/Collapse/example/edge_collapse_conserve_persistence_with_inf.cpp
@@ -0,0 +1,164 @@
+/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2020 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#include <gudhi/Flag_complex_sparse_matrix.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/Points_off_io.h>
+#include <gudhi/graph_simplicial_complex.h>
+
+#include<utility>  // for std::pair
+#include<vector>
+#include<tuple>
+#include <limits>  // for numeric_limits<>
+
+// Types definition
+
+using Simplex_tree = Gudhi::Simplex_tree<>;
+using Filtration_value = Simplex_tree::Filtration_value;
+using Vertex_handle = Simplex_tree::Vertex_handle;
+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 Persistence_pair = std::tuple<int, Filtration_value, Filtration_value>;
+/*
+ * Compare two intervals by dimension, then by length.
+ */
+struct cmp_intervals_by_dim_then_length {
+  explicit cmp_intervals_by_dim_then_length(Simplex_tree * sc)
+      : sc_(sc) { }
+
+  template<typename Persistent_interval>
+  bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) {
+    if (sc_->dimension(get < 0 > (p1)) == sc_->dimension(get < 0 > (p2)))
+      return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1))
+              > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2)));
+    else
+      return (sc_->dimension(get < 0 > (p1)) > sc_->dimension(get < 0 > (p2)));
+  }
+  Simplex_tree* sc_;
+};
+
+std::vector<Persistence_pair> get_persistence_pairs(Simplex_tree& st, int ambient_dim) {
+  std::vector<Persistence_pair> ppairs;
+  std::cout << "expansion" << std::endl;
+  st.expansion(ambient_dim);
+  
+  // Sort the simplices in the order of the filtration
+  std::cout << "initialize_filtration" << std::endl;
+  st.initialize_filtration();
+  std::cout << "Persistent_cohomology" << std::endl;
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(st);
+  // initializes the coefficient field for homology - must be a prime number
+  int p = 11;
+  std::cout << "init_coefficients" << std::endl;
+  pcoh.init_coefficients(p);
+
+  // Default min_interval_length = 0.
+  std::cout << "compute_persistent_cohomology" << std::endl;
+  pcoh.compute_persistent_cohomology();
+  // Custom sort and output persistence
+  cmp_intervals_by_dim_then_length cmp(&st);
+  auto persistent_pairs = pcoh.get_persistent_pairs();
+  std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp);
+  for (auto pair : persistent_pairs) {
+    ppairs.push_back({st.dimension(get<0>(pair)),
+                      st.filtration(get<0>(pair)),
+                      st.filtration(get<1>(pair)) });
+  }
+  return ppairs;
+}
+
+int main(int argc, char* argv[]) {
+  if (argc != 3) {
+    std::cerr << "This program requires an OFF file and minimal threshold value as parameter\n";
+    std::cerr << "For instance: ./Edge_collapse_conserve_persistence ../../data/points/tore3D_300.off 1.\n";
+    exit(-1);  // ----- >>
+  }
+
+  std::string off_file_points {argv[1]};
+  double threshold {atof(argv[2])};
+
+  Gudhi::Points_off_reader<Point> off_reader(off_file_points);
+  if (!off_reader.is_valid()) {
+    std::cerr << "Unable to read file " << off_file_points << "\n";
+    exit(-1);  // ----- >>
+  }
+
+  Vector_of_points point_vector = off_reader.get_point_cloud();
+  if (point_vector.size() <= 0) {
+    std::cerr << "Empty point cloud." << std::endl;
+    exit(-1);  // ----- >>
+  }
+
+  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;
+    exit(-1);
+  }
+
+  int ambient_dim = point_vector[0].size();
+
+  // ***** Simplex tree from a flag complex built after collapse *****
+  Flag_complex_sparse_matrix mat_filt_edge_coll(proximity_graph);
+
+  Simplex_tree stree_from_collapse;
+  for (Vertex_handle i = 0; i < point_vector.size(); i++) {
+    // insert the 2 vertices with a 0. filtration value just like a Rips
+    stree_from_collapse.insert_simplex({i}, 0.);
+    for (Vertex_handle j = 0; j < i; j++) {
+      stree_from_collapse.insert_simplex({i, j}, std::numeric_limits<Filtration_value>::infinity());
+    }
+  }
+  mat_filt_edge_coll.filtered_edge_collapse(
+    [&stree_from_collapse](const std::vector<Vertex_handle>& edge, Filtration_value filtration) {
+        // insert the edge
+        stree_from_collapse.assign_filtration(stree_from_collapse.find(edge), filtration);
+      });
+
+  std::vector<Persistence_pair> ppairs_from_collapse = get_persistence_pairs(stree_from_collapse, ambient_dim);
+
+  // ***** Simplex tree from the complete flag complex *****
+  Simplex_tree stree_wo_collapse;
+  stree_wo_collapse.insert_graph(proximity_graph);
+
+  std::vector<Persistence_pair> ppairs_wo_collapse = get_persistence_pairs(stree_wo_collapse, ambient_dim);
+
+  if (ppairs_wo_collapse.size() != ppairs_from_collapse.size()) {
+    std::cerr << "Number of persistence pairs with    collapse is " << ppairs_from_collapse.size() << std::endl;
+    std::cerr << "Number of persistence pairs without collapse is " << ppairs_wo_collapse.size()   << std::endl;
+    exit(-1);
+  }
+
+  int return_value = 0;
+  auto ppwoc_ptr = ppairs_wo_collapse.begin();
+  for (auto ppfc: ppairs_from_collapse) {
+    if (ppfc != *ppwoc_ptr) {
+      return_value++;
+      std::cerr << "Without collapse: "
+                << std::get<0>(*ppwoc_ptr) << " " << std::get<1>(*ppwoc_ptr) << " " << std::get<2>(*ppwoc_ptr)
+                << " - With collapse: "
+                << std::get<0>(ppfc) << " " << std::get<1>(ppfc) << " " << std::get<2>(ppfc) << std::endl;
+    }
+    ppwoc_ptr++;
+  }
+  return return_value;
+}