Hasse_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):       Clément Maria
 *
 *    Copyright (C) 2014 Inria
 *
 *    Modification(s):
 *      - YYYY/MM Author: Description of the modification
 */
 
#ifndef HASSE_COMPLEX_H_
#define HASSE_COMPLEX_H_
 
#include <gudhi/allocator.h>
 
#include <boost/iterator/counting_iterator.hpp>
 
#include <algorithm>
#include <utility>  // for pair
#include <vector>
#include <limits>  // for infinity value
 
#ifdef GUDHI_USE_TBB
#include <tbb/parallel_for.h>
#endif
 
namespace Gudhi {
 
template < class HasseCpx >
struct Hasse_simplex {
  // Complex_ds must verify that cpx->key(sh) is the order of sh in the filtration
 
  template< class Complex_ds >
  Hasse_simplex(Complex_ds & cpx
                , typename Complex_ds::Simplex_handle sh)
      : filtration_(cpx.filtration(sh))
      , boundary_() {
    boundary_.reserve(cpx.dimension(sh) + 1);
    for (auto b_sh : cpx.boundary_simplex_range(sh)) {
      boundary_.push_back(cpx.key(b_sh));
    }
  }
 
  Hasse_simplex(typename HasseCpx::Simplex_key key
                , typename HasseCpx::Filtration_value fil
                , std::vector<typename HasseCpx::Simplex_handle> const& boundary)
      : key_(key)
      , filtration_(fil)
      , boundary_(boundary) { }
 
  typename HasseCpx::Simplex_key key_;
  typename HasseCpx::Filtration_value filtration_;
  std::vector<typename HasseCpx::Simplex_handle> boundary_;
};
 
template < typename FiltrationValue = double
, typename SimplexKey = int
, typename VertexHandle = int
>
class Hasse_complex {
 public:
  typedef Hasse_simplex<Hasse_complex> Hasse_simp;
  typedef FiltrationValue Filtration_value;
  typedef SimplexKey Simplex_key;
  typedef int Simplex_handle;  // index in vector complex_
 
  typedef boost::counting_iterator< Simplex_handle > Filtration_simplex_iterator;
  typedef boost::iterator_range<Filtration_simplex_iterator> Filtration_simplex_range;
 
  typedef typename std::vector< Simplex_handle >::iterator Boundary_simplex_iterator;
  typedef boost::iterator_range<Boundary_simplex_iterator> Boundary_simplex_range;
 
  typedef typename std::vector< Simplex_handle >::iterator Skeleton_simplex_iterator;
  typedef boost::iterator_range< Skeleton_simplex_iterator > Skeleton_simplex_range;
 
  /*  only dimension 0 skeleton_simplex_range(...) */
  Skeleton_simplex_range skeleton_simplex_range(int dim = 0) {
    if (dim != 0) {
      std::cerr << "Dimension must be 0 \n";
    }
    return Skeleton_simplex_range(vertices_.begin(), vertices_.end());
  }
 
  template < class Complex_ds >
  Hasse_complex(Complex_ds & cpx)
      : complex_(cpx.num_simplices())
      , vertices_()
      , num_vertices_()
      , dim_max_(cpx.dimension()) {
    int size = complex_.size();
#ifdef GUDHI_USE_TBB
    tbb::parallel_for(0, size, [&](int idx){new (&complex_[idx]) Hasse_simp(cpx, cpx.simplex(idx));});
    for (int idx=0; idx < size; ++idx)
      if (complex_[idx].boundary_.empty())
        vertices_.push_back(idx);
#else
    for (int idx=0; idx < size; ++idx) {
      new (&complex_[idx]) Hasse_simp(cpx, cpx.simplex(idx));
      if (complex_[idx].boundary_.empty())
        vertices_.push_back(idx);
    }
#endif
  }
 
  Hasse_complex()
      : complex_()
      , vertices_()
      , num_vertices_(0)
      , dim_max_(-1) { }
 
  size_t num_simplices() {
    return complex_.size();
  }
 
  Filtration_simplex_range filtration_simplex_range() {
    return Filtration_simplex_range(Filtration_simplex_iterator(0)
                                    , Filtration_simplex_iterator(complex_.size()));
  }
 
  Simplex_key key(Simplex_handle sh) {
    return complex_[sh].key_;
  }
 
  Simplex_key null_key() {
    return -1;
  }
 
  Simplex_handle simplex(Simplex_key key) {
    if (key == null_key()) return null_simplex();
    return key;
  }
 
  Simplex_handle null_simplex() {
    return -1;
  }
 
  Filtration_value filtration(Simplex_handle sh) {
    if (sh == null_simplex()) {
      return std::numeric_limits<Filtration_value>::infinity();
    }
    return complex_[sh].filtration_;
  }
 
  int dimension(Simplex_handle sh) {
    if (complex_[sh].boundary_.empty()) return 0;
    return complex_[sh].boundary_.size() - 1;
  }
 
  int dimension() {
    return dim_max_;
  }
 
  std::pair<Simplex_handle, Simplex_handle> endpoints(Simplex_handle sh) {
    return std::pair<Simplex_handle, Simplex_handle>(complex_[sh].boundary_[0]
                                                     , complex_[sh].boundary_[1]);
  }
 
  void assign_key(Simplex_handle sh, Simplex_key key) {
    complex_[sh].key_ = key;
  }
 
  Boundary_simplex_range boundary_simplex_range(Simplex_handle sh) {
    return Boundary_simplex_range(complex_[sh].boundary_.begin()
                                  , complex_[sh].boundary_.end());
  }
 
  void display_simplex(Simplex_handle sh) {
    std::clog << dimension(sh) << "  ";
    for (auto sh_b : boundary_simplex_range(sh)) std::clog << sh_b << " ";
    std::clog << "  " << filtration(sh) << "         key=" << key(sh);
  }
 
  void initialize_filtration() {
    // Setting the keys is done by pcoh, Simplex_tree doesn't do it either.
#if 0
    Simplex_key key = 0;
    for (auto & h_simp : complex_)
      h_simp.key_ = key++;
#endif
  }
 
  std::vector< Hasse_simp, Gudhi::no_init_allocator<Hasse_simp> > complex_;
  std::vector<Simplex_handle> vertices_;
  size_t num_vertices_;
  int dim_max_;
};
 
template< typename T1, typename T2, typename T3 >
std::istream& operator>>(std::istream & is
                         , Hasse_complex< T1, T2, T3 > & hcpx) {
  assert(hcpx.num_simplices() == 0);
 
  size_t num_simp;
  is >> num_simp;
  hcpx.complex_.reserve(num_simp);
 
  std::vector< typename Hasse_complex<T1, T2, T3>::Simplex_key > boundary;
  typename Hasse_complex<T1, T2, T3>::Filtration_value fil;
  typename Hasse_complex<T1, T2, T3>::Filtration_value max_fil = 0;
  int max_dim = -1;
  int key = 0;
  // read all simplices in the file as a list of vertices
  while (read_hasse_simplex(is, boundary, fil)) {
    // insert every simplex in the simplex tree
    hcpx.complex_.emplace_back(key, fil, boundary);
 
    if (max_dim < hcpx.dimension(key)) {
      max_dim = hcpx.dimension(key);
    }
    if (hcpx.dimension(key) == 0) {
      hcpx.vertices_.push_back(key);
    }
    if (max_fil < fil) {
      max_fil = fil;
    }
 
    ++key;
    boundary.clear();
  }
 
  hcpx.dim_max_ = max_dim;
 
  return is;
}
 
}  // namespace Gudhi
 
#endif  // HASSE_COMPLEX_H_