8545078607

Committed 03 Apr 2024 08:26PM UTC coverage: 96.51% (+0.09%) from 96.417%

Build # 8545078607

Build Type

Pull #166

github

Committed by

web-flow

Commit Message

Merge ed46e50d4 into 3f2d63fae

Pull Request Pull Request #166: Better homology

Run Details

245 of 250 new or added lines in 5 files covered. (98.0%)

3 existing lines in 3 files now uncovered.

1770 of 1834 relevant lines covered (96.51%)

3302629.22 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

99.11

/src/filtrations/rips.jl

"""
    AbstractRipsFiltration{I<:Signed, T} <: AbstractFiltration{I, T}

An abstract Vietoris-Rips filtration. Its subtypes can only overload
[`adjacency_matrix`](@ref) and get default implementations for the rest of the filtration
interface.

# Example

```jldoctest
julia> struct MyRips <: Ripserer.AbstractRipsFiltration{Int, Float16} end

julia> Ripserer.adjacency_matrix(::MyRips) = [0 1 1; 1 0 1; 1 1 0]

julia> ripserer(MyRips())
2-element Vector{PersistenceDiagrams.PersistenceDiagram}:
 3-element 0-dimensional PersistenceDiagram
 0-element 1-dimensional PersistenceDiagram

```
"""
abstract type AbstractRipsFiltration{I<:Signed,T} <: AbstractFiltration{I,T} end

nv(rips::AbstractRipsFiltration) = size(adjacency_matrix(rips), 1)

function births(rips::AbstractRipsFiltration)
    adj = adjacency_matrix(rips)
    return view(adj, diagind(adj))
end

simplex_type(::Type{<:AbstractRipsFiltration{I,T}}, D) where {I,T} = Simplex{D,T,I}

function edges(rips::AbstractRipsFiltration)
    if issparse(adjacency_matrix(rips))
        _sparse_edges(rips)
    else
        _full_edges(rips)
    end
end

function _full_edges(rips::AbstractRipsFiltration)
    result = edge_type(rips)[]
    @inbounds for u in 1:nv(rips), v in (u + 1):nv(rips)
        sx = unsafe_simplex(rips, Val(1), (v, u))
        !isnothing(sx) && push!(result, sx)
    end
    return result
end

function _sparse_edges(rips::AbstractRipsFiltration)
    adj = adjacency_matrix(rips)
    result = edge_type(rips)[]
    rows = rowvals(adj)
    for u in 1:nv(rips)
        for i in nzrange(adj, u)
            v = rows[i]
            if v > u
                sx = unsafe_simplex(rips, Val(1), (v, u))
                !isnothing(sx) && push!(result, sx)
            end
        end
    end
    return result
end

@inline @propagate_inbounds function unsafe_simplex(
    ::Type{S}, rips::AbstractRipsFiltration{I,T}, vertices
) where {I,T,S<:Simplex}
    if dim(S) == 0
        return S(vertices[1], births(rips)[vertices[1]])
    else
        adj = adjacency_matrix(rips)
        n = length(vertices)
        diameter = typemin(T)
        for i in 1:n, j in (i + 1):n
            e = adj[vertices[j], vertices[i]]
            (iszero(e) || e > threshold(rips)) && return nothing
            diameter = ifelse(e > diameter, e, diameter)
        end
        return S(index(vertices), diameter)
    end
end

@inline @propagate_inbounds function unsafe_cofacet(
    ::Type{S}, rips::AbstractRipsFiltration{I,T}, simplex, cofacet_vertices, new_vertex
) where {I,T,S<:Simplex}
    diameter = birth(simplex)
    adj = adjacency_matrix(rips)
    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 new_vertex and v.
        e = adj[new_vertex, v]
        (iszero(e) || e > threshold(rips)) && return nothing
        diameter = ifelse(e > diameter, e, diameter)
    end
    return S(index(cofacet_vertices), diameter)
end

# This definition is used in SparseCoboundary
@inline @propagate_inbounds function unsafe_cofacet(
    ::Type{S}, rips::AbstractRipsFiltration{I}, simplex, cofacet_vertices, _, new_edges
) where {I,S<:Simplex}
    diameter = birth(simplex)
    for e in new_edges
        e > threshold(rips) && return nothing
        diameter = ifelse(e > diameter, e, diameter)
    end
    return S(index(cofacet_vertices), diameter)
