Added the possibility to limit dimension

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

Former-commit-id: 16ec8706ffb981ff892960468cc4b809e61612e2

7 files changed, 57 insertions, 27 deletions

diff --git a/src/Witness_complex/example/example_strong_witness_persistence.cpp b/src/Witness_complex/example/example_strong_witness_persistence.cpp
index bf6b4011..4e06867d 100644
--- a/src/Witness_complex/example/example_strong_witness_persistence.cpp
+++ b/src/Witness_complex/example/example_strong_witness_persistence.cpp
@@ -53,17 +53,18 @@ void program_options(int argc, char * argv[]
                      , std::string & filediag
                      , Filtration_value & max_squared_alpha
                      , int & p
+                     , int & dim_max
                      , Filtration_value & min_persistence);
 
 int main(int argc, char * argv[]) {
   std::string file_name;
   std::string filediag;
   Filtration_value max_squared_alpha;
-  int p, nbL;
+  int p, nbL, lim_d;
   Filtration_value min_persistence;
   SimplexTree simplex_tree;
 
-  program_options(argc, argv, nbL, file_name, filediag, max_squared_alpha, p, min_persistence);
+  program_options(argc, argv, nbL, file_name, filediag, max_squared_alpha, p, lim_d, min_persistence);
 
   // Extract the points from the file file_name
   Point_vector witnesses, landmarks;
@@ -119,6 +120,7 @@ void program_options(int argc, char * argv[]
                      , std::string & filediag
                      , Filtration_value & max_squared_alpha
                      , int & p
+                     , int & dim_max
                      , Filtration_value & min_persistence) {
   
   namespace po = boost::program_options;
@@ -141,7 +143,9 @@ void program_options(int argc, char * argv[]
       ("field-charac,p", po::value<int>(&p)->default_value(11),
        "Characteristic p of the coefficient field Z/pZ for computing homology.")
       ("min-persistence,m", po::value<Filtration_value>(&min_persistence)->default_value(0),
-       "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals");
+       "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals")
+      ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(std::numeric_limits<int>::max()),
+       "Maximal dimension of the strong witness complex we want to compute.");
  
   po::positional_options_description pos;
   pos.add("input-file", 1);
diff --git a/src/Witness_complex/example/example_witness_complex_off.cpp b/src/Witness_complex/example/example_witness_complex_off.cpp
index 849a0a40..37646faf 100644
--- a/src/Witness_complex/example/example_witness_complex_off.cpp
+++ b/src/Witness_complex/example/example_witness_complex_off.cpp
@@ -45,12 +45,12 @@ typedef std::vector< Point_d > Point_vector;
 int main(int argc, char * const argv[]) {
   if (argc != 4) {
     std::cerr << "Usage: " << argv[0]
-        << " path_to_point_file nbL alpha^2\n";
+        << " path_to_point_file number_of_landmarks max_squared_alpha limit_dimension\n";
     return 0;
   }
 
   std::string file_name = argv[1];
-  int nbL = atoi(argv[2]);
+  int nbL = atoi(argv[2]), lim_dim = atoi(argv[4]);
   double alpha2 = atof(argv[3]);
   clock_t start, end;
   Gudhi::Simplex_tree<> simplex_tree;
@@ -77,7 +77,7 @@ int main(int argc, char * const argv[]) {
                                   point_vector.begin(),
                                   point_vector.end());
 
-  witness_complex.create_complex(simplex_tree, alpha2);
+  witness_complex.create_complex(simplex_tree, alpha2, lim_dim);
   end = clock();
   std::cout << "Witness complex took "
       << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n";
diff --git a/src/Witness_complex/example/example_witness_complex_persistence.cpp b/src/Witness_complex/example/example_witness_complex_persistence.cpp
index 9b06b504..a4388047 100644
--- a/src/Witness_complex/example/example_witness_complex_persistence.cpp
+++ b/src/Witness_complex/example/example_witness_complex_persistence.cpp
@@ -53,17 +53,18 @@ void program_options(int argc, char * argv[]
                      , std::string & filediag
                      , Filtration_value & max_squared_alpha
                      , int & p
+                     , int & dim_max
                      , Filtration_value & min_persistence);
 
 int main(int argc, char * argv[]) {
   std::string file_name;
   std::string filediag;
   Filtration_value max_squared_alpha;
-  int p, nbL;
+  int p, nbL, lim_d;
   Filtration_value min_persistence;
   SimplexTree simplex_tree;
 
-  program_options(argc, argv, nbL, file_name, filediag, max_squared_alpha, p, min_persistence);
+  program_options(argc, argv, nbL, file_name, filediag, max_squared_alpha, p, lim_d, min_persistence);
 
   // Extract the points from the file file_name
   Point_vector witnesses, landmarks;
@@ -119,6 +120,7 @@ void program_options(int argc, char * argv[]
                      , std::string & filediag
                      , Filtration_value & max_squared_alpha
                      , int & p
+                     , int & dim_max
                      , Filtration_value & min_persistence) {
   
   namespace po = boost::program_options;
@@ -141,8 +143,10 @@ void program_options(int argc, char * argv[]
       ("field-charac,p", po::value<int>(&p)->default_value(11),
        "Characteristic p of the coefficient field Z/pZ for computing homology.")
       ("min-persistence,m", po::value<Filtration_value>(&min_persistence)->default_value(0),
-       "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals");
- 
+       "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals")
+      ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(std::numeric_limits<int>::max()),
+       "Maximal dimension of the weak witness complex we want to compute.");
+  
   po::positional_options_description pos;
   pos.add("input-file", 1);
 
diff --git a/src/Witness_complex/include/gudhi/Active_witness/Active_witness.h b/src/Witness_complex/include/gudhi/Active_witness/Active_witness.h
index 9ae41a69..e52410e4 100644
--- a/src/Witness_complex/include/gudhi/Active_witness/Active_witness.h
+++ b/src/Witness_complex/include/gudhi/Active_witness/Active_witness.h
@@ -31,6 +31,9 @@ namespace Gudhi {
 
 namespace witness_complex {
 
+  /** \brief Class representing a list of nearest neighbors to a given witness.
+   *  \detail Every element is a pair of a landmark identifier and the squared distance to it.
+   */
 template< typename Id_distance_pair,
           typename INS_range >
 class Active_witness {
diff --git a/src/Witness_complex/include/gudhi/Active_witness/Active_witness_iterator.h b/src/Witness_complex/include/gudhi/Active_witness/Active_witness_iterator.h
index 5d4f3d75..9c96f7e8 100644
--- a/src/Witness_complex/include/gudhi/Active_witness/Active_witness_iterator.h
+++ b/src/Witness_complex/include/gudhi/Active_witness/Active_witness_iterator.h
@@ -31,6 +31,13 @@ namespace Gudhi {
 
 namespace witness_complex {
 
+  /** \brief Iterator in the nearest landmark list.
+   *  \detail After the iterator reaches the end of the list,
+   *          the list is augmented by a (nearest landmark, distance) pair if possible.
+   *          If all the landmarks are present in the list, iterator returns the specific end value
+   *          of the corresponding 'Active_witness' object.
+   */
+  
 template< typename Active_witness,
           typename Id_distance_pair,
           typename INS_iterator >
diff --git a/src/Witness_complex/include/gudhi/Strong_witness_complex.h b/src/Witness_complex/include/gudhi/Strong_witness_complex.h
index 3a862d71..e125d307 100644
--- a/src/Witness_complex/include/gudhi/Strong_witness_complex.h
+++ b/src/Witness_complex/include/gudhi/Strong_witness_complex.h
@@ -138,19 +138,26 @@ private:
   /** \brief Outputs the (weak) witness complex with 
    *         squared relaxation parameter 'max_alpha_square'
    *         to simplicial complex 'complex'.
+   *         The parameter 'limit_dimension' represents the maximal dimension of the simplicial complex
+   *         (default value = no limit).
    */
   template < typename SimplicialComplexForWitness >
   bool create_complex(SimplicialComplexForWitness& complex,
-                      FT  max_alpha_square)       
+                      FT  max_alpha_square,
+                      Landmark_id limit_dimension = std::numeric_limits<Landmark_id>::max())       
   {
-    unsigned nbL = landmarks_.size();
-    unsigned complex_dim = 0;
+    std::size_t nbL = landmarks_.size();
+    Landmark_id complex_dim = 0;
     if (complex.num_vertices() > 0) {
-      std::cerr << "Witness complex cannot create complex - complex is not empty.\n";
+      std::cerr << "Strong witness complex cannot create complex - complex is not empty.\n";
       return false;
     }
     if (max_alpha_square < 0) {
-      std::cerr << "Witness complex cannot create complex - squared relaxation parameter must be non-negative.\n";
+      std::cerr << "Strong witness complex cannot create complex - squared relaxation parameter must be non-negative.\n";
+      return false;
+    }
+    if (limit_dimension < 0) {
+      std::cerr << "Strong witness complex cannot create complex - limit dimension must be non-negative.\n";
       return false;
     }
     typeVectorVertex vv;
@@ -165,13 +172,13 @@ private:
       ActiveWitness aw(landmark_tree_.query_incremental_nearest_neighbors(w));
       typeVectorVertex simplex;
       typename ActiveWitness::iterator aw_it = aw.begin();
-      float lim_d2 = aw.begin()->second + max_alpha_square;
-      while (aw_it != aw.end() && aw_it->second < lim_d2) {
+      float lim_dist2 = aw.begin()->second + max_alpha_square;
+      while ((Landmark_id)simplex.size() <= limit_dimension + 1 && aw_it != aw.end() && aw_it->second < lim_dist2) {
         simplex.push_back(aw_it->first);
         complex.insert_simplex_and_subfaces(simplex, aw_it->second - aw.begin()->second);
         aw_it++;
       }
-      if (simplex.size() - 1 > complex_dim)
+      if ((Landmark_id)simplex.size() - 1 > complex_dim)
         complex_dim = simplex.size() - 1;
     }
     complex.set_dimension(complex_dim);
diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h
index 7d8fce86..2a89306d 100644
--- a/src/Witness_complex/include/gudhi/Witness_complex.h
+++ b/src/Witness_complex/include/gudhi/Witness_complex.h
@@ -129,7 +129,8 @@ private:
    */
   template < typename SimplicialComplexForWitness >
   bool create_complex(SimplicialComplexForWitness& complex,
-                      FT  max_alpha_square)       
+                      FT  max_alpha_square,
+                      Landmark_id limit_dimension = std::numeric_limits<Landmark_id>::max())       
   {
     std::size_t nbL = landmarks_.size();
     if (complex.num_vertices() > 0) {
@@ -140,6 +141,10 @@ private:
       std::cerr << "Witness complex cannot create complex - squared relaxation parameter must be non-negative.\n";
       return false;
     }
+    if (limit_dimension < 0) {
+      std::cerr << "Witness complex cannot create complex - limit dimension must be non-negative.\n";
+      return false;
+    }
     typeVectorVertex vv;
     ActiveWitnessList active_witnesses;// = new ActiveWitnessList();
     for (unsigned i = 0; i != nbL; ++i) {
@@ -150,13 +155,13 @@ private:
       complex.insert_simplex(vv, Filtration_value(0.0));
       /* TODO Error if not inserted : normally no need here though*/
     }
-    unsigned k = 1; /* current dimension in iterative construction */
+    Landmark_id k = 1; /* current dimension in iterative construction */
     for (auto w: witnesses_)
       active_witnesses.push_back(ActiveWitness(landmark_tree_.query_incremental_nearest_neighbors(w)));
     ActiveWitness aw_copy(active_witnesses.front());
-    while (!active_witnesses.empty() && k < nbL ) {
+    while (!active_witnesses.empty() && k <= limit_dimension ) {
       typename ActiveWitnessList::iterator aw_it = active_witnesses.begin();
-      std::vector<int> simplex;
+      std::vector<Landmark_id> simplex;
       simplex.reserve(k+1);
       while (aw_it != active_witnesses.end()) {
         bool ok = add_all_faces_of_dimension(k,
@@ -191,7 +196,7 @@ private:
                                   double alpha2,
                                   double norelax_dist2,
                                   typename ActiveWitness::iterator curr_l,
-                                  std::vector<int>& simplex,
+                                  std::vector<Landmark_id>& simplex,
                                   SimplicialComplexForWitness& sc,
                                   typename ActiveWitness::iterator end)
   {
@@ -251,17 +256,17 @@ private:
    *  inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id
    */
   template < typename SimplicialComplexForWitness >
-  bool all_faces_in(std::vector<int>& simplex,
+  bool all_faces_in(typeVectorVertex& simplex,
                     double* filtration_value,
                     SimplicialComplexForWitness& sc)
   {
     typedef typename SimplicialComplexForWitness::Simplex_handle Simplex_handle;
 
-    std::vector< int > facet;
-    for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it)
+    typeVectorVertex facet;
+    for (typename typeVectorVertex::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it)
       {
         facet.clear();
-        for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it)
+        for (typename typeVectorVertex::iterator it = simplex.begin(); it != simplex.end(); ++it)
           if (it != not_it)
             facet.push_back(*it);
         Simplex_handle facet_sh = sc.find(facet);