From 64e1aba1a08cf200b62d3fb7c5e0b4f0b4c2ed3d Mon Sep 17 00:00:00 2001 From: etienneINSA Date: Fri, 24 Jun 2022 17:22:00 +0200 Subject: [PATCH 1/7] faster implementation of erdos_renyi --- src/SimpleGraphs/generators/randgraphs.jl | 29 +++++++++++++++++------ 1 file changed, 22 insertions(+), 7 deletions(-) diff --git a/src/SimpleGraphs/generators/randgraphs.jl b/src/SimpleGraphs/generators/randgraphs.jl index 9e879e032..e1f2fec1b 100644 --- a/src/SimpleGraphs/generators/randgraphs.jl +++ b/src/SimpleGraphs/generators/randgraphs.jl @@ -160,13 +160,28 @@ function erdos_renyi( seed::Union{Nothing,Integer}=nothing, ) p >= 1 && return is_directed ? complete_digraph(n) : complete_graph(n) - m = is_directed ? n * (n - 1) : div(n * (n - 1), 2) - ne = randbn(m, p; rng=rng, seed=seed) - return if is_directed - SimpleDiGraph(n, ne; rng=rng, seed=seed) - else - SimpleGraph(n, ne; rng=rng, seed=seed) + rng = getRNG(seed) + m = 0 + fadj = [Int[] for u in 1:n] + if is_directed + badj = [Int[] for u in 1:n] + end + for u in 1:n + start = is_directed ? 1 : u+1 + for v in start:n + v == u && continue + if rand(rng) < p + m += 1 + push!(fadj[u], v) + if is_directed + push!(badj[v], u) + else + push!(fadj[v], u) + end + end + end end + return is_directed ? SimpleDiGraph(m, fadj, badj) : SimpleGraph(m, fadj) end """ @@ -289,7 +304,7 @@ end watts_strogatz(n, k, β) Return a [Watts-Strogatz](https://en.wikipedia.org/wiki/Watts_and_Strogatz_model) -small world random graph with `n` vertices, each with expected degree `k` +small world random graph with `n` vertices, each with expected degree `k` (or `k - 1` if `k` is odd). Edges are randomized per the model based on probability `β`. The algorithm proceeds as follows. First, a perfect 1-lattice is constructed, From 6caf058286825129c9cafa2cd9f92d3b84cca8fb Mon Sep 17 00:00:00 2001 From: etienneINSA Date: Mon, 12 Feb 2024 17:42:49 +0100 Subject: [PATCH 2/7] erdos-renyi with fast sampling --- src/SimpleGraphs/generators/randgraphs.jl | 1870 ++++++++++++--------- 1 file changed, 1082 insertions(+), 788 deletions(-) diff --git a/src/SimpleGraphs/generators/randgraphs.jl b/src/SimpleGraphs/generators/randgraphs.jl index e1f2fec1b..d008c7716 100644 --- a/src/SimpleGraphs/generators/randgraphs.jl +++ b/src/SimpleGraphs/generators/randgraphs.jl @@ -1,9 +1,116 @@ -using Random: randperm, shuffle! +using Random: randperm, shuffle!, randexp using Statistics: mean using Graphs: sample! +using Base: OneTo + +""" + _d_sample!(buffer, N, n) + +Sample `n` sorted values from `1:N`. +It follows Witter's implementation from "An Efficient Algorithm for Sequential Random Sampling" +""" +function _d_sample!(rng, buffer, N, n) + current_index = 1 + current_sample = 0 + nreal = n + ninv = 1.0 / nreal + Nreal = N + Vprime = exp(log(rand(rng)) * ninv) + qu1 = -n + 1 + N + qu1real = -nreal + 1.0 + Nreal + negalphainv = -13 + threshold = -negalphainv * n + while (n > 1) && (threshold < N) + local negSreal + nmin1inv = 1.0 / (-1.0 + nreal) + while true + local X + while true + X = Nreal * (- Vprime + 1.0) + S = trunc(Int, X) + S < qu1 && break + Vprime = exp(log(rand(rng))*ninv) + end + U = rand(rng) + negSreal = -S + y1 = exp(log(U * Nreal/qu1real) * nmin1inv) + Vprime = y1 * (-X / Nreal + 1.0) * (qu1real / (negSreal + qu1real)) + Vprime <= 1.0 && break + y2 = 1.0 + top = -1.0 + Nreal + if - 1 + n > S + bottom = -nreal + Nreal + limit = -S + N + else + bottom = -1.0 + negSreal + Nreal + limit = qu1 + end + for t in (-1 + N):limit + y2 = (y2 * top) / bottom + top -= 1.0 + bottom -= 1.0 + end + if Nreal/(-X + Nreal) >= y1 * exp(log(y2) * nmin1inv) + Vprime = exp(log(rand(rng)) * nmin1inv) + break + end + Vprime = exp(log(rand(rng)) * ninv) + end + # Skip over the next S records and select the following one for the sample + current_sample += S + 1 + buffer[current_index] = current_sample + current_index += 1 + N = -S + (-1 + N) + Nreal = negSreal+ (-1.0 + Nreal) + n -= 1 + nreal += -1.0 + ninv = nmin1inv + qu1 = -S + qu1 + qu1real += negSreal + threshold += negalphainv + end + if n > 1 + # then Use Method A to finish the sampling + _a_sample!(rng, buffer, N, n, current_index, current_sample) + else # Special case n = 1 + S = trunc(Int, N * Vprime) + # Skip over the next S records and select the following one for the sample + current_sample += S + 1 + buffer[current_index] = current_sample + current_index += 1 + end +end + +function _a_sample!(rng, buffer, N, n, current_index, current_sample) + top = N - n + Nreal = N + while n >= 2 + V = rand(rng) + S = 0 + quot = top / Nreal + while quot > V + S += 1 + top -= 1.0 + Nreal -= 1.0 + quot = (quot * top) / Nreal + end + # Skip over the next S records and select the following one for the sample + current_sample += S + 1 + buffer[current_index] = current_sample + current_index += 1 + Nreal -= 1.0 + n -= 1 + end + # Special case n = 1 + S = trunc(Int, round(Int, Nreal) * rand(rng)) + # Skip over the next S records and select the following one for the sample + current_sample += S + 1 + buffer[current_index] = current_sample + current_index += 1 +end """ - SimpleGraph{T}(nv, ne; rng=nothing, seed=nothing) + SimpleGraph{T}(nv, ne; self_loops = false, rng=nothing, seed=nothing) Construct a random `SimpleGraph{T}` with `nv` vertices and `ne` edges. The graph is sampled uniformly from all such graphs. @@ -22,38 +129,72 @@ julia> SimpleGraph(5, 7) ``` """ function SimpleGraph{T}( - nv::Integer, - ne::Integer; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, -) where {T<:Integer} - tnv = T(nv) - maxe = div(Int(nv) * (nv - 1), 2) - @assert(ne <= maxe, "Maximum number of edges for this graph is $maxe") - rng = rng_from_rng_or_seed(rng, seed) - ne > div((2 * maxe), 3) && return complement(SimpleGraph(tnv, maxe - ne; rng=rng)) - - g = SimpleGraph(tnv) - - while g.ne < ne - source = rand(rng, one(T):tnv) - dest = rand(rng, one(T):tnv) - source != dest && add_edge!(g, source, dest) - end - return g + nv::Integer, + ne::Integer; + self_loops = false, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) where {T <: Integer} + maxe = self_loops ? div(Int(nv) * (nv + 1), 2) : div(Int(nv) * (nv - 1), 2) + + @assert(ne <= maxe, "Maximum number of edges for this graph is $maxe") + + remaining_edges(u) = maxe - (self_loops ? div(Int(nv - u) * (nv - u + 1), 2) : div(Int(nv - u) * (nv - u - 1), 2)) + rng = rng_from_rng_or_seed(rng, seed) + + edge_sample = Array{Int}(undef, max(ne, nv)) + _d_sample!(rng, edge_sample, maxe, ne) + fadjlist = Vector{Vector{T}}(undef, nv) + + first_index = 1 + current_index = 1 + indeg = zeros(T, nv) + + @inbounds for u in 1:nv + while current_index <= ne && remaining_edges(u) >= edge_sample[current_index] + current_index += 1 + end + indeg_u = indeg[u] + list_u = Vector{T}(undef, indeg_u + current_index-first_index) + for j in first_index:current_index-1 + v = edge_sample[j]-remaining_edges(u-1) + u + if self_loops + v -= 1 + end + list_u[j - first_index + 1 + indeg_u] = v + indeg[v] += 1 + end + fadjlist[u] = list_u + first_index = current_index + end + + insert_at = resize!(edge_sample, 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{T}(ne, fadjlist) end function SimpleGraph( - nv::T, - ne::Integer; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, -) where {T<:Integer} - return SimpleGraph{T}(nv, ne; rng=rng, seed=seed) + nv::T, + ne::Integer; + self_loops = false, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) where {T <: Integer} + return SimpleGraph{T}(nv, ne; self_loops = self_loops, rng = rng, seed = seed) end """ - SimpleDiGraph{T}(nv, ne; rng=nothing, seed=nothing) + SimpleDiGraph{T}(nv, ne; self_loops = false, rng=nothing, seed=nothing) Construct a random `SimpleDiGraph{T}` with `nv` vertices and `ne` edges. The graph is sampled uniformly from all such graphs. @@ -72,68 +213,210 @@ julia> SimpleDiGraph(5, 7) ``` """ function SimpleDiGraph{T}( - nv::Integer, - ne::Integer; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, -) where {T<:Integer} - tnv = T(nv) - maxe = Int(nv) * (nv - 1) - @assert(ne <= maxe, "Maximum number of edges for this graph is $maxe") - rng = rng_from_rng_or_seed(rng, seed) - ne > div((2 * maxe), 3) && return complement(SimpleDiGraph{T}(tnv, maxe - ne; rng=rng)) - g = SimpleDiGraph(tnv) - while g.ne < ne - source = rand(rng, one(T):tnv) - dest = rand(rng, one(T):tnv) - source != dest && add_edge!(g, source, dest) + nv::Integer, + ne::Integer; + self_loops = false, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) where {T <: Integer} + maxd = self_loops ? Int(nv) : nv - 1 + maxe = self_loops ? (Int(nv) * nv) : (Int(nv) * (nv - 1)) + + @assert(ne <= maxe, "Maximum number of edges for this graph is $maxe") + rng = rng_from_rng_or_seed(rng, seed) + + edge_sample = Array{Int}(undef, ne) + _d_sample!(rng, edge_sample, maxe, ne) + fadjlist = Vector{Vector{T}}(undef, nv) + badjlist = Vector{Vector{T}}(undef, nv) + + first_index = 1 + current_index = 1 + + @inbounds for u in 1:nv + while current_index <= ne && maxd*u >= edge_sample[current_index] + current_index += 1 + end + list_u = Vector{T}(undef, current_index-first_index) + for j in first_index:current_index-1 + v = edge_sample[j] - maxd*(u - 1) + if !self_loops && v >= u + v += 1 + end + list_u[j-first_index+1] = v + end + fadjlist[u] = list_u + first_index = current_index + end + + indeg = zeros(T, nv) + @inbounds for u in OneTo(nv) + for v in fadjlist[u] + indeg[v] += 1 + end + end + + @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 g + return SimpleDiGraph(ne, fadjlist, badjlist) end function SimpleDiGraph( - nv::T, - ne::Integer; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, -) where {T<:Integer} - return SimpleDiGraph{Int}(nv, ne; rng=rng, seed=seed) + nv::T, + ne::Integer; + self_loops = false, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) where {T <: Integer} + return SimpleDiGraph{Int}(nv, ne; self_loops = self_loops, rng = rng, seed = seed) end -""" - randbn(n, p; rng=nothing, seed=nothing) -Return a binomially-distributed random number with parameters `n` and `p` and optional `seed`. +# like randsubseq!, but without push! since our buffer is always sufficiently large +# we might change if randsubseq! gets faster in newer versions of julia +@inline function _randsubseq!(rng, buffer, sample_from, p) + buffer_len = 0 + if p >= 0.75 + @inbounds for v in sample_from + # for v in sample_from + if rand(rng) < p + buffer_len += 1 + buffer[buffer_len] = v + end + end + else + L = -1 / log1p(-p) + sample_from_len = length(sample_from) + i = 0 + @inbounds while true + # while true + i += floor(Int, randexp(rng) * L) + 1 + i > sample_from_len && return buffer_len + buffer_len += 1 + buffer[buffer_len] = sample_from[i] + end + end + return buffer_len +end -### References -- "Non-Uniform Random Variate Generation," Luc Devroye, p. 522. Retrieved via http://www.eirene.de/Devroye.pdf. -- http://stackoverflow.com/questions/23561551/a-efficient-binomial-random-number-generator-code-in-java -""" -function randbn( - n::Integer, - p::Real; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, -) - rng = rng_from_rng_or_seed(rng, seed) - log_q = log(1.0 - p) - x = 0 - sum = 0.0 - for _ in 1:n - sum += log(rand(rng)) / (n - x) - sum < log_q && break - x += 1 - end - return x +function _erdos_renyi_directed( + nv::T, + p::Union{Rational, AbstractFloat, AbstractIrrational}; + self_loops = false, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) where {T <: Integer} + rng = rng_from_rng_or_seed(rng, seed) + + @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) + buffer_len = _randsubseq!(rng, buffer, sample_from, p) + fadjlist[u] = copy(buffer[1:buffer_len]) + ne += buffer_len + + if !self_loops && u < nv + 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) + 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, 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, + p::Union{Rational, AbstractFloat, AbstractIrrational}; + self_loops = false, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) where {T <: Integer} + rng = rng_from_rng_or_seed(rng, seed) + + @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 = _randsubseq!(rng, 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 """ - erdos_renyi(n, p) + erdos_renyi(n, p) Create an [Erdős–Rényi](http://en.wikipedia.org/wiki/Erdős–Rényi_model) random graph with `n` vertices. Edges are added between pairs of vertices with probability `p`. ### Optional Arguments +- `self_loops=false`: if true, self loops will also be sampled - `is_directed=false`: if true, return a directed graph. - `rng=nothing`: set the Random Number Generator. - `seed=nothing`: set the RNG seed. @@ -153,44 +436,28 @@ julia> erdos_renyi(10, 0.5, is_directed=true, seed=123) ``` """ function erdos_renyi( - n::Integer, - p::Real; - is_directed=false, - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + n::Integer, + p::Real; + is_directed = false, + self_loops = false, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - p >= 1 && return is_directed ? complete_digraph(n) : complete_graph(n) - rng = getRNG(seed) - m = 0 - fadj = [Int[] for u in 1:n] - if is_directed - badj = [Int[] for u in 1:n] - end - for u in 1:n - start = is_directed ? 1 : u+1 - for v in start:n - v == u && continue - if rand(rng) < p - m += 1 - push!(fadj[u], v) - if is_directed - push!(badj[v], u) - else - push!(fadj[v], u) - end - end - end - end - return is_directed ? SimpleDiGraph(m, fadj, badj) : SimpleGraph(m, fadj) + if is_directed + return _erdos_renyi_directed(n, p; self_loops = self_loops, rng = rng, seed = seed) + else + return _erdos_renyi_undirected(n, p; self_loops = self_loops, rng = rng, seed = seed) + end end """ - erdos_renyi(n, ne) + erdos_renyi(n, ne) Create an [Erdős–Rényi](http://en.wikipedia.org/wiki/Erdős–Rényi_model) random graph with `n` vertices and `ne` edges. ### Optional Arguments +- `self_loops=false`: if true, self loops will also be sampled - `is_directed=false`: if true, return a directed graph. - `rng=nothing`: set the Random Number Generator. - `seed=nothing`: set the RNG seed. @@ -207,21 +474,22 @@ julia> erdos_renyi(10, 30, is_directed=true, seed=123) ``` """ function erdos_renyi( - n::Integer, - ne::Integer; - is_directed=false, - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + n::Integer, + ne::Integer; + is_directed = false, + self_loops = false, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - return if is_directed - SimpleDiGraph(n, ne; rng=rng, seed=seed) - else - SimpleGraph(n, ne; rng=rng, seed=seed) - end + return if is_directed + SimpleDiGraph(n, ne; self_loops = self_loops, rng = rng, seed = seed) + else + SimpleGraph(n, ne; self_loops = self_loops, rng = rng, seed = seed) + end end """ - expected_degree_graph(ω) + expected_degree_graph(ω) Given a vector of expected degrees `ω` indexed by vertex, create a random undirected graph in which vertices `i` and `j` are connected with probability `ω[i]*ω[j]/sum(ω)`. @@ -258,50 +526,50 @@ julia> print(degree(g)) ``` """ function expected_degree_graph( - ω::Vector{T}; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, -) where {T<:Real} - g = SimpleGraph(length(ω)) - return expected_degree_graph!(g, ω; rng=rng, seed=seed) + ω::Vector{T}; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) where {T <: Real} + g = SimpleGraph(length(ω)) + return expected_degree_graph!(g, ω; rng = rng, seed = seed) end function expected_degree_graph!( - g::SimpleGraph, - ω::Vector{T}; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, -) where {T<:Real} - n = length(ω) - @assert all(zero(T) .<= ω .<= n - one(T)) "Elements of ω needs to be at least 0 and at most n-1" - - π = sortperm(ω; rev=true) - rng = rng_from_rng_or_seed(rng, seed) - - S = sum(ω) - - for u in 1:(n - 1) - v = u + 1 - p = min(ω[π[u]] * ω[π[v]] / S, one(T)) - while v <= n && p > zero(p) - if p != one(T) - v += floor(Int, log(rand(rng)) / log(one(T) - p)) - end - if v <= n - q = min(ω[π[u]] * ω[π[v]] / S, one(T)) - if rand(rng) < q / p - add_edge!(g, π[u], π[v]) - end - p = q - v += 1 - end - end - end - return g + g::SimpleGraph, + ω::Vector{T}; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) where {T <: Real} + n = length(ω) + @assert all(zero(T) .<= ω .<= n - one(T)) "Elements of ω needs to be at least 0 and at most n-1" + + π = sortperm(ω; rev = true) + rng = rng_from_rng_or_seed(rng, seed) + + S = sum(ω) + + for u in 1:(n-1) + v = u + 1 + p = min(ω[π[u]] * ω[π[v]] / S, one(T)) + while v <= n && p > zero(p) + if p != one(T) + v += floor(Int, log(rand(rng)) / log(one(T) - p)) + end + if v <= n + q = min(ω[π[u]] * ω[π[v]] / S, one(T)) + if rand(rng) < q / p + add_edge!(g, π[u], π[v]) + end + p = q + v += 1 + end + end + end + return g end """ - watts_strogatz(n, k, β) + watts_strogatz(n, k, β) Return a [Watts-Strogatz](https://en.wikipedia.org/wiki/Watts_and_Strogatz_model) small world random graph with `n` vertices, each with expected degree `k` @@ -345,66 +613,66 @@ julia> watts_strogatz(Int8(10), 4, 0.8, is_directed=true, seed=123) - Small Worlds, Duncan J. watts. [https://en.wikipedia.org/wiki/Special:BookSources?isbn=978-0691005416](https://en.wikipedia.org/wiki/Special:BookSources?isbn=978-0691005416) """ function watts_strogatz( - n::Integer, - k::Integer, - β::Real; - is_directed::Bool=false, - remove_edges::Bool=true, - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + n::Integer, + k::Integer, + β::Real; + is_directed::Bool = false, + remove_edges::Bool = true, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - @assert k < n - - # If we have n - 1 neighbors (exactly k/2 on each side), then the graph is - # necessarily complete. No need to run the Watts-Strogatz procedure: - if k == n - 1 && iseven(k) - return is_directed ? complete_digraph(n) : complete_graph(n) - end - - g = is_directed ? SimpleDiGraph(n) : SimpleGraph(n) - - # The ith next vertex, in clockwise order. - # (Reduce to zero-based indexing, so the modulo works, by subtracting 1 - # before and adding 1 after.) - @inline target(s, i) = ((s + i - 1) % n) + 1 - - # Phase 1: For each step size i, add an edge from each vertex s to the ith - # next vertex, in clockwise order. - - for i in 1:div(k, 2), s in 1:n - add_edge!(g, s, target(s, i)) - end - - # Phase 2: For each step size i and each vertex s, consider the edge to the - # ith next vertex, in clockwise order. With probability β, delete the edge - # and rewire it to any (valid) target, chosen uniformly at random. - - rng = rng_from_rng_or_seed(rng, seed) - for i in 1:div(k, 2), s in 1:n - - # We only rewire with a probability β, and we only worry about rewiring - # if there is some vertex not connected to s; otherwise, the only valid - # rewiring is to reconnect to the ith next vertex, and there is no work - # to do. - (rand(rng) < β && degree(g, s) < n - 1) || continue - - t = target(s, i) - - while true - d = rand(rng, 1:n) # Tentative new target - d == s && continue # Self-loops prohibited - d == t && break # Rewired to original target - if add_edge!(g, s, d) # Was this valid (i.e., unconnected)? - remove_edges && rem_edge!(g, s, t) # True rewiring: Delete original edge - break # We found a valid target - end - end - end - return g + @assert k < n + + # If we have n - 1 neighbors (exactly k/2 on each side), then the graph is + # necessarily complete. No need to run the Watts-Strogatz procedure: + if k == n - 1 && iseven(k) + return is_directed ? complete_digraph(n) : complete_graph(n) + end + + g = is_directed ? SimpleDiGraph(n) : SimpleGraph(n) + + # The ith next vertex, in clockwise order. + # (Reduce to zero-based indexing, so the modulo works, by subtracting 1 + # before and adding 1 after.) + @inline target(s, i) = ((s + i - 1) % n) + 1 + + # Phase 1: For each step size i, add an edge from each vertex s to the ith + # next vertex, in clockwise order. + + for i in 1:div(k, 2), s in 1:n + add_edge!(g, s, target(s, i)) + end + + # Phase 2: For each step size i and each vertex s, consider the edge to the + # ith next vertex, in clockwise order. With probability β, delete the edge + # and rewire it to any (valid) target, chosen uniformly at random. + + rng = rng_from_rng_or_seed(rng, seed) + for i in 1:div(k, 2), s in 1:n + + # We only rewire with a probability β, and we only worry about rewiring + # if there is some vertex not connected to s; otherwise, the only valid + # rewiring is to reconnect to the ith next vertex, and there is no work + # to do. + (rand(rng) < β && degree(g, s) < n - 1) || continue + + t = target(s, i) + + while true + d = rand(rng, 1:n) # Tentative new target + d == s && continue # Self-loops prohibited + d == t && break # Rewired to original target + if add_edge!(g, s, d) # Was this valid (i.e., unconnected)? + remove_edges && rem_edge!(g, s, t) # True rewiring: Delete original edge + break # We found a valid target + end + end + end + return g end """ - newman_watts_strogatz(n, k, β) + newman_watts_strogatz(n, k, β) Return a Newman-Watts-Strogatz small world random graph with `n` vertices, each with expected degree `k(1 + β)` (or `(k - 1)(1 + β)` if `k` is odd). Edges are @@ -451,73 +719,73 @@ julia> newman_watts_strogatz(Int8(10), 4, 0.8, is_directed=true, seed=123) - Scaling and percolation in the small-world network model, M. E. J. Newman, Duncan J. Watts. [https://doi.org/10.1103/PhysRevE.60.7332](https://doi.org/10.1103/PhysRevE.60.7332) """ function newman_watts_strogatz( - n::Integer, - k::Integer, - β::Real; - is_directed::Bool=false, - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + n::Integer, + k::Integer, + β::Real; + is_directed::Bool = false, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - return watts_strogatz( - n, k, β; is_directed=is_directed, remove_edges=false, rng=rng, seed=seed - ) + return watts_strogatz( + n, k, β; is_directed = is_directed, remove_edges = false, rng = rng, seed = seed, + ) end -function _suitable(edges::Set{SimpleEdge{T}}, potential_edges::Dict{T,T}) where {T<:Integer} - isempty(potential_edges) && return true - list = keys(potential_edges) - for s1 in list, s2 in list - s1 >= s2 && continue - (SimpleEdge(s1, s2) ∉ edges) && return true - end - return false +function _suitable(edges::Set{SimpleEdge{T}}, potential_edges::Dict{T, T}) where {T <: Integer} + isempty(potential_edges) && return true + list = keys(potential_edges) + for s1 in list, s2 in list + s1 >= s2 && continue + (SimpleEdge(s1, s2) ∉ edges) && return true + end + return false end _try_creation(n::Integer, k::Integer, rng::AbstractRNG) = _try_creation(n, fill(k, n), rng) -function _try_creation(n::T, k::Vector{T}, rng::AbstractRNG) where {T<:Integer} - edges = Set{SimpleEdge{T}}() - m = 0 - stubs = zeros(T, sum(k)) - for i in one(T):n - for j in one(T):k[i] - m += 1 - stubs[m] = i - end - end - # stubs = vcat([fill(i, k[i]) for i = 1:n]...) # slower - - while !isempty(stubs) - potential_edges = Dict{T,T}() - shuffle!(rng, stubs) - for i in 1:2:length(stubs) - s1, s2 = stubs[i:(i + 1)] - if (s1 > s2) - s1, s2 = s2, s1 - end - e = SimpleEdge(s1, s2) - if s1 != s2 && ∉(e, edges) - push!(edges, e) - else - potential_edges[s1] = get(potential_edges, s1, 0) + 1 - potential_edges[s2] = get(potential_edges, s2, 0) + 1 - end - end - - if !_suitable(edges, potential_edges) - return Set{SimpleEdge{T}}() - end - - stubs = Vector{T}() - for (e, ct) in potential_edges - append!(stubs, fill(e, ct)) - end - end - return edges +function _try_creation(n::T, k::Vector{T}, rng::AbstractRNG) where {T <: Integer} + edges = Set{SimpleEdge{T}}() + m = 0 + stubs = zeros(T, sum(k)) + for i in one(T):n + for j in one(T):k[i] + m += 1 + stubs[m] = i + end + end + # stubs = vcat([fill(i, k[i]) for i = 1:n]...) # slower + + while !isempty(stubs) + potential_edges = Dict{T, T}() + shuffle!(rng, stubs) + for i in 1:2:length(stubs) + s1, s2 = stubs[i:(i+1)] + if (s1 > s2) + s1, s2 = s2, s1 + end + e = SimpleEdge(s1, s2) + if s1 != s2 && ∉(e, edges) + push!(edges, e) + else + potential_edges[s1] = get(potential_edges, s1, 0) + 1 + potential_edges[s2] = get(potential_edges, s2, 0) + 1 + end + end + + if !_suitable(edges, potential_edges) + return Set{SimpleEdge{T}}() + end + + stubs = Vector{T}() + for (e, ct) in potential_edges + append!(stubs, fill(e, ct)) + end + end + return edges end """ - barabasi_albert(n, k) + barabasi_albert(n, k) Create a [Barabási–Albert model](https://en.wikipedia.org/wiki/Barab%C3%A1si%E2%80%93Albert_model) random graph with `n` vertices. It is grown by adding new vertices to an initial @@ -544,7 +812,7 @@ julia> barabasi_albert(100, Int8(10), is_directed=true, complete=true, seed=123) barabasi_albert(n::Integer, k::Integer; keyargs...) = barabasi_albert(n, k, k; keyargs...) """ - barabasi_albert(n::Integer, n0::Integer, k::Integer) + barabasi_albert(n::Integer, n0::Integer, k::Integer) Create a [Barabási–Albert model](https://en.wikipedia.org/wiki/Barab%C3%A1si%E2%80%93Albert_model) random graph with `n` vertices. It is grown by adding new vertices to an initial @@ -570,26 +838,26 @@ julia> barabasi_albert(100, Int8(10), 3, is_directed=true, seed=123) ``` """ function barabasi_albert( - n::Integer, - n0::Integer, - k::Integer; - is_directed::Bool=false, - complete::Bool=false, - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + n::Integer, + n0::Integer, + k::Integer; + is_directed::Bool = false, + complete::Bool = false, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - if complete - g = is_directed ? complete_digraph(n0) : complete_graph(n0) - else - g = is_directed ? SimpleDiGraph(n0) : SimpleGraph(n0) - end - - barabasi_albert!(g, n, k; rng=rng, seed=seed) - return g + if complete + g = is_directed ? complete_digraph(n0) : complete_graph(n0) + else + g = is_directed ? SimpleDiGraph(n0) : SimpleGraph(n0) + end + + barabasi_albert!(g, n, k; rng = rng, seed = seed) + return g end """ - barabasi_albert!(g::AbstractGraph, n::Integer, k::Integer) + barabasi_albert!(g::AbstractGraph, n::Integer, k::Integer) Create a [Barabási–Albert model](https://en.wikipedia.org/wiki/Barab%C3%A1si%E2%80%93Albert_model) random graph with `n` vertices. It is grown by adding new vertices to an initial @@ -613,78 +881,78 @@ julia> g ``` """ function barabasi_albert!( - g::AbstractGraph, - n::Integer, - k::Integer; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + g::AbstractGraph, + n::Integer, + k::Integer; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - n0 = nv(g) - 1 <= k <= n0 <= n || - throw(ArgumentError("Barabási-Albert model requires 1 <= k <= nv(g) <= n")) - n0 == n && return g - - # seed random number generator - rng = rng_from_rng_or_seed(rng, seed) - - # add missing vertices - sizehint!(g.fadjlist, n) - add_vertices!(g, n - n0) - - # if initial graph doesn't contain any edges - # expand it by one vertex and add k edges from this additional node - if ne(g) == 0 - # expand initial graph - n0 += one(n0) - - # add edges to k existing vertices - for target in sample!(rng, collect(1:(n0 - 1)), k) - add_edge!(g, n0, target) - end - end - - # vector of weighted vertices (each node is repeated once for each adjacent edge) - weightedVs = Vector{Int}(undef, 2 * (n - n0) * k + 2 * ne(g)) - - # initialize vector of weighted vertices - offset = 0 - for e in edges(g) - weightedVs[offset += 1] = src(e) - weightedVs[offset += 1] = dst(e) - end - - # array to record if a node is picked - picked = fill(false, n) - - # vector of targets - targets = Vector{Int}(undef, k) - - for source in (n0 + 1):n - # choose k targets from the existing vertices - # pick uniformly from weightedVs (preferential attachment) - i = 0 - while i < k - target = weightedVs[rand(rng, 1:offset)] - if !picked[target] - targets[i += 1] = target - picked[target] = true - end - end - - # add edges to k targets - for target in targets - add_edge!(g, source, target) - - weightedVs[offset += 1] = source - weightedVs[offset += 1] = target - picked[target] = false - end - end - return g + n0 = nv(g) + 1 <= k <= n0 <= n || + throw(ArgumentError("Barabási-Albert model requires 1 <= k <= nv(g) <= n")) + n0 == n && return g + + # seed random number generator + rng = rng_from_rng_or_seed(rng, seed) + + # add missing vertices + sizehint!(g.fadjlist, n) + add_vertices!(g, n - n0) + + # if initial graph doesn't contain any edges + # expand it by one vertex and add k edges from this additional node + if ne(g) == 0 + # expand initial graph + n0 += one(n0) + + # add edges to k existing vertices + for target in sample!(rng, collect(1:(n0-1)), k) + add_edge!(g, n0, target) + end + end + + # vector of weighted vertices (each node is repeated once for each adjacent edge) + weightedVs = Vector{Int}(undef, 2 * (n - n0) * k + 2 * ne(g)) + + # initialize vector of weighted vertices + offset = 0 + for e in edges(g) + weightedVs[offset+=1] = src(e) + weightedVs[offset+=1] = dst(e) + end + + # array to record if a node is picked + picked = fill(false, n) + + # vector of targets + targets = Vector{Int}(undef, k) + + for source in (n0+1):n + # choose k targets from the existing vertices + # pick uniformly from weightedVs (preferential attachment) + i = 0 + while i < k + target = weightedVs[rand(rng, 1:offset)] + if !picked[target] + targets[i+=1] = target + picked[target] = true + end + end + + # add edges to k targets + for target in targets + add_edge!(g, source, target) + + weightedVs[offset+=1] = source + weightedVs[offset+=1] = target + picked[target] = false + end + end + return g end """ - static_fitness_model(m, fitness) + static_fitness_model(m, fitness) Generate a random graph with ``|fitness|`` vertices and `m` edges, in which the probability of the existence of ``Edge_{ij}`` is proportional @@ -729,33 +997,33 @@ julia> edges(g) |> collect ``` """ function static_fitness_model( - m::Integer, - fitness::Vector{T}; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, -) where {T<:Real} - m < 0 && throw(ArgumentError("number of edges must be positive")) - n = length(fitness) - m == 0 && return SimpleGraph(n) - nvs = 0 - for f in fitness - # sanity check for the fitness - f < zero(T) && throw(ArgumentError("fitness scores must be non-negative")) - f > zero(T) && (nvs += 1) - end - # avoid getting into an infinite loop when too many edges are requested - max_no_of_edges = div(nvs * (nvs - 1), 2) - m > max_no_of_edges && - throw(ArgumentError("too many edges requested ($m > $max_no_of_edges)")) - # calculate the cumulative fitness scores - cum_fitness = cumsum(fitness) - g = SimpleGraph(n) - _create_static_fitness_graph!(g, m, cum_fitness, cum_fitness, rng, seed) - return g + m::Integer, + fitness::Vector{T}; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) where {T <: Real} + m < 0 && throw(ArgumentError("number of edges must be positive")) + n = length(fitness) + m == 0 && return SimpleGraph(n) + nvs = 0 + for f in fitness + # sanity check for the fitness + f < zero(T) && throw(ArgumentError("fitness scores must be non-negative")) + f > zero(T) && (nvs += 1) + end + # avoid getting into an infinite loop when too many edges are requested + max_no_of_edges = div(nvs * (nvs - 1), 2) + m > max_no_of_edges && + throw(ArgumentError("too many edges requested ($m > $max_no_of_edges)")) + # calculate the cumulative fitness scores + cum_fitness = cumsum(fitness) + g = SimpleGraph(n) + _create_static_fitness_graph!(g, m, cum_fitness, cum_fitness, rng, seed) + return g end """ - static_fitness_model(m, fitness_out, fitness_in) + static_fitness_model(m, fitness_out, fitness_in) Generate a random directed graph with ``|fitness\\_out + fitness\\_in|`` vertices and `m` edges, in which the probability of the existence of ``Edge_{ij}`` is proportional with @@ -789,62 +1057,62 @@ julia> edges(g) |> collect ``` """ function static_fitness_model( - m::Integer, - fitness_out::Vector{T}, - fitness_in::Vector{S}; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, -) where {T<:Real} where {S<:Real} - m < 0 && throw(ArgumentError("number of edges must be positive")) - n = length(fitness_out) - length(fitness_in) != n && - throw(ArgumentError("fitness_in must have the same size as fitness_out")) - m == 0 && return SimpleDiGraph(n) - # avoid getting into an infinite loop when too many edges are requested - noutvs = ninvs = nvs = 0 - @inbounds for i in 1:n - # sanity check for the fitness - (fitness_out[i] < zero(T) || fitness_in[i] < zero(S)) && - error("fitness scores must be non-negative") # TODO 0.7: change to DomainError? - fitness_out[i] > zero(T) && (noutvs += 1) - fitness_in[i] > zero(S) && (ninvs += 1) - (fitness_out[i] > zero(T) && fitness_in[i] > zero(S)) && (nvs += 1) - end - max_no_of_edges = noutvs * ninvs - nvs - m > max_no_of_edges && - throw(ArgumentError("too many edges requested ($m > $max_no_of_edges)")) - # calculate the cumulative fitness scores - cum_fitness_out = cumsum(fitness_out) - cum_fitness_in = cumsum(fitness_in) - g = SimpleDiGraph(n) - _create_static_fitness_graph!(g, m, cum_fitness_out, cum_fitness_in, rng, seed) - return g + m::Integer, + fitness_out::Vector{T}, + fitness_in::Vector{S}; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) where {T <: Real} where {S <: Real} + m < 0 && throw(ArgumentError("number of edges must be positive")) + n = length(fitness_out) + length(fitness_in) != n && + throw(ArgumentError("fitness_in must have the same size as fitness_out")) + m == 0 && return SimpleDiGraph(n) + # avoid getting into an infinite loop when too many edges are requested + noutvs = ninvs = nvs = 0 + @inbounds for i in 1:n + # sanity check for the fitness + (fitness_out[i] < zero(T) || fitness_in[i] < zero(S)) && + error("fitness scores must be non-negative") # TODO 0.7: change to DomainError? + fitness_out[i] > zero(T) && (noutvs += 1) + fitness_in[i] > zero(S) && (ninvs += 1) + (fitness_out[i] > zero(T) && fitness_in[i] > zero(S)) && (nvs += 1) + end + max_no_of_edges = noutvs * ninvs - nvs + m > max_no_of_edges && + throw(ArgumentError("too many edges requested ($m > $max_no_of_edges)")) + # calculate the cumulative fitness scores + cum_fitness_out = cumsum(fitness_out) + cum_fitness_in = cumsum(fitness_in) + g = SimpleDiGraph(n) + _create_static_fitness_graph!(g, m, cum_fitness_out, cum_fitness_in, rng, seed) + return g end function _create_static_fitness_graph!( - g::AbstractGraph, - m::Integer, - cum_fitness_out::Vector{T}, - cum_fitness_in::Vector{S}, - rng::Union{Nothing,AbstractRNG}, - seed::Union{Nothing,Integer}, -) where {T<:Real} where {S<:Real} - max_out = cum_fitness_out[end] - max_in = cum_fitness_in[end] - rng = rng_from_rng_or_seed(rng, seed) - while m > 0 - source = searchsortedfirst(cum_fitness_out, rand(rng) * max_out) - target = searchsortedfirst(cum_fitness_in, rand(rng) * max_in) - # skip if loop edge - (source == target) && continue - # is there already an edge? If so, try again - add_edge!(g, source, target) || continue - m -= one(m) - end + g::AbstractGraph, + m::Integer, + cum_fitness_out::Vector{T}, + cum_fitness_in::Vector{S}, + rng::Union{Nothing, AbstractRNG}, + seed::Union{Nothing, Integer}, +) where {T <: Real} where {S <: Real} + max_out = cum_fitness_out[end] + max_in = cum_fitness_in[end] + rng = rng_from_rng_or_seed(rng, seed) + while m > 0 + source = searchsortedfirst(cum_fitness_out, rand(rng) * max_out) + target = searchsortedfirst(cum_fitness_in, rand(rng) * max_in) + # skip if loop edge + (source == target) && continue + # is there already an edge? If so, try again + add_edge!(g, source, target) || continue + m -= one(m) + end end """ - static_scale_free(n, m, α) + static_scale_free(n, m, α) Generate a random graph with `n` vertices, `m` edges and expected power-law degree distribution with exponent `α`. @@ -864,21 +1132,21 @@ Time complexity is ``\\mathcal{O}(|V| + |E| log |E|)``. - Cho YS, Kim JS, Park J, Kahng B, Kim D: Percolation transitions in scale-free networks under the Achlioptas process. Phys Rev Lett 103:135702, 2009. """ function static_scale_free( - n::Integer, - m::Integer, - α::Real; - finite_size_correction::Bool=true, - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + n::Integer, + m::Integer, + α::Real; + finite_size_correction::Bool = true, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - n < 0 && throw(ArgumentError("number of vertices must be positive")) - α < 2 && throw(ArgumentError("out-degree exponent must be >= 2")) - fitness = _construct_fitness(n, α, finite_size_correction) - return static_fitness_model(m, fitness; rng=rng, seed=seed) + n < 0 && throw(ArgumentError("number of vertices must be positive")) + α < 2 && throw(ArgumentError("out-degree exponent must be >= 2")) + fitness = _construct_fitness(n, α, finite_size_correction) + return static_fitness_model(m, fitness; rng = rng, seed = seed) end """ - static_scale_free(n, m, α_out, α_in) + static_scale_free(n, m, α_out, α_in) Generate a random graph with `n` vertices, `m` edges and expected power-law degree distribution with exponent `α_out` for outbound edges and `α_in` for @@ -899,43 +1167,43 @@ Time complexity is ``\\mathcal{O}(|V| + |E| log |E|)``. - Cho YS, Kim JS, Park J, Kahng B, Kim D: Percolation transitions in scale-free networks under the Achlioptas process. Phys Rev Lett 103:135702, 2009. """ function static_scale_free( - n::Integer, - m::Integer, - α_out::Real, - α_in::Float64; - finite_size_correction::Bool=true, - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + n::Integer, + m::Integer, + α_out::Real, + α_in::Float64; + finite_size_correction::Bool = true, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - n < 0 && throw(ArgumentError("number of vertices must be positive")) - α_out < 2 && throw(ArgumentError("out-degree exponent must be >= 2")) - α_in < 2 && throw(ArgumentError("out-degree exponent must be >= 2")) - # construct the fitness - fitness_out = _construct_fitness(n, α_out, finite_size_correction) - fitness_in = _construct_fitness(n, α_in, finite_size_correction) - # eliminate correlation - shuffle!(fitness_in) - return static_fitness_model(m, fitness_out, fitness_in; rng=rng, seed=seed) + n < 0 && throw(ArgumentError("number of vertices must be positive")) + α_out < 2 && throw(ArgumentError("out-degree exponent must be >= 2")) + α_in < 2 && throw(ArgumentError("out-degree exponent must be >= 2")) + # construct the fitness + fitness_out = _construct_fitness(n, α_out, finite_size_correction) + fitness_in = _construct_fitness(n, α_in, finite_size_correction) + # eliminate correlation + shuffle!(fitness_in) + return static_fitness_model(m, fitness_out, fitness_in; rng = rng, seed = seed) end function _construct_fitness(n::Integer, α::Real, finite_size_correction::Bool) - α = -1 / (α - 1) - fitness = zeros(n) - j = float(n) - if finite_size_correction && α < -0.5 - # See the Cho et al paper, first page first column + footnote 7 - j += n^(1 + 1 / 2 * α) * (10 * sqrt(2) * (1 + α))^(-1 / α) - 1 - end - j = max(j, n) - @inbounds for i in 1:n - fitness[i] = j^α - j -= 1 - end - return fitness + α = -1 / (α - 1) + fitness = zeros(n) + j = float(n) + if finite_size_correction && α < -0.5 + # See the Cho et al paper, first page first column + footnote 7 + j += n^(1 + 1 / 2 * α) * (10 * sqrt(2) * (1 + α))^(-1 / α) - 1 + end + j = max(j, n) + @inbounds for i in 1:n + fitness[i] = j^α + j -= 1 + end + return fitness end """ - random_regular_graph(n, k) + random_regular_graph(n, k) Create a random undirected [regular graph](https://en.wikipedia.org/wiki/Regular_graph) with `n` vertices, @@ -953,36 +1221,36 @@ Allocates an array of `nk` `Int`s, and . For ``k > \\frac{n}{2}``, generates a g ``n-k-1`` and returns its complement. """ function random_regular_graph( - n::Integer, - k::Integer; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + n::Integer, + k::Integer; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - !iseven(n * k) && throw(ArgumentError("n * k must be even")) - !(0 <= k < n) && throw(ArgumentError("the 0 <= k < n inequality must be satisfied")) - if k == 0 - return SimpleGraph(n) - end - rng = rng_from_rng_or_seed(rng, seed) - if (k > n / 2) && iseven(n * (n - k - 1)) - return complement(random_regular_graph(n, n - k - 1; rng=rng)) - end - - edges = _try_creation(n, k, rng) - while isempty(edges) - edges = _try_creation(n, k, rng) - end - - g = SimpleGraph(n) - for edge in edges - add_edge!(g, edge) - end - - return g + !iseven(n * k) && throw(ArgumentError("n * k must be even")) + !(0 <= k < n) && throw(ArgumentError("the 0 <= k < n inequality must be satisfied")) + if k == 0 + return SimpleGraph(n) + end + rng = rng_from_rng_or_seed(rng, seed) + if (k > n / 2) && iseven(n * (n - k - 1)) + return complement(random_regular_graph(n, n - k - 1; rng = rng)) + end + + edges = _try_creation(n, k, rng) + while isempty(edges) + edges = _try_creation(n, k, rng) + end + + g = SimpleGraph(n) + for edge in edges + add_edge!(g, edge) + end + + return g end """ - random_configuration_model(n, ks) + random_configuration_model(n, ks) Create a random undirected graph according to the [configuration model] (http://tuvalu.santafe.edu/~aaronc/courses/5352/fall2013/csci5352_2013_L11.pdf) @@ -1000,35 +1268,35 @@ Time complexity is approximately ``\\mathcal{O}(n \\bar{k}^2)``. Allocates an array of ``n \\bar{k}`` `Int`s. """ function random_configuration_model( - n::Integer, - k::Array{T}; - check_graphical::Bool=false, - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, -) where {T<:Integer} - n != length(k) && throw(ArgumentError("a degree sequence of length n must be provided")) - m = sum(k) - !iseven(m) && throw(ArgumentError("sum(k) must be even")) - !all(0 .<= k .< n) && - throw(ArgumentError("the 0 <= k[i] < n inequality must be satisfied")) - if check_graphical - isgraphical(k) || throw(ArgumentError("degree sequence must be graphical")) - end - rng = rng_from_rng_or_seed(rng, seed) - - edges = _try_creation(n, k, rng) - while m > 0 && isempty(edges) - edges = _try_creation(n, k, rng) - end - - g = SimpleGraphFromIterator(edges) - if nv(g) < n - add_vertices!(g, n - nv(g)) - end - return g + n::Integer, + k::Array{T}; + check_graphical::Bool = false, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) where {T <: Integer} + n != length(k) && throw(ArgumentError("a degree sequence of length n must be provided")) + m = sum(k) + !iseven(m) && throw(ArgumentError("sum(k) must be even")) + !all(0 .<= k .< n) && + throw(ArgumentError("the 0 <= k[i] < n inequality must be satisfied")) + if check_graphical + isgraphical(k) || throw(ArgumentError("degree sequence must be graphical")) + end + rng = rng_from_rng_or_seed(rng, seed) + + edges = _try_creation(n, k, rng) + while m > 0 && isempty(edges) + edges = _try_creation(n, k, rng) + end + + g = SimpleGraphFromIterator(edges) + if nv(g) < n + add_vertices!(g, n - nv(g)) + end + return g end """ - uniform_tree(n) + uniform_tree(n) Generates a random labelled tree, drawn uniformly at random over the ``n^{n-2}`` such trees. A uniform word of length `n-2` over the alphabet `1:n` is generated (Prüfer sequence) then decoded. See also the `prufer_decode` function and [this page on Prüfer codes](https://en.wikipedia.org/wiki/Pr%C3%BCfer_sequence). @@ -1043,16 +1311,16 @@ julia> uniform_tree(10) {10, 9} undirected simple Int64 graph ``` """ -function uniform_tree(n::Integer; rng::Union{Nothing,AbstractRNG}=nothing) - n <= 1 && return Graph(n) - n == 2 && return path_graph(n) - rng = rng_from_rng_or_seed(rng, nothing) - random_code = rand(rng, Base.OneTo(n), n - 2) - return prufer_decode(random_code) +function uniform_tree(n::Integer; rng::Union{Nothing, AbstractRNG} = nothing) + n <= 1 && return Graph(n) + n == 2 && return path_graph(n) + rng = rng_from_rng_or_seed(rng, nothing) + random_code = rand(rng, Base.OneTo(n), n - 2) + return prufer_decode(random_code) end """ - random_regular_digraph(n, k) + random_regular_digraph(n, k) Create a random directed [regular graph](https://en.wikipedia.org/wiki/Regular_graph) with `n` vertices, each with degree `k`. @@ -1067,40 +1335,40 @@ Allocates an ``n × n`` sparse matrix of boolean as an adjacency matrix and uses that to generate the directed graph. """ function random_regular_digraph( - n::Integer, - k::Integer; - dir::Symbol=:out, - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + n::Integer, + k::Integer; + dir::Symbol = :out, + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - # TODO remove the function sample from StatsBase for one allowing the use - # of a local rng - !(0 <= k < n) && throw(ArgumentError("the 0 <= k < n inequality must be satisfied")) - - if k == 0 - return SimpleDiGraph(n) - end - rng = rng_from_rng_or_seed(rng, seed) - if (k > n / 2) && iseven(n * (n - k - 1)) - return complement(random_regular_digraph(n, n - k - 1; dir=dir, rng=rng)) - end - cs = collect(2:n) - i = 1 - I = Vector{Int}(undef, n * k) - J = Vector{Int}(undef, n * k) - V = fill(true, n * k) - for r in 1:n - l = ((r - 1) * k + 1):(r * k) - I[l] .= r - J[l] = sample!(rng, cs, k; exclude=r) - end - - m = dir == :out ? sparse(I, J, V, n, n) : sparse(I, J, V, n, n)' - return SimpleDiGraph(m) + # TODO remove the function sample from StatsBase for one allowing the use + # of a local rng + !(0 <= k < n) && throw(ArgumentError("the 0 <= k < n inequality must be satisfied")) + + if k == 0 + return SimpleDiGraph(n) + end + rng = rng_from_rng_or_seed(rng, seed) + if (k > n / 2) && iseven(n * (n - k - 1)) + return complement(random_regular_digraph(n, n - k - 1; dir = dir, rng = rng)) + end + cs = collect(2:n) + i = 1 + I = Vector{Int}(undef, n * k) + J = Vector{Int}(undef, n * k) + V = fill(true, n * k) + for r in 1:n + l = ((r-1)*k+1):(r*k) + I[l] .= r + J[l] = sample!(rng, cs, k; exclude = r) + end + + m = dir == :out ? sparse(I, J, V, n, n) : sparse(I, J, V, n, n)' + return SimpleDiGraph(m) end """ - random_tournament_digraph(n) + random_tournament_digraph(n) Create a random directed [tournament graph] (https://en.wikipedia.org/wiki/Tournament_%28graph_theory%29) @@ -1122,28 +1390,55 @@ julia> random_tournament_digraph(Int8(10), seed=123) ``` """ function random_tournament_digraph( - n::Integer; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + n::Integer; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - rng = rng_from_rng_or_seed(rng, seed) - g = SimpleDiGraph(n) + rng = rng_from_rng_or_seed(rng, seed) + g = SimpleDiGraph(n) - for i in 1:n, j in (i + 1):n - rand(rng, Bool) ? add_edge!(g, SimpleEdge(i, j)) : add_edge!(g, SimpleEdge(j, i)) - end + for i in 1:n, j in (i+1):n + rand(rng, Bool) ? add_edge!(g, SimpleEdge(i, j)) : add_edge!(g, SimpleEdge(j, i)) + end - return g + return g end """ - stochastic_block_model(c, n) + randbn(n, p; rng=nothing, seed=nothing) + +Return a binomially-distributed random number with parameters `n` and `p` and optional `seed`. + +### References +- "Non-Uniform Random Variate Generation," Luc Devroye, p. 522. Retrieved via http://www.eirene.de/Devroye.pdf. +- http://stackoverflow.com/questions/23561551/a-efficient-binomial-random-number-generator-code-in-java +""" +function randbn( + n::Integer, + p::Real; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) + rng = rng_from_rng_or_seed(rng, seed) + log_q = log(1.0 - p) + x = 0 + sum = 0.0 + for _ in 1:n + sum += log(rand(rng)) / (n - x) + sum < log_q && break + x += 1 + end + return x +end + +""" + stochastic_block_model(c, n) Return a Graph generated according to the Stochastic Block Model (SBM). `c[a,b]` : Mean number of neighbors of a vertex in block `a` belonging to block `b`. - Only the upper triangular part is considered, since the lower triangular is - determined by ``c[b,a] = c[a,b] * \\frac{n[a]}{n[b]}``. + Only the upper triangular part is considered, since the lower triangular is + determined by ``c[b,a] = c[a,b] * \\frac{n[a]}{n[b]}``. `n[a]` : Number of vertices in block `a` ### Optional Arguments @@ -1154,68 +1449,68 @@ For a dynamic version of the SBM see the [`StochasticBlockModel`](@ref) type and related functions. """ function stochastic_block_model( - c::Matrix{T}, - n::Vector{U}; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, -) where {T<:Real} where {U<:Integer} - size(c, 1) == size(c, 2) == length(n) || - throw(ArgumentError("matrix-vector size mismatch")) - - # init dsfmt generator with a fixed seed - rng = rng_from_rng_or_seed(rng, seed) - N = sum(n) - K = length(n) - nedg = zeros(Int, K, K) - g = SimpleGraph(N) - cum = [sum(n[1:a]) for a in 0:K] - for a in 1:K - ra = (cum[a] + 1):cum[a + 1] - for b in a:K - ((a == b) && !(c[a, b] <= n[b] - 1)) || - ((a != b) && !(c[a, b] <= n[b])) && error( - "Mean degree cannot be greater than available neighbors in the block.", - ) # TODO 0.7: turn into some other error? - - m = a == b ? div(n[a] * (n[a] - 1), 2) : n[a] * n[b] - p = a == b ? n[a] * c[a, b] / (2m) : n[a] * c[a, b] / m - nedg = randbn(m, p; rng=rng) - rb = (cum[b] + 1):cum[b + 1] - i = 0 - while i < nedg - source = rand(rng, ra) - dest = rand(rng, rb) - if source != dest - if add_edge!(g, source, dest) - i += 1 - end - end - end - end - end - return g + c::Matrix{T}, + n::Vector{U}; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) where {T <: Real} where {U <: Integer} + size(c, 1) == size(c, 2) == length(n) || + throw(ArgumentError("matrix-vector size mismatch")) + + # init dsfmt generator with a fixed seed + rng = rng_from_rng_or_seed(rng, seed) + N = sum(n) + K = length(n) + nedg = zeros(Int, K, K) + g = SimpleGraph(N) + cum = [sum(n[1:a]) for a in 0:K] + for a in 1:K + ra = (cum[a]+1):cum[a+1] + for b in a:K + ((a == b) && !(c[a, b] <= n[b] - 1)) || + ((a != b) && !(c[a, b] <= n[b])) && error( + "Mean degree cannot be greater than available neighbors in the block.", + ) # TODO 0.7: turn into some other error? + + m = a == b ? div(n[a] * (n[a] - 1), 2) : n[a] * n[b] + p = a == b ? n[a] * c[a, b] / (2m) : n[a] * c[a, b] / m + nedg = randbn(m, p; rng = rng) + rb = (cum[b]+1):cum[b+1] + i = 0 + while i < nedg + source = rand(rng, ra) + dest = rand(rng, rb) + if source != dest + if add_edge!(g, source, dest) + i += 1 + end + end + end + end + end + return g end """ - stochastic_block_model(cint, cext, n) + stochastic_block_model(cint, cext, n) Return a Graph generated according to the Stochastic Block Model (SBM), sampling from an SBM with ``c_{a,a}=cint``, and ``c_{a,b}=cext``. """ function stochastic_block_model( - cint::T, - cext::T, - n::Vector{U}; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, -) where {T<:Real} where {U<:Integer} - K = length(n) - c = [ifelse(a == b, cint, cext) for a in 1:K, b in 1:K] - return stochastic_block_model(c, n; rng=rng, seed=seed) + cint::T, + cext::T, + n::Vector{U}; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) where {T <: Real} where {U <: Integer} + K = length(n) + c = [ifelse(a == b, cint, cext) for a in 1:K, b in 1:K] + return stochastic_block_model(c, n; rng = rng, seed = seed) end """ - StochasticBlockModel{T,P} + StochasticBlockModel{T,P} A type capturing the parameters of the SBM. Each vertex is assigned to a block and the probability of edge `(i,j)` @@ -1231,71 +1526,71 @@ block `k` and any vertex in block `l`. Graphs are generated by taking random ``i,j ∈ V`` and flipping a coin with probability `affinities[nodemap[i],nodemap[j]]`. """ -mutable struct StochasticBlockModel{T<:Integer,P<:Real} - n::T - nodemap::Array{T} - affinities::Matrix{P} +mutable struct StochasticBlockModel{T <: Integer, P <: Real} + n::T + nodemap::Array{T} + affinities::Matrix{P} end function ==(sbm::StochasticBlockModel, other::StochasticBlockModel) - return (sbm.n == other.n) && - (sbm.nodemap == other.nodemap) && - (sbm.affinities == other.affinities) + return (sbm.n == other.n) && + (sbm.nodemap == other.nodemap) && + (sbm.affinities == other.affinities) end # A constructor for StochasticBlockModel that uses the sizes of the blocks # and the affinity matrix. This construction implies that consecutive # vertices will be in the same blocks, except for the block boundaries. function StochasticBlockModel(sizes::AbstractVector, affinities::AbstractMatrix) - csum = cumsum(sizes) - j = 1 - nodemap = zeros(Int, csum[end]) - for i in 1:csum[end] - if i > csum[j] - j += 1 - end - nodemap[i] = j - end - return StochasticBlockModel(csum[end], nodemap, affinities) + csum = cumsum(sizes) + j = 1 + nodemap = zeros(Int, csum[end]) + for i in 1:csum[end] + if i > csum[j] + j += 1 + end + nodemap[i] = j + end + return StochasticBlockModel(csum[end], nodemap, affinities) end ### TODO: This documentation needs work. sbromberger 20170326 """ - sbmaffinity(internalp, externalp, sizes) + sbmaffinity(internalp, externalp, sizes) Produce the sbm affinity matrix with internal probabilities `internalp` and external probabilities `externalp`. """ function sbmaffinity( - internalp::Vector{T}, externalp::Real, sizes::Vector{U} -) where {T<:Real} where {U<:Integer} - numblocks = length(sizes) - numblocks == length(internalp) || - throw(ArgumentError("Inconsistent input dimensions: internalp, sizes")) - B = diagm(0 => internalp) + externalp * (ones(numblocks, numblocks) - I) - return B + internalp::Vector{T}, externalp::Real, sizes::Vector{U}, +) where {T <: Real} where {U <: Integer} + numblocks = length(sizes) + numblocks == length(internalp) || + throw(ArgumentError("Inconsistent input dimensions: internalp, sizes")) + B = diagm(0 => internalp) + externalp * (ones(numblocks, numblocks) - I) + return B end function StochasticBlockModel( - internalp::Real, externalp::Real, size::Integer, numblocks::Integer + internalp::Real, externalp::Real, size::Integer, numblocks::Integer, ) - sizes = fill(size, numblocks) - B = sbmaffinity(fill(internalp, numblocks), externalp, sizes) - return StochasticBlockModel(sizes, B) + sizes = fill(size, numblocks) + B = sbmaffinity(fill(internalp, numblocks), externalp, sizes) + return StochasticBlockModel(sizes, B) end function StochasticBlockModel( - internalp::Vector{T}, externalp::Real, sizes::Vector{U} -) where {T<:Real,U<:Integer} - B = sbmaffinity(internalp, externalp, sizes) - return StochasticBlockModel(sizes, B) + internalp::Vector{T}, externalp::Real, sizes::Vector{U}, +) where {T <: Real, U <: Integer} + B = sbmaffinity(internalp, externalp, sizes) + return StochasticBlockModel(sizes, B) end const biclique = ones(2, 2) - Matrix{Float64}(I, 2, 2) # TODO: this documentation needs work. sbromberger 20170326 """ - nearbipartiteaffinity(sizes, between, intra) + nearbipartiteaffinity(sizes, between, intra) Construct the affinity matrix for a near bipartite SBM. `between` is the affinity between the two parts of each bipartite community. @@ -1306,66 +1601,66 @@ Each half is connected as a random bipartite graph with probability `intra` The blocks are connected with probability `between`. """ function nearbipartiteaffinity( - sizes::AbstractVector{T}, between::Real, intra::Real -) where {T<:Integer} - numblocks = div(length(sizes), 2) - return kron(between * Matrix{Float64}(I, numblocks, numblocks), biclique) + - Matrix{Float64}(I, 2 * numblocks, 2 * numblocks) * intra + sizes::AbstractVector{T}, between::Real, intra::Real, +) where {T <: Integer} + numblocks = div(length(sizes), 2) + return kron(between * Matrix{Float64}(I, numblocks, numblocks), biclique) + + Matrix{Float64}(I, 2 * numblocks, 2 * numblocks) * intra end # Return a generator for edges from a stochastic block model near-bipartite graph. function nearbipartiteaffinity( - sizes::Vector{T}, between::Real, inter::Real, noise::Real -) where {T<:Integer} - return nearbipartiteaffinity(sizes, between, inter) .+ noise + sizes::Vector{T}, between::Real, inter::Real, noise::Real, +) where {T <: Integer} + return nearbipartiteaffinity(sizes, between, inter) .+ noise end function nearbipartiteSBM(sizes, between, inter, noise) - return StochasticBlockModel(sizes, nearbipartiteaffinity(sizes, between, inter, noise)) + return StochasticBlockModel(sizes, nearbipartiteaffinity(sizes, between, inter, noise)) end """ - random_pair(rng, n) + random_pair(rng, n) Generate a stream of random pairs in `1:n` using random number generator `RNG`. """ function random_pair(rng::AbstractRNG, n::Integer) - f(ch) = begin - while true - put!(ch, SimpleEdge(rand(rng, 1:n), rand(rng, 1:n))) - end - end - return f + f(ch) = begin + while true + put!(ch, SimpleEdge(rand(rng, 1:n), rand(rng, 1:n))) + end + end + return f end """ - make_edgestream(sbm; rng=nothing, seed=nothing) + make_edgestream(sbm; rng=nothing, seed=nothing) Take an infinite sample from the Stochastic Block Model `sbm`. Pass to `Graph(nvg, neg, edgestream)` to get a Graph object based on `sbm`. """ function make_edgestream( - sbm::StochasticBlockModel; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + sbm::StochasticBlockModel; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - rng = rng_from_rng_or_seed(rng, seed) - pairs = Channel(random_pair(rng, sbm.n); ctype=SimpleEdge, csize=32) - edges(ch) = begin - for e in pairs - i, j = Tuple(e) - i == j && continue - p = sbm.affinities[sbm.nodemap[i], sbm.nodemap[j]] - if rand(rng) < p - put!(ch, e) - end - end - end - return Channel(edges; ctype=SimpleEdge, csize=32) + rng = rng_from_rng_or_seed(rng, seed) + pairs = Channel(random_pair(rng, sbm.n); ctype = SimpleEdge, csize = 32) + edges(ch) = begin + for e in pairs + i, j = Tuple(e) + i == j && continue + p = sbm.affinities[sbm.nodemap[i], sbm.nodemap[j]] + if rand(rng) < p + put!(ch, e) + end + end + end + return Channel(edges; ctype = SimpleEdge, csize = 32) end """ - SimpleGraph{T}(nv, ne, edgestream::Channel) + SimpleGraph{T}(nv, ne, edgestream::Channel) Construct a `SimpleGraph{T}` with `nv` vertices and `ne` edges from `edgestream`. Can result in less than `ne` edges if the channel `edgestream` is closed prematurely. @@ -1373,58 +1668,57 @@ Duplicate edges are only counted once. The element type is the type of `nv`. """ function SimpleGraph(nvg::Integer, neg::Integer, edgestream::Channel) - g = SimpleGraph(nvg) - # println(g) - for e in edgestream - add_edge!(g, e) - ne(g) >= neg && break - end - return g + g = SimpleGraph(nvg) + for e in edgestream + add_edge!(g, e) + ne(g) >= neg && break + end + return g end """ - SimpleGraph{T}(nv, ne, smb::StochasticBlockModel) + SimpleGraph{T}(nv, ne, smb::StochasticBlockModel) Construct a random `SimpleGraph{T}` with `nv` vertices and `ne` edges. The graph is sampled according to the stochastic block model `smb`. The element type is the type of `nv`. """ function SimpleGraph( - nvg::Integer, - neg::Integer, - sbm::StochasticBlockModel; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + nvg::Integer, + neg::Integer, + sbm::StochasticBlockModel; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - return SimpleGraph(nvg, neg, make_edgestream(sbm; rng=rng, seed=seed)) + return SimpleGraph(nvg, neg, make_edgestream(sbm; rng = rng, seed = seed)) end # TODO: this documentation needs work. sbromberger 20170326 """ - blockcounts(sbm, A) + blockcounts(sbm, A) Count the number of edges that go between each block. """ function blockcounts(sbm::StochasticBlockModel, A::AbstractMatrix) - I = collect(1:(sbm.n)) - J = [sbm.nodemap[i] for i in 1:(sbm.n)] - V = ones(sbm.n) - Q = sparse(I, J, V) - return (Q'A) * Q + I = collect(1:(sbm.n)) + J = [sbm.nodemap[i] for i in 1:(sbm.n)] + V = ones(sbm.n) + Q = sparse(I, J, V) + return (Q'A) * Q end function blockcounts(sbm::StochasticBlockModel, g::AbstractGraph) - return blockcounts(sbm, adjacency_matrix(g)) + return blockcounts(sbm, adjacency_matrix(g)) end -function blockfractions(sbm::StochasticBlockModel, g::Union{AbstractGraph,AbstractMatrix}) - bc = blockcounts(sbm, g) - bp = bc ./ sum(bc) - return bp +function blockfractions(sbm::StochasticBlockModel, g::Union{AbstractGraph, AbstractMatrix}) + bc = blockcounts(sbm, g) + bp = bc ./ sum(bc) + return bp end """ - kronecker(SCALE, edgefactor, A=0.57, B=0.19, C=0.19; rng=nothing, seed=nothing) + kronecker(SCALE, edgefactor, A=0.57, B=0.19, C=0.19; rng=nothing, seed=nothing) Generate a directed [Kronecker graph](https://en.wikipedia.org/wiki/Kronecker_graph) with the default Graph500 parameters. @@ -1434,43 +1728,43 @@ References - http://www.graph500.org/specifications#alg:generator """ function kronecker( - SCALE, - edgefactor, - A=0.57, - B=0.19, - C=0.19; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + SCALE, + edgefactor, + A = 0.57, + B = 0.19, + C = 0.19; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - N = 2^SCALE - M = edgefactor * N - ij = ones(Int, M, 2) - ab = A + B - c_norm = C / (1 - (A + B)) - a_norm = A / (A + B) - rng = rng_from_rng_or_seed(rng, seed) - - for ib in 1:SCALE - ii_bit = rand(rng, M) .> (ab) # bitarray - jj_bit = rand(rng, M) .> (c_norm .* (ii_bit) + a_norm .* .!(ii_bit)) - ij .+= 2^(ib - 1) .* (hcat(ii_bit, jj_bit)) - end - - p = randperm(rng, N) - ij = p[ij] - - p = randperm(rng, M) - ij = ij[p, :] - - g = SimpleDiGraph(N) - for (s, d) in zip(@view(ij[:, 1]), @view(ij[:, 2])) - add_edge!(g, s, d) - end - return g + N = 2^SCALE + M = edgefactor * N + ij = ones(Int, M, 2) + ab = A + B + c_norm = C / (1 - (A + B)) + a_norm = A / (A + B) + rng = rng_from_rng_or_seed(rng, seed) + + for ib in 1:SCALE + ii_bit = rand(rng, M) .> (ab) # bitarray + jj_bit = rand(rng, M) .> (c_norm .* (ii_bit) + a_norm .* .!(ii_bit)) + ij .+= 2^(ib - 1) .* (hcat(ii_bit, jj_bit)) + end + + p = randperm(rng, N) + ij = p[ij] + + p = randperm(rng, M) + ij = ij[p, :] + + g = SimpleDiGraph(N) + for (s, d) in zip(@view(ij[:, 1]), @view(ij[:, 2])) + add_edge!(g, s, d) + end + return g end """ - dorogovtsev_mendes(n) + dorogovtsev_mendes(n) Generate a random `n` vertex graph by the Dorogovtsev-Mendes method (with `n \\ge 3`). @@ -1497,37 +1791,37 @@ julia> dorogovtsev_mendes(11, seed=123) ``` """ function dorogovtsev_mendes( - n::Integer; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, + n::Integer; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, ) - n < 3 && throw(DomainError("n=$n must be at least 3")) - rng = rng_from_rng_or_seed(rng, seed) - g = cycle_graph(3) - - for iteration in 1:(n - 3) - chosenedge = rand(rng, 1:(2 * ne(g))) # undirected so each edge is listed twice in adjlist - u, v = -1, -1 - for i in 1:nv(g) - edgelist = outneighbors(g, i) - if chosenedge > length(edgelist) - chosenedge -= length(edgelist) - else - u = i - v = edgelist[chosenedge] - break - end - end - - add_vertex!(g) - add_edge!(g, nv(g), u) - add_edge!(g, nv(g), v) - end - return g + n < 3 && throw(DomainError("n=$n must be at least 3")) + rng = rng_from_rng_or_seed(rng, seed) + g = cycle_graph(3) + + for iteration in 1:(n-3) + chosenedge = rand(rng, 1:(2*ne(g))) # undirected so each edge is listed twice in adjlist + u, v = -1, -1 + for i in 1:nv(g) + edgelist = outneighbors(g, i) + if chosenedge > length(edgelist) + chosenedge -= length(edgelist) + else + u = i + v = edgelist[chosenedge] + break + end + end + + add_vertex!(g) + add_edge!(g, nv(g), u) + add_edge!(g, nv(g), v) + end + return g end """ - random_orientation_dag(g) + random_orientation_dag(g) Generate a random oriented acyclical digraph. The function takes in a simple graph and a random number generator as an argument. The probability of each @@ -1548,20 +1842,20 @@ julia> random_orientation_dag(star_graph(Int8(10)), seed=123) ``` """ function random_orientation_dag( - g::SimpleGraph{T}; - rng::Union{Nothing,AbstractRNG}=nothing, - seed::Union{Nothing,Integer}=nothing, -) where {T<:Integer} - nvg = length(g.fadjlist) - rng = rng_from_rng_or_seed(rng, seed) - order = randperm(rng, nvg) - g2 = SimpleDiGraph(nv(g)) - @inbounds for i in vertices(g) - for j in outneighbors(g, i) - if order[i] < order[j] - add_edge!(g2, i, j) - end - end - end - return g2 + g::SimpleGraph{T}; + rng::Union{Nothing, AbstractRNG} = nothing, + seed::Union{Nothing, Integer} = nothing, +) where {T <: Integer} + nvg = length(g.fadjlist) + rng = rng_from_rng_or_seed(rng, seed) + order = randperm(rng, nvg) + g2 = SimpleDiGraph(nv(g)) + @inbounds for i in vertices(g) + for j in outneighbors(g, i) + if order[i] < order[j] + add_edge!(g2, i, j) + end + end + end + return g2 end From faf90690b692e899119745b65006dd3f16852a7d Mon Sep 17 00:00:00 2001 From: etienneINSA Date: Thu, 15 Feb 2024 16:04:11 +0100 Subject: [PATCH 3/7] This line is actually useful --- src/SimpleGraphs/generators/randgraphs.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/SimpleGraphs/generators/randgraphs.jl b/src/SimpleGraphs/generators/randgraphs.jl index d008c7716..f44bb0963 100644 --- a/src/SimpleGraphs/generators/randgraphs.jl +++ b/src/SimpleGraphs/generators/randgraphs.jl @@ -285,7 +285,6 @@ end buffer_len = 0 if p >= 0.75 @inbounds for v in sample_from - # for v in sample_from if rand(rng) < p buffer_len += 1 buffer[buffer_len] = v @@ -296,9 +295,9 @@ end sample_from_len = length(sample_from) i = 0 @inbounds while true - # while true - i += floor(Int, randexp(rng) * L) + 1 - i > sample_from_len && return buffer_len + s = randexp(rng) * L + i + s > sample_from_len && return buffer_len # compare before ceil to avoid overflow + i += ceil(Int, s) buffer_len += 1 buffer[buffer_len] = sample_from[i] end From 9c3b4a404cd0069702d2a3c596fc11efd2415ffb Mon Sep 17 00:00:00 2001 From: etienneINSA Date: Fri, 3 May 2024 19:11:18 +0200 Subject: [PATCH 4/7] add tests --- src/SimpleGraphs/generators/randgraphs.jl | 1792 ++++++++++---------- test/simplegraphs/generators/randgraphs.jl | 27 + 2 files changed, 927 insertions(+), 892 deletions(-) diff --git a/src/SimpleGraphs/generators/randgraphs.jl b/src/SimpleGraphs/generators/randgraphs.jl index f44bb0963..b93b37055 100644 --- a/src/SimpleGraphs/generators/randgraphs.jl +++ b/src/SimpleGraphs/generators/randgraphs.jl @@ -26,22 +26,22 @@ function _d_sample!(rng, buffer, N, n) while true local X while true - X = Nreal * (- Vprime + 1.0) + X = Nreal * (-Vprime + 1.0) S = trunc(Int, X) S < qu1 && break - Vprime = exp(log(rand(rng))*ninv) + Vprime = exp(log(rand(rng)) * ninv) end U = rand(rng) negSreal = -S - y1 = exp(log(U * Nreal/qu1real) * nmin1inv) + y1 = exp(log(U * Nreal / qu1real) * nmin1inv) Vprime = y1 * (-X / Nreal + 1.0) * (qu1real / (negSreal + qu1real)) Vprime <= 1.0 && break y2 = 1.0 top = -1.0 + Nreal - if - 1 + n > S + if -1 + n > S bottom = -nreal + Nreal - limit = -S + N - else + limit = -S + N + else bottom = -1.0 + negSreal + Nreal limit = qu1 end @@ -50,7 +50,7 @@ function _d_sample!(rng, buffer, N, n) top -= 1.0 bottom -= 1.0 end - if Nreal/(-X + Nreal) >= y1 * exp(log(y2) * nmin1inv) + if Nreal / (-X + Nreal) >= y1 * exp(log(y2) * nmin1inv) Vprime = exp(log(rand(rng)) * nmin1inv) break end @@ -61,7 +61,7 @@ function _d_sample!(rng, buffer, N, n) buffer[current_index] = current_sample current_index += 1 N = -S + (-1 + N) - Nreal = negSreal+ (-1.0 + Nreal) + Nreal = negSreal + (-1.0 + Nreal) n -= 1 nreal += -1.0 ninv = nmin1inv @@ -69,7 +69,7 @@ function _d_sample!(rng, buffer, N, n) qu1real += negSreal threshold += negalphainv end - if n > 1 + if n > 1 # then Use Method A to finish the sampling _a_sample!(rng, buffer, N, n, current_index, current_sample) else # Special case n = 1 @@ -106,7 +106,7 @@ function _a_sample!(rng, buffer, N, n, current_index, current_sample) # Skip over the next S records and select the following one for the sample current_sample += S + 1 buffer[current_index] = current_sample - current_index += 1 + return current_index += 1 end """ @@ -129,68 +129,76 @@ julia> SimpleGraph(5, 7) ``` """ function SimpleGraph{T}( - nv::Integer, - ne::Integer; - self_loops = false, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, -) where {T <: Integer} - maxe = self_loops ? div(Int(nv) * (nv + 1), 2) : div(Int(nv) * (nv - 1), 2) - - @assert(ne <= maxe, "Maximum number of edges for this graph is $maxe") - - remaining_edges(u) = maxe - (self_loops ? div(Int(nv - u) * (nv - u + 1), 2) : div(Int(nv - u) * (nv - u - 1), 2)) - rng = rng_from_rng_or_seed(rng, seed) - - edge_sample = Array{Int}(undef, max(ne, nv)) - _d_sample!(rng, edge_sample, maxe, ne) + nv::Integer, + ne::Integer; + self_loops=false, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, +) where {T<:Integer} + maxe = self_loops ? div(Int(nv) * (nv + 1), 2) : div(Int(nv) * (nv - 1), 2) + + @assert(ne <= maxe, "Maximum number of edges for this graph is $maxe") + + function remaining_edges(u) + return maxe - ( + if self_loops + div(Int(nv - u) * (nv - u + 1), 2) + else + div(Int(nv - u) * (nv - u - 1), 2) + end + ) + end + rng = rng_from_rng_or_seed(rng, seed) + + edge_sample = Array{Int}(undef, max(ne, nv)) + _d_sample!(rng, edge_sample, maxe, ne) fadjlist = Vector{Vector{T}}(undef, nv) - first_index = 1 - current_index = 1 - indeg = zeros(T, nv) - - @inbounds for u in 1:nv - while current_index <= ne && remaining_edges(u) >= edge_sample[current_index] - current_index += 1 - end - indeg_u = indeg[u] - list_u = Vector{T}(undef, indeg_u + current_index-first_index) - for j in first_index:current_index-1 - v = edge_sample[j]-remaining_edges(u-1) + u - if self_loops - v -= 1 - end - list_u[j - first_index + 1 + indeg_u] = v - indeg[v] += 1 - end - fadjlist[u] = list_u - first_index = current_index - end - - insert_at = resize!(edge_sample, 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{T}(ne, fadjlist) + first_index = 1 + current_index = 1 + indeg = zeros(T, nv) + + @inbounds for u in 1:nv + while current_index <= ne && remaining_edges(u) >= edge_sample[current_index] + current_index += 1 + end + indeg_u = indeg[u] + list_u = Vector{T}(undef, indeg_u + current_index - first_index) + for j in first_index:(current_index - 1) + v = edge_sample[j] - remaining_edges(u - 1) + u + if self_loops + v -= 1 + end + list_u[j - first_index + 1 + indeg_u] = v + indeg[v] += 1 + end + fadjlist[u] = list_u + first_index = current_index + end + + insert_at = resize!(edge_sample, 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{T}(ne, fadjlist) end function SimpleGraph( - nv::T, - ne::Integer; - self_loops = false, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, -) where {T <: Integer} - return SimpleGraph{T}(nv, ne; self_loops = self_loops, rng = rng, seed = seed) + nv::T, + ne::Integer; + self_loops=false, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, +) where {T<:Integer} + return SimpleGraph{T}(nv, ne; self_loops=self_loops, rng=rng, seed=seed) end """ @@ -213,41 +221,41 @@ julia> SimpleDiGraph(5, 7) ``` """ function SimpleDiGraph{T}( - nv::Integer, - ne::Integer; - self_loops = false, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, -) where {T <: Integer} - maxd = self_loops ? Int(nv) : nv - 1 - maxe = self_loops ? (Int(nv) * nv) : (Int(nv) * (nv - 1)) - - @assert(ne <= maxe, "Maximum number of edges for this graph is $maxe") - rng = rng_from_rng_or_seed(rng, seed) - - edge_sample = Array{Int}(undef, ne) - _d_sample!(rng, edge_sample, maxe, ne) + nv::Integer, + ne::Integer; + self_loops=false, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, +) where {T<:Integer} + maxd = self_loops ? Int(nv) : nv - 1 + maxe = self_loops ? (Int(nv) * nv) : (Int(nv) * (nv - 1)) + + @assert(ne <= maxe, "Maximum number of edges for this graph is $maxe") + rng = rng_from_rng_or_seed(rng, seed) + + edge_sample = Array{Int}(undef, ne) + _d_sample!(rng, edge_sample, maxe, ne) fadjlist = Vector{Vector{T}}(undef, nv) badjlist = Vector{Vector{T}}(undef, nv) - first_index = 1 - current_index = 1 - - @inbounds for u in 1:nv - while current_index <= ne && maxd*u >= edge_sample[current_index] - current_index += 1 - end - list_u = Vector{T}(undef, current_index-first_index) - for j in first_index:current_index-1 - v = edge_sample[j] - maxd*(u - 1) - if !self_loops && v >= u - v += 1 - end - list_u[j-first_index+1] = v - end - fadjlist[u] = list_u - first_index = current_index - end + first_index = 1 + current_index = 1 + + @inbounds for u in 1:nv + while current_index <= ne && maxd * u >= edge_sample[current_index] + current_index += 1 + end + list_u = Vector{T}(undef, current_index - first_index) + for j in first_index:(current_index - 1) + v = edge_sample[j] - maxd * (u - 1) + if !self_loops && v >= u + v += 1 + end + list_u[j - first_index + 1] = v + end + fadjlist[u] = list_u + first_index = current_index + end indeg = zeros(T, nv) @inbounds for u in OneTo(nv) @@ -258,7 +266,7 @@ function SimpleDiGraph{T}( @inbounds for v in OneTo(nv) badjlist[v] = Vector{T}(undef, indeg[v]) - end + end @inbounds for u in OneTo(nv) for v in fadjlist[u] @@ -270,141 +278,141 @@ function SimpleDiGraph{T}( end function SimpleDiGraph( - nv::T, - ne::Integer; - self_loops = false, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, -) where {T <: Integer} - return SimpleDiGraph{Int}(nv, ne; self_loops = self_loops, rng = rng, seed = seed) + nv::T, + ne::Integer; + self_loops=false, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, +) where {T<:Integer} + return SimpleDiGraph{Int}(nv, ne; self_loops=self_loops, rng=rng, seed=seed) end # like randsubseq!, but without push! since our buffer is always sufficiently large # we might change if randsubseq! gets faster in newer versions of julia @inline function _randsubseq!(rng, buffer, sample_from, p) - buffer_len = 0 - if p >= 0.75 - @inbounds for v in sample_from - if rand(rng) < p - buffer_len += 1 - buffer[buffer_len] = v - end - end - else - L = -1 / log1p(-p) - sample_from_len = length(sample_from) - i = 0 - @inbounds while true - s = randexp(rng) * L - i + s > sample_from_len && return buffer_len # compare before ceil to avoid overflow - i += ceil(Int, s) - buffer_len += 1 - buffer[buffer_len] = sample_from[i] - end - end - return buffer_len + buffer_len = 0 + if p >= 0.75 + @inbounds for v in sample_from + if rand(rng) < p + buffer_len += 1 + buffer[buffer_len] = v + end + end + else + L = -1 / log1p(-p) + sample_from_len = length(sample_from) + i = 0 + @inbounds while true + s = randexp(rng) * L + i + s > sample_from_len && return buffer_len # compare before ceil to avoid overflow + i += ceil(Int, s) + buffer_len += 1 + buffer[buffer_len] = sample_from[i] + end + end + return buffer_len end function _erdos_renyi_directed( - nv::T, - p::Union{Rational, AbstractFloat, AbstractIrrational}; - self_loops = false, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, -) where {T <: Integer} - rng = rng_from_rng_or_seed(rng, seed) - - @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) - buffer_len = _randsubseq!(rng, buffer, sample_from, p) - fadjlist[u] = copy(buffer[1:buffer_len]) - ne += buffer_len - - if !self_loops && u < nv - 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) - 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, 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) + nv::T, + p::Union{Rational,AbstractFloat,AbstractIrrational}; + self_loops=false, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, +) where {T<:Integer} + rng = rng_from_rng_or_seed(rng, seed) + + @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) + buffer_len = _randsubseq!(rng, buffer, sample_from, p) + fadjlist[u] = copy(buffer[1:buffer_len]) + ne += buffer_len + + if !self_loops && u < nv + 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) + 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, 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, - p::Union{Rational, AbstractFloat, AbstractIrrational}; - self_loops = false, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, -) where {T <: Integer} - rng = rng_from_rng_or_seed(rng, seed) + nv::T, + p::Union{Rational,AbstractFloat,AbstractIrrational}; + self_loops=false, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, +) where {T<:Integer} + rng = rng_from_rng_or_seed(rng, seed) - @assert 0 <= p <= 1 + @assert 0 <= p <= 1 - fadjlist = Vector{Vector{T}}(undef, nv) + fadjlist = Vector{Vector{T}}(undef, nv) - ne = 0 - buffer = Vector{T}(undef, nv) - indeg = zeros(T, 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) + @inbounds for u in OneTo(nv) + sample_from = self_loops ? (u:nv) : ((u + one(T)):nv) - buffer_len = _randsubseq!(rng, buffer, sample_from, p) + buffer_len = _randsubseq!(rng, buffer, sample_from, p) - indeg_u = indeg[u] - list_u = Vector{T}(undef, indeg_u + buffer_len) + 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 + for i in OneTo(buffer_len) + v = buffer[i] + list_u[i + indeg_u] = v + indeg[v] += 1 + end - fadjlist[u] = list_u + fadjlist[u] = list_u - ne += buffer_len - end + ne += buffer_len + end - insert_at = resize!(buffer, nv) - fill!(insert_at, one(T)) + 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 + @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) + return SimpleGraph(ne, fadjlist) end """ @@ -435,18 +443,18 @@ julia> erdos_renyi(10, 0.5, is_directed=true, seed=123) ``` """ function erdos_renyi( - n::Integer, - p::Real; - is_directed = false, - self_loops = false, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + n::Integer, + p::Real; + is_directed=false, + self_loops=false, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - if is_directed - return _erdos_renyi_directed(n, p; self_loops = self_loops, rng = rng, seed = seed) - else - return _erdos_renyi_undirected(n, p; self_loops = self_loops, rng = rng, seed = seed) - end + if is_directed + return _erdos_renyi_directed(n, p; self_loops=self_loops, rng=rng, seed=seed) + else + return _erdos_renyi_undirected(n, p; self_loops=self_loops, rng=rng, seed=seed) + end end """ @@ -473,18 +481,18 @@ julia> erdos_renyi(10, 30, is_directed=true, seed=123) ``` """ function erdos_renyi( - n::Integer, - ne::Integer; - is_directed = false, - self_loops = false, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + n::Integer, + ne::Integer; + is_directed=false, + self_loops=false, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - return if is_directed - SimpleDiGraph(n, ne; self_loops = self_loops, rng = rng, seed = seed) - else - SimpleGraph(n, ne; self_loops = self_loops, rng = rng, seed = seed) - end + return if is_directed + SimpleDiGraph(n, ne; self_loops=self_loops, rng=rng, seed=seed) + else + SimpleGraph(n, ne; self_loops=self_loops, rng=rng, seed=seed) + end end """ @@ -525,46 +533,46 @@ julia> print(degree(g)) ``` """ function expected_degree_graph( - ω::Vector{T}; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, -) where {T <: Real} - g = SimpleGraph(length(ω)) - return expected_degree_graph!(g, ω; rng = rng, seed = seed) + ω::Vector{T}; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, +) where {T<:Real} + g = SimpleGraph(length(ω)) + return expected_degree_graph!(g, ω; rng=rng, seed=seed) end function expected_degree_graph!( - g::SimpleGraph, - ω::Vector{T}; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, -) where {T <: Real} - n = length(ω) - @assert all(zero(T) .<= ω .<= n - one(T)) "Elements of ω needs to be at least 0 and at most n-1" - - π = sortperm(ω; rev = true) - rng = rng_from_rng_or_seed(rng, seed) - - S = sum(ω) - - for u in 1:(n-1) - v = u + 1 - p = min(ω[π[u]] * ω[π[v]] / S, one(T)) - while v <= n && p > zero(p) - if p != one(T) - v += floor(Int, log(rand(rng)) / log(one(T) - p)) - end - if v <= n - q = min(ω[π[u]] * ω[π[v]] / S, one(T)) - if rand(rng) < q / p - add_edge!(g, π[u], π[v]) - end - p = q - v += 1 - end - end - end - return g + g::SimpleGraph, + ω::Vector{T}; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, +) where {T<:Real} + n = length(ω) + @assert all(zero(T) .<= ω .<= n - one(T)) "Elements of ω needs to be at least 0 and at most n-1" + + π = sortperm(ω; rev=true) + rng = rng_from_rng_or_seed(rng, seed) + + S = sum(ω) + + for u in 1:(n - 1) + v = u + 1 + p = min(ω[π[u]] * ω[π[v]] / S, one(T)) + while v <= n && p > zero(p) + if p != one(T) + v += floor(Int, log(rand(rng)) / log(one(T) - p)) + end + if v <= n + q = min(ω[π[u]] * ω[π[v]] / S, one(T)) + if rand(rng) < q / p + add_edge!(g, π[u], π[v]) + end + p = q + v += 1 + end + end + end + return g end """ @@ -612,62 +620,62 @@ julia> watts_strogatz(Int8(10), 4, 0.8, is_directed=true, seed=123) - Small Worlds, Duncan J. watts. [https://en.wikipedia.org/wiki/Special:BookSources?isbn=978-0691005416](https://en.wikipedia.org/wiki/Special:BookSources?isbn=978-0691005416) """ function watts_strogatz( - n::Integer, - k::Integer, - β::Real; - is_directed::Bool = false, - remove_edges::Bool = true, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + n::Integer, + k::Integer, + β::Real; + is_directed::Bool=false, + remove_edges::Bool=true, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - @assert k < n - - # If we have n - 1 neighbors (exactly k/2 on each side), then the graph is - # necessarily complete. No need to run the Watts-Strogatz procedure: - if k == n - 1 && iseven(k) - return is_directed ? complete_digraph(n) : complete_graph(n) - end - - g = is_directed ? SimpleDiGraph(n) : SimpleGraph(n) - - # The ith next vertex, in clockwise order. - # (Reduce to zero-based indexing, so the modulo works, by subtracting 1 - # before and adding 1 after.) - @inline target(s, i) = ((s + i - 1) % n) + 1 - - # Phase 1: For each step size i, add an edge from each vertex s to the ith - # next vertex, in clockwise order. - - for i in 1:div(k, 2), s in 1:n - add_edge!(g, s, target(s, i)) - end - - # Phase 2: For each step size i and each vertex s, consider the edge to the - # ith next vertex, in clockwise order. With probability β, delete the edge - # and rewire it to any (valid) target, chosen uniformly at random. - - rng = rng_from_rng_or_seed(rng, seed) - for i in 1:div(k, 2), s in 1:n - - # We only rewire with a probability β, and we only worry about rewiring - # if there is some vertex not connected to s; otherwise, the only valid - # rewiring is to reconnect to the ith next vertex, and there is no work - # to do. - (rand(rng) < β && degree(g, s) < n - 1) || continue - - t = target(s, i) - - while true - d = rand(rng, 1:n) # Tentative new target - d == s && continue # Self-loops prohibited - d == t && break # Rewired to original target - if add_edge!(g, s, d) # Was this valid (i.e., unconnected)? - remove_edges && rem_edge!(g, s, t) # True rewiring: Delete original edge - break # We found a valid target - end - end - end - return g + @assert k < n + + # If we have n - 1 neighbors (exactly k/2 on each side), then the graph is + # necessarily complete. No need to run the Watts-Strogatz procedure: + if k == n - 1 && iseven(k) + return is_directed ? complete_digraph(n) : complete_graph(n) + end + + g = is_directed ? SimpleDiGraph(n) : SimpleGraph(n) + + # The ith next vertex, in clockwise order. + # (Reduce to zero-based indexing, so the modulo works, by subtracting 1 + # before and adding 1 after.) + @inline target(s, i) = ((s + i - 1) % n) + 1 + + # Phase 1: For each step size i, add an edge from each vertex s to the ith + # next vertex, in clockwise order. + + for i in 1:div(k, 2), s in 1:n + add_edge!(g, s, target(s, i)) + end + + # Phase 2: For each step size i and each vertex s, consider the edge to the + # ith next vertex, in clockwise order. With probability β, delete the edge + # and rewire it to any (valid) target, chosen uniformly at random. + + rng = rng_from_rng_or_seed(rng, seed) + for i in 1:div(k, 2), s in 1:n + + # We only rewire with a probability β, and we only worry about rewiring + # if there is some vertex not connected to s; otherwise, the only valid + # rewiring is to reconnect to the ith next vertex, and there is no work + # to do. + (rand(rng) < β && degree(g, s) < n - 1) || continue + + t = target(s, i) + + while true + d = rand(rng, 1:n) # Tentative new target + d == s && continue # Self-loops prohibited + d == t && break # Rewired to original target + if add_edge!(g, s, d) # Was this valid (i.e., unconnected)? + remove_edges && rem_edge!(g, s, t) # True rewiring: Delete original edge + break # We found a valid target + end + end + end + return g end """ @@ -718,69 +726,69 @@ julia> newman_watts_strogatz(Int8(10), 4, 0.8, is_directed=true, seed=123) - Scaling and percolation in the small-world network model, M. E. J. Newman, Duncan J. Watts. [https://doi.org/10.1103/PhysRevE.60.7332](https://doi.org/10.1103/PhysRevE.60.7332) """ function newman_watts_strogatz( - n::Integer, - k::Integer, - β::Real; - is_directed::Bool = false, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + n::Integer, + k::Integer, + β::Real; + is_directed::Bool=false, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - return watts_strogatz( - n, k, β; is_directed = is_directed, remove_edges = false, rng = rng, seed = seed, - ) + return watts_strogatz( + n, k, β; is_directed=is_directed, remove_edges=false, rng=rng, seed=seed + ) end -function _suitable(edges::Set{SimpleEdge{T}}, potential_edges::Dict{T, T}) where {T <: Integer} - isempty(potential_edges) && return true - list = keys(potential_edges) - for s1 in list, s2 in list - s1 >= s2 && continue - (SimpleEdge(s1, s2) ∉ edges) && return true - end - return false +function _suitable(edges::Set{SimpleEdge{T}}, potential_edges::Dict{T,T}) where {T<:Integer} + isempty(potential_edges) && return true + list = keys(potential_edges) + for s1 in list, s2 in list + s1 >= s2 && continue + (SimpleEdge(s1, s2) ∉ edges) && return true + end + return false end _try_creation(n::Integer, k::Integer, rng::AbstractRNG) = _try_creation(n, fill(k, n), rng) -function _try_creation(n::T, k::Vector{T}, rng::AbstractRNG) where {T <: Integer} - edges = Set{SimpleEdge{T}}() - m = 0 - stubs = zeros(T, sum(k)) - for i in one(T):n - for j in one(T):k[i] - m += 1 - stubs[m] = i - end - end - # stubs = vcat([fill(i, k[i]) for i = 1:n]...) # slower - - while !isempty(stubs) - potential_edges = Dict{T, T}() - shuffle!(rng, stubs) - for i in 1:2:length(stubs) - s1, s2 = stubs[i:(i+1)] - if (s1 > s2) - s1, s2 = s2, s1 - end - e = SimpleEdge(s1, s2) - if s1 != s2 && ∉(e, edges) - push!(edges, e) - else - potential_edges[s1] = get(potential_edges, s1, 0) + 1 - potential_edges[s2] = get(potential_edges, s2, 0) + 1 - end - end - - if !_suitable(edges, potential_edges) - return Set{SimpleEdge{T}}() - end - - stubs = Vector{T}() - for (e, ct) in potential_edges - append!(stubs, fill(e, ct)) - end - end - return edges +function _try_creation(n::T, k::Vector{T}, rng::AbstractRNG) where {T<:Integer} + edges = Set{SimpleEdge{T}}() + m = 0 + stubs = zeros(T, sum(k)) + for i in one(T):n + for j in one(T):k[i] + m += 1 + stubs[m] = i + end + end + # stubs = vcat([fill(i, k[i]) for i = 1:n]...) # slower + + while !isempty(stubs) + potential_edges = Dict{T,T}() + shuffle!(rng, stubs) + for i in 1:2:length(stubs) + s1, s2 = stubs[i:(i + 1)] + if (s1 > s2) + s1, s2 = s2, s1 + end + e = SimpleEdge(s1, s2) + if s1 != s2 && ∉(e, edges) + push!(edges, e) + else + potential_edges[s1] = get(potential_edges, s1, 0) + 1 + potential_edges[s2] = get(potential_edges, s2, 0) + 1 + end + end + + if !_suitable(edges, potential_edges) + return Set{SimpleEdge{T}}() + end + + stubs = Vector{T}() + for (e, ct) in potential_edges + append!(stubs, fill(e, ct)) + end + end + return edges end """ @@ -837,22 +845,22 @@ julia> barabasi_albert(100, Int8(10), 3, is_directed=true, seed=123) ``` """ function barabasi_albert( - n::Integer, - n0::Integer, - k::Integer; - is_directed::Bool = false, - complete::Bool = false, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + n::Integer, + n0::Integer, + k::Integer; + is_directed::Bool=false, + complete::Bool=false, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - if complete - g = is_directed ? complete_digraph(n0) : complete_graph(n0) - else - g = is_directed ? SimpleDiGraph(n0) : SimpleGraph(n0) - end - - barabasi_albert!(g, n, k; rng = rng, seed = seed) - return g + if complete + g = is_directed ? complete_digraph(n0) : complete_graph(n0) + else + g = is_directed ? SimpleDiGraph(n0) : SimpleGraph(n0) + end + + barabasi_albert!(g, n, k; rng=rng, seed=seed) + return g end """ @@ -880,74 +888,74 @@ julia> g ``` """ function barabasi_albert!( - g::AbstractGraph, - n::Integer, - k::Integer; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + g::AbstractGraph, + n::Integer, + k::Integer; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - n0 = nv(g) - 1 <= k <= n0 <= n || - throw(ArgumentError("Barabási-Albert model requires 1 <= k <= nv(g) <= n")) - n0 == n && return g - - # seed random number generator - rng = rng_from_rng_or_seed(rng, seed) - - # add missing vertices - sizehint!(g.fadjlist, n) - add_vertices!(g, n - n0) - - # if initial graph doesn't contain any edges - # expand it by one vertex and add k edges from this additional node - if ne(g) == 0 - # expand initial graph - n0 += one(n0) - - # add edges to k existing vertices - for target in sample!(rng, collect(1:(n0-1)), k) - add_edge!(g, n0, target) - end - end - - # vector of weighted vertices (each node is repeated once for each adjacent edge) - weightedVs = Vector{Int}(undef, 2 * (n - n0) * k + 2 * ne(g)) - - # initialize vector of weighted vertices - offset = 0 - for e in edges(g) - weightedVs[offset+=1] = src(e) - weightedVs[offset+=1] = dst(e) - end - - # array to record if a node is picked - picked = fill(false, n) - - # vector of targets - targets = Vector{Int}(undef, k) - - for source in (n0+1):n - # choose k targets from the existing vertices - # pick uniformly from weightedVs (preferential attachment) - i = 0 - while i < k - target = weightedVs[rand(rng, 1:offset)] - if !picked[target] - targets[i+=1] = target - picked[target] = true - end - end - - # add edges to k targets - for target in targets - add_edge!(g, source, target) - - weightedVs[offset+=1] = source - weightedVs[offset+=1] = target - picked[target] = false - end - end - return g + n0 = nv(g) + 1 <= k <= n0 <= n || + throw(ArgumentError("Barabási-Albert model requires 1 <= k <= nv(g) <= n")) + n0 == n && return g + + # seed random number generator + rng = rng_from_rng_or_seed(rng, seed) + + # add missing vertices + sizehint!(g.fadjlist, n) + add_vertices!(g, n - n0) + + # if initial graph doesn't contain any edges + # expand it by one vertex and add k edges from this additional node + if ne(g) == 0 + # expand initial graph + n0 += one(n0) + + # add edges to k existing vertices + for target in sample!(rng, collect(1:(n0 - 1)), k) + add_edge!(g, n0, target) + end + end + + # vector of weighted vertices (each node is repeated once for each adjacent edge) + weightedVs = Vector{Int}(undef, 2 * (n - n0) * k + 2 * ne(g)) + + # initialize vector of weighted vertices + offset = 0 + for e in edges(g) + weightedVs[offset += 1] = src(e) + weightedVs[offset += 1] = dst(e) + end + + # array to record if a node is picked + picked = fill(false, n) + + # vector of targets + targets = Vector{Int}(undef, k) + + for source in (n0 + 1):n + # choose k targets from the existing vertices + # pick uniformly from weightedVs (preferential attachment) + i = 0 + while i < k + target = weightedVs[rand(rng, 1:offset)] + if !picked[target] + targets[i += 1] = target + picked[target] = true + end + end + + # add edges to k targets + for target in targets + add_edge!(g, source, target) + + weightedVs[offset += 1] = source + weightedVs[offset += 1] = target + picked[target] = false + end + end + return g end """ @@ -996,29 +1004,29 @@ julia> edges(g) |> collect ``` """ function static_fitness_model( - m::Integer, - fitness::Vector{T}; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, -) where {T <: Real} - m < 0 && throw(ArgumentError("number of edges must be positive")) - n = length(fitness) - m == 0 && return SimpleGraph(n) - nvs = 0 - for f in fitness - # sanity check for the fitness - f < zero(T) && throw(ArgumentError("fitness scores must be non-negative")) - f > zero(T) && (nvs += 1) - end - # avoid getting into an infinite loop when too many edges are requested - max_no_of_edges = div(nvs * (nvs - 1), 2) - m > max_no_of_edges && - throw(ArgumentError("too many edges requested ($m > $max_no_of_edges)")) - # calculate the cumulative fitness scores - cum_fitness = cumsum(fitness) - g = SimpleGraph(n) - _create_static_fitness_graph!(g, m, cum_fitness, cum_fitness, rng, seed) - return g + m::Integer, + fitness::Vector{T}; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, +) where {T<:Real} + m < 0 && throw(ArgumentError("number of edges must be positive")) + n = length(fitness) + m == 0 && return SimpleGraph(n) + nvs = 0 + for f in fitness + # sanity check for the fitness + f < zero(T) && throw(ArgumentError("fitness scores must be non-negative")) + f > zero(T) && (nvs += 1) + end + # avoid getting into an infinite loop when too many edges are requested + max_no_of_edges = div(nvs * (nvs - 1), 2) + m > max_no_of_edges && + throw(ArgumentError("too many edges requested ($m > $max_no_of_edges)")) + # calculate the cumulative fitness scores + cum_fitness = cumsum(fitness) + g = SimpleGraph(n) + _create_static_fitness_graph!(g, m, cum_fitness, cum_fitness, rng, seed) + return g end """ @@ -1056,58 +1064,58 @@ julia> edges(g) |> collect ``` """ function static_fitness_model( - m::Integer, - fitness_out::Vector{T}, - fitness_in::Vector{S}; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, -) where {T <: Real} where {S <: Real} - m < 0 && throw(ArgumentError("number of edges must be positive")) - n = length(fitness_out) - length(fitness_in) != n && - throw(ArgumentError("fitness_in must have the same size as fitness_out")) - m == 0 && return SimpleDiGraph(n) - # avoid getting into an infinite loop when too many edges are requested - noutvs = ninvs = nvs = 0 - @inbounds for i in 1:n - # sanity check for the fitness - (fitness_out[i] < zero(T) || fitness_in[i] < zero(S)) && - error("fitness scores must be non-negative") # TODO 0.7: change to DomainError? - fitness_out[i] > zero(T) && (noutvs += 1) - fitness_in[i] > zero(S) && (ninvs += 1) - (fitness_out[i] > zero(T) && fitness_in[i] > zero(S)) && (nvs += 1) - end - max_no_of_edges = noutvs * ninvs - nvs - m > max_no_of_edges && - throw(ArgumentError("too many edges requested ($m > $max_no_of_edges)")) - # calculate the cumulative fitness scores - cum_fitness_out = cumsum(fitness_out) - cum_fitness_in = cumsum(fitness_in) - g = SimpleDiGraph(n) - _create_static_fitness_graph!(g, m, cum_fitness_out, cum_fitness_in, rng, seed) - return g + m::Integer, + fitness_out::Vector{T}, + fitness_in::Vector{S}; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, +) where {T<:Real} where {S<:Real} + m < 0 && throw(ArgumentError("number of edges must be positive")) + n = length(fitness_out) + length(fitness_in) != n && + throw(ArgumentError("fitness_in must have the same size as fitness_out")) + m == 0 && return SimpleDiGraph(n) + # avoid getting into an infinite loop when too many edges are requested + noutvs = ninvs = nvs = 0 + @inbounds for i in 1:n + # sanity check for the fitness + (fitness_out[i] < zero(T) || fitness_in[i] < zero(S)) && + error("fitness scores must be non-negative") # TODO 0.7: change to DomainError? + fitness_out[i] > zero(T) && (noutvs += 1) + fitness_in[i] > zero(S) && (ninvs += 1) + (fitness_out[i] > zero(T) && fitness_in[i] > zero(S)) && (nvs += 1) + end + max_no_of_edges = noutvs * ninvs - nvs + m > max_no_of_edges && + throw(ArgumentError("too many edges requested ($m > $max_no_of_edges)")) + # calculate the cumulative fitness scores + cum_fitness_out = cumsum(fitness_out) + cum_fitness_in = cumsum(fitness_in) + g = SimpleDiGraph(n) + _create_static_fitness_graph!(g, m, cum_fitness_out, cum_fitness_in, rng, seed) + return g end function _create_static_fitness_graph!( - g::AbstractGraph, - m::Integer, - cum_fitness_out::Vector{T}, - cum_fitness_in::Vector{S}, - rng::Union{Nothing, AbstractRNG}, - seed::Union{Nothing, Integer}, -) where {T <: Real} where {S <: Real} - max_out = cum_fitness_out[end] - max_in = cum_fitness_in[end] - rng = rng_from_rng_or_seed(rng, seed) - while m > 0 - source = searchsortedfirst(cum_fitness_out, rand(rng) * max_out) - target = searchsortedfirst(cum_fitness_in, rand(rng) * max_in) - # skip if loop edge - (source == target) && continue - # is there already an edge? If so, try again - add_edge!(g, source, target) || continue - m -= one(m) - end + g::AbstractGraph, + m::Integer, + cum_fitness_out::Vector{T}, + cum_fitness_in::Vector{S}, + rng::Union{Nothing,AbstractRNG}, + seed::Union{Nothing,Integer}, +) where {T<:Real} where {S<:Real} + max_out = cum_fitness_out[end] + max_in = cum_fitness_in[end] + rng = rng_from_rng_or_seed(rng, seed) + while m > 0 + source = searchsortedfirst(cum_fitness_out, rand(rng) * max_out) + target = searchsortedfirst(cum_fitness_in, rand(rng) * max_in) + # skip if loop edge + (source == target) && continue + # is there already an edge? If so, try again + add_edge!(g, source, target) || continue + m -= one(m) + end end """ @@ -1131,17 +1139,17 @@ Time complexity is ``\\mathcal{O}(|V| + |E| log |E|)``. - Cho YS, Kim JS, Park J, Kahng B, Kim D: Percolation transitions in scale-free networks under the Achlioptas process. Phys Rev Lett 103:135702, 2009. """ function static_scale_free( - n::Integer, - m::Integer, - α::Real; - finite_size_correction::Bool = true, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + n::Integer, + m::Integer, + α::Real; + finite_size_correction::Bool=true, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - n < 0 && throw(ArgumentError("number of vertices must be positive")) - α < 2 && throw(ArgumentError("out-degree exponent must be >= 2")) - fitness = _construct_fitness(n, α, finite_size_correction) - return static_fitness_model(m, fitness; rng = rng, seed = seed) + n < 0 && throw(ArgumentError("number of vertices must be positive")) + α < 2 && throw(ArgumentError("out-degree exponent must be >= 2")) + fitness = _construct_fitness(n, α, finite_size_correction) + return static_fitness_model(m, fitness; rng=rng, seed=seed) end """ @@ -1166,39 +1174,39 @@ Time complexity is ``\\mathcal{O}(|V| + |E| log |E|)``. - Cho YS, Kim JS, Park J, Kahng B, Kim D: Percolation transitions in scale-free networks under the Achlioptas process. Phys Rev Lett 103:135702, 2009. """ function static_scale_free( - n::Integer, - m::Integer, - α_out::Real, - α_in::Float64; - finite_size_correction::Bool = true, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + n::Integer, + m::Integer, + α_out::Real, + α_in::Float64; + finite_size_correction::Bool=true, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - n < 0 && throw(ArgumentError("number of vertices must be positive")) - α_out < 2 && throw(ArgumentError("out-degree exponent must be >= 2")) - α_in < 2 && throw(ArgumentError("out-degree exponent must be >= 2")) - # construct the fitness - fitness_out = _construct_fitness(n, α_out, finite_size_correction) - fitness_in = _construct_fitness(n, α_in, finite_size_correction) - # eliminate correlation - shuffle!(fitness_in) - return static_fitness_model(m, fitness_out, fitness_in; rng = rng, seed = seed) + n < 0 && throw(ArgumentError("number of vertices must be positive")) + α_out < 2 && throw(ArgumentError("out-degree exponent must be >= 2")) + α_in < 2 && throw(ArgumentError("out-degree exponent must be >= 2")) + # construct the fitness + fitness_out = _construct_fitness(n, α_out, finite_size_correction) + fitness_in = _construct_fitness(n, α_in, finite_size_correction) + # eliminate correlation + shuffle!(fitness_in) + return static_fitness_model(m, fitness_out, fitness_in; rng=rng, seed=seed) end function _construct_fitness(n::Integer, α::Real, finite_size_correction::Bool) - α = -1 / (α - 1) - fitness = zeros(n) - j = float(n) - if finite_size_correction && α < -0.5 - # See the Cho et al paper, first page first column + footnote 7 - j += n^(1 + 1 / 2 * α) * (10 * sqrt(2) * (1 + α))^(-1 / α) - 1 - end - j = max(j, n) - @inbounds for i in 1:n - fitness[i] = j^α - j -= 1 - end - return fitness + α = -1 / (α - 1) + fitness = zeros(n) + j = float(n) + if finite_size_correction && α < -0.5 + # See the Cho et al paper, first page first column + footnote 7 + j += n^(1 + 1 / 2 * α) * (10 * sqrt(2) * (1 + α))^(-1 / α) - 1 + end + j = max(j, n) + @inbounds for i in 1:n + fitness[i] = j^α + j -= 1 + end + return fitness end """ @@ -1220,32 +1228,32 @@ Allocates an array of `nk` `Int`s, and . For ``k > \\frac{n}{2}``, generates a g ``n-k-1`` and returns its complement. """ function random_regular_graph( - n::Integer, - k::Integer; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + n::Integer, + k::Integer; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - !iseven(n * k) && throw(ArgumentError("n * k must be even")) - !(0 <= k < n) && throw(ArgumentError("the 0 <= k < n inequality must be satisfied")) - if k == 0 - return SimpleGraph(n) - end - rng = rng_from_rng_or_seed(rng, seed) - if (k > n / 2) && iseven(n * (n - k - 1)) - return complement(random_regular_graph(n, n - k - 1; rng = rng)) - end - - edges = _try_creation(n, k, rng) - while isempty(edges) - edges = _try_creation(n, k, rng) - end - - g = SimpleGraph(n) - for edge in edges - add_edge!(g, edge) - end - - return g + !iseven(n * k) && throw(ArgumentError("n * k must be even")) + !(0 <= k < n) && throw(ArgumentError("the 0 <= k < n inequality must be satisfied")) + if k == 0 + return SimpleGraph(n) + end + rng = rng_from_rng_or_seed(rng, seed) + if (k > n / 2) && iseven(n * (n - k - 1)) + return complement(random_regular_graph(n, n - k - 1; rng=rng)) + end + + edges = _try_creation(n, k, rng) + while isempty(edges) + edges = _try_creation(n, k, rng) + end + + g = SimpleGraph(n) + for edge in edges + add_edge!(g, edge) + end + + return g end """ @@ -1267,32 +1275,32 @@ Time complexity is approximately ``\\mathcal{O}(n \\bar{k}^2)``. Allocates an array of ``n \\bar{k}`` `Int`s. """ function random_configuration_model( - n::Integer, - k::Array{T}; - check_graphical::Bool = false, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, -) where {T <: Integer} - n != length(k) && throw(ArgumentError("a degree sequence of length n must be provided")) - m = sum(k) - !iseven(m) && throw(ArgumentError("sum(k) must be even")) - !all(0 .<= k .< n) && - throw(ArgumentError("the 0 <= k[i] < n inequality must be satisfied")) - if check_graphical - isgraphical(k) || throw(ArgumentError("degree sequence must be graphical")) - end - rng = rng_from_rng_or_seed(rng, seed) - - edges = _try_creation(n, k, rng) - while m > 0 && isempty(edges) - edges = _try_creation(n, k, rng) - end - - g = SimpleGraphFromIterator(edges) - if nv(g) < n - add_vertices!(g, n - nv(g)) - end - return g + n::Integer, + k::Array{T}; + check_graphical::Bool=false, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, +) where {T<:Integer} + n != length(k) && throw(ArgumentError("a degree sequence of length n must be provided")) + m = sum(k) + !iseven(m) && throw(ArgumentError("sum(k) must be even")) + !all(0 .<= k .< n) && + throw(ArgumentError("the 0 <= k[i] < n inequality must be satisfied")) + if check_graphical + isgraphical(k) || throw(ArgumentError("degree sequence must be graphical")) + end + rng = rng_from_rng_or_seed(rng, seed) + + edges = _try_creation(n, k, rng) + while m > 0 && isempty(edges) + edges = _try_creation(n, k, rng) + end + + g = SimpleGraphFromIterator(edges) + if nv(g) < n + add_vertices!(g, n - nv(g)) + end + return g end """ uniform_tree(n) @@ -1310,12 +1318,12 @@ julia> uniform_tree(10) {10, 9} undirected simple Int64 graph ``` """ -function uniform_tree(n::Integer; rng::Union{Nothing, AbstractRNG} = nothing) - n <= 1 && return Graph(n) - n == 2 && return path_graph(n) - rng = rng_from_rng_or_seed(rng, nothing) - random_code = rand(rng, Base.OneTo(n), n - 2) - return prufer_decode(random_code) +function uniform_tree(n::Integer; rng::Union{Nothing,AbstractRNG}=nothing) + n <= 1 && return Graph(n) + n == 2 && return path_graph(n) + rng = rng_from_rng_or_seed(rng, nothing) + random_code = rand(rng, Base.OneTo(n), n - 2) + return prufer_decode(random_code) end """ @@ -1334,36 +1342,36 @@ Allocates an ``n × n`` sparse matrix of boolean as an adjacency matrix and uses that to generate the directed graph. """ function random_regular_digraph( - n::Integer, - k::Integer; - dir::Symbol = :out, - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + n::Integer, + k::Integer; + dir::Symbol=:out, + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - # TODO remove the function sample from StatsBase for one allowing the use - # of a local rng - !(0 <= k < n) && throw(ArgumentError("the 0 <= k < n inequality must be satisfied")) - - if k == 0 - return SimpleDiGraph(n) - end - rng = rng_from_rng_or_seed(rng, seed) - if (k > n / 2) && iseven(n * (n - k - 1)) - return complement(random_regular_digraph(n, n - k - 1; dir = dir, rng = rng)) - end - cs = collect(2:n) - i = 1 - I = Vector{Int}(undef, n * k) - J = Vector{Int}(undef, n * k) - V = fill(true, n * k) - for r in 1:n - l = ((r-1)*k+1):(r*k) - I[l] .= r - J[l] = sample!(rng, cs, k; exclude = r) - end - - m = dir == :out ? sparse(I, J, V, n, n) : sparse(I, J, V, n, n)' - return SimpleDiGraph(m) + # TODO remove the function sample from StatsBase for one allowing the use + # of a local rng + !(0 <= k < n) && throw(ArgumentError("the 0 <= k < n inequality must be satisfied")) + + if k == 0 + return SimpleDiGraph(n) + end + rng = rng_from_rng_or_seed(rng, seed) + if (k > n / 2) && iseven(n * (n - k - 1)) + return complement(random_regular_digraph(n, n - k - 1; dir=dir, rng=rng)) + end + cs = collect(2:n) + i = 1 + I = Vector{Int}(undef, n * k) + J = Vector{Int}(undef, n * k) + V = fill(true, n * k) + for r in 1:n + l = ((r - 1) * k + 1):(r * k) + I[l] .= r + J[l] = sample!(rng, cs, k; exclude=r) + end + + m = dir == :out ? sparse(I, J, V, n, n) : sparse(I, J, V, n, n)' + return SimpleDiGraph(m) end """ @@ -1389,18 +1397,18 @@ julia> random_tournament_digraph(Int8(10), seed=123) ``` """ function random_tournament_digraph( - n::Integer; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + n::Integer; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - rng = rng_from_rng_or_seed(rng, seed) - g = SimpleDiGraph(n) + rng = rng_from_rng_or_seed(rng, seed) + g = SimpleDiGraph(n) - for i in 1:n, j in (i+1):n - rand(rng, Bool) ? add_edge!(g, SimpleEdge(i, j)) : add_edge!(g, SimpleEdge(j, i)) - end + for i in 1:n, j in (i + 1):n + rand(rng, Bool) ? add_edge!(g, SimpleEdge(i, j)) : add_edge!(g, SimpleEdge(j, i)) + end - return g + return g end """ @@ -1413,21 +1421,21 @@ Return a binomially-distributed random number with parameters `n` and `p` and op - http://stackoverflow.com/questions/23561551/a-efficient-binomial-random-number-generator-code-in-java """ function randbn( - n::Integer, - p::Real; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + n::Integer, + p::Real; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - rng = rng_from_rng_or_seed(rng, seed) - log_q = log(1.0 - p) - x = 0 - sum = 0.0 - for _ in 1:n - sum += log(rand(rng)) / (n - x) - sum < log_q && break - x += 1 - end - return x + rng = rng_from_rng_or_seed(rng, seed) + log_q = log(1.0 - p) + x = 0 + sum = 0.0 + for _ in 1:n + sum += log(rand(rng)) / (n - x) + sum < log_q && break + x += 1 + end + return x end """ @@ -1448,46 +1456,46 @@ For a dynamic version of the SBM see the [`StochasticBlockModel`](@ref) type and related functions. """ function stochastic_block_model( - c::Matrix{T}, - n::Vector{U}; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, -) where {T <: Real} where {U <: Integer} - size(c, 1) == size(c, 2) == length(n) || - throw(ArgumentError("matrix-vector size mismatch")) - - # init dsfmt generator with a fixed seed - rng = rng_from_rng_or_seed(rng, seed) - N = sum(n) - K = length(n) - nedg = zeros(Int, K, K) - g = SimpleGraph(N) - cum = [sum(n[1:a]) for a in 0:K] - for a in 1:K - ra = (cum[a]+1):cum[a+1] - for b in a:K - ((a == b) && !(c[a, b] <= n[b] - 1)) || - ((a != b) && !(c[a, b] <= n[b])) && error( - "Mean degree cannot be greater than available neighbors in the block.", - ) # TODO 0.7: turn into some other error? - - m = a == b ? div(n[a] * (n[a] - 1), 2) : n[a] * n[b] - p = a == b ? n[a] * c[a, b] / (2m) : n[a] * c[a, b] / m - nedg = randbn(m, p; rng = rng) - rb = (cum[b]+1):cum[b+1] - i = 0 - while i < nedg - source = rand(rng, ra) - dest = rand(rng, rb) - if source != dest - if add_edge!(g, source, dest) - i += 1 - end - end - end - end - end - return g + c::Matrix{T}, + n::Vector{U}; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, +) where {T<:Real} where {U<:Integer} + size(c, 1) == size(c, 2) == length(n) || + throw(ArgumentError("matrix-vector size mismatch")) + + # init dsfmt generator with a fixed seed + rng = rng_from_rng_or_seed(rng, seed) + N = sum(n) + K = length(n) + nedg = zeros(Int, K, K) + g = SimpleGraph(N) + cum = [sum(n[1:a]) for a in 0:K] + for a in 1:K + ra = (cum[a] + 1):cum[a + 1] + for b in a:K + ((a == b) && !(c[a, b] <= n[b] - 1)) || + ((a != b) && !(c[a, b] <= n[b])) && error( + "Mean degree cannot be greater than available neighbors in the block.", + ) # TODO 0.7: turn into some other error? + + m = a == b ? div(n[a] * (n[a] - 1), 2) : n[a] * n[b] + p = a == b ? n[a] * c[a, b] / (2m) : n[a] * c[a, b] / m + nedg = randbn(m, p; rng=rng) + rb = (cum[b] + 1):cum[b + 1] + i = 0 + while i < nedg + source = rand(rng, ra) + dest = rand(rng, rb) + if source != dest + if add_edge!(g, source, dest) + i += 1 + end + end + end + end + end + return g end """ @@ -1497,15 +1505,15 @@ Return a Graph generated according to the Stochastic Block Model (SBM), sampling from an SBM with ``c_{a,a}=cint``, and ``c_{a,b}=cext``. """ function stochastic_block_model( - cint::T, - cext::T, - n::Vector{U}; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, -) where {T <: Real} where {U <: Integer} - K = length(n) - c = [ifelse(a == b, cint, cext) for a in 1:K, b in 1:K] - return stochastic_block_model(c, n; rng = rng, seed = seed) + cint::T, + cext::T, + n::Vector{U}; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, +) where {T<:Real} where {U<:Integer} + K = length(n) + c = [ifelse(a == b, cint, cext) for a in 1:K, b in 1:K] + return stochastic_block_model(c, n; rng=rng, seed=seed) end """ @@ -1525,32 +1533,32 @@ block `k` and any vertex in block `l`. Graphs are generated by taking random ``i,j ∈ V`` and flipping a coin with probability `affinities[nodemap[i],nodemap[j]]`. """ -mutable struct StochasticBlockModel{T <: Integer, P <: Real} - n::T - nodemap::Array{T} - affinities::Matrix{P} +mutable struct StochasticBlockModel{T<:Integer,P<:Real} + n::T + nodemap::Array{T} + affinities::Matrix{P} end function ==(sbm::StochasticBlockModel, other::StochasticBlockModel) - return (sbm.n == other.n) && - (sbm.nodemap == other.nodemap) && - (sbm.affinities == other.affinities) + return (sbm.n == other.n) && + (sbm.nodemap == other.nodemap) && + (sbm.affinities == other.affinities) end # A constructor for StochasticBlockModel that uses the sizes of the blocks # and the affinity matrix. This construction implies that consecutive # vertices will be in the same blocks, except for the block boundaries. function StochasticBlockModel(sizes::AbstractVector, affinities::AbstractMatrix) - csum = cumsum(sizes) - j = 1 - nodemap = zeros(Int, csum[end]) - for i in 1:csum[end] - if i > csum[j] - j += 1 - end - nodemap[i] = j - end - return StochasticBlockModel(csum[end], nodemap, affinities) + csum = cumsum(sizes) + j = 1 + nodemap = zeros(Int, csum[end]) + for i in 1:csum[end] + if i > csum[j] + j += 1 + end + nodemap[i] = j + end + return StochasticBlockModel(csum[end], nodemap, affinities) end ### TODO: This documentation needs work. sbromberger 20170326 @@ -1561,28 +1569,28 @@ Produce the sbm affinity matrix with internal probabilities `internalp` and external probabilities `externalp`. """ function sbmaffinity( - internalp::Vector{T}, externalp::Real, sizes::Vector{U}, -) where {T <: Real} where {U <: Integer} - numblocks = length(sizes) - numblocks == length(internalp) || - throw(ArgumentError("Inconsistent input dimensions: internalp, sizes")) - B = diagm(0 => internalp) + externalp * (ones(numblocks, numblocks) - I) - return B + internalp::Vector{T}, externalp::Real, sizes::Vector{U} +) where {T<:Real} where {U<:Integer} + numblocks = length(sizes) + numblocks == length(internalp) || + throw(ArgumentError("Inconsistent input dimensions: internalp, sizes")) + B = diagm(0 => internalp) + externalp * (ones(numblocks, numblocks) - I) + return B end function StochasticBlockModel( - internalp::Real, externalp::Real, size::Integer, numblocks::Integer, + internalp::Real, externalp::Real, size::Integer, numblocks::Integer ) - sizes = fill(size, numblocks) - B = sbmaffinity(fill(internalp, numblocks), externalp, sizes) - return StochasticBlockModel(sizes, B) + sizes = fill(size, numblocks) + B = sbmaffinity(fill(internalp, numblocks), externalp, sizes) + return StochasticBlockModel(sizes, B) end function StochasticBlockModel( - internalp::Vector{T}, externalp::Real, sizes::Vector{U}, -) where {T <: Real, U <: Integer} - B = sbmaffinity(internalp, externalp, sizes) - return StochasticBlockModel(sizes, B) + internalp::Vector{T}, externalp::Real, sizes::Vector{U} +) where {T<:Real,U<:Integer} + B = sbmaffinity(internalp, externalp, sizes) + return StochasticBlockModel(sizes, B) end const biclique = ones(2, 2) - Matrix{Float64}(I, 2, 2) @@ -1600,22 +1608,22 @@ Each half is connected as a random bipartite graph with probability `intra` The blocks are connected with probability `between`. """ function nearbipartiteaffinity( - sizes::AbstractVector{T}, between::Real, intra::Real, -) where {T <: Integer} - numblocks = div(length(sizes), 2) - return kron(between * Matrix{Float64}(I, numblocks, numblocks), biclique) + - Matrix{Float64}(I, 2 * numblocks, 2 * numblocks) * intra + sizes::AbstractVector{T}, between::Real, intra::Real +) where {T<:Integer} + numblocks = div(length(sizes), 2) + return kron(between * Matrix{Float64}(I, numblocks, numblocks), biclique) + + Matrix{Float64}(I, 2 * numblocks, 2 * numblocks) * intra end # Return a generator for edges from a stochastic block model near-bipartite graph. function nearbipartiteaffinity( - sizes::Vector{T}, between::Real, inter::Real, noise::Real, -) where {T <: Integer} - return nearbipartiteaffinity(sizes, between, inter) .+ noise + sizes::Vector{T}, between::Real, inter::Real, noise::Real +) where {T<:Integer} + return nearbipartiteaffinity(sizes, between, inter) .+ noise end function nearbipartiteSBM(sizes, between, inter, noise) - return StochasticBlockModel(sizes, nearbipartiteaffinity(sizes, between, inter, noise)) + return StochasticBlockModel(sizes, nearbipartiteaffinity(sizes, between, inter, noise)) end """ @@ -1624,12 +1632,12 @@ end Generate a stream of random pairs in `1:n` using random number generator `RNG`. """ function random_pair(rng::AbstractRNG, n::Integer) - f(ch) = begin - while true - put!(ch, SimpleEdge(rand(rng, 1:n), rand(rng, 1:n))) - end - end - return f + f(ch) = begin + while true + put!(ch, SimpleEdge(rand(rng, 1:n), rand(rng, 1:n))) + end + end + return f end """ @@ -1639,23 +1647,23 @@ Take an infinite sample from the Stochastic Block Model `sbm`. Pass to `Graph(nvg, neg, edgestream)` to get a Graph object based on `sbm`. """ function make_edgestream( - sbm::StochasticBlockModel; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + sbm::StochasticBlockModel; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - rng = rng_from_rng_or_seed(rng, seed) - pairs = Channel(random_pair(rng, sbm.n); ctype = SimpleEdge, csize = 32) - edges(ch) = begin - for e in pairs - i, j = Tuple(e) - i == j && continue - p = sbm.affinities[sbm.nodemap[i], sbm.nodemap[j]] - if rand(rng) < p - put!(ch, e) - end - end - end - return Channel(edges; ctype = SimpleEdge, csize = 32) + rng = rng_from_rng_or_seed(rng, seed) + pairs = Channel(random_pair(rng, sbm.n); ctype=SimpleEdge, csize=32) + edges(ch) = begin + for e in pairs + i, j = Tuple(e) + i == j && continue + p = sbm.affinities[sbm.nodemap[i], sbm.nodemap[j]] + if rand(rng) < p + put!(ch, e) + end + end + end + return Channel(edges; ctype=SimpleEdge, csize=32) end """ @@ -1667,12 +1675,12 @@ Duplicate edges are only counted once. The element type is the type of `nv`. """ function SimpleGraph(nvg::Integer, neg::Integer, edgestream::Channel) - g = SimpleGraph(nvg) - for e in edgestream - add_edge!(g, e) - ne(g) >= neg && break - end - return g + g = SimpleGraph(nvg) + for e in edgestream + add_edge!(g, e) + ne(g) >= neg && break + end + return g end """ @@ -1683,13 +1691,13 @@ The graph is sampled according to the stochastic block model `smb`. The element type is the type of `nv`. """ function SimpleGraph( - nvg::Integer, - neg::Integer, - sbm::StochasticBlockModel; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + nvg::Integer, + neg::Integer, + sbm::StochasticBlockModel; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - return SimpleGraph(nvg, neg, make_edgestream(sbm; rng = rng, seed = seed)) + return SimpleGraph(nvg, neg, make_edgestream(sbm; rng=rng, seed=seed)) end # TODO: this documentation needs work. sbromberger 20170326 @@ -1699,21 +1707,21 @@ end Count the number of edges that go between each block. """ function blockcounts(sbm::StochasticBlockModel, A::AbstractMatrix) - I = collect(1:(sbm.n)) - J = [sbm.nodemap[i] for i in 1:(sbm.n)] - V = ones(sbm.n) - Q = sparse(I, J, V) - return (Q'A) * Q + I = collect(1:(sbm.n)) + J = [sbm.nodemap[i] for i in 1:(sbm.n)] + V = ones(sbm.n) + Q = sparse(I, J, V) + return (Q'A) * Q end function blockcounts(sbm::StochasticBlockModel, g::AbstractGraph) - return blockcounts(sbm, adjacency_matrix(g)) + return blockcounts(sbm, adjacency_matrix(g)) end -function blockfractions(sbm::StochasticBlockModel, g::Union{AbstractGraph, AbstractMatrix}) - bc = blockcounts(sbm, g) - bp = bc ./ sum(bc) - return bp +function blockfractions(sbm::StochasticBlockModel, g::Union{AbstractGraph,AbstractMatrix}) + bc = blockcounts(sbm, g) + bp = bc ./ sum(bc) + return bp end """ @@ -1727,39 +1735,39 @@ References - http://www.graph500.org/specifications#alg:generator """ function kronecker( - SCALE, - edgefactor, - A = 0.57, - B = 0.19, - C = 0.19; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + SCALE, + edgefactor, + A=0.57, + B=0.19, + C=0.19; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - N = 2^SCALE - M = edgefactor * N - ij = ones(Int, M, 2) - ab = A + B - c_norm = C / (1 - (A + B)) - a_norm = A / (A + B) - rng = rng_from_rng_or_seed(rng, seed) - - for ib in 1:SCALE - ii_bit = rand(rng, M) .> (ab) # bitarray - jj_bit = rand(rng, M) .> (c_norm .* (ii_bit) + a_norm .* .!(ii_bit)) - ij .+= 2^(ib - 1) .* (hcat(ii_bit, jj_bit)) - end - - p = randperm(rng, N) - ij = p[ij] - - p = randperm(rng, M) - ij = ij[p, :] - - g = SimpleDiGraph(N) - for (s, d) in zip(@view(ij[:, 1]), @view(ij[:, 2])) - add_edge!(g, s, d) - end - return g + N = 2^SCALE + M = edgefactor * N + ij = ones(Int, M, 2) + ab = A + B + c_norm = C / (1 - (A + B)) + a_norm = A / (A + B) + rng = rng_from_rng_or_seed(rng, seed) + + for ib in 1:SCALE + ii_bit = rand(rng, M) .> (ab) # bitarray + jj_bit = rand(rng, M) .> (c_norm .* (ii_bit) + a_norm .* .!(ii_bit)) + ij .+= 2^(ib - 1) .* (hcat(ii_bit, jj_bit)) + end + + p = randperm(rng, N) + ij = p[ij] + + p = randperm(rng, M) + ij = ij[p, :] + + g = SimpleDiGraph(N) + for (s, d) in zip(@view(ij[:, 1]), @view(ij[:, 2])) + add_edge!(g, s, d) + end + return g end """ @@ -1790,33 +1798,33 @@ julia> dorogovtsev_mendes(11, seed=123) ``` """ function dorogovtsev_mendes( - n::Integer; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, + n::Integer; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, ) - n < 3 && throw(DomainError("n=$n must be at least 3")) - rng = rng_from_rng_or_seed(rng, seed) - g = cycle_graph(3) - - for iteration in 1:(n-3) - chosenedge = rand(rng, 1:(2*ne(g))) # undirected so each edge is listed twice in adjlist - u, v = -1, -1 - for i in 1:nv(g) - edgelist = outneighbors(g, i) - if chosenedge > length(edgelist) - chosenedge -= length(edgelist) - else - u = i - v = edgelist[chosenedge] - break - end - end - - add_vertex!(g) - add_edge!(g, nv(g), u) - add_edge!(g, nv(g), v) - end - return g + n < 3 && throw(DomainError("n=$n must be at least 3")) + rng = rng_from_rng_or_seed(rng, seed) + g = cycle_graph(3) + + for iteration in 1:(n - 3) + chosenedge = rand(rng, 1:(2 * ne(g))) # undirected so each edge is listed twice in adjlist + u, v = -1, -1 + for i in 1:nv(g) + edgelist = outneighbors(g, i) + if chosenedge > length(edgelist) + chosenedge -= length(edgelist) + else + u = i + v = edgelist[chosenedge] + break + end + end + + add_vertex!(g) + add_edge!(g, nv(g), u) + add_edge!(g, nv(g), v) + end + return g end """ @@ -1841,20 +1849,20 @@ julia> random_orientation_dag(star_graph(Int8(10)), seed=123) ``` """ function random_orientation_dag( - g::SimpleGraph{T}; - rng::Union{Nothing, AbstractRNG} = nothing, - seed::Union{Nothing, Integer} = nothing, -) where {T <: Integer} - nvg = length(g.fadjlist) - rng = rng_from_rng_or_seed(rng, seed) - order = randperm(rng, nvg) - g2 = SimpleDiGraph(nv(g)) - @inbounds for i in vertices(g) - for j in outneighbors(g, i) - if order[i] < order[j] - add_edge!(g2, i, j) - end - end - end - return g2 + g::SimpleGraph{T}; + rng::Union{Nothing,AbstractRNG}=nothing, + seed::Union{Nothing,Integer}=nothing, +) where {T<:Integer} + nvg = length(g.fadjlist) + rng = rng_from_rng_or_seed(rng, seed) + order = randperm(rng, nvg) + g2 = SimpleDiGraph(nv(g)) + @inbounds for i in vertices(g) + for j in outneighbors(g, i) + if order[i] < order[j] + add_edge!(g2, i, j) + end + end + end + return g2 end diff --git a/test/simplegraphs/generators/randgraphs.jl b/test/simplegraphs/generators/randgraphs.jl index 8004b7cc8..2113c7eb5 100644 --- a/test/simplegraphs/generators/randgraphs.jl +++ b/test/simplegraphs/generators/randgraphs.jl @@ -4,12 +4,30 @@ @testset "(Int, Int)" begin r1 = SimpleGraph(10, 20; rng=rng) r2 = SimpleDiGraph(5, 10; rng=rng) + r3 = SimpleGraph(10, 20; self_loops=true, rng=rng) + r4 = SimpleDiGraph(5, 10; self_loops=true, rng=rng) + @test isvalid_simplegraph(r1) + @test isvalid_simplegraph(r2) + @test isvalid_simplegraph(r3) + @test isvalid_simplegraph(r4) + @test has_self_loops(r1) == false + @test has_self_loops(r2) == false + @test is_directed(r1) == false + @test is_directed(r2) == true + @test is_directed(r3) == false + @test is_directed(r4) == true @test nv(r1) == 10 @test ne(r1) == 20 @test nv(r2) == 5 @test ne(r2) == 10 + @test nv(r3) == 10 + @test ne(r3) == 20 + @test nv(r4) == 5 + @test ne(r4) == 10 @test eltype(r1) == Int @test eltype(r2) == Int + @test eltype(r3) == Int + @test eltype(r4) == Int @test SimpleGraph(10, 20; rng=StableRNG(3)) == SimpleGraph(10, 20; rng=StableRNG(3)) @test SimpleGraph(10, 40; rng=StableRNG(3)) == SimpleGraph(10, 40; rng=StableRNG(3)) @@ -40,8 +58,17 @@ er = erdos_renyi(10, 0.5; rng=rng) @test nv(er) == 10 @test is_directed(er) == false + @test has_self_loops(er) == false er = erdos_renyi(10, 0.5; is_directed=true, rng=rng) @test nv(er) == 10 + @test has_self_loops(er) == false + @test is_directed(er) == true + + er = erdos_renyi(10, 0.5; has_self_loops=true, rng=rng) + @test nv(er) == 10 + @test is_directed(er) == false + er = erdos_renyi(10, 0.5; is_directed=true, has_self_loops=true, rng=rng) + @test nv(er) == 10 @test is_directed(er) == true er = erdos_renyi(10, 0.5; rng=StableRNG(17)) From 505634582dd380f0223062f7a4f68612007baddd Mon Sep 17 00:00:00 2001 From: etienneINSA Date: Mon, 6 May 2024 10:48:24 +0200 Subject: [PATCH 5/7] fix tests --- src/SimpleGraphs/generators/randgraphs.jl | 17 ++++++++++----- src/SimpleGraphs/simpledigraph.jl | 4 ++-- src/SimpleGraphs/simplegraph.jl | 4 ++-- test/simplegraphs/generators/randgraphs.jl | 24 ++++++++++++++++++++-- 4 files changed, 38 insertions(+), 11 deletions(-) diff --git a/src/SimpleGraphs/generators/randgraphs.jl b/src/SimpleGraphs/generators/randgraphs.jl index b93b37055..426afe8db 100644 --- a/src/SimpleGraphs/generators/randgraphs.jl +++ b/src/SimpleGraphs/generators/randgraphs.jl @@ -135,6 +135,9 @@ function SimpleGraph{T}( rng::Union{Nothing,AbstractRNG}=nothing, seed::Union{Nothing,Integer}=nothing, ) where {T<:Integer} + + ne == 0 && return Graph{T}(nv) + maxe = self_loops ? div(Int(nv) * (nv + 1), 2) : div(Int(nv) * (nv - 1), 2) @assert(ne <= maxe, "Maximum number of edges for this graph is $maxe") @@ -227,6 +230,9 @@ function SimpleDiGraph{T}( rng::Union{Nothing,AbstractRNG}=nothing, seed::Union{Nothing,Integer}=nothing, ) where {T<:Integer} + + ne == 0 && return DiGraph{T}(nv) + maxd = self_loops ? Int(nv) : nv - 1 maxe = self_loops ? (Int(nv) * nv) : (Int(nv) * (nv - 1)) @@ -322,8 +328,6 @@ function _erdos_renyi_directed( ) where {T<:Integer} rng = rng_from_rng_or_seed(rng, seed) - @assert 0 <= p <= 1 - fadjlist = Vector{Vector{T}}(undef, nv) ne = 0 @@ -371,10 +375,9 @@ function _erdos_renyi_undirected( rng::Union{Nothing,AbstractRNG}=nothing, seed::Union{Nothing,Integer}=nothing, ) where {T<:Integer} + rng = rng_from_rng_or_seed(rng, seed) - @assert 0 <= p <= 1 - fadjlist = Vector{Vector{T}}(undef, nv) ne = 0 @@ -439,7 +442,7 @@ julia> erdos_renyi(10, 0.5) julia> using Graphs julia> erdos_renyi(10, 0.5, is_directed=true, seed=123) -{10, 49} directed simple Int64 graph +{10, 47} directed simple Int64 graph ``` """ function erdos_renyi( @@ -450,6 +453,10 @@ function erdos_renyi( rng::Union{Nothing,AbstractRNG}=nothing, seed::Union{Nothing,Integer}=nothing, ) + p >= 1 && return is_directed ? complete_digraph(n) : complete_graph(n) + + @assert p >= 0 + if is_directed return _erdos_renyi_directed(n, p; self_loops=self_loops, rng=rng, seed=seed) else diff --git a/src/SimpleGraphs/simpledigraph.jl b/src/SimpleGraphs/simpledigraph.jl index 4b7605a1a..27a845aee 100644 --- a/src/SimpleGraphs/simpledigraph.jl +++ b/src/SimpleGraphs/simpledigraph.jl @@ -11,7 +11,7 @@ mutable struct SimpleDiGraph{T<:Integer} <: AbstractSimpleGraph{T} badjlist::Vector{Vector{T}} # [dst]: (src, src, src) function SimpleDiGraph{T}( - ne::Int, fadjlist::Vector{Vector{T}}, badjlist::Vector{Vector{T}} + ne::Integer, fadjlist::Vector{Vector{T}}, badjlist::Vector{Vector{T}} ) where {T} throw_if_invalid_eltype(T) return new(ne, fadjlist, badjlist) @@ -19,7 +19,7 @@ mutable struct SimpleDiGraph{T<:Integer} <: AbstractSimpleGraph{T} end function SimpleDiGraph( - ne::Int, fadjlist::Vector{Vector{T}}, badjlist::Vector{Vector{T}} + ne::Integer, fadjlist::Vector{Vector{T}}, badjlist::Vector{Vector{T}} ) where {T} return SimpleDiGraph{T}(ne, fadjlist, badjlist) end diff --git a/src/SimpleGraphs/simplegraph.jl b/src/SimpleGraphs/simplegraph.jl index b896c33ae..cb022051c 100644 --- a/src/SimpleGraphs/simplegraph.jl +++ b/src/SimpleGraphs/simplegraph.jl @@ -9,13 +9,13 @@ mutable struct SimpleGraph{T<:Integer} <: AbstractSimpleGraph{T} ne::Int fadjlist::Vector{Vector{T}} # [src]: (dst, dst, dst) - function SimpleGraph{T}(ne::Int, fadjlist::Vector{Vector{T}}) where {T} + function SimpleGraph{T}(ne::Integer, fadjlist::Vector{Vector{T}}) where {T} throw_if_invalid_eltype(T) return new{T}(ne, fadjlist) end end -function SimpleGraph(ne, fadjlist::Vector{Vector{T}}) where {T} +function SimpleGraph(ne::Integer, fadjlist::Vector{Vector{T}}) where {T} return SimpleGraph{T}(ne, fadjlist) end diff --git a/test/simplegraphs/generators/randgraphs.jl b/test/simplegraphs/generators/randgraphs.jl index 2113c7eb5..f0bce87d7 100644 --- a/test/simplegraphs/generators/randgraphs.jl +++ b/test/simplegraphs/generators/randgraphs.jl @@ -29,6 +29,16 @@ @test eltype(r3) == Int @test eltype(r4) == Int + er = Graph(0, 0; rng=rng) + @test er == Graph() + er = DiGraph(0, 0; rng=rng) + @test er == DiGraph() + + er = Graph(3, 0; rng=rng) + @test er == Graph(3) + er = DiGraph(3, 0; rng=rng) + @test er == DiGraph(3) + @test SimpleGraph(10, 20; rng=StableRNG(3)) == SimpleGraph(10, 20; rng=StableRNG(3)) @test SimpleGraph(10, 40; rng=StableRNG(3)) == SimpleGraph(10, 40; rng=StableRNG(3)) @test SimpleDiGraph(10, 20; rng=StableRNG(3)) == @@ -64,10 +74,10 @@ @test has_self_loops(er) == false @test is_directed(er) == true - er = erdos_renyi(10, 0.5; has_self_loops=true, rng=rng) + er = erdos_renyi(10, 0.5; self_loops=true, rng=rng) @test nv(er) == 10 @test is_directed(er) == false - er = erdos_renyi(10, 0.5; is_directed=true, has_self_loops=true, rng=rng) + er = erdos_renyi(10, 0.5; is_directed=true, self_loops=true, rng=rng) @test nv(er) == 10 @test is_directed(er) == true @@ -80,6 +90,16 @@ @test erdos_renyi(5, 2.1; rng=rng) == complete_graph(5) @test erdos_renyi(5, 2.1; is_directed=true, rng=rng) == complete_digraph(5) + er = erdos_renyi(0, 0.0; rng=rng) + @test er == Graph() + er = erdos_renyi(0, 0.0; is_directed=true, rng=rng) + @test er == DiGraph() + + er = erdos_renyi(3, 0.0; rng=rng) + @test er == Graph(3) + er = erdos_renyi(3, 0.0; is_directed=true, rng=rng) + @test er == DiGraph(3) + # issue #173 er = erdos_renyi(4, 6; seed=1) @test nv(er) == 4 From 28e6c92006baaecf2b43472f94b023a5ba93aa54 Mon Sep 17 00:00:00 2001 From: etienneINSA Date: Mon, 6 May 2024 10:48:55 +0200 Subject: [PATCH 6/7] use randexp --- src/SimpleGraphs/generators/randgraphs.jl | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/src/SimpleGraphs/generators/randgraphs.jl b/src/SimpleGraphs/generators/randgraphs.jl index 426afe8db..bc4299e6c 100644 --- a/src/SimpleGraphs/generators/randgraphs.jl +++ b/src/SimpleGraphs/generators/randgraphs.jl @@ -15,7 +15,7 @@ function _d_sample!(rng, buffer, N, n) nreal = n ninv = 1.0 / nreal Nreal = N - Vprime = exp(log(rand(rng)) * ninv) + Vprime = exp(-randexp(rng) * ninv) qu1 = -n + 1 + N qu1real = -nreal + 1.0 + Nreal negalphainv = -13 @@ -29,7 +29,7 @@ function _d_sample!(rng, buffer, N, n) X = Nreal * (-Vprime + 1.0) S = trunc(Int, X) S < qu1 && break - Vprime = exp(log(rand(rng)) * ninv) + Vprime = exp(-randexp(rng) * ninv) end U = rand(rng) negSreal = -S @@ -51,10 +51,10 @@ function _d_sample!(rng, buffer, N, n) bottom -= 1.0 end if Nreal / (-X + Nreal) >= y1 * exp(log(y2) * nmin1inv) - Vprime = exp(log(rand(rng)) * nmin1inv) + Vprime = exp(-randexp(rng) * nmin1inv) break end - Vprime = exp(log(rand(rng)) * ninv) + Vprime = exp(-randexp(rng) * ninv) end # Skip over the next S records and select the following one for the sample current_sample += S + 1 @@ -135,7 +135,6 @@ function SimpleGraph{T}( rng::Union{Nothing,AbstractRNG}=nothing, seed::Union{Nothing,Integer}=nothing, ) where {T<:Integer} - ne == 0 && return Graph{T}(nv) maxe = self_loops ? div(Int(nv) * (nv + 1), 2) : div(Int(nv) * (nv - 1), 2) @@ -230,7 +229,6 @@ function SimpleDiGraph{T}( rng::Union{Nothing,AbstractRNG}=nothing, seed::Union{Nothing,Integer}=nothing, ) where {T<:Integer} - ne == 0 && return DiGraph{T}(nv) maxd = self_loops ? Int(nv) : nv - 1 @@ -375,7 +373,6 @@ function _erdos_renyi_undirected( rng::Union{Nothing,AbstractRNG}=nothing, seed::Union{Nothing,Integer}=nothing, ) where {T<:Integer} - rng = rng_from_rng_or_seed(rng, seed) fadjlist = Vector{Vector{T}}(undef, nv) @@ -454,7 +451,7 @@ function erdos_renyi( seed::Union{Nothing,Integer}=nothing, ) p >= 1 && return is_directed ? complete_digraph(n) : complete_graph(n) - + @assert p >= 0 if is_directed @@ -567,7 +564,7 @@ function expected_degree_graph!( p = min(ω[π[u]] * ω[π[v]] / S, one(T)) while v <= n && p > zero(p) if p != one(T) - v += floor(Int, log(rand(rng)) / log(one(T) - p)) + v += floor(Int, -randexp(rng) / log(one(T) - p)) end if v <= n q = min(ω[π[u]] * ω[π[v]] / S, one(T)) @@ -1438,7 +1435,7 @@ function randbn( x = 0 sum = 0.0 for _ in 1:n - sum += log(rand(rng)) / (n - x) + sum += -randexp(rng) / (n - x) sum < log_q && break x += 1 end From 6e22b60d7835010de17520743438d4454e4d3f99 Mon Sep 17 00:00:00 2001 From: etienneINSA Date: Mon, 6 May 2024 12:08:38 +0200 Subject: [PATCH 7/7] try to please codecov --- test/simplegraphs/generators/randgraphs.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/simplegraphs/generators/randgraphs.jl b/test/simplegraphs/generators/randgraphs.jl index de74a672b..7f6ca0c10 100644 --- a/test/simplegraphs/generators/randgraphs.jl +++ b/test/simplegraphs/generators/randgraphs.jl @@ -66,18 +66,22 @@ @testset "Erdös-Renyí" begin er = erdos_renyi(10, 0.5; rng=rng) + @test isvalid_simplegraph(er) @test nv(er) == 10 @test is_directed(er) == false @test has_self_loops(er) == false er = erdos_renyi(10, 0.5; is_directed=true, rng=rng) + @test isvalid_simplegraph(er) @test nv(er) == 10 @test has_self_loops(er) == false @test is_directed(er) == true er = erdos_renyi(10, 0.5; self_loops=true, rng=rng) + @test isvalid_simplegraph(er) @test nv(er) == 10 @test is_directed(er) == false er = erdos_renyi(10, 0.5; is_directed=true, self_loops=true, rng=rng) + @test isvalid_simplegraph(er) @test nv(er) == 10 @test is_directed(er) == true @@ -85,6 +89,10 @@ @test nv(er) == 10 @test is_directed(er) == false + er = erdos_renyi(10, 0.8; rng=StableRNG(17)) + @test nv(er) == 10 + @test is_directed(er) == false + @test erdos_renyi(5, 1.0; rng=rng) == complete_graph(5) @test erdos_renyi(5, 1.0; is_directed=true, rng=rng) == complete_digraph(5) @test erdos_renyi(5, 2.1; rng=rng) == complete_graph(5)