Simplex_tree.h

/*    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):       Clément Maria
 *
 *    Copyright (C) 2014  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 SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_H_
#define SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_H_

#include <gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h>
#include <gudhi/Simplex_tree/Simplex_tree_siblings.h>
#include <gudhi/Simplex_tree/Simplex_tree_iterators.h>
#include <gudhi/Simplex_tree/indexing_tag.h>

#include <boost/container/flat_map.hpp>
#include <boost/iterator/transform_iterator.hpp>
#include <boost/graph/adjacency_list.hpp>

#include <algorithm>
#include <utility>
#include <vector>

namespace Gudhi {
  
template<typename IndexingTag = linear_indexing_tag,
    typename FiltrationValue = double, typename SimplexKey = int  // must be a signed integer type
    , typename VertexHandle = int  // must be a signed integer type, int convertible to it
//         , bool ContiguousVertexHandles = true   //true is Vertex_handles are exactly the set [0;n)
>
class Simplex_tree {
 public:
  typedef IndexingTag Indexing_tag;
  typedef FiltrationValue Filtration_value;
  typedef SimplexKey Simplex_key;
  typedef VertexHandle Vertex_handle;

  /* Type of node in the simplex tree. */
  typedef Simplex_tree_node_explicit_storage<Simplex_tree> Node;
  /* Type of dictionary Vertex_handle -> Node for traversing the simplex tree. */
  typedef typename boost::container::flat_map<Vertex_handle, Node> Dictionary;

  friend class Simplex_tree_node_explicit_storage< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >;
  friend class Simplex_tree_siblings< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle>, Dictionary>;
  friend class Simplex_tree_simplex_vertex_iterator< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >;
  friend class Simplex_tree_boundary_simplex_iterator< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >;
  friend class Simplex_tree_complex_simplex_iterator< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >;
  friend class Simplex_tree_skeleton_simplex_iterator< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >;

  /* \brief Set of nodes sharing a same parent in the simplex tree. */
  /* \brief Set of nodes sharing a same parent in the simplex tree. */
  typedef Simplex_tree_siblings<Simplex_tree, Dictionary> Siblings;

 public:
  typedef typename Dictionary::iterator Simplex_handle;

 private:
  typedef typename Dictionary::iterator Dictionary_it;
  typedef typename Dictionary_it::value_type Dit_value_t;

  struct return_first {
    Vertex_handle operator()(const Dit_value_t& p_sh) const {
      return p_sh.first;
    }
  };

 public:
  typedef boost::transform_iterator<return_first, Dictionary_it> Complex_vertex_iterator;
  typedef boost::iterator_range<Complex_vertex_iterator> Complex_vertex_range;
  typedef Simplex_tree_simplex_vertex_iterator<Simplex_tree> Simplex_vertex_iterator;
  typedef boost::iterator_range<Simplex_vertex_iterator> Simplex_vertex_range;
  typedef Simplex_tree_boundary_simplex_iterator<Simplex_tree> Boundary_simplex_iterator;
  typedef boost::iterator_range<Boundary_simplex_iterator> Boundary_simplex_range;
  typedef Simplex_tree_complex_simplex_iterator<Simplex_tree> Complex_simplex_iterator;
  typedef boost::iterator_range<Complex_simplex_iterator> Complex_simplex_range;
  typedef Simplex_tree_skeleton_simplex_iterator<Simplex_tree> Skeleton_simplex_iterator;
  typedef boost::iterator_range<Skeleton_simplex_iterator> Skeleton_simplex_range;
  typedef typename std::vector<Simplex_handle>::iterator Filtration_simplex_iterator;
  typedef boost::iterator_range<Filtration_simplex_iterator> Filtration_simplex_range;

  /* @} */  // end name range and iterator types
  Complex_vertex_range complex_vertex_range() {
    return Complex_vertex_range(
        boost::make_transform_iterator(root_.members_.begin(), return_first()),
        boost::make_transform_iterator(root_.members_.end(), return_first()));
  }

  Complex_simplex_range complex_simplex_range() {
    return Complex_simplex_range(Complex_simplex_iterator(this),
                                 Complex_simplex_iterator());
  }

  Skeleton_simplex_range skeleton_simplex_range(int dim) {
    return Skeleton_simplex_range(Skeleton_simplex_iterator(this, dim),
                                  Skeleton_simplex_iterator());
  }

  Filtration_simplex_range filtration_simplex_range(linear_indexing_tag) {
    if (filtration_vect_.empty()) {
      initialize_filtration();
    }
    return Filtration_simplex_range(filtration_vect_.begin(),
                                    filtration_vect_.end());
  }

  Filtration_simplex_range filtration_simplex_range() {
    return filtration_simplex_range(Indexing_tag());
  }
  Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) {
    return Simplex_vertex_range(Simplex_vertex_iterator(this, sh),
                                Simplex_vertex_iterator(this));
  }

  Boundary_simplex_range boundary_simplex_range(Simplex_handle sh) {
    return Boundary_simplex_range(Boundary_simplex_iterator(this, sh),
                                  Boundary_simplex_iterator(this));
  }
  // end range and iterator methods
  Simplex_tree()
      : null_vertex_(-1),
        threshold_(0),
        num_simplices_(0),
        root_(NULL, null_vertex_),
        filtration_vect_(),
        dimension_(-1) {
  }

  ~Simplex_tree() {
    for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh) {
      if (has_children(sh)) {
        rec_delete(sh->second.children());
      }
    }
  }  // end constructor/destructor
 private:
  void rec_delete(Siblings * sib) {
    for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) {
      if (has_children(sh)) {
        rec_delete(sh->second.children());
      }
    }
    delete sib;
  }

 public:
  Simplex_key key(Simplex_handle sh) {
    return sh->second.key();
  }
  Simplex_handle simplex(Simplex_key key) {
    return filtration_vect_[key];
  }
  Filtration_value filtration(Simplex_handle sh) {
    if (sh != null_simplex()) {
      return sh->second.filtration();
    } else {
      return INFINITY;
    }  // filtration(); }
  }
  Filtration_value filtration() {
    return threshold_;
  }
  Simplex_handle null_simplex() {
    return Dictionary_it(NULL);
  }
  Simplex_key null_key() {
    return -1;
  }
  Vertex_handle null_vertex() {
    return null_vertex_;
  }
  size_t num_vertices() {
    return root_.members_.size();
  }
  const unsigned int& num_simplices() const {
    return num_simplices_;
  }

  int dimension(Simplex_handle sh) {
    Siblings * curr_sib = self_siblings(sh);
    int dim = 0;
    while (curr_sib != NULL) {
      ++dim;
      curr_sib = curr_sib->oncles();
    }
    return dim - 1;
  }
  int dimension() {
    return dimension_;
  }

  bool has_children(Simplex_handle sh) {
    return (sh->second.children()->parent() == sh->first);
  }

  template<class RandomAccessVertexRange>
  Simplex_handle find(const RandomAccessVertexRange & s) {
    if (s.begin() == s.end())
      std::cerr << "Empty simplex \n";

    sort(s.begin(), s.end());

    Siblings * tmp_sib = &root_;
    Dictionary_it tmp_dit;
    Vertex_handle last = s[s.size() - 1];
    for (auto v : s) {
      tmp_dit = tmp_sib->members_.find(v);
      if (tmp_dit == tmp_sib->members_.end()) {
        return null_simplex();
      }
      if (!has_children(tmp_dit) && v != last) {
        return null_simplex();
      }
      tmp_sib = tmp_dit->second.children();
    }
    return tmp_dit;
  }

  Simplex_handle find_vertex(Vertex_handle v) {
    return root_.members_.begin() + v;
  }
