Bitmap_cubical_complex.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_H_
#define BITMAP_CUBICAL_COMPLEX_H_

#include <gudhi/Bitmap_cubical_complex_base.h>
#include <gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h>

#ifdef GUDHI_USE_TBB
#include <tbb/parallel_sort.h>
#endif

#include <limits>
#include <utility>    // for pair<>
#include <algorithm>  // for sort
#include <vector>
#include <numeric>  // for iota
#include <cstddef>

namespace Gudhi {

namespace cubical_complex {

// global variable, was used just for debugging.
const bool globalDbg = false;

template <typename T>
class is_before_in_filtration;

template <typename T>
class Bitmap_cubical_complex : public T {
 public:
  //*********************************************//
  // Typedefs and typenames
  //*********************************************//
  typedef std::size_t Simplex_key;
  typedef typename T::filtration_type Filtration_value;
  typedef Simplex_key Simplex_handle;

  //*********************************************//
  // Constructors
  //*********************************************//
  // Over here we need to define various input types. I am proposing the following ones:
  // Perseus style
  // TODO(PD) H5 files?
  // TODO(PD) binary files with little endiangs / big endians ?
  // TODO(PD) constructor from a vector of elements of a type T. ?

  Bitmap_cubical_complex(const char* perseus_style_file)
      : T(perseus_style_file), key_associated_to_simplex(this->total_number_of_cells + 1) {
    if (globalDbg) {
      std::clog << "Bitmap_cubical_complex( const char* perseus_style_file )\n";
    }
    for (std::size_t i = 0; i != this->total_number_of_cells; ++i) {
      this->key_associated_to_simplex[i] = i;
    }
    // we initialize this only once, in each constructor, when the bitmap is constructed.
    // If the user decide to change some elements of the bitmap, then this procedure need
    // to be called again.
    this->initialize_simplex_associated_to_key();
  }

  Bitmap_cubical_complex(const std::vector<unsigned>& dimensions,
                         const std::vector<Filtration_value>& top_dimensional_cells)
      : T(dimensions, top_dimensional_cells), key_associated_to_simplex(this->total_number_of_cells + 1) {
    for (std::size_t i = 0; i != this->total_number_of_cells; ++i) {
      this->key_associated_to_simplex[i] = i;
    }
    // we initialize this only once, in each constructor, when the bitmap is constructed.
    // If the user decide to change some elements of the bitmap, then this procedure need
    // to be called again.
    this->initialize_simplex_associated_to_key();
  }

  Bitmap_cubical_complex(const std::vector<unsigned>& dimensions,
                         const std::vector<Filtration_value>& top_dimensional_cells,
                         std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed)
      : T(dimensions, top_dimensional_cells, directions_in_which_periodic_b_cond_are_to_be_imposed),
        key_associated_to_simplex(this->total_number_of_cells + 1) {
    for (std::size_t i = 0; i != this->total_number_of_cells; ++i) {
      this->key_associated_to_simplex[i] = i;
    }
    // we initialize this only once, in each constructor, when the bitmap is constructed.
    // If the user decide to change some elements of the bitmap, then this procedure need
    // to be called again.
    this->initialize_simplex_associated_to_key();
  }

  virtual ~Bitmap_cubical_complex() {}

  //*********************************************//
  // Other 'easy' functions
  //*********************************************//

  std::size_t num_simplices() const { return this->total_number_of_cells; }

  static Simplex_handle null_simplex() {
    if (globalDbg) {
      std::clog << "Simplex_handle null_simplex()\n";
    }
    return std::numeric_limits<Simplex_handle>::max();
  }

  inline std::size_t dimension() const { return this->sizes.size(); }

  inline unsigned dimension(Simplex_handle sh) const {
    if (globalDbg) {
      std::clog << "unsigned dimension(const Simplex_handle& sh)\n";
    }
    if (sh != null_simplex()) return this->get_dimension_of_a_cell(sh);
    return -1;
  }

  Filtration_value filtration(Simplex_handle sh) {
    if (globalDbg) {
      std::clog << "Filtration_value filtration(const Simplex_handle& sh)\n";
    }
    // Returns the filtration value of a simplex.
    if (sh != null_simplex()) return this->data[sh];
    return std::numeric_limits<Filtration_value>::infinity();
  }

