1 files changed, 214 insertions, 0 deletions
diff --git a/example/Persistent_cohomology/performance_rips_persistence.cpp b/example/Persistent_cohomology/performance_rips_persistence.cpp
new file mode 100644
index 00000000..b4d282ac
--- /dev/null
+++ b/example/Persistent_cohomology/performance_rips_persistence.cpp
@@ -0,0 +1,214 @@
+/*    This file is part of the Gudhi Library. The Gudhi library 
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++ 
+ *    library for computational topology.
+ *
+ *    Author(s):       Clément Maria
+ *
+ *    Copyright (C) 2014  INRIA Sophia Antipolis-Méditerranée (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <gudhi/reader_utils.h>
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Persistent_cohomology/Multi_field.h>
+#include <gudhi/Hasse_complex.h>
+
+#include <chrono>
+#include <string>
+#include <vector>
+
+using namespace Gudhi;
+using namespace Gudhi::persistent_cohomology;
+
+/* Compute the persistent homology of the complex cpx with coefficients in Z/pZ. */
+template< typename FilteredComplex>
+void timing_persistence(FilteredComplex & cpx
+                        , int p);
+
+/* Compute multi-field persistent homology of the complex cpx with coefficients in
+ * Z/rZ for all prime number r in [p;q].*/
+template< typename FilteredComplex>
+void timing_persistence(FilteredComplex & cpx
+                        , int p
+                        , int q);
+
+/* Timings for the computation of persistent homology with different 
+ * representations of a Rips complex and different coefficient fields. The 
+ * Rips complex is built on a set of 10000 points sampling a Klein bottle embedded 
+ * in dimension 5.
+ * We represent complexes with a simplex tree and 
+ * with a Hasse diagram. The Hasse diagram represents explicitly all 
+ * codimension 1 incidence relations in the complex, and hence leads to 
+ * a faster computation of persistence because boundaries are precomputed. 
+ * Hovewer, the simplex tree may be constructed directly from a point cloud and
+ * is more compact.
+ * We compute persistent homology with coefficient fields Z/2Z and Z/1223Z.
+ * We present also timings for the computation of multi-field persistent 
+ * homology in all fields Z/rZ for r prime between 2 and 1223.
+ */
+int main(int argc, char * argv[]) {
+  std::chrono::time_point<std::chrono::system_clock> start, end;
+  int elapsed_sec;
+  {
+
+  std::string filepoints = "../../../data/points/Kl.txt";
+  Filtration_value threshold = 0.27;
+  int dim_max = 3;
+  int p = 2;
+  int q = 1223;
+
+  // Extract the points from the file filepoints
+  typedef std::vector<double> Point_t;
+  std::vector< Point_t > points;
+  read_points(filepoints, points);
+
+  // Compute the proximity graph of the points
+  start = std::chrono::system_clock::now();
+  Graph_t prox_graph = compute_proximity_graph(points, threshold
+                                               , euclidean_distance<Point_t>);
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "Compute Rips graph in " << elapsed_sec << " ms.\n";
+
+  // Construct the Rips complex in a Simplex Tree
+  Simplex_tree<Simplex_tree_options_fast_persistence> st;
+  start = std::chrono::system_clock::now();
+
+  // insert the proximity graph in the simplex tree
+  st.insert_graph(prox_graph);
+  // expand the graph until dimension dim_max
+  st.expansion(dim_max);
+
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "Compute Rips complex in " << elapsed_sec << " ms.\n";
+  std::cout << "  - dimension           = " << st.dimension() << std::endl;
+  std::cout << "  - number of simplices = " << st.num_simplices() << std::endl;
+
+  // Sort the simplices in the order of the filtration
+  start = std::chrono::system_clock::now();
+  st.initialize_filtration();
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "Order the simplices of the filtration in " << elapsed_sec << " ms.\n";
+
+  // Copy the keys inside the simplices
+  start = std::chrono::system_clock::now();
+  {
+    int count = 0;
+    for (auto sh : st.filtration_simplex_range())
+      st.assign_key(sh, count++);
+  }
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "Copied the keys inside the simplices in " << elapsed_sec << " ms.\n";
+
+  // Convert the simplex tree into a hasse diagram
+  start = std::chrono::system_clock::now();
+  Hasse_complex<> hcpx(st);
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "Convert the simplex tree into a Hasse diagram in " << elapsed_sec << " ms.\n";
+
+
+  std::cout << "Timings when using a simplex tree: \n";
+  timing_persistence(st, p);
+  timing_persistence(st, q);
+  timing_persistence(st, p, q);
+
+  std::cout << "Timings when using a Hasse complex: \n";
+  timing_persistence(hcpx, p);
+  timing_persistence(hcpx, q);
+  timing_persistence(hcpx, p, q);
+
+  start = std::chrono::system_clock::now();
+  }
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "Running the complex destructors in " << elapsed_sec << " ms.\n";
+  return 0;
+}
+
+template< typename FilteredComplex>
+void
+timing_persistence(FilteredComplex & cpx
+                   , int p) {
+  std::chrono::time_point<std::chrono::system_clock> start, end;
+  int elapsed_sec;
+  {
+  start = std::chrono::system_clock::now();
+  Persistent_cohomology< FilteredComplex, Field_Zp > pcoh(cpx);
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Initialize pcoh in " << elapsed_sec << " ms.\n";
+  // initializes the coefficient field for homology
+  start = std::chrono::system_clock::now();
+  pcoh.init_coefficients(p);
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Initialize the coefficient field in " << elapsed_sec << " ms.\n";
+
+  start = std::chrono::system_clock::now();
+
+  pcoh.compute_persistent_cohomology(INFINITY);
+
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Compute persistent homology in Z/" << p << "Z in " << elapsed_sec << " ms.\n";
+  start = std::chrono::system_clock::now();
+  }
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Run the persistence destructors in " << elapsed_sec << " ms.\n";
+}
+
+template< typename FilteredComplex>
+void
+timing_persistence(FilteredComplex & cpx
+                   , int p
+                   , int q) {
+  std::chrono::time_point<std::chrono::system_clock> start, end;
+  int elapsed_sec;
+  {
+  start = std::chrono::system_clock::now();
+  Persistent_cohomology< FilteredComplex, Multi_field > pcoh(cpx);
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Initialize pcoh in " << elapsed_sec << " ms.\n";
+  // initializes the coefficient field for homology
+  start = std::chrono::system_clock::now();
+  pcoh.init_coefficients(p, q);
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Initialize the coefficient field in " << elapsed_sec << " ms.\n";
+  // compute persistent homology, disgarding persistent features of life shorter than min_persistence
+
+  start = std::chrono::system_clock::now();
+
+  pcoh.compute_persistent_cohomology(INFINITY);
+
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Compute multi-field persistent homology in all coefficient fields Z/pZ "
+      << "with p in [" << p << ";" << q << "] in " << elapsed_sec << " ms.\n";
+  start = std::chrono::system_clock::now();
+  }
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Run the persistence destructors in " << elapsed_sec << " ms.\n";
+}