Refactoring with `StatPersistence`

[jamesotto852]

I've just pushed to the landscape-refactor branch. Its current state reflects all the work I've done to refactor the computations performed by the various geoms that represent persistence data into StatPersistence. There are some outstanding issues in the refactor, I'm going to enumerate them in this thread.

@corybrunson I'd be happy to schedule some time to discuss this live, if you think that would be helpful?

[corybrunson]

@jamesotto852 fabulous. I'll wait for the enumerated issues, then experiment a bit, and after that i'll definitely be keen to talk through everything in real time.

[jamesotto852]

Just to refresh, the intent of this refactoring has been to ensure logic implementing the conversion of point cloud data to a persistence homology is only implemented in a single place: StatPersistence. This separates the TDA computations from computations which convert the persistence homology to the particular representation that is to be plotted -- the conversion to a plot-ready representation happens in the relevant Geom, currently one of GeomPersistence, GeomBarcode, or GeomLandscape.

This lets us simplify the handling of aesthetics. Each of the Stat implementations require a dataset aesthetic and do not accept start and end which would correspond to a precomputed persistence homology. Instead, a precomputed homology can be used with geom_*(stat = "identity") which requires specification of start and end. This is a small breach of convention -- start and end as "aesthetics" used with stat = "identity" feels odd, there is still transformation happening but it's just happening in the Geom* methods.¹

Footnotes

1. Now that I'm writing this out, I'm realizing this could be implemented differently, with two stats that accept 1.) point cloud and 2.) persistence homology data, ultimately outputting the same data with x, y, and other columns needed by the Geom plot methods. This may be worth investigating at some point, I wonder if some of the issues I'll be posting here might be fixed by this strategy. ↩

[jamesotto852]

Basic Examples

Here is a "complete" set of basic usage of the 6 Geom and Stat implementations. Below, I generate and visualize basic point cloud data on S¹:

library("ggplot2"); theme_set(theme_bw(14))
library("patchwork")
library("tibble")
library("ggtda")
#> Loading required package: rlang
library("ripserr")

# Simulate point cloud data -----------------------------------------------
# Unit circle (Betti-1 number = 1)
set.seed(1); angles <- runif(100, 0, 2 * pi)
circle2d <- cbind(x = cos(angles), y = sin(angles)) 

circle2d |> 
  as_tibble() |> 
  ggplot(aes(x, y)) +
  geom_point() +
  coord_fixed()

Using StatPersistence with each Geom*

Below, I visualize the barcode, persistence diagram, and landscape of the unit circle data. These rely on the same methods in StatPersistence to convert the point cloud data to the persistence homology.

# Using StatPersistent to plot from point cloud data ----------------------
df_pointcloud <- tibble::tibble(data = list(circle2d))

p_barcode <-
  ggplot(df_pointcloud, aes(dataset = data, colour = after_stat(dimension))) +
  geom_barcode() 

p_persistence <-
  ggplot(df_pointcloud, aes(dataset = data, colour = after_stat(dimension))) +
  geom_persistence() +
  geom_abline(slope = 1) +
  coord_fixed(xlim = c(0, 2), ylim = c(0, 2))

p_landscape <-
  ggplot(df_pointcloud, aes(dataset = data, colour = after_stat(dimension))) +
  geom_landscape() 

p_barcode / p_persistence / p_landscape

Using StatIdentity with each Geom*

Now we recreate the above visuals, feeding in manually computed persistence data.

# Precompute PH and plot with stat = "identity" ---------------------------
circ.phom <- as.data.frame(vietoris_rips(circle2d))
circ.phom$dimension <- ordered(circ.phom$dimension, levels = c(0, 1))

p_barcode_identity <- 
  ggplot(circ.phom, aes(start = birth, end = death, colour = dimension)) +
  geom_barcode(stat = "identity")

p_persistence_identity <- 
  ggplot(circ.phom, aes(start = birth, end = death, colour = dimension)) +
  geom_persistence(stat = "identity") +
  geom_abline(slope = 1) +
  coord_fixed(xlim = c(0, 2), ylim = c(0, 2))

p_landscape_identity <- 
  dplyr::rename(circ.phom, "start" = "birth", "end" = "death") |>
  ggtda:::landscape_path() |> 
  ggplot(aes(x, y, colour = dimension)) +
  geom_landscape(stat = "identity")

p_barcode_identity / p_persistence_identity / p_landscape_identity

Notes

The above examples illustrate some of the issues with the current implementation.

1. Persistence diagrams show an extra point at (x = 0, y = Inf)
2. StatLandscape$compute_panel() has an additional step, it uses landscape_path() after StatPersistence$compute_panel(). Currently, this drops the dimension column from the persistence homology data.

