faster implementation of erdos_renyi

[etiennedeg]

Master :

julia> @benchmark erdos_renyi(500, 0.5)
BenchmarkTools.Trial: 115 samples with 1 evaluation.
 Range (min … max):  31.674 ms … 102.297 ms  ┊ GC (min … max): 0.00% … 14.32%
 Time  (median):     39.373 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   43.619 ms ±  13.058 ms  ┊ GC (mean ± σ):  1.75% ±  5.16%

     ▇▂█▂▆▁▂ ▄ 
  ▄█████████▃█▆▆█▁▆▄▃▆▃▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▃▁▁▃▁▁▁▁▁▁▄▁▁▁▄▁▃▁▁▁▁▃ ▃
  31.7 ms         Histogram: frequency by time         90.9 ms <

 Memory estimate: 3.69 MiB, allocs estimate: 2502.

julia> @benchmark erdos_renyi(5000, 0.5)
BenchmarkTools.Trial: 1 sample with 1 evaluation.
 Single result which took 17.436 s (1.51% GC) to evaluate,
 with a memory estimate of 263.33 MiB, over 35003 allocations.

julia> @benchmark erdos_renyi(500, is_directed = true, 0.5)
BenchmarkTools.Trial: 82 samples with 1 evaluation.
 Range (min … max):  58.535 ms … 67.722 ms  ┊ GC (min … max): 0.00% … 5.64%
 Time  (median):     61.109 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   61.466 ms ±  1.972 ms  ┊ GC (mean ± σ):  0.98% ± 2.18%

         █▂▂  ▅ ▅  █▂▂    ▂      ▂
  ████▁▁███████▅█▅▁████▁██████▅▅▅█▅▁▅▁█▁▁▅█▅▅▅▁▅▅▁▁▁▁▁▅▁▁▁▁▁▅ ▁
  58.5 ms         Histogram: frequency by time        66.9 ms <

 Memory estimate: 7.38 MiB, allocs estimate: 5003.

julia> @benchmark erdos_renyi(5000, is_directed = true, 0.5)
BenchmarkTools.Trial: 1 sample with 1 evaluation.
 Single result which took 37.086 s (0.62% GC) to evaluate,
 with a memory estimate of 526.66 MiB, over 70005 allocations.

This commit :

julia> @benchmark erdos_renyi(500, 0.5)
BenchmarkTools.Trial: 847 samples with 1 evaluation.
 Range (min … max):  3.717 ms … 12.779 ms  ┊ GC (min … max): 0.00% … 55.33%
 Time  (median):     5.758 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   5.883 ms ±  1.727 ms  ┊ GC (mean ± σ):  8.63% ± 15.70%

  ▂             █▆▁▁        
  █▆▄▃▃▃▃▃▃▃▄▅▆▆████▄▃▃▃▁▁▂▁▂▂▂▁▂▂▂▂▂▁▂▁▁▁▁▂▁▃▃▂▂▃▂▃▂▃▂▃▂▂▂▂ ▃
  3.72 ms        Histogram: frequency by time        11.7 ms <

 Memory estimate: 3.69 MiB, allocs estimate: 2502.

julia> @benchmark erdos_renyi(5000, 0.5)
BenchmarkTools.Trial: 7 samples with 1 evaluation.
 Range (min … max):  519.211 ms …    1.668 s  ┊ GC (min … max):  6.73% … 67.66%
 Time  (median):     586.219 ms               ┊ GC (median):    11.76%
 Time  (mean ± σ):   752.332 ms ± 411.881 ms  ┊ GC (mean ± σ):  29.61% ± 23.51%

  █▁ ▁    ▁  ▁                                                ▁  
  ██▁█▁▁▁▁█▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  519 ms           Histogram: frequency by time          1.67 s <

 Memory estimate: 263.33 MiB, allocs estimate: 35003.

julia> @benchmark erdos_renyi(500, is_directed = true, 0.5)
BenchmarkTools.Trial: 576 samples with 1 evaluation.
 Range (min … max):  7.551 ms … 13.974 ms  ┊ GC (min … max): 0.00% … 23.43%
 Time  (median):     8.205 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   8.665 ms ±  1.101 ms  ┊ GC (mean ± σ):  5.53% ±  9.83%

  ▃▃▂▂▂▂ ▅█▇▄▃▁▂    ▁                    ▁  ▁▁▁▁ ▁
  ███████████████▇▇▇█▅▆▄▁▄▆▅▆▁▁▁▇▁▁▁▁▄▄▁▇███████▄█▆▅▇▇▆▄▅▄▁▄ ▇
  7.55 ms      Histogram: log(frequency) by time     11.7 ms <

 Memory estimate: 7.38 MiB, allocs estimate: 5003.

julia> @benchmark erdos_renyi(5000, is_directed = true, 0.5)
BenchmarkTools.Trial: 4 samples with 1 evaluation.
 Range (min … max):  1.143 s …    1.453 s  ┊ GC (min … max):  0.00% … 22.82%
 Time  (median):     1.376 s               ┊ GC (median):    21.56%
 Time  (mean ± σ):   1.337 s ± 135.752 ms  ┊ GC (mean ± σ):  17.84% ± 11.62%

  █                                     █         █        █  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁█ ▁
  1.14 s         Histogram: frequency by time         1.45 s <

 Memory estimate: 526.66 MiB, allocs estimate: 70005.

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 99.23372% with 2 lines in your changes are missing coverage. Please review.

Project coverage is 97.36%. Comparing base (9c8efd9) to head (6e22b60).

Files	Patch %	Lines
src/SimpleGraphs/generators/randgraphs.jl	99.22%	2 Missing :warning:

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #150      +/-   ##
==========================================
+ Coverage   97.31%   97.36%   +0.05%     
==========================================
  Files         120      120              
  Lines        6953     7175     +222     
==========================================
+ Hits         6766     6986     +220     
- Misses        187      189       +2     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[simonschoelly]

The problem this benchmark is, that you only benchmark if for p=0.5, which is rather unrealistic for larger graphs.

The implementation we currently have in Graph.jl, can certainly speed up, and I think in some edges cases, there might even be error. But there are better algorithms than in this PR. I implemented something 4 years ago, but then got sidetracked, so I am not sure if it is correct, but I will attach that code, so you can take it or maybe use it for creating a better implementation.

My memory is very shaky here, so there might be some errors:

• For the (n, p) graph (where each edge is there with probability p), one iterates over all potential edges. If an edge (u, v) is the last included one, one uses the geometric distribution to directly jump to the next edge '(u', v')`.
• For the (n, m) graph, (where we select one graph with exactly m edges uniformely), one uses something called Vitters Method D, described here: https://www.ittc.ku.edu/~jsv/Papers/Vit87.RandomSampling.pdf

And here is the code I once wrote, no guarantee on correctness:

using LightGraphs
using Random
using Base: OneTo


function erdos_renyi_directed(nv::T, m::Integer; self_loops=false) where {T <: Integer}
    nvi = Int(nv)

    max_possible_edges = self_loops ? (nvi * nvi) : (nvi * (nvi - 1))
    @assert 0 <= m <= max_possible_edges

    fadjlist = Vector{Vector{T}}(undef, nvi)
    taken = 0
    remaining = max_possible_edges
    u = 1
    v = 1
    buffer = Vector{T}(undef, nv)
    buffer_len = 0

    @inbounds while u <= nv
        v += ifelse(!self_loops && u == v, 1, 0)
        
        if v > nv
            fadjlist[u] = buffer[OneTo(buffer_len)]
            u += 1
            v = 1
            buffer_len = 0
        end

        if rand() * remaining < (m - taken)
            buffer_len += 1
            buffer[buffer_len] = v
            taken += 1
        end

        v += 1
        remaining -= 1
    end

    # reuse the memory allocated for  buffer
    indeg = buffer
    fill!(indeg, zero(T))
    @inbounds for u in OneTo(nv)
        for v in fadjlist[u]
            indeg[v] += 1
        end
    end

    badjlist = Vector{Vector{T}}(undef, nvi)
    @inbounds for v in OneTo(nv)
        badjlist[v] = Vector{T}(undef, indeg[v])
    end 

    @inbounds for u in OneTo(nv)
        for v in fadjlist[u]
            badjlist[v][end - indeg[v] + 1] = u
            indeg[v] -= 1
        end
    end

    return SimpleDiGraph(m, fadjlist, badjlist)
end



