Rewrite for the Weak Witness Complex

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

Former-commit-id: 00d6d906598893ab551157c00020ce33fe19cb5b

2 files changed, 253 insertions, 78 deletions

diff --git a/src/Witness_complex/include/gudhi/Construct_closest_landmark_table.h b/src/Witness_complex/include/gudhi/Construct_closest_landmark_table.h
new file mode 100644
index 00000000..ef711c34
--- /dev/null
+++ b/src/Witness_complex/include/gudhi/Construct_closest_landmark_table.h
@@ -0,0 +1,90 @@
+/*    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):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2015  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/>.
+ */
+
+#ifndef CONSTRUCT_CLOSEST_LANDMARK_TABLE_H_
+#define CONSTRUCT_CLOSEST_LANDMARK_TABLE_H_
+
+#include <boost/range/size.hpp>
+
+#include <queue>  // for priority_queue<>
+#include <utility>  // for pair<>
+#include <iterator>
+#include <vector>
+#include <set>
+
+namespace Gudhi {
+
+namespace witness_complex {
+
+  /**
+   *  \ingroup witness_complex
+   * \brief Construct the closest landmark tables for all witnesses.
+   *  \details Output a table 'knn', each line of which represents a witness and
+   *   consists of landmarks sorted by
+   *   euclidean distance from the corresponding witness.
+   *
+   *  The type WitnessContainer is a random access range and
+   *  the type LandmarkContainer is a range.
+   *  The type KNearestNeighbors can be seen as 
+   *  Witness_range<Closest_landmark_range<Vertex_handle>>, where
+   *  Witness_range and Closest_landmark_range are random access ranges and
+   *  Vertex_handle is the label type of a vertex in a simplicial complex.
+   *  Closest_landmark_range needs to have push_back operation.
+   */
+
+  template <typename WitnessContainer,
+            typename LandmarkContainer,
+            typename KNearestNeighbours>
+  void construct_closest_landmark_table(WitnessContainer const &points,
+                                        LandmarkContainer const &landmarks,
+                                        KNearestNeighbours &knn) {
+    int nbP = boost::size(points);
+    assert(nbP >= boost::size(landmarks));
+
+    int dim = boost::size(*std::begin(points));
+    typedef std::pair<double, int> dist_i;
+    typedef bool (*comp)(dist_i, dist_i);
+    knn = KNearestNeighbours(nbP);
+    for (int points_i = 0; points_i < nbP; points_i++) {
+      std::priority_queue<dist_i, std::vector<dist_i>, comp> l_heap([](dist_i j1, dist_i j2) {
+          return j1.first > j2.first;
+        });
+      typename LandmarkContainer::const_iterator landmarks_it;
+      int landmarks_i = 0;
+      for (landmarks_it = landmarks.begin(), landmarks_i = 0; landmarks_it != landmarks.end();
+           ++landmarks_it, landmarks_i++) {
+        dist_i dist = std::make_pair(euclidean_distance(points[points_i], *landmarks_it), landmarks_i);
+        l_heap.push(dist);
+      }
+      for (int i = 0; i < dim + 1; i++) {
+        dist_i dist = l_heap.top();
+        knn[points_i].push_back(dist.second);
+        l_heap.pop();
+      }
+    }
+  }
+
+}  // namespace witness_complex
+
+}  // namespace Gudhi
+
+#endif  // CONSTRUCT_CLOSEST_LANDMARK_TABLE_H_
diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h
index 5c6d087f..4c489e7a 100644
--- a/src/Witness_complex/include/gudhi/Witness_complex.h
+++ b/src/Witness_complex/include/gudhi/Witness_complex.h
@@ -31,6 +31,7 @@
 #include <boost/range/size.hpp>
 
 #include <gudhi/distance_functions.h>
+#include <gudhi/Kd_tree_search.h>
 
 #include <algorithm>
 #include <utility>
@@ -42,8 +43,10 @@
 #include <ctime>
 #include <iostream>
 
