| Title: | Topological Data Analysis: Simplicial Complex |
|---|---|
| Description: | Provides an implementation of simplicial complexes for Topological Data Analysis (TDA). The package includes functions to compute faces, boundary operators, Betti numbers, Euler characteristic, and to construct simplicial complexes, including Vietoris-Rips, Cech, Alpha, Delaunay, Witness, flood, and (via a Freudenthal triangulation) cubical complexes for grid and image data. It also implements persistent homology, from building filtrations (via a single build_filtration() entry point covering all of the above) to computing persistence diagrams, persistence landscapes, and Wasserstein/bottleneck distances between diagrams, with the aim of helping readers understand the core concepts of computational topology. Methods are based on standard references in persistent homology such as Zomorodian and Carlsson (2005) <doi:10.1007/s00454-004-1146-y>, Chazal and Michel (2021) <doi:10.3389/frai.2021.667963>, and Otter, Porter, Tillmann, Grindrod and Harrington (2017) <doi:10.1140/epjds/s13688-017-0109-5>. |
| Authors: | ChiChien Wang [aut, cre, trl] |
| Maintainer: | ChiChien Wang <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.2.1 |
| Built: | 2026-09-30 14:26:59 UTC |
| Source: | https://github.com/tda-r/simplicialcomplex |
Construct an Alpha Complex
AlphaComplex(points, epsilon = Inf)AlphaComplex(points, epsilon = Inf)
points |
A numeric matrix or data.frame with one point per row
(columns are coordinates); must be in "general position" (see
|
epsilon |
The alpha-complex scale. A simplex is included once its
alpha value is at most |
A list of class "alpha_complex" with elements:
A list of integer vectors, every included simplex (all dimensions, not just the maximal ones - membership in the alpha complex is not simply generated by cliques).
A numeric vector, the alpha value of each simplex in
simplices, in the same order.
Pass the result to as_filtration for the same
list(simplex=, t=) format used by build_filtration
and every persistence function, or just call
build_filtration(points, method = "Alpha", eps_max = epsilon).
points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) alpha_complex <- AlphaComplex(points, epsilon = 1)points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) alpha_complex <- AlphaComplex(points, epsilon = 1)
Convert a flood_complex object to a filtration list
as_filtration(fc)as_filtration(fc)
fc |
A |
A filtration list compatible with boundary_info.
The bottleneck distance between the points of two persistence diagrams in
a given homological dimension, allowing points to be matched to the
diagonal: the infimum, over all matchings, of the largest single
point-to-point distance (Cohen-Steiner, Edelsbrunner and Harer (2007),
"Stability of Persistence Diagrams") - the limit of
wasserstein_distance. Essential (death = Inf)
classes are kept rather than dropped - see
wasserstein_distance's Details for why and how. Computed
exactly by binary search over the (finitely many) candidate distance
values for the smallest one admitting a perfect matching, checked with a
standard augmenting-path bipartite matcher, on the same augmented
assignment problem wasserstein_distance uses (see
augmented_cost_matrix).
bottleneck_distance(df1, df2, dimension = 0, ground = c("Linf", "L2"))bottleneck_distance(df1, df2, dimension = 0, ground = c("Linf", "L2"))
df1, df2
|
Persistence diagram data frames with columns |
dimension |
Homological dimension to compare. |
ground |
Ground metric on the birth-death plane: |
A single non-negative number, the bottleneck distance.
df1 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 2)) df2 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 3)) bottleneck_distance(df1, df2, dimension = 0)df1 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 2)) df2 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 3)) bottleneck_distance(df1, df2, dimension = 0)
Compute the boundary operator for a simplicial complex
boundary(simplices, bound_dim)boundary(simplices, bound_dim)
simplices |
A list of simplices (each a numeric vector). |
bound_dim |
The dimension k of the boundary operator |
A sparse matrix representing .
simplices <- list(c(1, 2), c(3, 4), c(2, 1, 3), c(4, 2)) boundary(simplices, 0)simplices <- list(c(1, 2), c(3, 4), c(2, 1, 3), c(4, 2)) boundary(simplices, 0)
Get the boundary matrix and its reduction information in matrix form
boundary_info(filist, max_dimension = NULL)boundary_info(filist, max_dimension = NULL)
filist |
Filtration list, each element includes simplex and time. |
max_dimension |
Optional maximum homology dimension to eventually
report via |
A list containing the boundary matrix, the last boundary row, the
pivot owner for persistence extraction, and filist - the exact
(possibly max_dimension-restricted) filtration list the other
three elements were computed from, re-tagged with the same
"max_dimension" attribute so extract_persistence_pairs
can auto-detect it too. Always pass THIS filist back into
extract_persistence_pairs, not your original one - see its
Details for why.
points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) filtration <- build_filtration(points, method = "VR", eps_max = 1.2) res <- boundary_info(filtration)points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) filtration <- build_filtration(points, method = "VR", eps_max = 1.2) res <- boundary_info(filtration)
Build a cubical (grid) filtration from an image
build_cubical_filtration(image, superlevel = FALSE)build_cubical_filtration(image, superlevel = FALSE)
image |
A numeric matrix (grid of pixel/voxel values, e.g. a
grayscale image with values in |
superlevel |
If |
A filtration list: one list(simplex = <integer vector of
pixel ids>, t = <numeric>) per simplex, sorted by (t, dimension,
lexicographic order) - the same format as
build_filtration.
# a ring (value 1) around a hole (value 5) around a background (value 9): # the hole is born once the ring closes and dies once its center fills in image <- matrix(c( 9, 9, 9, 9, 9, 9, 1, 1, 1, 9, 9, 1, 5, 1, 9, 9, 1, 1, 1, 9, 9, 9, 9, 9, 9 ), nrow = 5, byrow = TRUE) filtration <- build_cubical_filtration(image) pairs <- persistence_pairs(filtration) pairs[pairs$dim == 1, ] # one H1 bar: birth = 1 (ring closes), death = 5 plot_persistence(pairs)# a ring (value 1) around a hole (value 5) around a background (value 9): # the hole is born once the ring closes and dies once its center fills in image <- matrix(c( 9, 9, 9, 9, 9, 9, 1, 1, 1, 9, 9, 1, 5, 1, 9, 9, 1, 1, 1, 9, 9, 9, 9, 9, 9 ), nrow = 5, byrow = TRUE) filtration <- build_cubical_filtration(image) pairs <- persistence_pairs(filtration) pairs[pairs$dim == 1, ] # one H1 bar: birth = 1 (ring closes), death = 5 plot_persistence(pairs)
Build a filtration from a point cloud, for any of five complex types
build_filtration( points, method = c("VR", "Delaunay", "Alpha", "Cech", "Witness"), eps_max = NULL, landmarks = NULL, nu = 1, max_dimension = NULL )build_filtration( points, method = c("VR", "Delaunay", "Alpha", "Cech", "Witness"), eps_max = NULL, landmarks = NULL, nu = 1, max_dimension = NULL )
points |
A numeric matrix or data.frame, one point per row. |
method |
One of |
eps_max |
Maximum scale. Required for |
landmarks |
For |
nu |
For |
max_dimension |
Optional maximum homology dimension you intend to compute persistence for. |
A filtration list, as described above. When max_dimension
is set, the list carries it as a "max_dimension" attribute (see
above); otherwise the attribute is absent, exactly as before this
parameter existed.
points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) vr_filtration <- build_filtration(points, method = "VR", eps_max = 1.2) alpha_filtration <- build_filtration(points, method = "Alpha") delaunay_filtration <- build_filtration(points, method = "Delaunay") cech_filtration <- build_filtration(points, method = "Cech", eps_max = 0.8) # cap at H1: persistence_pairs()/etc. need no max_dimension of their own vr_capped <- build_filtration(points, method = "VR", eps_max = 1.2, max_dimension = 1) pairs <- persistence_pairs(vr_capped)points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) vr_filtration <- build_filtration(points, method = "VR", eps_max = 1.2) alpha_filtration <- build_filtration(points, method = "Alpha") delaunay_filtration <- build_filtration(points, method = "Delaunay") cech_filtration <- build_filtration(points, method = "Cech", eps_max = 0.8) # cap at H1: persistence_pairs()/etc. need no max_dimension of their own vr_capped <- build_filtration(points, method = "VR", eps_max = 1.2, max_dimension = 1) pairs <- persistence_pairs(vr_capped)
Wraps flood_complex and returns the filtration in the same format as build_filtration: a list of
list(simplex = <integer vector>, t = <numeric>), sorted by (time, dimension, lexicographic order).
The result plugs directly into boundary_info() / extract_persistence_pairs() /
plot_persistence(), as well as into the faster flood_persistence.
build_flood_filtration(points, landmarks, ...)build_flood_filtration(points, landmarks, ...)
points |
A numeric matrix (N x d) point cloud. |
landmarks |
Number of FPS landmarks, or an explicit landmark matrix. |
... |
Passed on to |
A filtration list compatible with boundary_info.
Construct a Cech Complex
CechComplex(points, epsilon)CechComplex(points, epsilon)
points |
A numeric matrix or data.frame with one point per row (columns are coordinates). |
epsilon |
A non-negative numeric radius, used only to bound how many candidate simplices are generated (see Details) - not a promise that every returned simplex individually satisfies the Cech criterion at this scale. |
A list with:
An igraph object: the 1-skeleton of the candidate
graph (edges where distance ).
A list of integer vectors, each the vertex indices of
a maximal clique of network.
points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) cech_complex <- CechComplex(points, epsilon = 0.8)points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) cech_complex <- CechComplex(points, epsilon = 0.8)
Implements Algorithm 21 (CollectG): for every
critical edge (i.e. every edge not used as a DMVF
matching/tree edge - by construction these all have persistence greater
than dmvf$delta, so unlike the book's pseudocode no extra filter
is needed here), unions with the unique tree paths from
and from up to their respective roots.
collect_g(dmvf)collect_g(dmvf)
dmvf |
A result of |
A data.frame(u, v) of the (deduplicated, undirected) edges of the reconstructed graph.
Builds the same point cloud's filtration under several
build_filtration methods and summarizes, side by side, how
expensive each one was to build and what persistent homology it found -
useful for seeing directly how Cech's dimension blow-up, Alpha's
Delaunay-bounded dimension, VR's cheap-but-approximate cliques, and
Witness's landmark subsampling trade off against each other on the same
data (see Otter et al. (2017), Section 5.2, for the underlying
trade-offs).
compare_complexes( points, methods = c("VR", "Delaunay", "Alpha", "Cech", "Witness"), eps_max = NULL, landmarks = NULL, nu = 1 )compare_complexes( points, methods = c("VR", "Delaunay", "Alpha", "Cech", "Witness"), eps_max = NULL, landmarks = NULL, nu = 1 )
points |
A numeric matrix or data.frame, one point per row. |
methods |
Character vector of methods to compare, any subset of
|
eps_max |
Maximum scale, passed to |
landmarks |
Passed to |
nu |
Passed to |
A list with:
A data frame, one row per method, with the number of vertices and simplices, the highest simplex dimension reached, build and persistence-computation time in seconds, and the number of finite/essential persistence pairs found.
A named list of persistence diagram data frames (one per
method, from persistence_pairs), for further comparison
(e.g. with wasserstein_distance/
bottleneck_distance, or plot_persistence).
A named list of the raw filtration lists.
points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1, 0.5, 0.5), ncol = 2, byrow = TRUE) cmp <- compare_complexes(points, methods = c("VR", "Alpha", "Delaunay"), eps_max = 1.2) cmp$summarypoints <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1, 0.5, 0.5), ncol = 2, byrow = TRUE) cmp <- compare_complexes(points, methods = c("VR", "Alpha", "Delaunay"), eps_max = 1.2) cmp$summary
Picks a single build_filtration method, builds the
persistence diagram of each point cloud with it, and compares the two
diagrams with either wasserstein_distance or
bottleneck_distance. Useful for asking "how different are
these two datasets topologically", holding the complex construction fixed
so the comparison is apples-to-apples.
complex_distance( points1, points2, method = c("VR", "Delaunay", "Alpha", "Cech", "Witness"), distance = c("wasserstein", "bottleneck"), dimension = 0, p = 2, ground = NULL, eps_max = NULL, landmarks = NULL, nu = 1 )complex_distance( points1, points2, method = c("VR", "Delaunay", "Alpha", "Cech", "Witness"), distance = c("wasserstein", "bottleneck"), dimension = 0, p = 2, ground = NULL, eps_max = NULL, landmarks = NULL, nu = 1 )
points1, points2
|
Numeric matrices or data.frames, one point per row. |
method |
One of |
distance |
One of |
dimension |
Homological dimension to compare. |
p |
Wasserstein order, used only when |
ground |
Ground metric passed to the chosen distance function
( |
eps_max |
Maximum scale, passed to |
landmarks |
Passed to |
nu |
Passed to |
A list with:
The computed distance (a single number).
Echoed back for reference.
The two persistence diagrams (data frames) the distance was computed from.
set.seed(1) cloud_a <- matrix(rnorm(24), ncol = 2) cloud_b <- matrix(rnorm(24), ncol = 2) + 0.1 complex_distance(cloud_a, cloud_b, method = "VR", distance = "bottleneck", dimension = 0, eps_max = 0.6)set.seed(1) cloud_a <- matrix(rnorm(24), ncol = 2) cloud_b <- matrix(rnorm(24), ncol = 2) + 0.1 complex_distance(cloud_a, cloud_b, method = "VR", distance = "bottleneck", dimension = 0, eps_max = 0.6)
CROCKER = "Contour Realization Of Computed k-dimensional hole Evolution in the Rips complex".
Treats the -th Betti number as a function of two parameters at once.
crocker( point_clouds, dim, method = "VR", eps_max = NULL, n_eps = 50, max_dimension = NULL, n_cores = 1 )crocker( point_clouds, dim, method = "VR", eps_max = NULL, n_eps = 50, max_dimension = NULL, n_cores = 1 )
point_clouds |
A list of length |
dim |
The homology dimension |
method |
Complex type passed through to |
eps_max |
Maximum scale (epsilon) value. If |
n_eps |
Number of proximity values sampled uniformly from |
max_dimension |
Optional cap forwarded to |
n_cores |
Number of cores to use via |
An object of class "crocker": a list with
An n_eps x n numeric matrix; entry [j, i] is of frame at eps_grid[j].
The sampled scale values (length n_eps).
Frame indices, seq_along(point_clouds).
The homology dimension tracked.
A tidy data frame with columns t, epsilon, betti, ready for plot_crocker.
Construct a Delaunay Complex
DelaunayComplex(points)DelaunayComplex(points)
points |
A numeric matrix or data.frame with one point per row (columns are coordinates); must be in "general position". |
A list of class c("delaunay_complex", "alpha_complex");
see AlphaComplex for the element description.
points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) delaunay_complex <- DelaunayComplex(points)points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) delaunay_complex <- DelaunayComplex(points)
Computes the same optimal matching that
wasserstein_distance/bottleneck_distance
reduce to internally, and returns it decoded into a readable table
instead of just the distance - what plot_matching draws.
Every point of both diagrams appears in exactly one row: matched to a
point of the other diagram (type = "real-real"), or matched to the
diagonal, i.e. effectively unmatched (type = "x-diagonal" for a
point of df1, "y-diagonal" for a point of df2).
diagram_matching( df1, df2, dimension, distance = c("wasserstein", "bottleneck"), p = 2, ground = NULL )diagram_matching( df1, df2, dimension, distance = c("wasserstein", "bottleneck"), p = 2, ground = NULL )
df1, df2
|
Persistence diagram data frames with columns |
dimension |
Homological dimension to compare. |
distance |
Which distance's optimal matching to compute,
|
p |
Wasserstein order, used only when |
ground |
Ground metric on the birth-death plane, |
A list with distance (the matching's cost, equal to what
wasserstein_distance/bottleneck_distance
would return), matches (a data frame with columns
x_birth, x_death, x_essential, y_birth, y_death, y_essential,
type; a NA pair on one side means that row's point matched the
diagonal), and X, Y, essential_X, essential_Y
(the two diagrams' points after the same essential-pair capping
wasserstein_distance uses - see its Details - and which of
them were essential before capping).
df1 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 2)) df2 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 3)) diagram_matching(df1, df2, dimension = 0, distance = "wasserstein", p = 2)df1 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 2)) df2 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 3)) diagram_matching(df1, df2, dimension = 0, distance = "wasserstein", p = 2)
of a simplicial complexCompute the Euler characteristic of a simplicial complex
euler_characteristic(simplices, tol)euler_characteristic(simplices, tol)
simplices |
A list of simplices (each a numeric vector). |
tol |
Optional numerical tolerance to pass to |
The Euler characteristic is computed as:
where is the th Betti number, and is the highest dimension of any simplex in the complex.
Interpretation of values:
: Sphere-like surfaces
: Disk-like spaces
: Torus-like or circle-like spaces
: Surfaces with multiple handles or genus
An integer representing the Euler characteristic .
simplices <- list(c(1, 2), c(3, 4), c(2, 1, 3), c(4, 2)) euler_characteristic(simplices, tol=0.1)simplices <- list(c(1, 2), c(3, 4), c(2, 1, 3), c(4, 2)) euler_characteristic(simplices, tol=0.1)
This function extracts the persistence from combining the boundary matrix and its filtration
extract_persistence_pairs(filist, last_1, pivot_owner, max_dimension = NULL)extract_persistence_pairs(filist, last_1, pivot_owner, max_dimension = NULL)
filist |
Filtration list, each element includes simplex and time.
When |
last_1 |
The last 1 row index for each column in boundary matrix (after reduction). |
pivot_owner |
The column index owning the pivot row. |
max_dimension |
Optional maximum homology dimension to report. Set
this to the SAME value passed to |
boundary_info() and extract_persistence_pairs() are two
halves of one computation - last_1/pivot_owner only mean
anything relative to the exact filist boundary_info() used
internally. Whenever max_dimension is involved, always call as:
res <- boundary_info(filtration, max_dimension = k)
pairs <- extract_persistence_pairs(res$filist, res$last_1, res$pivot_owner,
max_dimension = k)
A length mismatch between filist and last_1/pivot_owner
is refused with an error rather than silently producing wrong pairs - see
persistence_pairs for a one-call alternative that cannot run
into this.
A data frame with columns: dimension, birth, and death.
points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) filtration <- build_filtration(points, method = "VR", eps_max = 1.2) res <- boundary_info(filtration) pairs <- extract_persistence_pairs(filtration, res$last_1, res$pivot_owner)points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) filtration <- build_filtration(points, method = "VR", eps_max = 1.2) res <- boundary_info(filtration) pairs <- extract_persistence_pairs(filtration, res$last_1, res$pivot_owner)
Generate all unique faces of a given dimension from simplices
faces(simplices, target_dim)faces(simplices, target_dim)
simplices |
A list of simplices (each a numeric vector). |
target_dim |
The target dimension |
The function generates all possible subsets (combinations) of each simplex, removes duplicates,
and filters them to only include those of length target_dim + 1.
For example, a 2-simplex c(1, 2, 3) has three 1-dimensional faces (edges):
c(1,2), c(1,3), and c(2,3), and three 0-dimensional faces (vertices):
1, 2, and 3.
A list of faces (each a numeric vector) of dimension target_dim.
simplices <- list(c(1, 2), c(3, 4), c(2, 1, 3), c(4, 2)) faces(simplices, target_dim=0)simplices <- list(c(1, 2), c(3, 4), c(2, 1, 3), c(4, 2)) faces(simplices, target_dim=0)
Builds the Flood complex of a point cloud: a Delaunay complex on a set of landmarks whose simplices are filtered by their covering radius with respect to the full point cloud ("flood time").
flood_complex( points, landmarks, max_dimension = NULL, points_per_edge = 30, backend = c("auto", "cpu", "torch"), batch_points = 2^22, delaunay = NULL )flood_complex( points, landmarks, max_dimension = NULL, points_per_edge = 30, backend = c("auto", "cpu", "torch"), batch_points = 2^22, delaunay = NULL )
points |
A numeric matrix (N x d) of witness points. |
landmarks |
Either an integer (number of FPS landmarks) or a numeric matrix (N_l x d) of explicit landmark coordinates. |
max_dimension |
Top dimension of the simplices. Defaults to d. |
points_per_edge |
Grid resolution per simplex edge (accuracy vs. speed trade-off). Defaults to 30, as in flooder. |
backend |
One of |
batch_points |
Maximum number of grid points processed per batch (bounds memory). Defaults to 2^22. |
delaunay |
Optional precomputed Delaunay triangulation of the landmarks: an m x (d+1) integer matrix of 1-based landmark indices. If NULL (default), computed via |
A list of class "flood_complex" with elements
simplices (list of integer vectors, landmark indices),
filtration (numeric vector of flood times),
landmarks (matrix), landmark_indices (or NULL).
Computes persistence pairs from a filtration list with the standard
column-reduction algorithm, but on sparse columns (integer index vectors
over GF(2)) instead of a dense matrix. Produces the same output format as
extract_persistence_pairs(filist, res$last_1, res$pivot_owner) while
scaling to the much larger complexes produced by build_flood_filtration.
flood_persistence(filist, max_dimension = NULL)flood_persistence(filist, max_dimension = NULL)
filist |
A filtration list (from |
max_dimension |
Optional maximum homology dimension to report. When set, one extra dimension is kept internally so dimension- |
A data frame with columns dim, birth, death.
Row-reduce (partial pivoting) a matrix
gauss_jordan_eliminate(M, tol = 1e-08)gauss_jordan_eliminate(M, tol = 1e-08)
M |
A numeric matrix. |
tol |
Pivot values smaller than this (in absolute value) are treated as zero. Defaults to |
Shared elimination step behind ker and im
(and, internally, the zigzag module's own linear solves): partial-pivoting
Gauss-Jordan elimination, stopping early once every row has a pivot.
A list with R (the row-reduced matrix) and pivots
(the column indices where a pivot was found, in row order).
Selects n_lms landmarks from a point cloud via (exact) Farthest-Point Sampling.
Equivalent to flooder::generate_landmarks (which uses an approximate bucket-FPS;
the exact version below gives the same qualitative coverage).
generate_landmarks(points, n_lms, start_idx = 1)generate_landmarks(points, n_lms, start_idx = 1)
points |
A numeric matrix (N x d) point cloud. |
n_lms |
Number of landmarks to sample (<= N). |
start_idx |
Index of the starting point. Defaults to 1 (flooder defaults to index 0, i.e. the same first point). |
A list with landmarks (n_lms x d matrix) and indices (row indices into points).
The classical graph-theory Laplacian - degree matrix minus adjacency matrix - built directly from a complex's 1-skeleton (vertices + edges).
graph_laplacian(simplices)graph_laplacian(simplices)
simplices |
A list of simplices; only the 0-simplices (vertices) and 1-simplices (edges) are used. |
A list with the Laplacian L, the degree matrix D, the adjacency matrix A, and the vertex basis (row/column order).
Computes .
hodge_laplacian(K_simplices, k)hodge_laplacian(K_simplices, k)
K_simplices |
A single complex (list of maximal simplices). |
k |
The dimension. |
A list with the full Laplacian, its down/up pieces, and the k-simplex basis (row/column order) everything is expressed in.
Compute a homology basis from a cycle basis and a boundary basis
homology(Z, B = matrix(numeric(0), nrow = nrow(Z), ncol = 0), tol = NULL)homology(Z, B = matrix(numeric(0), nrow = nrow(Z), ncol = 0), tol = NULL)
Z |
A matrix whose columns form a basis of the cycle space
|
B |
A matrix whose columns form a basis of the boundary space
|
tol |
Numerical tolerance passed to |
: two cycles
represent the same homology class exactly when they differ by a boundary.
This walks the columns of Z in order, greedily keeping any column
that increases the rank of the span accumulated so far (starting from
B's span) - i.e. any cycle that is not already a linear
combination of B and the cycles kept before it. The kept columns
are one representative chain per homology class.
betti_number computes directly from ranks
(rank-nullity), without ever materializing Z, B, or a
homology basis - that is cheaper when only the count is needed. Use
homology() when the actual representative cycles matter (e.g. to
visualize or track a specific hole), not to recompute a Betti number.
As with ker and im, the specific
representative chosen for each class depends on the order cycles in
Z are tested against the growing boundary span, and is not unique.
A matrix with nrow(Z) rows, one column per representative
of a basis of the quotient . If , the
result has 0 columns; ncol() of the result is the Betti
number .
Compute a basis for the image (column space) of a matrix
im(M, tol = 1e-08)im(M, tol = 1e-08)
M |
A numeric matrix (or an object coercible to one, e.g. a sparse
|
tol |
Pivoting tolerance passed to the Gauss-Jordan elimination used
internally. Defaults to |
, the
column space of M. The basis returned is a subset of M's own
columns, specifically, the pivot columns found by Gauss-Jordan
elimination rather than synthetic linear combinations, so each basis
vector is directly interpretable as one of the original columns of
M (e.g. the boundary of one specific simplex).
A matrix with nrow(M) rows, one column per basis vector of
. If the image is trivial (), the result
has 0 columns.
Compute a basis for the kernel (null space) of a matrix
ker(M, tol = 1e-08)ker(M, tol = 1e-08)
M |
A numeric matrix (or an object coercible to one, e.g. a sparse
|
tol |
Pivoting tolerance passed to the Gauss-Jordan elimination used
internally (see |
. The
basis is obtained from Gauss-Jordan elimination of : one basis
vector per free (non-pivot) column, in the usual parametric-solution
construction. If M has 0 rows (the zero map), every
standard basis vector of the domain is in the kernel, so the identity
matrix is returned.
As with any basis, the specific vectors returned are not unique, they
depend on the pivoting order of the elimination, only the number of
columns () is an invariant of M.
A matrix with ncol(M) rows, one column per basis vector of
. If the kernel is trivial (), the result has
0 columns.
build_cubical_filtration to any triangulation:
given the maximal simplices of a simplicial complex and a function defined at its vertices, builds the simplex-wise lower-star filtration .Generalizes build_cubical_filtration to any triangulation:
given the maximal simplices of a simplicial complex and a function defined at its vertices, builds the simplex-wise lower-star filtration .
lower_star_filtration(top_simplices, f)lower_star_filtration(top_simplices, f)
top_simplices |
A list of maximal simplices (each an integer vector of 1-based vertex ids). |
f |
A plain numeric vector giving the function value at each vertex; |
A filtration list: one list(simplex =integer vector of vertex ids, t = numeric) per simplex (every face of every maximal simplex),
sorted by (t, dimension, lexicographic order).
# two triangles sharing an edge, function increasing away from vertex 1 triangles <- list(c(1, 2, 3), c(2, 3, 4)) f <- c(0, 1, 1, 2) filtration <- lower_star_filtration(triangles, f) length(filtration) # 4 vertices + 5 edges + 2 triangles# two triangles sharing an edge, function increasing away from vertex 1 triangles <- list(c(1, 2, 3), c(2, 3, 4)) f <- c(0, 1, 1, 2) filtration <- lower_star_filtration(triangles, f) length(filtration) # 4 vertices + 5 edges + 2 triangles
Implements Algorithm 20 (MorseRecon):
given a triangulated domain and a density function rho that concentrates around a hidden geometric graph ,
computes the "mountain ridges" of - the 1-unstable manifolds of the discrete gradient field after cancelling
vertex-edge persistence pairs with persistence at most delta - as an approximation of .
morse_recon(top_simplices, rho, delta = 0, vertex_coords = NULL)morse_recon(top_simplices, rho, delta = 0, vertex_coords = NULL)
top_simplices |
List of maximal simplices (e.g. triangles) of the ambient 2-complex, as accepted by |
rho |
A plain numeric vector, the density value at each vertex. |
delta |
Persistence threshold used to cancel low-persistence vertex-edge pairs (noise); larger values denoise more aggressively. |
vertex_coords |
Optional n x 2 matrix of vertex coordinates, used only by |
An object of class "morse_recon": a list with
filtration |
the full lower-star filtration of the 2-complex |
dmvf |
the vertex-edge DMVF: |
graph_edges |
data.frame(u, v) - the edges of |
rho, delta, coords
|
the inputs, kept for plotting/inspection |
Implements Algorithm 19 (SimplePersDMVF) together
with the threshold- simplification: vertex-edge persistence pairs with persistence
at most delta are cancelled, leaving a discrete Morse vector field
where every non-critical vertex is matched with the tree edge
connecting it to its component's root, and every root is a critical vertex.
partial_pers_dmvf(filist1, delta = 0)partial_pers_dmvf(filist1, delta = 0)
filist1 |
A filtration list restricted to vertices and edges only (e.g. via |
delta |
Persistence threshold; vertex-edge pairs with persistence |
A list with components:
pers |
data.frame(u, v, t, persistence, type) for every edge, |
tree_edges |
the subset of |
critical_edges |
the rest of |
parent |
named integer vector, |
roots |
integer vector of critical (root) vertices |
delta |
the threshold used |
Converts a persistence diagram data frame (the direct output of
extract_persistence_pairs(), persistence_pairs(), or
flood_persistence()) into its persistence landscape: a collection of
continuous, piecewise-linear functions , obtained by
overlaying the "tent" function of every birth-death pair and, at each time
t, taking the th largest tent value (Bubenik, 2015; Chazal
and Michel, 2021, Section 5.4).
persistence_landscape( df, dimension = 0, k_max = NULL, resolution = 500, t_range = NULL )persistence_landscape( df, dimension = 0, k_max = NULL, resolution = 500, t_range = NULL )
df |
A persistence diagram data frame with columns |
dimension |
Homological dimension to extract the landscape for. |
k_max |
Number of landscape levels |
resolution |
Number of grid points used to discretize each landscape function. |
t_range |
Optional |
A data frame with columns t, k, value: the
discretized landscape functions, one row per (grid point, level) pair.
points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) filtration <- build_filtration(points, method = "VR", eps_max = 1.2) res <- boundary_info(filtration) pairs <- extract_persistence_pairs(filtration, res$last_1, res$pivot_owner) landscape <- persistence_landscape(pairs, dimension = 0)points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) filtration <- build_filtration(points, method = "VR", eps_max = 1.2) res <- boundary_info(filtration) pairs <- extract_persistence_pairs(filtration, res$last_1, res$pivot_owner) landscape <- persistence_landscape(pairs, dimension = 0)
Same standard column-reduction algorithm as
boundary_info() + extract_persistence_pairs(), but each
column is stored as a sorted integer vector over GF(2) instead of a row of
a dense n x n matrix. Memory drops from O(n^2) to O(total number of
non-zeros), which is what makes filtrations with thousands to hundreds of
thousands of simplices (e.g. Flood or large Vietoris-Rips complexes)
feasible. Output is identical to the dense pipeline.
persistence_pairs(filist, max_dimension = NULL)persistence_pairs(filist, max_dimension = NULL)
filist |
Filtration list, each element includes simplex and time. |
max_dimension |
Optional maximum homology dimension to report
( |
A data frame with columns: dim, birth, and death.
points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) filtration <- build_filtration(points, method = "VR", eps_max = 1.2) pairs <- persistence_pairs(filtration) pairs_h0_only <- persistence_pairs(filtration, max_dimension = 0) # or bake the cap into the filtration itself, and drop the argument here: capped <- build_filtration(points, method = "VR", eps_max = 1.2, max_dimension = 0) pairs_h0_only2 <- persistence_pairs(capped)points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) filtration <- build_filtration(points, method = "VR", eps_max = 1.2) pairs <- persistence_pairs(filtration) pairs_h0_only <- persistence_pairs(filtration, max_dimension = 0) # or bake the cap into the filtration itself, and drop the argument here: capped <- build_filtration(points, method = "VR", eps_max = 1.2, max_dimension = 0) pairs_h0_only2 <- persistence_pairs(capped)
Computes the persistent Laplacian .
persistent_laplacian(X_simplices, Y_simplices, q)persistent_laplacian(X_simplices, Y_simplices, q)
X_simplices, Y_simplices
|
Lists of (maximal) simplices, same format used everywhere else in the package. X must be a subcomplex of Y. |
q |
The dimension. |
A list with the full Laplacian, its upper/down pieces, and the q-simplex basis (in row/column order) everything is expressed in.
Plot a CROCKER matrix as a filled contour plot
plot_crocker(cr)plot_crocker(cr)
cr |
An object returned by |
A ggplot2 object: a tile plot of the Betti number over time and scale.
Plot a Persistence Landscape
plot_landscape(landscape_df)plot_landscape(landscape_df)
landscape_df |
Data frame from |
A ggplot2 object with one line per landscape level
.
points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) filtration <- build_filtration(points, method = "VR", eps_max = 1.2) res <- boundary_info(filtration) pairs <- extract_persistence_pairs(filtration, res$last_1, res$pivot_owner) landscape <- persistence_landscape(pairs, dimension = 0) plot_landscape(landscape)points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) filtration <- build_filtration(points, method = "VR", eps_max = 1.2) res <- boundary_info(filtration) pairs <- extract_persistence_pairs(filtration, res$last_1, res$pivot_owner) landscape <- persistence_landscape(pairs, dimension = 0) plot_landscape(landscape)
Visualizes the point correspondence that realizes the Wasserstein or
bottleneck distance between two persistence diagrams in a given
homological dimension: both diagrams' points, overlaid on the same axes,
joined by dashed lines to their matched partner - either a point of the
other diagram, or (for a point left unmatched) its own projection onto
the diagonal birth = death.
plot_matching( df1, df2, dimension, distance = c("wasserstein", "bottleneck"), p = 2, ground = NULL, labels = c("Diagram 1", "Diagram 2") )plot_matching( df1, df2, dimension, distance = c("wasserstein", "bottleneck"), p = 2, ground = NULL, labels = c("Diagram 1", "Diagram 2") )
df1, df2
|
Persistence diagram data frames with columns |
dimension |
Homological dimension to compare. |
distance |
Which distance's optimal matching to visualize,
|
p |
Wasserstein order, used only when |
ground |
Ground metric, |
labels |
Legend labels identifying |
Essential (death = Inf) points are capped exactly the way
wasserstein_distance does (see its Details) so they take
part in the matching instead of being dropped, and are drawn as triangles
at their capped height - marked "essential" in the legend - so they stay
visually distinguishable from genuine finite points that happen to reach
that height. Matches to the diagonal are always drawn to the point's
perpendicular projection ((birth+death)/2, (birth+death)/2)
regardless of ground, since that is the standard, readable way to
depict an unmatched point whichever ground metric produced the matching.
A ggplot2 object; the title reports the resulting distance.
df1 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 2)) df2 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 3)) plot_matching(df1, df2, dimension = 0, distance = "wasserstein", p = 2)df1 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 2)) df2 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 3)) plot_matching(df1, df2, dimension = 0, distance = "wasserstein", p = 2)
A 3D companion to plot_morse_recon: renders the input density rho as an actual terrain surface over the triangulated
domain (using the triangles stored in mr$filtration), and overlays the discrete Morse critical structure on it with the rgl package,
producing a mouse-rotatable 3D scene rather than a flat projection:
the root (critical) vertices of mr$dmvf, the surviving density peaks after persistence simplification - as labelled spheres ("M1", "M2", ...);
the endpoints of mr$dmvf$critical_edges, the surviving "saddle-like" edges of Section 10.4.1 as smaller spheres;
mr$graph_edges itself, i.e. the collected 1-unstable manifolds from collect_g, literally the vertex-to-
saddle connections along the reconstructed ridge lines as raised line segments on the surface.
plot_morse_landscape( mr, z_scale = NULL, show_critical_edges = TRUE, label_roots = TRUE, surface_col = c("#f7fbff", "#6baed6", "#08306b"), root_col = "red", saddle_col = "blue", edge_col = "black", point_radius = NULL, edge_lwd = 3, window_size = c(1400, 1000), title_cex = 1.4 )plot_morse_landscape( mr, z_scale = NULL, show_critical_edges = TRUE, label_roots = TRUE, surface_col = c("#f7fbff", "#6baed6", "#08306b"), root_col = "red", saddle_col = "blue", edge_col = "black", point_radius = NULL, edge_lwd = 3, window_size = c(1400, 1000), title_cex = 1.4 )
mr |
A |
z_scale |
Vertical exaggeration applied to |
show_critical_edges |
Logical; if |
label_roots |
Logical; if |
surface_col, root_col, saddle_col, edge_col
|
Colour for the density surface, the root markers, the saddle markers, and the |
point_radius, edge_lwd
|
Marker radius, auto-scaled from the domain size when |
window_size |
Length-2 |
title_cex |
Character expansion for the |
A interactive 3D rgl scene.
On macOS, opening a native rgl window requires XQuartz; Windows and Linux do not need it. To avoid a native window altogether, render the scene as a WebGL widget in the RStudio Viewer or a browser:
options(rgl.useNULL = TRUE) plot_morse_landscape(mr) rgl::rglwidget()
morse_recon
Draws the reconstructed graph (the union of 1-unstable manifolds of the surviving high-persistence critical edges)
over the input vertices, optionally coloured by the input density rho.
plot_morse_recon(mr, show_density = TRUE, point_size = 0.6, edge_size = 0.9)plot_morse_recon(mr, show_density = TRUE, point_size = 0.6, edge_size = 0.9)
mr |
A |
show_density |
Logical; if |
point_size, edge_size
|
Point/line sizes passed to ggplot2. |
A ggplot2 object.
Figure zooms into a small neighbourhood of the mesh (as most DMT papers illustrate the vector field, since drawing it over the whole domain is unreadable) and draws the actual computed structure on top of the local triangulation:
every triangulation edge with both endpoints in the window, in light grey, for context;
the surviving critical edges of mr$dmvf$critical_edges that fall inside the window, as thick red segments.
These are the "saddles" that did not get matched to a vertex;
for every non-critical vertex in the window whose DMVF match mr$dmvf$parent[[v]] is also in the window, a black arrow
from to the midpoint of its matched edge , the vertex-edge V-path arrows themselves.
plot_morse_vpath( mr, center, radius, vertex_size = 1.2, arrow_size = 0.12, critical_size = 1.3 )plot_morse_vpath( mr, center, radius, vertex_size = 1.2, arrow_size = 0.12, critical_size = 1.3 )
mr |
A |
center |
Length-2 numeric, an |
radius |
Euclidean radius (same units as |
vertex_size, arrow_size, critical_size
|
Point size, arrow-head size, and line width for the critical edges, respectively. |
A ggplot2 object.
Plot Persistence Diagram
plot_persistence(df)plot_persistence(df)
df |
Dataframe from plot_persistence. |
A ggplot2 object representing the persistence diagram.
points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) filtration <- build_filtration(points, method = "VR", eps_max = 1.2) res <- boundary_info(filtration) pairs <- extract_persistence_pairs(filtration, res$last_1, res$pivot_owner) plot_persistence(pairs)points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) filtration <- build_filtration(points, method = "VR", eps_max = 1.2) res <- boundary_info(filtration) pairs <- extract_persistence_pairs(filtration, res$last_1, res$pivot_owner) plot_persistence(pairs)
Drops every simplex of dimension greater than max_dimension from a
filtration list; everything else (order, ties, t values) is left
untouched.
restrict_filtration(filist, max_dimension)restrict_filtration(filist, max_dimension)
filist |
A filtration list, as produced by
|
max_dimension |
Maximum simplex dimension to keep ( |
A filtration list containing only the entries with dimension
<= max_dimension, in the same relative order.
points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) filtration <- build_filtration(points, method = "VR", eps_max = 1.2) restrict_filtration(filtration, max_dimension = 1) # vertices + edges onlypoints <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) filtration <- build_filtration(points, method = "VR", eps_max = 1.2) restrict_filtration(filtration, max_dimension = 1) # vertices + edges only
Construct a Vietoris–Rips Complex (1-skeleton + maximal simplices)
VietorisRipsComplex(points, epsilon)VietorisRipsComplex(points, epsilon)
points |
A numeric matrix or data.frame with one point per row (columns are coordinates). |
epsilon |
A positive numeric threshold; connect points with distance < |
The Vietoris–Rips complex at scale includes a simplex for every
finite set of points with pairwise distances < . This function
constructs the 1-skeleton (edges only) and then uses maximal cliques in that
graph as the maximal simplices.
A list with:
An igraph object representing the 1-skeleton.
A list of integer vectors, each the vertex indices of a maximal simplex.
points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) epsilon <- 1.5 vr_complex <- VietorisRipsComplex(points, epsilon)points <- matrix(c(0, 1, 1, 0, 0, 0, 1, 1), ncol = 2) epsilon <- 1.5 vr_complex <- VietorisRipsComplex(points, epsilon)
The -Wasserstein distance between the points of two persistence
diagrams in a given homological dimension, allowing points to be matched
to the diagonal (Cohen-Steiner, Edelsbrunner, Harer and Mileyko (2010),
"Lipschitz Functions Have -Stable Persistence"). Essential
(death = Inf) classes are kept rather than dropped - see Details.
Computed exactly via the Hungarian algorithm (clue::solve_LSAP) on
the augmented assignment problem described in
augmented_cost_matrix.
wasserstein_distance(df1, df2, dimension = 0, p = 2, ground = c("L2", "Linf"))wasserstein_distance(df1, df2, dimension = 0, p = 2, ground = c("L2", "Linf"))
df1, df2
|
Persistence diagram data frames with columns |
dimension |
Homological dimension to compare. |
p |
Wasserstein order ( |
ground |
Ground metric on the birth-death plane used to measure the
distance between two points (or a point and the diagonal): |
A single non-negative number, the -Wasserstein distance.
df1 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 2)) df2 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 3)) wasserstein_distance(df1, df2, dimension = 0, p = 2)df1 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 2)) df2 <- data.frame(dim = c(0, 0), birth = c(0, 0), death = c(1, 3)) wasserstein_distance(df1, df2, dimension = 0, p = 2)
Construct a Witness Complex
WitnessComplex(points, landmarks, epsilon, nu = 1)WitnessComplex(points, landmarks, epsilon, nu = 1)
points |
A numeric matrix or data.frame of witness points |
landmarks |
Either an integer (the number of landmarks to select from
|
epsilon |
A non-negative numeric scale. |
nu |
Landmark-distance parameter; defaults to |
A list with:
An igraph object: the 1-skeleton on the landmark
vertex set.
A list of integer vectors (indices into landmarks
/ landmark_indices), each the vertices of a maximal simplex.
The landmark coordinate matrix.
Row indices into points, or NULL
if landmarks was given as coordinates rather than indices.
A symmetric matrix giving, for every pair of
landmarks, the exact scale at which their edge is witnessed (used by
build_filtration to time every simplex exactly, rather
than only checking membership at the single scale epsilon).
points <- matrix(rnorm(200), ncol = 2) witness_complex <- WitnessComplex(points, landmarks = 15, epsilon = 0.5)points <- matrix(rnorm(200), ncol = 2) witness_complex <- WitnessComplex(points, landmarks = 15, epsilon = 0.5)
Internally, each step is expanded into a sequence of single simplex insertions (faces before cofaces)
or single simplex deletions (cofaces before faces), exactly as boundary_info/persistence_pairs
assume for ordinary filtrations. At every elementary insertion or deletion, the representative-cycle basis of each affected
homology dimension is updated directly via linear algebra (reusing boundary and faces for every
boundary-matrix computation), an insertion either creates a new cycle (birth in dimension ) or turns an existing
cycle into a boundary (death in dimension ); a deletion either destroys an existing cycle (death in dimension )
or frees a previously-trivial cycle from being a boundary (birth in dimension ), where is the
dimension of the simplex being inserted/deleted.
zigzag_persistence(complexes, max_dimension = NULL)zigzag_persistence(complexes, max_dimension = NULL)
complexes |
A list of length |
max_dimension |
Optional integer cap: dimensions above this are not tracked
(saves work for large complexes where only e.g. |
A data frame with columns dim, birth, death (integer
indices into 0, ..., n, i.e. into complexes), one row per bar. Every
bar uses a closed interval: death is the last index at which
the class is still present, and a class still alive at is reported with
death = n, since a zigzag filtration, unlike an ordinary one, has no
canonical "infinity" to extend to. This differs by one from
persistence_pairs's convention, where death is the index of the
killing simplex and the interval is half-open (birth <= i < death);
on a purely-growing filtration the two agree after death_here = death_there - 1.