Bitmap_cubical_complex_base.h

/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
 *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
 *    Author(s):       Pawel Dlotko
 *
 *    Copyright (C) 2015 Inria
 *
 *    Modification(s):
 *      - YYYY/MM Author: Description of the modification
 */
 
#ifndef BITMAP_CUBICAL_COMPLEX_BASE_H_
#define BITMAP_CUBICAL_COMPLEX_BASE_H_
 
#include <gudhi/Bitmap_cubical_complex/counter.h>
 
#include <boost/config.hpp>
 
#include <iostream>
#include <vector>
#include <string>
#include <fstream>
#include <algorithm>
#include <iterator>
#include <limits>
#include <utility>
#include <stdexcept>
#include <cstddef>
 
namespace Gudhi {
 
namespace cubical_complex {
 
template <typename T>
class Bitmap_cubical_complex_base {
 public:
  typedef T filtration_type;
 
  Bitmap_cubical_complex_base() : total_number_of_cells(0) {}
  Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes);
  Bitmap_cubical_complex_base(const char* perseus_style_file);
  Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, const std::vector<T>& top_dimensional_cells);
 
  virtual ~Bitmap_cubical_complex_base() {}
 
  virtual inline std::vector<std::size_t> get_boundary_of_a_cell(std::size_t cell) const;
  virtual inline std::vector<std::size_t> get_coboundary_of_a_cell(std::size_t cell) const;
 
  inline size_t get_top_dimensional_coface_of_a_cell(size_t splx);
 
  virtual int compute_incidence_between_cells(std::size_t coface, std::size_t face) const {
    // first get the counters for coface and face:
    std::vector<unsigned> coface_counter = this->compute_counter_for_given_cell(coface);
    std::vector<unsigned> face_counter = this->compute_counter_for_given_cell(face);
 
    // coface_counter and face_counter should agree at all positions except from one:
    int number_of_position_in_which_counters_do_not_agree = -1;
    std::size_t number_of_full_faces_that_comes_before = 0;
    for (std::size_t i = 0; i != coface_counter.size(); ++i) {
      if ((coface_counter[i] % 2 == 1) && (number_of_position_in_which_counters_do_not_agree == -1)) {
        ++number_of_full_faces_that_comes_before;
      }
      if (coface_counter[i] != face_counter[i]) {
        if (number_of_position_in_which_counters_do_not_agree != -1) {
          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.");
        }
        number_of_position_in_which_counters_do_not_agree = i;
      }
    }
 
    int incidence = 1;
    if (number_of_full_faces_that_comes_before % 2) incidence = -1;
    // if the face cell is on the right from coface cell:
    if (coface_counter[number_of_position_in_which_counters_do_not_agree] + 1 ==
        face_counter[number_of_position_in_which_counters_do_not_agree]) {
      incidence *= -1;
    }
 
    return incidence;
  }
 
  inline unsigned get_dimension_of_a_cell(std::size_t cell) const;
 
  inline T& get_cell_data(std::size_t cell);
 
  void impose_lower_star_filtration();  // assume that top dimensional cells are already set.
 
  inline unsigned dimension() const { return sizes.size(); }
 
  inline std::size_t size() const { return this->data.size(); }
 
  template <typename K>
  friend std::ostream& operator<<(std::ostream& os, const Bitmap_cubical_complex_base<K>& b);
 
  void put_data_to_bins(std::size_t number_of_bins);
 
  void put_data_to_bins(T diameter_of_bin);
 
  std::pair<T, T> min_max_filtration();
 
  // ITERATORS
 
  class All_cells_iterator : std::iterator<std::input_iterator_tag, T> {
   public:
    All_cells_iterator() { this->counter = 0; }
 
    All_cells_iterator operator++() {
      // first find first element of the counter that can be increased:
      ++this->counter;
      return *this;
    }
 
    All_cells_iterator operator++(int) {
      All_cells_iterator result = *this;
      ++(*this);
      return result;
    }
 
    All_cells_iterator& operator=(const All_cells_iterator& rhs) {
      this->counter = rhs.counter;
      return *this;
    }
 
    bool operator==(const All_cells_iterator& rhs) const {
      if (this->counter != rhs.counter) return false;
      return true;
    }
 
    bool operator!=(const All_cells_iterator& rhs) const { return !(*this == rhs); }
 