  static Simplex_key null_key() {
    if (globalDbg) {
      std::clog << "Simplex_key null_key()\n";
    }
    return std::numeric_limits<Simplex_handle>::max();
  }

  Simplex_key key(Simplex_handle sh) const {
    if (globalDbg) {
      std::clog << "Simplex_key key(const Simplex_handle& sh)\n";
    }
    if (sh != null_simplex()) {
      return this->key_associated_to_simplex[sh];
    }
    return this->null_key();
  }

  Simplex_handle simplex(Simplex_key key) {
    if (globalDbg) {
      std::clog << "Simplex_handle simplex(Simplex_key key)\n";
    }
    if (key != null_key()) {
      return this->simplex_associated_to_key[key];
    }
    return null_simplex();
  }

  void assign_key(Simplex_handle sh, Simplex_key key) {
    if (globalDbg) {
      std::clog << "void assign_key(Simplex_handle& sh, Simplex_key key)\n";
    }
    if (key == null_key()) return;
    this->key_associated_to_simplex[sh] = key;
    this->simplex_associated_to_key[key] = sh;
  }

  void initialize_simplex_associated_to_key();

  //*********************************************//
  // Iterators
  //*********************************************//

  typedef typename std::vector<Simplex_handle>::iterator Boundary_simplex_iterator;
  typedef typename std::vector<Simplex_handle> Boundary_simplex_range;

  class Filtration_simplex_range;

  class Filtration_simplex_iterator : std::iterator<std::input_iterator_tag, Simplex_handle> {
    // Iterator over all simplices of the complex in the order of the indexing scheme.
    // 'value_type' must be 'Simplex_handle'.
   public:
    Filtration_simplex_iterator(Bitmap_cubical_complex* b) : b(b), position(0) {}

    Filtration_simplex_iterator() : b(NULL), position(0) {}

    Filtration_simplex_iterator operator++() {
      if (globalDbg) {
        std::clog << "Filtration_simplex_iterator operator++\n";
      }
      ++this->position;
      return (*this);
    }

    Filtration_simplex_iterator operator++(int) {
      Filtration_simplex_iterator result = *this;
      ++(*this);
      return result;
    }

    Filtration_simplex_iterator& operator=(const Filtration_simplex_iterator& rhs) {
      if (globalDbg) {
        std::clog << "Filtration_simplex_iterator operator =\n";
      }
      this->b = rhs.b;
      this->position = rhs.position;
      return (*this);
    }

    bool operator==(const Filtration_simplex_iterator& rhs) const {
      if (globalDbg) {
        std::clog << "bool operator == ( const Filtration_simplex_iterator& rhs )\n";
      }
      return (this->position == rhs.position);
    }

    bool operator!=(const Filtration_simplex_iterator& rhs) const {
      if (globalDbg) {
        std::clog << "bool operator != ( const Filtration_simplex_iterator& rhs )\n";
      }
      return !(*this == rhs);
    }

    Simplex_handle operator*() {
      if (globalDbg) {
        std::clog << "Simplex_handle operator*()\n";
      }
      return this->b->simplex_associated_to_key[this->position];
    }

    friend class Filtration_simplex_range;

   private:
    Bitmap_cubical_complex<T>* b;
    std::size_t position;
  };

  class Filtration_simplex_range {
    // Range over the simplices of the complex in the order of the filtration.
    // .begin() and .end() return type Filtration_simplex_iterator.
   public:
    typedef Filtration_simplex_iterator const_iterator;
    typedef Filtration_simplex_iterator iterator;

    Filtration_simplex_range(Bitmap_cubical_complex<T>* b) : b(b) {}

    Filtration_simplex_iterator begin() {
      if (globalDbg) {
        std::clog << "Filtration_simplex_iterator begin() \n";
      }
      return Filtration_simplex_iterator(this->b);
    }

    Filtration_simplex_iterator end() {
      if (globalDbg) {
        std::clog << "Filtration_simplex_iterator end()\n";
      }
      Filtration_simplex_iterator it(this->b);
      it.position = this->b->simplex_associated_to_key.size();
      return it;
    }