end

function coboundary(
    rips::AbstractRipsFiltration, sx::AbstractSimplex, ::Val{A}=Val(true)
) where {A}
    if adjacency_matrix(rips) isa SparseMatrixCSC
        return SparseCoboundary{A}(rips, sx)
    else
        return Coboundary{A}(rips, sx)
    end
end

distance_matrix(rips::AbstractRipsFiltration) = adjacency_matrix(rips)

function _check_distance_matrix(dist::SparseMatrixCSC)
    n = LinearAlgebra.checksquare(dist)
    vals = nonzeros(dist)
    rows = rowvals(dist)
    for i in 1:n
        vertex_birth = dist[i, i]
        for j in nzrange(dist, i)
            i == rows[j] && continue
            edge_birth = vals[j]
            iszero(edge_birth) && throw(ArgumentError("zero edges in input matrix"))
            edge_birth < vertex_birth && throw(
                ArgumentError(
                    "edges with birth value lower than vertex births in input matrix"
                ),
            )
            edge_birth ≠ dist[i, rows[j]] &&
                throw(ArgumentError("input matrix not symmetric"))
        end
    end
end

function _check_distance_matrix(dist::AbstractMatrix)
    n = LinearAlgebra.checksquare(dist)
    for i in 1:n
        for j in (i + 1):n
            vertex_birth = max(dist[j, j], dist[i, i])
            edge_birth = dist[j, i]
            iszero(edge_birth) && throw(ArgumentError("zero edges in input matrix"))
            edge_birth < vertex_birth && throw(
                ArgumentError(
                    "edges with birth value lower than vertex births in input matrix"
                ),
            )
            edge_birth ≠ dist[i, j] && throw(ArgumentError("input matrix not symmetric"))
        end
    end
end

"""
    Rips{I, T} <: AbstractRipsFiltration{I, T}

This type represents a filtration of Vietoris-Rips complexes.

Diagonal items in the input matrix are treated as vertex birth times.

Zero values are not allowed due to how sparse matrices work in Julia. If you need zero birth
times, try offseting all values by a constant.

Threshold defaults to the radius of the input space. When using low `threshold`s, consider
using the `sparse=true` keyword argument. It will give the same result, but may be much
faster.

# Constructors

* `Rips(distance_matrix; threshold=nothing)`
* `Rips(points; metric=Euclidean(1e-12), threshold=nothing)`
* `Rips{I}(args...)`: `I` sets the size of integer used to represent simplices.

# Examples

```jldoctest
julia> data = [(sin(t), cos(t)) for t in range(0, 2π, length=101)][1:end-1];

julia> ripserer(Rips(data))
2-element Vector{PersistenceDiagrams.PersistenceDiagram}:
 100-element 0-dimensional PersistenceDiagram
 1-element 1-dimensional PersistenceDiagram

julia> ripserer(Rips(data, threshold=1.7))[2]
1-element 1-dimensional PersistenceDiagram:
 [0.0628, ∞)

julia> using Distances

julia> ripserer(Rips(data, metric=Cityblock()))[2]
1-element 1-dimensional PersistenceDiagram:
 [0.0888, 2.0)

```
"""
struct Rips{I,T,A<:AbstractMatrix{T}} <: AbstractRipsFiltration{I,T}
    adj::A
    threshold::T
end

function Rips{I}(
    dists::AbstractMatrix{T}; threshold=nothing, sparse=false, verbose=false
) where {I,T}
    _check_distance_matrix(dists)
    thresh = isnothing(threshold) ? radius(dists) : T(threshold)
    if sparse
        adj = SparseArrays.sparse(dists)
    else
        adj = copy(dists)
    end
    if issparse(adj)
        adj isa SparseMatrixCSC ||
            throw(ArgumentError("only SparseMatrixCSC sparse matrices supported"))

        # This is janky, but using the deprecated syntax prints depwarns in the docstrings.
        if VERSION ≥ v"1.9.0-DEV"
            SparseArrays.fkeep!((_, _, v) -> v ≤ thresh, adj)
        else
            SparseArrays.fkeep!(adj, (_, _, v) -> v ≤ thresh)
        end
    end
    return Rips{I,T,typeof(adj)}(adj, thresh)
