Skip to content

Commit

Permalink
fix: GudhiCohomology inf trick is not ready yet
Browse files Browse the repository at this point in the history
  • Loading branch information
@DavidLapous
DavidLapous committed Oct 20, 2024
1 parent 0459cff commit ca2417c
Show file tree
Hide file tree
Showing 3 changed files with 48 additions and 61 deletions.
5 changes: 3 additions & 2 deletions multipers/gudhi/mma_interface_coh.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,8 +136,9 @@ class Persistence_backend_cohomology {
dimension_type get_dimension(pos_index i) { return matrix_.dimension(matrix_.simplex(i)); }

Barcode get_barcode() {
Persistent_cohomology pcoh(matrix_);
pcoh.init_coefficients(2);
Persistent_cohomology pcoh(matrix_, true);
constexpr int coeff_field_characteristic = 11;
pcoh.init_coefficients(coeff_field_characteristic);
pcoh.compute_persistent_cohomology(0);

const auto &pairs = pcoh.get_persistent_pairs();
Expand Down
72 changes: 31 additions & 41 deletions multipers/multi_parameter_rank_invariant/persistence_slices.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,35 +7,34 @@

namespace Gudhi::multiparameter {

struct Simplex_tree_float { // smaller simplextrees
struct Simplex_tree_float { // smaller simplextrees
typedef linear_indexing_tag Indexing_tag;
typedef std::int32_t Vertex_handle;
typedef float Filtration_value;
typedef std::uint32_t Simplex_key;
static const bool store_key = true;
static const bool store_filtration = true;
static const bool contiguous_vertices =
false; // TODO OPTIMIZATION : maybe make the simplextree contiguous when
// calling grid_squeeze ?
static const bool contiguous_vertices = false; // TODO OPTIMIZATION : maybe make the simplextree contiguous when
// calling grid_squeeze ?
static const bool link_nodes_by_label = true;
static const bool stable_simplex_handles = false;
static const bool is_multi_parameter = false;
};

// using Simplex_tree_float = Simplex_tree_options_fast_persistence;

using Gudhi::multi_persistence::Box;
using Simplex_tree_std = Simplex_tree<Simplex_tree_float>;

using Barcode = std::vector<std::pair<Simplex_tree_std::Filtration_value,
Simplex_tree_std::Filtration_value>>;
using Barcode = std::vector<std::pair<Simplex_tree_std::Filtration_value, Simplex_tree_std::Filtration_value>>;

inline Barcode compute_dgm(Simplex_tree_std &st, int degree) {
st.initialize_filtration(true);
constexpr int coeff_field_characteristic = 11;
constexpr Simplex_tree_std::Filtration_value min_persistence = 0;
bool persistence_dim_max = st.dimension() == degree;
Gudhi::persistent_cohomology::Persistent_cohomology<
Simplex_tree_std, Gudhi::persistent_cohomology::Field_Zp>
pcoh(st, persistence_dim_max);
Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree_std, Gudhi::persistent_cohomology::Field_Zp> pcoh(
st, persistence_dim_max);
pcoh.init_coefficients(coeff_field_characteristic);
pcoh.compute_persistent_cohomology(min_persistence);
const auto &persistent_pairs = pcoh.intervals_in_dimension(degree);
Expand All @@ -46,9 +45,10 @@ inline Barcode compute_dgm(Simplex_tree_std &st, int degree) {
}

template <typename degree_type, class interface_std_like>
inline std::vector<Barcode>
compute_dgms(interface_std_like &st, const std::vector<degree_type> &degrees,
int num_collapses, int expansion_dim) {
inline std::vector<Barcode> compute_dgms(interface_std_like &st,
const std::vector<degree_type> &degrees,
int num_collapses,
int expansion_dim) {
std::vector<Barcode> out(degrees.size());
static_assert(!interface_std_like::Options::is_multi_parameter,
"Can only compute persistence for 1-parameter simplextrees.");
Expand All @@ -65,7 +65,7 @@ compute_dgms(interface_std_like &st, const std::vector<degree_type> &degrees,
st.expansion(expansion_dim);
}

st.initialize_filtration(true); // true is ignore_infinite_values
st.initialize_filtration(true); // true is ignore_infinite_values
constexpr int coeff_field_characteristic = 11;
constexpr typename interface_std_like::Filtration_value min_persistence = 0;

Expand All @@ -76,12 +76,12 @@ compute_dgms(interface_std_like &st, const std::vector<degree_type> &degrees,
break;
}
}
//
if constexpr (verbose){
//
if constexpr (verbose) {
std::cout << "Computing dgm of st:\n";
for (auto& sh : st.filtration_simplex_range()){
std::cout <<"dim: "<< st.dimension(sh) << " vertices: ";
for (auto v : st.simplex_vertex_range(sh)){
for (auto &sh : st.filtration_simplex_range()) {
std::cout << "dim: " << st.dimension(sh) << " vertices: ";
for (auto v : st.simplex_vertex_range(sh)) {
std::cout << v << " ";
}
std::cout << " filtration: ";
Expand All @@ -90,45 +90,37 @@ compute_dgms(interface_std_like &st, const std::vector<degree_type> &degrees,
return out;
}



Gudhi::persistent_cohomology::Persistent_cohomology<
interface_std_like, Gudhi::persistent_cohomology::Field_Zp>
pcoh(st, persistence_dim_max);
Gudhi::persistent_cohomology::Persistent_cohomology<interface_std_like, Gudhi::persistent_cohomology::Field_Zp> pcoh(
st, persistence_dim_max);
pcoh.init_coefficients(coeff_field_characteristic);
pcoh.compute_persistent_cohomology(min_persistence);
for (auto i = 0u; i < degrees.size(); i++) {
out[i] = pcoh.intervals_in_dimension(degrees[i]);
}
return out;
}

// small wrapper
template <typename degree_type, class interface_std_like>
inline std::vector<Barcode>
compute_dgms(interface_std_like &st, const std::vector<degree_type> &degrees,
int expand_collapse_dim = 0) {
if (expand_collapse_dim > 0)
return compute_dgms(st, degrees, 10, expand_collapse_dim);
inline std::vector<Barcode> compute_dgms(interface_std_like &st,
const std::vector<degree_type> &degrees,
int expand_collapse_dim = 0) {
if (expand_collapse_dim > 0) return compute_dgms(st, degrees, 10, expand_collapse_dim);
return compute_dgms(st, degrees, 0, 0);
}

// Adapted version from the Simplextree interface
template <class simplextree_like>
inline simplextree_like collapse_edges(simplextree_like &st,
int num_collapses) {
inline simplextree_like collapse_edges(simplextree_like &st, int num_collapses) {
using Filtered_edge = std::tuple<typename simplextree_like::Vertex_handle,
typename simplextree_like::Vertex_handle,
typename simplextree_like::Filtration_value>;
std::vector<Filtered_edge> edges;
for (auto sh : st.skeleton_simplex_range(1)) {
if (st.dimension(sh) == 1) {
const auto filtration = st.filtration(sh);
if (filtration ==
std::numeric_limits<
typename simplextree_like::Filtration_value>::infinity())
continue;
typename simplextree_like::Simplex_vertex_range rg =
st.simplex_vertex_range(sh);
if (filtration == std::numeric_limits<typename simplextree_like::Filtration_value>::infinity()) continue;
typename simplextree_like::Simplex_vertex_range rg = st.simplex_vertex_range(sh);
auto vit = rg.begin();
typename simplextree_like::Vertex_handle v = *vit;
typename simplextree_like::Vertex_handle w = *++vit;
Expand All @@ -138,14 +130,12 @@ inline simplextree_like collapse_edges(simplextree_like &st,
for (int iteration = 0; iteration < num_collapses; iteration++) {
auto current_size = edges.size();
edges = Gudhi::collapse::flag_complex_collapse_edges(std::move(edges));
if (edges.size() >= current_size)
break; // no need to do more
if (edges.size() >= current_size) break; // no need to do more
}
simplextree_like collapsed_stree;
// Copy the original 0-skeleton
for (auto sh : st.skeleton_simplex_range(0)) {
collapsed_stree.insert_simplex({*(st.simplex_vertex_range(sh).begin())},
st.filtration(sh));
collapsed_stree.insert_simplex({*(st.simplex_vertex_range(sh).begin())}, st.filtration(sh));
}
// Insert remaining edges
for (auto [x, y, filtration] : edges) {
Expand All @@ -154,4 +144,4 @@ inline simplextree_like collapse_edges(simplextree_like &st,
return collapsed_stree;
}

} // namespace Gudhi::multiparameter
} // namespace Gudhi::multiparameter
32 changes: 14 additions & 18 deletions multipers/multi_parameter_rank_invariant/rank_invariant.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -165,14 +165,14 @@ template <class PersBackend,
typename dtype,
typename index_type>
inline void compute_2d_rank_invariant_of_elbow(
interface::Truc<PersBackend, Structure, MultiFiltration> &slicer, // truc slicer
const tensor::static_tensor_view<dtype, index_type> &out, // assumes its a zero tensor
typename interface::Truc<PersBackend, Structure, MultiFiltration>::ThreadSafe &slicer, // truc slicer
const tensor::static_tensor_view<dtype, index_type> &out, // assumes its a zero tensor
const index_type I,
const index_type J,
const std::vector<index_type> &grid_shape,
const std::vector<index_type> &degrees,
std::vector<Index> &order_container, // constant size
std::vector<typename MultiFiltration::value_type> &one_persistence, // constant size
// std::vector<Index> &order_container, // constant size
// std::vector<typename MultiFiltration::value_type> &one_persistence, // constant size
const bool flip_death = false) {
using value_type = typename MultiFiltration::value_type;
const auto &filtrations_values = slicer.get_filtrations();
Expand All @@ -198,7 +198,7 @@ inline void compute_2d_rank_invariant_of_elbow(
filtration_in_slice = get_slice_rank_filtration(x, y, I, J);
}
if constexpr (verbose) std::cout << filtration_in_slice << ",";
one_persistence[i] = filtration_in_slice;
slicer.get_one_filtration()[i] = filtration_in_slice;
}
if constexpr (verbose) std::cout << "\b]" << std::endl;

Expand All @@ -211,7 +211,7 @@ inline void compute_2d_rank_invariant_of_elbow(
using bc_type = typename interface::Truc<PersBackend, Structure, MultiFiltration>::split_barcode;
bc_type barcodes;
if constexpr (PersBackend::is_vine) {
slicer.set_one_filtration(one_persistence);
// slicer.set_one_filtration(one_persistence);
if (I == 0 && J == 0) [[unlikely]] // this is dangerous, assumes it starts at 0 0
{
// TODO : This is a good optimization but needs a patch on
Expand All @@ -228,8 +228,8 @@ inline void compute_2d_rank_invariant_of_elbow(
}
barcodes = slicer.get_barcode();
} else {
PersBackend pers = slicer.compute_persistence_out(one_persistence, order_container);
barcodes = slicer.get_barcode(pers, one_persistence);
slicer.template compute_persistence<false>(); // PUT THAT TO TRUE...
barcodes = slicer.get_barcode();
}

// note one_pers not necesary when vine, but does the same computation
Expand Down Expand Up @@ -286,19 +286,15 @@ inline void compute_2d_rank_invariant(
std::cout << "Shape " << grid_shape[0] << " " << grid_shape[1] << " " << grid_shape[2] << " " << grid_shape[3]
<< " " << grid_shape[4] << std::endl;

using local_type = std::tuple<std::vector<Index>, // order
std::vector<typename MultiFiltration::value_type> // one
// filtration
>;
local_type local_template = {std::vector<Index>(slicer.num_generators()),
std::vector<typename MultiFiltration::value_type>(slicer.num_generators())};
tbb::enumerable_thread_specific<local_type> thread_locals(local_template);
using ThreadSafe = typename interface::Truc<PersBackend, Structure, MultiFiltration>::ThreadSafe;
ThreadSafe slicer_thread(slicer);
tbb::enumerable_thread_specific<ThreadSafe> thread_locals(slicer_thread);
tbb::parallel_for(0, X, [&](index_type I) {
tbb::parallel_for(0, Y, [&](index_type J) {
if constexpr (verbose) std::cout << "Computing elbow " << I << " " << J << "...";
auto &[order_container, one_filtration_container] = thread_locals.local();
compute_2d_rank_invariant_of_elbow(
slicer, out, I, J, grid_shape, degrees, order_container, one_filtration_container, flip_death);
ThreadSafe &slicer = thread_locals.local();
compute_2d_rank_invariant_of_elbow<PersBackend, Structure, MultiFiltration, dtype, index_type>(
slicer, out, I, J, grid_shape, degrees, flip_death);
if constexpr (verbose) std::cout << "Done!" << std::endl;
});
});
Expand Down

0 comments on commit ca2417c

Please sign in to comment.