//{ return root_.members_.find(v); }

  template<class RandomAccessVertexRange>
  std::pair<Simplex_handle, bool> insert(RandomAccessVertexRange & simplex,
                                         Filtration_value filtration) {
    if (simplex.empty()) {
      return std::pair<Simplex_handle, bool>(null_simplex(), true);
    }

    sort(simplex.begin(), simplex.end());  // must be sorted in increasing order

    Siblings * curr_sib = &root_;
    std::pair<Simplex_handle, bool> res_insert;
    typename RandomAccessVertexRange::iterator vi;
    for (vi = simplex.begin(); vi != simplex.end() - 1; ++vi) {
      res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration));
      if (!(has_children(res_insert.first))) {
        res_insert.first->second.assign_children(new Siblings(curr_sib, *vi));
      }
      curr_sib = res_insert.first->second.children();
    }
    res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration));
    if (!res_insert.second) {  // if already in the complex
      if (res_insert.first->second.filtration() > filtration) {  // if filtration value modified
        res_insert.first->second.assign_filtration(filtration);
        return res_insert;
      }
      return std::pair<Simplex_handle, bool>(null_simplex(), false);  // if filtration value unchanged
    }
    // otherwise the insertion has succeeded
    return res_insert;
  }

  void assign_key(Simplex_handle sh, Simplex_key key) {
    sh->second.assign_key(key);
  }

 public:
  std::pair<Simplex_handle, Simplex_handle> endpoints(Simplex_handle sh) {
    return std::pair<Simplex_handle, Simplex_handle>(
        root_.members_.find(sh->first),
        root_.members_.find(self_siblings(sh)->parent()));
  }

  Siblings * self_siblings(Simplex_handle sh) {
    if (sh->second.children()->parent() == sh->first)
      return sh->second.children()->oncles();
    else
      return sh->second.children();
  }

