`Custom` filtration gives an `OverflowError`

[lampretl]

I have the facets fct (maximal simplices) of the 8-dimensional real projective space RP8, obtained from SageMath. I wish to benchmark Ripser computing its homology (as a trivial filtration where all simplices are in level 0).

julia> typeof(fct), length(fct), extrema(v for σ∈fct for v∈σ)
(Vector{Vector{Int64}}, 1814400, (1, 511))

julia> flt = Ripserer.Custom([Tuple(σ)=>1 for σ∈fct]; threshold=1);
ERROR: OverflowError: simplex (511, 510, 368, 333, 328, 217, 214, 83, 36) overflows in UniformScaling{Bool}(true)! Try using a bigger index type in your filtration.
Stacktrace:
 [1] index_overflow_check(vertices::NTuple{9, Int64}, message::String)
   @ Ripserer ~/.julia/packages/Ripserer/r6Usj/src/base/abstractsimplex.jl:156
 [2] index_overflow_check
   @ ~/.julia/packages/Ripserer/r6Usj/src/base/abstractsimplex.jl:153 [inlined]
 [3] (Custom{Int64})(simplices::Vector{Pair{NTuple{9, Int64}, Int64}}; threshold::Int64, verbose::Bool)
   @ Ripserer ~/.julia/packages/Ripserer/r6Usj/src/filtrations/custom.jl:177
 [4] Custom
   @ ~/.julia/packages/Ripserer/r6Usj/src/filtrations/custom.jl:158 [inlined]
 [5] #Custom#65
   @ ~/.julia/packages/Ripserer/r6Usj/src/filtrations/custom.jl:181 [inlined]
 [6] top-level scope
   @ REPL[30]:1

I'm not sure how to fix this error.

[mtsch]

This happens because the index of the simplex is too large to fit into a 64-bit interger. You can change the index type used in a filtration by giving it a type argument as

flt = Ripserer.Custom{Int128}([Tuple(σ)=>1 for σ∈fct]; threshold=1);

If Int128 still errors, you can try some larger integer types from BitIntegers.jl and it should work, but will probably be a bit slower. If all else fails, you can also use the built-in BigInt, but that will probably be extremely slow.

[lampretl]

Thank you for a quick response! What is the index of a simplex? Something similar to a hash? I am asking because in my case, the number of vertices is very small (511), so the required type is quite small (UInt16).

[lampretl]

Hmm, if on RP4 I use Ripserer.Custom{UInt64}([Tuple(σ)=>1 for σ∈fct]; threshold=1), then Ripserer.ripserer(...) just hangs and doesn't finish:

[mtsch]

The index maps the vertices of a simplex to its combinatorial index, i.e. the position it would have in a list of all possible N-simplices ordered lexicographically. For example, the first few 2-simplices have the following indices

vertices <-> index
(3,2,1)  <-> 1
(4,2,1)  <-> 2
(4,3,1)  <-> 3
(4,3,2)  <-> 4
(5,2,1)  <-> 5
(5,3,1)  <-> 6
(5,3,2)  <-> 7
(5,4,1)  <-> 8
...

This is used for compression as it allows us to store a many-vertex simplex in a single integer. If you want to know more, the algorithms for converting indices to and from simplices is described in the original Ripser paper.
In Ripserer, the sign of the index is used internally to denote the simplex orientation, so the int type you use should definitely be signed. The fact that it can allow you to pass in an UInt64 is a bug.

Maybe that's why the code hangs as negating an UInt gives you something invalid. It's hard to say for sure though as I don't know what fct looks like.

Try it with Int128 first and let me know if it works!

[lampretl]

I can confirm that Ripserer.Custom{T} and Ripserer.ripserer work for T=Int128 (needs 12min in total for RP8). But already for T=Int120 from BitIntegers.jl it fails, even though the largest boundary matrix is cca 15M x 15M. So in this respect, this issue is solved. I would recommend @mtsch to improve the functionality, so that Ripserer.Custom would automatically detect what the smallest sufficient type is, and use that. And only allow types <:Signed, so the program won't hang indefinitely.

For anyone else interested, I did some digging and realized that Ripser uses a combinatorial number system, which in the Ripser article is paragraph Indexing simplices in the combinatorial number system below Proposition 2.1 on page 398. So Ripserer.jl is using something like:

function index_of_simplex(σ::AbstractVector{ti}, ::th, v0::Int=1) ::th  where {ti<:Integer, th<:Integer}
    return sum(binomial(th(v-v0),th(i)) for (i,v) ∈ enumerate(σ)) 
end
julia> using IterTools
julia> v0=1; n=5; k=3;  println("max =",binomial(n,k)-1); 
julia> for σ ∈ sort(collect(subsets(v0:v0+n-1,k)), by=reverse)  
          println("σ = $σ   index = $(index_of_simplex(σ,Int(0),v0))") end
max = 9
σ = [1, 2, 3]   index = 0
σ = [1, 2, 4]   index = 1
σ = [1, 3, 4]   index = 2
σ = [2, 3, 4]   index = 3
σ = [1, 2, 5]   index = 4
σ = [1, 3, 5]   index = 5
σ = [2, 3, 5]   index = 6
σ = [1, 4, 5]   index = 7
σ = [2, 4, 5]   index = 8
σ = [3, 4, 5]   index = 9

If the smallest vertex would be 0 instead of 1, like in C++, then we'd use v0=0 and we'd get the same indices. This is here https://github.com/mtsch/Ripserer.jl/blob/9583a3b7e64565f11f869155b346b675eb508740/src/base/abstractsimplex.jl#L138using, though I don't understand why it's done recursively, when a simple sum suffices.

What is still unclear to me: if vertices are between 1:n, and we have k-simplices, then the indices range between 0:binomial(n,k+1)-1. In the case of RP8, the indices range between 0:binomial(511,9)-1 = 6098989894499557054, but ``typemax(Int64) = 9223372036854775807` is larger. So why the overflow?

[mtsch]

I've added a quick fix that makes Custom only work with Signed. Picking the int size automatically is a good idea, but I'll have to think if there's a robust way to determine it. As you've noticed, the integer size needs to be quite a bit larger than the maximal index. Computing binomials and performing the binary search to find simplex vertices requires more bits than the index itself.