function erdos_renyi_directed(nv::T, p::Union{Rational, AbstractFloat, AbstractIrrational}; self_loops=false) where {T <: Integer}
    @assert 0 <= p <= 1 
    
    fadjlist = Vector{Vector{T}}(undef, nv)
    
    ne = 0
    buffer = Vector{T}(undef, nv)
    sample_from = self_loops ? collect(one(T):nv) : collect(T(2):nv)
    @inbounds for u in OneTo(nv)
        randsubseq!(buffer, sample_from, p)
        fadjlist[u] = copy(buffer)
        ne += length(buffer)
        
        if !self_loops
            sample_from[u] = u
        end
    end

    # reuse the memory allocated for  buffer
    # need to resize! because randsubseq! could have shrunken the buffer length
    indeg = resize!(buffer, nv)
    fill!(indeg, zero(T))
    @inbounds for u in OneTo(nv)
        for v in fadjlist[u]
            indeg[v] += 1
        end
    end

    badjlist = Vector{Vector{T}}(undef, nv)
    @inbounds for v in OneTo(nv)
        badjlist[v] = Vector{T}(undef, indeg[v])
    end 

    @inbounds for u in OneTo(nv)
        for v in fadjlist[u]
            badjlist[v][end - indeg[v] + 1] = u
            indeg[v] -= 1
        end
    end
    return SimpleDiGraph(ne, fadjlist, badjlist)
end

function erdos_renyi_undirected(nv::T, m::Integer; self_loops=false) where {T <: Integer}
    nvi = Int(nv)

    max_possible_edges = self_loops ? (nvi * (nvi + 1)) ÷ 2 : (nvi * (nvi - 1)) ÷ 2
    @assert 0 <= m <= max_possible_edges

    fadjlist = Vector{Vector{T}}(undef, nvi)
    taken = 0
    remaining = max_possible_edges
    u = 1
    v = self_loops ? one(T) : T(2)
    buffer = Vector{T}(undef, nv)
    buffer_len = 0
    indeg = zeros(T, nv)
    @inbounds while u <= nv
        if v > nv
            indeg_u = indeg[u]
            list_u = Vector{T}(undef, indeg_u + buffer_len)

            for i in OneTo(buffer_len)
                w = buffer[i]
                list_u[i + indeg_u] = w
                indeg[w] += 1
            end

            fadjlist[u] = list_u

            u += 1
            v = ifelse(self_loops, u, u + 1)
            buffer_len = 0
        end

        if rand() * remaining < (m - taken)
            buffer_len += 1
            buffer[buffer_len] = v
            taken += 1
        end

        v += 1
        remaining -= 1
    end

    insert_at = resize!(buffer, nv)
    fill!(insert_at, one(T))

    @inbounds for u in OneTo(nv)
        list_u = fadjlist[u]
        for i in (indeg[u] + 1):length(list_u)
            v = list_u[i]
            fadjlist[v][insert_at[v]] = u
            insert_at[v] += 1
        end
    end

    return SimpleGraph(m, fadjlist)
end




@inline function _sample_with_replacement!(buffer, sample_from, p, buffer_offset=0)
    buffer_len = 0
    if p >= 0.5
        @inbounds for v in sample_from
            if rand() < p
                buffer_len += 1
                buffer[buffer_len + buffer_offset] = v
            end
        end
    else
        L = -1 / log1p(-p)
        sample_from_len = length(sample_from)
        i = 0
        @inbounds while true
            i += floor(Int, randexp() * L) + 1
            i > sample_from_len && return buffer_len
            buffer_len += 1
            buffer[buffer_len + buffer_offset] = sample_from[i]
        end
    end
    return buffer_len
end


function erdos_renyi_undirected(nv::T, p::Union{Rational, AbstractFloat, AbstractIrrational}; self_loops=false) where {T <: Integer}
    @assert 0 <= p <= 1 

    fadjlist = Vector{Vector{T}}(undef, nv)

    ne = 0
    buffer = Vector{T}(undef, nv)
    indeg = zeros(T, nv)

    @inbounds for u in OneTo(nv)
        sample_from = self_loops ? (u:nv) : (u+one(T):nv)

        buffer_len = _sample_with_replacement!(buffer, sample_from, p)

        indeg_u = indeg[u]
        list_u = Vector{T}(undef, indeg_u + buffer_len)

        for i in OneTo(buffer_len)
            v = buffer[i] 
            list_u[i + indeg_u] = v
            indeg[v] += 1
        end

        fadjlist[u] = list_u
        
        ne += buffer_len
    end

    insert_at = resize!(buffer, nv)
    fill!(insert_at, one(T))

    @inbounds for u in OneTo(nv)
        list_u = fadjlist[u]
        for i in (indeg[u] + 1):length(list_u)
            v = list_u[i]
            fadjlist[v][insert_at[v]] = u
            insert_at[v] += 1
        end
    end

    return SimpleGraph(ne, fadjlist)
end

[etiennedeg]