// void display_simplex(Simplex_handle sh)
// {
//   std::cout << "   " << "[" << filtration(sh) << "] ";
//   for( auto vertex : simplex_vertex_range(sh) )
//     { std::cout << vertex << " "; }
// }

  // void print(Simplex_handle sh, std::ostream& os = std::cout)
  // {  for(auto v : simplex_vertex_range(sh)) {os << v << " ";}
  //    os << std::endl; }

 public:
  Siblings * root() {
    return &root_;
  }

 public:
  void set_filtration(Filtration_value fil) {
    threshold_ = fil;
  }
  void set_num_simplices(const unsigned int& num_simplices) {
    num_simplices_ = num_simplices;
  }
  void set_dimension(int dimension) {
    dimension_ = dimension;
  }

 public:
  void initialize_filtration() {
    filtration_vect_.clear();
    filtration_vect_.reserve(num_simplices());
    for (auto cpx_it = complex_simplex_range().begin();
        cpx_it != complex_simplex_range().end(); ++cpx_it) {
      filtration_vect_.push_back(*cpx_it);
    }

    // the stable sort ensures the ordering heuristic
    std::stable_sort(filtration_vect_.begin(), filtration_vect_.end(),
                     is_before_in_filtration(this));
  }

 private:
  bool reverse_lexicographic_order(Simplex_handle sh1, Simplex_handle sh2) {
    Simplex_vertex_range rg1 = simplex_vertex_range(sh1);
    Simplex_vertex_range rg2 = simplex_vertex_range(sh2);
    Simplex_vertex_iterator it1 = rg1.begin();
    Simplex_vertex_iterator it2 = rg2.begin();
    while (it1 != rg1.end() && it2 != rg2.end()) {
      if (*it1 == *it2) {
        ++it1;
        ++it2;
      } else {
        return *it1 < *it2;
      }
    }
    return ((it1 == rg1.end()) && (it2 != rg2.end()));
  }
  struct is_before_in_filtration {
    explicit is_before_in_filtration(Simplex_tree * st)
        : st_(st) {
    }

    bool operator()(const Simplex_handle sh1, const Simplex_handle sh2) const {
      if (st_->filtration(sh1) != st_->filtration(sh2)) {
        return st_->filtration(sh1) < st_->filtration(sh2);
      }

      return st_->reverse_lexicographic_order(sh1, sh2);  // is sh1 a proper subface of sh2
    }

    Simplex_tree * st_;
  };

 public:
  template<class OneSkeletonGraph>
  void insert_graph(const OneSkeletonGraph& skel_graph) {
    assert(num_simplices() == 0);  // the simplex tree must be empty

    if (boost::num_vertices(skel_graph) == 0) {
      return;
    }
    if (num_edges(skel_graph) == 0) {
      dimension_ = 0;
    } else {
      dimension_ = 1;
    }

    num_simplices_ = boost::num_vertices(skel_graph)
        + boost::num_edges(skel_graph);
    root_.members_.reserve(boost::num_vertices(skel_graph));

    typename boost::graph_traits<OneSkeletonGraph>::vertex_iterator v_it,
        v_it_end;
    for (tie(v_it, v_it_end) = boost::vertices(skel_graph); v_it != v_it_end;
        ++v_it) {
      root_.members_.emplace_hint(
          root_.members_.end(), *v_it,
          Node(&root_, boost::get(vertex_filtration_t(), skel_graph, *v_it)));
    }
    typename boost::graph_traits<OneSkeletonGraph>::edge_iterator e_it,
        e_it_end;
    for (tie(e_it, e_it_end) = boost::edges(skel_graph); e_it != e_it_end;
        ++e_it) {
      auto u = source(*e_it, skel_graph);
      auto v = target(*e_it, skel_graph);
      if (u < v) {  // count edges only once { std::swap(u,v); } // u < v
        auto sh = find_vertex(u);
        if (!has_children(sh)) {
          sh->second.assign_children(new Siblings(&root_, sh->first));
        }

        sh->second.children()->members().emplace(
            v,
            Node(sh->second.children(),
                 boost::get(edge_filtration_t(), skel_graph, *e_it)));
      }
    }
  }
  void expansion(int max_dim) {
    dimension_ = max_dim;
    for (Dictionary_it root_it = root_.members_.begin();
        root_it != root_.members_.end(); ++root_it) {
      if (has_children(root_it)) {
        siblings_expansion(root_it->second.children(), max_dim - 1);
      }
    }
    dimension_ = max_dim - dimension_;
  }

 private:
  void siblings_expansion(Siblings * siblings,  // must contain elements
      int k) {
    if (dimension_ > k) {
      dimension_ = k;
    }
    if (k == 0)
      return;
    Dictionary_it next = siblings->members().begin();
    ++next;

    static std::vector<std::pair<Vertex_handle, Node> > inter;  // static, not thread-safe.
    for (Dictionary_it s_h = siblings->members().begin();
        s_h != siblings->members().end(); ++s_h, ++next) {
      Simplex_handle root_sh = find_vertex(s_h->first);
      if (has_children(root_sh)) {
        intersection(
            inter,                      // output intersection
            next,                       // begin
            siblings->members().end(),  // end
            root_sh->second.children()->members().begin(),
            root_sh->second.children()->members().end(),
            s_h->second.filtration());
        if (inter.size() != 0) {
          this->num_simplices_ += inter.size();
          Siblings * new_sib = new Siblings(siblings,  // oncles
              s_h->first,  // parent
              inter);      // boost::container::ordered_unique_range_t
          inter.clear();
          s_h->second.assign_children(new_sib);
          siblings_expansion(new_sib, k - 1);
        } else {
          s_h->second.assign_children(siblings);  // ensure the children property
          inter.clear();
        }
      }
    }
  }
  void intersection(std::vector<std::pair<Vertex_handle, Node> >& intersection,
                    Dictionary_it begin1, Dictionary_it end1,
                    Dictionary_it begin2, Dictionary_it end2,
                    Filtration_value filtration) {
    if (begin1 == end1 || begin2 == end2)
      return;  // ----->>
    while (true) {
      if (begin1->first == begin2->first) {
        intersection.push_back(
            std::pair<Vertex_handle, Node>(
                begin1->first,
                Node(NULL, maximum(begin1->second.filtration(), begin2->second.filtration(), filtration))));
        ++begin1;
        ++begin2;
        if (begin1 == end1 || begin2 == end2)
          return;  // ----->>
      } else {
        if (begin1->first < begin2->first) {
          ++begin1;
          if (begin1 == end1)
            return;
        } else {
          ++begin2;
          if (begin2 == end2)
            return;  // ----->>
        }
      }
    }
  }
  Filtration_value maximum(Filtration_value a, Filtration_value b,
                           Filtration_value c) {
    Filtration_value max = (a < b) ? b : a;
    return ((max < c) ? c : max);
  }

 public:
  void print_hasse(std::ostream& os) {
    os << num_simplices() << " " << std::endl;
    for (auto sh : filtration_simplex_range(Indexing_tag())) {
      os << dimension(sh) << " ";
      for (auto b_sh : boundary_simplex_range(sh)) {
        os << key(b_sh) << " ";
      }
      os << filtration(sh) << " \n";
    }
  }

 private:
  Vertex_handle null_vertex_;
  Filtration_value threshold_;
  unsigned int num_simplices_;
  Siblings root_;
  std::vector<Simplex_handle> filtration_vect_;
  int dimension_;
};

