Skeleton_blocker_link_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):       David Salinas
 *
 *    Copyright (C) 2014 Inria
 *
 *    Modification(s):
 *      - YYYY/MM Author: Description of the modification
 */
 
#ifndef SKELETON_BLOCKER_LINK_COMPLEX_H_
#define SKELETON_BLOCKER_LINK_COMPLEX_H_
 
#include <gudhi/Skeleton_blocker_complex.h>
#include <gudhi/Debug_utils.h>
 
namespace Gudhi {
 
namespace skeleton_blocker {
 
template<class ComplexType> class Skeleton_blocker_sub_complex;
 
template<typename ComplexType>
class Skeleton_blocker_link_complex : public Skeleton_blocker_sub_complex<
ComplexType> {
  template<typename T> friend class Skeleton_blocker_link_superior;
  typedef typename ComplexType::Edge_handle Edge_handle;
 
  typedef typename ComplexType::boost_vertex_handle boost_vertex_handle;
 
 private:
  bool only_superior_vertices_;
 
 public:
  typedef typename ComplexType::Vertex_handle Vertex_handle;
  typedef typename ComplexType::Root_vertex_handle Root_vertex_handle;
 
  typedef typename ComplexType::Simplex Simplex;
  typedef typename ComplexType::Root_simplex_handle Root_simplex_handle;
 
  typedef typename ComplexType::Blocker_handle Blocker_handle;
 
  typedef typename ComplexType::Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_handle_iterator;
 
  explicit Skeleton_blocker_link_complex(bool only_superior_vertices = false)
      : only_superior_vertices_(only_superior_vertices) { }
 
  Skeleton_blocker_link_complex(const ComplexType & parent_complex,
                                const Simplex& alpha_parent_address,
                                bool only_superior_vertices = false,
                                bool only_vertices = false)
      : only_superior_vertices_(only_superior_vertices) {
    if (!alpha_parent_address.empty())
      build_link(parent_complex, alpha_parent_address, only_vertices);
  }
 
  Skeleton_blocker_link_complex(const ComplexType & parent_complex,
                                Vertex_handle a_parent_address,
                                bool only_superior_vertices = false)
      : only_superior_vertices_(only_superior_vertices) {
    Simplex alpha_simplex(a_parent_address);
    build_link(parent_complex, alpha_simplex);
  }
 
  Skeleton_blocker_link_complex(const ComplexType & parent_complex,
                                Edge_handle edge, bool only_superior_vertices =
                                false)
      : only_superior_vertices_(only_superior_vertices) {
    Simplex alpha_simplex(parent_complex.first_vertex(edge),
                          parent_complex.second_vertex(edge));
    build_link(parent_complex, alpha_simplex);
  }
 
 protected:
  void compute_link_vertices(const ComplexType & parent_complex,
                             const Simplex& alpha_parent_address,
                             bool only_superior_vertices,
                             bool is_alpha_blocker = false) {
    if (alpha_parent_address.dimension() == 0) {
      // for a vertex we know exactly the number of vertices of the link (and the size of the corresponding vector)
      // thus we call a specific function that will reserve a vector with appropriate size
      this->compute_link_vertices(parent_complex,
                                  alpha_parent_address.first_vertex(),
                                  only_superior_vertices_);
    } else {
      // we compute the intersection of neighbors of alpha and store it in link_vertices
      Simplex link_vertices_parent;
      parent_complex.add_neighbours(alpha_parent_address, link_vertices_parent,
                                    only_superior_vertices);
      // For all vertex 'v' in this intersection, we go through all its adjacent blockers.
      // If one blocker minus 'v' is included in alpha then the vertex is not in the link complex.
      for (auto v_parent : link_vertices_parent) {
        bool new_vertex = true;
        for (auto beta : parent_complex.const_blocker_range(v_parent)) {
          if (!is_alpha_blocker || *beta != alpha_parent_address) {
            new_vertex = !(alpha_parent_address.contains_difference(*beta,
                                                                   v_parent));
            if (!new_vertex)
              break;
          }
        }
        if (new_vertex)
          this->add_vertex(parent_complex.get_id(v_parent));
      }
    }
  }
 
  void compute_link_vertices(const ComplexType & parent_complex,
                             Vertex_handle alpha_parent_address,
                             bool only_superior_vertices) {
    // for a vertex we know exactly the number of vertices of the link (and the size of the corresponding vector
    this->skeleton.m_vertices.reserve(
                                      parent_complex.degree(alpha_parent_address));
 
    // For all vertex 'v' in this intersection, we go through all its adjacent blockers.
    // If one blocker minus 'v' is included in alpha then the vertex is not in the link complex.
    for (auto v_parent : parent_complex.vertex_range(alpha_parent_address)) {
      if (!only_superior_vertices
          || v_parent.vertex > alpha_parent_address.vertex)
        this->add_vertex(parent_complex.get_id(v_parent));
    }
  }
 
  void compute_link_edges(const ComplexType & parent_complex,
                          const Simplex& alpha_parent_address,
                          bool is_alpha_blocker = false) {
    if (this->num_vertices() <= 1)
      return;
 
    for (auto x_link = this->vertex_range().begin();
         x_link != this->vertex_range().end(); ++x_link) {
      for (auto y_link = x_link; ++y_link != this->vertex_range().end();) {
        Vertex_handle x_parent = *parent_complex.get_address(
                                                             this->get_id(*x_link));
        Vertex_handle y_parent = *parent_complex.get_address(
                                                             this->get_id(*y_link));
        if (parent_complex.contains_edge(x_parent, y_parent)) {
          // we check that there is no blocker subset of alpha passing through x and y
          bool new_edge = true;
          for (auto blocker_parent : parent_complex.const_blocker_range(
                                                                        x_parent)) {
            if (!is_alpha_blocker || *blocker_parent != alpha_parent_address) {
              if (blocker_parent->contains(y_parent)) {
                new_edge = !(alpha_parent_address.contains_difference(
                                                                     *blocker_parent, x_parent, y_parent));
                if (!new_edge)
                  break;
              }
            }
          }
          if (new_edge)
            this->add_edge_without_blockers(*x_link, *y_link);
        }
      }
    }
  }
 
  boost::optional<Vertex_handle> give_equivalent_vertex(const ComplexType & other_complex,
                                                        Vertex_handle address) const {
    Root_vertex_handle id((*this)[address].get_id());
    return other_complex.get_address(id);
  }
 
  /*
   * compute the blockers of the link if is_alpha_blocker is false.
   * Otherwise, alpha is a blocker, and the link is computed in the complex where
   * the blocker is anticollapsed.
   */
  void compute_link_blockers(const ComplexType & parent_complex,
                             const Simplex& alpha_parent,
                             bool is_alpha_blocker = false) {
    for (auto x_link : this->vertex_range()) {
      Vertex_handle x_parent = *this->give_equivalent_vertex(parent_complex,
                                                             x_link);
 
      for (auto blocker_parent : parent_complex.const_blocker_range(x_parent)) {
        if (!is_alpha_blocker || *blocker_parent != alpha_parent) {
          Simplex sigma_parent(*blocker_parent);
 
          sigma_parent.difference(alpha_parent);
 
          if (sigma_parent.dimension() >= 2
              && sigma_parent.first_vertex() == x_parent) {
            Root_simplex_handle sigma_id(parent_complex.get_id(sigma_parent));
            auto sigma_link = this->get_simplex_address(sigma_id);
            // ie if the vertices of sigma are vertices of the link
            if (sigma_link) {
              bool is_new_blocker = true;
              for (auto a : alpha_parent) {
                for (auto eta_parent : parent_complex.const_blocker_range(a)) {
                  if (!is_alpha_blocker || *eta_parent != alpha_parent) {
                    Simplex eta_minus_alpha(*eta_parent);
                    eta_minus_alpha.difference(alpha_parent);
                    if (eta_minus_alpha != sigma_parent
                        && sigma_parent.contains_difference(*eta_parent,
                                                            alpha_parent)) {
                      is_new_blocker = false;
                      break;
                    }
                  }
                }
                if (!is_new_blocker)
                  break;
              }
              if (is_new_blocker)
                this->add_blocker(new Simplex(*sigma_link));
            }
          }
        }
      }
    }
  }
 
 public:
  void build_link(const ComplexType & parent_complex,
                  const Simplex& alpha_parent_address,
                  bool is_alpha_blocker = false,
                  bool only_vertices = false) {
    assert(is_alpha_blocker || parent_complex.contains(alpha_parent_address));
    compute_link_vertices(parent_complex, alpha_parent_address, only_superior_vertices_);
    if (!only_vertices) {
      compute_link_edges(parent_complex, alpha_parent_address, is_alpha_blocker);
      compute_link_blockers(parent_complex, alpha_parent_address, is_alpha_blocker);
    }
  }
 
  friend void build_link_of_blocker(const ComplexType & parent_complex,
                                    Simplex& blocker,
                                    Skeleton_blocker_link_complex & result) {
    assert(blocker.dimension() >= 2);
    assert(parent_complex.contains_blocker(blocker));
    result.clear();
    result.build_link(parent_complex, blocker, true);
  }
};
 
}  // namespace skeleton_blocker
 
namespace skbl = skeleton_blocker;
 
}  // namespace Gudhi
 
#endif  // SKELETON_BLOCKER_LINK_COMPLEX_H_