All work up until now

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

Former-commit-id: a9379060f1edfbc419e259b2e207ebc608798803

7 files changed, 244 insertions, 129 deletions

diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt
index 80396392..d6a958b7 100644
--- a/src/Witness_complex/example/CMakeLists.txt
+++ b/src/Witness_complex/example/CMakeLists.txt
@@ -42,10 +42,12 @@ if(CGAL_FOUND)
       target_link_libraries(a0_persistence ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY})
       add_test( a0_persistence_10 ${CMAKE_CURRENT_BINARY_DIR}/a0_persistence 10)
       add_executable ( bench_rwit bench_rwit.cpp )
+      target_link_libraries(bench_rwit ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
                   
       add_executable ( relaxed_delaunay relaxed_delaunay.cpp )
       target_link_libraries(relaxed_delaunay ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY})
       add_test( relaxed_delaunay_10 ${CMAKE_CURRENT_BINARY_DIR}/relaxed_delaunay 10)   
+      add_executable ( generate_torus generate_torus.cpp )
     else()
       message(WARNING "Eigen3 not found. Version 3.1.0 is required for witness_complex_sphere example.")      
     endif()
diff --git a/src/Witness_complex/example/Landmark_choice_random_knn.h b/src/Witness_complex/example/Landmark_choice_random_knn.h
index 3797b4c5..fb91e116 100644
--- a/src/Witness_complex/example/Landmark_choice_random_knn.h
+++ b/src/Witness_complex/example/Landmark_choice_random_knn.h
@@ -131,7 +131,7 @@ typedef Neighbor_search::Tree Tree;
           distances[points_i].push_back(sqrt(search_it->second));
           knn[points_i].push_back((search_it++)->first);
         }
-      std::cout << "k = " << knn[points_i].size() << std::endl;
+      //std::cout << "k = " << knn[points_i].size() << std::endl;
     }
   }
 