// Print a Simplex_tree in os.
template<typename T1, typename T2, typename T3>
std::ostream& operator<<(std::ostream & os, Simplex_tree<T1, T2, T3> & st) {
  for (auto sh : st.filtration_simplex_range()) {
    os << st.dimension(sh) << " ";
    for (auto v : st.simplex_vertex_range(sh)) {
      os << v << " ";
    }
    os << st.filtration(sh) << "\n";  // TODO(VR): why adding the key ?? not read ?? << "     " << st.key(sh) << " \n";
  }
  return os;
}
template<typename T1, typename T2, typename T3>
std::istream& operator>>(std::istream & is, Simplex_tree<T1, T2, T3> & st) {
  // assert(st.num_simplices() == 0);

  std::vector<typename Simplex_tree<T1, T2, T3>::Vertex_handle> simplex;
  typename Simplex_tree<T1, T2, T3>::Filtration_value fil;
  typename Simplex_tree<T1, T2, T3>::Filtration_value max_fil = 0;
  int max_dim = -1;
  size_t num_simplices = 0;
  while (read_simplex(is, simplex, fil)) {  // read all simplices in the file as a list of vertices
    ++num_simplices;
    int dim = static_cast<int>(simplex.size() - 1);  // Warning : simplex_size needs to be casted in int - Can be 0
    if (max_dim < dim) {
      max_dim = dim;
    }
    if (max_fil < fil) {
      max_fil = fil;
    }
    st.insert(simplex, fil);  // insert every simplex in the simplex tree
    simplex.clear();
  }
  st.set_num_simplices(num_simplices);
  st.set_dimension(max_dim);
  st.set_filtration(max_fil);

  return is;
}
  // end defgroup simplex_tree
}  // namespace Gudhi

#endif  // SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_H_