    /*
     * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as
     * an argument of the following functions:
     * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell
     * boundary and coboundary and dimension
     * and in function get_cell_data to get a filtration of a cell.
     */
    std::size_t operator*() { return this->counter; }
    friend class Bitmap_cubical_complex_base;
 
   protected:
    std::size_t counter;
  };
 
  All_cells_iterator all_cells_iterator_begin() {
    All_cells_iterator a;
    return a;
  }
 
  All_cells_iterator all_cells_iterator_end() {
    All_cells_iterator a;
    a.counter = this->data.size();
    return a;
  }
 
  class All_cells_range {
   public:
    All_cells_range(Bitmap_cubical_complex_base* b) : b(b) {}
 
    All_cells_iterator begin() { return b->all_cells_iterator_begin(); }
 
    All_cells_iterator end() { return b->all_cells_iterator_end(); }
 
   private:
    Bitmap_cubical_complex_base<T>* b;
  };
 
  All_cells_range all_cells_range() { return All_cells_range(this); }
 
  typedef typename std::vector<std::size_t>::const_iterator Boundary_iterator;
  typedef typename std::vector<std::size_t> Boundary_range;
 
  Boundary_range boundary_range(std::size_t sh) { return this->get_boundary_of_a_cell(sh); }
 
  typedef typename std::vector<std::size_t>::const_iterator Coboundary_iterator;
  typedef typename std::vector<std::size_t> Coboundary_range;
 
  Coboundary_range coboundary_range(std::size_t sh) { return this->get_coboundary_of_a_cell(sh); }
 
  class Top_dimensional_cells_iterator : std::iterator<std::input_iterator_tag, T> {
   public:
    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 );
    }
 
    Top_dimensional_cells_iterator operator++() {
      // first find first element of the counter that can be increased:
      std::size_t dim = 0;
      while ((dim != this->b.dimension()) && (this->counter[dim] == this->b.sizes[dim] - 1)) ++dim;
 
      if (dim != this->b.dimension()) {
        ++this->counter[dim];
        for (std::size_t i = 0; i != dim; ++i) {
          this->counter[i] = 0;
        }
      } else {
        ++this->counter[0];
      }
      return *this;
    }
 
    Top_dimensional_cells_iterator operator++(int) {
      Top_dimensional_cells_iterator result = *this;
      ++(*this);
      return result;
    }
 
    Top_dimensional_cells_iterator& operator=(const Top_dimensional_cells_iterator& rhs) {
      this->counter = rhs.counter;
      this->b = rhs.b;
      return *this;
    }
 
    bool operator==(const Top_dimensional_cells_iterator& rhs) const {
      if (&this->b != &rhs.b) return false;
      if (this->counter.size() != rhs.counter.size()) return false;
      for (std::size_t i = 0; i != this->counter.size(); ++i) {
        if (this->counter[i] != rhs.counter[i]) return false;
      }
      return true;
    }
 
    bool operator!=(const Top_dimensional_cells_iterator& rhs) const { return !(*this == rhs); }
 
    /*
     * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as
     * an argument of the following functions:
     * get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell to get information about the cell
     * boundary and coboundary and dimension
     * and in function get_cell_data to get a filtration of a cell.
     */
    std::size_t operator*() { return this->compute_index_in_bitmap(); }
 
    std::size_t compute_index_in_bitmap() const {
      std::size_t index = 0;
      for (std::size_t i = 0; i != this->counter.size(); ++i) {
        index += (2 * this->counter[i] + 1) * this->b.multipliers[i];
      }
      return index;
    }
 
    void print_counter() const {
      for (std::size_t i = 0; i != this->counter.size(); ++i) {
        std::clog << this->counter[i] << " ";
      }
    }
    friend class Bitmap_cubical_complex_base;
 
   protected:
    std::vector<std::size_t> counter;
    Bitmap_cubical_complex_base& b;
  };
 
  Top_dimensional_cells_iterator top_dimensional_cells_iterator_begin() {
    Top_dimensional_cells_iterator a(*this);
    return a;
  }
 
  Top_dimensional_cells_iterator top_dimensional_cells_iterator_end() {
    Top_dimensional_cells_iterator a(*this);
    for (std::size_t i = 0; i != this->dimension(); ++i) {
      a.counter[i] = this->sizes[i] - 1;
    }
    a.counter[0]++;
    return a;
  }
 
  class Top_dimensional_cells_range {
   public:
    Top_dimensional_cells_range(Bitmap_cubical_complex_base* b) : b(b) {}
 
    Top_dimensional_cells_iterator begin() { return b->top_dimensional_cells_iterator_begin(); }
 
    Top_dimensional_cells_iterator end() { return b->top_dimensional_cells_iterator_end(); }
 
   private:
    Bitmap_cubical_complex_base<T>* b;
  };
 
  Top_dimensional_cells_range top_dimensional_cells_range() { return Top_dimensional_cells_range(this); }
 
  //****************************************************************************************************************//
  //****************************************************************************************************************//
  //****************************************************************************************************************//
  //****************************************************************************************************************//
 
  inline std::size_t number_cells() const { return this->total_number_of_cells; }
 
  //****************************************************************************************************************//
  //****************************************************************************************************************//
  //****************************************************************************************************************//
  //****************************************************************************************************************//
 
 protected:
  std::vector<unsigned> sizes;
  std::vector<unsigned> multipliers;
  std::vector<T> data;
  std::size_t total_number_of_cells;
 
  void set_up_containers(const std::vector<unsigned>& sizes) {
    unsigned multiplier = 1;
    for (std::size_t i = 0; i != sizes.size(); ++i) {
      this->sizes.push_back(sizes[i]);
      this->multipliers.push_back(multiplier);
      multiplier *= 2 * sizes[i] + 1;
    }
    this->data = std::vector<T>(multiplier, std::numeric_limits<T>::infinity());
    this->total_number_of_cells = multiplier;
  }
 
  std::size_t compute_position_in_bitmap(const std::vector<unsigned>& counter) {
    std::size_t position = 0;
    for (std::size_t i = 0; i != this->multipliers.size(); ++i) {
      position += this->multipliers[i] * counter[i];
    }
    return position;
  }
 
  std::vector<unsigned> compute_counter_for_given_cell(std::size_t cell) const {
    std::vector<unsigned> counter;
    counter.reserve(this->sizes.size());
    for (std::size_t dim = this->sizes.size(); dim != 0; --dim) {
      counter.push_back(cell / this->multipliers[dim - 1]);
      cell = cell % this->multipliers[dim - 1];
    }
    std::reverse(counter.begin(), counter.end());
    return counter;
  }
  void read_perseus_style_file(const char* perseus_style_file);
  void setup_bitmap_based_on_top_dimensional_cells_list(const std::vector<unsigned>& sizes_in_following_directions,
                                                        const std::vector<T>& top_dimensional_cells);
  Bitmap_cubical_complex_base(const char* perseus_style_file, std::vector<bool> directions);
  Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes, std::vector<bool> directions);
  Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, const std::vector<T>& top_dimensional_cells,
                              std::vector<bool> directions);
};
 