diff --git a/src/Witness_complex/example/Landmark_choice_sparsification.h b/src/Witness_complex/example/Landmark_choice_sparsification.h
index 1052b0c4..d5d5fba1 100644
--- a/src/Witness_complex/example/Landmark_choice_sparsification.h
+++ b/src/Witness_complex/example/Landmark_choice_sparsification.h
@@ -74,7 +74,7 @@ namespace witness_complex {
                                          FT mu_epsilon,
                                          Point_random_access_range & landmarks)
   {
-    int nbP = points.end() - points.begin();
+    unsigned nbP = points.end() - points.begin();
     assert(nbP >= nbL);
     CGAL::Random rand;
     // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety
@@ -87,7 +87,7 @@ namespace witness_complex {
                       typename Tree::Splitter(),
                       points_traits);
     
-    for (int points_i = 0; points_i < nbP; points_i++) {
+    for (unsigned points_i = 0; points_i < nbP; points_i++) {
       if (dropped_points[points_i])
         continue;
       Point_d & w = points[points_i];
@@ -98,7 +98,7 @@ namespace witness_complex {
         dropped_points[i] = true;
     }
 
-    for (int points_i = 0; points_i < nbP; points_i++) {
+    for (unsigned points_i = 0; points_i < nbP; points_i++) {
       if (dropped_points[points_i])
         landmarks.push_back(points[points_i]);
     }
diff --git a/src/Witness_complex/example/a0_persistence.cpp b/src/Witness_complex/example/a0_persistence.cpp
index 295c7115..422185ca 100644
--- a/src/Witness_complex/example/a0_persistence.cpp
+++ b/src/Witness_complex/example/a0_persistence.cpp
@@ -138,6 +138,12 @@ void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha2, int limD, FT
   std::cout << "Witness complex for " << nbL << " landmarks took "
             << time << " s. \n";
   std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n";
+  int streedim = 0;
+  for (auto s: simplex_tree.complex_simplex_range())
+    if (simplex_tree.dimension(s) > streedim)
+      streedim = simplex_tree.dimension(s);
+  std::cout << "Dimension of the simplicial complex is " << streedim << std::endl;
+
   //std::cout << simplex_tree << "\n";
   
   // Compute the persistence diagram of the complex
@@ -348,7 +354,7 @@ int experiment1 (int argc, char * const argv[])
   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";
-
+  
   Simplex_tree<> simplex_tree;
   std::vector<std::vector< int > > knn;
   std::vector<std::vector< FT > > distances;
diff --git a/src/Witness_complex/example/bench_rwit.cpp b/src/Witness_complex/example/bench_rwit.cpp
index f83f9f42..2d3a009c 100644
--- a/src/Witness_complex/example/bench_rwit.cpp
+++ b/src/Witness_complex/example/bench_rwit.cpp
@@ -26,6 +26,7 @@
 #include <gudhi/Dim_lists.h>
 #include <gudhi/reader_utils.h>
 #include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Good_links.h>
 #include "Landmark_choice_random_knn.h"
 #include "Landmark_choice_sparsification.h"
 
@@ -44,6 +45,8 @@
 #include <boost/iterator/counting_iterator.hpp>
 #include <boost/range/iterator_range.hpp>
 
+#include <boost/program_options.hpp>
+
 #include "generators.h"
 #include "output.h"
 #include "output_tikz.h"
@@ -53,14 +56,96 @@ using namespace Gudhi::witness_complex;
 using namespace Gudhi::persistent_cohomology;
 
 typedef std::vector<Point_d> Point_Vector;
-typedef Simplex_tree<Simplex_tree_options_fast_persistence> STree;
-//typedef Simplex_tree<> STree;
+//typedef Simplex_tree<Simplex_tree_options_fast_persistence> STree;
+typedef Simplex_tree<> STree;
 typedef A0_complex< STree> A0Complex;
 typedef STree::Simplex_handle Simplex_handle;
 
 typedef A0_complex< STree > SRWit;
 typedef Relaxed_witness_complex< STree > WRWit;
 
+/** Program options ***************************************************************
+***********************************************************************************
+***********************************************************************************
+***********************************************************************************
+**********************************************************************************/
+
+
+void program_options(int argc, char * const argv[]
+                     , int & experiment_number
+                     , std::string & filepoints
+                     , std::string & landmark_file
+                     , std::string & experiment_name
+                     , int & nbL
+                     , double & alpha2_s
+                     , double & alpha2_w
+                     , double & mu_epsilon
+                     , int & dim_max
+                     , std::vector<int> & desired_homology
+                     , double & min_persistence) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()
+    ("option", po::value<int>(& experiment_number),
+     "Experiment id.")
+    ("input-file", po::value<std::string>(&filepoints),
+       "Name of file containing a point set. Format is one point per line:   X1 ... Xd ");
+  
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()
+      ("help,h", "produce help message")
+      ("output-file,o", po::value<std::string>(&experiment_name)->default_value("witness"),
+       "The prefix of all the output files. Default is 'witness'")
+      ("landmarks,L", po::value<int>(&nbL)->default_value(0),
+       "Number of landmarks.")
+      ( "landmark-file,l", po::value<std::string>(&landmark_file),
+        "Name of a fike containing landmarks")
+      ("alpha2_s,A", po::value<double>(&alpha2_s)->default_value(0),
+       "Relaxation parameter for the strong complex.")
+      ("alpha2_w,a", po::value<double>(&alpha2_w)->default_value(0),
+       "Relaxation parameter for the weak complex.")
+      ("mu_epsilon,e", po::value<double>(&mu_epsilon)->default_value(0),
+       "Sparsification parameter.")
+      ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
+       "Maximal dimension of the Witness complex we want to compute.")
+    
+    ("homology,H", po::value<std::vector<int>>(&desired_homology)->multitoken(),
+    "The desired Betti numbers.")
+      ("min-persistence,m", po::value<Filtration_value>(&min_persistence),
+       "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals");
+
+  po::positional_options_description pos;
+  pos.add("option", 1);
+  pos.add("input-file", 2);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).
+            options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/3Z \n";
+    std::cout << "of a Strong relaxed witness complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
+
+
+
+
+
 /**
  * \brief Customized version of read_points
  * which takes into account a possible nbP first line
@@ -175,6 +260,7 @@ void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha2, int limD, FT
   Gudhi::witness_complex::Dim_lists<STree> simplices(simplex_tree, limD); 
   
   // std::vector<Simplex_handle> simplices;
+  std::cout << "Starting collapses...\n";
   simplices.collapse();
   simplices.output_simplices();
 
@@ -204,6 +290,21 @@ void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha2, int limD, FT
     chi += 1-2*(simplex_tree.dimension(sh)%2);
   std::cout << "Euler characteristic is " << chi << std::endl;
   write_witness_mesh(point_vector, landmarks_ind, collapsed_tree, collapsed_tree.complex_simplex_range(), false, true, mesh_filename+"_after_collapse.mesh");
+  Gudhi::Good_links<STree> gl(collapsed_tree);
+  if (gl.complex_is_pseudomanifold())
+    std::cout << "Collapsed complex is a pseudomanifold.\n";
+  else
+    std::cout << "Collapsed complex is NOT a pseudomanifold.\n";
+  bool good = true;
+  for (auto v: collapsed_tree.complex_vertex_range())
+    if (!gl.has_good_link(v)) {
+      std::cout << "Bad link around " << v << std::endl;
+      good = false;
+    }
+  if (good)
+    std::cout << "All links are good.\n";
+  else
+    std::cout << "There are bad links.\n";
 }
 
 void rips_experiment(Point_Vector & points, double threshold, int dim_max)
@@ -267,83 +368,6 @@ int experiment0 (int argc, char * const argv[])
   return 0;
 }
 
-int experiment1 (int argc, char * const argv[])
-{
-  if (argc != 8) {
-    std::cerr << "Usage: " << argv[0]
-              << " 1 file_name nbL alpha mu_epsilon limD experiment_name\n";
-    return 0;
-  }
-  /*
-    boost::filesystem::path p;
-    for (; argc > 2; --argc, ++argv)
-    p /= argv[1];
-  */
-  
-  std::string file_name = argv[2];
-  int nbL = atoi(argv[3]), limD = atoi(argv[6]);
-  double alpha2 = atof(argv[4]), mu_epsilon = atof(argv[5]);
-  std::string experiment_name = argv[7];
-
-  // 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";
-
-  STree simplex_tree;
-  std::vector<std::vector< int > > knn;
-  std::vector<std::vector< FT > > distances;
-  std::vector<Point_d> landmarks;
-  Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
-  Gudhi::witness_complex::build_distance_matrix(point_vector,   // aka witnesses
-                                                landmarks,  // aka landmarks
-                                                alpha2,
-                                                limD,
-                                                knn,
-                                                distances);
-  WRWit      rw(distances,
-                knn,
-                simplex_tree,
-                nbL,
-                0,
-                limD);
-  std::vector<int> landmarks_ind(nbL); 
-  for (unsigned i = 0; i != distances.size(); ++i) {
-    if (distances[i][0] == 0)
-      landmarks_ind[knn[i][0]] = i;
-  }
-
-  write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplex_tree.complex_simplex_range(), false, true, experiment_name+"_witness_complex.mesh");
-
-  int chi = 0;
-  for (auto sh: simplex_tree.complex_simplex_range())
-    chi += 1-2*(simplex_tree.dimension(sh)%2);
-  std::cout << "Euler characteristic is " << chi << std::endl;
-
-  rw_experiment(point_vector, nbL, alpha2, limD, mu_epsilon, experiment_name);
-  
-  // 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;
-  //   }
-  // }
-  return 0;
-}
 
 /********************************************************************************************
  * Length of the good interval experiment
@@ -380,17 +404,25 @@ double good_interval_length(const std::vector<int> & desired_homology, STree & s
   std::string line;
   std::vector<Pers_endpoint> pers_endpoints;
   while (getline(in_stream, line)) {
-    int p, dim;
+    unsigned p, dim;
     double alpha_start, alpha_end;
     std::istringstream iss(line);
     iss >> p >> dim >> alpha_start >> alpha_end;
+    if (iss.fail())
+      alpha_end = alpha2;
+    //std::cout << p << " " << dim << " " << alpha_start << " " << alpha_end << "\n";
+    //if (dim < desired_homology.size()+1)
     if (alpha_start != alpha_end) {
-      if (alpha_end < alpha_start)
-        alpha_end = alpha2;
+    //   if (alpha_end < alpha_start)
+    //     alpha_end = alpha2;
       pers_endpoints.push_back(Pers_endpoint(alpha_start, true, dim));
       pers_endpoints.push_back(Pers_endpoint(alpha_end, false, dim));
+      std::cout << p << " " << dim << " " << alpha_start << " " << alpha_end << "\n";
     }
   }
+  std::cout << "desired_homology.size() = " << desired_homology.size() << "\n";
+  for (auto nd: desired_homology)
+    std::cout << nd << std::endl;
   std::cout << "Pers_endpoints.size = " << pers_endpoints.size() << std::endl;
   in_stream.close();
   std::sort(pers_endpoints.begin(),
@@ -404,14 +436,21 @@ double good_interval_length(const std::vector<int> & desired_homology, STree & s
     std::cout << p.alpha << " " << p.dim << " " << p.start << "\n";
   }
   */