Oh ok, I see the point of using the binomial distribution. The biggest problem here is the use of add_edge!, which is very inefficient. I will check your code to see if it gives improvement when I will have time.

[simonschoelly]

I mean, if we have some benchmarks, that show that this implementation is also faster for very small and very large p such as p = k/n and p = (n-k)/n (for some small integer k) then we could already merge this.

[abraunst]

The following code implements the good strategy above for (n,p) graphs. If there's interest I can prepare a PR.

function trianglemap(r)
    j = floor(Int, (1+sqrt(8r-7))/2) + 1
    i = r - binomial(j-1,2)
    i, j
end

function nondiagmap(r,n)
    j = div(r-1, n-1)
    i = r - j*(n-1)
    j += 1
    i += (i >= j)
    i, j
end

function erdos_renyi(
    n::Integer,
    p::Real;
    is_directed=false,
    rng::Union{Nothing,AbstractRNG}=nothing,
    seed::Union{Nothing,Integer}=nothing,
)
    rng = rng_from_rng_or_seed(rng, seed)
    m = is_directed ? n*(n-1) : binomial(n,2)
    R = randsubseq(rng, 1:m, min(1.0, p))
    g = if is_directed
        SimpleDiGraphFromIterator(Edge(nondiagmap(r,n)...) for r in R)
    else
        SimpleGraphFromIterator(Edge(trianglemap(r)...) for r in R)
    end
    add_vertices!(g, n - nv(g))
    return g
end

To compare I've used

for N in (500, 5000), is_directed in (true, false), p in (10/N, 0.5, 1-10/N)
    @btime erdos_renyi($N, $p; is_directed=$is_directed)
end

That gives

634.683 μs (2601 allocations: 387.77 KiB)
  15.538 ms (5023 allocations: 8.35 MiB)
  20.433 ms (5023 allocations: 9.27 MiB)
  407.949 μs (1285 allocations: 187.91 KiB)
  9.689 ms (2517 allocations: 4.18 MiB)
  14.804 ms (2517 allocations: 4.64 MiB)
  7.411 ms (26507 allocations: 4.08 MiB)
  1.679 s (70029 allocations: 622.30 MiB)
  2.491 s (80029 allocations: 1.35 GiB)
  4.633 ms (13244 allocations: 2.04 MiB)
  1.066 s (35020 allocations: 311.18 MiB)
  1.666 s (40020 allocations: 691.04 MiB)

which seems uniformly better than master:

1.230 ms (2637 allocations: 356.81 KiB)
  68.167 ms (5003 allocations: 7.38 MiB)
  47.454 ms (7675 allocations: 7.74 MiB)
  583.637 μs (1307 allocations: 174.48 KiB)
  30.522 ms (2502 allocations: 3.69 MiB)
  22.270 ms (3823 allocations: 3.86 MiB)
  15.997 ms (26913 allocations: 3.66 MiB)
  42.191 s (70005 allocations: 526.66 MiB)
  10.789 s (107180 allocations: 1.17 GiB)
  7.240 ms (13430 allocations: 1.82 MiB)
  18.297 s (35003 allocations: 263.33 MiB)
  4.420 s (53510 allocations: 597.51 MiB)

[abraunst]

I went ahead and opened https://github.com/JuliaGraphs/Graphs.jl/pull/212

[etiennedeg]

I added the implementation of @simonschoelly for G(n, p) graphs, and I implemented Vitter's algorithm for G(n, m) graphs. As was initiated in Simon's implementation, we can now allow self loops generation by passing the self_loops argument.

It misses some tests. I think a necessary test is to ensure the generated graphs have a coherent structure. For test about assessing the distribution, I still don't know the best way to test for it.

[abraunst]

I added the implementation of @simonschoelly for G(n, p) graphs, and I implemented Vitter's algorithm for G(n, m) graphs. As was initiated in Simon's implementation, we can now allow self loops generation by passing the self_loops argument.

It misses some tests. I think a necessary test is to ensure the generated graphs have a coherent structure. For test about assessing the distribution, I still don't know the best way to test for it.

Great! Note that there are implementations of the various algorithms by Vitter, including this one, in StatsBase. Best.

[etiennedeg]

Thanks, I will take a look. Maybe better to not add a dependency here, the algorithm is sufficiently simple to be contained in Graphs.jl

[etiennedeg]

@gdalle This should be ready for review.

I'm unable to reach this last branch for codecov

Edit: This is triggered non-deterministically and very rarely, I don't think it makes sense to try to cover it in the tests