1 files changed, 60 insertions, 23 deletions

diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
index 0d6299d2..2bf62f9b 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
@@ -13,6 +13,8 @@
 
 #include <gudhi/Bitmap_cubical_complex/counter.h>
 
+#include <boost/config.hpp>
+
 #include <iostream>
 #include <vector>
 #include <string>
@@ -41,7 +43,7 @@ namespace cubical_complex {
  * Each cell is represented by a single
  * bit (in case of black and white bitmaps, or by a single element of a type T
  * (here T is a filtration type of a bitmap, typically a double).
- * All the informations needed for homology and
+ * All the information needed for homology and
  * persistent homology computations (like dimension of a cell, boundary and
  * coboundary elements of a cell, are then obtained from the
  * position of the element in C.
@@ -110,6 +112,16 @@ class Bitmap_cubical_complex_base {
   virtual inline std::vector<std::size_t> get_coboundary_of_a_cell(std::size_t cell) const;
 
   /**
+   * This function finds a top-dimensional cell that is incident to the input cell and has
+   * the same filtration value. In case several cells are suitable, an arbitrary one is
+   * returned. Note that the input parameter can be a cell of any dimension (vertex, edge, etc).
+   * On the other hand, the output is always indicating the position of
+   * a top-dimensional cube in the data structure.
+   * \pre The filtration values are assigned as per `impose_lower_star_filtration()`.
+   **/
+  inline size_t get_top_dimensional_coface_of_a_cell(size_t splx);
+
+  /**
   * This procedure compute incidence numbers between cubes. For a cube \f$A\f$ of
   * dimension n and a cube \f$B \subset A\f$ of dimension n-1, an incidence
   * between \f$A\f$ and \f$B\f$ is the integer with which \f$B\f$ appears in the boundary of \f$A\f$.
@@ -142,7 +154,7 @@ class Bitmap_cubical_complex_base {
       }
       if (coface_counter[i] != face_counter[i]) {
         if (number_of_position_in_which_counters_do_not_agree != -1) {
-          std::cout << "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.\n";
+          std::cerr << "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.\n";
           throw std::logic_error(
               "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.");
         }
@@ -197,7 +209,7 @@ class Bitmap_cubical_complex_base {
   /**
    * Returns number of all cubes in the data structure.
    **/
-  inline unsigned size() const { return this->data.size(); }
+  inline std::size_t size() const { return this->data.size(); }
 
   /**
    * Writing to stream operator. By using it we get the values T of cells in order in which they are stored in the
@@ -239,8 +251,14 @@ class Bitmap_cubical_complex_base {
    * @brief Iterator through all cells in the complex (in order they appear in the structure -- i.e.
    * in lexicographical order).
    **/
-  class All_cells_iterator : std::iterator<std::input_iterator_tag, T> {
+  class All_cells_iterator {
    public:
+    typedef std::input_iterator_tag iterator_category;
+    typedef std::size_t value_type;
+    typedef std::ptrdiff_t difference_type;
+    typedef value_type* pointer;
+    typedef value_type reference;
+
     All_cells_iterator() { this->counter = 0; }
 
     All_cells_iterator operator++() {
@@ -343,8 +361,14 @@ class Bitmap_cubical_complex_base {
    * @brief Iterator through top dimensional cells of the complex. The cells appear in order they are stored
    * in the structure (i.e. in lexicographical order)
    **/
-  class Top_dimensional_cells_iterator : std::iterator<std::input_iterator_tag, T> {
+  class Top_dimensional_cells_iterator {
    public:
+    typedef std::input_iterator_tag iterator_category;
+    typedef std::size_t value_type;
+    typedef std::ptrdiff_t difference_type;
+    typedef value_type* pointer;
+    typedef value_type reference;
+
     Top_dimensional_cells_iterator(Bitmap_cubical_complex_base& b) : b(b) {
       this->counter = std::vector<std::size_t>(b.dimension());
       // std::fill( this->counter.begin() , this->counter.end() , 0 );
@@ -408,7 +432,7 @@ class Bitmap_cubical_complex_base {
 
     void print_counter() const {
       for (std::size_t i = 0; i != this->counter.size(); ++i) {
-        std::cout << this->counter[i] << " ";
+        std::clog << this->counter[i] << " ";
       }
     }
     friend class Bitmap_cubical_complex_base;
@@ -521,11 +545,11 @@ void Bitmap_cubical_complex_base<T>::put_data_to_bins(std::size_t number_of_bins
   // now put the data into the appropriate bins:
   for (std::size_t i = 0; i != this->data.size(); ++i) {
     if (dbg) {
-      std::cerr << "Before binning : " << this->data[i] << std::endl;
+      std::clog << "Before binning : " << this->data[i] << std::endl;
     }
     this->data[i] = min_max.first + dx * (this->data[i] - min_max.first) / number_of_bins;
     if (dbg) {
-      std::cerr << "After binning : " << this->data[i] << std::endl;
+      std::clog << "After binning : " << this->data[i] << std::endl;
     }
   }
 }
@@ -539,11 +563,11 @@ void Bitmap_cubical_complex_base<T>::put_data_to_bins(T diameter_of_bin) {
   // now put the data into the appropriate bins:
   for (std::size_t i = 0; i != this->data.size(); ++i) {
     if (dbg) {
-      std::cerr << "Before binning : " << this->data[i] << std::endl;
+      std::clog << "Before binning : " << this->data[i] << std::endl;
     }
     this->data[i] = min_max.first + diameter_of_bin * (this->data[i] - min_max.first) / number_of_bins;
     if (dbg) {
-      std::cerr << "After binning : " << this->data[i] << std::endl;
+      std::clog << "After binning : " << this->data[i] << std::endl;
     }
   }
 }
@@ -603,6 +627,19 @@ void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells
 }
 
 template <typename T>
+size_t Bitmap_cubical_complex_base<T>::get_top_dimensional_coface_of_a_cell(size_t splx) {
+  if (this->get_dimension_of_a_cell(splx) == this->dimension()){return splx;}
+  else{
+    for (auto v : this->get_coboundary_of_a_cell(splx)){
+      if(this->get_cell_data(v) == this->get_cell_data(splx)){
+        return this->get_top_dimensional_coface_of_a_cell(v);
+      }
+    }
+  }
+  BOOST_UNREACHABLE_RETURN(-2);
+}
+
+template <typename T>
 Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes_in_following_directions,
                                                             const std::vector<T>& top_dimensional_cells) {
   this->setup_bitmap_based_on_top_dimensional_cells_list(sizes_in_following_directions, top_dimensional_cells);
@@ -617,7 +654,7 @@ void Bitmap_cubical_complex_base<T>::read_perseus_style_file(const char* perseus
   inFiltration >> dimensionOfData;
 
   if (dbg) {
-    std::cerr << "dimensionOfData : " << dimensionOfData << std::endl;
+    std::clog << "dimensionOfData : " << dimensionOfData << std::endl;
   }
 
   std::vector<unsigned> sizes;
@@ -630,7 +667,7 @@ void Bitmap_cubical_complex_base<T>::read_perseus_style_file(const char* perseus
     sizes.push_back(size_in_this_dimension);
     dimensions *= size_in_this_dimension;
     if (dbg) {
-      std::cerr << "size_in_this_dimension : " << size_in_this_dimension << std::endl;
+      std::clog << "size_in_this_dimension : " << size_in_this_dimension << std::endl;
     }
   }
   this->set_up_containers(sizes);
@@ -651,7 +688,7 @@ void Bitmap_cubical_complex_base<T>::read_perseus_style_file(const char* perseus
       }
 
       if (dbg) {
-        std::cerr << "Cell of an index : " << it.compute_index_in_bitmap()
+        std::clog << "Cell of an index : " << it.compute_index_in_bitmap()
                   << " and dimension: " << this->get_dimension_of_a_cell(it.compute_index_in_bitmap())
                   << " get the value : " << filtrationLevel << std::endl;
       }
@@ -754,20 +791,20 @@ std::vector<std::size_t> Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cel
 template <typename T>
 unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell(std::size_t cell) const {
   bool dbg = false;
-  if (dbg) std::cerr << "\n\n\n Computing position o a cell of an index : " << cell << std::endl;
+  if (dbg) std::clog << "\n\n\n Computing position o a cell of an index : " << cell << std::endl;
   unsigned dimension = 0;
   for (std::size_t i = this->multipliers.size(); i != 0; --i) {
     unsigned position = cell / this->multipliers[i - 1];
 
     if (dbg) {
-      std::cerr << "i-1 :" << i - 1 << std::endl;
-      std::cerr << "cell : " << cell << std::endl;
-      std::cerr << "position : " << position << std::endl;
-      std::cerr << "multipliers[" << i - 1 << "] = " << this->multipliers[i - 1] << std::endl;
+      std::clog << "i-1 :" << i - 1 << std::endl;
+      std::clog << "cell : " << cell << std::endl;
+      std::clog << "position : " << position << std::endl;
+      std::clog << "multipliers[" << i - 1 << "] = " << this->multipliers[i - 1] << std::endl;
     }
 
     if (position % 2 == 1) {
-      if (dbg) std::cerr << "Nonzero length in this direction \n";
+      if (dbg) std::clog << "Nonzero length in this direction \n";
       dimension++;
     }
     cell = cell % this->multipliers[i - 1];
@@ -803,9 +840,9 @@ void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration() {
 
   while (indices_to_consider.size()) {
     if (dbg) {
-      std::cerr << "indices_to_consider in this iteration \n";
+      std::clog << "indices_to_consider in this iteration \n";
       for (std::size_t i = 0; i != indices_to_consider.size(); ++i) {
-        std::cout << indices_to_consider[i] << "  ";
+        std::clog << indices_to_consider[i] << "  ";
       }
     }
     std::vector<std::size_t> new_indices_to_consider;
@@ -813,14 +850,14 @@ void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration() {
       std::vector<std::size_t> bd = this->get_boundary_of_a_cell(indices_to_consider[i]);
       for (std::size_t boundaryIt = 0; boundaryIt != bd.size(); ++boundaryIt) {
         if (dbg) {
-          std::cerr << "filtration of a cell : " << bd[boundaryIt] << " is : " << this->data[bd[boundaryIt]]
+          std::clog << "filtration of a cell : " << bd[boundaryIt] << " is : " << this->data[bd[boundaryIt]]
                     << " while of a cell: " << indices_to_consider[i] << " is: " << this->data[indices_to_consider[i]]
                     << std::endl;
         }
         if (this->data[bd[boundaryIt]] > this->data[indices_to_consider[i]]) {
           this->data[bd[boundaryIt]] = this->data[indices_to_consider[i]];
           if (dbg) {
-            std::cerr << "Setting the value of a cell : " << bd[boundaryIt]
+            std::clog << "Setting the value of a cell : " << bd[boundaryIt]
                       << " to : " << this->data[indices_to_consider[i]] << std::endl;
           }
         }