template <typename T>
void Bitmap_cubical_complex_base<T>::put_data_to_bins(std::size_t number_of_bins) {
  bool dbg = false;
 
  std::pair<T, T> min_max = this->min_max_filtration();
  T dx = (min_max.second - min_max.first) / (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::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::clog << "After binning : " << this->data[i] << std::endl;
    }
  }
}
 
template <typename T>
void Bitmap_cubical_complex_base<T>::put_data_to_bins(T diameter_of_bin) {
  bool dbg = false;
  std::pair<T, T> min_max = this->min_max_filtration();
 
  std::size_t number_of_bins = (min_max.second - min_max.first) / 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::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::clog << "After binning : " << this->data[i] << std::endl;
    }
  }
}
 
template <typename T>
std::pair<T, T> Bitmap_cubical_complex_base<T>::min_max_filtration() {
  std::pair<T, T> min_max(std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity());
  for (std::size_t i = 0; i != this->data.size(); ++i) {
    if (this->data[i] < min_max.first) min_max.first = this->data[i];
    if (this->data[i] > min_max.second) min_max.second = this->data[i];
  }
  return min_max;
}
 
template <typename K>
std::ostream& operator<<(std::ostream& out, const Bitmap_cubical_complex_base<K>& b) {
  for (typename Bitmap_cubical_complex_base<K>::all_cells_const_iterator it = b.all_cells_const_begin();
       it != b.all_cells_const_end(); ++it) {
    out << *it << " ";
  }
  return out;
}
 
template <typename T>
Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes) {
  this->set_up_containers(sizes);
}
 