   private:
    Bitmap_cubical_complex<T>* b;
  };

  //*********************************************//
  // Methods to access iterators from the container:

  Boundary_simplex_range boundary_simplex_range(Simplex_handle sh) { return this->get_boundary_of_a_cell(sh); }

  Filtration_simplex_range filtration_simplex_range() {
    if (globalDbg) {
      std::clog << "Filtration_simplex_range filtration_simplex_range()\n";
    }
    // Returns a range over the simplices of the complex in the order of the filtration
    return Filtration_simplex_range(this);
  }
  //*********************************************//

  //*********************************************//
  // Elements which are in Gudhi now, but I (and in all the cases I asked also Marc) do not understand why they are
  // there.
  // TODO(PD) the file IndexingTag.h in the Gudhi library contains an empty structure, so
  // I understand that this is something that was planned (for simplicial maps?)
  // but was never finished. The only idea I have here is to use the same empty structure from
  // IndexingTag.h file, but only if the compiler needs it. If the compiler
  // do not need it, then I would rather not add here elements which I do not understand.
  // typedef Indexing_tag

  std::pair<Simplex_handle, Simplex_handle> endpoints(Simplex_handle sh) {
    std::vector<std::size_t> bdry = this->get_boundary_of_a_cell(sh);
    if (globalDbg) {
      std::clog << "std::pair<Simplex_handle, Simplex_handle> endpoints( Simplex_handle sh )\n";
      std::clog << "bdry.size() : " << bdry.size() << "\n";
    }
    // this method returns two first elements from the boundary of sh.
    if (bdry.size() < 2)
      throw(
          "Error in endpoints in Bitmap_cubical_complex class. The cell have less than two elements in the "
          "boundary.");
    return std::make_pair(bdry[0], bdry[1]);
  }

  class Skeleton_simplex_range;

  class Skeleton_simplex_iterator : std::iterator<std::input_iterator_tag, Simplex_handle> {
    // Iterator over all simplices of the complex in the order of the indexing scheme.
    // 'value_type' must be 'Simplex_handle'.
   public:
    Skeleton_simplex_iterator(Bitmap_cubical_complex* b, std::size_t d) : b(b), dimension(d) {
      if (globalDbg) {
        std::clog << "Skeleton_simplex_iterator ( Bitmap_cubical_complex* b , std::size_t d )\n";
      }
      // find the position of the first simplex of a dimension d
      this->position = 0;
      while ((this->position != b->data.size()) &&
             (this->b->get_dimension_of_a_cell(this->position) != this->dimension)) {
        ++this->position;
      }
    }

    Skeleton_simplex_iterator() : b(NULL), position(0), dimension(0) {}

    Skeleton_simplex_iterator operator++() {
      if (globalDbg) {
        std::clog << "Skeleton_simplex_iterator operator++()\n";
      }
      // increment the position as long as you did not get to the next element of the dimension dimension.
      ++this->position;
      while ((this->position != this->b->data.size()) &&
             (this->b->get_dimension_of_a_cell(this->position) != this->dimension)) {
        ++this->position;
      }
      return (*this);
    }

    Skeleton_simplex_iterator operator++(int) {
      Skeleton_simplex_iterator result = *this;
      ++(*this);
      return result;
    }

    Skeleton_simplex_iterator& operator=(const Skeleton_simplex_iterator& rhs) {
      if (globalDbg) {
        std::clog << "Skeleton_simplex_iterator operator =\n";
      }
      this->b = rhs.b;
      this->position = rhs.position;
      this->dimension = rhs.dimension;
      return (*this);
    }

    bool operator==(const Skeleton_simplex_iterator& rhs) const {
      if (globalDbg) {
        std::clog << "bool operator ==\n";
      }
      return (this->position == rhs.position);
    }

    bool operator!=(const Skeleton_simplex_iterator& rhs) const {
      if (globalDbg) {
        std::clog << "bool operator != ( const Skeleton_simplex_iterator& rhs )\n";
      }
      return !(*this == rhs);
    }

    Simplex_handle operator*() {
      if (globalDbg) {
        std::clog << "Simplex_handle operator*() \n";
      }
      return this->position;
    }