^{Created on 2025-09-20 with reprex v2.1.1}

[jamesotto852]

GeomLandscape and StatLandscape

Organizing the steps involved in transforming point cloud data to the landscape representation has been challenging. Right now this is how StatLandscape$compute_panel() is defined:

  compute_panel = function(
    data, scales,
    order_by = c("persistence", "start"),
    decreasing = FALSE,
    filtration = "Rips",
    diameter_max = NULL, radius_max = NULL, dimension_max = 1L,
    field_order = 2L,
    engine = NULL,
    n_levels = Inf
  ) { 
    
    # First, get persistence representation
    data <- StatPersistence$compute_panel(
      data, scales,
      order_by,
      decreasing,
      filtration,
      diameter_max,
      radius_max,
      dimension_max,
      field_order,
      engine
    )
    
    # Next, transform to landscape (path) representation
    # TODO -- need to join these results back into all the extra columns in `data`
    landscape_path(data, n_levels)
  }

See that, instead of purely relying on StatPersistence$compute_panel(), I am using landscape_path() to further manipulate the persistence homology data (this helper function is simply refactoring the previous StatPersistence$compute_group() method).

landscape_path() bug:

Right now, landscape_path() is not preserving the dimension column from the output of StatPersistence$compute_panel(). It's been a while since I've looked at how landscape_path() works, I'm not sure how exactly to identify the dimensions of the frontiers correctly -- I feel like {phutil} may be able to help here?

GeomLandscape's aesthetics

Because landscape_path() converts the persistence homology to the frontiers represented with (x,y) coordinates, GeomLandscape accepts the x and y aesthetics, this is different from GeomBarcode and GeomPersistence. I don't think this is a problem, but it's worth pointing out.

Also, we'll need to export landscape_path() so that GeomLandscape can be used with StatIdentity to visualize pre-computed persistence homology data, this code excerpt from above illustrates the point:

p_landscape_identity <- 
  dplyr::rename(circ.phom, "start" = "birth", "end" = "death") |>
  ggtda:::landscape_path() |> 
  ggplot(aes(x, y, colour = dimension)) +
  geom_landscape(stat = "identity")

Alternatively, we could move the use of landscape_path() to GeomLandscape's draw method. This feels like an odd choice as the logic in landscape_path() feels like it belongs in a compute method, but it may resolve some of the issues and result in a more consistent API. This is something I plan to look into, please let me know if you have thoughts.

[corybrunson]

@jamesotto852 this is very helpful. Some initial thoughts:

This separates the TDA computations from computations which convert the persistence homology to the particular representation that is to be plotted -- the conversion to a plot-ready representation happens in the relevant Geom, currently one of GeomPersistence, GeomBarcode, or GeomLandscape.

I was resistant to this early on, but i've come around: It simplifies the internals while setting only very basic and intuitive expectations of the user. I leant too far toward flexibility, i think, and perhaps tried to be too clever. In the same spirit, though, i wonder whether we really need StatBarcode as an alias of StatPersistence; having geom_barcode() default to stat = "persistence" strikes me as more helpful to the user, by making clear that both, or all three (see below), geoms receive the same transformed data.

Alternatively, we could move the use of landscape_path() to GeomLandscape's draw method.

I also think this makes more sense than exporting landscape_path(). Maybe the function could be exported anyway, but i find it difficult to imagine wanting to combine geom_landscape() with any other geom_*() that takes x,y; one might want to annotate(), but that must be done manually. So i don't think it's necessary to give the user direct access to the x,y transformation of the start,end data obtained from StatPersistence (or passed directly). If a user really wants, they can use the function, exported or not (:::), together with geom_path() or geom_point() or whatever they want to do other than plot a standard landscape. Anyway, this would also obviate the need for StatLandscape, further simplifying the API—and, something i still find valuable, reminding users that not all stats and geoms come in neat cognominal pairs.

Regarding your numbered notes:

1. I think this is appropriate. Some existing software also does this, and the way {ggplot2} locates "points at infinity" makes them more evidently beyond the plotting window than other visualizations tend to. We might consider defaulting the "segments extending to infinity" to be depicted with arrowheads, as some other software does; or maybe that's a possibility better raised with {ggplot2} itself.
2. It actually drops the dimension column, not just drops features with infinite (births or) deaths? In that case, it certainly needs to be fixed. I can take a close look soon.

[jamesotto852]

I've just pushed 66016f1 which addresses a lot of these points:

• Removed aliases of StatPersistence
  I think I agree that we don't need these basic aliases. Originally I wrote up that way so we could have different default behavior depending on the ultimate method of visualizing persistence. {ggdensity} has an example of this kind of implementation, see StatHdr and StatHdrLines. But, now that I've done the refactoring I think it is safe to just reference the same StatPersistence for each geom -- if needed, we could always re-introduce the different stats in a non-breaking way in the future. I went ahead and removed the stat_barcode() and stat_landscape() layer functions, if you think we should include them I can add them back but I figured it's cleaner if we don't include them -- I agree with your perspective on not always offering geom + stat pairs.

• Moved landscape_path() to GeomLandscape$draw_group()
  This was necessary to "collapse" the stats into StatPersistence. Now each of the Geom* representations of persistence homology accept the same aesthetics (dataset and start + end). This also fixed the bug where the dimensions were being dropped (now that this happens at the group-level everything works as originally intended). I do think we should export this function, I've found that users really appreciate having access to the "pipeline" that transforms data into the values being plotted.

[corybrunson]

@jamesotto852 do i understand correctly that landscape-refactor is a descendant of StatPersistence-Refactor? It looks that way from the commit history but i thought i'd check.

I agree with exporting landscape_path(), but i (irresolutely) suggest giving it a cols argument or something similar (e.g. pairs) to which the user is expected to pass c(x,y) orc(start,end), and which is passed c(start,end) internally. For example, a user might want to

landscape_path(data = pd, cols = c(birth, death), n_levels = 3)

instead of

names(pd)[2:3] <- c("start", "end")
landscape_path(data = pd, n_levels = 3)

What is the state of the unit tests? Do you have a mental list of unwritten needs?

Conditional on the above, i'd approve / second a merge.

[jamesotto852]

That's correct, it is a descendant, StatPersistence-Refactor can safely be deleted at this point.

I totally agree with the cols argument for landscape_path().

Here are my thoughts on the work that remains for this branch:

1. The default plotting window for geom_landscape() is not determined correctly
   I'm not sure how to fix this. I can't find good documentation on how to "force" the plot limits -- because we're determining the x and y values to plot in the $draw*() methods the plot doesn't render correctly.

   In each of the Geom* for the 3 visualizations of persistence homologies I've tried to fix this with the $setup_data() method, creating "faux" x and y columns from the start and end columns created by StatPersistence
   (note: these x and y columns are only used to determine the plotting window, they are totally ignored in the draw_*() methods). This works very well in GeomBarcode$setup_data() and GeomPersistence$setup_data() but I don't have a good way to do it in GeomLandscape$setup_data(). Do you have thoughts on either how to either derive better x and y values during setup or another strategy to fix the plotting window?

   Here's a basic illustration of the problem, see that I'm using the default plotting window in the top row and manually setting the window in the bottom row (except in the case of geom_barcode():

library("ggplot2"); theme_set(theme_bw(14))
library("patchwork")
library("ggtda")
#> Loading required package: rlang

# Simulate point cloud data -----------------------------------------------
# Unit circle (Betti-1 number = 1)
set.seed(1); angles <- runif(100, 0, 2 * pi)
circle2d <- cbind(x = cos(angles), y = sin(angles)) 
df_pointcloud <- tibble::tibble(data = list(circle2d))

p_barcode <-
  ggplot(df_pointcloud, aes(dataset = data, colour = after_stat(dimension))) +
  geom_barcode() 

p_persistence <-
  ggplot(df_pointcloud, aes(dataset = data, colour = after_stat(dimension))) +
  geom_persistence() +
  geom_abline(slope = 1) 

p_landscape <-
  ggplot(df_pointcloud, aes(dataset = data, colour = after_stat(dimension))) +
  geom_landscape()


p_barcode_lims <-
  ggplot(df_pointcloud, aes(dataset = data, colour = after_stat(dimension))) +
  geom_barcode() 

p_persistence_lims <-
  ggplot(df_pointcloud, aes(dataset = data, colour = after_stat(dimension))) +
  geom_persistence() +
  geom_abline(slope = 1) +
  coord_fixed(xlim = c(0, 2), ylim = c(0, 2))

p_landscape_lims <-
  ggplot(df_pointcloud, aes(dataset = data, colour = after_stat(dimension))) +
  geom_landscape() +
  coord_fixed(xlim = c(0, 2), ylim = c(0, 2))

p_barcode      + p_persistence      + p_landscape +
p_barcode_lims + p_persistence_lims + p_landscape_lims +
  plot_layout(ncol = 3, guides = "collect")

^{Created on 2025-09-21 with reprex v2.1.1}

2. Dealing with points at infinity in GeomPersistence
   I see that what I thought was a bug was the intended functionality. I'm fine with this default behavior, I do think having a computed variable that could be mapped to would be a nice way to let users customize how these points are represented (e.g. one could map the shape of these points to arrows to achieve what you described).

3. Function-level documentation
   I haven't been keeping up with the {roxygen2} documentation, once we're totally aligned on the implementation and it's finalized I will go through and make sure everything is up-to-date.

4. Long-form documentation
   This refactoring has enabled some really cool, powerful ways to compute and plot persistent homologies, especially when faceting is leveraged. I have some examples that I would like to include as a vignette (and possibly in the README). This would also be an opportunity to fully document how the dataset aesthetic is to be used.

5. Unit Testing
   I think I need to do a overhaul of the current unit testing. I have a decent understanding of where to start, I may need your help to be sure I'm appropriately testing the high-dimensional space of parameters governing the way StatPersistence computes the persistence homology. Similarly, I may need help to cover all the common ways persistence is visualized.

[corybrunson]

1. This is not a problem with the refactor; it occurs in main as well. I don't know whether we've discussed this explicitly, but my outlook has been that, in wrapper functions like those of {TDAstats}, it makes sense to incorporate coord_fixed() into the ggplot output, but when building plots from scratch as here, it's OK to require the extra step. Though i agree that it "feels wrong". {ggplot2} was recently overhauled, so i'm definitely not sure there's not an appropriate trick for defaulting to a different Coord within a layer call like geom_persistence(), similar to how a Position can be defaulted to. Sounds like a good question for the extenders group discussion!
2. This is another one i'd like to ask the extenders group about. So far as i can tell, no such behavior is supported.
3. Copy.
4. Please! Though we should also carefully review the existing vignettes (and examples) under the refactor.
5. Feel free to tag me on any tests i should review. Some of these parameters might be renamed ahead of submission to CRAN, e.g. radius_max and diameter_max should probably be changed from *_max to max_*. (Note: The arguments should agree with those of {tdarec}; if the ones we adopt here differ from those there, then those there should be soft-deprecated in favor of those adopted here.) I also wonder if we should replace order_by with the order aesthetic, which could default to interaction(after_stat(persistence), start) (or similar if that's buggy).

[jamesotto852]

1. I've read through the ggplot2 4.0.0 announcement a few days ago, I don't remember anything that stands out as immediately relevant but it's definitely possible there could be new functionality to help here. I agree that this would be a good point to discuss with the extender group!

   If we can't find a fix, I agree that requiring the extra step is fine, for now. I didn't show it, but the problem also comes up with different diagram choices. I wonder if there's a way we could trigger a warning if the plotting window is incorrect?

2. I think it would actually be really simple -- at the end of StatPersistence$compute_panel() we would just need to add a column to data classifying the type of feature (i.e. infinite end) which could be mapped to in the aes() call with after_stat().

3. .

4. Definitely! I think the current examples + vignettes "Visualizing persistent homology data" and "Illustrating topological constructions on point clouds" are really great and just need minor updating to keep up with changes to the API. I'm primarily thinking about updating "Topological analysis of grouped data and list-columns", leveraging some of the improvements to showcase the nice integration into the {ggplot2} tools, not so much focused on walking through TDA theory.

5. Will do! I definitely agree that we should match the {tdarec} arguments, I can make a first pass to align everything. I need to review the use of the order aesthetic but the change you propose seems reasonable.

   Note, {ggplot2} v4.0.0 is deprecating (or heralding the deprecating) some of the functionality we use in testing. I think I know what needs to change, I'll tag you where it's relevant.

---

Here's the prioritization I would propose, taking into account dependencies in the work to be done:

1. Standardize parameters to align with {tdarec},
2. Investigate order aesthetic
3. Update function-level documentation and examples
4. Round out unit testing
5. Update relevant code in README, "Visualizing persistent homology data" and "Illustrating topological constructions on point clouds"
6. Major update to "Topological analysis of grouped data and list-columns"

A second list of quick fixes which don't warrant inclusion in the prioritization are:

• Fixing the plotting windows
• Adding the computed variable for points at infinity.

---

What do you think, am I missing any work here or should things be moved around? I feel like I can make good progress on these items this week, I'm really hoping we can get this PR wrapped up by next weekend.

[corybrunson]

The prioritization makes sense! Nothing comes to mind, but of course i'll walk through everything closer to the merge. I'll pose those questions to the discussion group soon unless you want to.

[jamesotto852]

Sounds good! I'm happy for you to lead the discussion with the group, although I'd love to be there to participate. Do you have a scheduled meeting in mind?

[corybrunson]

Not yet; i'll just ask the questions on the forum. Though i think we should volunteer to present {ggtda} at a future virtual meetup. (Is that what you meant by a scheduled meeting?)

[jamesotto852]

I agree that we should volunteer to present, maybe right before submitting to CRAN in case we get feedback we'd like to implement? I wasn't thinking about discussing these items then, though — I'd rather do it sooner. It makes sense to ask on the forum to find the right date.