Improved the benching program a bit

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/relaxed-witness@1030 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: dc3d13fa785bee42a574134519b05e98b934d6f5

1 files changed, 204 insertions, 29 deletions

diff --git a/src/Witness_complex/example/relaxed_witness_persistence.cpp b/src/Witness_complex/example/relaxed_witness_persistence.cpp
index 09fe0d31..b8e4866f 100644
--- a/src/Witness_complex/example/relaxed_witness_persistence.cpp
+++ b/src/Witness_complex/example/relaxed_witness_persistence.cpp
@@ -50,41 +50,63 @@ using namespace Gudhi::persistent_cohomology;
 typedef std::vector<Point_d> Point_Vector;
 typedef Relaxed_witness_complex< Simplex_tree<> > RelaxedWitnessComplex;
 
-int main (int argc, char * const argv[])
-{
-  if (argc != 6) {
-      std::cerr << "Usage: " << argv[0]
-                << " nbP nbL dim alpha limD\n";
-      return 0;
+/**
+ * \brief Customized version of read_points
+ * which takes into account a possible nbP first line
+ *
+ */
+inline void
+read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) {
+  std::ifstream in_file(file_name.c_str(), std::ios::in);
+  if (!in_file.is_open()) {
+    std::cerr << "Unable to open file " << file_name << std::endl;
+    return;
+  }
+  std::string line;
+  double x;
+  while (getline(in_file, line)) {
+    std::vector< double > point;
+    std::istringstream iss(line);
+    while (iss >> x) {
+      point.push_back(x);
     }
-  /*
-  boost::filesystem::path p;
-  for (; argc > 2; --argc, ++argv)
-    p /= argv[1];
-  */
-  
-  int nbP       = atoi(argv[1]);
-  int nbL       = atoi(argv[2]);
-  int dim       = atoi(argv[3]);
-  double alpha  = atof(argv[4]);
-  int limD      = atoi(argv[5]);
-  //Construct the Simplex Tree
-  clock_t start, end;
+    if (point.size() != 1)
+      points.push_back(point);
+  }
+  in_file.close();
+}
 
-  // Construct the Simplex Tree
-  Simplex_tree<> simplex_tree;
+void output_experiment_information(char * const file_name)
+{
+    std::cout << "Enter a valid experiment number. Usage: "
+              << file_name << " exp_no options\n";
+    std::cout << "Experiment description:\n"
+              << "0 nbP nbL dim alpha limD: "
+              << "Build persistence diagram on relaxed witness complex "
+              << "built from a point cloud on (dim-1)-dimensional sphere "
+              << "consisting of nbP witnesses and nbL landmarks. "
+              << "The maximal relaxation is alpha and the limit on simplicial complex "
+              << "dimension is limD.\n";
+    std::cout << "1 file_name nbL alpha limD: "
+              << "Build persistence diagram on relaxed witness complex "
+              << "build from a point cloud stored in a file and nbL landmarks. "
+              << "The maximal relaxation is alpha and the limit on simplicial complex dimension is limD\n";
+}
 
-  // Read the point file
-  Point_Vector point_vector;
-  generate_points_sphere(point_vector, nbP, dim);
-  std::cout << "Successfully generated " << point_vector.size() << " points.\n";
-  std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n";
+void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha, int limD)
+{
+  clock_t start, end;
+  Simplex_tree<> simplex_tree;
 
   // Choose landmarks
   std::vector<std::vector< int > > knn;
   std::vector<std::vector< FT > > distances;
+  start = clock();
   Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha, limD, knn, distances);
-  
+  end = clock();
+  double time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
+  std::cout << "Choice of " << nbL << " landmarks took "
+            << time << " s. \n";
   // Compute witness complex
   start = clock();
   RelaxedWitnessComplex(distances,
@@ -94,15 +116,168 @@ int main (int argc, char * const argv[])
                         alpha,
                         limD);
   end = clock();
-  double time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
+  time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
   std::cout << "Witness complex for " << nbL << " landmarks took "
             << time << " s. \n";
+  std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n";
   // std::cout << simplex_tree << "\n";
   
   // Compute the persistence diagram of the complex
   persistent_cohomology::Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(simplex_tree, false);
   int p = 3;
   pcoh.init_coefficients( p ); //initilizes the coefficient field for homology
-  pcoh.compute_persistent_cohomology( alpha/10 );
+  start = clock();
+  pcoh.compute_persistent_cohomology( alpha/5 );
+  end = clock();
+  time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
+  std::cout << "Persistence diagram took "
+            << time << " s. \n";
   pcoh.output_diagram();
 }