template <typename T>
void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells_list(
    const std::vector<unsigned>& sizes_in_following_directions, const std::vector<T>& top_dimensional_cells) {
  this->set_up_containers(sizes_in_following_directions);
 
  std::size_t number_of_top_dimensional_elements = 1;
  for (std::size_t i = 0; i != sizes_in_following_directions.size(); ++i) {
    number_of_top_dimensional_elements *= sizes_in_following_directions[i];
  }
  if (number_of_top_dimensional_elements != top_dimensional_cells.size()) {
    std::cerr << "Error in constructor Bitmap_cubical_complex_base ( std::vector<std::size_t> "
              << "sizes_in_following_directions, std::vector<T> top_dimensional_cells ). Number of top dimensional "
              << "elements that follow from sizes_in_following_directions vector is different than the size of "
              << "top_dimensional_cells vector."
              << std::endl;
    throw(
        "Error in constructor Bitmap_cubical_complex_base( std::vector<std::size_t> sizes_in_following_directions,"
        "std::vector<T> top_dimensional_cells ). Number of top dimensional elements that follow from "
        "sizes_in_following_directions vector is different than the size of top_dimensional_cells vector.");
  }
 
  Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);
  std::size_t index = 0;
  for (it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) {
    this->get_cell_data(*it) = top_dimensional_cells[index];
    ++index;
  }
  this->impose_lower_star_filtration();
}
 
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);
}
 
template <typename T>
void Bitmap_cubical_complex_base<T>::read_perseus_style_file(const char* perseus_style_file) {
  bool dbg = false;
  std::ifstream inFiltration;
  inFiltration.open(perseus_style_file);
  unsigned dimensionOfData;
  inFiltration >> dimensionOfData;
 
  if (dbg) {
    std::clog << "dimensionOfData : " << dimensionOfData << std::endl;
  }
 
  std::vector<unsigned> sizes;
  sizes.reserve(dimensionOfData);
  // all dimensions multiplied
  std::size_t dimensions = 1;
  for (std::size_t i = 0; i != dimensionOfData; ++i) {
    unsigned size_in_this_dimension;
    inFiltration >> size_in_this_dimension;
    sizes.push_back(size_in_this_dimension);
    dimensions *= size_in_this_dimension;
    if (dbg) {
      std::clog << "size_in_this_dimension : " << size_in_this_dimension << std::endl;
    }
  }
  this->set_up_containers(sizes);
 
  Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);
  it = this->top_dimensional_cells_iterator_begin();
 
  T filtrationLevel = 0.;
  std::size_t filtration_counter = 0;
  while (!inFiltration.eof()) {
    std::string line;
    getline(inFiltration, line);
    if (line.length() != 0) {
      int n = sscanf(line.c_str(), "%lf", &filtrationLevel);
      if (n != 1) {
        std::string perseus_error("Bad Perseus file format. This line is incorrect : " + line);
        throw std::ios_base::failure(perseus_error.c_str());
      }
 
      if (dbg) {
        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;
      }
      this->get_cell_data(*it) = filtrationLevel;
      ++it;
      ++filtration_counter;
    }
  }
 
  if (filtration_counter != dimensions) {
    std::string perseus_error("Bad Perseus file format. Read " + std::to_string(filtration_counter) + " expected " + \
      std::to_string(dimensions) + " values");
    throw std::ios_base::failure(perseus_error.c_str());
  }
 
  inFiltration.close();
  this->impose_lower_star_filtration();
}
 
template <typename T>
Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const char* perseus_style_file,
                                                            std::vector<bool> directions) {
  // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary
  // conditions.
  // It ignores the last parameter of the function.
  this->read_perseus_style_file(perseus_style_file);
}
 
template <typename T>
Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes,
                                                            std::vector<bool> directions) {
  // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary
  // conditions.
  // It ignores the last parameter of the function.
  this->set_up_containers(sizes);
}
 
template <typename T>
Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions,
                                                            const std::vector<T>& top_dimensional_cells,
                                                            std::vector<bool> directions) {
  // this constructor is here just for compatibility with a class that creates cubical complexes with periodic boundary
  // conditions.
  // It ignores the last parameter of the function.
  this->setup_bitmap_based_on_top_dimensional_cells_list(dimensions, top_dimensional_cells);
}
 
template <typename T>
Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const char* perseus_style_file) {
  this->read_perseus_style_file(perseus_style_file);
}
 
