Fixed a major bug in the Active_witness_iterator

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

Former-commit-id: d4f474b44835b11e65a152f9d3fe87cc32b527f5

5 files changed, 116 insertions, 39 deletions

diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp
index 310fb1be..fe7f628f 100644
--- a/src/Witness_complex/example/witness_complex_sphere.cpp
+++ b/src/Witness_complex/example/witness_complex_sphere.cpp
@@ -85,7 +85,7 @@ int main(int argc, char * const argv[]) {
                                     landmarks.end(),
                                     point_vector.begin(),
                                     point_vector.end());
-    witness_complex.create_complex(simplex_tree, 0);
+    witness_complex.create_complex(simplex_tree, 0.1);
     end = clock();
     double time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
     std::cout << "Witness complex for " << number_of_landmarks << " landmarks took "
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 7b169784..9ae41a69 100644
--- a/src/Witness_complex/include/gudhi/Active_witness/Active_witness.h
+++ b/src/Witness_complex/include/gudhi/Active_witness/Active_witness.h
@@ -34,23 +34,26 @@ namespace witness_complex {
 template< typename Id_distance_pair,
           typename INS_range >
 class Active_witness {
-public:
+public:  
   typedef Active_witness<Id_distance_pair, INS_range> ActiveWitness;
   typedef typename INS_range::iterator INS_iterator;
   typedef Active_witness_iterator< ActiveWitness, Id_distance_pair, INS_iterator > iterator;
   typedef typename std::list<Id_distance_pair> Table;
-  
+
+  Table end_element_table_ = {Id_distance_pair(-1,0)};
+  typename Table::iterator end_pointer = end_element_table_.begin();
+
   Table nearest_landmark_table_;
   INS_range    search_range_;
   INS_iterator iterator_last_;
   INS_iterator iterator_end_;
 
   Active_witness(INS_range search_range)
-    : search_range_(search_range), iterator_end_(search_range.end()), iterator_last_(search_range.begin())
+    : search_range_(search_range), iterator_last_(search_range.begin()), iterator_end_(search_range.end())
   {
     nearest_landmark_table_.push_back(*iterator_last_);
   }
-
+ 
   iterator begin()
   {
     return iterator(this, nearest_landmark_table_.begin());
@@ -60,9 +63,6 @@ public:
   {
     return iterator(this);
   }
-
-  Table end_element_table_ = {Id_distance_pair(-1,0)};
-  typename Table::iterator end_pointer = end_element_table_.begin();
 };
 
 }
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 658405f6..5d4f3d75 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
@@ -64,7 +64,7 @@ public:
   
 private :
 
-  Id_distance_pair dereference() const
+  Id_distance_pair& dereference() const
   {
     return *lh_;
   }
@@ -76,9 +76,17 @@ private :
   
   void increment()
   {
+    // if neighbor search is at its end, check if lh_++ is end
+    if (aw_->iterator_last_ == aw_->iterator_end_) {
+      if (lh_++ == aw_->nearest_landmark_table_.end()) {
+        lh_ = aw_->end_pointer;
+        return;
+      }
+      return;
+    }
     // if the id of the current landmark is the same as the last one
     if (lh_->first == aw_->iterator_last_->first) {
-      // if the next iterator is end, lh_it = nullptr
+      // if the next iterator is end, lh_it = end pointer
       if (++(aw_->iterator_last_) == aw_->iterator_end_) {
         lh_ = aw_->end_pointer;
         return;
diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h
index 158f92ec..e939de34 100644
--- a/src/Witness_complex/include/gudhi/Witness_complex.h
+++ b/src/Witness_complex/include/gudhi/Witness_complex.h
@@ -84,9 +84,9 @@ private:
   typedef FT Filtration_value;
 
   
-  typedef std::size_t Witness_id;
-  typedef std::size_t Landmark_id;
-  typedef std::pair<std::size_t, FT> Id_distance_pair;
+  typedef std::ptrdiff_t Witness_id;
+  typedef std::ptrdiff_t Landmark_id;
+  typedef std::pair<Landmark_id, FT> Id_distance_pair;
   typedef Active_witness<Id_distance_pair, Nearest_landmark_range> ActiveWitness;
   typedef std::list< ActiveWitness > ActiveWitnessList;
   typedef std::vector< Landmark_id > typeVectorVertex;
@@ -155,7 +155,7 @@ private:
     }
     typeVectorVertex vv;
     ActiveWitnessList active_witnesses;// = new ActiveWitnessList();
-    for (auto i = 0; i != nbL; ++i) {
+    for (unsigned i = 0; i != nbL; ++i) {
       // initial fill of 0-dimensional simplices
       // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore
       //counter++;
@@ -166,10 +166,12 @@ private:
     unsigned 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 ) {
       typename ActiveWitnessList::iterator aw_it = active_witnesses.begin();
       while (aw_it != active_witnesses.end()) {
         std::vector<int> simplex;
+        //simplex.reserve(k+1);
         bool ok = add_all_faces_of_dimension(k,
                                              max_alpha_square,
                                              std::numeric_limits<double>::infinity(),
@@ -178,10 +180,13 @@ private:
                                              complex,
                                              aw_it->end());
         if (!ok)
-          active_witnesses.erase(aw_it++); //First increase the iterator and then erase the previous element
+          //{aw_it++;}
+        active_witnesses.erase(aw_it++); //First increase the iterator and then erase the previous element
         else
           aw_it++;
-      }
+      } 
+      std::cout << "Active witnesses after dim=" << k << " is finished: " << active_witnesses.size() << "\n";
+      std::cout << complex << "\n";
       k++;
     }
     return true;
@@ -190,7 +195,6 @@ private:
   //@}
 
  private:
-
   /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness
    * Iterator is needed to know until how far we can take landmarks to form simplexes
    * simplex is the prefix of the simplexes to insert
@@ -222,12 +226,13 @@ private:
                                                       sc,
                                                       end) || will_be_active;
         }
+        assert(!simplex.empty());
         simplex.pop_back();
         // If norelax_dist is infinity, change to first omitted distance
         if (l_it->second <= norelax_dist2)
           norelax_dist2 = l_it->second;
         typename ActiveWitness::iterator next_it = l_it;
-        will_be_active = add_all_faces_of_dimension(dim-1,
+        will_be_active = add_all_faces_of_dimension(dim,
                                                     alpha2,
                                                     norelax_dist2,
                                                     ++next_it,
@@ -240,12 +245,16 @@ private:
         simplex.push_back(l_it->first);
         double filtration_value = 0;
         // if norelax_dist is infinite, relaxation is 0.
+        //std::cout << "landmark_id=" << l_it->first << " distance=" << l_it->second << "\n";
+        // std::size_t landmark_id = l_it->first;
+        // double distance = l_it->second;
         if (l_it->second > norelax_dist2) 
           filtration_value = l_it->second - norelax_dist2; 
         if (all_faces_in(simplex, &filtration_value, sc)) {
           will_be_active = true;
           sc.insert_simplex(simplex, filtration_value);
         }
+        assert(!simplex.empty());
         simplex.pop_back();
         // If norelax_dist is infinity, change to first omitted distance
         if (l_it->second < norelax_dist2)
@@ -301,9 +310,23 @@ private:
   //   }
   //   return (fv_it == fvr.end());
   // }
-  
+
   
  public:
+  template < typename SimplicialComplexForWitness >
+  void print_complex(SimplicialComplexForWitness& complex)
+  {
+    std::cout << complex << "\n";
+  }
+
+  template < typename Container >
+  void print_container(Container& container)
+  {
+    for (auto l: container)
+      std::cout << l << ", ";
+    std::cout << "\n";
+  }
+
   // /*
   //  *  \brief Verification if every simplex in the complex is witnessed by witnesses in knn.
   //  *  \param print_output =true will print the witnesses for each simplex
diff --git a/src/Witness_complex/test/simple_witness_complex.cpp b/src/Witness_complex/test/simple_witness_complex.cpp
index 03df78ee..30b168c1 100644
--- a/src/Witness_complex/test/simple_witness_complex.cpp
+++ b/src/Witness_complex/test/simple_witness_complex.cpp
@@ -25,6 +25,8 @@
 #include <boost/test/unit_test.hpp>
 #include <boost/mpl/list.hpp>
 
+#include <CGAL/Epick_d.h>
+
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Witness_complex.h>
 
@@ -34,26 +36,70 @@
 
 typedef Gudhi::Simplex_tree<> Simplex_tree;
 typedef std::vector< Vertex_handle > typeVectorVertex;
-typedef Gudhi::witness_complex::Witness_complex<Simplex_tree> WitnessComplex;
+typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> Kernel;
+typedef typename Kernel::FT FT;
+typedef typename Kernel::Point_d Point_d;
+typedef Gudhi::witness_complex::Witness_complex<Kernel> WitnessComplex;
+
+/* All landmarks and witnesses are taken on the grid in the following manner.
+   LWLWL  2W4W7
+   WW.WW  WW.WW
+   L...L  1...6
+   WW.WW  WW.WW
+   LWLWL  0W3W5
+
+   Witness complex consists of 8 vertices, 12 edges and 4 triangles
+ */
 
 BOOST_AUTO_TEST_CASE(simple_witness_complex) {
-  Simplex_tree complex;
-  std::vector< typeVectorVertex > knn;
-
-  knn.push_back({1, 0, 5, 2, 6, 3, 4});
-  knn.push_back({2, 6, 4, 5, 0, 1, 3});
-  knn.push_back({3, 4, 2, 1, 5, 6, 0});
-  knn.push_back({4, 2, 1, 3, 5, 6, 0});
-  knn.push_back({5, 1, 6, 0, 2, 3, 4});
-  knn.push_back({6, 0, 5, 2, 1, 3, 4});
-  knn.push_back({0, 5, 6, 1, 2, 3, 4});
-  knn.push_back({2, 6, 4, 5, 3, 1, 0});
-  knn.push_back({1, 2, 5, 4, 3, 6, 0});
-  knn.push_back({3, 4, 0, 6, 5, 1, 2});
-  knn.push_back({5, 0, 1, 3, 6, 2, 4});
-  knn.push_back({5, 6, 1, 0, 2, 3, 4});
-  knn.push_back({1, 6, 0, 5, 2, 3, 4});
-  WitnessComplex witnessComplex(knn, 7, 7, complex);
-
-  BOOST_CHECK(witnessComplex.is_witness_complex(knn, false));
+  Simplex_tree complex, relaxed_complex;
+
+  std::vector<Point_d> witnesses, landmarks;
+
+  landmarks.push_back(Point_d(std::vector<FT>{-2,-2}));
+  landmarks.push_back(Point_d(std::vector<FT>{-2, 0}));
+  landmarks.push_back(Point_d(std::vector<FT>{-2, 2}));
+  landmarks.push_back(Point_d(std::vector<FT>{ 0,-2}));
+  landmarks.push_back(Point_d(std::vector<FT>{ 0, 2}));
+  landmarks.push_back(Point_d(std::vector<FT>{ 2,-2}));
+  landmarks.push_back(Point_d(std::vector<FT>{ 2, 0}));
+  landmarks.push_back(Point_d(std::vector<FT>{ 2, 2}));
+  witnesses.push_back(Point_d(std::vector<FT>{-2,-1}));
+  witnesses.push_back(Point_d(std::vector<FT>{-2, 1}));
+  witnesses.push_back(Point_d(std::vector<FT>{-1,-2}));
+  witnesses.push_back(Point_d(std::vector<FT>{-1,-1}));
+  witnesses.push_back(Point_d(std::vector<FT>{-1, 1}));
+  witnesses.push_back(Point_d(std::vector<FT>{-1, 2}));
+  witnesses.push_back(Point_d(std::vector<FT>{ 1,-2}));
+  witnesses.push_back(Point_d(std::vector<FT>{ 1,-1}));
+  witnesses.push_back(Point_d(std::vector<FT>{ 1, 1}));
+  witnesses.push_back(Point_d(std::vector<FT>{ 1, 2}));
+  witnesses.push_back(Point_d(std::vector<FT>{ 2,-1}));
+  witnesses.push_back(Point_d(std::vector<FT>{ 2, 1}));
+
+  // landmarks.push_back(Point_d(std::vector<FT>{1,0}));
+  // landmarks.push_back(Point_d(std::vector<FT>{0,0}));
+  // landmarks.push_back(Point_d(std::vector<FT>{0,1}));
+  // witnesses.push_back(Point_d(std::vector<FT>{1,0}));
+  
+  WitnessComplex witness_complex(landmarks.begin(),
+                                 landmarks.end(),
+                                 witnesses.begin(),
+                                 witnesses.end());
+  // witness_complex.create_complex(complex, 0);
+
+  //std::cout << complex << "\n";
+
+  // BOOST_CHECK(complex.num_simplices() == 24);
+
+  witness_complex.create_complex(relaxed_complex, 8.01);
+  
+  std::cout << "Num_simplices: " << relaxed_complex.num_simplices() << "\n";
+  std::cout << relaxed_complex << "\n";
+
+  //BOOST_CHECK(relaxed_complex.num_simplices() == 24);
+  
+  //BOOST_CHECK(witnessComplex.is_witness_complex(knn, false));
+
+  
 }