### oxicuda-tda TODO

Topological Data Analysis for OxiCUDA — a pure Rust toolkit for persistent homology,
Vietoris-Rips filtrations, persistence diagrams, the Mapper algorithm, witness
complexes, and standard TDA metrics (Betti numbers, persistent entropy, persistence
landscapes). Part of [OxiCUDA](https://github.com/cool-japan/oxicuda) (Vol.49).

(C) 2026 COOLJAPAN OU (Team KitaSan)

## Implementation Status

**Actual: 53 source files** — implements the standard ELZ-2002 persistence reduction
over Z_2 (plus a host-side chunk-parallel reducer that is bit-identical to it), exact
bottleneck and Wasserstein-1 distances, landscape distance, Mapper with single-linkage
clustering, lazy-witness complex with maxmin landmark sampling, two-parameter persistence
(rank invariant / Hilbert function), Stable Mapper bootstrap, and 7 + 3 PTX kernels
× 6 SM versions.

### Completed [x]

#### Core Infrastructure
- [x] `error.rs` — `TdaError` enum (13 variants: `EmptyPointCloud`, `DimensionMismatch`,
  `InvalidSimplex`, `ClosureViolation`, `EmptyComplex`, `FiltrationNotSorted`,
  `ReductionFailed`, `NanFiltrationValue`, `InvalidCoverParameter`,
  `LandmarkSelectionFailed`, `MatchingFailed`, `ParameterOutOfRange`,
  `DimensionTooLarge`) plus `TdaResult<T>`
- [x] `handle.rs` — `SmVersion`, `LcgRng` (MMIX 64-bit LCG, bit-32 boolean trick),
  `TdaHandle` central entry point
- [x] `ptx_kernels.rs` — 7 kernels × 6 SM versions (sm_70/75/80/86/89/90):
  `pairwise_dist`, `filtration_sort`, `boundary_reduce`, `diagram_match`, `witness_dist`,
  `betti_count`, `mapper_cluster`; all generated by string concatenation

#### Simplicial Complexes & Filtration
- [x] `complex/simplex.rs` — `Simplex { vertices: Vec<usize> }` (always sorted):
  `dim()`, `faces()`, `boundary()` returning (±1, face) pairs (verifies ∂² = 0),
  `contains_face`
- [x] `complex/complex.rs` — `SimplicialComplex` with `add_simplex`,
  `add_simplex_with_closure`, automatic closure verification, indexed by dimension
- [x] `complex/filtration.rs` — `FilteredSimplex { simplex, value }` and `Filtration`:
  `vietoris_rips(dist, n_pts, max_radius, max_dim)` (0-simplices at value 0;
  edges / triangles / higher faces at the max pairwise distance among vertices),
  `vietoris_rips_from_points`, `sublevel_set`

#### Distance Matrices
- [x] `distance/pairwise.rs` — `pairwise_euclidean_sq`, `pairwise_euclidean`,
  `pairwise_manhattan`, `knn_graph`, `points_to_distance_matrix`

#### Persistent Homology
- [x] `homology/boundary.rs` — `BoundaryMatrix` with sparse Z_2 columns
  (`Vec<Vec<usize>>`): `from_filtration`, `low()`, `add_cols()` (XOR), `is_zero()`
- [x] `homology/reduction.rs` — ELZ-2002 standard column-reduction algorithm with
  pivot tracking via `HashMap<row, col>` over Z_2
- [x] `homology/persistent.rs` — `PersistencePair { dim, birth, death: Option<f64> }`
  and `extract_persistence_pairs` joining a reduced boundary matrix with its filtration

#### Diagrams, Barcodes, Distances
- [x] `persistence/diagram.rs` — `PersistenceDiagram { pairs, dimension }`: finite /
  essential filtering, `from_pairs_by_dim`
- [x] `persistence/barcode.rs` — `Bar { birth, death, dim }`, `Barcode`: `from_diagram`,
  `lifetimes`, `count_significant`, `n_finite`
- [x] `persistence/distance.rs` — `bottleneck_distance` (binary search + bipartite
  matching, exact), `wasserstein_1` (Hungarian O(n³)), diagonal projection distance

#### Mapper Algorithm
- [x] `mapper/mapper.rs` — `MapperGraph { nodes, edges }`: overlapping interval cover,
  Union-Find single-linkage clustering per interval, edges = clusters sharing a point;
  `betti_1` = edges − nodes + components; BFS `connected_components`

#### Witness Complex
- [x] `witness/witness.rs` — `maxmin_landmarks` (greedy farthest-point sampling),
  `lazy_witness_complex` with parameter R and 0th-nearest-neighbour distance m_w offset

#### Topological Metrics
- [x] `metrics/metrics.rs` — `betti_numbers`, `persistent_entropy`
  H = -Σ (l_i / L) · log(l_i / L), `persistence_landscape` L_k(t) = k-th largest tent
  function, `landscape_distance` (L2), `total_persistence`, `count_components`

#### Tests & Benchmarks
- [x] `e2e_tests.rs` — 19 cross-module tests: ∂² = 0 sanity, Vietoris-Rips at r = 0
  and r = ∞, filtration sorted-by-value, 3-point H_1 = 1, boundary-matrix reduction,
  birth ≤ death invariant, 4-point square topology, bottleneck self-distance 0,
  Wasserstein triangle inequality, Mapper on circular point cloud (loop detection),
  maxmin landmarks separation, witness vertices, Betti numbers from diagram,
  entropy non-negativity, landscape positivity, essential-class handling,
  PTX × 6 SM versions
- [x] `benches/tda_ops.rs` — Criterion suite: 7 PTX kernels × 4 SM versions plus a
  Vietoris-Rips 20-point 2D end-to-end bench

### Future Enhancements [ ]

#### P0 — Correctness & Foundation
- [x] `complex/cech.rs` — Exact Čech filtration via minimum enclosing ball / smallest
  circle (Welzl). Currently only Vietoris-Rips is provided.
  (complex/cech.rs -- Čech filtration with minimum-enclosing-ball radius per simplex via deterministic Badoiu-Clarkson MEB)
- [x] `complex/alpha.rs` — Alpha filtration via Delaunay triangulation
  (planar + 3D) for sparser but topologically equivalent complexes
- [x] `homology/cohomology.rs` — Cohomology reduction (transpose of boundary matrix)
  giving identical pairs but often a faster reduction path
- [x] `homology/clearing.rs` — Clearing optimisation (Chen-Kerber 2011): skip columns
  whose pivots are already known from a paired (apparent) simplex (implemented inside
  `cohomology.rs` as the `use_clearing` flag in `CohomologyConfig`)
- [x] `homology/twist.rs` — Twist reduction (Chen-Kerber): process columns by
  decreasing dimension to enable additional clearing
  (homology/twist.rs -- Chen-Kerber 2011; standard mod-2 persistence reduction processed by ascending dimension with clearing of paired columns from next dimension; same pairs as standard reduction)

#### P1 — Algorithmic Coverage
- [x] `homology/cohomology_z.rs` — Persistence over Z (with rank trick) for orientable
  manifold detection — current implementation is Z_2 only
  (homology/cohomology_z.rs -- de Silva 2011 cohomology duality; reverse-order reduction of the boundary matrix equivalent to cohomology pairs identical to standard persistent homology; per-pair generator columns)
- [x] `persistence/wasserstein_p.rs` — Wasserstein-p distance for general p > 1
  (current implementation is p = 1)
- [x] `persistence/landscape_p.rs` — L^p landscape distance for p ≠ 2
  (persistence/landscape_p.rs -- Bubenik 2015 persistence landscapes (tent functions, k-th largest) + Lᵖ norm/distance via grid trapezoid)
- [x] `mapper/ball_mapper.rs` — Ball Mapper variant: clusters as balls of fixed radius
  around landmark points (Dlotko 2019)
- [x] `mapper/multiscale.rs` — Multiscale Mapper with cover refinement and persistence
  of the Mapper graph
- [x] `witness/strong_witness.rs` — Strong witness complex (de Silva & Carlsson 2004)
  as a stricter alternative to the lazy variant

#### P1 — Vectorisation & Representations
- [x] `persistence/persistence_image.rs` — Persistence images (Adams et al. 2017):
  vectorisation of diagrams via weighted Gaussian smoothing for downstream ML
- [x] `persistence/persistence_silhouette` — Persistence silhouettes (Chazal et al.):
  weighted average of landscape functions (implemented inside `persistence_image.rs`
  as `persistence_silhouette`)
- [x] `vector/betti_curve.rs` — Betti curves as a vector representation of barcodes

#### P2 — Advanced
- [x] `homology/zigzag.rs` — Zigzag persistence (Carlsson & de Silva 2010) for
  diagrams that allow both insertion and deletion of simplices
- [x] `homology/multi_parameter.rs` — Multi-parameter persistence (RIVET interface,
  Hilbert function, fibered barcodes)
  (homology/multi_parameter.rs -- Carlsson-Zomorodian 2009 two-parameter persistence:
  `BiFiltration`/`BigradedSimplex` with closure+bigraded-monotonicity validation,
  `HilbertFunction` (bigraded Betti grid), `MultiParameterPersistence::rank_invariant`
  via two-step ELZ reduction; function-Rips bi-filtration constructor; already present)
- [x] `distance/kernel.rs` — Persistence-based kernels: scale-space kernel
  (Reininghaus 2015), sliced Wasserstein kernel (Carrière 2017)
- [x] `mapper/stable.rs` — Stable Mapper (Carrière & Oudot 2018) with statistical
  bootstrap confidence intervals
  (mapper/stable.rs -- bootstrap resampling of the point cloud, comparing each replicate's
  (β₀, β₁) Mapper signature to the reference; reports stability/component/loop confidence
  driven by the deterministic LcgRng; already present)
- [x] `benches/algo_bench.rs` — Extended algorithm benches on standard datasets
  (torus / sphere / circle parametric surfaces)
  (benches/algo_bench.rs -- full host-side persistence pipeline (VR → boundary → reduce →
  pairs) on S¹/T²/S² deterministic parametric clouds, sequential-vs-chunk-parallel reducer,
  bottleneck/Wasserstein-1 diagram distances, plus gpu_reduction PTX codegen timing.
  Stanford Dragon / MNIST left out: those need bundled external asset files, not CPU-derivable)
- [x] `persistence/tropical.rs` — Tropical geometry encoding (Monod 2019): tropical coordinates of persistence diagrams as Γ-type statistics; vectorisation stable under reordering and permutation; `TropicalCoordinates`
- [x] `homology/gpu_reduction.rs` — GPU persistence reduction (Otter 2017 spectral): parallel column-reduction with chunk-based pivot lookups in shared memory; PTX kernel for batched boundary-matrix reduction; sm_80+ warp-level synchronisation
  (homology/gpu_reduction.rs -- host-side `chunked_parallel_reduce` (bulk-synchronous
  round-based reduction with a HashMap pivot-lookup chunk; output bit-identical to the
  sequential ELZ reducer, cross-checked on random Rips) + `GpuReductionPlan` tiling +
  `batched_column_reduce_ptx` with cp.async owner-chunk staging on sm_80+. Launching the
  PTX on a real device still requires GPU hardware.)

#### P2 — GPU / Architecture-Specific
- [x] PTX kernel for batched parallel column-reduction with chunk-based pivot lookups
  (homology/gpu_reduction.rs -- `batched_column_reduce_ptx`, one synchronous round; emits
  cp.async owner-chunk prefetch on sm_80+; structurally validated across 7 SM versions)
- [x] PTX kernel for parallel Vietoris-Rips edge enumeration with distance thresholding
  (homology/gpu_reduction.rs -- `vietoris_rips_edges_ptx`, strict-upper-triangle thread map,
  per-edge Euclidean length + threshold flag for host stream-compaction)
- [x] PTX kernel for Wasserstein matching via auction algorithm parallel rounds
  (homology/gpu_reduction.rs -- `wasserstein_auction_ptx`, Bertsekas best/second-best net-value
  scan per unassigned person with ε-bid emission)
- [ ] On-device validation of the above three kernels on a real GPU (requires GPU hardware:
  needs PTX/SASS execution under an NVIDIA driver — codegen + host reference are done above)

## Dependencies

| Dependency | Purpose | Pure Rust? |
|------------|---------|------------|
| thiserror | Error derive macros | Yes |
| criterion (dev) | Benchmark harness | Yes |

No CUDA SDK, no `nvcc`, no link-time `libcuda`, no external linear-algebra or
distance-matrix dependencies. Distance and clustering kernels are implemented from
scratch.

## Quality Status

- Warnings: 0 (clippy clean, `--all-features --all-targets -D warnings`)
- Tests: 398 passing (host-side; incl. 2 worked-example reproductions — textbook two-circles
  persistence and Mapper circular-topology recovery); PTX kernel strings validated per SM version
- `unwrap()` calls in production code: 0
- `unsafe` code: 0
- macOS: compiles; GPU integration paths return `UnsupportedPlatform` at runtime
- Linux + NVIDIA driver 525+: GPU paths exercised behind `#[cfg(feature = "gpu-tests")]`

## Performance Targets

| Workload | Size | Target |
|----------|------|--------|
| Pairwise Euclidean distance | n = 4096 in dim ≤ 64 | bandwidth-bound (≥ 90% theoretical) |
| Vietoris-Rips filtration enumeration | n = 1024, max dim = 2 | edges + triangles ≤ 50 ms on sm_80 |
| Boundary matrix column reduction | ~50K simplices | within 2× CPU baseline (PHAT / Ripser) |
| Bottleneck distance | diagrams with ≤ 1000 pairs | exact within 10 ms |
| Wasserstein-1 distance | diagrams with ≤ 500 pairs | exact (Hungarian O(n³)) within 100 ms |
| Mapper construction | n = 4096, cover = 25, overlap = 0.3 | ≤ 100 ms total |

The persistence reduction is currently a CPU-side ELZ implementation; the
`boundary_reduce` PTX kernel is a stub awaiting integration with a GPU column-reduction
pass (chunk-based pivot lookup with shared-memory hashing).

## Notes

- The PTX kernels are deliberately written as string concatenation rather than `format!`
  to keep generated register names stable across SM versions.
- `BoundaryMatrix` stores each column as `Vec<usize>` of row indices (Z_2): sparse
  enough for typical filtrations and compatible with XOR-style column addition.
- Both `bottleneck_distance` and `wasserstein_1` use the standard reference algorithms
  (binary-searched bipartite matching and the Hungarian method respectively) so they
  match every TDA library that follows Edelsbrunner-Harer.
- `Filtration` always stores simplices sorted by `(value, dim, lex)` so the boundary
  matrix can be built directly without resorting.

---

## Architecture-Specific Deepening

### Volta / Turing (sm_70, sm_75)
- [x] PTX `pairwise_dist` uses 32-wide warp reductions for sum-of-squares
- [ ] Verified on T4 / V100 against CPU reference for Vietoris-Rips on 1K points

### Ampere (sm_80, sm_86)
- [x] PTX kernels emit `cp.async` for the candidate-pair buffer prefetch
- [ ] `boundary_reduce` 3-stage software pipeline benchmarked vs. naive load on A100

### Ada (sm_89)
- [x] Reuses Ampere code path with `cp.async` prefetch
- [ ] Investigate L2 persistence policy for the pivot-row hash table during reduction

### Hopper (sm_90)
- [x] PTX strings emit `wgmma` warp-group MMA where kernels are amenable
- [ ] TMA-based asynchronous load of distance-matrix tiles in `pairwise_dist`
- [ ] Distributed shared-memory multi-CTA column reduction for very large complexes

---

## Deepening Opportunities

### Verification Gaps
- [x] ∂² = 0 sanity is enforced on every `Simplex` boundary computation
- [x] Vietoris-Rips correctness at edge cases r = 0 (vertices only) and r = ∞
  (full complex) is covered in `e2e_tests.rs`
- [x] Bottleneck self-distance = 0 and triangle inequality for Wasserstein-1 are
  validated end-to-end
- [ ] Persistence diagram cross-validation against Ripser / Gudhi reference outputs
  on the same point clouds
- [x] Mapper graph stability under bootstrap resampling (Carrière & Oudot 2018)
  (mapper/stable.rs -- `stable_mapper` bootstrap with (β₀, β₁) signature confidence scores;
  tested on circle / two-blob / single-interval cases)

### Implementation Deepening
- [x] `Filtration::vietoris_rips` correctly assigns simplex birth value as the maximum
  pairwise distance among its vertices
- [x] `BoundaryMatrix::add_cols` is XOR over `Vec<usize>` row-index lists (sparse Z_2)
- [x] `bottleneck_distance` uses binary search over candidate distances with bipartite
  matching at each candidate
- [x] Implement clearing + twist optimisations to bring reduction performance closer
  to PHAT / Ripser (homology/twist.rs -- Chen-Kerber 2011 twist+clearing; homology/cohomology.rs
  `use_clearing` flag)
- [x] Implement cohomology reduction as the typically faster path for "thin"
  filtrations (homology/cohomology.rs `persistent_cohomology` + homology/cohomology_z.rs
  reverse-order duality reduction)

### Documentation Gaps
- [x] Each public type carries a doc comment summarising its semantics
- [x] Worked example reproducing the textbook "two circles" persistence diagram
  (worked_examples.rs/two_circles_persistence_recovers_two_loops -- two disjoint noisy
  unit circles → Vietoris-Rips persistent homology recovers 2 persistent H0 components
  and exactly 2 long-lived H1 loops, ≈1.0 persistence each, clearly above the < 0.7
  H0 sampling-noise floor; no spurious H1 noise bars)
- [x] Worked example showing Mapper recovering a circular topology from noisy 2D points
  (worked_examples.rs/mapper_recovers_circle_topology -- 40 noisy circle points,
  x-coordinate filter, Mapper nerve is a single cycle: 1 connected component, β₁ = 1,
  |E| == |V|, every node degree 2)