template <typename T>
std::vector<std::size_t> Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell(std::size_t cell) const {
  std::vector<std::size_t> boundary_elements;
 
  // Speed traded of for memory. Check if it is better in practice.
  boundary_elements.reserve(this->dimension() * 2);
 
  std::size_t sum_of_dimensions = 0;
  std::size_t cell1 = cell;
  for (std::size_t i = this->multipliers.size(); i != 0; --i) {
    unsigned position = cell1 / this->multipliers[i - 1];
    if (position % 2 == 1) {
      if (sum_of_dimensions % 2) {
        boundary_elements.push_back(cell + this->multipliers[i - 1]);
        boundary_elements.push_back(cell - this->multipliers[i - 1]);
      } else {
        boundary_elements.push_back(cell - this->multipliers[i - 1]);
        boundary_elements.push_back(cell + this->multipliers[i - 1]);
      }
      ++sum_of_dimensions;
    }
    cell1 = cell1 % this->multipliers[i - 1];
  }
 
  return boundary_elements;
}
 
template <typename T>
std::vector<std::size_t> Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(std::size_t cell) const {
  std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell);
  std::vector<std::size_t> coboundary_elements;
  std::size_t cell1 = cell;
  for (std::size_t i = this->multipliers.size(); i != 0; --i) {
    unsigned position = cell1 / this->multipliers[i - 1];
    if (position % 2 == 0) {
      if ((cell > this->multipliers[i - 1]) && (counter[i - 1] != 0)) {
        coboundary_elements.push_back(cell - this->multipliers[i - 1]);
      }
      if ((cell + this->multipliers[i - 1] < this->data.size()) && (counter[i - 1] != 2 * this->sizes[i - 1])) {
        coboundary_elements.push_back(cell + this->multipliers[i - 1]);
      }
    }
    cell1 = cell1 % this->multipliers[i - 1];
  }
  return coboundary_elements;
}
 
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::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::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::clog << "Nonzero length in this direction \n";
      dimension++;
    }
    cell = cell % this->multipliers[i - 1];
  }
  return dimension;
}
 
template <typename T>
inline T& Bitmap_cubical_complex_base<T>::get_cell_data(std::size_t cell) {
  return this->data[cell];
}
 
template <typename T>
void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration() {
  bool dbg = false;
 
  // this vector will be used to check which elements have already been taken care of in imposing lower star filtration
  std::vector<bool> is_this_cell_considered(this->data.size(), false);
 
  std::size_t size_to_reserve = 1;
  for (std::size_t i = 0; i != this->multipliers.size(); ++i) {
    size_to_reserve *= (std::size_t)((this->multipliers[i] - 1) / 2);
  }
 
  std::vector<std::size_t> indices_to_consider;
  indices_to_consider.reserve(size_to_reserve);
  // we assume here that we already have a filtration on the top dimensional cells and
  // we have to extend it to lower ones.
  typename Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this);
  for (it = this->top_dimensional_cells_iterator_begin(); it != this->top_dimensional_cells_iterator_end(); ++it) {
    indices_to_consider.push_back(it.compute_index_in_bitmap());
  }
 
  while (indices_to_consider.size()) {
    if (dbg) {
      std::clog << "indices_to_consider in this iteration \n";
      for (std::size_t i = 0; i != indices_to_consider.size(); ++i) {
        std::clog << indices_to_consider[i] << "  ";
      }
    }
    std::vector<std::size_t> new_indices_to_consider;
    for (std::size_t i = 0; i != indices_to_consider.size(); ++i) {
      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::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::clog << "Setting the value of a cell : " << bd[boundaryIt]
                      << " to : " << this->data[indices_to_consider[i]] << std::endl;
          }
        }
        if (is_this_cell_considered[bd[boundaryIt]] == false) {
          new_indices_to_consider.push_back(bd[boundaryIt]);
          is_this_cell_considered[bd[boundaryIt]] = true;
        }
      }
    }
    indices_to_consider.swap(new_indices_to_consider);
  }
}
 
template <typename T>
bool compareFirstElementsOfTuples(const std::pair<std::pair<T, std::size_t>, char>& first,
                                  const std::pair<std::pair<T, std::size_t>, char>& second) {
  if (first.first.first < second.first.first) {
    return true;
  } else {
    if (first.first.first > second.first.first) {
      return false;
    }
    // in this case first.first.first == second.first.first, so we need to compare dimensions
    return first.second < second.second;
  }
}
 
}  // namespace cubical_complex
 
namespace Cubical_complex = cubical_complex;
 
}  // namespace Gudhi
 
#endif  // BITMAP_CUBICAL_COMPLEX_BASE_H_