implement persistence engines in a companion package

[astamm]

In the spirit of https://github.com/tdaverse/ggtda/blob/main/R/persistence-engines.R, we could implement persistence engine dispatch that will be used by several downstream packages such as ggtda or tdarec.

We need these to take as input either a point cloud or a function.

[corybrunson]

It's not exactly the same question, but this issue on the scopes of the step functions groups data structures on which the same PH calculations can be performed and lists the implementations of the filtrations that could be used.

I'd like to have a spreadsheet (or tabular data set) somewhere that lines up all the feasible combinations of mathematical object (e.g. point cloud), data structure (e.g. distance matrix), filtration (e.g. Vietoris–Rips), and implementation (e.g. TDA::diagRips(library = "GUDHI")). It might be made an invisible object in {phutil} that guides engine deployment. Does that make sense to you?

[astamm]

It makes a lot of sense yes. I might be needing your help though as you seem to have a better view on this than I do. In mathematical objects, I would have either point clouds or sub-levels of a function. For data structure, only distance matrix? As filtrations, Vietoris-Rips, Cech, Alpha, Cubical, Delaunay, Witness, Cover, Tangential.

Note that you call them filtrations, but sometimes there are referred to as complexes and I am not sure I understand this very well. See GUDHI doc for Delaunay for instance:

Delaunay complex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation. The Simplicial complex filtration values can be computed with different filtrations (Delaunay complex, Delaunay Čech complex or Alpha complex)

When the simplicial complex filtration values are computed, it has the same persistent homology as the Čech complex, while being significantly smaller.

Could it be that these (e.g. Rips) are complexes and then you build a filtration upon them using a data structure like the simplex tree (but not necessarily this one although this is the standard)?

Anyway, for {phutil} v0.0.1 that could be submitted today, I will just have the persistence class with the modifications you suggested and maybe the extended one and the *_set() ones.

[corybrunson]

Oh! I'm sorry, we haven't discussed this in detail. A simplicial filtration is a nested sequence of simplicial complexes, e.g. $\varnothing = K_0 \subset K_1 \subset K_2 \subset \cdots \subset K_{n-1} \subset K_n = K$ is a filtration of $K$. Given a point cloud $X$, the Vietoris–Rips filtration comprises the Vietoris–Rips complexes the $V_r$, which includes a 1-simplex between every pair of vertices within distance $r$ and every higher-dimensional simplex whose edges (1-simplices) are present. If $r < s$, then $V_r \subset V_s$, so the Vietoris–Rips complexes form a filtration indexed by the distance $r$. The persistence of a feature is, roughly speaking, the range of $r$ values for which the feature exists in $K_r$—but this requires some disambiguation that is governed by the "elder rule" and handled by the algebra.

My limited understanding is that the alpha and Čech filtrations are different but homologically equivalent (possibly only under specific circumstances), but that the former are much more efficient to compute, for example.

It is possible to construct PH from the actual simplicial filtration directly—i think GUDHI can do this—but if you're starting from a point cloud then it's much more efficient to do the linear algebra directly from the distance matrix, which is what Ripser does.

As for mathematical objects versus data structures, certainly it's not a 1-to-1 correspondence. A point cloud could be represented as a distance matrix (class "dist" in R, e.g. eurodist), a coordinate matrix (like the output of TDA::circleUnif()), or even a weighted "igraph" object where missing edges indicate infinite distance. Meanwhile, in addition to a point cloud, a matrix could represent what i'm calling a lattice (e.g. volcano) suitable for cubical persistent homology.

Here's my understanding of the significance of these different features for the user:

• mathematical object: the theoretical framework used to analyze the data; ties the data to its meaning, context, and how it was collected
• data structure: how the mathematical object is stored and accessed, basically its type or class; constrains the available computational engines
• filtration: the choice of persistent homology calculation to perform; may be informed by the meaning of the data, e.g. a principled reason for preferring the Čech to the Vietoris–Rips filtration for a specific problem
• implementation: the software used to compute the PH using the desired filtration; the user might trust one package over another, or prefer the metadata one provides, or only care about speed but know that a certain algorithm is faster for data with their characteristics

