Efficient Construction of a k-Critical Bifiltration in SimplexTreeMulti

[odinhg]

Firstly, thank you for developing such a great tool for working with multiparameter persistence.

I am constructing a $k$-critical bifiltration using the multipers.SimplexTreeMulti class by inserting each simplex multiple times with different minimal filtration values. I follow the approach described in the documentation, using the .insert() method of the simplex tree.

My question is: Is there a more efficient way to achieve this? Any guidance on best practices for handling this kind of filtration efficiently would be greatly appreciated.

[DavidLapous]

Hey Odin, thanks for your interest in this library !

Inserting simplices one by one is indeed very slow. This needs a python call for each simplex.
I didn't think a lot about this interface, since I didn't see a lot of practical use case with $k$-critical filtrations in the wild; I'm happy to see that I was wrong :)
I thought about redirecting you to the insert_batch method, but it seems that I disabled it for the multicritical case; although it makes sense since inputs are not arrays anymore :/

Currently, the fastest option currently I think is using a "hidden" Slicer constructor (this is a class that can handle any filtration):

import multipers as mp

options = {"is_kcritical":True}
HiddenSlicer_type = mp.Slicer(return_type_only=True, **options)
sparse_boundary_matrix = ... # cf example below
simplices_dimensions = ...   # np.ndarray[ndim=1, dtype=int]
simplicies_filtrations = ... # [[f_1(splx), ..., f_k(splx) ] for splx]
HiddenSlicer_type(
    sparse_boundary_matrix,
    simplices_dimensions,
    simplicies_filtrations
)

To give you some example of these inputs:

st = mp.SimplexTreeMulti(kcritical=True)
st.insert([0,1,2], [0,0])
st.remove_maximal_simplex([0,1,2])
st.insert([0,1,2], [1,0])
st.insert([0,1,2], [0,1])
s = mp.Slicer(st, **options) # same filtration, different c++ class
B = s.get_boundaries()  # ((), (), (), (0, 1), (0, 2), (1, 2), (3, 4, 5))
F = s.get_filtrations() # [f for splx,f in st] up to some reordering
D = s.get_dimensions()  # ordered dims : [0,0,0,1,1,1,2]
assert s == SomeSlicer_type(B,D,F)

Another option is writing your filtration to an scc2020 file and read it afterwards with multipers with, e.g., mp.Slicer(file). This is the file format used by rivet, mpfree,etc.

Note that I also noticed that the function mp.slicer.to_simplextree is only working with 1-critical simplices, but the Slicer class should still be able to compute every invariants. I'll have to fix all of this. This is fixed in 234597e.

I agree that this solution is not good, so I'll have to think about a better interface.
What I have in mind currently:

1. re-enable insert_batch for k-critical simplextrees, but I don't know if this is very practical.
   with interface insert_batch(simplices, filtrations), where filtration is a num_simplex x k x num_parameters array (non-useful values can be filled with $\infty$ s). I did a first implementation of this in e8f50c1.
2. Do a helper for insert with multiple filtrations for a given simplex (but it wont be significantly faster)
3. Do you have a specific filtration in mind (e.g., Implement standard multiparameter persistence multifiltration #7 ) ?
   In that case, I think it will be immensely faster to implement it in multipers. I'd be happy to help you with that.

What do you think about this ? I'm open to anny suggestion so, if you have another idea in mind, --don't hesitate ;)

[odinhg]

With your recent fix to the insert_batch method and the clarifying example, everything seems to work perfectly. Thank you again!

Currently, I'm implementing it in Python. If you're interested, here's the repository where I'm experimenting: https://github.com/odinhg/core_multipers.

The implementation only involves constructing the Delaunay complex, computing miniball radii and distances to the $k$ nearest neighbors. I think the Python version might be sufficient for running our experiments, but I also want to develop a C++ version. And having it part of multipers would be amazing. If you're okay with it, I can email you once I have a C++ version working.

[DavidLapous]

This looks great ! Thank you Odin :)
For the implementation in multipers:

• I don't really know where to put it in the library. I think I need to do a filtration module to make it cleaner and put every filtration in this : mp.filtration.DelaunayCore, mp.filtration.RipsFunction, mp.filtration.Cubical, mp.filtration.FunctionDelaunay; what do you think about this ?
• I quickly checked that the multipers computations works with your filtration; the only issue is that mpfree doesn't seem to support multicritical filtrations yet (even though the scc2020 format supports it), so it will be hard to handle larger datasets.
• So in particular, I think that your current python implementation is OK for the moment since it won't be the bottleneck for large datasets, and we can optimize it afterward

You can already do a PR with the python version, and add the mp.filtration.DelaunayCore function if that's fine for you. I'll add the other filtrations there in the meantime.

For a faster implementation, I think the easiest route is the following.

• [in python] compute the alpha complex, knn distance if it's fast enough
• [multipers] translate it to a multipers simplextree
• [in cython/cpp] take the simplextreemulti's pointer from python, to get a multipers simplextree (with the correct simplices). I'll have to help you with that.
• [cpp] fix its filtration (for loop on the simplices).
  Useful functions for that:
  • simplextree.complex_simplex_range() (iterator of simplex handle)
  • simplextree.filtration_mutable(simplex_handle) to get the multicritical filtration value of the simplex
  • the multicritical filtration value class lives there gudhi/gudhi/Multi_critical_filtration.h (in particular, the add_generator will be useful here)

[odinhg]

Thanks for the reply, David!

I don't really know where to put it in the library. I think I need to do a filtration module to make it cleaner and put every filtration in this : mp.filtration.DelaunayCore, mp.filtration.RipsFunction, mp.filtration.Cubical, mp.filtration.FunctionDelaunay; what do you think about this ?

This sounds like a great way of organizing it! Very user-friendly.

You can already do a PR with the python version, and add the mp.filtration.DelaunayCore function if that's fine for you. I'll add the other filtrations there in the meantime.

Perfect! I will try to get the Python version working and submit a PR when I am happy with it.

For a faster implementation, I think the easiest route is the following.

[in python] compute the alpha complex, knn distance if it's fast enough
[multipers] translate it to a multipers simplextree
[in cython/cpp] take the simplextreemulti's pointer from python, to get a multipers simplextree (with the correct simplices). I'll have to help you with that.
[cpp] fix its filtration (for loop on the simplices).
Useful functions for that:
simplextree.complex_simplex_range() (iterator of simplex handle)
simplextree.filtration_mutable(simplex_handle) to get the multicritical filtration value of the simplex
the multicritical filtration value class lives there gudhi/gudhi/Multi_critical_filtration.h (in particular, the add_generator will be useful here)

This also sounds like a good long-term plan. Is there any reason not to do everything in C++ and directly provide the simplextreemulti with the entire filtration? I will try to get the Python version working first, and then I will come back to you for help with the Cython/C++ part.

[DavidLapous]

This sounds like a great way of organizing it! Very user-friendly.
Perfect! I will try to get the Python version working and submit a PR when I am happy with it.

Great ! many thanks :) I'll try to do the filtration module shortly then

Is there any reason not to do everything in C++ and directly provide the simplextreemulti with the entire filtration?

You can do everything in c++ directly, but this can lead to more work/maintenance, and I'm not really sure that it will be faster. Also, having more python code is usually easier to deal with, as it's more modular (more and smaller functions are exposed to the interface).
Tell me if I'm wrong but from what I see from your python code, the only thing that can be (significantly) speed up by doing everything in c++, is avoiding a (cpp) simplextree copy, by re-implementing the delaunay triangulation for the simplextreemulti.

These are my reasons, but as you are doing this part of the code, feel free to pick what is easier for you / what you prefer !

[odinhg]

I was mostly curious, but you make convincing points! I'm not familiar with Cython or combining C++ and Python in the same project, so I trust you to know the best. I also agree in that there's little to gain in terms of speed from doing everything in C++ for this. Thank you for clarifying, and have a good day!