    friend class Skeleton_simplex_range;

   private:
    Bitmap_cubical_complex<T>* b;
    std::size_t position;
    unsigned dimension;
  };

  class Skeleton_simplex_range {
    // Range over the simplices of the complex in the order of the filtration.
    // .begin() and .end() return type Filtration_simplex_iterator.
   public:
    typedef Skeleton_simplex_iterator const_iterator;
    typedef Skeleton_simplex_iterator iterator;

    Skeleton_simplex_range(Bitmap_cubical_complex<T>* b, unsigned dimension) : b(b), dimension(dimension) {}

    Skeleton_simplex_iterator begin() {
      if (globalDbg) {
        std::clog << "Skeleton_simplex_iterator begin()\n";
      }
      return Skeleton_simplex_iterator(this->b, this->dimension);
    }

    Skeleton_simplex_iterator end() {
      if (globalDbg) {
        std::clog << "Skeleton_simplex_iterator end()\n";
      }
      Skeleton_simplex_iterator it(this->b, this->dimension);
      it.position = this->b->data.size();
      return it;
    }

   private:
    Bitmap_cubical_complex<T>* b;
    unsigned dimension;
  };

  Skeleton_simplex_range skeleton_simplex_range(unsigned dimension) {
    if (globalDbg) {
      std::clog << "Skeleton_simplex_range skeleton_simplex_range( unsigned dimension )\n";
    }
    return Skeleton_simplex_range(this, dimension);
  }

  friend class is_before_in_filtration<T>;

 protected:
  std::vector<std::size_t> key_associated_to_simplex;
  std::vector<std::size_t> simplex_associated_to_key;
};  // Bitmap_cubical_complex

template <typename T>
void Bitmap_cubical_complex<T>::initialize_simplex_associated_to_key() {
  if (globalDbg) {
    std::clog << "void Bitmap_cubical_complex<T>::initialize_elements_ordered_according_to_filtration() \n";
  }
  this->simplex_associated_to_key = std::vector<std::size_t>(this->data.size());
  std::iota(std::begin(simplex_associated_to_key), std::end(simplex_associated_to_key), 0);
#ifdef GUDHI_USE_TBB
  tbb::parallel_sort(simplex_associated_to_key.begin(), simplex_associated_to_key.end(),
                     is_before_in_filtration<T>(this));
#else
  std::sort(simplex_associated_to_key.begin(), simplex_associated_to_key.end(), is_before_in_filtration<T>(this));
#endif

  // we still need to deal here with a key_associated_to_simplex:
  for (std::size_t i = 0; i != simplex_associated_to_key.size(); ++i) {
    this->key_associated_to_simplex[simplex_associated_to_key[i]] = i;
  }
}

template <typename T>
class is_before_in_filtration {
 public:
  explicit is_before_in_filtration(Bitmap_cubical_complex<T>* CC) : CC_(CC) {}

  bool operator()(const typename Bitmap_cubical_complex<T>::Simplex_handle& sh1,
                  const typename Bitmap_cubical_complex<T>::Simplex_handle& sh2) const {
    // Not using st_->filtration(sh1) because it uselessly tests for null_simplex.
    typedef typename T::filtration_type Filtration_value;
    Filtration_value fil1 = CC_->data[sh1];
    Filtration_value fil2 = CC_->data[sh2];
    if (fil1 != fil2) {
      return fil1 < fil2;
    }
    // in this case they are on the same filtration level, so the dimension decide.
    std::size_t dim1 = CC_->get_dimension_of_a_cell(sh1);
    std::size_t dim2 = CC_->get_dimension_of_a_cell(sh2);
    if (dim1 != dim2) {
      return dim1 < dim2;
    }
    // in this case both filtration and dimensions of the considered cubes are the same. To have stable sort, we simply
    // compare their positions in the bitmap:
    return sh1 < sh2;
  }

 protected:
  Bitmap_cubical_complex<T>* CC_;
};

}  // namespace cubical_complex

namespace Cubical_complex = cubical_complex;

}  // namespace Gudhi

#endif  // BITMAP_CUBICAL_COMPLEX_H_