end
function Rips{I}(points::AbstractVector; metric=Euclidean(1e-12), kwargs...) where {I}
    if !allunique(points)
        @warn "points not unique"
        points = unique(points)
    end
    return Rips{I}(distances(points, metric); kwargs...)
end
function Rips(dist_or_points; kwargs...)
    return Rips{Int}(dist_or_points; kwargs...)
end

function Base.show(io::IO, rips::Rips{I,T}) where {I,T}
    return print(io, "Rips{$I, $T}(nv=$(nv(rips)), sparse=$(issparse(rips.adj)))")
end

threshold(rips::Rips) = rips.threshold
adjacency_matrix(rips::Rips) = rips.adj

1	"""
2	AbstractRipsFiltration{I<:Signed, T} <: AbstractFiltration{I, T}
3
4	An abstract Vietoris-Rips filtration. Its subtypes can only overload
5	[`adjacency_matrix`](@ref) and get default implementations for the rest of the filtration
6	interface.
7
8	# Example
9
10	```jldoctest
11	julia> struct MyRips <: Ripserer.AbstractRipsFiltration{Int, Float16} end
12
13	julia> Ripserer.adjacency_matrix(::MyRips) = [0 1 1; 1 0 1; 1 1 0]
14
15	julia> ripserer(MyRips())
16	2-element Vector{PersistenceDiagrams.PersistenceDiagram}:
17	3-element 0-dimensional PersistenceDiagram
18	0-element 1-dimensional PersistenceDiagram
19
20	```
21	"""
22	abstract type AbstractRipsFiltration{I<:Signed,T} <: AbstractFiltration{I,T} end
23
24	nv(rips::AbstractRipsFiltration) = size(adjacency_matrix(rips), 1)	1,566,498✔
25
26	function births(rips::AbstractRipsFiltration)	49,328✔
27	adj = adjacency_matrix(rips)	49,328✔
28	return view(adj, diagind(adj))	49,328✔
29	end
30
31	simplex_type(::Type{<:AbstractRipsFiltration{I,T}}, D) where {I,T} = Simplex{D,T,I}	43,274,929✔
32
33	function edges(rips::AbstractRipsFiltration)	314✔
34	if issparse(adjacency_matrix(rips))	314✔
35	_sparse_edges(rips)	79✔
36	else
37	_full_edges(rips)	235✔
38	end
39	end
40
41	function _full_edges(rips::AbstractRipsFiltration)	235✔
42	result = edge_type(rips)[]	235✔
43	@inbounds for u in 1:nv(rips), v in (u + 1):nv(rips)	470✔
44	sx = unsafe_simplex(rips, Val(1), (v, u))	493,354✔
45	!isnothing(sx) && push!(result, sx)	785,559✔
46	end	501,531✔
47	return result	235✔
48	end
49
50	function _sparse_edges(rips::AbstractRipsFiltration)	79✔
51	adj = adjacency_matrix(rips)	79✔
52	result = edge_type(rips)[]	79✔
53	rows = rowvals(adj)	79✔
54	for u in 1:nv(rips)	79✔
55	for i in nzrange(adj, u)	2,691✔
56	v = rows[i]	52,021✔
57	if v > u	52,021✔
58	sx = unsafe_simplex(rips, Val(1), (v, u))	25,805✔
59	!isnothing(sx) && push!(result, sx)	51,432✔
60	end
61	end	101,751✔
62	end	4,903✔
63	return result	79✔
64	end
65
66	@inline @propagate_inbounds function unsafe_simplex(	778,850✔
67	::Type{S}, rips::AbstractRipsFiltration{I,T}, vertices
68	) where {I,T,S<:Simplex}
69	if dim(S) == 0	778,850✔
70	return S(vertices[1], births(rips)[vertices[1]])	51,321✔
71	else
72	adj = adjacency_matrix(rips)	729,747✔
73	n = length(vertices)	729,747✔
74	diameter = typemin(T)	729,747✔
75	for i in 1:n, j in (i + 1):n	1,258,167✔
76	e = adj[vertices[j], vertices[i]]	1,240,385✔
77	(iszero(e) \|\| e > threshold(rips)) && return nothing	2,428,124✔
78	diameter = ifelse(e > diameter, e, diameter)	1,012,735✔
79	end	1,759,055✔
80	return S(index(vertices), diameter)	528,420✔
81	end
82	end
83
84	@inline @propagate_inbounds function unsafe_cofacet(	17,206,623✔
85	::Type{S}, rips::AbstractRipsFiltration{I,T}, simplex, cofacet_vertices, new_vertex
86	) where {I,T,S<:Simplex}
87	diameter = birth(simplex)	17,206,623✔
88	adj = adjacency_matrix(rips)	17,206,623✔
89	for v in cofacet_vertices	17,206,623✔
90	v == new_vertex && continue	46,460,930✔
91	# Even though this looks like a tight loop, v changes way more often than us, so
92	# this is the faster order of indexing by new_vertex and v.
93	e = adj[new_vertex, v]	30,758,364✔
94	(iszero(e) \|\| e > threshold(rips)) && return nothing	30,758,364✔
95	diameter = ifelse(e > diameter, e, diameter)	25,282,547✔
96	end	52,715,919✔
97	return S(index(cofacet_vertices), diameter)	11,730,806✔
98	end
99
100	# This definition is used in SparseCoboundary
101	@inline @propagate_inbounds function unsafe_cofacet(	3,511,618✔
102	::Type{S}, rips::AbstractRipsFiltration{I}, simplex, cofacet_vertices, _, new_edges
103	) where {I,S<:Simplex}
104	diameter = birth(simplex)	3,511,618✔
105	for e in new_edges	3,511,618✔
106	e > threshold(rips) && return nothing	10,617,285✔
107	diameter = ifelse(e > diameter, e, diameter)	10,617,235✔
108	end	17,722,902✔
109	return S(index(cofacet_vertices), diameter)	3,511,568✔
110	end
111
112	function coboundary(	2,139,698✔
113	rips::AbstractRipsFiltration, sx::AbstractSimplex, ::Val{A}=Val(true)
114	) where {A}
115	if adjacency_matrix(rips) isa SparseMatrixCSC	2,139,698✔
116	return SparseCoboundary{A}(rips, sx)	530,784✔
117	else
118	return Coboundary{A}(rips, sx)	591,424✔
119	end
120	end
121
122	distance_matrix(rips::AbstractRipsFiltration) = adjacency_matrix(rips)	47✔
123
124	function _check_distance_matrix(dist::SparseMatrixCSC)	10✔
125	n = LinearAlgebra.checksquare(dist)	10✔
126	vals = nonzeros(dist)	10✔
127	rows = rowvals(dist)	10✔
128	for i in 1:n	10✔
129	vertex_birth = dist[i, i]	920✔
130	for j in nzrange(dist, i)	460✔
131	i == rows[j] && continue	40,501✔
132	edge_birth = vals[j]	40,090✔
133	iszero(edge_birth) && throw(ArgumentError("zero edges in input matrix"))	40,090✔
134	edge_birth < vertex_birth && throw(	40,090✔
135	ArgumentError(
136	"edges with birth value lower than vertex births in input matrix"
137	),
138	)
139	edge_birth ≠ dist[i, rows[j]] &&	40,090✔
140	throw(ArgumentError("input matrix not symmetric"))
141	end	80,542✔
142	end	460✔
143	end
144
145	function _check_distance_matrix(dist::AbstractMatrix)	210✔
146	n = LinearAlgebra.checksquare(dist)	212✔
147	for i in 1:n	208✔
148	for j in (i + 1):n	8,357✔
149	vertex_birth = max(dist[j, j], dist[i, i])	921,370✔
150	edge_birth = dist[j, i]	921,370✔
151	iszero(edge_birth) && throw(ArgumentError("zero edges in input matrix"))	921,370✔
152	edge_birth < vertex_birth && throw(	921,370✔
153	ArgumentError(
154	"edges with birth value lower than vertex births in input matrix"
155	),
156	)
157	edge_birth ≠ dist[i, j] && throw(ArgumentError("input matrix not symmetric"))	921,368✔
158	end	1,834,787✔
159	end	8,149✔
160	end
161
162	"""
163	Rips{I, T} <: AbstractRipsFiltration{I, T}
164
165	This type represents a filtration of Vietoris-Rips complexes.
166
167	Diagonal items in the input matrix are treated as vertex birth times.
168
169	Zero values are not allowed due to how sparse matrices work in Julia. If you need zero birth
170	times, try offseting all values by a constant.
171
172	Threshold defaults to the radius of the input space. When using low `threshold`s, consider
173	using the `sparse=true` keyword argument. It will give the same result, but may be much
174	faster.
175
176	# Constructors
177
178	* `Rips(distance_matrix; threshold=nothing)`
179	* `Rips(points; metric=Euclidean(1e-12), threshold=nothing)`
180	* `Rips{I}(args...)`: `I` sets the size of integer used to represent simplices.
181
182	# Examples
183
184	```jldoctest
185	julia> data = [(sin(t), cos(t)) for t in range(0, 2π, length=101)][1:end-1];
186
187	julia> ripserer(Rips(data))
188	2-element Vector{PersistenceDiagrams.PersistenceDiagram}:
189	100-element 0-dimensional PersistenceDiagram
190	1-element 1-dimensional PersistenceDiagram
191
192	julia> ripserer(Rips(data, threshold=1.7))[2]
193	1-element 1-dimensional PersistenceDiagram:
194	[0.0628, ∞)
195
196	julia> using Distances
197
198	julia> ripserer(Rips(data, metric=Cityblock()))[2]
199	1-element 1-dimensional PersistenceDiagram:
200	[0.0888, 2.0)
201
202	```
203	"""
204	struct Rips{I,T,A<:AbstractMatrix{T}} <: AbstractRipsFiltration{I,T}
205	adj::A	214✔
206	threshold::T
207	end
208
209	function Rips{I}(	440✔
210	dists::AbstractMatrix{T}; threshold=nothing, sparse=false, verbose=false
211	) where {I,T}
212	_check_distance_matrix(dists)	220✔
213	thresh = isnothing(threshold) ? radius(dists) : T(threshold)	354✔
214	if sparse	214✔
215	adj = SparseArrays.sparse(dists)	9✔
216	else
217	adj = copy(dists)	205✔
218	end
219	if issparse(adj)	223✔
220	adj isa SparseMatrixCSC \|\|	19✔
221	throw(ArgumentError("only SparseMatrixCSC sparse matrices supported"))
222
223	# This is janky, but using the deprecated syntax prints depwarns in the docstrings.
224	if VERSION ≥ v"1.9.0-DEV"	19✔
225	SparseArrays.fkeep!((_, _, v) -> v ≤ thresh, adj)	90,350✔
226	else
NEW 227	SparseArrays.fkeep!(adj, (_, _, v) -> v ≤ thresh)	×
228	end
229	end
230	return Rips{I,T,typeof(adj)}(adj, thresh)	214✔
231	end
232	function Rips{I}(points::AbstractVector; metric=Euclidean(1e-12), kwargs...) where {I}	150✔
233	if !allunique(points)	119✔
234	@warn "points not unique"	2✔
235	points = unique(points)	2✔
236	end
237	return Rips{I}(distances(points, metric); kwargs...)	75✔
238	end
239	function Rips(dist_or_points; kwargs...)	386✔
240	return Rips{Int}(dist_or_points; kwargs...)	193✔
241	end
242
243	function Base.show(io::IO, rips::Rips{I,T}) where {I,T}	147✔
244	return print(io, "Rips{$I, $T}(nv=$(nv(rips)), sparse=$(issparse(rips.adj)))")	147✔
245	end
246
247	threshold(rips::Rips) = rips.threshold	39,599,655✔
248	adjacency_matrix(rips::Rips) = rips.adj	23,382,590✔

mtsch / Ripserer.jl / 8545078607

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous