diff --git a/Project.toml b/Project.toml
index 053d72e5..90350a00 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.8.2"
+version = "0.9.0"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
@@ -24,7 +24,7 @@ DataStructures = "0.17"
 Distances = "0.8, 0.9"
 Hungarian = "0.6"
 IterTools = "1"
-PersistenceDiagrams = "0.4"
+PersistenceDiagrams = "0.5"
 ProgressMeter = "1"
 RecipesBase = "1"
 StaticArrays = "0.12"
diff --git a/docs/src/examples/image_sublevel.jl b/docs/src/examples/image_sublevel.jl
index 6c44bd28..f114c2ab 100644
--- a/docs/src/examples/image_sublevel.jl
+++ b/docs/src/examples/image_sublevel.jl
@@ -33,7 +33,7 @@ heatmap(blackhole_grayscale, aspect_ratio=1, size=(800, 800))
 #     be careful with computing representatives because collecting them for thousands of
 #     persistence intervals can take a while.
 
-result_min = ripserer(Cubical(blackhole_grayscale))
+result_min = ripserer(Cubical{128}(blackhole_grayscale))
 
 # We plot the diagram with the `persistence` argument. This plots persistence vs birth
 # instead of death vs birth.
@@ -44,7 +44,7 @@ plot(result_min, persistence=true)
 # remove the noise by supplying the `cutoff` argument to `ripserer`. This removes all
 # intervals with persistence strictly smaller than cutoff.
 
-result_min = ripserer(Cubical(blackhole_grayscale), cutoff=0.1)
+result_min = ripserer(Cubical{128}(blackhole_grayscale), cutoff=0.1)
 
 #
 
@@ -52,7 +52,7 @@ plot(result_min, persistence=true)
 
 # Now that we know we have a smaller number of intervals, we can compute the
 # representatives.
-result_min = ripserer(Cubical(blackhole_grayscale), cutoff=0.1, reps=true)
+result_min = ripserer(Cubical{128}(blackhole_grayscale), cutoff=0.1, reps=true)
 
 # We separate the finite and infinite intervals. The finite one corresponds to the local
 # minimum in the middle of the image and the infinite one corresponds to the global minimum.
