mtsch / ripserer.jl Goto Github PK
View Code? Open in Web Editor NEWFlexible and efficient persistent homology computation.
Home Page: https://mtsch.github.io/Ripserer.jl/dev/
License: MIT License
Flexible and efficient persistent homology computation.
Home Page: https://mtsch.github.io/Ripserer.jl/dev/
License: MIT License
Hi, thank you for this awesome package. Right now if I pass a SparseCSCMatrix object to Ripserer directly it gives an error saying that some edge lengths are 0. However if I first try to specify a dense matrix (filling the 0s with a large number) and then supply it with a threshold, I run into memory issues constructing the dense matrix.
Would it be possible to allow SparseCSCMatrix inputs directly, where off-diagonal zeros are considered as above the threshold by default?
Is it possible to implement reconstruct_cycle
for more than 2-dimension? Let's say in 3D case, the function will compute all the 2-simplexes that "minimally" surround a 3-dimensional hole.
Also the cycles.jl
code has no comment, it's very difficult to understand what is going on, specially in the _find_cycle()
function. For example,
best_path = edgetype(g)[] # what does this line do?
best_sx = first(g.removed) # what does this line do?
Is it possible to add some algorithmic descriptions for the functions in the cycles.jl
file? It will be really helpful.
Thanks!! I really appreciate your work.
Hi, Tom from Codecov here.
We noticed that you are using Codecov with fairly high frequency, and we’re so excited to see that! However, because you are not using our app, you may have experienced issues with uploading reports or viewing coverage information. This is due to rate-limiting issues from GitHub.
In order to prevent any future outages, we ask that you move over to our GitHub app integration.
The process is extremely simple and shouldn’t require more than a few clicks, and you should not expect any downtime. By moving to our app, you will no longer need an admin or separate account to manage the relationship with GitHub as the team bot.
Let me know if you have any questions, or if I can help at all with this process.
I'm not sure how hard this would be or if it is even a new idea. One of the things I'm interested in is doing persistent homology on a Cubical complex where we identify the top and bottom of the matrix and the right and left sides, essentially wrapping around like a torus. I haven't dug into the code to see what the right approach would be. Maybe making a TorusComplex object, and that somehow tells ripserer what to identify? Or should there just be a flag to ripserer when you submit a Cubical object? Thanks!
Thank you @mtsch for this great package. Could you please add a section with that comparison you shared on Discourse?
Hi,
Thanks for this awesome package.
I'm trying to use it on some fairly large set of data points (~10k data points of 5-10 dimensions). This is already after substantial dimensionality reduction and downsampling, but it still computationally very demanding to run the algorithm on.
Would you have any advice on things to try to speed it up? Would downsampling (e.g. dropping X% of the data points) would be a valid approach?
Thank you,
Federico
I get the following error when trying to use the reconstruct cycles function. I have tried manually upgrading the Graphs library to see if that helps but it still gives the same error. I can provide the underlying data if needed. Thank you for including this functionality - I don't see it available in any other packages!
BoundsError: attempt to access Tuple{Int64} at index [2]
Stacktrace:
[1] getindex(t::Tuple, i::Int64)
@ Base ./tuple.jl:29
[2] getindex
@ ~/.julia/packages/StaticArrays/J9itA/src/SArray.jl:62 [inlined]
[3] _next_common
@ ~/.julia/packages/Ripserer/RPZgP/src/base/abstractsimplex.jl:232 [inlined]
[4] iterate
@ ~/.julia/packages/Ripserer/RPZgP/src/base/abstractsimplex.jl:256 [inlined]
[5] iterate
@ ~/.julia/packages/Ripserer/RPZgP/src/base/abstractsimplex.jl:252 [inlined]
[6] outneighbors(g::OneSkeleton{Simplex{1, Float64, Int64}, Rips{Int64, Float64, SparseMatrixCSC{Float64, Int64}}, Simplex{1, Float64, Int64}, SparseMatrixCSC{Float64, Int64}}, u::Int64)
@ Main ./In[108]:63
[7] a_star_impl!(g::OneSkeleton{Simplex{1, Float64, Int64}, Rips{Int64, Float64, SparseMatrixCSC{Float64, Int64}}, Simplex{1, Float64, Int64}, SparseMatrixCSC{Float64, Int64}}, goal::Int64, open_set::DataStructures.PriorityQueue{Int64, Float64, Base.Order.ForwardOrdering}, closed_set::Vector{Bool}, g_score::Vector{Float64}, came_from::Vector{Int64}, distmx::SparseMatrixCSC{Float64, Int64}, heuristic::var"#38#39"{Int64, SparseMatrixCSC{Float64, Int64}}, edgetype_to_return::Type{Graphs.SimpleGraphs.SimpleEdge{Int64}})
@ Graphs ~/.julia/packages/Graphs/7SMZs/src/shortestpaths/astar.jl:43
[8] a_star(g::OneSkeleton{Simplex{1, Float64, Int64}, Rips{Int64, Float64, SparseMatrixCSC{Float64, Int64}}, Simplex{1, Float64, Int64}, SparseMatrixCSC{Float64, Int64}}, s::Int64, t::Int64, distmx::SparseMatrixCSC{Float64, Int64}, heuristic::Function, edgetype_to_return::Type{Graphs.SimpleGraphs.SimpleEdge{Int64}})
@ Graphs ~/.julia/packages/Graphs/7SMZs/src/shortestpaths/astar.jl:94
[9] a_star
@ ~/.julia/packages/Graphs/7SMZs/src/shortestpaths/astar.jl:81 [inlined]
[10] _find_cycle(g::OneSkeleton{Simplex{1, Float64, Int64}, Rips{Int64, Float64, SparseMatrixCSC{Float64, Int64}}, Simplex{1, Float64, Int64}, SparseMatrixCSC{Float64, Int64}}, dists::SparseMatrixCSC{Float64, Int64})
@ Main ./In[108]:115
[11] reconstruct_cycle_2(filtration::Rips{Int64, Float64, SparseMatrixCSC{Float64, Int64}}, interval::PersistenceDiagrams.PersistenceInterval, r::Simplex{1, Float64, Int64}; distances::SparseMatrixCSC{Float64, Int64})
@ Main ./In[108]:179
[12] reconstruct_cycle_2(filtration::Rips{Int64, Float64, SparseMatrixCSC{Float64, Int64}}, interval::PersistenceDiagrams.PersistenceInterval, r::Simplex{1, Float64, Int64})
@ Main ./In[108]:158
[13] reconstruct_cycle_2(filtration::Rips{Int64, Float64, SparseMatrixCSC{Float64, Int64}}, interval::PersistenceDiagrams.PersistenceInterval)
@ Main ./In[108]:158
[14] top-level scope
@ In[109]:2
A promote_rule is not supposed to non-terminating, as that violates the specification of promote_type.
Look at the code here:
Ripserer.jl/src/base/primefield.jl
Line 84 in 62a369c
There should be an additional rule to resolve the ambiguity there. That looks like:
Base.promote_rule(::Type{Mod{M}}, ::Type{<:Mod}) where {M} = Union{}
Additionally, I noticed that you appear to have assumed that all Integers
are an Int
, which is not always true. You can make your current rule more accurate by first checking it the Integer is directly promotable to Int:
Base.promote_rule(::Type{Mod{M}}, ::Type{N}) where {M, N<:Integer} = Base.promote_type(N, Int) === Int ? Mod{M} : Union{}
This issue is used to trigger TagBot; feel free to unsubscribe.
If you haven't already, you should update your TagBot.yml
to include issue comment triggers.
Please see this post on Discourse for instructions and more details.
Hi!
I've been wondering if it's possible to add multithreading or GPU support to the computation of persistence intervals in a straightforward way. This could be particularly desirable for very large filtrations. So far, the only package I've found that supports this is Giotto-tda. Do you think this could be of interest to a wider audience?
I would be happy to help if you believe it's worthwhile.
Hello. I am still relatively new to computing with simplicial complexes, but I use Julia all the time. I was wondering if Ripserer is able to compute the discrete morse theory--or essentially the collapses of a simplicial complex? I believe Eirene does this internally, for its computation of persistent homology, but not sure about Ripserer. If you or anyone happens to know of a julia package that computes the collapses for a complex, please let me know.
Hi Matija,
I've been trying to see how to add the feature of keeping track which component merged with whom. I have some preliminary (messy code) that does this. I guess this would only work for dim=0 homology since only there these merges make sense.
1.) Do you perhaps have any suggestion on how/what to improve on this?
2.) Would it be worth to add support in Ripserer for this, as well as constructing/plotting the corresponding merge tree?
I haven't found anything that supports this out there so perhaps could be useful
The code below is a simple modification of the zeroth_intervals
function
function zeroth_intervals_with_merges(filtration, cutoff::T = T(0); verbose=true) where T <: AbstractFloat
F = Mod
reps = false
V = simplex_type(filtration, 0)
CE = chain_element_type(V, F)
dset = DisjointSetsWithBirth(vertices(filtration), births(filtration))
intervals = PersistenceInterval[]
to_skip = simplex_type(filtration, 1)[]
to_reduce = simplex_type(filtration, 1)[]
simplices = sort!(edges(filtration))
merged_vertices = Dict()
persistence_merged_vertices = Dict()
if verbose
progbar = Progress(
length(simplices) + nv(filtration); desc="Computing 0d intervals... "
)
end
for edge in simplices
u, v = vertices(edge)
i = find_root!(dset, u)
j = find_root!(dset, v)
if i ≠ j
# According to the elder rule, the vertex with the higher birth will die first.
last_vertex = birth(dset, i) > birth(dset, j) ? i : j
int = interval(dset, filtration, last_vertex, edge, cutoff, reps)
btime_last, bvertex_last = birth(dset, last_vertex)
if !isnothing(int)
push!(intervals, int)
if last_vertex == i
btime_first, bvertex_first = birth(dset, j)
else
btime_first, bvertex_first = birth(dset, i)
end
merged_vertices[simplex(filtration, Val(0), (bvertex_last,))] = simplex(filtration, Val(0), (bvertex_first,))
persistence_merged_vertices[simplex(filtration, Val(0), (bvertex_last,))] = persistence(int)
end
union!(dset, i, j)
push!(to_skip, edge)
else
push!(to_reduce, edge)
end
verbose && next!(progbar; showvalues=((:n_intervals, length(intervals)),))
end
for v in vertices(filtration)
if find_root!(dset, v) == v && !isnothing(simplex(filtration, Val(0), (v,)))
int = interval(dset, filtration, v, nothing, cutoff, reps)
push!(intervals, int)
end
verbose && next!(progbar; showvalues=((:n_intervals, length(intervals)),))
end
reverse!(to_reduce)
thresh = T(threshold(filtration))
diagram = PersistenceDiagram(
sort!(intervals; by=persistence);
threshold=thresh,
dim=0,
field=F,
filtration=filtration,
)
return merged_vertices, Ripserer.postprocess_diagram(filtration, diagram), persistence_merged_vertices
end
I can also share the code to construct and plot the corresponding merge tree if you think it's worthwhile or necessary
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.