Implement standard multiparameter persistence multifiltration

[DavidLapous]

• Function-simplextree bifiltrations, e.g., Degree-Rips. Done, but I'm not satisfied with the interface
• Multicover bifiltration, using Rhomboid tiling
• Delaunay Bifiltrations
• Sparse Multicover bifiltration
• Tell me

[odinhg]

I think the rhomboid tiling bifiltration should be doable if it's okay to call an external binary. Of course this would result in some file I/O overhead, which is not optimal. And the implementation [1] currently only supports point clouds in $\mathbb{R}^d$ for $d=2,3$. What remains for me to figure out is how to build it with a newer version of CGAL so it can use the same toolchain as mpfree currently uses when installing it with multipers. Currently, I have been using old Nix packages to get the versions correct.

If this is an okay (preliminary) solution for multipers, I can probably look into it.

[1] https://github.com/geoo89/rhomboidtiling

[DavidLapous]

Thank you for your help Odin. I actually looked at doing this a couple of years ago, but I didn't integrate it to multipers because of this old CGAL version issue. IIRC, I had it working in a conda environment up to some patches, but I didn't go much further than this.

If you want to take a look into it, it'd be very helpful ! I'm fine with the I/O overhead, as it will definitely not be the bottleneck here.
Furthermore, the input file can be created efficiently with numpy, and rhomboidtiling supports firep files, which are scc files up to some sed-ing, and can thus efficiently be loaded into multipers, with the Slicer._build_from_scc_file method.
For the integration in multipers, doing this would be very similar to how function_delaunay is handled currently:

• an interface function in

  multipers/multipers/filtrations/filtrations.py

  Lines 111 to 174 in ce14eed

  	def DelaunayLowerstar(
  	points: ArrayLike,
  	function: ArrayLike,
  	*,
  	distance_matrix: Optional[ArrayLike] = None,
  	threshold_radius: Optional[float] = None,
  	reduce_degree: int = -1,
  	vineyard: Optional[bool] = None,
  	dtype=np.float64,
  	verbose: bool = False,
  	clear: bool = True,
  	flagify: bool = False,
  	):
  	"""
  	Computes the Function Delaunay bifiltration. Similar to RipsLowerstar, but most suited for low-dimensional euclidean data.
  	See [Delaunay bifiltrations of functions on point clouds, Alonso et al] https://doi.org/10.1137/1.9781611977912.173
  	Input:
  	- points or distance_matrix: ArrayLike
  	- function : ArrayLike of shape (num_data, )
  	- threshold_radius: max edge length of the rips. Defaults at min(max(distance_matrix, axis=1)).
  	"""
  	from multipers.slicer import from_function_delaunay
  	if flagify and reduce_degree >= 0:
  	raise ValueError(
  	"Got {reduce_degree=} and {flagify=}. Cannot flagify with reduce degree."
  	)
  	assert distance_matrix is None, "Delaunay cannot be built from distance matrices"
  	if threshold_radius is not None:
  	raise NotImplementedError("Delaunay with threshold not implemented yet.")
  	api = api_from_tensors(points, function)
  	if not flagify and (api.has_grad(points) or api.has_grad(function)):
  	warn("Cannot keep points gradient unless using `flagify=True`.")
  	points = api.astensor(points)
  	function = api.astensor(function).squeeze()
  	assert (
  	function.ndim == 1
  	), "Delaunay Lowerstar is only compatible with 1 additional parameter."
  	slicer = from_function_delaunay(
  	api.asnumpy(points),
  	api.asnumpy(function),
  	degree=reduce_degree,
  	vineyard=vineyard,
  	dtype=dtype,
  	verbose=verbose,
  	clear=clear,
  	)
  	if reduce_degree >= 0:
  	# Force resolution to avoid confusion with hilbert.
  	slicer = slicer.minpres(degree=reduce_degree, force=True)
  	if flagify:
  	from multipers.slicer import to_simplextree
  	slicer = to_simplextree(slicer)
  	slicer.flagify(2)
  	if api.has_grad(points) or api.has_grad(function):
  	distances = api.cdist(points, points) / 2
  	grid = compute_grid([distances.ravel(), function])
  	slicer = slicer.grid_squeeze(grid)
  	slicer = slicer._clean_filtration_grid()
  	return slicer

• an io function in

  multipers/multipers/io.pyx

  Lines 473 to 515 in ce14eed

  	def function_delaunay_presentation_to_slicer(
  	slicer,
  	point_cloud:np.ndarray,
  	function_values:np.ndarray,
  	id:Optional[str] = None,
  	bool clear:bool = True,
  	bool verbose:bool=False,
  	int degree = -1,
  	bool multi_chunk = False,
  	):
  	"""
  	Computes a function delaunay presentation, and returns it as a slicer.
  	slicer: empty slicer to fill
  	points : (num_pts, n) float array
  	grades : (num_pts,) float array
  	degree (opt) : if given, computes a minimal presentation of this homological degree first
  	clear:bool, removes temporary files if true
  	degree: computes minimal presentation of this degree if given
  	verbose : bool
  	"""
  	if id is None:
  	id = str(threading.get_native_id())
  	global input_path, output_path, pathes
  	backend = "function_delaunay"
  	_init_external_softwares(requires=[backend])
  	to_write = np.concatenate([point_cloud, function_values.reshape(-1,1)], axis=1)
  	np.savetxt(input_path+id,to_write,delimiter=' ')
  	verbose_arg = "> /dev/null 2>&1" if not verbose else ""
  	degree_arg = f"--minpres {degree}" if degree >= 0 else ""
  	multi_chunk_arg = "--multi-chunk" if multi_chunk else ""
  	if os.path.exists(output_path + id):
  	os.remove(output_path+ id)
  	command = f"{pathes[backend]} {degree_arg} {multi_chunk_arg} {input_path+id} {output_path+id} {verbose_arg} --no-delaunay-compare"
  	if verbose:
  	print(command)
  	os.system(command)
  	slicer._build_from_scc_file(path=output_path+id, shift_dimension=-1 if degree <= 0 else degree-1 )
  	if clear:
  	clear_io(output_path + id, input_path + id)

Thank you for looking into it ! and feel free to ping me if you need help somewhere.

On another topic, such filtration (i.e., multicritical filtrations of the form float + integer as the degree-rips / your core-delaunay / multicover, etc.) will “soon” benefit from some nice performance and interface improvements when #47 will be finished (in particular with the Degree_rips_filtration data-structure).

[odinhg]

Great, thanks @DavidLapous

I now have a fork (link) of rhomboidtiling which compiles in my multipers conda env at least. I added the option for saving scc files so we don't have to do any sed-ing, keeping it a bit cleaner.

I don't think rhomboidtiling includes mpfree. Would it be possible to use io.scc_reduce_from_str_to_slicer (instead of Slicer._build_from_scc_file) on the scc file outputted by rhomboidtiling?

[DavidLapous]

Thank you Odin, the fork is a good idea indeed.

I don't think rhomboidtiling includes mpfree. Would it be possible to use io.scc_reduce_from_str_to_slicer (instead of Slicer._build_from_scc_file) on the scc file outputted by rhomboidtiling?

I need to implement the multi-critical to 1-critical thing before, but yes, that'd be interesting to also have a reduce_degree flag here as well to avoid some copies.

[odinhg]

I need to implement the multi-critical to 1-critical thing before, but yes, that'd be interesting to also have a reduce_degree flag here as well to avoid some copies.

Thanks for the reply, David. I might be mistaken here, but Michael Lesnick remarks in his Notes on Multiparameter Persistence (page 119) that the rhomboid bifiltration is 1-critical, despite the (weakly equivalent) multicover bifiltration being multi-critical.

What I have been doing in some experiments is something like this:

subprocess.run(["rhomboidtiling", str(point_cloud_file), str(output_file), str(dim), str(k_max), method, str(homology_dimension)], capture_output=True, text=True)
slicer = mp.Slicer()
mp.io.scc_reduce_from_str_to_slicer(output_file, slicer, backend="mpfree", dimension=homology_dimension, verbose=False)

[DavidLapous]

Thanks for the reply, David. I might be mistaken here, but Michael Lesnick remarks in his Notes on Multiparameter Persistence (page 119) that the rhomboid bifiltration is 1-critical, despite the (weakly equivalent) multicover bifiltration being multi-critical.

Yep you're right, mb.

What I have been doing in some experiments is something like this [...]

Thanks for your fork, it works on my end (although this generates enormous filtrations).