BiCM

Maximum entropy model for the Undirected Bipartite Configuration Model (BiCM).

The object holds the maximum likelihood parameters of the model (θ), the expected bi-adjacency matrix (Ĝ), and the variance for the elements of the adjacency matrix (σ).

BiCM(G::T; precision::N=Float64, kwargs...) where {T<:Graphs.AbstractGraph, N<:Real}
BiCM(;d⊥::Vector{T}, d⊤::Vector{T}, precision::Type{<:AbstractFloat}=Float64, kwargs...)

Constructor function for the BiCM type. The graph you provide should be bipartite

By default and dependng on the graph type T, the definition of degree from $Graphs.jl$ is applied. If you want to use a different definition of degrees, you can pass vectors of degrees sequences as keyword arguments (d⊥, d⊤). If you want to generate a model directly from degree sequences without an underlying graph , you can simply pass the degree sequences as arguments (d⊥, d⊤). If you want to work from an adjacency matrix, or edge list, you can use the graph constructors from the $JuliaGraphs$ ecosystem.

Zero degree nodes have a zero probability of being connected to other nodes, so they are skipped in the computation of the model.

Examples

# generating a model from a graph
julia> G = corporateclub();

julia> model =  BiCM(G)
BiCM{Graphs.SimpleGraphs.SimpleGraph{Int64}, Float64} (25 + 15 vertices, 6 + 6 unique degrees, 0.30 compression ratio)

# generating a model directly from a degree sequence
julia> model = model = BiCM(d⊥=[1,1,2,2,2,3,3,1,1,2], d⊤=[3,4,5,2,5,6,6,1,1,2])
BiCM{Nothing, Float64} (10 + 10 vertices, 3 + 6 unique degrees, 0.45 compression ratio)

# generating a model directly from a degree sequence with a different precision
julia> model = model = BiCM(d⊥=[1,1,2,2,2,3,3,1,1,2], d⊤=[3,4,5,2,5,6,6,1,1,2], precision=Float32)
BiCM{Nothing, Float32} (10 + 10 vertices, 3 + 6 unique degrees, 0.45 compression ratio)

# generating a model from an adjacency matrix
julia> A = [0 0 0 1 0;0 0 0 1 0;0 0 0 0 1;1 1 0 0 0;0 0 1 0 0];

julia> G = MaxEntropyGraphs.Graphs.SimpleGraph(A);

julia> @assert MaxEntropyGraphs.Graphs.is_bipartite(G); # check if the graph is bipartite

julia> model = BiCM(G) # generating the model
BiCM{Graphs.SimpleGraphs.SimpleGraph{Int64}, Float64} (3 + 2 vertices, 1 + 2 unique degrees, 0.60 compression ratio)

# generating a model from a biadjacency matrix
julia> biadjacency = [1 0;1 0; 0 1];

julia> N⊥,N⊤ = size(biadjacency); # layer dimensions