-namespace Gudhi {
+namespace gss = Gudhi::spatial_searching;
 
+namespace Gudhi {
+  
 namespace witness_complex {
 
 // /*
@@ -55,9 +58,15 @@ namespace witness_complex {
 template< class SimplicialComplex,
           class Kernel_ >
 class Witness_complex {
-  typedef Kernel_        K;
-  typedef K::Point_d     Point_d;
-  
+  typedef Kernel_                                                 K;
+  typedef typename K::Point_d                                     Point_d;
+  typedef typename K::FT                                          FT;
+  typedef std::vector<Point_d>                                    Point_range;
+  typedef gss::Kd_tree_search<Kernel_, Point_range>               Kd_tree;
+  typedef typename Kd_tree::INS_range                             Nearest_landmark_range;
+  typedef typename std::vector<Nearest_landmark_range>            Nearest_landmark_table;
+  typedef typename Nearest_landmark_table::const_iterator         Nearest_landmark_table_iterator;
+  //typedef std::vector<std::pair<unsigned,FT>>                     Nearest_landmarks;
   
  private:
   struct Active_witness {
@@ -86,9 +95,10 @@ class Witness_complex {
 
  private:
   int nbL_;  // Number of landmarks
-  SimplicialComplex& sc_;  // Simplicial complex
-  std::vector<Point_d> witnesses_, landmarks_; 
-  
+  SimplicialComplex&             sc_;  // Simplicial complex
+  Point_range                         witnesses_, landmarks_;
+  Kd_tree                             landmark_tree_;
+  std::vector<Nearest_landmark_range> nearest_landmarks_;
   
  public:
   /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -118,52 +128,53 @@ class Witness_complex {
   Witness_complex(InputIteratorLandmarks landmarks_first,
                   InputIteratorLandmarks landmarks_last,
                   InputIteratorWitnesses witnesses_first,
-                  InputIteratorWitnesses witnesses_last)
-    : witnesses_(witnesses_first, witnesses_last), landmarks_(landmarks_first, landmarks_last)
+                  InputIteratorWitnesses witnesses_last,
+                  SimplicialComplex& sc)
+    : witnesses_(witnesses_first, witnesses_last), landmarks_(landmarks_first, landmarks_last), landmark_tree_(landmarks_), sc_(sc)
   {    
+    for (Point_d w : witnesses_)
+      nearest_landmarks_.push_back(landmark_tree_.query_incremental_nearest_neighbors(w));
   }
 
+  Point_d get_point( std::size_t vertex ) const
+  {
+    return landmarks_[vertex];
+  }
   
-  
-  template< typename KNearestNeighbors >
-  Witness_complex(KNearestNeighbors const & knn,
-                  int nbL,
-                  int dim,
-                  SimplicialComplex & sc) : nbL_(nbL), sc_(sc) {
-    // Construction of the active witness list
-    int nbW = boost::size(knn);
+  template < typename SimplicialComplexForWitness >
+  bool create_complex(SimplicialComplexForWitness& complex,
+                      FT  max_alpha_square) const
+  {
+    unsigned nbW = witnesses_.size(),
+             nbL = landmarks_.size();
     typeVectorVertex vv;
-    int counter = 0;
-    /* The list of still useful witnesses
-     * it will diminuish in the course of iterations
-     */
-    ActiveWitnessList active_w;  // = new ActiveWitnessList();
-    for (Vertex_handle i = 0; i != nbL_; ++i) {
+    ActiveWitnessList active_w;// = new ActiveWitnessList();
+    for (int 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++;
+      //counter++;
       vv = {i};
-      sc_.insert_simplex(vv);
-      // TODO(SK) Error if not inserted : normally no need here though
+      complex.insert_simplex(vv, Filtration_value(0.0));
+      /* TODO Error if not inserted : normally no need here though*/
     }
-    int k = 1; /* current dimension in iterative construction */
-    for (int i = 0; i != nbW; ++i)
+    unsigned k = 1; /* current dimension in iterative construction */
+    for (int i=0; i != nbW; ++i)
       active_w.push_back(i);
-    while (!active_w.empty() && k < dim) {
-      typename ActiveWitnessList::iterator it = active_w.begin();
-      while (it != active_w.end()) {
-        typeVectorVertex simplex_vector;
-        /* THE INSERTION: Checking if all the subfaces are in the simplex tree*/
-        bool ok = all_faces_in(knn, *it, k);
-        if (ok) {
-          for (int i = 0; i != k + 1; ++i)
-            simplex_vector.push_back(knn[*it][i]);
-          sc_.insert_simplex(simplex_vector);
-          // TODO(SK) Error if not inserted : normally no need here though
-          ++it;
-        } else {
-          active_w.erase(it++);  // First increase the iterator and then erase the previous element
-        }
+    while (!active_w.empty() && k < nbL ) {
+      typename ActiveWitnessList::iterator aw_it = active_w.begin();
+      while (aw_it != active_w.end()) {
+        std::vector<int> simplex;
+        bool ok = add_all_faces_of_dimension(k,
+                                             max_alpha_square,
+                                             std::numeric_limits<double>::infinity(),
+                                             nearest_landmarks_[*aw_it].begin(),
+                                             simplex,
+                                             complex,
+                                             nearest_landmarks_[*aw_it].end());
+        if (!ok)
+          active_w.erase(aw_it++); //First increase the iterator and then erase the previous element
+        else
+          aw_it++;
       }
       k++;
     }
@@ -172,40 +183,114 @@ class Witness_complex {
   //@}
 
  private:
-  /* \brief Check if the facets of the k-dimensional simplex witnessed 
-   *  by witness witness_id are already in the complex.
-   *  inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id
+
+  /* \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
+   * The output value indicates if the witness rests active or not
    */
-  template <typename KNearestNeighbors>
-  bool all_faces_in(KNearestNeighbors const &knn, int witness_id, int k) {
-    std::vector< Vertex_handle > facet;
-    // CHECK ALL THE FACETS
-    for (int i = 0; i != k + 1; ++i) {
-      facet = {};
-      for (int j = 0; j != k + 1; ++j) {
-        if (j != i) {
-          facet.push_back(knn[witness_id][j]);
+  template < typename SimplicialComplexForWitness >
+  bool add_all_faces_of_dimension(int dim,
+                                  double alpha2,
+                                  double norelax_dist2,
+                                  Nearest_landmark_table_iterator curr_l,
+                                  std::vector<int>& simplex,
+                                  SimplicialComplexForWitness& sc,
+                                  Nearest_landmark_table_iterator end)
+  {
+    if (curr_l == end)
+      return false;
+    bool will_be_active = false;
+    Nearest_landmark_table_iterator l_it = curr_l;
+    if (dim > 0)
+      for (; l_it->second - alpha2 <= norelax_dist2 && l_it != end; ++l_it) {
+        simplex.push_back(l_it->first);
+        if (sc.find(simplex) != sc.null_simplex())
+          will_be_active = add_all_faces_of_dimension(dim-1,
+                                                      alpha2,
+                                                      norelax_dist2,
+                                                      l_it+1,
+                                                      simplex,
+                                                      sc,
+                                                      end) || will_be_active;
+        simplex.pop_back();
+        // If norelax_dist is infinity, change to first omitted distance
+        if (l_it->second <= norelax_dist2)
+          norelax_dist2 = l_it->second;
+        will_be_active = add_all_faces_of_dimension(dim-1,
+                                                    alpha2,
+                                                    norelax_dist2,
+                                                    l_it+1,
+                                                    simplex,
+                                                    sc,
+                                                    end) || will_be_active;
+      } 
+    else if (dim == 0)
+      for (; l_it->second - alpha2 <= norelax_dist2 && l_it != end; ++l_it) {
+        simplex.push_back(l_it->second);
+        double filtration_value = 0;
+        // if norelax_dist is infinite, relaxation is 0.
+        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);
         }
-      }  // endfor
-      if (sc_.find(facet) == sc_.null_simplex())
-        return false;
-    }  // endfor
-    return true;
+        simplex.pop_back();
+        // If norelax_dist is infinity, change to first omitted distance
+        if (l_it->second < norelax_dist2)
+          norelax_dist2 = l_it->second;
+      } 
+    return will_be_active;
   }
-
-  template <typename T>
-  static void print_vector(const std::vector<T>& v) {
-    std::cout << "[";
-    if (!v.empty()) {
-      std::cout << *(v.begin());
-      for (auto it = v.begin() + 1; it != v.end(); ++it) {
-        std::cout << ",";
-        std::cout << *it;
+  
+  /** \brief Check if the facets of the k-dimensional simplex witnessed 
+   *  by witness witness_id are already in the complex.
+   *  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,
+                    double* filtration_value,
+                    SimplicialComplexForWitness& sc)
+  {
+    std::vector< int > facet;
+    for (std::vector<int>::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)
+          if (it != not_it)
+            facet.push_back(*it);
+        Simplex_handle facet_sh = sc.find(facet);
+        if (facet_sh == sc.null_simplex())
+          return false;
+        else if (sc.filtration(facet_sh) > *filtration_value)
+          *filtration_value = sc.filtration(facet_sh);
       }
-    }
-    std::cout << "]";
+    return true;
   }
 
+  // bool is_face(Simplex_handle face, Simplex_handle coface)
+  // {
+  //   // vertex range is sorted in decreasing order
+  //   auto fvr = sc.simplex_vertex_range(face);
+  //   auto cfvr = sc.simplex_vertex_range(coface);
+  //   auto fv_it = fvr.begin();
+  //   auto cfv_it = cfvr.begin();
+  //   while (fv_it != fvr.end() && cfv_it != cfvr.end()) {
+  //     if (*fv_it < *cfv_it)
+  //       ++cfv_it;
+  //     else if (*fv_it == *cfv_it) {
+  //       ++cfv_it;
+  //       ++fv_it;
+  //     }
+  //     else
+  //       return false;
+      
+  //   }
+  //   return (fv_it == fvr.end());
+  // }
+  
+  
  public:
   // /*
   //  *  \brief Verification if every simplex in the complex is witnessed by witnesses in knn.
@@ -269,13 +354,13 @@ class Witness_complex {
    *
    *  Landmarks are supposed to be in [0,nbL_-1]
    */
-  template <class KNearestNeighbors, class SimplicialComplexForWitness>
-  void witness_complex(KNearestNeighbors const & knn,
-                       int nbL,
-                       int dim,
-                       SimplicialComplexForWitness & sc) {
-    Witness_complex<SimplicialComplexForWitness>(knn, nbL, dim, sc);
-  }
+  // template <class KNearestNeighbors, class SimplicialComplexForWitness>
+  // void witness_complex(KNearestNeighbors const & knn,
+  //                      int nbL,
+  //                      int dim,
+  //                      SimplicialComplexForWitness & sc) {
+  //   Witness_complex<SimplicialComplexForWitness>(knn, nbL, dim, sc);
+  // }
 
 }  // namespace witness_complex