Skeleton_blocker_contractor.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):
 *      - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL
 *      - YYYY/MM Author: Description of the modification
 */
 
#ifndef SKELETON_BLOCKER_CONTRACTOR_H_
#define SKELETON_BLOCKER_CONTRACTOR_H_
 
#include <gudhi/Contraction/Edge_profile.h>
#include <gudhi/Contraction/policies/Cost_policy.h>
#include <gudhi/Contraction/policies/Edge_length_cost.h>
#include <gudhi/Contraction/policies/Placement_policy.h>
#include <gudhi/Contraction/policies/First_vertex_placement.h>
#include <gudhi/Contraction/policies/Valid_contraction_policy.h>
#include <gudhi/Contraction/policies/Dummy_valid_contraction.h>  // xxx remove
#include <gudhi/Contraction/policies/Link_condition_valid_contraction.h>  // xxx remove
#include <gudhi/Contraction/policies/Contraction_visitor.h>
 
#include <gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h>
#include <gudhi/Debug_utils.h>
 
// todo remove the queue to be independent from cgal
#include <CGAL/Modifiable_priority_queue.h>
#include <CGAL/version.h>  // for CGAL_VERSION_NR
 
#include <boost/scoped_array.hpp>
#include <boost/scoped_ptr.hpp>
 
#include <memory>
#include <cassert>
#include <list>
#include <utility>  // for pair
#include <vector>
 
// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
#if CGAL_VERSION_NR < 1041101000
# error Skeleton_blocker_contractor is only available for CGAL >= 4.11
#endif
 
namespace Gudhi {
 
namespace contraction {
 
template <class Profile>
Placement_policy<Profile>* make_first_vertex_placement() {
  return new First_vertex_placement<Profile>();
}
 
template <class Profile>
Valid_contraction_policy<Profile>* make_link_valid_contraction() {
  return new Link_condition_valid_contraction<Profile>();
}
 
template <class Profile>
class Contraction_visitor_remove_popable : public Contraction_visitor<Profile> {
 public:
  typedef typename Profile::Point Point;
  typedef typename Profile::Complex Complex;
  typedef typename Complex::Vertex_handle Vertex_handle;
 
