edge_collapse_conserve_persistence.cpp
/* 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): Vincent Rouvreau
*
* Copyright (C) 2020 Inria
*
* Modification(s):
* - YYYY/MM Author: Description of the modification
*/
#include <gudhi/Flag_complex_edge_collapser.h>
#include <gudhi/Simplex_tree.h>
#include <gudhi/Persistent_cohomology.h>
#include <gudhi/Points_off_io.h>
#include <boost/range/adaptor/transformed.hpp>
#include<utility> // for std::pair
#include<vector>
#include<tuple>
// Types definition
using Point = std::vector<Filtration_value>;
using Vector_of_points = std::vector<Point>;
using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>;
using Persistence_interval = std::tuple<int, Filtration_value, Filtration_value>;
/*
* Compare two intervals by dimension, then by length.
*/
struct cmp_intervals_by_length {
explicit cmp_intervals_by_length(Simplex_tree * sc)
: sc_(sc) { }
template<typename Persistent_interval>
bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) {
return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1))
> sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2)));
}
};
std::vector<Persistence_interval> get_persistence_intervals(Simplex_tree& st, int ambient_dim) {
std::vector<Persistence_interval> persistence_intervals;
st.expansion(ambient_dim);
// Sort the simplices in the order of the filtration
// Compute the persistence diagram of the complex
// initializes the coefficient field for homology - must be a prime number
int p = 11;
pcoh.init_coefficients(p);
// Default min_interval_length = 0.
pcoh.compute_persistent_cohomology();
// Custom sort and output persistence
cmp_intervals_by_length cmp(&st);
auto persistent_pairs = pcoh.get_persistent_pairs();
std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp);
for (auto pair : persistent_pairs) {
persistence_intervals.emplace_back(st.dimension(get<0>(pair)),
st.filtration(get<0>(pair)),
st.filtration(get<1>(pair)));
}
return persistence_intervals;
}
int main(int argc, char* argv[]) {
if (argc != 3) {
std::cerr << "This program requires an OFF file and minimal threshold value as parameter\n";
std::cerr << "For instance: ./Edge_collapse_conserve_persistence ../../data/points/tore3D_300.off 1.\n";
exit(-1); // ----- >>
}
std::string off_file_points {argv[1]};
double threshold {atof(argv[2])};
Gudhi::Points_off_reader<Point> off_reader(off_file_points);
if (!off_reader.is_valid()) {
std::cerr << "Unable to read file " << off_file_points << "\n";
exit(-1); // ----- >>
}
Vector_of_points point_vector = off_reader.get_point_cloud();
if (point_vector.size() <= 0) {
std::cerr << "Empty point cloud." << std::endl;
exit(-1); // ----- >>
}
Proximity_graph proximity_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(),
threshold,
if (num_edges(proximity_graph) <= 0) {
std::cerr << "Total number of edges is zero." << std::endl;
exit(-1);
}
int ambient_dim = point_vector[0].size();
// ***** Simplex tree from a flag complex built after collapse *****
boost::adaptors::transform(edges(proximity_graph), [&](auto&&edge){
return std::make_tuple(static_cast<Vertex_handle>(source(edge, proximity_graph)),
static_cast<Vertex_handle>(target(edge, proximity_graph)),
get(Gudhi::edge_filtration_t(), proximity_graph, edge));
})
);
Simplex_tree stree_from_collapse;
for (Vertex_handle vertex = 0; static_cast<std::size_t>(vertex) < point_vector.size(); vertex++) {
// insert the vertex with a 0. filtration value just like a Rips
stree_from_collapse.insert_simplex({vertex}, 0.);
}
for (auto remaining_edge : remaining_edges) {
stree_from_collapse.insert_simplex({std::get<0>(remaining_edge), std::get<1>(remaining_edge)},
std::get<2>(remaining_edge));
}
std::vector<Persistence_interval> persistence_intervals_from_collapse = get_persistence_intervals(stree_from_collapse, ambient_dim);
// ***** Simplex tree from the complete flag complex *****
Simplex_tree stree_wo_collapse;
stree_wo_collapse.insert_graph(proximity_graph);
std::vector<Persistence_interval> persistence_intervals_wo_collapse = get_persistence_intervals(stree_wo_collapse, ambient_dim);
// ***** Comparison *****
if (persistence_intervals_wo_collapse.size() != persistence_intervals_from_collapse.size()) {
std::cerr << "Number of persistence pairs with collapse is " << persistence_intervals_from_collapse.size() << std::endl;
std::cerr << "Number of persistence pairs without collapse is " << persistence_intervals_wo_collapse.size() << std::endl;
exit(-1);
}
int return_value = 0;
auto ppwoc_ptr = persistence_intervals_wo_collapse.begin();
for (auto ppfc: persistence_intervals_from_collapse) {
if (ppfc != *ppwoc_ptr) {
return_value++;
std::cerr << "Without collapse: "
<< std::get<0>(*ppwoc_ptr) << " " << std::get<1>(*ppwoc_ptr) << " " << std::get<2>(*ppwoc_ptr)
<< " - With collapse: "
<< std::get<0>(ppfc) << " " << std::get<1>(ppfc) << " " << std::get<2>(ppfc) << std::endl;
}
ppwoc_ptr++;
}
return return_value;
}
Compute the Euclidean distance between two Points given by a range of coordinates....
Definition: distance_functions.h:32
OFF file reader implementation in order to read points from an OFF file.
Definition: Points_off_io.h:122
Options::Filtration_value Filtration_value
Type for the value of the filtration function.
Definition: Simplex_tree.h:86
std::pair< Simplex_handle, bool > insert_simplex(const InputVertexRange &simplex, Filtration_value filtration=0)
Insert a simplex, represented by a range of Vertex_handles, in the simplicial complex.
Definition: Simplex_tree.h:778
static Filtration_value filtration(Simplex_handle sh)
Returns the filtration value of a simplex.
Definition: Simplex_tree.h:535
void initialize_filtration()
Initializes the filtration cache, i.e. sorts the simplices according to their order in the filtration...
Definition: Simplex_tree.h:912
Options::Vertex_handle Vertex_handle
Type for the vertex handle.
Definition: Simplex_tree.h:94
void expansion(int max_dim)
Expands the Simplex_tree containing only its one skeleton until dimension max_dim.
Definition: Simplex_tree.h:1170
int dimension(Simplex_handle sh)
Returns the dimension of a simplex.
Definition: Simplex_tree.h:600
void insert_graph(const OneSkeletonGraph &skel_graph)
Inserts a 1-skeleton in an empty Simplex_tree.
Definition: Simplex_tree.h:1106
Structure representing the coefficient field .
Definition: Field_Zp.h:27
Computes the persistent cohomology of a filtered complex.
Definition: Persistent_cohomology.h:52
Global distance functions.
Graph simplicial complex methods.
typename boost::adjacency_list< boost::vecS, boost::vecS, boost::directedS, boost::property< vertex_filtration_t, typename SimplicialComplexForProximityGraph::Filtration_value >, boost::property< edge_filtration_t, typename SimplicialComplexForProximityGraph::Filtration_value > > Proximity_graph
Proximity_graph contains the vertices and edges with their filtration values in order to store the re...
Definition: graph_simplicial_complex.h:45
auto flag_complex_collapse_edges(const FilteredEdgeRange &edges)
Implicitly constructs a flag complex from edges as an input, collapses edges while preserving the per...
Definition: Flag_complex_edge_collapser.h:329
Value type for a filtration function on a cell complex.
Definition: FiltrationValue.h:20
Handle type for the vertices of a cell complex.
Definition: VertexHandle.h:15