New Way of Getting Simplices From Filtrations (#60)

* new way of getting simplices from filtrations

* sparse filtration performance improvement

.travis.yml

@@ -5,7 +5,7 @@ os:
 
 julia:
   - 1.0
-  - 1.4
+  - 1
   - nightly
 
 notifications:

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"

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

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)
 ```

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
 

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.

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