[Port] Have added Wilson's algorithm and loop erased random walks #24
Conversation
Codecov Report
@@ Coverage Diff @@
## master #24 +/- ##
==========================================
- Coverage 99.45% 99.44% -0.02%
==========================================
Files 106 107 +1
Lines 5554 5602 +48
==========================================
+ Hits 5524 5571 +47
- Misses 30 31 +1 |
|
@nassarhuda, are you familiar with Wilson's algorithm? Could you take a look here? |
| function wilson_rst end | ||
| @traitfn function wilson_rst(g::AG::(!IsDirected), | ||
| distmx::AbstractMatrix{T}=weights(g); | ||
| seed::Int=-1, |
There was a problem hiding this comment.
Do we really need a seed, when we can pass in a random generator?
| tree = wilson_rst(g) | ||
| probability_of_tree = prod([weights(g)[src(e), dst(e)] for e in tree]) | ||
| probability_of_tree /= det(laplacian_matrix(g)[2:nv(g), 2:nv(g)])) |
There was a problem hiding this comment.
Can we use a julia or jldoctest block here?
| start1 = rand(rng, 1:nv(g)) | ||
| start2 = rand(rng, 1:nv(g)) |
There was a problem hiding this comment.
| start1 = rand(rng, 1:nv(g)) | |
| start2 = rand(rng, 1:nv(g)) | |
| start1 = rand(rng, vertices(g)) | |
| start2 = rand(rng, vertices(g)) |
might make more sense, although at least for SimpleGraph, this will not be of type Int but of type eltype(g)
| seed::Int=-1, | ||
| rng::AbstractRNG=GLOBAL_RNG |
There was a problem hiding this comment.
Seem here, do we need seed?
| f::Union{Set{V},Vector{V}}=Set{Integer}(), | ||
| seed::Int=-1, | ||
| rng::AbstractRNG=GLOBAL_RNG | ||
| )::Vector{Int} where {T <: Real, U, AG <: AbstractGraph{U}, V <: Integer} |
There was a problem hiding this comment.
If I remember ::SomeType after a function does actually force convert to be used on the result of that function.
|
@etiennedeg @gjh6 I think this is something useful to have. Given that this PR has been lying dormant for quite a while, there might be some things that you would do differently? Then maybe it is better, if you would refactor this before me or someone else continues with reviewing. |
| rng = getRNG(seed) | ||
| end | ||
|
|
||
| visited = Vector{Integer}(undef, 1) |
There was a problem hiding this comment.
Vector should be typed with eltype(g) (passed by dynamic dispatch)
| i += 1 | ||
| end | ||
|
|
||
| length(f) == 0 || visited_view[cur_pos] in f || throw(ErrorException("termiating set was not reached")) |
There was a problem hiding this comment.
| length(f) == 0 || visited_view[cur_pos] in f || throw(ErrorException("termiating set was not reached")) | |
| length(f) == 0 || visited_view[cur_pos] in f || throw(ErrorException("terminating set was not reached")) |
| break | ||
| end | ||
| nbrs = neighbors(g, cur) | ||
| length(nbrs) == 0 && break |
There was a problem hiding this comment.
Can we do this check only for the first vertex? Graph is undirected so if we move to an adjacent neighbor, this vertex must have at least one neighbor ?
| end | ||
| nbrs = neighbors(g, cur) | ||
| length(nbrs) == 0 && break | ||
| wght = [distmx[cur, n] for n in nbrs] |
There was a problem hiding this comment.
| wght = [distmx[cur, n] for n in nbrs] | |
| wght = @view distmx[cur, nbrs] |
Maybe this is better?
| if v in visited_view | ||
| cur_pos = indexin(v, visited_view)[1] | ||
| visited_view = view(visited, 1:cur_pos) |
There was a problem hiding this comment.
| if v in visited_view | |
| cur_pos = indexin(v, visited_view)[1] | |
| visited_view = view(visited, 1:cur_pos) | |
| new_cur_pos = findfirst(v, visited_view) | |
| if !isnothing(new_cur_pos) | |
| cur_pos = new_cur_pos | |
| visited_view = view(visited, 1:cur_pos) | |
| end |
| end | ||
|
|
||
| visited = Vector{Integer}(undef, 1) | ||
| visited_view = view(visited, 1:1) |
There was a problem hiding this comment.
Is it necessary to carry a view? The only place where it is used is line 164-165 (and only in one line with my proposed change). The indexing of visited_view is the same as visited because it always start at the first index of visited, so in line 152, 156 and 178, we can replace visited_view by visited.
| end | ||
|
|
||
| visited_vertices = Set(walk) | ||
| unvisited_vertices = setdiff(Set([i for i = 1:nv(g)]), visited_vertices) |
There was a problem hiding this comment.
It is more efficient to just flag visited vertices (with a linear pass over walk) than to do complex set operations (which I think are not even linear...). We just maintain a BitVector is_visitedof size nv(g).
In loop_erased_randomwalk, e in f becomes is_visited[e], which is also more efficient.
|
|
||
| walk = loop_erased_randomwalk(g, start1, distmx=distmx, f=[start2], rng=rng) | ||
|
|
||
| tree = SimpleGraph(nv(g)) |
There was a problem hiding this comment.
tree is only used to together with add_edge!, which is a costly operation. It is more efficient to directly store the unique edges.
| end | ||
|
|
||
| visited_vertices = Set(walk) | ||
| unvisited_vertices = setdiff(Set([i for i = 1:nv(g)]), visited_vertices) |
There was a problem hiding this comment.
| unvisited_vertices = setdiff(Set([i for i = 1:nv(g)]), visited_vertices) | |
| unvisited_vertices = setdiff(Set(vertices(g)), visited_vertices) |
this is a port of #1576