GUDHI / gudhi-devel

The GUDHI library is a generic open source C++ library, with a Python interface, for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.
https://gudhi.inria.fr/
MIT License
258 stars 66 forks source link

Betti numbers computation on 0-isolated simplices crashes #1118

Closed VincentRouvreau closed 1 month ago

VincentRouvreau commented 2 months ago
  Simplex_tree st;

  st.insert_simplex_and_subfaces({0});
  st.insert_simplex_and_subfaces({1});
  st.insert_simplex_and_subfaces({2});
  st.insert_simplex_and_subfaces({3});
  st.insert_simplex_and_subfaces({4});

  // Sort the simplices in the order of the filtration
  st.initialize_filtration();

  // Class for homology computation
  St_persistence pcoh(st);

  // Initialize the coefficient field Z/11Z for homology
  pcoh.init_coefficients(3);

  // Compute the persistence diagram of the complex
  pcoh.compute_persistent_cohomology();

  std::vector<int> bns = pcoh.betti_numbers();
  // -> unknown location(0): fatal error: in "betti_numbers_isolated_zero_simplices": memory access violation at address: 0x0: no mapping at fault address
mglisse commented 2 months ago

The code of betti_numbers looks like

    std::vector<int> betti_numbers(std::max(dim_max_, 0));

    for (auto pair : persistent_pairs_) {
      // Count never ended persistence intervals
      if (cpx_->null_simplex() == get<1>(pair)) {
        // Increment corresponding betti number
        betti_numbers[cpx_->dimension(get<0>(pair))] += 1;
      }
    }

In this example, dim_max_ is 0, so we create a vector of size 0, but persistent_pairs_ is not empty, it contains infinite intervals for the vertices. Note that persistence_dim_max is false. If I set it to true, we create a vector of the right size, and things work. Arguably, the bug is that we computed H0 even though the combination of dimension() and persistence_dim_max tells us not to. A completely different testcase:

import gudhi
st=gudhi.SimplexTree()
st.insert([0])
st.insert([1])
st.persistence()

[(0, (0.0, inf)), (0, (0.0, inf))]

but I think according to the documentation it should return [].