  void on_contracted(const Profile &profile, boost::optional< Point > placement) override {
    profile.complex().remove_all_popable_blockers(profile.v0_handle());
  }
};
 
template <class Profile>
Contraction_visitor<Profile>* make_remove_popable_blockers_visitor() {
  return new Contraction_visitor_remove_popable<Profile>();
}
 
template<class GeometricSimplifiableComplex, class EdgeProfile = Edge_profile<GeometricSimplifiableComplex>>
class Skeleton_blocker_contractor : private skeleton_blocker::Dummy_complex_visitor<
typename GeometricSimplifiableComplex::Vertex_handle> {
  GeometricSimplifiableComplex& complex_;
 
 public:
  typedef typename GeometricSimplifiableComplex::Graph_vertex Graph_vertex;
  typedef typename GeometricSimplifiableComplex::Vertex_handle Vertex_handle;
  typedef typename GeometricSimplifiableComplex::Simplex Simplex;
  typedef typename GeometricSimplifiableComplex::Root_vertex_handle Root_vertex_handle;
  typedef typename GeometricSimplifiableComplex::Graph_edge Graph_edge;
  typedef typename GeometricSimplifiableComplex::Edge_handle Edge_handle;
  typedef typename GeometricSimplifiableComplex::Point Point;
 
  typedef EdgeProfile Profile;
 
 
  typedef Cost_policy<Profile> Cost_policy_;
  typedef Placement_policy<Profile> Placement_policy_;
  typedef Valid_contraction_policy<Profile> Valid_contraction_policy_;
  typedef Contraction_visitor<EdgeProfile> Contraction_visitor_;
  typedef Edge_profile_factory<EdgeProfile> Edge_profile_factory_;
 
 
 
  typedef boost::optional<double> Cost_type;
  typedef boost::optional<Point> Placement_type;
 
  typedef size_t size_type;
 
  typedef Skeleton_blocker_contractor Self;
 
 private:
  struct Compare_id {
    Compare_id() : algorithm_(0) { }
 
    Compare_id(Self const* aAlgorithm) : algorithm_(aAlgorithm) { }
 
    bool operator()(Edge_handle a, Edge_handle b) const {
      return algorithm_->get_undirected_edge_id(a) < algorithm_->get_undirected_edge_id(b);
    }
 
    Self const* algorithm_;
  };
 
  struct Compare_cost {
    Compare_cost() : algorithm_(0) { }
 
    Compare_cost(Self const* aAlgorithm) : algorithm_(aAlgorithm) { }
 
    bool operator()(Edge_handle a, Edge_handle b) const {
      // NOTE: A cost is an optional<> value.
      // Absent optionals are ordered first; that is, "none < T" and "T > none" for any defined T != none.
      // In consequence, edges with undefined costs will be promoted to the top of the priority queue and popped out
      // first.
      return algorithm_->get_data(a).cost() < algorithm_->get_data(b).cost();
    }
 
    Self const* algorithm_;
  };
 
  struct Undirected_edge_id : boost::put_get_helper<size_type, Undirected_edge_id> {
    typedef boost::readable_property_map_tag category;
    typedef size_type value_type;
    typedef size_type reference;
    typedef Edge_handle key_type;
 
    Undirected_edge_id() : algorithm_(0) { }
 
    Undirected_edge_id(Self const* aAlgorithm) : algorithm_(aAlgorithm) { }
 
    size_type operator[](Edge_handle e) const {
      return algorithm_->get_undirected_edge_id(e);
    }
 
    Self const* algorithm_;
  };
 
#if CGAL_VERSION_NR < 1050500000
  typedef CGAL::Modifiable_priority_queue<Edge_handle, Compare_cost, Undirected_edge_id> PQ;
#else
  typedef CGAL::Modifiable_priority_queue<Edge_handle, Compare_cost, Undirected_edge_id, CGAL::CGAL_BOOST_PENDING_RELAXED_HEAP> PQ;
#endif
  
  typedef bool pq_handle;
 
 
  // An Edge_data is associated with EVERY edge in the complex (collapsible or not).
  // It relates the edge with the PQ-handle needed to update the priority queue
  // It also relates the edge with a policy-based cache
 
  class Edge_data {
   public:
    Edge_data() : PQHandle_(), cost_() { }
 
    Cost_type const& cost() const {
      return cost_;
    }
 
    Cost_type & cost() {
      return cost_;
    }
 
    pq_handle PQ_handle() const {
      return PQHandle_;
    }
 
    bool is_in_PQ() const {
      return PQHandle_ != false;
    }
 
    void set_PQ_handle(pq_handle h) {
      PQHandle_ = h;
    }
 
    void reset_PQ_handle() {
      PQHandle_ = false;
    }
 
   private:
    pq_handle PQHandle_;
    Cost_type cost_;
  };
  typedef Edge_data* Edge_data_ptr;
  typedef boost::scoped_array<Edge_data> Edge_data_array;
 
  int get_undirected_edge_id(Edge_handle edge) const {
    return complex_[edge].index();
  }
 
  const Edge_data& get_data(Edge_handle edge) const {
    return edge_data_array_[get_undirected_edge_id(edge)];
  }
 
  Edge_data& get_data(Edge_handle edge) {
    return edge_data_array_[get_undirected_edge_id(edge)];
  }
 
  Cost_type get_cost(const Profile & profile) const {
    return (*cost_policy_)(profile, get_placement(profile));
  }
 
  Profile create_profile(Edge_handle edge) const {
    if (edge_profile_factory_)
      return edge_profile_factory_->make_profile(complex_, edge);
    else
      return Profile(complex_, edge);
  }
 
  void insert_in_PQ(Edge_handle edge, Edge_data& data) {
    heap_PQ_->push(edge);
    data.set_PQ_handle(true);
    ++current_num_edges_heap_;
  }
 
  void update_in_PQ(Edge_handle edge, Edge_data& data) {
#if CGAL_VERSION_NR < 1050500000
    data.set_PQ_handle(heap_PQ_->update(edge, data.PQ_handle()));
#else
    heap_PQ_->update(edge);
#endif
  }
 
  void remove_from_PQ(Edge_handle edge, Edge_data& data) {
    heap_PQ_->erase(edge);
    data.set_PQ_handle(false);
    --current_num_edges_heap_;
  }
 
  boost::optional<Edge_handle> pop_from_PQ() {
    boost::optional<Edge_handle> edge = heap_PQ_->extract_top();
    if (edge)
      get_data(*edge).reset_PQ_handle();
    return edge;
  }
 
 private:
  void collect_edges() {
    //
    // Loop over all the edges in the complex in the heap
    //
    size_type size = complex_.num_edges();
    DBG("Collecting edges ...");
    DBGMSG("num edges :", size);
 
    edge_data_array_.reset(new Edge_data[size]);
 
    heap_PQ_.reset(new PQ(size, Compare_cost(this), Undirected_edge_id(this)));
 
    std::size_t id = 0;
 
    // xxx do a parallel for
    for (auto edge : complex_.edge_range()) {
      complex_[edge].index() = id++;
      Profile const& profile = create_profile(edge);
      Edge_data& data = get_data(edge);
      data.cost() = get_cost(profile);
      ++initial_num_edges_heap_;
      insert_in_PQ(edge, data);
      if (contraction_visitor_) contraction_visitor_->on_collected(profile, data.cost());
    }
 
    DBG("Edges collected.");
  }
 
  bool should_stop(double lCost, const Profile &profile) const {
    return false;
  }
 
  boost::optional<Point> get_placement(const Profile& profile) const {
    return (*placement_policy_)(profile);
  }
 
  bool is_contraction_valid(Profile const& profile, Placement_type placement) const {
    if (!valid_contraction_policy_) return true;
    return (*valid_contraction_policy_)(profile, placement);
  }
 
 
 public:
  void contract_edges(int num_max_contractions = -1) {
    DBG("\n\nContract edges");
    int num_contraction = 0;
 
    bool unspecified_num_contractions = (num_max_contractions == -1);
    //
    // Pops and processes each edge from the PQ
    //
    boost::optional<Edge_handle> edge;
    while ((edge = pop_from_PQ()) && ((num_contraction < num_max_contractions) || (unspecified_num_contractions))) {
      Profile const& profile = create_profile(*edge);
      Cost_type cost(get_data(*edge).cost());
      if (contraction_visitor_) contraction_visitor_->on_selected(profile, cost, 0, 0);
 
      DBGMSG("\n\n---- Pop edge - num vertices :", complex_.num_vertices());
 
      if (cost) {
        DBGMSG("sqrt(cost):", std::sqrt(*cost));
        if (should_stop(*cost, profile)) {
          if (contraction_visitor_) contraction_visitor_->on_stop_condition_reached();
          DBG("should_stop");
          break;
        }
        Placement_type placement = get_placement(profile);
        if (is_contraction_valid(profile, placement) && placement) {
          DBG("contraction_valid");
          contract_edge(profile, placement);
          ++num_contraction;
        } else {
          DBG("contraction not valid");
          if (contraction_visitor_) contraction_visitor_->on_non_valid(profile);
        }
      } else {
        DBG("uncomputable cost");
      }
    }
    if (contraction_visitor_) contraction_visitor_->on_stop_condition_reached();
  }
 
  bool is_in_heap(Edge_handle edge) const {
    if (heap_PQ_->empty()) {
      return false;
    } else {
      return edge_data_array_[get_undirected_edge_id(edge)].is_in_PQ();
    }
  }
 
  bool is_heap_empty() const {
    return heap_PQ_->empty();
  }
 
  boost::optional<std::pair<Edge_handle, Placement_type > > top_edge() {
    boost::optional<std::pair<Edge_handle, Placement_type > > res;
 
    if (!heap_PQ_->empty()) {
      auto edge = heap_PQ_->top();
      Profile const& profile = create_profile(edge);
      Placement_type placement = get_placement(profile);
      res = std::make_pair(edge, placement);
      DBGMSG("top edge:", complex_[edge]);
    }
    return res;
  }
 
  Skeleton_blocker_contractor(GeometricSimplifiableComplex& complex)
      : complex_(complex),
      cost_policy_(new Edge_length_cost<Profile>),
      placement_policy_(new First_vertex_placement<Profile>),
      valid_contraction_policy_(new Link_condition_valid_contraction<Profile>),
      contraction_visitor_(new Contraction_visitor_()),
      edge_profile_factory_(0),
      initial_num_edges_heap_(0),
      current_num_edges_heap_(0) {
    complex_.set_visitor(this);
    if (contraction_visitor_) contraction_visitor_->on_started(complex);
    collect_edges();
  }
 
  Skeleton_blocker_contractor(GeometricSimplifiableComplex& complex,
                              Cost_policy_ *cost_policy,
                              Placement_policy_ * placement_policy = new First_vertex_placement<Profile>,
                              Valid_contraction_policy_ * valid_contraction_policy =
                              new Link_condition_valid_contraction<Profile>,
                              Contraction_visitor_* contraction_visitor = new Contraction_visitor_(),
                              Edge_profile_factory_* edge_profile_factory = NULL) :
      complex_(complex),
      cost_policy_(cost_policy),
      placement_policy_(placement_policy),
      valid_contraction_policy_(valid_contraction_policy),
      contraction_visitor_(contraction_visitor),
      edge_profile_factory_(edge_profile_factory),
      initial_num_edges_heap_(0),
      current_num_edges_heap_(0) {
    complex_.set_visitor(this);
    if (contraction_visitor) contraction_visitor->on_started(complex);
    collect_edges();
  }
 
  ~Skeleton_blocker_contractor() {
    complex_.set_visitor(0);
  }
 
 private:
  void contract_edge(const Profile& profile, Placement_type placement) {
    if (contraction_visitor_) contraction_visitor_->on_contracting(profile, placement);
 
    assert(complex_.contains_vertex(profile.v0_handle()));
    assert(complex_.contains_vertex(profile.v1_handle()));
    assert(placement);
 
    profile.complex().point(profile.v0_handle()) = *placement;
 
    // remark : this is not necessary since v1 will be deactivated
    // profile.complex().point(profile.v1_handle()) = *placement;
 
    complex_.contract_edge(profile.v0_handle(), profile.v1_handle());
 
    assert(complex_.contains_vertex(profile.v0_handle()));
    assert(!complex_.contains_vertex(profile.v1_handle()));
 
    update_changed_edges();
 
    // the visitor could do something as complex_.remove_popable_blockers();
    if (contraction_visitor_) contraction_visitor_->on_contracted(profile, placement);
  }
 
 private:
  // every time the visitor's method on_changed_edge is called, it adds an
  // edge to changed_edges_
  std::vector< Edge_handle > changed_edges_;
 
  inline void on_changed_edge(Vertex_handle a, Vertex_handle b) override {
    boost::optional<Edge_handle> ab(complex_[std::make_pair(a, b)]);
    assert(ab);
    changed_edges_.push_back(*ab);
  }
 
  void update_changed_edges() {
    // xxx do a parallel for
    DBG("update edges");
 
    // sequential loop
    for (auto ab : changed_edges_) {
      // 1-get the Edge_handle corresponding to ab
      // 2-change the data in mEdgeArray[ab.id()]
      // 3-update the heap
      Edge_data& data = get_data(ab);
      Profile const& profile = create_profile(ab);
      data.cost() = get_cost(profile);
      if (data.is_in_PQ()) {
        update_in_PQ(ab, data);
      } else {
        insert_in_PQ(ab, data);
      }
    }
    changed_edges_.clear();
  }
 
 
 private:
  void on_remove_edge(Vertex_handle a, Vertex_handle b) override {
    boost::optional<Edge_handle> ab((complex_[std::make_pair(a, b)]));
    assert(ab);
    Edge_data& lData = get_data(*ab);
    if (lData.is_in_PQ()) {
      remove_from_PQ(*ab, lData);
    }
  }
 
 private:
  void on_swaped_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x) override {
    boost::optional<Edge_handle> ax(complex_[std::make_pair(a, x)]);
    boost::optional<Edge_handle> bx(complex_[std::make_pair(b, x)]);
    assert(ax && bx);
    complex_[*ax].index() = complex_[*bx].index();
  }
 
 private:
  void on_delete_blocker(const Simplex * blocker) override {
    // we go for all pairs xy that belongs to the blocker
    // note that such pairs xy are necessarily edges of the complex
    // by definition of a blocker
 
    // todo uniqument utile pour la link condition
    // laisser a l'utilisateur ? boolean update_heap_on_removed_blocker?
    Simplex blocker_copy(*blocker);
    for (auto x = blocker_copy.begin(); x != blocker_copy.end(); ++x) {
      for (auto y = x; ++y != blocker_copy.end();) {
        auto edge_descr(complex_[std::make_pair(*x, *y)]);
        assert(edge_descr);
        Edge_data& data = get_data(*edge_descr);
        Profile const& profile = create_profile(*edge_descr);
        data.cost() = get_cost(profile);
 
        // If the edge is already in the heap
        // its priority has not changed.
        // If the edge is not present, we reinsert it
        // remark : we could also reinsert the edge
        // only if it is valid
        if (!data.is_in_PQ()) {
          insert_in_PQ(*edge_descr, data);
        }
      }
    }
  }
 
 
 private:
  std::shared_ptr<Cost_policy_> cost_policy_;
  std::shared_ptr<Placement_policy_> placement_policy_;
  std::shared_ptr<Valid_contraction_policy_> valid_contraction_policy_;
  std::shared_ptr<Contraction_visitor_> contraction_visitor_;
 
  // in case the user wants to do something special when the edge profile
  // are created (for instance add some info)
  std::shared_ptr<Edge_profile_factory_> edge_profile_factory_;
  Edge_data_array edge_data_array_;
 
  boost::scoped_ptr<PQ> heap_PQ_;
  int initial_num_edges_heap_;
  int current_num_edges_heap_;
};
 
}  // namespace contraction
 
}  // namespace Gudhi
 
#endif  // SKELETON_BLOCKER_CONTRACTOR_H_