Fix terrible performance with custom filtrations (#171)

* Fix terrible performance with custom filtrations * Fix type params, prevent unnecessary overflows in high-dimensional simplices * fixup aqua tests * update changelog * Fix aqua and formatting * Update doctests * Update test/aqua.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update test/doctests.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * don't fix in doctest * revert doctests * cleanup Project.toml * remove future * Add assertion to Mod * Fixup docs * Format --------- Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
mtsch · Jun 4, 2024 · 8b296a0 · 8b296a0 · mtsch · Jun 4, 2024
1 parent 3f2d63f
commit 8b296a0
Show file tree

Hide file tree

Showing 13 changed files with 52 additions and 246 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,3 +1,7 @@
+# v0.16.13
+
+Fix performance issues with high-dimensional `Custom` filtrations.
+
 # v0.16.9
 
 Replace LightGraps with Graphs.

diff --git a/Project.toml b/Project.toml
@@ -1,13 +1,12 @@
 name = "Ripserer"
 uuid = "aa79e827-bd0b-42a8-9f10-2b302677a641"
 authors = ["mtsch <[email protected]>"]
-version = "0.16.12"
+version = "0.16.13"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
-Future = "9fa8497b-333b-5362-9e8d-4d0656e87820"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -21,17 +20,26 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
 
 [compat]
+Aqua = "0.8"
 Compat = "^3.10, 4"
 DataStructures = "0.17, 0.18"
 Distances = "0.8, 0.9, 0.10"
+Documenter = "1"
 Graphs = "1"
 IterTools = "1"
+LinearAlgebra = "1"
+MLJBase = "1"
 MLJModelInterface = "^0.3.5, 0.4, 1"
 MiniQhull = "0.2, 0.3, 0.4"
 PersistenceDiagrams = "0.9"
 ProgressMeter = "1"
+Random = "1"
 RecipesBase = "1"
+SafeTestsets = "0.1"
+SparseArrays = "1"
 StaticArrays = "0.12, 1"
+Suppressor = "0.2"
+Test = "1"
 TupleTools = "1"
 julia = "1.6"
 

diff --git a/docs/make.jl b/docs/make.jl
@@ -28,12 +28,8 @@ makedocs(;
         "Home" => "index.md",
         "Usage Guide" => "generated/basics.md",
         "API" => "api.md",
-        "Examples" => [
-            "generated/stability.md",
-            "generated/cocycles.md",
-            "generated/cubical.md",
-            "generated/malaria.md",
-        ],
+        "Examples" =>
+            ["generated/stability.md", "generated/cocycles.md", "generated/cubical.md"],
         "Benchmarks" => "benchmarks.md",
         "Related Julia Packages" => "related-work.md",
         "Acknowledgements and References" => "references.md",

diff --git a/docs/src/examples/malaria.jl b/docs/src/examples/malaria.jl
diff --git a/src/base/abstractcell.jl b/src/base/abstractcell.jl
@@ -13,9 +13,8 @@ is determined by the following.
 * The `birth` of type `T` determines when the cell enters the filtration. Note that two
   cells with the same `index` should also have the same `birth`.
 
-* `vertices` should return the cell's vertices as a tuple. For 0-cells (vertices), the
-  `vertices` are also used to index into a filtration's
-  [`vertices`](vertices(::AbstractFiltration)).
+* `vertices` should return the cell's vertices as a tuple. For 0-cells, these are also used
+  to index into a filtration's [`vertices`](@ref vertices(::AbstractFiltration)).
 
 * The `sign` determens its orientation. Note that `cell == -cell`.
 

diff --git a/src/base/abstractfiltration.jl b/src/base/abstractfiltration.jl
@@ -162,10 +162,7 @@ julia> vertices(Rips([0 1 1; 1 0 1; 1 1 0]))
 Base.OneTo(3)
 
 julia> vertices(Cubical([0 1 1; 1 0 1; 1 1 0]))
-3×3 CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}:
- CartesianIndex(1, 1)  CartesianIndex(1, 2)  CartesianIndex(1, 3)
- CartesianIndex(2, 1)  CartesianIndex(2, 2)  CartesianIndex(2, 3)
- CartesianIndex(3, 1)  CartesianIndex(3, 2)  CartesianIndex(3, 3)
+CartesianIndices((3, 3))
 
 ```
 """

diff --git a/src/base/abstractsimplex.jl b/src/base/abstractsimplex.jl
@@ -65,7 +65,11 @@ Use binary search to find index of first vertex in `(K - 1)`-dimensional simplex
 """
 function _first_vertex(index::I, ::Val{K}) where {I,K}
     lo = I(K - 1)
-    hi = I(K + 100)
+    if K < 3
+        hi = I(K + 100)
+    else
+        hi = I(K)
+    end
     while _binomial(hi, Val(K)) ≤ index
         lo = hi
         hi <<= 0x01

diff --git a/src/base/chain.jl b/src/base/chain.jl
@@ -22,7 +22,7 @@ julia> data = [(rand(), rand(), rand()) for _ in 1:100];
 julia> result = ripserer(data; reps=true, modulus=7);
 
 julia> chain = result[end][end].representative
-28-element Chain{Ripserer.Mod{7},Ripserer.Simplex{1, Float64, Int64}}:
+28-element Chain{Mod{7},Simplex{1, Float64, Int64}}:
  +Simplex{1}((68, 54), 0.25178097927369875) => 6 mod 7
  +Simplex{1}((54, 46), 0.2575262844682746) => 1 mod 7
  +Simplex{1}((88, 56), 0.25936896586973557) => 1 mod 7

diff --git a/src/base/primefield.jl b/src/base/primefield.jl
@@ -27,7 +27,12 @@ is_prime(::Any) = false
 Like `mod`, but with prime `M`.
 """
 function mod_prime(i, ::Val{M}) where {M}
-    is_prime(M) || throw(DomainError(M, "modulus must be a prime number"))
+    if !is_prime(M)
+        throw(DomainError(M, "modulus must be a prime number"))
+    end
+    if !(1 < M ≤ 3037000499)
+        throw(DomainError(M, "modulus must be smaller than √(typemax(Int))"))
+    end
     i = i % M
     return i + ifelse(signbit(i), M, 0)
 end

diff --git a/src/extra/circularcoordinates.jl b/src/extra/circularcoordinates.jl
@@ -192,7 +192,7 @@ CircularCoordinates(
   radius=(0.9510565162951535,),
   n_landmarks=10,
   partition=linear,
-  metric=Distances.Euclidean(0.0),
+  metric=Euclidean(0.0),
 )
 
 julia> summary(cc(data))

diff --git a/src/filtrations/alpha.jl b/src/filtrations/alpha.jl
@@ -180,13 +180,13 @@ julia> length(Ripserer.edges(rips))
 
 julia> sort(ripserer(alpha)[2], by=persistence)[end]
 [0.375, 2.01) with:
- birth_simplex: Ripserer.Simplex{1, Float64, Int64}
- death_simplex: Ripserer.Simplex{2, Float64, Int64}
+ birth_simplex: Simplex{1, Float64, Int64}
+ death_simplex: Simplex{2, Float64, Int64}
 
 julia> sort(ripserer(rips)[2], by=persistence)[end]
 [0.375, 2.01) with:
- birth_simplex: Ripserer.Simplex{1, Float64, Int64}
- death_simplex: Ripserer.Simplex{2, Float64, Int64}
+ birth_simplex: Simplex{1, Float64, Int64}
+ death_simplex: Simplex{2, Float64, Int64}
 ```
 """
 struct Alpha{I,P<:SVector} <: AbstractCustomFiltration{I,Float64}

diff --git a/src/filtrations/custom.jl b/src/filtrations/custom.jl
@@ -100,8 +100,7 @@ struct Custom{I,T} <: AbstractCustomFiltration{I,T}
     threshold::T
 end
 
-@inline insert_simplex!(::Vector{Dict{I,T}}, ::Tuple{}, _, _) where {I,T} = nothing
-@inline function insert_simplex!(
+function insert_simplex!(
     dicts::Vector{Dict{I,T}}, vertices::NTuple{N,I}, birth, threshold
 ) where {N,T,I}
     if birth > threshold
@@ -111,11 +110,19 @@ end
         dim = N - 1
         d_dict = dicts[dim + 1]
         d_dict[idx] = min(birth, get(d_dict, idx, typemax(T)))
-        for vs in IterTools.subsets(vertices, Val(N - 1))
-            # Some extra work is being done here, but this tends to be type stable.
-            insert_simplex!(dicts, vs, birth, threshold)
-        end
+        _insert_simplex_facets!(dicts, vertices, birth, Val(N - 1))
+    end
+end
+@inline _insert_simplex_facets!(_, _, _, ::Val{0}) = nothing
+@inline function _insert_simplex_facets!(dicts, vertices, birth, ::Val{N}) where {N}
+    dim = N - 1
+    d_dict = dicts[dim + 1]
+    T = valtype(d_dict)
+    for vs in IterTools.subsets(vertices, Val(N))
+        idx = index(vs)
+        d_dict[idx] = min(birth, get(d_dict, idx, typemax(T)))
     end
+    return _insert_simplex_facets!(dicts, vertices, birth, Val(N - 1))
 end
 
 function _adjacency_matrix(dicts)