From d2961fa1fcb8d9bb8be2e82c33a6e4e4c4866a7c Mon Sep 17 00:00:00 2001
From: mtsch <matijacufar@gmail.com>
Date: Fri, 26 Jun 2020 18:51:49 +1200
Subject: [PATCH] New Way of Getting Simplices From Filtrations (#60)

* new way of getting simplices from filtrations

* sparse filtration performance improvement
---
 .travis.yml                 |   2 +-
 Project.toml                |   2 +-
 docs/src/api/filtrations.md |  14 +-
 docs/src/api/simplices.md   |   8 -
 src/Ripserer.jl             |   2 +-
 src/abstractfiltration.jl   |  70 +++++---
 src/abstractsimplex.jl      | 157 +++++++++++-------
 src/cubical.jl              | 177 +++++++++++---------
 src/main.jl                 |  30 ++--
 src/primefield.jl           |   2 +-
 src/reductionmatrix.jl      |  65 ++++----
 src/ripsfiltration.jl       | 310 +++++++++++++++++++++++++-----------
 src/simplex.jl              |  95 +----------
 src/zerodimensional.jl      |  12 +-
 test/cubical.jl             |  54 +++----
 test/reductionmatrix.jl     |  23 ++-
 test/ripserer.jl            |  26 ++-
 test/ripsfiltration.jl      |  88 +++-------
 test/simplex.jl             |  68 ++++----
 test/testhelpers.jl         |   9 +-
 20 files changed, 663 insertions(+), 551 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index 110c7ed3..95b3ca60 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -5,7 +5,7 @@ os:
 
 julia:
   - 1.0
-  - 1.4
+  - 1
   - nightly
 
 notifications:
diff --git a/Project.toml b/Project.toml
index 62dc9e87..476f65d6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Ripserer"
 uuid = "aa79e827-bd0b-42a8-9f10-2b302677a641"
 authors = ["mtsch <matijacufar@gmail.com>"]
-version = "0.11.0"
+version = "0.12.0"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
diff --git a/docs/src/api/filtrations.md b/docs/src/api/filtrations.md
index d6a0ac0c..8ee92dce 100644
--- a/docs/src/api/filtrations.md
+++ b/docs/src/api/filtrations.md
@@ -13,23 +13,27 @@ Ripserer.edges(::Ripserer.AbstractFiltration)
 ```
 
 ```@docs
-Ripserer.diam(::Ripserer.AbstractFiltration, ::Any)
+Ripserer.birth(::Ripserer.AbstractFiltration, ::Any)
 ```
 
 ```@docs
-Ripserer.diam(::Ripserer.AbstractFiltration, ::Ripserer.AbstractSimplex, ::Any, ::Any)
+Ripserer.threshold(::Ripserer.AbstractFiltration)
 ```
 
 ```@docs
-Ripserer.birth(::Ripserer.AbstractFiltration, ::Any)
+Ripserer.simplex_type
 ```
 
 ```@docs
-Ripserer.threshold(::Ripserer.AbstractFiltration)
+Ripserer.simplex
 ```
 
 ```@docs
-Ripserer.simplex_type
+Ripserer.unsafe_simplex
+```
+
+```@docs
+Ripserer.unsafe_cofacet
 ```
 
 ```@docs
diff --git a/docs/src/api/simplices.md b/docs/src/api/simplices.md
index 8510e2d1..91f6502e 100644
--- a/docs/src/api/simplices.md
+++ b/docs/src/api/simplices.md
@@ -20,14 +20,6 @@ Base.sign(::Ripserer.AbstractSimplex)
 Base.:-(::Ripserer.AbstractSimplex)
 ```
 
-```@docs
-Ripserer.coface_type(::Ripserer.AbstractSimplex)
-```
-
-```@docs
-Ripserer.face_type(::Ripserer.AbstractSimplex)
-```
-
 ```@docs
 vertices(::Ripserer.AbstractSimplex)
 ```
diff --git a/src/Ripserer.jl b/src/Ripserer.jl
index c22bbbd8..49bf827a 100644
--- a/src/Ripserer.jl
+++ b/src/Ripserer.jl
@@ -24,7 +24,7 @@ using StaticArrays
 using TupleTools
 
 # reexporting basics these makes Ripserer usable without having to import another package.
-import PersistenceDiagrams: birth, threshold
+import PersistenceDiagrams: birth, threshold, dim
 export
     birth, death, persistence, representative, birth_simplex, death_simplex, barcode
 
diff --git a/src/abstractfiltration.jl b/src/abstractfiltration.jl
index 245fe5e0..43ec49af 100644
--- a/src/abstractfiltration.jl
+++ b/src/abstractfiltration.jl
@@ -8,14 +8,13 @@ simplices.
 
 * [`n_vertices(::AbstractFiltration)`](@ref)
 * [`edges(::AbstractFiltration)`](@ref)
-* [`diam(::AbstractFiltration, vertices)`](@ref)
-* [`diam(::AbstractFiltration, ::AbstractSimplex, ::Any, ::Any)`](@ref) - only used when
-  [`simplex_type`](@ref) is an [`IndexedSimplex`](@ref).
 * [`simplex_type(::AbstractFiltration, dim)`](@ref)
-* [`birth(::AbstractFiltration, v)`](@ref) - optional, defaults to returning `zero(T)`.
-* [`threshold(::AbstractFiltration)`](@ref) - optional, defaults to returning `Inf`.
-* [`postprocess_interval(::AbstractFiltration, ::Any)`](@ref) - optional
-  postprocessing function that is applied to each interval in resulting persistence diagram.
+* [`simplex(::AbstractFiltration, ::Val{dim}, vertices, sign)`](@ref)
+* [`unsafe_simplex(::AbstractFiltration, ::Val{dim}, vertices, sign)`](@ref)
+* [`unsafe_cofacet`](@ref)`(::AbstractFiltration, simplex, vertices, vertex[, edges, sign])`
+* [`birth(::AbstractFiltration, v)`](@ref)
+* [`threshold(::AbstractFiltration)`](@ref)
+* [`postprocess_interval(::AbstractFiltration, ::Any)`](@ref)
 """
 abstract type AbstractFiltration end
 
@@ -43,28 +42,59 @@ n_vertices(::AbstractFiltration)
 """
     edges(filtration::AbstractFiltration)
 
-Get edges in distance matrix in `filtration`, sorted by decresing length and increasing
-combinatorial index. Edges should be of type [`simplex_type`](@ref)(filtration, 1)`.
+Get edges (1-simplices) in `filtration`. Edges should be of type
+[`simplex_type`](@ref)`(filtration, 1)`.
 """
 edges(::AbstractFiltration)
 
 """
-    diam(::AbstractFiltration, vertices)
+     simplex(::AbstractFiltration, ::Val{D}, vertices, sign=1)
 
-Get the diameter of a simplex with the vertex set `vertices`. Should return `missing` if
-`vertices` do not form a valid simplex.
+Return `D`-simplex constructed from `vertices` with sign equal to `sign`. Return `nothing`
+if simplex is not in filtration. This function is safe to call with vertices that are out of
+order. Default implementation sorts `vertices` and calls [`unsafe_simplex`](@ref).
 """
-diam(::AbstractFiltration, ::Any)
+function simplex(flt::AbstractFiltration, ::Val{D}, vertices, sign=1) where D
+    vxs = TupleTools.sort(Tuple(vertices), rev=true)
+    if allunique(vxs) && all(x -> x > 0, vxs)
+        return unsafe_simplex(flt, Val(D), vxs, sign)
+    else
+        throw(ArgumentError("invalid vertices $(vertices)"))
+    end
+end
+
+"""
+    unsafe_simplex(::AbstractFiltration, ::Val{D}, vertices, sign=1)
+
+Return `D`-simplex constructed from `vertices` with sign equal to `sign`. Return `nothing`
+if simplex is not in filtration. The unsafe in the name implies that it's up to the caller
+to ensure vertices are sorted and unique.
+"""
+unsafe_simplex(::AbstractFiltration, ::Val, vertices, sign)
 
 """
-    diam(::AbstractFiltration, simplex, vertices, new_vertex)
+    unsafe_cofacet(filtration, simplex, cofacet_vertices, new_vertex[, edges, sign=1])
+
+Return cofacet of `simplex` with vertices equal to `cofacet_vertices`. `new_vertex` is the
+vertex that was added to construct the cofacet. In the case of sparse rips filtrations, an
+additional argument `edges` is used. `edges` is a vector that contains the weights on edges
+connecting the new vertex to old vertices.
 
-Get the diameter of coface of a `Simplex` that is formed by adding `new_vertex` to
-`vertices`. Should return `missing` if new simplex is not valid.
+The unsafe in the name implies that it's up to the caller to ensure vertices are sorted and
+unique.
 
-This functions is used with the [`coboundary`](@ref) function for [`IndexedSimplex`](@ref)
+Default implementation uses [`unsafe_simplex`](@ref).
 """
-diam(::AbstractFiltration, ::AbstractSimplex, ::Any, ::Any)
+function unsafe_cofacet(
+    flt::AbstractFiltration, ::AbstractSimplex{D}, vertices, v, sign=1
+) where D
+    return unsafe_simplex(flt, Val(D + 1), vertices, sign)
+end
+function unsafe_cofacet(
+    flt::AbstractFiltration, ::AbstractSimplex{D}, vertices, v, edges::SVector, sign=1
+) where D
+    return unsafe_simplex(flt, Val(D + 1), vertices, sign)
+end
 
 """
     birth(::AbstractFiltration, v)
@@ -77,12 +107,12 @@ birth(::AbstractFiltration, _) = false # false is a strong zero.
     threshold(::AbstractFiltration)
 
 Get the threshold of filtration. This is the maximum diameter a simplex in the filtration
-can have. Used only for placing the infinity line in plotting. Defaults to `missing`.
+can have. Defaults to `Inf`.
 """
 threshold(::AbstractFiltration) = Inf
 
 """
-    postprocess_diagram(::AbstractFiltration, interval)
+    postprocess_interval(::AbstractFiltration, interval)
 
 This function is called on each resulting persistence interval. The default implementation
 does nothing.
diff --git a/src/abstractsimplex.jl b/src/abstractsimplex.jl
index 047203a0..50c5cedf 100644
--- a/src/abstractsimplex.jl
+++ b/src/abstractsimplex.jl
@@ -16,13 +16,35 @@ actually needed for the main algorithm.
 * [`Base.:-(::AbstractSimplex)`](@ref)
 * `Base.isless(::AbstractSimplex, ::AbstractSimplex)`
 * [`vertices(::AbstractSimplex)`](@ref)
-* [`coface_type(::AbstractSimplex)`](@ref)
 * [`coboundary(::Any, ::AbstractSimplex)`](@ref)
-* [`face_type(::AbstractSimplex)`](@ref) - only required for homology.
-* [`boundary(::Any, ::AbstractSimplex)`](@ref) - only required for homology.
+* [`boundary(::Any, ::AbstractSimplex)`](@ref)
+* `length(::Type{AbstractSimplex})` - only if `D`-simplexs does not have `D + 1` vertices.
 """
 abstract type AbstractSimplex{D, T, I} <: AbstractVector{I} end
 
+Base.:(==)(::AbstractSimplex, ::AbstractSimplex) = false
+function Base.:(==)(sx1::A, sx2::A) where A<:AbstractSimplex
+    return diam(sx1) == diam(sx2) && vertices(sx1) == vertices(sx2)
+end
+function Base.isequal(sx1::A, sx2::A) where A<:AbstractSimplex
+    return diam(sx1) == diam(sx2) && vertices(sx1) == vertices(sx2)
+end
+Base.hash(sx::AbstractSimplex, h::UInt64) = hash(vertices(sx), hash(diam(sx), h))
+
+Base.iterate(sx::AbstractSimplex) = iterate(vertices(sx))
+Base.getindex(sx::AbstractSimplex, i) = vertices(sx)[i]
+Base.firstindex(::AbstractSimplex) = 1
+Base.lastindex(sx::AbstractSimplex) = length(sx)
+Base.size(sx::AbstractSimplex) = (length(sx),)
+
+Base.length(sx::AbstractSimplex) = length(typeof(sx))
+Base.length(::Type{<:AbstractSimplex{D}}) where D = D + 1
+
+function Base.show(io::IO, sx::AbstractSimplex{D}) where D
+    print(io, nameof(typeof(sx)), "{", D, "}(",
+          sign(sx) == 1 ? :+ : :-, vertices(sx), ", ", diam(sx), ")")
+end
+
 """
     diam(simplex::AbstractSimplex)
 
@@ -70,49 +92,6 @@ Reverse the simplex orientation.
 Base.:-(::AbstractSimplex)
 Base.:+(sx::AbstractSimplex) = sx
 
-Base.:(==)(::AbstractSimplex, ::AbstractSimplex) = false
-function Base.:(==)(sx1::A, sx2::A) where A<:AbstractSimplex
-    return diam(sx1) == diam(sx2) && vertices(sx1) == vertices(sx2)
-end
-function Base.isequal(sx1::A, sx2::A) where A<:AbstractSimplex
-    return diam(sx1) == diam(sx2) && vertices(sx1) == vertices(sx2)
-end
-Base.hash(sx::AbstractSimplex, h::UInt64) = hash(vertices(sx), hash(diam(sx), h))
-
-"""
-    coface_type(::AbstractSimplex)
-    coface_type(::Type{<:AbstractSimplex})
-
-Get the type of simplex's coface. For a `D`-dimensional simplex, this is usually its
-`D+1`-dimensional counterpart. Only the method for the type needs to be implemented.
-
-```jldoctest
-coface_type(Simplex{2}((3, 2, 1), 3.2))
-
-# output
-
-Simplex{3, Float64, Int}
-```
-"""
-coface_type(sx::AbstractSimplex) = coface_type(typeof(sx))
-
-"""
-    face_type(::AbstractSimplex)
-    face_type(::Type{<:AbstractSimplex})
-
-Get the type of a simplex's face. For a `D`-dimensional simplex, this is usually its
-`D-1`-dimensional counterpart. Only the method for the type needs to be implemented.
-
-```jldoctest
-face_type(Simplex{2}((3, 2, 1), 3.2))
-
-# output
-
-Simplex{1, Float64, Int}
-```
-"""
-face_type(sx::AbstractSimplex) = face_type(typeof(sx))
-
 """
     dim(::AbstractSimplex)
     dim(::Type{<:AbstractSimplex})
@@ -135,8 +114,7 @@ Base.abs(sx::AbstractSimplex) = sign(sx) == 1 ? sx : -sx
 """
     vertices(simplex::AbstractSimplex{dim})
 
-Get the vertices of `simplex`. Returns `NTuple{dim+1, Int}`. In the algorithm, only the
-method for 2-simplices is actually used.
+Get the vertices of `simplex`. Returns `SVector{length(simplex), Int}`.
 
 ```jldoctest
 vertices(Simplex{2}((3, 2, 1), 3.2))
@@ -153,12 +131,13 @@ vertices(Simplex{2}((3, 2, 1), 3.2))
 vertices(::AbstractSimplex)
 
 """
-    coboundary(filtration, simplex[, Val{all_cofaces}])
+    coboundary(filtration, simplex[, Val{all_cofacets}])
 
 Iterate over the coboundary of `simplex`. Use the `filtration` to determine the diameters
-and validity of cofaces. Iterates values of the type [`coface_type`](@ref)`(simplex)`. If
-`all_cofaces` is `false`, only return cofaces with vertices added to the beginning of vertex
-list.
+and validity of cofacets. If `all_cofacets` is `false`, only return cofaces with vertices
+added to the beginning of vertex list.
+
+Comes with a default implementation.
 
 ```jldoctest coboundary
 filtration = Rips([0 1 1 1; 1 0 1 1; 1 1 0 1; 1 1 1 0])
@@ -184,13 +163,17 @@ end
   +[4, 3, 1]
 ```
 """
-coboundary(::Any, ::AbstractSimplex)
+function coboundary(filtration, simplex::AbstractSimplex, ::Val{A}=Val(true)) where A
+    return Coboundary{A}(filtration, simplex)
+end
 
 """
-    boundary(filtration, simplex[, Val{all_cofaces}])
+    boundary(filtration, simplex[, Val{all_cofacets}])
 
 Iterate over the boundary of `simplex`. Use the `filtration` to determine the diameters
-and validity of cofaces. Iterates values of the type [`face_type`](@ref)`(simplex)`.
+and validity of cofacets.
+
+Comes with a default implementation.
 
 ```jldoctest boundary
 filtration = Rips([0 1 1 1; 1 0 1 1; 1 1 0 1; 1 1 1 0])
@@ -206,4 +189,66 @@ Simplex{2}(-[4, 1], 1)
 Simplex{2}(+[4, 2], 1)
 ```
 """
-boundary(::Any, ::AbstractSimplex)
+boundary(filtration, simplex::AbstractSimplex) = Boundary(filtration, simplex)
+
+# (co)boundaries ========================================================================= #
+struct Coboundary{A, D, I, F, S}
+    filtration::F
+    simplex::S
+    vertices::NTuple{D, I}
+
+    function Coboundary{A}(
+        filtration::F, simplex::AbstractSimplex{D, <:Any, I}
+    ) where {A, D, I, F}
+        S = typeof(simplex)
+        return new{A, D + 1, I, F, S}(filtration, simplex, Tuple(vertices(simplex)))
+    end
+end
+
+function Base.iterate(
+    ci::Coboundary{A, D, I}, (v, k)=(I(n_vertices(ci.filtration) + 1), D),
+) where {A, D, I}
+    @inbounds while true
+        v -= one(I)
+        while k > 0 && v == ci.vertices[end + 1 - k]
+            A || return nothing
+            v -= one(I)
+            k -= 1
+        end
+        v > 0 || return nothing
+        sign = ifelse(iseven(k), one(I), -one(I))
+        new_vertices = TupleTools.insertafter(ci.vertices, D - k, (v,))
+        sx = unsafe_cofacet(ci.filtration, ci.simplex, new_vertices, v, sign)
+        if !isnothing(sx)
+            _sx::simplex_type(ci.filtration, D) = sx
+            return _sx, (v, k)
+        end
+    end
+end
+
+struct Boundary{D, I, F, S}
+    filtration::F
+    simplex::S
+    vertices::NTuple{D, I}
+
+    function Boundary(
+        filtration::F, simplex::AbstractSimplex{D, <:Any, I}
+    ) where {D, I, F}
+        S = typeof(simplex)
+        return new{D + 1, I, F, S}(filtration, simplex, Tuple(vertices(simplex)))
+    end
+end
+
+function Base.iterate(bi::Boundary{D, I}, k=1) where {D, I}
+    while k ≤ D
+        facet_vertices = TupleTools.deleteat(bi.vertices, k)
+        k += 1
+        sign = ifelse(iseven(k), one(I), -one(I))
+        sx = unsafe_simplex(bi.filtration, Val(D - 2), facet_vertices, sign)
+        if !isnothing(sx)
+            _sx::simplex_type(bi.filtration, D - 2) = sx
+            return _sx, k
+        end
+    end
+    return nothing
+end
diff --git a/src/cubical.jl b/src/cubical.jl
index 2e2c5d47..53e99712 100644
--- a/src/cubical.jl
+++ b/src/cubical.jl
@@ -36,27 +36,15 @@ end
     end
 end
 
-function Base.show(io::IO, ::MIME"text/plain", cube::Cubelet{D, T, I}) where {D, T, I}
-    print(io, D, "-dim Cubelet", (index(cube), diam(cube)))
-    print(io, ":\n  $(sign(cube) == 1 ? '+' : '-')$(vertices(cube))")
-end
-
-function Base.show(io::IO, cube::Cubelet{D, M}) where {D, M}
-    print(io, "Cubelet{$D}($(sign(cube) == 1 ? '+' : '-')$(vertices(cube)), $(diam(cube)))")
-end
-
 index(cube::Cubelet) = cube.index
 diam(cube::Cubelet) = cube.diam
-coface_type(::Type{Cubelet{D, T, I}}) where {D, T, I} = Cubelet{D + 1, T, I}
-face_type(::Type{Cubelet{D, T, I}}) where {D, T, I} = Cubelet{D - 1, T, I}
 
 # @generated used because inference doesn't work well for 2^D
 @generated function vertices(cube::Cubelet{D, <:Any, I}) where {D, I}
     return :(vertices(index(cube), Val($(2^D))))
 end
 
-Base.lastindex(::Cubelet{D}) where D = 2^D
-Base.size(::Cubelet{D}) where D = (2^D,)
+Base.length(::Type{<:Cubelet{D}}) where D = 1 << D
 
 """
     Cubical{I<:Signed, T} <: AbstractFiltration
@@ -87,53 +75,79 @@ birth(cf::Cubical, i) = cf.data[i]
 simplex_type(::Cubical{I, T}, dim) where {I, T} = Cubelet{dim, T, I}
 
 dim(cf::Cubical) = length(size(cf.data))
-Base.CartesianIndices(cf::Cubical) = CartesianIndices(cf.data)
-Base.LinearIndices(cf::Cubical) = LinearIndices(cf.data)
 
-function diam(cf::Cubical{<:Any, T}, vertices) where {T}
-    res = typemin(T)
+function to_linear(cf::Cubical{I}, vertices) where I
+    indices = LinearIndices(cf.data)
+    return map(v -> I(get(indices, v, 0)), vertices)
+end
+function to_cartesian(cf::Cubical, vertices)
+    indices = CartesianIndices(cf.data)
+    return map(v -> indices[v], vertices)
+end
+
+function unsafe_simplex(
+    cf::Cubical{I, T}, ::Val{D}, vertices, sign=one(I)
+) where {I, T, D}
+    diam = typemin(T)
     for v in vertices
-        if isnothing(get(cf.data, v, nothing))
-            return missing
+        d = get(cf.data, v, missing)
+        if ismissing(d)
+            return nothing
         else
-            res = max(res, cf.data[v])
+            _d::T = d
+            diam = ifelse(_d > diam, _d, diam)
         end
     end
-    return res
+    return simplex_type(cf, D)(sign * index(vertices), diam)
+end
+
+# Check if linear indices u and v are adjacent in array arr.
+function is_adjacent(arr, u, v)
+    if isassigned(arr, u) && isassigned(arr, v)
+        diff = Tuple(CartesianIndices(arr)[u] - CartesianIndices(arr)[v])
+        ones = abs.(diff) .== 1
+        zeros = diff .== 0
+        return count(ones) == 1 && count(zeros) == length(diff) - 1
+    else
+        return false
+    end
 end
 
 function edges(cf::Cubical{I, <:Any}) where {I}
-    E = edge_type(cf)
-    result = E[]
-    for u_lin in eachindex(cf.data)
-        u_car = CartesianIndices(cf)[u_lin]
+    result = edge_type(cf)[]
+    for u in eachindex(cf.data)
         for d in 1:dim(cf)
-            v_car = u_car + CartesianIndex{dim(cf)}(ntuple(i -> i == d ? 1 : 0, dim(cf)))
-            if v_car in CartesianIndices(cf)
-                v_lin = LinearIndices(cf)[v_car]
-                push!(result, E(
-                    index((I(v_lin), I(u_lin))), max(cf.data[u_lin], cf.data[v_lin])
-                ))
+            if d == 1
+                v = u + 1
+            else
+                v = u + prod(size(cf.data, i) for i in 1:d - 1)
+            end
+            if is_adjacent(cf.data, u, v)
+                sx = unsafe_simplex(cf, Val(1), (v, u), 1)
+                if !isnothing(sx)
+                    _sx::edge_type(cf) = sx
+                    push!(result, _sx)
+                end
             end
         end
     end
-    return sort!(result)
+    return result
 end
 
 # coboundary ============================================================================= #
-struct CubeletCoboundary{A, I, N, C<:Cubelet, F<:Cubical{I}, K}
+struct CubeletCoboundary{A, D, N, I, F<:Cubical{I}, K}
     filtration::F
-    cubelet::C
-    vertices::SVector{K, CartesianIndex{N}}
+    vertices::NTuple{K, CartesianIndex{N}}
 end
 
 function CubeletCoboundary{A}(
     filtration::F, cubelet::C
-) where {A, I, D, F<:Cubical{I}, C<:Cubelet{D, <:Any, I}}
-    K = length(cubelet)
+) where {A, I, D, F<:Cubical{I}, C<:Cubelet{D}}
+    K = length(C)
     N = dim(filtration)
-    vxs = map(v -> CartesianIndices(filtration)[v], vertices(cubelet))
-    return CubeletCoboundary{A, I, N, C, F, K}(filtration, cubelet, vxs)
+    return CubeletCoboundary{A, D, N, I, F, K}(
+        filtration, to_cartesian(filtration, Tuple(vertices(cubelet)))
+    )
 end
 
 function coboundary(filtration::Cubical, cubelet::Cubelet)
@@ -151,40 +165,51 @@ function all_equal_in_dim(dim, vertices)
     return true
 end
 
+# Type safe way to get half of tuple a-la TupleTools.
+function second_half(tup::NTuple{N}) where N
+    return _half(tup, Val(N÷2), TupleTools.unsafe_tail)
+end
+function first_half(tup::NTuple{N}) where N
+    return _half(tup, Val(N÷2), TupleTools.unsafe_front)
+end
+@inline function _half(tup, ::Val{N}, f) where N
+    if N == 0
+        return tup
+    else
+        return _half(f(tup), Val(N-1), f)
+    end
+end
+
 function Base.iterate(
-    cc::CubeletCoboundary{A, I, N, C}, (dim, dir)=(one(I), one(I))
-) where {A, I, N, C}
-    # If not all indices in a given dimension are equal, we can't create a coface by
-    # expanding in that direction.
+    cc::CubeletCoboundary{A, D, N, I}, (dim, dir)=(one(I), one(I))
+) where {A, D, N, I}
     while true
+        # Idea: expand cube by adding vertices that have ±1 added to one of the dimensions
+        # where all old vertices have the same value.
         while dim ≤ N && !all_equal_in_dim(dim, cc.vertices)
             dim += one(I)
         end
-        if dim > N
-            break
-        end
-
-        diff = CartesianIndex{N}(ntuple(isequal(dim), N) .* dir)
-        new_vertices = cc.vertices .+ Ref(diff)
-        diameter = max(diam(cc.cubelet), diam(cc.filtration, new_vertices))
+        dim > N && return nothing
 
+        diff = CartesianIndex{N}(ntuple(isequal(dim), Val(N)) .* dir)
+        new_vertices = TupleTools.sort(to_linear(
+            cc.filtration, TupleTools.vcat(cc.vertices, cc.vertices .+ Ref(diff))
+        ), rev=true)
         dim += I(dir == -1)
         dir *= -one(I)
-
-        if !ismissing(diameter)
-            all_vertices = sort(map(v -> I(LinearIndices(cc.filtration)[v]),
-                                    vcat(cc.vertices, new_vertices)), rev=true)
-            if A || all_vertices[end÷2+1:end] == LinearIndices(cc.filtration)[cc.vertices]
-                return coface_type(C)(-dir * index(all_vertices), diameter), (dim, dir)
+        if A || second_half(new_vertices) == to_linear(cc.filtration, cc.vertices)
+            sx = unsafe_simplex(cc.filtration, Val(D + 1), new_vertices, -dir)
+            if !isnothing(sx)
+                _sx::simplex_type(cc.filtration, D + 1) = sx
+                return _sx, (dim, dir)
             end
         end
     end
 end
 
-struct CubeletBoundary{D, I, C<:Cubelet, N, F<:Cubical{I}, K}
+struct CubeletBoundary{D, I, N, F<:Cubical{I}, K}
     filtration::F
-    cubelet::C
-    vertices::SVector{K, CartesianIndex{N}}
+    vertices::NTuple{K, CartesianIndex{N}}
 end
 
 function CubeletBoundary(
@@ -192,34 +217,34 @@ function CubeletBoundary(
 ) where {D, I, F<:Cubical{I}, C<:Cubelet{D}}
     K = length(cubelet)
     N = dim(filtration)
-    vxs = map(v -> CartesianIndices(filtration)[v], vertices(cubelet))
-    return CubeletBoundary{D, I, C, N, F, K}(filtration, cubelet, vxs)
+    return CubeletBoundary{D, I, N, F, K}(
+        filtration, to_cartesian(filtration, Tuple(vertices(cubelet)))
+    )
 end
 
 boundary(filtration::Cubical, cubelet::Cubelet) = CubeletBoundary(filtration, cubelet)
 
-function first_half(vec::SVector{K, T}) where {K, T}
-    SVector{K÷2, T}(Tuple(vec)[1:K÷2])
-end
-function second_half(vec::SVector{K, T}) where {K, T}
-    SVector{K÷2, T}(Tuple(vec)[K÷2+1:K])
-end
-
-# Idea: split the cube in each dimension, returning two halves depending on dir.
-function Base.iterate(cb::CubeletBoundary{D, I, C}, (dim, dir)=(1, 1)) where {D, I, C}
+# Idea: split the cube in each dimension, returning one of two halves depending on dir.
+function Base.iterate(cb::CubeletBoundary{D, I}, (dim, dir)=(1, 1)) where {D, I}
     if dim > D
         return nothing
     else
-        lin_vertices = map(v -> I(LinearIndices(cb.filtration)[v]),
-                           sort(cb.vertices, by=v -> v[dim], rev=true))
+        lin_vertices = to_linear(
+            cb.filtration, TupleTools.sort(cb.vertices, by=v -> v[dim], rev=true)
+        )
         if dir == 1
             new_vertices = first_half(lin_vertices)
-            diameter = diam(cb.filtration, new_vertices)
-            return face_type(C)(new_vertices, diameter), (dim, -dir)
+            sign = one(I)
+            state = (dim, -dir)
         else
             new_vertices = second_half(lin_vertices)
-            diameter = diam(cb.filtration, new_vertices)
-            return -face_type(C)(new_vertices, diameter), (dim + 1, 1)
+            sign = -one(I)
+            state = (dim + 1, 1)
+        end
+        sx = simplex(cb.filtration, Val(D - 1), new_vertices, sign)
+        if !isnothing(sx)
+            _sx::simplex_type(cb.filtration, D - 1) = sx
+            return _sx, state
         end
     end
 end
diff --git a/src/main.jl b/src/main.jl
index 1b729d12..7f28b828 100644
--- a/src/main.jl
+++ b/src/main.jl
@@ -1,28 +1,23 @@
 """
     overflows(filtration::AbstractFiltration, dim_max)
+    overflows(::Type{AbstractSimplex}, n_vertices, field)
 
 Check if `dim`-dimensional simplex with largest index on `n_vertices(filtration)` can safely
-be constructed. Throw `OverflowError` error if it can't.
+be constructed and if it can be represented as a chain element.
 """
 function overflows(flt::AbstractFiltration, dim_max, field)
     return overflows(simplex_type(flt, dim_max + 1), n_vertices(flt), field)
 end
 
-"""
-    overflows(::Type{AbstractSimplex}, n_vertices, field)
-
-Check if simplex with largest index on `n_vertices` can safely be constructed with
-overflow. Throw `OverflowError` error if it would overflow.
-"""
 overflows(::Type{<:AbstractSimplex}, ::Any, ::Any) = false
 
 function overflows(S::Type{<:IndexedSimplex{<:Any, T, I}}, n_vertices, field) where {T, I}
-    len = length(S(1, oneunit(T)))
-    len > n_vertices && throw(ArgumentError("$S has more than $(n_vertices) vertices."))
+    length(S) > n_vertices && throw(ArgumentError("$S has more than $(n_vertices) vertices."))
+    # Idea: calculate index in I and BigInt and compare if they are always the same.
     acc_int = zero(I)
     acc_big = big(0)
     i = n_vertices
-    for k in len:-1:1
+    for k in length(S):-1:1
         acc_int += small_binomial(I(i), Val(k))
         acc_big += binomial(big(i), big(k))
         acc_int ≠ acc_big && return true
@@ -98,19 +93,16 @@ function ripserer(
     cohomology=true,
 )
     if overflows(filtration, dim_max, field_type)
+        S = simplex_type(filtration, dim_max + 1)
         throw(OverflowError(
-            "Constructing some of the $dim_max-simplices in this filtration would overflow. " *
+            "$S on $(n_vertices(filtration)) vertices overflows." *
             "Try using a larger index type or a smaller `dim_max`."
         ))
     end
-    if cohomology
-        return _cohomology(
-            filtration, cutoff, progress, field_type, Val(dim_max), Val(reps)
-        )
+    return if cohomology
+        _cohomology(filtration, cutoff, progress, field_type, Val(dim_max), Val(reps))
     else
-        return _homology(
-            filtration, cutoff, progress, field_type, Val(dim_max), Val(reps)
-        )
+        _homology(filtration, cutoff, progress, field_type, Val(dim_max), Val(reps))
     end
 end
 
@@ -133,7 +125,6 @@ function _cohomology(
                 matrix = next_matrix(matrix, progress)
             end
         end
-
         return result
     end
 end
@@ -165,6 +156,5 @@ function _homology(
             matrix = next_matrix(matrix, progress)
         end
     end
-
     return result
 end
diff --git a/src/primefield.jl b/src/primefield.jl
index 8a5bb7f9..a4cb9f57 100644
--- a/src/primefield.jl
+++ b/src/primefield.jl
@@ -19,7 +19,7 @@ Return `true` if `n` is a prime number.
     end
 end
 # The following causes an error if you try to use anything other than Int as a modulus.
-@pure is_prime(::Any) = false
+is_prime(::Any) = false
 
 """
     mod_prime(i, ::Val{M})
diff --git a/src/reductionmatrix.jl b/src/reductionmatrix.jl
index 67cfa0b3..fb74c256 100644
--- a/src/reductionmatrix.jl
+++ b/src/reductionmatrix.jl
@@ -206,31 +206,32 @@ end
 # TODO: the following code is pretty messy. There must be a way to handle all this with a
 # bit less type magic.
 struct ReductionMatrix{
-    Cohomology, Field, Filtration, Simplex, SimplexElem, Face, FaceElem, O<:Base.Ordering
+    Cohomology, T, Filtration, Simplex, SimplexElem, Facet, FacetElem, O<:Base.Ordering
 }
     filtration::Filtration
-    reduced::ReducedMatrix{Face, SimplexElem, O}
-    working_boundary::WorkingBoundary{FaceElem, O}
+    reduced::ReducedMatrix{Facet, SimplexElem, O}
+    working_boundary::WorkingBoundary{FacetElem, O}
     columns_to_reduce::Vector{Simplex}
     columns_to_skip::Vector{Simplex}
 end
 
-function ReductionMatrix{Co, Field}(
+function ReductionMatrix{C, T}(
     filtration::Filtration, columns_to_reduce, columns_to_skip
-) where {Co, Field, Filtration}
+) where {C, T, Filtration}
 
     Simplex = eltype(columns_to_reduce)
-    Face = Co ? coface_type(Simplex) : face_type(Simplex)
-    ordering = Co ? Base.Order.Forward : Base.Order.Reverse
+    facet_dim = dim(Simplex) + (C ? 1 : -1)
+    Facet = simplex_type(filtration, facet_dim)
+    ordering = C ? Base.Order.Forward : Base.Order.Reverse
     O = typeof(ordering)
-    SimplexElem = chain_element_type(Simplex, Field)
-    FaceElem = chain_element_type(Face, Field)
+    SimplexElem = chain_element_type(Simplex, T)
+    FacetElem = chain_element_type(Facet, T)
 
-    reduced = ReducedMatrix{Face, SimplexElem}(ordering)
+    reduced = ReducedMatrix{Facet, SimplexElem}(ordering)
     sizehint!(reduced, length(columns_to_reduce))
-    working_boundary = WorkingBoundary{FaceElem}(ordering)
+    working_boundary = WorkingBoundary{FacetElem}(ordering)
 
-    return ReductionMatrix{Co, Field, Filtration, Simplex, SimplexElem, Face, FaceElem, O}(
+    return ReductionMatrix{C, T, Filtration, Simplex, SimplexElem, Facet, FacetElem, O}(
         filtration,
         reduced,
         working_boundary,
@@ -252,11 +253,11 @@ function co_boundary(matrix::ReductionMatrix{false}, simplex::AbstractSimplex)
     return boundary(matrix.filtration, simplex)
 end
 
-is_cohomology(::ReductionMatrix{Co}) where Co = Co
+is_cohomology(::ReductionMatrix{C}) where C = C
 field_type(::ReductionMatrix{<:Any, F}) where F = F
 simplex_type(::ReductionMatrix{<:Any, <:Any, <:Any, S}) where S = S
 simplex_element(::ReductionMatrix{<:Any, <:Any, <:Any, <:Any, E}) where E = E
-face_element(::ReductionMatrix{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, E}) where E = E
+facet_element(::ReductionMatrix{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, E}) where E = E
 dim(::ReductionMatrix{true, <:Any, <:Any, S}) where S = dim(S)
 dim(::ReductionMatrix{false, <:Any, <:Any, S}) where S = dim(S) - 1
 
@@ -264,16 +265,16 @@ function initialize_boundary!(matrix::ReductionMatrix, column_simplex)
     # TODO: can emergent pairs be safely be enabled for all kinds of complexes?
     # An argument that disables it could be added.
     empty!(matrix.working_boundary)
-    for face in co_boundary(matrix, column_simplex)
-        # Checking this on every face helps if more than one face has the same diameter.
+    for facet in co_boundary(matrix, column_simplex)
+        # Checking this on every facet helps if more than one facet has the same diameter.
         if (is_cohomology(matrix) &&
-            diam(face) == diam(column_simplex) &&
-            !haskey(matrix.reduced, face)
+            diam(facet) == diam(column_simplex) &&
+            !haskey(matrix.reduced, facet)
             )
             empty!(matrix.working_boundary)
-            return face_element(matrix)(face)
+            return facet_element(matrix)(facet)
         end
-        nonheap_push!(matrix.working_boundary, face)
+        nonheap_push!(matrix.working_boundary, facet)
     end
     if isempty(matrix.working_boundary)
         return nothing
@@ -286,10 +287,10 @@ end
 function add!(matrix::ReductionMatrix, column, factor)
     record!(matrix.reduced, column, factor)
     for element in column
-        for face in co_boundary(matrix, simplex(element))
+        for facet in co_boundary(matrix, simplex(element))
             push!(
                 matrix.working_boundary,
-                face_element(matrix)(face, coefficient(element) * factor)
+                facet_element(matrix)(facet, coefficient(element) * factor)
             )
         end
     end
@@ -405,11 +406,11 @@ simplex_name(::Type{<:Simplex{2}}) = "triangles"
 simplex_name(::Type{<:Simplex{3}}) = "tetrahedra"
 simplex_name(::Type{<:AbstractSimplex{D}}) where D = "$D-simplices"
 
-function next_matrix(matrix::ReductionMatrix{Co, Field}, progress) where {Co, Field}
-    C = coface_type(simplex_type(matrix))
+function next_matrix(matrix::ReductionMatrix, progress)
+    C = simplex_type(matrix.filtration, dim(simplex_type(matrix)) + 1)
     new_to_reduce = C[]
     new_to_skip = C[]
-    Co && sizehint!(new_to_skip, length(matrix.reduced))
+    is_cohomology(matrix) && sizehint!(new_to_skip, length(matrix.reduced))
 
     if progress
         progbar = Progress(
@@ -418,12 +419,12 @@ function next_matrix(matrix::ReductionMatrix{Co, Field}, progress) where {Co, Fi
         )
     end
     for simplex in Iterators.flatten((matrix.columns_to_reduce, matrix.columns_to_skip))
-        for coface in coboundary(matrix.filtration, simplex, Val(false))
+        for cofacet in coboundary(matrix.filtration, simplex, Val(false))
             # Clearing optimization only enabled for cohomology.
-            if Co && haskey(matrix.reduced, coface)
-                push!(new_to_skip, abs(coface))
+            if is_cohomology(matrix) && haskey(matrix.reduced, cofacet)
+                push!(new_to_skip, abs(cofacet))
             else
-                push!(new_to_reduce, abs(coface))
+                push!(new_to_reduce, abs(cofacet))
             end
         end
         progress && next!(progbar)
@@ -432,8 +433,10 @@ function next_matrix(matrix::ReductionMatrix{Co, Field}, progress) where {Co, Fi
         stderr, "Assembled $(length(new_to_reduce)) $(simplex_name(C)). Sorting... ",
         color=:green
     )
-    sort!(new_to_reduce, rev=Co)
+    sort!(new_to_reduce, rev=is_cohomology(matrix))
     progress && printstyled(stderr, "done.\n", color=:green)
 
-    return ReductionMatrix{Co, Field}(matrix.filtration, new_to_reduce, new_to_skip)
+    return ReductionMatrix{is_cohomology(matrix), field_type(matrix)}(
+        matrix.filtration, new_to_reduce, new_to_skip
+    )
 end
diff --git a/src/ripsfiltration.jl b/src/ripsfiltration.jl
index 2ff37d93..c245aa55 100644
--- a/src/ripsfiltration.jl
+++ b/src/ripsfiltration.jl
@@ -1,110 +1,146 @@
-# distance matrix stuff ================================================================== #
 """
-    edges(dist::AbstractMatrix, thresh, E)
-
-Return sorted edges of type `E` in distance matrix with length lower than `thresh`.
-"""
-function edges(dist::AbstractMatrix, thresh, ::Type{E}) where E
-    @boundscheck begin
-        size(dist, 1) == size(dist, 2) && size(dist, 1) > 1 ||
-            throw(ArgumentError("invalid input matrix"))
-    end
-
-    n = size(dist, 1)
-    res = E[]
-    @inbounds for j in 1:n, i in j+1:n
-        l = dist[i, j]
-        l ≤ thresh && push!(res, E((i, j), l))
-    end
-    return sort!(res)
-end
+    AbstractRipsFiltration{I<:Signed, T} <: AbstractFiltration
 
-function edges(dist::AbstractSparseMatrix, thresh, ::Type{E}) where E
-    res = E[]
-    I, J, V = findnz(dist)
-    for (i, j, l) in zip(I, J, V)
-        i > j || continue
-        l ≤ thresh && push!(res, E((i, j), l))
-    end
-    return sort!(res)
-end
+An abstract Vietoris-Rips filtration. Its subtypes can overload [`dist`](@ref) and get the
+following default implementations.
 
+* [`edges`](@ref)
+* [`simplex_type`](@ref).
+* [`simplex`](@ref).
+* [`unsafe_simplex`](@ref).
+* [`unsafe_cofacet`](@ref).
 """
-    distances(metric, points[, births])
+abstract type AbstractRipsFiltration{I<:Signed, T} <: AbstractFiltration end
 
-Return distance matrix calculated from `points` with `metric`. If given, add birth times
-from `births` to the diagonal.
 """
-function distances(metric, points)
-    isempty(points) && throw(ArgumentError("`points` must be nonempty"))
-
-    dim = length(first(points))
-    T = eltype(first(points))
-    dists = pairwise(metric, reshape(reinterpret(T, points), (dim, length(points))), dims=2)
-    return dists
-end
+    dist(::AbstractRipsFiltration, u, v)
+    dist(::AbstractRipsFiltration)
 
-# rips filtrations ======================================================================= #
+Return the distance between vertices `u` and `v`. If the distance is somehow invalid, it may
+return `missing` instead. If `u` and `v` are not given, return the distance matrix.
 """
-    AbstractRipsFiltration{I<:Signed, T} <: AbstractFiltration
-
-An abstract Vietoris-Rips filtration. Its subtypes can overload
-[`dist`](@ref)`(::AbstractRipsFiltration{T}, u, v)::Union{T, Missing}` instead of
-[`diam`](@ref).
+dist(::AbstractRipsFiltration, ::Any, ::Any)
 
-[`diam`](@ref)`(::AbstractRipsFiltration, ...)` then defaults to maximum [`dist`] among
-vertices.
+simplex_type(::AbstractRipsFiltration{I, T}, dim) where {I, T} = Simplex{dim, T, I}
 
-Comes with default implementations of [`edges`](@ref) and [`simplex_type`](@ref).
-"""
-abstract type AbstractRipsFiltration{I<:Signed, T} <: AbstractFiltration end
+function unsafe_simplex(
+    rips::AbstractRipsFiltration{I}, ::Val{0}, vertex::NTuple{1}, sign=1
+) where I
+    v = only(vertex)
+    return simplex_type(rips, 0)(I(sign) * v, birth(rips, v))
+end
 
-@propagate_inbounds function diam(rips::AbstractRipsFiltration{<:Any, T}, vertices) where T
+@inline @propagate_inbounds function unsafe_simplex(
+    rips::AbstractRipsFiltration{I, T}, ::Val{D}, vertices, sign=1
+) where {I, T, D}
     n = length(vertices)
-    res = typemin(T)
+    diameter = typemin(T)
     for i in 1:n, j in i+1:n
         d = dist(rips, vertices[j], vertices[i])
-        ismissing(d) && return missing
-        res = ifelse(d < res, res, d)
+        if ismissing(d) || d > threshold(rips)
+            return nothing
+        else
+            _d::T = d
+            diameter = ifelse(_d > diameter, _d, diameter)
+        end
     end
-    return ifelse(res > threshold(rips), missing, res)
+    return simplex_type(rips, D)(I(sign) * index(vertices), diameter)
 end
 
-@propagate_inbounds function diam(rips::AbstractRipsFiltration, sx::AbstractSimplex, us, v)
-    res = diam(sx)
-    for u in us
+@inline @propagate_inbounds function unsafe_cofacet(
+    rips::AbstractRipsFiltration{I, T},
+    simplex::IndexedSimplex{D},
+    cofacet_vertices,
+    new_vertex,
+    sign=1,
+) where {I, T, D}
+    diameter = diam(simplex)
+    for v in cofacet_vertices
+        v == new_vertex && continue
         # Even though this looks like a tight loop, v changes way more often than us, so
-        # this is the faster order of indexing by u and v.
-        d = dist(rips, v, u)
+        # this is the faster order of indexing by new_vertex and v.
+        d = dist(rips, new_vertex, v)
         if ismissing(d) || d > threshold(rips)
-            return missing
+            return nothing
         else
-            res = ifelse(res > d, res, d)
+            _d::T = d
+            diameter = ifelse(_d > diameter, _d, diameter)
         end
     end
-    return res
+    return simplex_type(rips, D + 1)(I(sign) * index(cofacet_vertices), diameter)
+end
+
+@inline @propagate_inbounds function unsafe_cofacet(
+    rips::AbstractRipsFiltration{I},
+    sx::IndexedSimplex{D},
+    cofacet_vertices,
+    ::Any,
+    new_edges::SVector,
+    sign=1,
+) where {I, D}
+    new_diam = diam(sx)
+    for i in 1:D + 1
+        e = new_edges[i]
+        e > threshold(rips) && return nothing
+        new_diam = ifelse(e > new_diam, e, new_diam)
+    end
+    new_index = index(cofacet_vertices)
+    return simplex_type(rips, D + 1)(I(sign) * new_index, new_diam)
+end
+
+function edges(rips::AbstractRipsFiltration)
+    if issparse(dist(rips))
+        _sparse_edges(rips)
+    else
+        _full_edges(rips)
+    end
 end
 
-edges(rips::AbstractRipsFiltration) = edges(dist(rips), threshold(rips), edge_type(rips))
-simplex_type(rips::AbstractRipsFiltration{I, T}, dim) where {I, T} = Simplex{dim, T, I}
+function _full_edges(rips::AbstractRipsFiltration)
+    result = edge_type(rips)[]
+    @inbounds for u in 1:size(dist(rips), 1), v in u+1:size(dist(rips), 1)
+        sx = unsafe_simplex(rips, Val(1), (v, u), 1)
+        !isnothing(sx) && push!(result, sx)
+    end
+    return result
+end
+
+function _sparse_edges(rips::AbstractRipsFiltration)
+    result = edge_type(rips)[]
+    rows = rowvals(rips.dist)
+    vals = nonzeros(rips.dist)
+    for u in 1:size(rips.dist, 1)
+        for i in nzrange(dist(rips), u)
+            v = rows[i]
+            if v > u
+                sx = unsafe_simplex(rips, Val(1), (v, u), 1)
+                !isnothing(sx) && push!(result, sx)
+            end
+        end
+    end
+    return result
+end
 
 """
-    dist(::AbstractRipsFiltration, u, v)
-    dist(::AbstractRipsFiltration)
+    distances(metric, points[, births])
 
-Return the distance between vertices `u` and `v`. If the distance is higher than the
-threshold, return `missing` instead. If `u` and `v` are not given, return the distance
-matrix.
+Return distance matrix calculated from `points` with `metric`.
 """
-dist(::AbstractRipsFiltration, ::Any, ::Any)
+function distances(metric, points)
+    isempty(points) && throw(ArgumentError("`points` must be nonempty"))
+
+    dim = length(first(points))
+    T = eltype(first(points))
+    dists = pairwise(metric, reshape(reinterpret(T, points), (dim, length(points))), dims=2)
+    return dists
+end
 
-# rips filtration ======================================================================== #
 """
     default_rips_threshold(dists)
 
-The default threshold is equal to the radius of the input space. At this threshold, all
-vertices are connected to a vertex `x` and the homology becomes trivial. If any of the
-distances is negative, default threshold defaults to `typemax(eltype(dists))`.
+The default threshold is equal to the radius of the input space. At this threshold, there
+exists a vertex ``v`` such that all vertices are connected to it and the homology becomes
+trivial.
 """
 function default_rips_threshold(dists::AbstractMatrix{T}) where T
     return minimum(maximum(abs, dists[:, i]) for i in 1:size(dists, 1))
@@ -116,10 +152,12 @@ end
 This type represents a filtration of Vietoris-Rips complexes.
 Diagonal items are treated as vertex birth times.
 
+Threshold defaults to radius of input space.
+
 # Constructors
 
-* `Rips(distance_matrix; threshold=default_rips_threshold(dist))`
-* `Rips(points; metric=Euclidean(), threshold=default_rips_threshold(dist))`
+* `Rips(distance_matrix; threshold=nothing)`
+* `Rips(points; metric=Euclidean(), threshold=nothing)`
 * `Rips{I}(args...)`: `I` sets the size of integer used to represent simplices.
 """
 struct Rips{I, T, A<:AbstractMatrix{T}} <: AbstractRipsFiltration{I, T}
@@ -131,7 +169,7 @@ function Rips{I}(
     dist::AbstractMatrix{T}; threshold=nothing
 ) where {I, T}
     issymmetric(dist) ||
-        throw(ArgumentError("`dist` must be a distance matrix"))
+        throw(ArgumentError("`dist` must be symmetric"))
     !issparse(dist) ||
         throw(ArgumentError("`dist` is sparse. Use `SparseRips` instead"))
     thresh = isnothing(threshold) ? default_rips_threshold(dist) : T(threshold)
@@ -145,7 +183,7 @@ function Rips(dist; kwargs...)
 end
 
 @propagate_inbounds function dist(rips::Rips{<:Any, T}, i, j) where T
-    return ifelse(i == j, zero(T), rips.dist[i, j])
+    return rips.dist[i, j]
 end
 dist(rips::Rips) = rips.dist
 
@@ -157,7 +195,7 @@ birth(rips::Rips, i) = rips.dist[i, i]
 """
     SparseRips{I, T} <: AbstractRipsFiltration{T, Simplex}
 
-This type represents a filtration of Vietoris-Rips complexes.
+This type represents a sparse filtration of Vietoris-Rips complexes.
 The distance matrix will be converted to a sparse matrix with all values greater than
 threshold deleted. Off-diagonal zeros in the matrix are treated as `missing`. Diagonal items
 are treated as vertex birth times.
@@ -168,38 +206,120 @@ are treated as vertex birth times.
 * `SparseRips(distance_matrix; threshold=nothing)`: `I` sets the integer size used to
   represent simplices.
 """
-struct SparseRips{I, T, A<:AbstractSparseMatrix{T}} <: AbstractRipsFiltration{I, T}
-    dist::A
+struct SparseRips{I, T} <: AbstractRipsFiltration{I, T}
+    dist::SparseMatrixCSC{T, Int}
     threshold::T
 end
 
 function SparseRips{I}(
-    dist::AbstractSparseMatrix{T}; threshold=nothing
+    dist::AbstractMatrix{T}; threshold=nothing
 ) where {I, T}
-    issymmetric(dist) || throw(ArgumentError("`dist` must be a distance matrix"))
+    issymmetric(dist) || throw(ArgumentError("`dist` must be symmetric"))
     if isnothing(threshold)
-        thresh = maximum(dist)
-        dists = sparse(dist)
-    else
-        thresh = threshold
-        dists = SparseArrays.fkeep!(copy(dist), (_, _, v) -> v ≤ threshold)
+        threshold = issparse(dist) ? maximum(dist) : default_rips_threshold(dist)
     end
-    return SparseRips{I, T, typeof(dists)}(dists, thresh)
+    dists = SparseArrays.fkeep!(SparseMatrixCSC(dist), (_, _, v) -> v ≤ threshold)
+
+    return SparseRips{I, T}(dists, threshold)
 end
 function SparseRips(dist; threshold=nothing)
     return SparseRips{Int}(dist; threshold=threshold)
 end
-function SparseRips{I}(dist; threshold=nothing) where I
-    return SparseRips{I}(sparse(dist); threshold=threshold)
-end
 
 @propagate_inbounds function dist(rips::SparseRips{<:Any, T}, i, j) where T
     res = rips.dist[i, j]
-    return ifelse(i == j, zero(T), ifelse(iszero(res), missing, res))
+    return ifelse(i == j, res, ifelse(iszero(res), missing, res))
 end
 dist(rips::SparseRips) = rips.dist
 
 n_vertices(rips::SparseRips) = size(rips.dist, 1)
 threshold(rips::SparseRips) = rips.threshold
 birth(rips::SparseRips, i) = rips.dist[i, i]
-simplex_type(rips::SparseRips{I, T}, dim) where {I, T} = Simplex{dim, T, I}
+simplex_type(::SparseRips{I, T}, dim) where {I, T} = Simplex{dim, T, I}
+
+# This is the coboundary used when distance matrix in AbstractRipsFiltration is sparse.
+function coboundary(
+    rips::AbstractRipsFiltration, sx::AbstractSimplex, ::Val{A}=Val(true)
+) where A
+    if dist(rips) isa SparseMatrixCSC
+        return SparseCoboundary{A}(rips, sx)
+    else
+        return Coboundary{A}(rips, sx)
+    end
+end
+
+struct SparseCoboundary{A, D, I, F, S}
+    filtration::F
+    simplex::S
+    vertices::SVector{D, I}
+    ptrs_begin::SVector{D, I}
+    ptrs_end::SVector{D, I}
+
+    function SparseCoboundary{A}(
+        filtration::F, simplex::AbstractSimplex{D, <:Any, I}
+    ) where {A, D, I, F}
+        verts = vertices(simplex)
+        colptr = dist(filtration).colptr
+        ptrs_begin = colptr[verts .+ 1]
+        ptrs_end = colptr[verts]
+        return new{A, D + 1, I, F, typeof(simplex)}(
+            filtration, simplex, verts, ptrs_begin, ptrs_end
+        )
+    end
+end
+
+@propagate_inbounds @inline function next_common(
+    ptrs::SVector{D}, ptrs_end::SVector{D}, rowval
+) where D
+    # could also indicate when m is equal to one of the points
+    ptrs = ptrs .- 1
+    for i in 1:D
+        ptrs[i] < ptrs_end[i] && return zero(ptrs), 0
+    end
+    m = rowval[ptrs[2]]
+    i = 1
+    while true
+        ptrs_i = ptrs[i]
+        row = rowval[ptrs_i]
+        while row > m
+            ptrs -= SVector(ntuple(isequal(i), Val(D)))
+            ptrs_i -= 1
+            ptrs_i < ptrs_end[i] && return zero(ptrs), zero(eltype(rowval))
+            row = rowval[ptrs_i]
+        end
+        i = ifelse(row == m, i + 1, 1)
+        i > D && return ptrs, m
+        m = row
+    end
+end
+
+function Base.iterate(
+    it::SparseCoboundary{A, D, I}, (ptrs, k)=(it.ptrs_begin, D)
+) where {A, D, I}
+    rowval = dist(it.filtration).rowval
+    nzval = dist(it.filtration).nzval
+    @inbounds while true
+        ptrs, v = next_common(ptrs, it.ptrs_end, rowval)
+        if iszero(first(ptrs))
+            return nothing
+        elseif k > 0 && v == it.vertices[end + 1 - k]
+            k -= 1
+        else
+            while k > 0 && v < it.vertices[end + 1 - k]
+                k -= 1
+            end
+            !A && k ≠ D && return nothing
+
+            sign = ifelse(iseven(k), 1, -1)
+            new_vertices = insert(it.vertices, D - k + 1, v)
+            new_edges = nzval[ptrs]
+            sx = unsafe_cofacet(
+                it.filtration, it.simplex, new_vertices, v, new_edges, sign
+            )
+            if !isnothing(sx)
+                _sx::simplex_type(it.filtration, D) = sx
+                return _sx, (ptrs, k)
+            end
+        end
+    end
+end
diff --git a/src/simplex.jl b/src/simplex.jl
index cd1d09a7..32290c55 100644
--- a/src/simplex.jl
+++ b/src/simplex.jl
@@ -14,11 +14,15 @@ By defining the [`index`](@ref), a default implementation of `sign`, `isless`,
 * `IndexedSimplex{D[, T, I]}(index::I, diam::T)` - constructor.
 * `IndexedSimplex{D[, T, I]}(vertices::NTuple{D+1, I}, diam::T)` - constructor.
 * [`diam(::AbstractSimplex)`](@ref)
-* [`coface_type(::AbstractSimplex)`](@ref)
 * [`index(::IndexedSimplex)`](@ref)
 """
 abstract type IndexedSimplex{D, T, I<:Integer} <: AbstractSimplex{D, T, I} end
 
+function Base.show(io::IO, ::MIME"text/plain", sx::IndexedSimplex{D}) where {D}
+    print(io, D, "-dim ", nameof(typeof(sx)), (index(sx), diam(sx)))
+    print(io, ":\n  $(sign(sx) == 1 ? '+' : '-')$(vertices(sx))")
+end
+
 """
     index(simplex::IndexedSimplex)
 
@@ -39,6 +43,8 @@ Base.sign(sx::IndexedSimplex) = sign(index(sx))
 Base.abs(sx::S) where S<:IndexedSimplex = S(abs(index(sx)), diam(sx))
 Base.:-(sx::S) where S<:IndexedSimplex = S(-index(sx), diam(sx))
 
+Base.length(::Type{<:IndexedSimplex{D}}) where D = D + 1
+
 function Base.isless(sx1::S, sx2::S) where S<:IndexedSimplex
     return ifelse(
         diam(sx1) ≠ diam(sx2),
@@ -57,12 +63,6 @@ function Base.hash(sx::IndexedSimplex, h::UInt64)
     return hash(index(sx), hash(diam(sx), h))
 end
 
-Base.eltype(sx::IndexedSimplex{<:Any, <:Any, I}) where I = I
-Base.getindex(sx::IndexedSimplex, i) = vertices(sx)[i]
-Base.firstindex(sx::IndexedSimplex) = 1
-Base.lastindex(sx::IndexedSimplex{D}) where D = D + 1
-Base.size(sx::IndexedSimplex{D}) where D = (D + 1,)
-
 # vertices and indices =================================================================== #
 """
     small_binomial(n, ::Val{k})
@@ -179,76 +179,6 @@ where ``i_k`` are the simplex vertex indices.
 end
 
 index(vertex::Integer) = vertex
-index(vertex::CartesianIndex) = index(vertex.I)
-
-# (co)boundaries ========================================================================= #
-struct IndexedCobounary{all_cofaces, D, I, F, S<:IndexedSimplex}
-    filtration::F
-    simplex::S
-    vertices::NTuple{D, I}
-
-    function IndexedCobounary{A}(
-        filtration::F, simplex::S
-    ) where {A, D, I, F, S<:IndexedSimplex{D, <:Any, I}}
-        return new{A, D + 1, I, F, S}(filtration, simplex, Tuple(vertices(simplex)))
-    end
-end
-
-function coboundary(filtration, simplex::IndexedSimplex, ::Val{A}=Val(true)) where A
-    return IndexedCobounary{A}(filtration, simplex)
-end
-
-function Base.iterate(
-    ci::IndexedCobounary{all_cofaces, D, I}, (v, k)=(I(n_vertices(ci.filtration) + 1), D),
-) where {all_cofaces, D, I}
-    @inbounds while true
-        v -= one(I)
-        while k > 0 && v == ci.vertices[end + 1 - k]
-            all_cofaces || return nothing
-            v -= one(I)
-            k -= 1
-        end
-        v > 0 || return nothing
-        diameter = diam(ci.filtration, ci.simplex, ci.vertices, v)
-        if !ismissing(diameter)
-            sign = ifelse(iseven(k), one(I), -one(I))
-            new_index = index(TupleTools.insertafter(ci.vertices, D - k, (v,))) * sign
-
-            return coface_type(ci.simplex)(new_index, diameter), (v, k)
-        end
-    end
-end
-
-struct IndexedBoundary{D, I, F, S<:IndexedSimplex}
-    filtration::F
-    simplex::S
-    vertices::NTuple{D, I}
-
-    function IndexedBoundary(
-        filtration::F, simplex::S
-    ) where {D, I, F, S<:IndexedSimplex{D, <:Any, I}}
-        return new{D + 1, I, F, S}(filtration, simplex, Tuple(vertices(simplex)))
-    end
-end
-
-function boundary(filtration, simplex::IndexedSimplex)
-    return IndexedBoundary(filtration, simplex)
-end
-
-function Base.iterate(bi::IndexedBoundary{D, I}, k=1) where {D, I}
-    while k ≤ D
-        face_vertices = TupleTools.deleteat(bi.vertices, k)
-        diameter = diam(bi.filtration, face_vertices)
-        k += 1
-        if !ismissing(diameter)
-            sign = ifelse(iseven(k), one(I), -one(I))
-            new_index = index(face_vertices) * sign
-
-            return face_type(bi.simplex)(new_index, diameter), k
-        end
-    end
-    return nothing
-end
 
 """
     Simplex{D, T, I} <: IndexedSimplex{D, T, I}
@@ -302,17 +232,6 @@ function Simplex{D, T, I}(vertices, diam) where {D, T, I<:Integer}
     return Simplex{D, T, I}(index(vertices_svec), T(diam))
 end
 
-function Base.show(io::IO, ::MIME"text/plain", sx::Simplex{D, T, I}) where {D, T, I}
-    print(io, D, "-dim Simplex", (index(sx), diam(sx)))
-    print(io, ":\n  $(sign(sx) == 1 ? '+' : '-')$(vertices(sx))")
-end
-
-function Base.show(io::IO, sx::Simplex{D, M}) where {D, M}
-    print(io, "Simplex{$D}($(sign(sx) == 1 ? '+' : '-')$(vertices(sx)), $(diam(sx)))")
-end
-
 # Interface implementation =============================================================== #
 index(sx::Simplex) = sx.index
 diam(sx::Simplex) = sx.diam
-coface_type(::Type{<:Simplex{D, T, I}}) where {D, T, I} = Simplex{D+1, T, I}
-face_type(::Type{<:Simplex{D, T, I}}) where {D, T, I} = Simplex{D-1, T, I}
diff --git a/src/zerodimensional.jl b/src/zerodimensional.jl
index 4ba05d62..9ca39b9e 100644
--- a/src/zerodimensional.jl
+++ b/src/zerodimensional.jl
@@ -71,14 +71,12 @@ end
 
 function interval(
     ::Type{R}, dset::DisjointSetsWithBirth, filtration, vertex, edge, cutoff
-) where {V, R<:RepresentativeInterval{PersistenceInterval, V}}
+) where R<:RepresentativeInterval
     birth, death = birth_death(dset, vertex, edge)
     if death - birth > cutoff
-        rep = map(find_leaves!(dset, vertex)) do u
-            V((u,), Ripserer.birth(filtration, u))
-        end
-        birth_vertex = findfirst(r -> diam(r) == birth, rep)
-        return R(PersistenceInterval(birth, death), V((birth_vertex,), birth), edge, rep)
+        rep = simplex.(Ref(filtration), Val(0), tuple.(find_leaves!(dset, vertex)))
+        birth_simplex = simplex(filtration, Val(0), (findfirst(r -> diam(r) == birth, rep),))
+        return R(PersistenceInterval(birth, death), birth_simplex, edge, rep)
     else
         return nothing
     end
@@ -122,7 +120,7 @@ function zeroth_intervals(
     end
     to_skip = edge_type(filtration)[]
     to_reduce = edge_type(filtration)[]
-    simplices = edges(filtration)
+    simplices = sort!(edges(filtration))
     if progress
         progbar = Progress(length(simplices), desc="Computing 0d intervals... ")
     end
diff --git a/test/cubical.jl b/test/cubical.jl
index 9787abae..2e2171ab 100644
--- a/test/cubical.jl
+++ b/test/cubical.jl
@@ -3,8 +3,7 @@ using Ripserer
 using Compat
 using StaticArrays
 
-using Ripserer: all_equal_in_dim, boundary, coboundary, face_type, coface_type,
-    edges, n_vertices, index
+using Ripserer: all_equal_in_dim, boundary, coboundary, edges, n_vertices, index
 
 data1d = cos.(range(0, 4π, length=1000))
 
@@ -23,8 +22,8 @@ using ..TestHelpers: test_indexed_simplex_interface, test_filtration_interface
 @testset "Cubelet" begin
     test_indexed_simplex_interface(Cubelet, D -> 2^D)
 
-    @testset "Vertices, index." begin
-        @testset "Basics." begin
+    @testset "Vertices, index" begin
+        @testset "Basics" begin
             @test vertices(Cubelet{2}(1, rand())) === SVector{4}(4, 3, 2, 1)
             @test vertices(Cubelet{2}(Int128(2), rand())) === SVector{4, Int128}(5, 3, 2, 1)
             @test vertices(Cubelet{1}(3, rand())) === SVector{2}(3, 2)
@@ -32,7 +31,7 @@ using ..TestHelpers: test_indexed_simplex_interface, test_filtration_interface
             @test vertices(Cubelet{3}(5, rand())) === SVector{8}(9, 8, 7, 6, 4, 3, 2, 1)
         end
 
-        @testset "Index to vertices and back." begin
+        @testset "Index to vertices and back" begin
             for i in 1:20
                 cube = Cubelet{5}(i, rand())
                 @test index(vertices(cube)) == i
@@ -42,7 +41,7 @@ using ..TestHelpers: test_indexed_simplex_interface, test_filtration_interface
             end
         end
     end
-    @testset "Printing." begin
+    @testset "Printing" begin
         @test sprint(print, Cubelet{1}(1, 1)) == "Cubelet{1}(+[2, 1], 1)"
         @test sprint(print, Cubelet{2}(-1, 1)) == "Cubelet{2}(-[4, 3, 2, 1], 1)"
 
@@ -59,37 +58,34 @@ end
 @testset "Cubical" begin
     test_filtration_interface(Cubical, (data1d, data2d, data3d))
 
-    @testset "n_vertices, indices, birth, diam" begin
+    @testset "n_vertices, birth" begin
         for data in (data1d, data2d, data3d)
             filtration = Cubical(data)
 
             @test n_vertices(filtration) == length(data)
-            @test CartesianIndices(filtration) == CartesianIndices(data)
-            @test LinearIndices(filtration) == LinearIndices(data)
             @test birth(filtration, 10) == data[10]
-            @test diam(filtration, (10, 9)) == max(data[10], data[9])
         end
     end
 end
 
-@testset "Coboundary." begin
+@testset "Coboundary" begin
     @testset "all_equal_in_dim" begin
         @test all_equal_in_dim(1, [(1, 1), (1, 2), (1, 3), (1, 4)])
         @test !all_equal_in_dim(2, [(1, 1), (1, 1), (1, 2), (1, 2)])
     end
 
-    @testset "Edges have no coboundary in 1d." begin
+    @testset "Edges have no coboundary in 1d" begin
         cob = Cubelet{2, Float64, Int}[]
         flt = Cubical(data1d)
-        cub = Cubelet{1}((3, 2), diam(flt, (3, 2)))
+        cub = simplex(flt, Val(1), (3, 2))
         for c in coboundary(flt, cub)
             push!(cob, c)
         end
         @test isempty(cob)
     end
 
-    @testset "Coboundary of edges in 2d." begin
-        cob = Cubelet{2, Int, Int}[]
+    @testset "Coboundary of edges in 2d" begin
+        cob = []
         flt = Cubical(data2d)
         cub = Cubelet{1}((10, 9), 1)
         for c in coboundary(flt, cub)
@@ -101,7 +97,7 @@ end
         @test cob[2] == -Cubelet{2}((10, 9, 3, 2), 2)
         @test sign(cob[2]) == -1
 
-        cocob = coface_type(eltype(cob))[]
+        cocob = []
         for c in coboundary(flt, cob[1])
             push!(cocob, c)
         end
@@ -111,14 +107,14 @@ end
         @test isempty(cocob)
     end
 
-    @testset "Dimaters of 2-cubelet coboundary in 2×2×2 3d data." begin
+    @testset "Diameters of 2-cubelet coboundary in 2×2×2 3d data" begin
         data = zeros(2, 2, 2)
         data[:, 1, 1] .= 1
         data[:, 2, 1] .= 2
         data[:, 1, 2] .= 3
         data[:, 2, 2] .= 4
 
-        @testset "Diameter increasing, cofaces are supersets of original." begin
+        @testset "Diameter increasing, cofacets are supersets of original" begin
             cob = Cubelet{2, Float64, Int}[]
             for cube in edges(Cubical(data))
                 for c in coboundary(Cubical(data), cube)
@@ -127,7 +123,7 @@ end
                 end
             end
         end
-        @testset "Coboundary with all_cofaces=false." begin
+        @testset "Coboundary with all_cofacets=false" begin
             cob = Cubelet{2, Float64, Int}[]
             for cube in edges(Cubical(data))
                 for c in coboundary(Cubical(data), cube, Val(false))
@@ -137,7 +133,7 @@ end
             @test length(cob) == 6
             @test allunique(cob)
         end
-        @testset "All 2-cubelets have the same coboundary." begin
+        @testset "All 2-cubelets have the same coboundary" begin
             for cube in [
                 Cubelet{2}((1, 2, 3, 4), 2.0),
                 Cubelet{2}((1, 2, 5, 6), 4.0),
@@ -146,13 +142,13 @@ end
                 Cubelet{2}((3, 4, 7, 8), 4.0),
                 Cubelet{2}((5, 6, 7, 8), 4.0),
             ]
-                coface = only(coboundary(Cubical(data), cube))
-                @test abs(coface) == Cubelet{3}((8, 7, 6, 5, 4, 3, 2, 1), 4.0)
+                cofacet = only(coboundary(Cubical(data), cube))
+                @test abs(cofacet) == Cubelet{3}((8, 7, 6, 5, 4, 3, 2, 1), 4.0)
             end
         end
     end
-    @testset "Counting cofaces in 3d." begin
-        cob = Cubelet{2, Int, Int}[]
+    @testset "Counting cofacets in 3d" begin
+        cob = []
         flt = Cubical(data3d)
         cub = Cubelet{1}((14, 13), 1)
         for c in coboundary(flt, cub)
@@ -160,7 +156,7 @@ end
         end
         @test length(cob) == 3
 
-        cocob = coface_type(eltype(cob))[]
+        cocob = []
         for c in coboundary(flt, cob[1])
             push!(cocob, c)
         end
@@ -168,8 +164,8 @@ end
     end
 end
 
-@testset "Boundary." begin
-    @testset "Boundary of an edge in 1d." begin
+@testset "Boundary" begin
+    @testset "Boundary of an edge in 1d" begin
         bnd = Cubelet{0, Float64, Int}[]
         flt = Cubical(data3d)
         cub = Cubelet{1}((3, 2), 1.0)
@@ -178,7 +174,7 @@ end
         end
         @test bnd == [Cubelet{0}((3,), 1.0), -Cubelet{0}((2,), 1.0)]
     end
-    @testset "Boundary of 2-cubelet in 2d." begin
+    @testset "Boundary of 2-cubelet in 2d" begin
         bnd = Cubelet{1, Int, Int}[]
         flt = Cubical(data2d)
         cub = Cubelet{2}((10, 9, 3, 2), 2)
@@ -188,7 +184,7 @@ end
         @test bnd == [Cubelet{1}((10, 3), 2), -Cubelet{1}((9, 2), 2),
                       Cubelet{1}((10, 9), 1), -Cubelet{1}((3, 2), 2)]
     end
-    @testset "Boundary of 3-cubelet in 3d." begin
+    @testset "Boundary of 3-cubelet in 3d" begin
         bnd = Cubelet{2, Float64, Int}[]
         flt = Cubical(data3d)
         cub = Cubelet{3}((112, 111, 102, 101, 12, 11, 2, 1), 1.0)
diff --git a/test/reductionmatrix.jl b/test/reductionmatrix.jl
index 2be58a58..114b702f 100644
--- a/test/reductionmatrix.jl
+++ b/test/reductionmatrix.jl
@@ -2,15 +2,24 @@ using Ripserer
 
 using Random
 
-using Ripserer: chain_element_type, coefficient, coface_type, face_type, index
+using Ripserer: chain_element_type, coefficient, index
 
 using Ripserer: ReducedMatrix, record!, commit!, discard!
 using Ripserer: WorkingBoundary, nonheap_push!, get_pivot!, repair!
-using Ripserer: ReductionMatrix, simplex_type, simplex_element, face_element, next_matrix
+using Ripserer: ReductionMatrix, simplex_type, simplex_element, facet_element, next_matrix
+
+cofacet_type(::Type{<:A}) where {D, T, I, A<:Cubelet{D, T, I}} =
+    Cubelet{D + 1, T, I}
+cofacet_type(::Type{<:A}) where {D, T, I, A<:Simplex{D, T, I}} =
+    Simplex{D + 1, T, I}
+facet_type(::Type{<:A}) where {D, T, I, A<:Cubelet{D, T, I}} =
+    Cubelet{D - 1, T, I}
+facet_type(::Type{<:A}) where {D, T, I, A<:Simplex{D, T, I}} =
+    Simplex{D - 1, T, I}
 
 @testset "ReducedMatrix" begin
     for S in (Simplex{2, Int, Int}, Cubelet{1, Int, Int}), T in (Mod{3}, Rational{Int})
-        C = coface_type(S)
+        C = cofacet_type(S)
         CE = chain_element_type(C, T)
         SE = chain_element_type(S, T)
 
@@ -21,7 +30,7 @@ using Ripserer: ReductionMatrix, simplex_type, simplex_element, face_element, ne
         rev = Base.Order.Reverse
 
         @testset "ReducedMatrix with simplex type $S and field type $T" begin
-            @testset "a fresh ReducedMatrix is empty even if you commit nothing." begin
+            @testset "a fresh ReducedMatrix is empty even if you commit nothing" begin
                 matrix = ReducedMatrix{C, SE}(fwd)
                 commit!(matrix, columns[1], T(2))
 
@@ -231,7 +240,7 @@ end
     end
 end
 
-@testset "ReducedMatrix" begin
+@testset "ReductionMatrix" begin
     for Co in (true, false),
         S in (Simplex{2, Int, Int}, Cubelet{1, Int, Int}),
         T in (Mod{3}, Rational{Int})
@@ -243,7 +252,7 @@ end
         co_str = Co ? "Coh" : "H"
 
         SE = chain_element_type(S, T)
-        F = Co ? coface_type(S) : face_type(S)
+        F = Co ? cofacet_type(S) : facet_type(S)
         FE = chain_element_type(F, T)
 
         @testset "$(co_str)omology, $S, $T" begin
@@ -257,7 +266,7 @@ end
 
                 @test simplex_type(matrix) == S
                 @test simplex_element(matrix) == SE
-                @test face_element(matrix) == FE
+                @test facet_element(matrix) == FE
                 @test dim(matrix) == (Co ? dim(S) : dim(S) - 1)
             end
 
diff --git a/test/ripserer.jl b/test/ripserer.jl
index d69d70fb..3cc7927b 100644
--- a/test/ripserer.jl
+++ b/test/ripserer.jl
@@ -100,6 +100,24 @@ end
         @test d1_3 == d1_331 == d1_r == []
         @test d2_3 == d2_331 == d1_r == []
     end
+    @testset "Equal to Rips" begin
+        for thresh in (nothing, 1, 0.5, 0.126)
+            data = rand_torus(100)
+            r_res = ripserer(data, threshold=thresh, dim_max=2)
+            s_res_1 = ripserer(SparseRips(data, threshold=thresh), dim_max=2)
+
+            # Add zeros to diagonal. Adding ones first actually changes the structure of
+            # the matrix.
+            data2 = sparse(data)
+            for i in axes(data2, 1)
+                data2[i, i] = 1
+                data2[i, i] = 0
+            end
+            s_res_2 = ripserer(data2, threshold=thresh, dim_max=2)
+
+            @test r_res == s_res_1 == s_res_2
+        end
+    end
 end
 
 @testset "Representatives" begin
@@ -301,8 +319,12 @@ Ripserer.n_vertices(::CustomRips) = 10
 end
 
 struct CustomFiltration <: Ripserer.AbstractFiltration end
-Ripserer.diam(::CustomFiltration, ::Simplex, _, _) = 1
-Ripserer.diam(::CustomFiltration, _) = 1
+function Ripserer.unsafe_simplex(::CustomFiltration, ::Val{0}, (v,), sign=1)
+    return Simplex{0}(sign * v, 0)
+end
+function Ripserer.unsafe_simplex(::CustomFiltration, ::Val{D}, vertices, sign=1) where D
+    return Simplex{D}(sign * Ripserer.index(vertices), 1)
+end
 Ripserer.n_vertices(::CustomFiltration) = 10
 Ripserer.simplex_type(::CustomFiltration, D) = Simplex{D, Int, Int}
 Ripserer.edges(::CustomFiltration) = Simplex{1}.(10:-1:1, 1)
diff --git a/test/ripsfiltration.jl b/test/ripsfiltration.jl
index 9081b964..de914687 100644
--- a/test/ripsfiltration.jl
+++ b/test/ripsfiltration.jl
@@ -1,5 +1,5 @@
 using Ripserer
-using Ripserer: distances, vertex_type, edge_type, dist, edges, n_vertices
+using Ripserer: distances, vertex_type, edge_type, dist, edges, n_vertices, unsafe_simplex
 
 using Compat
 using StaticArrays
@@ -7,59 +7,6 @@ using StaticArrays
 using ..TestHelpers: test_filtration_interface
 include("data.jl")
 
-@testset "Edges for matrix inputs..." begin
-    @testset "in a dense matrix." begin
-        dist = [0 1 2 3;
-                1 0 4 5;
-                2 4 0 4;
-                3 5 4 0]
-        res = edges(dist, 4, Simplex{1, Int, Int})
-
-        @test eltype(res) === Simplex{1, Int, Int}
-        @test length(res) == 5
-        @test res[1] == Simplex{1}(1, 1)
-        @test res[2] == Simplex{1}(2, 2)
-        @test res[3] == Simplex{1}(4, 3)
-        @test res[4] == Simplex{1}(6, 4)
-        @test res[5] == Simplex{1}(3, 4)
-    end
-    @testset "in a sparse matrix." begin
-        dist = sparse([0 1 0 3;
-                       1 0 4 5;
-                       0 4 0 4;
-                       3 5 4 0])
-        res = edges(dist, Inf, Simplex{1, Int, Int128})
-
-        @test eltype(res) === Simplex{1, Int, Int128}
-        @test length(res) == 5
-        @test res[1] === Simplex{1, Int, Int128}(1, 1)
-        @test res[2] === Simplex{1, Int, Int128}(4, 3)
-        @test res[3] === Simplex{1, Int, Int128}(6, 4)
-        @test res[4] === Simplex{1, Int, Int128}(3, 4)
-        @test res[5] === Simplex{1, Int, Int128}(5, 5)
-    end
-    @testset "number of edges in a dense matrix." begin
-        dist = rand_dist_matrix(100)
-        n_edges = binomial(size(dist, 1), 2)
-        res = edges(dist, Inf, Simplex{1, Float64, Int})
-
-        @test eltype(res) ≡ Simplex{1, Float64, Int}
-        @test length(res) == n_edges
-        @test issorted(res, by=diam)
-    end
-    @testset "number of edges in a sparse matrix." begin
-        for _ in 1:10
-            dist = rand_dist_matrix(100, 0.5)
-            n_edges = nnz(dist) ÷ 2
-            res = edges(dist, Inf, Simplex{1, Float64, Int128})
-
-            @test eltype(res) ≡ Simplex{1, Float64, Int128}
-            @test length(res) == n_edges
-            @test issorted(res, by=diam)
-        end
-    end
-end
-
 @testset "distances" begin
     for points in (
         [(0, 0), (0, 1), (1, 1), (1, 0)],
@@ -74,7 +21,7 @@ for Filtration in (Rips, SparseRips)
     @testset "$(string(Filtration))" begin
         test_filtration_interface(Filtration, (icosahedron, torus(100)))
 
-        @testset "With no threshold." begin
+        @testset "With no threshold" begin
             filtration = Filtration(Float64[0 1 2 9; 1 0 3 9; 2 3 0 4; 9 9 4 0])
 
             @test n_vertices(filtration) == 4
@@ -82,25 +29,29 @@ for Filtration in (Rips, SparseRips)
             @test dist(filtration, 1, 2) == 1.0
             @test dist(filtration, 1, 3) == 2.0
             @test dist(filtration, 3, 2) == 3.0
-            @test diam(filtration, Simplex{1}(3, 3), [3, 2], 1) == 3.0
-            @test diam(filtration, (4, 2, 1)) === (filtration isa Rips ? missing : 9.0)
 
-            @test threshold(filtration) == (filtration isa Rips ? 4.0 : 9.0)
+            @test threshold(filtration) == 4.0
             @test vertex_type(filtration) === Simplex{0, Float64, Int}
             @test edge_type(filtration) === Simplex{1, Float64, Int}
             @test birth(filtration, 1) == 0
             @test birth(filtration, 4) == 0
+
+            @test dist(filtration) == filtration.dist
+
+            @test unsafe_simplex(filtration, Val(0), (1,)) === Simplex{0}(1, 0.0)
+            @test unsafe_simplex(filtration, Val(1), (2, 1), -1) === -Simplex{1}(1, 1.0)
+            @test simplex(filtration, Val(2), (1, 3, 2)) === Simplex{2}(1, 3.0)
+            @test_throws ArgumentError simplex(filtration, Val(2), (1, 1, 2))
         end
-        @testset "With threshold and index type." begin
+        @testset "With threshold and index type" begin
             filtration = Filtration{Int128}([1 1 2; 1 2 3; 2 3 2]; threshold=2)
 
             @test n_vertices(filtration) == 3
             @test threshold(filtration) == 2
-            @test dist(filtration, 3, 3) == 0
+            @test dist(filtration, 3, 3) == 2
             @test dist(filtration, 1, 2) == 1
             @test dist(filtration, 1, 3) == 2
-            @test dist(filtration, 3, 2) ≡ (issparse(filtration.dist) ? missing : 3)
-            @test ismissing(diam(filtration, Simplex{2}(1, 1), [1, 2], 3))
+            @test dist(filtration, 3, 2) === (issparse(filtration.dist) ? missing : 3)
             @test threshold(filtration) == 2
             @test vertex_type(filtration) === Simplex{0, Int, Int128}
             @test edge_type(filtration) === Simplex{1, Int, Int128}
@@ -108,11 +59,18 @@ for Filtration in (Rips, SparseRips)
             @test birth(filtration, 1) == 1
             @test birth(filtration, 2) == 2
             @test birth(filtration, 3) == 2
+
+            @test dist(filtration) == filtration.dist
+
+            @test unsafe_simplex(filtration, Val(0), (1,)) === Simplex{0, Int, Int128}(1, 1)
+            @test unsafe_simplex(filtration, Val(2), (3, 2, 1)) === nothing
+            @test simplex(filtration, Val(1), (1, 2), -1) === -Simplex{1, Int, Int128}(1, 1)
+            @test_throws ArgumentError simplex(filtration, Val(2), (0, 1, 2))
         end
     end
 end
 
-@testset "Rips points constructor." begin
+@testset "Rips points constructor" begin
     filtration = Rips(torus_points(9))
     @test all(x -> x > 0, dist(filtration, i, j) for i in 1:10 for j in i+1:9)
     @test eltype(edges(filtration)) === Simplex{1, Float64, Int}
@@ -123,13 +81,13 @@ end
 end
 
 @testset "Errors" begin
-    @testset "Non-distance matrices throw an error." begin
+    @testset "Non-distance matrices throw an error" begin
         @test_throws ArgumentError Rips([1 2 3; 4 5 6; 7 8 9])
         @test_throws ArgumentError Rips(zeros(3, 2))
         @test_throws ArgumentError SparseRips([1 2 3; 4 5 6; 7 8 9])
         @test_throws ArgumentError SparseRips(zeros(3, 2))
     end
-    @testset "Constructing Rips filtration with sparse matrix not allowed." begin
+    @testset "Constructing Rips filtration with sparse matrix not allowed" begin
         @test_throws ArgumentError Rips(sparse([0 1 1; 1 0 1; 1 1 0]))
     end
 end
diff --git a/test/simplex.jl b/test/simplex.jl
index 727df416..1cf31c28 100644
--- a/test/simplex.jl
+++ b/test/simplex.jl
@@ -1,20 +1,28 @@
 using Ripserer
 using StaticArrays
 
-using Ripserer: small_binomial, boundary, coboundary, face_type, coface_type, index
+using Ripserer: small_binomial, boundary, coboundary, index
 using ..TestHelpers: test_indexed_simplex_interface
 
-struct FakeFiltration end
-Ripserer.diam(::FakeFiltration, args...) =
-    1
-Ripserer.n_vertices(::FakeFiltration) =
-    20
+struct FakeFiltration<:Ripserer.AbstractFiltration end
+function Ripserer.unsafe_simplex(::FakeFiltration, ::Val{D}, vertices, sign=1) where D
+    return Simplex{D, Int, Int}(sign * index(vertices), 1)
+end
+Ripserer.n_vertices(::FakeFiltration) = 20
+Ripserer.simplex_type(::FakeFiltration, D) = Simplex{D, Int, Int}
 
-struct FakeFiltrationWithThreshold end
-Ripserer.diam(::FakeFiltrationWithThreshold, _, _, v) =
-    v ≤ 10 ? 1 : missing
-Ripserer.n_vertices(::FakeFiltrationWithThreshold) =
-    20
+struct FakeFiltrationWithThreshold<:Ripserer.AbstractFiltration end
+function Ripserer.unsafe_simplex(
+    ::FakeFiltrationWithThreshold, ::Val{D}, vertices, sign=1
+) where D
+    if maximum(vertices) > 10
+        return nothing
+    else
+        return Simplex{D, Int, Int}(sign * index(vertices), 1)
+    end
+end
+Ripserer.n_vertices(::FakeFiltrationWithThreshold) = 20
+Ripserer.simplex_type(::FakeFiltrationWithThreshold, D) = Simplex{D, Int, Int}
 
 @testset "Binomials" begin
     @test all(binomial(n, k) == small_binomial(n, Val(k)) for n in 0:1000 for k in 0:7)
@@ -54,38 +62,38 @@ end
 
     end
     @testset "coboundary" begin
-        @testset "number of cofaces" begin
+        @testset "number of cofacets" begin
             for dim in 1:10
                 simplex = Simplex{dim}(10, 1)
                 simplex_vxs = vertices(simplex)
-                cob = coface_type(simplex)[]
-                for coface in coboundary(FakeFiltration(), simplex)
-                    push!(cob, coface)
+                cob = []
+                for cofacet in coboundary(FakeFiltration(), simplex)
+                    push!(cob, cofacet)
                 end
                 @test length(cob) == 20 - dim - 1
                 @test all(isequal(1), diam.(cob))
                 @test all(issubset(simplex_vxs, vertices(c)) for c in cob)
             end
         end
-        @testset "number of cofaces, all_cofaces=false" begin
+        @testset "number of cofacets, all_cofacets=false" begin
             for dim in 1:5
                 cob = Simplex{dim+1, Int, Int}[]
                 for idx in binomial(20, dim+1):-1:1
                     simplex = Simplex{dim}(idx, 1)
-                    for coface in coboundary(FakeFiltration(), simplex, Val(false))
-                        push!(cob, coface)
+                    for cofacet in coboundary(FakeFiltration(), simplex, Val(false))
+                        push!(cob, cofacet)
                     end
                 end
                 @test sort(index.(abs.(cob))) == 1:binomial(20, dim+2)
             end
         end
-        @testset "number of cofaces, thresholding" begin
+        @testset "number of cofacets, thresholding" begin
             for dim in 1:8 # at 9, simplex with index 10 becomes invalid.
                 simplex = Simplex{dim}(10, 1)
                 simplex_vxs = vertices(simplex)
-                cob = coface_type(simplex)[]
-                for coface in coboundary(FakeFiltrationWithThreshold(), simplex)
-                    push!(cob, coface)
+                cob = []
+                for cofacet in coboundary(FakeFiltrationWithThreshold(), simplex)
+                    push!(cob, cofacet)
                 end
                 @test length(cob) == 10 - dim - 1
                 @test all(isequal(1), diam.(cob))
@@ -100,11 +108,11 @@ end
                 @static if VERSION ≥ v"1.1.0"
                     @test begin @inferred Union{
                         Nothing,
-                        Tuple{coface_type(simplex), Tuple{Int, Int}},
+                        Tuple{Simplex{dim + 1, Int, Int}, Tuple{Int, Int}},
                     } iterate(cobiter); true end
                     @test begin @inferred Union{
                         Nothing,
-                        Tuple{coface_type(simplex), Tuple{Int, Int}},
+                        Tuple{Simplex{dim + 1, Int, Int}, Tuple{Int, Int}},
                     } iterate(cobiter, (13, dim+1)); true end
                 end
             end
@@ -112,13 +120,13 @@ end
     end
 
     @testset "boundary" begin
-        @testset "number of faces" begin
+        @testset "number of facets" begin
             for dim in 1:10
                 simplex = Simplex{dim}(10, 1)
                 simplex_vxs = vertices(simplex)
-                bnd = face_type(simplex)[]
-                for face in boundary(FakeFiltration(), simplex)
-                    push!(bnd, face)
+                bnd = []
+                for facet in boundary(FakeFiltration(), simplex)
+                    push!(bnd, facet)
                 end
                 @test length(bnd) == dim + 1
                 @test all(isequal(1), diam.(bnd))
@@ -133,11 +141,11 @@ end
                 @static if VERSION ≥ v"1.1.0"
                     @test begin @inferred Union{
                         Nothing,
-                        Tuple{face_type(simplex), Int},
+                        Tuple{Simplex{dim - 1, Int, Int}, Int},
                     } iterate(biter); true end
                     @test begin @inferred Union{
                         Nothing,
-                        Tuple{face_type(simplex), Int},
+                        Tuple{Simplex{dim - 1, Int, Int}, Int},
                     } iterate(biter, 1); true end
                 end
             end
diff --git a/test/testhelpers.jl b/test/testhelpers.jl
index 9a5c2f49..6ad9af37 100644
--- a/test/testhelpers.jl
+++ b/test/testhelpers.jl
@@ -3,7 +3,7 @@ using Test
 
 using Ripserer
 using Ripserer: vertex_type, edges, edge_type,
-    AbstractSimplex, IndexedSimplex, face_type, coface_type, index
+    AbstractSimplex, IndexedSimplex, index
 
 """
     test_indexed_simplex_interface(S, n_vertices)
@@ -45,11 +45,6 @@ function test_indexed_simplex_interface(S, n_vertices)
                 @test sign(+S{D}(i, d)) == sign(i)
                 @test sign(-S{D}(i, d)) == -sign(i)
 
-                @test coface_type(S{D}(i, d)) ≡ S{D + 1, T, I}
-                @test coface_type(typeof(S{D}(i, d))) ≡ S{D + 1, T, I}
-                @test face_type(S{D}(i, d)) ≡ S{D - 1, T, I}
-                @test face_type(typeof(S{D}(i, d))) ≡ S{D - 1, T, I}
-
                 @test dim(S{D}(i, d)) == D
                 @test abs(S{D}(i, d)) == S{D}(abs(i), d)
                 @test abs(-S{D}(i, d)) == abs(S{D}(i, d))
@@ -111,8 +106,6 @@ function test_filtration_interface(Filtration, datasets)
             @test diam(sx) isa T
             @test birth(flt, 1) isa Union{T, Bool}
 
-            @test diam(flt, (4, 3, 2, 1)) isa Union{T, Missing}
-
             @test begin @inferred edges(flt); true end
         end
     end