DBCM

Maximum entropy model for the Directed Binary Configuration Model (UBCM).

The object holds the maximum likelihood parameters of the model (θ) and optionally the expected adjacency matrix (G), and the variance for the elements of the adjacency matrix (σ). All settings and other metadata are stored in the status field.

DBCM(G::T; precision::N=Float64, kwargs...) where {T<:Graphs.AbstractGraph, N<:Real}
DBCM(;d_out::Vector{T}, d_in::Vector{T}, precision::Type{<:AbstractFloat}=Float64, kwargs...)

Constructor function for the DBCM type.

By default and dependng on the graph type T, the definition of in- and outdegree 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_out, d_in). 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_out, d_in). If you want to work from an adjacency matrix, or edge list, you can use the graph constructors from the $JuliaGraphs$ ecosystem.

Examples

# generating a model from a graph
julia> G = MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques())
{16, 111} directed simple Int64 graph
julia> model = DBCM(G)
DBCM{Graphs.SimpleGraphs.SimpleDiGraph{Int64}, Float64} (16 vertices, 15 unique degree pairs, 0.94 compression ratio)

# generating a model directly from a degree sequence
julia> model = DBCM(d_out=MaxEntropyGraphs.Graphs.outdegree(G), d_in=MaxEntropyGraphs.Graphs.indegree(G))
DBCM{Nothing, Float64} (16 vertices, 15 unique degree pairs, 0.94 compression ratio)

# generating a model directly from a degree sequence with a different precision
julia>  model = DBCM(d_out=MaxEntropyGraphs.Graphs.outdegree(G), d_in=MaxEntropyGraphs.Graphs.indegree(G), precision=Float32)
DBCM{Nothing, Float32} (16 vertices, 15 unique degree pairs, 0.94 compression ratio)

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

julia> G = MaxEntropyGraphs.Graphs.SimpleDiGraph(A)
{3, 5} directed simple Int64 graph
julia> model = DBCM(G)
DBCM{Graphs.SimpleGraphs.SimpleDiGraph{Int64}, Float64} (3 vertices, 3 unique degree pairs, 1.00 compression ratio)

# generating a model from an edge list
julia> E = [(1,2),(1,3),(2,3)];

julia> edgelist = [MaxEntropyGraphs.Graphs.Edge(x,y) for (x,y) in E];

julia> G = MaxEntropyGraphs.Graphs.SimpleDiGraphFromIterator(edgelist)
{3, 3} directed simple Int64 graph
julia> model = DBCM(G)
DBCM{Graphs.SimpleGraphs.SimpleDiGraph{Int64}, Float64} (3 vertices, 3 unique degree pairs, 1.00 compression ratio)

See also Graphs.outdegree, Graphs.indegree.

solve_model!(m::DBCM)

Compute the likelihood maximising parameters of the DBCM 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), :degreesminor, :random, :uniform, or :chunglu.
• 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 = DBCM(MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques()));

julia> solve_model!(model);

# using analytical gradient and degrees minor initial guess
julia> solve_model!(model, method=:BFGS, analytical_gradient=true, initial=:degrees_minor)
(DBCM{Graphs.SimpleGraphs.SimpleDiGraph{Int64}, Float64} (16 vertices, 15 unique degree pairs, 0.94 compression ratio), retcode: Success
u: [3.118482950362848, 2.2567400402511617, 2.2467332710940333, 0.8596258292464105, 0.4957550197436504, 0.3427782029923598, 0.126564995232929, -0.3127732185244699, -0.3967757456352901, -0.43450987676209596  …  -0.5626916621021604, 1.223396713832784, 0.10977479732876981, -1.0367565290851806, -2.0427364999923148, -0.650376357149203, -1.5165614611776657, 0.7532475835319463, 0.39856890694767605, -0.6704522097652438]
Final objective value:     120.15942408828177
)

initial_guess(m::DBCM, method::Symbol=:degrees)

Compute an initial guess for the maximum likelihood parameters of the DBCM 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_{out}); -\log(d_{in})]$
• :degrees_minor: the initial guess is computed using the degrees of the graph and the number of edges, i.e. $\theta = [-\log(d_{out}/(\sqrt{E} + 1)); -\log(d_{in}//(\sqrt{E} + 1) )]$
• :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_{out}/(2E)); -\log(d_{in}/(2E))]$