+
+void rips_experiment(Point_Vector & points, double threshold, int dim_max)
+{
+  typedef std::vector<double> Point_t;
+  typedef Simplex_tree<Simplex_tree_options_fast_persistence> ST;
+  clock_t start, end;
+  ST st;
+
+  // Compute the proximity graph of the points
+  start = clock();
+  Graph_t prox_graph = compute_proximity_graph(points, threshold
+                                               , euclidean_distance<Point_t>);
+  // Construct the Rips complex in a Simplex Tree
+  // 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 = clock();
+  
+  double time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
+  std::cout << "Rips complex took "
+            << time << " s. \n";
+  std::cout << "The complex contains " << st.num_simplices() << " simplices \n";
+  //std::cout << "   and has dimension " << st.dimension() << " \n";
+
+  // Sort the simplices in the order of the filtration
+  st.initialize_filtration();
+
+  // Compute the persistence diagram of the complex
+  persistent_cohomology::Persistent_cohomology<ST, Field_Zp > pcoh(st);
+  // initializes the coefficient field for homology
+  int p = 3;
+  double min_persistence = threshold/5;
+  pcoh.init_coefficients(p);
+  pcoh.compute_persistent_cohomology(min_persistence);
+  pcoh.output_diagram();
+}
+
+int experiment0 (int argc, char * const argv[])
+{
+  if (argc != 7) {
+    std::cerr << "Usage: " << argv[0]
+              << " 0 nbP nbL dim alpha limD\n";
+    return 0;
+  }
+  /*
+    boost::filesystem::path p;
+    for (; argc > 2; --argc, ++argv)
+    p /= argv[1];
+  */
+  
+  int nbP       = atoi(argv[2]);
+  int nbL       = atoi(argv[3]);
+  int dim       = atoi(argv[4]);
+  double alpha  = atof(argv[5]);
+  int limD      = atoi(argv[6]);
+
+  // Read the point file
+  Point_Vector point_vector;
+  generate_points_sphere(point_vector, nbP, dim);
+  std::cout << "Successfully generated " << point_vector.size() << " points.\n";
+  std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n";
+
+  rw_experiment(point_vector, nbL, alpha, limD);
+}
+
+int experiment1 (int argc, char * const argv[])
+{
+  if (argc != 3) {
+    std::cerr << "Usage: " << argv[0]
+              << " 1 file_name\n";
+    return 0;
+  }
+  /*
+    boost::filesystem::path p;
+    for (; argc > 2; --argc, ++argv)
+    p /= argv[1];
+  */
+  
+  std::string file_name = argv[2];
+
+  // Read the point file
+  Point_Vector point_vector;
+  read_points_cust(file_name, point_vector);
+  std::cout << "The file contains " << point_vector.size() << " points.\n";
+  std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n";
+
+  bool ok = false;
+  int nbL, limD;
+  double alpha;
+  while (!ok) {
+    std::cout << "Relaxed witness complex: parameters nbL, alpha, limD.\n";
+    std::cout << "Enter nbL: ";
+    std::cin >> nbL;
+    std::cout << "Enter alpha: ";
+    std::cin >> alpha;
+    std::cout << "Enter limD: ";
+    std::cin >> limD;
+    std::cout << "Start relaxed witness complex...\n";
+    rw_experiment(point_vector, nbL, alpha, limD);
+    std::cout << "Is the result correct? [y/n]: ";
+    char answer;
+    std::cin >> answer;
+    switch (answer) {
+    case 'n':
+      ok = false; break;
+    default :
+      ok = true; break;
+    }
+  }
+  ok = false;
+  while (!ok) {
+    std::cout << "Rips complex: parameters threshold, limD.\n";
+    std::cout << "Enter threshold: ";
+    std::cin >> alpha;
+    std::cout << "Enter limD: ";
+    std::cin >> limD;
+    std::cout << "Start Rips complex...\n";
+    rips_experiment(point_vector, alpha, limD);
+    std::cout << "Is the result correct? [y/n]: ";
+    char answer;
+    std::cin >> answer;
+    switch (answer) {
+    case 'n':
+      ok = false; break;
+    default :
+      ok = true; break;
+    }
+  }
+}
+
+int main (int argc, char * const argv[])
+{
+  if (argc == 1) {
+    output_experiment_information(argv[0]);
+    return 1;
+  }
+  switch (atoi(argv[1])) {
+  case 0 :
+    return experiment0(argc, argv);
+  case 1 :
+    return experiment1(argc, argv);
+  default :
+    output_experiment_information(argv[0]);
+    return 1;
+  }
+}