-  std::vector<int> current_homology(nbL-1,0);
-  current_homology[0] = 1; // for the compulsary "0 0 inf" entry
+  std::vector<int> current_homology(desired_homology.size(),0);
+  //current_homology[0] = 1; // for the compulsary "0 0 inf" entry
   double good_start = 0, good_end = 0;
   double sum_intervals = 0;
   int num_pieces = 0;
   bool interval_in_process = (desired_homology == current_homology);
   for (auto p: pers_endpoints) {
     /*
+      std::cout << p.alpha << " " << p.dim << " ";
+    if (p.start)
+      std::cout << "s\n";
+    else
+      std::cout << "e\n";
+    */
+    /*
     std::cout << "Treating " << p.alpha << " " << p.dim << " " << p.start
               << " [";
     for (int v: current_homology)
@@ -450,66 +489,114 @@ void run_comparison(std::vector<std::vector< int > > const & knn,
                     std::vector<std::vector< FT > > const & distances,
                     Point_Vector & points,
                     unsigned nbL,
-                    double alpha2,
+                    unsigned limD,
+                    double alpha2_s,
+                    double alpha2_w,
                     std::vector<int>& desired_homology)
 {
   clock_t start, end;
   STree simplex_tree;
 
-  alpha2 = 0.55;
-  std::cout << "alpha2 = " << alpha2 << "\n";
+  //std::cout << "alpha2 = " << alpha2_s << "\n";
   start = clock();
   SRWit srwit(distances,
               knn,
               simplex_tree,
               nbL,
-              alpha2,
-              nbL-1);
+              alpha2_s,
+              limD);
   end = clock();
   std::cout << "SRWit.size = " << simplex_tree.num_simplices() << std::endl;
-  simplex_tree.set_dimension(nbL-1);
+  simplex_tree.set_dimension(desired_homology.size());
 
   std::cout << "Good homology interval length for SRWit is "
-                << good_interval_length(desired_homology, simplex_tree, alpha2) << "\n";
+                << good_interval_length(desired_homology, simplex_tree, alpha2_s) << "\n";
   std::cout << "Time: " << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n";
+  int chi = 0;
+  for (auto sh: simplex_tree.complex_simplex_range())
+    chi += 1-2*(simplex_tree.dimension(sh)%2);
+  std::cout << "Euler characteristic is " << chi << std::endl;
+
 
   STree simplex_tree2;
-  alpha2 = 0.35;
-  std::cout << "alpha2 = " << alpha2 << "\n";
+  std::cout << "alpha2 = " << alpha2_w << "\n";
   start = clock();
   WRWit wrwit(distances,
               knn,
               simplex_tree2,
               nbL,
-              alpha2,
-              nbL-1);
+              alpha2_w,
+              limD);
   end = clock();
   std::cout << "WRWit.size = " << simplex_tree2.num_simplices() << std::endl;
-  simplex_tree.set_dimension(nbL-1);
+  simplex_tree2.set_dimension(nbL-1);
   
   std::cout << "Good homology interval length for WRWit is "
-            << good_interval_length(desired_homology, simplex_tree2, alpha2) << "\n";
+            << good_interval_length(desired_homology, simplex_tree2, alpha2_w) << "\n";
   std::cout << "Time: " << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n";
+  chi = 0;
+  for (auto sh: simplex_tree2.complex_simplex_range())
+    chi += 1-2*(simplex_tree2.dimension(sh)%2);
+  std::cout << "Euler characteristic is " << chi << std::endl;
+
+  //write_witness_mesh(points, landmarks_ind, simplex_tree2, simplex_tree2.complex_simplex_range(), false, true, "wrwit.mesh");
 
-  STree simplex_tree3;
-  std::cout << "Enter alpha2 for Rips: ";
-  std::cin >> alpha2;
-  std::cout << "\n";
-  std::cout << "points.size() = " << points.size() << "\n";
-  start = clock();
-  rips(points,
-       0.2, 
-       nbL-1,
-       simplex_tree3);
-  end = clock();
-  std::cout << "Rips.size = " << simplex_tree3.num_simplices() << std::endl;
-  simplex_tree.set_dimension(nbL-1);
   
-  std::cout << "Good homology interval length for WRWit is "
-            << good_interval_length(desired_homology, simplex_tree3, alpha2) << "\n";
-  std::cout << "Time: " << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n";
 }
 
+int experiment1 (int argc, char * const argv[])
+{
+  /*
+    boost::filesystem::path p;
+    for (; argc > 2; --argc, ++argv)
+    p /= argv[1];
+  */
+  
+  // std::string file_name = argv[2];
+  // int nbL = atoi(argv[3]), limD = atoi(argv[6]);
+  // double alpha2 = atof(argv[4]), mu_epsilon = atof(argv[5]);
+  // std::string experiment_name = argv[7];
+
+  int option = 1;
+  std::string file_name, landmark_file;
+  int nbL = 0, limD;
+  double alpha2_s, alpha2_w, mu_epsilon, min_pers;
+  std::string experiment_name;
+  std::vector<int> desired_homology = {1};
+  std::vector<Point_d> landmarks;
+  
+  program_options(argc, argv, option, file_name, landmark_file, experiment_name, nbL, alpha2_s, alpha2_w, mu_epsilon, limD, desired_homology, min_pers);
+
+  // 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";
+  //std::cout << "Limit dimension for the complexes is " << limD << ".\n";
+
+  if (landmark_file == "")
+    Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
+    //Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha2_s, limD, knn, distances);
+  else
+    read_points_cust(landmark_file, landmarks);
+  nbL = landmarks.size();
+  STree simplex_tree;
+  std::vector<std::vector< int > > knn;
+  std::vector<std::vector< FT > > distances;
+
+  //Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
+  Gudhi::witness_complex::build_distance_matrix(point_vector,   // aka witnesses
+                                                landmarks,  // aka landmarks
+                                                alpha2_s,
+                                                limD,
+                                                knn,
+                                                distances);
+
+  run_comparison(knn, distances, point_vector, nbL, limD, alpha2_s, alpha2_w, desired_homology);
+  return 0;
+}
+
+
 int experiment2(int argc, char * const argv[])
 {
   for (unsigned d = 3; d < 4; d++) {
@@ -524,7 +611,6 @@ int experiment2(int argc, char * const argv[])
     case 4: alpha2 = 1.4; break;
     default: alpha2 = 1.4; break;
     }
-    std::cout << "alpha2 = " << alpha2 << "\n";
     unsigned nbL = 20;
     std::vector<int> desired_homology(nbL-1,0);
     desired_homology[0] = 1; desired_homology[d] = 1;
@@ -536,19 +622,22 @@ int experiment2(int argc, char * const argv[])
       std::cout << "Running test S"<< d <<", |W|=" << nbW << ", |L|=" << nbL << std::endl;
       generate_points_sphere(point_vector, i*1000, d+1);
       std::vector<Point_d> landmarks;
-      Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
 
+      Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
+      
       std::vector<std::vector< int > > knn;
       std::vector<std::vector< FT > > distances;
-    
+
+      
       std::cout << "|L| after sparsification: " << landmarks.size() << "\n";
+      
       Gudhi::witness_complex::build_distance_matrix(point_vector,   // aka witnesses
                                                     landmarks,  // aka landmarks
                                                     alpha2,
                                                     nbL-1,
                                                     knn,
                                                     distances);
-      run_comparison(knn, distances, point_vector, nbL, alpha2, desired_homology);
+      run_comparison(knn, distances, point_vector, nbL, nbL-1, alpha2, alpha2, desired_homology);
     }
   }
   /*
@@ -588,7 +677,8 @@ int experiment2(int argc, char * const argv[])
 
 int experiment3(int argc, char * const argv[])
 {
-  // COLLAPSES EXPERIMENT
+  // Both witnesses and landmarks are given as input
+  
   
   return 0;
 }
@@ -609,6 +699,9 @@ int main (int argc, char * const argv[])
   case 2 :
     return experiment2(argc, argv);
     break;
+  case 3 :
+    return experiment3(argc, argv);
+    break;
   default :
     output_experiment_information(argv[0]);
     return 1;
diff --git a/src/Witness_complex/example/generators.h b/src/Witness_complex/example/generators.h
index 23915546..731a52b0 100644
--- a/src/Witness_complex/example/generators.h
+++ b/src/Witness_complex/example/generators.h
@@ -25,10 +25,12 @@
 
 #include <CGAL/Epick_d.h>
 #include <CGAL/point_generators_d.h>
+#include <CGAL/Random.h>
 
 #include <fstream>
 #include <string>
 #include <vector>
+#include <cmath>
 
 typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K;
 typedef K::FT FT;
@@ -144,12 +146,21 @@ void generate_points_sphere(Point_Vector& W, int nbP, int dim) {
     W.push_back(*rp++);
 }
 
-/** \brief Generate nbP points on a d-torus
+/** \brief Generate nbP points on a (flat) d-torus embedded in R^{2d}
  *
  */
 void generate_points_torus(Point_Vector& W, int nbP, int dim) {
-  CGAL::Random_points_on_sphere_d<Point_d> rp(2, 1);
-  
+  CGAL::Random rand;
+  const double pi = std::acos(-1);
+  for (int i = 0; i < nbP; i++) {
+    std::vector<FT> point;
+    for (int j = 0; j < dim; j++) {
+      double alpha = rand.uniform_real((double)0, 2*pi);
+      point.push_back(sin(alpha));
+      point.push_back(cos(alpha));
+    }
+    W.push_back(Point_d(point));
+  }
 }
 
 #endif  // EXAMPLE_WITNESS_COMPLEX_GENERATORS_H_
diff --git a/src/Witness_complex/example/output.h b/src/Witness_complex/example/output.h
index 0e74387a..d3f534af 100644
--- a/src/Witness_complex/example/output.h
+++ b/src/Witness_complex/example/output.h
@@ -73,9 +73,12 @@ void write_edges(std::string file_name, STree& witness_complex, Point_Vector& la
   ofs.close();
 }
 
-/** \brief Write triangles (tetrahedra in 3d) of a witness
- *  complex in a file, compatible with medit.
- *  l_is_v = landmark is vertex
+/** \brief Write triangles (tetrahedra in 3d) of a simplicial complex in a file, compatible with medit.
+ *  `landmarks_ind` represents the set of landmark indices in W  
+ *  `st` is the Simplex_tree to be visualized,
+ *  `shr` is the Simplex_handle_range of simplices in `st` to be visualized
+ *  `is2d` should be true if the simplicial complex is 2d, false if 3d
+ *  `l_is_v` = landmark is vertex
  */
 template <typename SimplexHandleRange,
           typename STree >