It would take some work to figure out the right way to store the information about what packages and options are available to construct what filtrations for what data structures encoding what mathematical objects; it would not at all be a simple merging of crosswalks, e.g. from objects to structures, from structures to filtrations, and from filtrations to implementations (which is what i just spent a half-hour trying to do). Maybe we could define a class of "engines", analogous to Docker containers, that contain all this as metadata and can output syntactically valid expressions as well as actually compute the PH for a compatible input object?

[astamm]

Nice detailed answer, thanks! It makes me think that this is exactly how I should structure the metadata field of the 'persistence' class where implementation could itself be a list with the engine (e.g. TDA::ripsDiag()) and the parameters used when calling the engine. What do you think?

[corybrunson]

I'd have to see it, but it sounds reasonable; anyway, it sounds like the right approach, even if we have to workshop it.

[corybrunson]

I've put a draft spreadsheet and illustrative function in the engines branch. Maybe the metadata could be populated from this table, but my guess is not—it'll require more delicate documentation, e.g. alphaShapeDiag() requires a coordinate matrix of exactly 3 dimensions and we will probably want to pre-process input before passing it to gridDiag().

Certainly this feature is for a second release!

[corybrunson]

The tidy design principles book has a nice concise chapter on solutions to the engine problem. Ours will necessarily be more complicated but (i think) could be done without any additional machinery. It would rely on dependencies ({ggtda} and {tdarec}) to handle the exported functions correctly, but there'll be some amount of that no matter how we design engine dispatch in {phutil}.

[corybrunson]

@astamm this plan poses a problem we missed earlier: Should {ripserr} (and perhaps other engines) depend on {phutil} for the persistence class, or should {phutil} depend on {ripserr} and other engines for engine deployment?

I see a few ways out of this dilemma:

1. Create a new package for engine deployment rather than using {phutil}. I don't like this idea because it would result in two small but essential packages rather than one.
2. Have engines like {ripserr} return persistent homology using the diagram class from {TDA}, which is structurally similar to the persistence class. I like this idea even less because it introduces an at least implicit dependency on {TDA}.

Do you see a better way?

Maybe (1) will make sense as more shared utilities are identified to fill out the packages. For a similar case, {interplex} is quite small, as it currently serves only to convert between simplicial complex data structures; but it could be enhanced with engine deployment for constructing simplicial comlpexes (rather than filtrations and persistent homology), which would introduce no new dependencies but support {ggtda} and perhaps other future packages under the hood. Similarly, a new package for PH engine deployment might find other uses based on the same dependencies—maybe validation and benchmarking in addition to class conversion.

[astamm]

I am not sure I have the big picture as much as you do @corybrunson. But it seems to me that by creating the TDAverse we implicitly opted in for creating (and maintaining) small packages designed to do well one particular task and to work well together, in line with the tidyverse principles.

Following these principles would suggest to go for your first option. Even more extreme, we could have a separate package for the distances that are currently in {phutil}. That way we would have:

• {phutil}: for the 'persistence' class and related classes (similar to {tibble});
• another package for engine deployment;
• another package for distance computations.

[corybrunson]

@astamm maybe i've developed a more specific outlook in the years since i connected with Raoul, but the big picture is for all of us to shape as we will. : )

But thanks a lot for your thoughts. I'm iffy about separating {phutil} into the class and the distances, since the distances will likely eventually be used by several packages while they depend on none—just right for a common dependency! But let's at least plan to put engine deployment in a separate package. If OK with you, i'll initialize it in a new repo and migrate this issue.

[corybrunson]

If we want to keep up the naming theme with {phutil}, some options for the engine deployment are {phlourish} (as to flourish a tool), {phlagon} (by analogy to container, e.g. Docker), {phacile} (making things easier), and {phetch} (as to fetch the appropriate instrument). : )

[astamm]

I like most of them :-)

[corybrunson]

I vote for {phetch}, since it evokes the get*() naming convention in both base R and Tidyverse.

[corybrunson]

@astamm i know i volunteered to start this package, but you have a stronger sense for the (R6?) objects that would handle the engines. Would you be able to create a template of one, for both of us (and others) to build upon? (You could do it in a new repo for the eventual package or just in a script that i could incorporate into the package repo.)