julia> adjacency = [zeros(Int, N⊥,N⊥) biadjacency; biadjacency' zeros(Int,N⊤,N⊤)];

julia> G = MaxEntropyGraphs.Graphs.SimpleGraph(adjacency); # generate graph

julia> model = BiCM(G) # generate model
BiCM{Graphs.SimpleGraphs.SimpleGraph{Int64}, Float64} (3 + 2 vertices, 1 + 2 unique degrees, 0.60 compression ratio)

# generating a model from an edge list
julia> edges = MaxEntropyGraphs.Graphs.SimpleEdge.([(1,4); (2,4); (3,5)]);

julia> G = MaxEntropyGraphs.Graphs.SimpleGraph(edges); # generate graph

julia> model = BiCM(G) # generate model
BiCM{Graphs.SimpleGraphs.SimpleGraph{Int64}, Float64} (3 + 2 vertices, 1 + 2 unique degrees, 0.60 compression ratio)

solve_model!(m::BiCM)

Compute the likelihood maximising parameters of the BiCM model m.

Arguments

• method::Symbol: solution method to use, can be :fixedpoint (default), or :NelderMead, :BFGS, :LBFGS and :Newton.
• initial::Symbol: initial guess for the parameters $\Theta$, can be :degrees (default), :random, :uniform, or :chung_lu.
• maxiters::Int: maximum number of iterations for the solver (defaults to 1000).
• verbose::Bool: set to show log messages (defaults to false).
• ftol::Real: function tolerance for convergence with the fixedpoint method (defaults to 1e-8).
• abstol::Union{Number, Nothing}: absolute function tolerance for convergence with the other methods (defaults to nothing).
• reltol::Union{Number, Nothing}: relative function tolerance for convergence with the other methods (defaults to nothing).
• AD_method::Symbol: autodiff method to use, can be any of :AutoZygote, :AutoReverseDiff, :AutoForwardDiff and :AutoFiniteDiff. Performance depends on the size of the problem (defaults to :AutoZygote),
• analytical_gradient::Bool: set the use the analytical gradient instead of the one generated with autodiff (defaults to false)

Examples

# default use
julia> model = BiCM(corporateclub());

julia> solve_model!(model);

# using analytical gradient and uniform initial guess
julia> solve_model!(model, method=:BFGS, analytical_gradient=true, initial=:uniform)
(BiCM{Graphs.SimpleGraphs.SimpleGraph{Int64}, Float64} (25 + 15 vertices, 6 + 6 unique degrees, 0.30 compression ratio), retcode: Success
u: [1.449571644621672, 0.8231752829683303, 0.34755085972479766, -0.04834480708852856, -0.3984299800917503, -0.7223268299919358, 1.6090554004279671, 1.2614196476197532, 0.9762560461922147, 0.11406188481061938, -0.24499004480426345, -2.2646067641037333]
Final objective value:     171.15095803718134
)

initial_guess(m::BiCM; method::Symbol=:degrees)

Compute an initial guess for the maximum likelihood parameters of the BiCM model m using the method method.

The methods available are:

• :degrees (default): the initial guess is computed using the degrees of the graph, i.e. $\theta = [-\log(d_{ot}); -\log(d_{ op})]$
• :random: the initial guess is computed using random values between 0 and 1, i.e. $\theta_{i} = -\log(r_{i})$ where $r_{i} \sim U(0,1)$
• :uniform: the initial guess is uniformily set to 0.5, i.e. $\theta_{i} = -\log(0.5)$
• :chung_lu: the initial guess is computed using the degrees of the graph and the number of edges, i.e. $\theta = [-\log(d_{ot}/(2E)); -\log(d_{ op}/(2E))]$

Examples

julia> model = BiCM(corporateclub());

julia> initial_guess(model, method=:random);

julia> initial_guess(model, method=:uniform);

julia> initial_guess(model, method=:chung_lu);

julia> initial_guess(model);

rand(m::BiCM; precomputed::Bool=false)

Generate a random graph from the BiCM model m.

Arguments:

• precomputed::Bool: if true, the precomputed expected biadjacency matrix (m.Ĝ) is used to generate the random graph, otherwise the maximum likelihood parameters are used to generate the random graph on the fly. For larger networks, it is recommended to not precompute the expected adjacency matrix to limit memory pressure.

Note: The generated graph will also be bipartite and respect the layer membership of the original graph used to define the model.

rand(m::BiCM, n::Int; precomputed::Bool=false)

Generate `n` random graphs from the BiCM model `m`. If multithreading is available, the graphs are generated in parallel.

Arguments:

• precomputed::Bool: if true, the precomputed expected biadjacency matrix (m.Ĝ) is used to generate the random graph, otherwise the maximum likelihood parameters are used to generate the random graph on the fly. For larger networks, it is recommended to not precompute the expected adjacency matrix to limit memory pressure.

Note: The generated graph will also be bipartite and respect the layer membership of the original graph used to define the model.

AIC(m::BiCM)

Compute the Akaike Information Criterion (AIC) for the BiCM model m. The parameters of the models most be computed beforehand. If the number of empirical observations becomes too small with respect to the number of parameters, you will get a warning. In that case, the corrected AIC (AICc) should be used instead.

Examples

julia> model = BiCM(corporateclub());

julia> solve_model!(model);

julia> AIC(model);
[...]

See also AICc, L_BiCM_reduced.

AICc(m::BiCM)

Compute the corrected Akaike Information Criterion (AICc) for the BiCM model m. The parameters of the models most be computed beforehand.

Examples

julia> model = BiCM(corporateclub());

julia> solve_model!(model);

julia> AICc(model)
432.12227535579956

See also AIC, L_BiCM_reduced.

BIC(m::BiCM)

Compute the Bayesian Information Criterion (BIC) for the BiCM model m. The parameters of the models most be computed beforehand. BIC is believed to be more restrictive than AIC, as the former favors models with a lower number of parameters than those favored by the latter.

Examples

julia> model = BiCM(corporateclub());

julia> solve_model!(model);

julia> BIC(model)
579.3789571131789

See also AIC, L_BiCM_reduced.

Return the reduced number of nodes in the UBCM network

L_BiCM_reduced(θ::Vector, k⊥::Vector, k⊤::Vector, f⊥::Vector, f⊤::Vector, nz⊥::UnitRange, nz⊤::UnitRange, n⊥ᵣ::Int)

Compute the log-likelihood of the reduced BiCM model using the exponential formulation in order to maintain convexity.

Arguments

• θ: the maximum likelihood parameters of the model ([α; β])
• k⊥: the reduced degree sequence of the ⊥ layer
• k⊤: the reduced degree sequence of the ⊤ layer
• f⊥: the frequency of each degree in the ⊥ layer
• f⊤: the frequency of each degree in the ⊤ layer
• nz⊥: the indices of non-zero elements in the reduced ⊥ layer degree sequence
• nz⊤: the indices of non-zero elements in the reduced ⊤ layer degree sequence
• n⊥ᵣ: the number unique values in the reduced ⊥ layer degree sequence

The function returns the log-likelihood of the reduced model. For the optimisation, this function will be used to generate an anonymous function associated with a specific model.

Examples

# Generic use:
julia> k⊥ = [1, 2, 3, 4];

julia> k⊤  = [1, 2, 4];

julia> f⊥  = [1; 3; 1; 1];

julia> f⊤  = [4; 2; 1];

julia> nz⊥ = 1:length(k⊥);

julia> nz⊤ = 1:length(k⊤);

julia> n⊥ᵣ = length(k⊥);

julia> θ   = collect(range(0.1, step=0.1, length=length(k⊥) + length(k⊤)));

julia> L_BiCM_reduced(θ, k⊥, k⊤, f⊥, f⊤, nz⊥, nz⊤, n⊥ᵣ)
-26.7741690720244

# Use with DBCM model:
julia> G = corporateclub();

julia> model = BiCM(G);

julia> model_fun = θ -> L_BiCM_reduced(θ, model.d⊥ᵣ, model.d⊤ᵣ, model.f⊥, model.f⊤, model.d⊥ᵣ_nz, model.d⊤ᵣ_nz, model.status[:d⊥_unique]);

julia> model_fun(ones(size(model.θᵣ)))
-237.5980041411147

L_BiCM_reduced(m::BiCM)

Return the log-likelihood of the BiCM model m based on the computed maximum likelihood parameters.

Examples

# Use with DBCM model:
julia> G = corporateclub();

julia> model = BiCM(G);

julia> solve_model!(model);

julia> L_BiCM_reduced(model);

See also L_BiCM_reduced(::Vector, ::Vector, ::Vector, ::Vector, ::Vector, ::UnitRange, ::UnitRange, ::Int)

∇L_BiCM_reduced!(∇L::AbstractVector, θ::AbstractVector, k⊥::Vector, k⊤::Vector, f⊥::Vector, f⊤::Vector,  nz⊥::UnitRange{T}, nz⊤::UnitRange{T}, x::AbstractVector, y::AbstractVector, n⊥::Int) where {T<:Signed}

Compute the gradient of the log-likelihood of the reduced DBCM model using the exponential formulation in order to maintain convexity.

For the optimisation, this function will be used togenerate an anonymous function associated with a specific model. The function will update pre-allocated vectors (∇L,x and y) for speed. The gradient is non-allocating.

Arguments

• ∇L: the gradient of the log-likelihood of the reduced model
• θ: the maximum likelihood parameters of the model ([α; β])
• k⊥: the reduced degree sequence of the ⊥ layer
• k⊤: the reduced degree sequence of the ⊤ layer
• f⊥: the frequency of each degree in the ⊥ layer
• f⊤: the frequency of each degree in the ⊤ layer
• nz⊥: the indices of non-zero elements in the reduced ⊥ layer degree sequence
• nz⊤: the indices of non-zero elements in the reduced ⊤ layer degree sequence
• x: the exponentiated maximum likelihood parameters of the model ( xᵢ = exp(-αᵢ) )
• y: the exponentiated maximum likelihood parameters of the model ( yᵢ = exp(-βᵢ) )
• n⊥: the number unique values in the reduced ⊥ layer degree sequence

Examples

# Explicit use with BiCM model:
julia> G = corporateclub();

julia> model = BiCM(G);

julia> ∇L = zeros(Real, length(model.θᵣ));

julia> x  = zeros(Real, length(model.xᵣ));

julia> y  = zeros(Real, length(model.yᵣ));

julia> ∇model_fun! = θ -> ∇L_BiCM_reduced!(∇L, θ, model.d⊥ᵣ, model.d⊤ᵣ, model.f⊥, model.f⊤, model.d⊥ᵣ_nz, model.d⊤ᵣ_nz, x, y, model.status[:d⊥_unique]);

julia> ∇model_fun!(model.θᵣ);

∇L_BiCM_reduced_minus!(args...)

Compute minus the gradient of the log-likelihood of the reduced BiCM model using the exponential formulation in order to maintain convexity. Used for optimisation in a non-allocating manner.

See also ∇L_BiCM_reduced!

BiCM_reduced_iter!(θ::AbstractVector, k⊥::Vector, k⊤::Vector, f⊥::Vector, f⊤::Vector, nz⊥::UnitRange{T}, nz⊤::UnitRange{T}, x::AbstractVector, y::AbstractVector, G::AbstractVector, n⊥::Int) where {T<:Signed}

Compute the next fixed-point iteration for the BiCM model using the exponential formulation in order to maintain convexity. The function is non-allocating and will update pre-allocated vectors (θ, x, y and G) for speed.

Arguments

• θ: the maximum likelihood parameters of the model ([α; β])
• k⊥: the reduced degree sequence of the ⊥ layer
• k⊤: the reduced degree sequence of the ⊤ layer
• f⊥: the frequency of each degree in the ⊥ layer
• f⊤: the frequency of each degree in the ⊤ layer
• nz⊥: the indices of non-zero elements in the reduced ⊥ layer degree sequence
• nz⊤: the indices of non-zero elements in the reduced ⊤ layer degree sequence
• x: the exponentiated maximum likelihood parameters of the model ( xᵢ = exp(-αᵢ) )
• y: the exponentiated maximum likelihood parameters of the model ( yᵢ = exp(-βᵢ) )
• G: buffer for computations
• n⊥: the number unique values in the reduced ⊥ layer degree sequence

Examples

# Use with BiCM model:
julia> G = corporateclub();

julia> model = BiCM(G);

julia> G = zeros(eltype(model.θᵣ), length(model.θᵣ));

julia> x = zeros(eltype(model.θᵣ), length(model.xᵣ));

julia> y = zeros(eltype(model.θᵣ), length(model.yᵣ));


julia> BiCM_FP! = θ -> BiCM_reduced_iter!(θ, model.d⊥ᵣ, model.d⊤ᵣ, model.f⊥, model.f⊤, model.d⊥ᵣ_nz, model.d⊤ᵣ_nz, x, y, G, model.status[:d⊥_unique]);

julia> BiCM_FP!(model.θᵣ);

set_xᵣ!(m::BiCM)

Set the value of xᵣ to exp(-αᵣ) for the BiCM model m

set_yᵣ!(m::BiCM)

Set the value of yᵣ to exp(-βᵣ) for the BiCM model m

Ĝ(m::BiCM)

Compute the expected biadjacency matrix for the BiCM model m

set_Ĝ!(m::BiCM)

Set the expected biadjacency matrix for the BiCM model m

precision(m::BiCM)

Determine the compute precision of the BiCM model m.

Examples

julia> model = BiCM(corporateclub());

julia> MaxEntropyGraphs.precision(model)
Float64

julia> model = BiCM(corporateclub(), precision=Float32);

julia> MaxEntropyGraphs.precision(model)
Float32

A(m::BiCM,i::Int,j::Int)

Return the expected value of the biadjacency matrix for the BiCM model m at the node pair (i,j).

❗ For perfomance reasons, the function does not check:

• if the node pair is valid.
• if the parameters of the model have been computed.

f_BiCM(x::T)

Helper function for the BiCM model to compute the expected value of the biadjacency matrix. The function computes the expression x / (1 + x). As an argument you need to pass the product of the maximum likelihood parameters xᵣ[i] * yᵣ[j] from a BiCM model.

biadjacency_matrix(G::Graphs.SimpleGraph; skipchecks::Bool=false)

Return the biadjacency matrix of the bipartite graph G.

If the graph is not bipartite, an error is thrown. If the adjacency matrix can be written as: A = [O B; B' O] where O is the null matrix, then the returned biadjacency matrix is B and and B' is the transposed biadjacency matrix.

Arguments

• skipchecks: if true, skip the check for the graph being bipartite.

Examples

julia> A = [0 0 0 0 1 0 0;
0 0 0 0 1 1 0;
0 0 0 0 0 0 1;
0 0 0 0 0 1 1;
1 1 0 0 0 0 0;
0 1 0 1 0 0 0;
0 0 1 1 0 0 0];

julia> G = MaxEntropyGraphs.Graphs.SimpleGraph(A);

julia> Array(biadjacency_matrix(A))
4x3 Matrix{Int64}:
 1  0  0
 1  1  0
 0  0  1
 0  1  1

project(G::Graphs.SimpleGraph;  membership::Vector=Graphs.bipartite_map(G), 
                                bottom::Vector=findall(membership .== 1), 
                                top::Vector=findall(membership .== 2); 
                                layer::Symbol=:bottom)

Project the bipartite graph G onto one of its layers and return the projected graph.

Arguments

• membership: the bipartite mapping of the graphs. This can be computed using Graphs.bipartite_map(G).
• bottom: the nodes in the bottom layer. This can be computed using findall(membership .== 1).
• top: the nodes in the top layer. This can be computed using findall(membership .== 2).
• layer: the layer can be specified by passing layer=:bottom or layer=:top. Layer membership is determined by the bipartite map of the graph.
• method: the method used to compute the adjacency matrix of the projected graph. This can be :simple or :weighted. Both methods compute the product of the biadjacency matrix with its transposed, but the :weighted method uses the weights of the edges in the projected graph.

Examples

julia> using Graphs

julia> G = SimpleGraph(5); add_edge!(G, 1, 4); add_edge!(G, 2, 4); add_edge!(G, 3, 4); add_edge!(G, 3, 5);

julia> project(G, layer=:bottom)
{3, 3} undirected simple Int64 graph

julia> project(G, layer=:top)
{2, 1} undirected simple Int64 graph

project(B::T; layer::Symbol=:bottom, method::Symbol=:simple) where {T<:AbstractMatrix}

Project the biadjacency matrix B onto one of its layers.

Arguments

• layer: the layer can be specified by passing layer=:bottom or layer=:top. Layer membership is determined by the bipartite map of the graph.
• method: the method used to compute the adjacency matrix of the projected graph. This can be :simple or :weighted. Both methods compute the product of the biadjacency matrix with its transposed, but the :weighted method uses the weights of the edges in the projected graph.

Examples

julia> B = [0 0 0 1 1; 0 0 0 1 1; 0 0 0 1 0];

julia> project(B, layer=:bottom)
3×3 Matrix{Bool}:
 0  1  1
 1  0  1
 1  1  0

julia> project(B, layer=:top)
5×5 Matrix{Bool}:
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  1
 0  0  0  1  0

julia> B = [0 0 0 1 1; 0 0 0 1 1; 0 0 0 1 0];

julia> project(B, layer=:bottom, method=:weighted)
3×3 Matrix{Int64}:
 0  2  1
 2  0  1
 1  1  0

julia> project(B, layer=:top, method=:weighted)
5×5 Matrix{Int64}:
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  2
 0  0  0  2  0

project(m::BiCM;   α::Float64=0.05, layer::Symbol=:bottom, precomputed::Bool=true, 
                    distribution::Symbol=:Poisson, adjustment::PValueAdjustment=BenjaminiHochberg(),
                        multithreaded::Bool=false)

Obtain the statistically validated projected graph of the BiCM model m onto the layer layer using the V-motifs and the significance level α combined with the p-value adjustment method adjustment.

Arguments

• α: the significance level.
• layer: the layer can be specified by passing layer=:bottom or layer=:top.
• precomputed: if true, the expected values of the biadjacency matrix are used, otherwise the parameters are computed from the model parameters.
• distribution: the distribution used to compute the p-values. This can be :Poisson or :PoissonBinomial.
• adjustment: the method used to adjust the p-values for multiple testing. This can be any of the methods in the PValueAdjustment type (see MultipleTesting.jl). By default, the Benjamini-Hochberg method is used.
• multithreaded: if true, the p-values are computed using multithreading.

Examples

julia> model = BiCM(corporateclub());

julia> solve_model!(model);

julia> project(model, layer=:bottom)
{25, 0} undirected simple Int64 graph

julia> project(model, layer=:top)
{15, 0} undirected simple Int64 graph