@@ -113,7 +113,7 @@ scatter!(result_min[2][2], blackhole_grayscale,
 
 # Let's repeat what we just did, but with the image negated.
 
-result_max = ripserer(Cubical(-blackhole_grayscale), cutoff=0.1, reps=true)
+result_max = ripserer(Cubical{128}(-blackhole_grayscale), cutoff=0.1, reps=true)
 plot(result_max)
 
 heatmap(blackhole_grayscale, aspect_ratio=1, size=(800, 800))
diff --git a/src/Ripserer.jl b/src/Ripserer.jl
index acbb3913..e422fd8f 100644
--- a/src/Ripserer.jl
+++ b/src/Ripserer.jl
@@ -25,17 +25,19 @@ using TupleTools
 
 # reexporting basics these makes Ripserer usable without having to import another package.
 import PersistenceDiagrams: birth, threshold
-export PersistenceDiagram, PersistenceInterval
-export birth, death, persistence, representative, barcode
-
-export Mod
-export AbstractSimplex, diam, coface_type, face_type, vertices, coboundary, boundary, dim
-export IndexedSimplex, index, Simplex
-export AbstractFiltration, n_vertices, edges, threshold
-export AbstractFlagFiltration, Rips, SparseRips
-export Cubelet, Cubical
-export simplex, coefficient
-export ripserer
+export
+    PersistenceDiagram, PersistenceInterval, RepresentativeInterval,
+    birth, death, persistence, representative, birth_simplex, death_simplex, barcode
+
+export
+    Mod,
+    AbstractSimplex, diam, coface_type, face_type, vertices, coboundary, boundary, dim,
+    IndexedSimplex, index, Simplex,
+    AbstractFiltration, n_vertices, edges, threshold,
+    AbstractFlagFiltration, Rips, SparseRips,
+    Cubelet, Cubical,
+    simplex, coefficient,
+    ripserer
 
 include("primefield.jl")
 
diff --git a/src/chainelement.jl b/src/chainelement.jl
index 46edbcd2..9d602a9c 100644
--- a/src/chainelement.jl
+++ b/src/chainelement.jl
@@ -59,6 +59,10 @@ function Base.show(io::IO, ce::AbstractChainElement)
     print(io, simplex(ce), " => ", coefficient(ce))
 end
 
+function Base.convert(::Type{C}, simplex::S) where {S, C<:AbstractChainElement{S}}
+    C(simplex)
+end
+
 """
     chain_element_type(simplex, coefficient)
 
diff --git a/src/cubical.jl b/src/cubical.jl
index 7a624762..dea019ed 100644
--- a/src/cubical.jl
+++ b/src/cubical.jl
@@ -154,9 +154,7 @@ function Base.iterate(
 ) where {A, N, I, C}
     # If not all indices in a given dimension are equal, we can't create a coface by
     # expanding in that direction.
-    diameter = missing
-    new_vertices = cc.vertices
-    while ismissing(diameter)
+    while true
         while dim ≤ N && !all_equal_in_dim(dim, cc.vertices)
             dim += one(I)
         end
@@ -170,14 +168,14 @@ function Base.iterate(
 
         dim += I(dir == -1)
         dir *= -one(I)
-    end
 
-    if ismissing(diameter)
-        return nothing
-    else
-        all_vertices = sort(map(v -> I(LinearIndices(cc.filtration)[v]),
-                                vcat(cc.vertices, new_vertices)), rev=true)
-        return coface_type(C)(-dir * index(all_vertices), diameter), (dim, dir)
+        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)
+            end
+        end
     end
 end
 
diff --git a/src/reductionmatrix.jl b/src/reductionmatrix.jl
index f34b3de4..236d5a16 100644
--- a/src/reductionmatrix.jl
+++ b/src/reductionmatrix.jl
@@ -309,42 +309,50 @@ function reduce_column!(matrix::ReductionMatrix, column_simplex)
     return pivot
 end
 
-function birth_death(matrix::ReductionMatrix{true}, simplex, pivot)
-    return diam(simplex), isnothing(pivot) ? Inf : diam(pivot)
+function birth_death(matrix::ReductionMatrix{true}, column, pivot)
+    return diam(column), isnothing(pivot) ? Inf : diam(pivot)
 end
-function birth_death(matrix::ReductionMatrix{false}, simplex, pivot)
-    return isnothing(pivot) ? Inf : diam(pivot), diam(simplex)
+function birth_death(matrix::ReductionMatrix{false}, column, pivot)
+    return isnothing(pivot) ? Inf : diam(pivot), diam(column)
 end
 
 # Interval with no representative, (co)homology.
 function interval(
-    ::Type{PersistenceInterval{Nothing}}, matrix, simplex, pivot, cutoff
+    ::Type{PersistenceInterval}, matrix::ReductionMatrix, column, pivot, cutoff
 )
-    birth, death = birth_death(matrix, simplex, pivot)
+    birth, death = birth_death(matrix, column, pivot)
     return death - birth > cutoff ? PersistenceInterval(birth, death) : nothing
 end
 # With representative, cohomology.
 function interval(
-    ::Type{PersistenceInterval{V}}, matrix::ReductionMatrix{true}, simplex, pivot, cutoff
-) where V<:AbstractVector
-    birth, death = birth_death(matrix, simplex, pivot)
+    ::Type{R}, matrix::ReductionMatrix{true}, column, pivot, cutoff
+) where R<:RepresentativeInterval
+    birth, death = birth_death(matrix, column, pivot)
 
     if !isfinite(death)
-        return PersistenceInterval(birth, death, V())
+        return R(
+            PersistenceInterval(birth, death), column, nothing, eltype(matrix.reduced)[]
+        )
     elseif death - birth > cutoff
-        return PersistenceInterval(birth, death, collect(matrix.reduced[pivot]))
+        return R(
+            PersistenceInterval(birth, death),
+            column, simplex(pivot), collect(matrix.reduced[pivot])
+        )
     else
         return nothing
     end
 end
 # With representative, homology.
 function interval(
-    ::Type{PersistenceInterval{V}}, matrix::ReductionMatrix{false}, simplex, pivot, cutoff
-) where V<:AbstractVector
-    birth, death = birth_death(matrix, simplex, pivot)
+    ::Type{R}, matrix::ReductionMatrix{false}, column, pivot, cutoff
+) where R<:RepresentativeInterval
+    birth, death = birth_death(matrix, column, pivot)
 
     if death - birth > cutoff
-        return PersistenceInterval(birth, death, move!(matrix.working_boundary))
+        return R(
+            PersistenceInterval(birth, death),
+            simplex(pivot), column, move!(matrix.working_boundary)
+        )
     else
         return nothing
     end
@@ -355,9 +363,16 @@ function compute_intervals!(
 ) where {Co, reps}
     if reps
         representative_type = Co ? simplex_element(matrix) : face_element(matrix)
-        intervals = PersistenceInterval{Vector{representative_type}}[]
+        critical_birth_type = simplex_type(matrix.filtration, dim(matrix))
+        critical_death_type = simplex_type(matrix.filtration, dim(matrix) + 1)
+        intervals = RepresentativeInterval{
+            PersistenceInterval,
+            critical_birth_type,
+            Union{critical_death_type, Nothing},
+            Vector{representative_type},
+        }[]
     else
-        intervals = PersistenceInterval{Nothing}[]
+        intervals = PersistenceInterval[]
     end
     if progress
         progbar = Progress(
diff --git a/src/simplexrecipes.jl b/src/simplexrecipes.jl
index b43b4cb0..7c4e9cdc 100644
--- a/src/simplexrecipes.jl
+++ b/src/simplexrecipes.jl
@@ -1,9 +1,3 @@
-# simplex plots ========================================================================== #
-const SxVector{D} = AbstractVector{<:AbstractSimplex{D}}
-const IntervalWithRep{D} = PersistenceInterval{
-    <:AbstractVector{<:Pair{<:AbstractSimplex{D}, <:Any}}
-}
-
 """
     index_data(indices, args...)
 
@@ -35,7 +29,7 @@ end
 
 """
     plottable(sx::AbstractSimplex, args...)
-    plottable(sx::PersistenceInterval, args...)
+    plottable(sx::RepresentativeInterval, args...)
     plottable(sx::Vector{<:AbstractSimplex}, args...)
 
 Turn `sx` into a series that can be plotted. `args...`, should contain the data which tells
@@ -44,23 +38,26 @@ the plot where to place simplex vertices.
 function plottable(sx::AbstractSimplex, args...)
     return plottable([sx], args...)
 end
-function plottable(int::PersistenceInterval, args...)
+function plottable(int::RepresentativeInterval, args...)
     return plottable(simplex.(representative(int)), args...)
 end
-function plottable(int::PersistenceInterval{Nothing}, args...)
-    error("interval has no representative. Run `ripserer` with `representatives=true`")
+function plottable(rep::AbstractVector{<:AbstractChainElement}, args...)
+    return plottable(simplex.(rep), args...)
+end
+function plottable(element::AbstractChainElement, args...)
+    return plottable(simplex(element), args...)
 end
-function plottable(sxs::SxVector{0}, args...)
+function plottable(sxs::AbstractVector{<:AbstractSimplex{0}}, args...)
     indices = only.(vertices.(sxs))
     return index_data(indices, args...), [:seriestype => :scatter], 0
 end
-function plottable(sxs::SxVector{1}, args...)
+function plottable(sxs::AbstractVector{<:AbstractSimplex{1}}, args...)
     indices = mapreduce(vcat, vertices.(sxs)) do (u, v)
         [u, v, 0]
     end
     return index_data(indices, args...), [:seriestype => :path], 1
 end
-function plottable(sxs::SxVector{D}, args...) where D
+function plottable(sxs::AbstractVector{<:AbstractSimplex{D}}, args...) where D
     indices = mapreduce(vcat, vertices.(sxs)) do vs
         idxs = Int[]
         for (u, v, w) in subsets(vs, Val(3))
@@ -75,16 +72,32 @@ end
 function apply_threshold(sx::AbstractSimplex, thresh, thresh_strict)
     return diam(sx) ≤ thresh && diam(sx) < thresh_strict ? sx : nothing
 end
-function apply_threshold(sxs::SxVector, thresh, thresh_strict)
+function apply_threshold(sxs::AbstractVector{<:AbstractSimplex}, thresh, thresh_strict)
     return filter(sx -> diam(sx) ≤ thresh && diam(sx) < thresh_strict, sxs)
 end
-function apply_threshold(int::PersistenceInterval, thresh, thresh_strict)
-    reps = filter(sx -> diam(sx) ≤ thresh && diam(sx) < thresh_strict, representative(int))
-    return PersistenceInterval(birth(int), death(int), reps)
+function apply_threshold(els::AbstractVector{<:AbstractChainElement}, thresh, thresh_strict)
+    return apply_threshold(simplex.(els), thresh, thresh_strict)
+end
+function apply_threshold(int::PersistenceDiagrams.AbstractInterval, thresh, thresh_strict)
+    error("interval has no representative. Run `ripserer` with `representatives=true`")
+end
+function apply_threshold(int::RepresentativeInterval, thresh, thresh_strict)
+    return apply_threshold(representative(int), thresh, thresh_strict)
 end
 
-@recipe function f(sx::Union{AbstractSimplex, SxVector, PersistenceInterval}, args...;
-                   threshold=Inf, threshold_strict=Inf)
+@recipe function f(
+    sx::Union{
+        AbstractSimplex,
+        AbstractChainElement,
+        AbstractVector{<:AbstractSimplex},
+        AbstractVector{<:AbstractChainElement},
+        PersistenceDiagrams.AbstractInterval,
+        RepresentativeInterval,
+    },
+    args...;
+    threshold=Inf,
+    threshold_strict=Inf
+)
     sx = apply_threshold(sx, threshold, threshold_strict)
     isnothing(sx) && return ()
     series, attrs, D = plottable(sx, args...)
diff --git a/src/zerodimensional.jl b/src/zerodimensional.jl
index 5d12c73f..bff33c16 100644
--- a/src/zerodimensional.jl
+++ b/src/zerodimensional.jl
@@ -65,16 +65,31 @@ end
 
 birth(dset::DisjointSetsWithBirth, i) = dset.births[i]
 
-"""
-    zeroth_representative(dset, vertex, reps, CE, V)
+function birth_death(dset::DisjointSetsWithBirth, vertex, edge)
+    return birth(dset, vertex), isnothing(edge) ? Inf : diam(edge)
+end
 
-Collect the zero dimensional representative of `vertex`.
-"""
-function zeroth_representative(filtration, dset, vertex, reps, CE, V)
-    if reps
-        return map(find_leaves!(dset, vertex)) do u
-            CE(V((u,), birth(filtration, u)))
+function interval(
+    ::Type{R}, dset::DisjointSetsWithBirth, filtration, vertex, edge, cutoff
+) where {V, R<:RepresentativeInterval{PersistenceInterval, V}}
+    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)
+    else
+        return nothing
+    end
+end
+
+function interval(
+    ::Type{PersistenceInterval}, dset::DisjointSetsWithBirth, filtration, vertex, edge, cutoff
+)
+    birth, death = birth_death(dset, vertex, edge)
+    if death - birth > cutoff
+        return PersistenceInterval(birth, death)
     else
         return nothing
     end
@@ -97,9 +112,14 @@ function zeroth_intervals(
     CE = chain_element_type(V, F)
     dset = DisjointSetsWithBirth([birth(filtration, v) for v in 1:n_vertices(filtration)])
     if reps
-        intervals = PersistenceInterval{Vector{CE}}[]
+        intervals = RepresentativeInterval{
+            PersistenceInterval,
+            vertex_type(filtration),
+            Union{edge_type(filtration), Nothing},
+            Vector{CE},
+        }[]
     else
-        intervals = PersistenceInterval{Nothing}[]
+        intervals = PersistenceInterval[]
     end
     to_skip = edge_type(filtration)[]
     to_reduce = edge_type(filtration)[]
@@ -107,30 +127,28 @@ function zeroth_intervals(
     if progress
         progbar = Progress(length(simplices), desc="Computing 0d intervals... ")
     end
-    for sx in simplices
-        u, v = index.(vertices(sx))
+    for edge in simplices
+        u, v = index.(vertices(edge))
         i = find_root!(dset, u)
         j = find_root!(dset, v)
         if i ≠ j
             # According to the elder rule, the vertex with the lower birth will fall
             # into a later interval.
-            dead = birth(dset, i) > birth(dset, j) ? i : j
-            if diam(sx) - birth(dset, dead) > cutoff
-                representative = zeroth_representative(filtration, dset, dead, reps, CE, V)
-                interval = PersistenceInterval(birth(dset, dead), diam(sx), representative)
-                push!(intervals, interval)
+            v = birth(dset, i) > birth(dset, j) ? i : j
+            int = interval(eltype(intervals), dset, filtration, v, edge, cutoff)
+            if !isnothing(int)
+                push!(intervals, int)
             end
             union!(dset, i, j)
-            push!(to_skip, sx)
+            push!(to_skip, edge)
         else
-            push!(to_reduce, sx)
+            push!(to_reduce, edge)
         end
         progress && next!(progbar)
     end
     for v in 1:n_vertices(filtration)
         if find_root!(dset, v) == v
-            representative = zeroth_representative(filtration, dset, v, reps, CE, V)
-            push!(intervals, PersistenceInterval(birth(dset, v), Inf, representative))
+            push!(intervals, interval(eltype(intervals), dset, filtration, v, nothing, 0))
         end
     end
     reverse!(to_reduce)
diff --git a/test/cubical.jl b/test/cubical.jl
index 45f877a0..fe6144a8 100644
--- a/test/cubical.jl
+++ b/test/cubical.jl
@@ -1,4 +1,6 @@
 using Ripserer
+
+using Compat
 using StaticArrays
 
 using Ripserer: all_equal_in_dim
@@ -104,7 +106,47 @@ end
         @test isempty(cocob)
     end
 
-    @testset "3d" begin
+    @testset "3d values" begin
+        data = zeros(2, 2, 2)
+        data[:, 1, 1] .= 1
+        data[:, 2, 1] .= 2
+        data[:, 1, 2] .= 3
+        data[:, 2, 2] .= 4
+
+        @testset "Cubelet{1} coboundary" begin
+            cob = Cubelet{2, Float64, Int}[]
+            for cube in edges(Cubical(data))
+                for c in coboundary(Cubical(data), cube)
+                    @test issubset(cube, c)
+                    @test diam(cube) ≤ diam(c)
+                end
+            end
+        end
+        @testset "Cubelet{1} coboundary, all=false" begin
+            cob = Cubelet{2, Float64, Int}[]
+            for cube in edges(Cubical(data))
+                for c in coboundary(Cubical(data), cube, Val(false))
+                    push!(cob, c)
+                end
+            end
+            @test length(cob) == 6
+            @test allunique(cob)
+        end
+        @testset "Cubelet{2} coboundary" begin
+            for cube in [
+                Cubelet{2}((1, 2, 3, 4), 2.0),
+                Cubelet{2}((1, 2, 5, 6), 4.0),
+                Cubelet{2}((1, 3, 5, 7), 3.0),
+                Cubelet{2}((2, 4, 6, 8), 4.0),
+                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)
+            end
+        end
+    end
+    @testset "3d counts" begin
         cob = Cubelet{2, Int, Int}[]
         flt = Cubical(data3d)
         cub = Cubelet{1}((14, 13), 1)
diff --git a/test/plotting.jl b/test/plotting.jl
index 5b6018d3..b7450eb9 100644
--- a/test/plotting.jl
+++ b/test/plotting.jl
@@ -43,13 +43,14 @@ series(args...; kwargs...) =
         data = collect(1:100)
         for dim in 0:3
             sx = Simplex{dim}(1, 1)
+            fsx = Simplex{dim + 1}(1, 1)
             @test length(series(sx, data)) == 1
             @test length(series([sx], data)) == 1
             @test length(series([sx], data)) == 1
-            @test length(series(PersistenceInterval(
-                1.0, 1.0, [ChainElement{typeof(sx), typeof(1//1)}(sx, 1//1)]
+            @test length(series(RepresentativeInterval(
+                1.0, 1.0, sx, fsx, [ChainElement{typeof(sx), typeof(1//1)}(sx, 1//1)]
             ), data)) == 1
-            @test_throws ErrorException series(PersistenceInterval(1.0, 1.0, nothing), data)
+            @test_throws ErrorException series(PersistenceInterval(1.0, 1.0), data)
         end
     end
 end
diff --git a/test/ripserer.jl b/test/ripserer.jl
index 127c6835..f7133a9d 100644
--- a/test/ripserer.jl
+++ b/test/ripserer.jl
@@ -170,55 +170,90 @@ end
     end
 end
 
-@testset "reps" begin
-    _, d1, d2 = ripserer(projective_plane, dim_max=2, reps=true)
-
-    @test simplex.(representative(only(d1))) == [
-        Simplex{1}((11, 10), 1),
-        Simplex{1}((10, 7), 1),
-        Simplex{1}((10, 6), 1),
-        Simplex{1}((8, 1), 1),
-        Simplex{1}((7, 3), 1),
-        Simplex{1}((7, 1), 1),
-        Simplex{1}((6, 2), 1),
-        Simplex{1}((5, 1), 1),
-        Simplex{1}((2, 1), 1),
-    ]
-    @test coefficient.(representative(only(d1))) == fill(Mod{2}(1), 9)
-    @test simplex.(representative(only(d2))) == [Simplex{2}((6, 2, 1), 1)]
-    @test coefficient.(representative(only(d2))) == [Mod{2}(1)]
-end
-
-@testset "reps - types" begin
-    d0, d1, d2, d3 = ripserer(cycle, dim_max=3, reps=true)
-    @test eltype(d0) <: PersistenceInterval{
-        <:Vector{<:PackedElement{Simplex{0, Int, Int}, Mod{2}}}}
-    @test eltype(d1) <: PersistenceInterval{
-        <:Vector{<:PackedElement{Simplex{1, Int, Int}, Mod{2}}}}
-    @test eltype(d2) <: PersistenceInterval{
-        <:Vector{<:PackedElement{Simplex{2, Int, Int}, Mod{2}}}}
-    @test eltype(d3) <: PersistenceInterval{
-        <:Vector{<:PackedElement{Simplex{3, Int, Int}, Mod{2}}}}
-
-    d0, d1, d2, d3 = ripserer(cycle, dim_max=3, reps=true,
-                              field_type=Rational{Int})
-    @test eltype(d0) ≡ PersistenceInterval{
-        Vector{ChainElement{Simplex{0, Int, Int}, Rational{Int}}}}
-    @test eltype(d1) ≡ PersistenceInterval{
-        Vector{ChainElement{Simplex{1, Int, Int}, Rational{Int}}}}
-    @test eltype(d2) ≡ PersistenceInterval{
-        Vector{ChainElement{Simplex{2, Int, Int}, Rational{Int}}}}
-    @test eltype(d3) ≡ PersistenceInterval{
-        Vector{ChainElement{Simplex{3, Int, Int}, Rational{Int}}}}
-end
-
-@testset "Representative of infinite interval." begin
-    _, d1 = ripserer(
-        cycle, dim_max=1, reps=true, threshold=1, field_type=Rational{Int}
-    )
-    @test representative(only(d1)) == ChainElement{
-        Simplex{0, Int, Int}, Rational{Int}
-    }[]
+@testset "Representatives" begin
+    @testset "example from ripser" begin
+        _, d1, d2 = ripserer(projective_plane, dim_max=2, reps=true)
+
+        @test simplex.(representative(only(d1))) == [
+            Simplex{1}((11, 10), 1),
+            Simplex{1}((10, 7), 1),
+            Simplex{1}((10, 6), 1),
+            Simplex{1}((8, 1), 1),
+            Simplex{1}((7, 3), 1),
+            Simplex{1}((7, 1), 1),
+            Simplex{1}((6, 2), 1),
+            Simplex{1}((5, 1), 1),
+            Simplex{1}((2, 1), 1),
+        ]
+        @test coefficient.(representative(only(d1))) == fill(Mod{2}(1), 9)
+        @test simplex.(representative(only(d2))) == [Simplex{2}((6, 2, 1), 1)]
+        @test coefficient.(representative(only(d2))) == [Mod{2}(1)]
+    end
+    @testset "types" begin
+        d0, d1, d2, d3 = ripserer(cycle, dim_max=3, reps=true)
+        @test eltype(d0) <: RepresentativeInterval{
+            PersistenceInterval,
+            Simplex{0, Int, Int},
+            Union{Nothing, Simplex{1, Int, Int}},
+            <:Vector{<:PackedElement{Simplex{0, Int, Int}, Mod{2}}}}
+        @test eltype(d1) <: RepresentativeInterval{
+            PersistenceInterval,
+            Simplex{1, Int, Int},
+            Union{Nothing, Simplex{2, Int, Int}},
+            <:Vector{<:PackedElement{Simplex{1, Int, Int}, Mod{2}}}}
+        @test eltype(d2) <: RepresentativeInterval{
+            PersistenceInterval,
+            Simplex{2, Int, Int},
+            Union{Nothing, Simplex{3, Int, Int}},
+            <:Vector{<:PackedElement{Simplex{2, Int, Int}, Mod{2}}}}
+        @test eltype(d3) <: RepresentativeInterval{
+            PersistenceInterval,
+            Simplex{3, Int, Int},
+            Union{Nothing, Simplex{4, Int, Int}},
+            <:Vector{<:PackedElement{Simplex{3, Int, Int}, Mod{2}}}}
+
+        d0, d1, d2, d3 = ripserer(cycle, dim_max=3, reps=true,
+                                  field_type=Rational{Int})
+        @test eltype(d0) ≡ RepresentativeInterval{
+            PersistenceInterval,
+            Simplex{0, Int, Int},
+            Union{Nothing, Simplex{1, Int, Int}},
+            Vector{ChainElement{Simplex{0, Int, Int}, Rational{Int}}}}
+        @test eltype(d1) ≡ RepresentativeInterval{
+            PersistenceInterval,
+            Simplex{1, Int, Int},
+            Union{Nothing, Simplex{2, Int, Int}},
+            Vector{ChainElement{Simplex{1, Int, Int}, Rational{Int}}}}
+        @test eltype(d2) ≡ RepresentativeInterval{
+            PersistenceInterval,
+            Simplex{2, Int, Int},
+            Union{Nothing, Simplex{3, Int, Int}},
+            Vector{ChainElement{Simplex{2, Int, Int}, Rational{Int}}}}
+        @test eltype(d3) ≡ RepresentativeInterval{
+            PersistenceInterval,
+            Simplex{3, Int, Int},
+            Union{Nothing, Simplex{4, Int, Int}},
+            Vector{ChainElement{Simplex{3, Int, Int}, Rational{Int}}}}
+    end
+    @testset "infinite interval" begin
+        _, d1 = ripserer(
+            cycle, dim_max=1, reps=true, threshold=1, field_type=Rational{Int}
+        )
+        @test representative(only(d1)) == ChainElement{
+            Simplex{0, Int, Int}, Rational{Int}
+        }[]
+    end
+    @testset "critical simplices" begin
+        t = torus(100)
+        result = ripserer(t, reps=true, threshold=0.5)
+        for diag in result
+            @test birth.(diag) == diam.(birth_simplex.(diag))
+            finite = filter(isfinite, diag)
+            @test death.(finite) == diam.(death_simplex.(finite))
+            infinite = filter(!isfinite, diag)
+            @test all(isnothing, death_simplex.(infinite))
+        end
+    end
 end
 
 @testset "Sublevel using SparseRips" begin
@@ -246,35 +281,66 @@ end
     @test sort(maxs) == [1.0, 2.0]
 end
 
-@testset "2D image sublevel using Cubical" begin
-    data = [0 0 0 0 0;
-            0 2 2 2 0;
-            0 2 1 2 0;
-            0 2 2 2 0;
-            0 0 0 0 0]
+@testset "Cubical" begin
+    @testset "2D image" begin
+        data = [0 0 0 0 0;
+                0 2 2 2 0;
+                0 2 1 2 0;
+                0 2 2 2 0;
+                0 0 0 0 0]
 
-    d0, d1, d2 = ripserer(Cubical(data), reps=true, dim_max=2)
+        d0, d1, d2 = ripserer(Cubical(data), reps=true, dim_max=2)
 
-    @test d0 == [(0, Inf), (1, 2)]
-    @test d1 == [(0, 2)]
-    @test isempty(d2)
+        @test d0 == [(0, Inf), (1, 2)]
+        @test d1 == [(0, 2)]
+        @test isempty(d2)
 
-    @test vertices.(representative(d0[1])) == [SVector(i) for i in 1:length(data)]
-    @test vertices(only(representative(d0[2]))) == SVector(13)
+        @test vertices.(representative(d0[1])) == [SVector(i) for i in 1:length(data)]
+        @test vertices(only(representative(d0[2]))) == SVector(13)
+    end
+    @testset "3D image" begin
+        # Cube with hole in the middle.
+        data = zeros(5, 5, 5)
+        data[2, 2:4, 2:4] .= 1
+        data[3, :, :] .= [0 0 0 0 0; 0 1 1 1 0; 0 1 0 1 0; 0 1 1 1 0; 0 0 0 0 0]
+        data[4, 2:4, 2:4] .= 1
+
+        d0, d1, d2 = ripserer(Cubical(data), dim_max=2)
+
+        @test d0 == [(0, 1.0), (0, Inf)]
+        @test isempty(d1)
+        @test d2 == [(0, 1)]
+    end
 end
 
 @testset "Persistent homology." begin
-    res_hom = ripserer(cycle, cohomology=false, dim_max=3)
-    res_coh = ripserer(cycle, cohomology=true, dim_max=3)
-
-    @test res_hom == res_coh
-
-    res_hom = ripserer(cycle, cohomology=false, reps=true, dim_max=3)
-    @test vertices.(simplex.(representative(res_hom[2][1]))) == sort!(vcat(
-        [SVector(i+1, i) for i in 1:17], [SVector(18, 1)]
-    ))
+    @testset "same diagrams" begin
+        res_hom = ripserer(cycle, cohomology=false, dim_max=3)
+        res_coh = ripserer(cycle, cohomology=true, dim_max=3)
 
-    @test_broken ripserer(cycle, cohomology=false, threshold=2)[2][1] == (1.0, Inf)
+        @test res_hom == res_coh
+    end
+    @testset "same critical simplices" begin
+        # Add some noise because critical simplices might be different if values are exactly
+        # the same.
+        cyc = cycle .+ 0.01 .* rand_dist_matrix(18)
+        res_hom = ripserer(cyc, cohomology=false, dim_max=3, reps=true)
+        res_coh = ripserer(cyc, cohomology=true, dim_max=3, reps=true)
+
+        for i in 1:4
+            @test birth_simplex.(res_hom[i]) == birth_simplex.(res_coh[i])
+            @test death_simplex.(res_hom[i]) == death_simplex.(res_coh[i])
+        end
+    end
+    @testset "representative cycle" begin
+        res_hom = ripserer(cycle, cohomology=false, reps=true, dim_max=3)
+        @test vertices.(simplex.(representative(res_hom[2][1]))) == sort!(vcat(
+            [SVector(i+1, i) for i in 1:17], [SVector(18, 1)]
+        ))
+    end
+    @testset "infinite intervals" begin
+        @test_broken ripserer(cycle, cohomology=false, threshold=2)[2][1] == (1.0, Inf)
+    end
 end
 
 @testset "Only print to stderr and only when progress is enabled." begin