Examples

julia> model = DBCM(MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques()));

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

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

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

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

julia> initial_guess(model);

rand(m::DBCM; precomputed=false)

Generate a random graph from the DBCM model m.

Arguments:

• precomputed::Bool: if true, the precomputed expected adjacency 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.

Examples

# generate a DBCM model macaques network
julia> G = MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques());

julia> model = DBCM(G); 

julia> solve_model!(model); # compute the maximum likelihood parameters

julia> typeof(rand(model))
Graphs.SimpleGraphs.SimpleDiGraph{Int64}

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

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

Arguments

• precomputed::Bool: if true, the precomputed expected adjacency 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.

Examples

# generate a DBCM model macaques network
julia> G = MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques());

julia> model = DBCM(G); 

julia> solve_model!(model); # compute the maximum likelihood parameters

julia> typeof(rand(model, 10))
Vector{SimpleDiGraph{Int64}} (alias for Array{Graphs.SimpleGraphs.SimpleDiGraph{Int64}, 1})

AIC(m::DBCM)

Compute the Akaike Information Criterion (AIC) for the DBCM 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 = DBCM(MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques()));

julia> solve_model!(model);

julia> AIC(model);
┌ Warning: The number of observations is small with respect to the number of parameters (n/k < 40). Consider using the corrected AIC (AICc) instead.
[...]

See also AICc, L_DBCM_reduced.

AICc(m::DBCM)

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

Examples

julia> model = DBCM(MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques()));

julia> solve_model!(model);

julia> AICc(model)
314.5217467272881

See also AIC, L_DBCM_reduced.

BIC(m::DBCM)

Compute the Bayesian Information Criterion (BIC) for the DBCM 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 = DBCM(MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques()));

julia> solve_model!(model);

julia> BIC(model)
415.69929372350714

See also AIC, L_DBCM_reduced.

Return the reduced number of nodes in the UBCM network

L_DBCM_reduced(θ::Vector, k_out::Vector, k_in::Vector, F::Vector, nz_out::Vector, nz_in::Vector, n::Int=length(k_out))

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

Arguments

- `θ`: the maximum likelihood parameters of the model ([α; β])
- `k_out`: the reduced outdegree sequence
- `k_in`: the reduced indegree sequence
- `F`: the frequency of each pair in the degree sequence
- `nz_out`: the indices of non-zero elements in the reduced outdegree sequence
- `nz_in`: the indices of non-zero elements in the reduced indegree sequence
- `n`: the number of nodes in the reduced model

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_out  = [1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 5];

julia> k_in   = [2, 3, 4, 1, 3, 5, 2, 4, 1, 2, 4, 0, 4];

julia> F      = [2, 2, 1, 1, 1, 2, 3, 1, 1, 2, 2, 1, 1];

julia> θ      = collect(range(0.1, step=0.1, length=length(k_out)));

julia> nz_out = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13];

julia> nz_in  = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13];

julia> n      = length(k_out);

julia> L_DBCM_reduced(θ, k_out, k_in, F, nz_out, nz_in, n);

# Use with DBCM model:
julia> G = MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques());

julia> model = DBCM(G);

julia> model_fun = θ -> L_DBCM_reduced(θ, model.dᵣ_out, model.dᵣ_in, model.f, model.dᵣ_out_nz, model.dᵣ_in_nz, model.status[:d_unique]);

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

L_DBCM_reduced(m::DBCM)

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

Examples

# Use with DBCM model:
julia> G = MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques());

julia> model = DBCM(G);

julia> solve_model!(model);

julia> L_DBCM_reduced(model);

See also L_DBCM_reduced(::Vector, ::Vector, ::Vector, ::Vector, ::Vector, ::Vector)

∇L_DBCM_reduced!(∇L::AbstractVector, θ::AbstractVector, k_out::AbstractVector, k_in::AbstractVector, F::AbstractVector, nz_out::Vector, nz_in::Vector, x::AbstractVector, y::AbstractVector,n::Int)

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_out: the reduced outdegree sequence
• k_in: the reduced indegree sequence
• F: the frequency of each pair in the degree sequence
• nz_out: the indices of non-zero elements in the reduced outdegree sequence
• nz_in: the indices of non-zero elements in the reduced indegree 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 of nodes in the reduced model

Examples

# Explicit use with DBCM model:
julia> G = MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques());

julia> model = DBCM(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_DBCM_reduced!(∇L, θ, model.dᵣ_out, model.dᵣ_in, model.f, model.dᵣ_out_nz, model.dᵣ_in_nz, x, y, model.status[:d_unique]);

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

# Use within optimisation.jl framework:
julia> fun = (θ,p) -> -L_DBCM_reduced(θ, model.dᵣ_out, model.dᵣ_in, model.f, model.dᵣ_out_nz, model.dᵣ_in_nz, model.status[:d_unique]);

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

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

julia> ∇fun! = (∇L, θ, p) -> ∇L_DBCM_reduced!(∇L, θ, model.dᵣ_out, model.dᵣ_in, model.f, model.dᵣ_out_nz, model.dᵣ_in_nz, x, y, model.status[:d_unique]);

julia> θ₀ = initial_guess(model); # initial condition

julia> foo = MaxEntropyGraphs.Optimization.OptimizationFunction(fun, grad=∇fun!); # define target function 

julia> prob  = MaxEntropyGraphs.Optimization.OptimizationProblem(foo, θ₀); # define the optimisation problem

julia> method = MaxEntropyGraphs.OptimizationOptimJL.LBFGS(); # set the optimisation method

julia> MaxEntropyGraphs.Optimization.solve(prob, method); # solve it

∇L_DBCM_reduced_minus!(args...)

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

See also ∇L_DBCM_reduced!

DBCM_reduced_iter!(θ::AbstractVector, k_out::AbstractVector, k_in::AbstractVector, F::AbstractVector, nz_out::Vector, nz_in::Vector,x::AbstractVector, y::AbstractVector, G::AbstractVector, H::AbstractVector, n::Int)

Compute the next fixed-point iteration for the DBCM 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_out: the reduced outdegree sequence
• k_in: the reduced indegree sequence
• F: the frequency of each pair in the degree sequence
• nz_out: the indices of non-zero elements in the reduced outdegree sequence
• nz_in: the indices of non-zero elements in the reduced indegree 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 of nodes in the reduced model

Examples

# Use with DBCM model:
julia> G = MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques());

julia> model = DBCM(G);

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

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

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

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

julia> DBCM_FP! = θ -> DBCM_reduced_iter!(θ, model.dᵣ_out, model.dᵣ_in, model.f, model.dᵣ_out_nz, model.dᵣ_in_nz, x, y, G, model.status[:d_unique]);

julia> DBCM_FP!(model.θᵣ);

set_xᵣ!(m::DBCM)

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

set_yᵣ!(m::DBCM)

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

Ĝ(m::DBCM)

Compute the expected adjacency matrix for the DBCM model m

set_Ĝ!(m::DBCM)

Set the expected adjacency matrix for the DBCM model m

σˣ(m::DBCM{T,N}) where {T,N}

Compute the standard deviation for the elements of the adjacency matrix for the DBCM model m.

Note: read as "sigma star"

set_σ!(m::DBCM)

Set the standard deviation for the elements of the adjacency matrix for the DBCM model m

precision(m::DBCM)

Determine the compute precision of the DBCM model m.

Examples

julia> model = DBCM(MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques()));

julia> MaxEntropyGraphs.precision(model)
Float64

julia> model = DBCM(MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques()), precision=Float32);

julia> MaxEntropyGraphs.precision(model)
Float32

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

Return the expected value of the adjacency matrix for the DBCM 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_DBCM(x::T)

Helper function for the DBCM model to compute the expected value of the adjacency 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 DBCM model.

σₓ(m::DBCM, X::function)

Compute the standard deviation of metric X for the DBCM model m.

This requires that both the expected values (m.Ĝ) and standard deviations (m.σ) are computed for m.