RBCM

Maximum entropy model for the Reciprocal Binary Configuration Model (RBCM).

The model constrains, for every node, the non-reciprocated out-degree $k^{→}_i = \sum_{j≠i} a_{ij}(1-a_{ji})$, the non-reciprocated in-degree $k^{←}_i = \sum_{j≠i} a_{ji}(1-a_{ij})$ and the reciprocated degree $k^{↔}_i = \sum_{j≠i} a_{ij}a_{ji}$ (Squartini & Garlaschelli (2011), App. C.1, there called the reciprocal configuration model RCM; also used in Di Vece et al. (2023) and the NuMeTriS package).

The object holds the maximum likelihood parameters of the model (θᵣ = [α; β; γ]), and optionally the expected adjacency matrix (Ĝ) and the variance of its elements (σ). All settings and other metadata are stored in the status field.

Note: within a dyad the entries $a_{ij}$ and $a_{ji}$ are not independent under the RBCM (unlike the DBCM): their covariance is $p^{↔}_{ij} - ⟨a_{ij}⟩⟨a_{ji}⟩ ≠ 0$. Metrics that combine both directions of a dyad (e.g. motifs) must therefore be computed from the dyadic probabilities, not from Ĝ alone; the package's σₓ and motif methods for the RBCM take care of this.

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

Constructor function for the RBCM type.

By default the three reciprocal degree sequences are computed from the directed graph G (see nonreciprocated_outdegree, nonreciprocated_indegree and reciprocated_degree). You can also pass the degree sequences directly as keyword arguments (d_out = k→, d_in = k←, d_rec = k↔), either alongside the graph or without one.

Examples

# generating a model from a graph
julia> G = MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques());

julia> model = RBCM(G)
RBCM{Graphs.SimpleGraphs.SimpleDiGraph{Int64}, Float64} (16 vertices, 15 unique degree triples, 0.94 compression ratio)

# generating a model directly from degree sequences
julia> model = RBCM(d_out=nonreciprocated_outdegree(G), d_in=nonreciprocated_indegree(G), d_rec=reciprocated_degree(G))
RBCM{Nothing, Float64} (16 vertices, 15 unique degree triples, 0.94 compression ratio)

solve_model!(m::RBCM)

Compute the likelihood maximising parameters of the RBCM model m.

Arguments

• method::Symbol: solution method to use, can be :fixedpoint (default), or :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::Union{Real, Nothing}: tolerance for the fixedpoint method (defaults to nothing, i.e. 1e-8). It bounds the fixed-point increment $\|G(\theta) - \theta\|_\infty$ in parameter space; it is not the constraint residual, and it is ignored by every other method. Use constraint_residual to measure how well the expected degrees actually match the observed ones.
• 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).
• g_tol::Union{Number, Nothing}: gradient tolerance for the gradient-based methods (maps to Optim's g_abstol); set e.g. 1e-5 to stop before over-converging (defaults to nothing, i.e. Optim's tight default). The gradient of this model is its constraint residual (up to the degree multiplicities), but g_abstol is a stopping criterion rather than a guarantee: Optim can also stop on its function or parameter convergence checks and report success without the gradient ever reaching g_tol. Verify what was actually achieved with constraint_residual.
• 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 = RBCM(MaxEntropyGraphs.Graphs.SimpleDiGraph(rhesus_macaques()));

julia> solve_model!(model);

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

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

The methods available are:

• :degrees (default): the initial guess is computed using the degree sequences, i.e. $\theta = [-\log(k^{→}); -\log(k^{←}); -\log(k^{↔})]$
• :degrees_minor: the initial guess is computed using the degree sequences and the number of edges, i.e. $\theta_{ch} = -\log(k_{ch}/(\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 degree sequences and the number of edges, i.e. $\theta_{ch} = -\log(k_{ch}/(2E))$

Entries whose constraint is zero get an infinite parameter (their fitness is exactly zero); the solution methods handle these separately.

Examples

julia> model = RBCM(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::RBCM; precomputed::Bool=false)

Generate a random graph from the RBCM model m. Each dyad (i,j) is sampled independently from its four-state distribution (single link i→j, single link j→i, reciprocated pair, absent) using a single uniform draw against the cumulative probabilities (p⭢, p⭠, p⭤).

Arguments:

• precomputed::Bool: only false is supported for the RBCM. The expected adjacency matrix Ĝ does not characterise the dyadic joint distribution (the states of aᵢⱼ and aⱼᵢ are correlated), so sampling from Ĝ would not reproduce the reciprocated degree sequence. An ArgumentError is thrown when precomputed=true.

Examples

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

julia> model = RBCM(G);

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

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

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

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

Arguments:

• precomputed::Bool: only false is supported for the RBCM (see rand(::RBCM)).

Examples

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

julia> model = RBCM(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::RBCM)

Compute the Akaike Information Criterion (AIC) for the RBCM 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.

The number of parameters is 3N and the number of observations N(N-1) (the ordered node pairs), identical to the DBCM convention, so AIC values of a DBCM and an RBCM fitted to the same network are directly comparable.

Examples

julia> model = RBCM(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_RBCM_reduced.

AICc(m::RBCM)

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

Examples

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

julia> solve_model!(model);

julia> AICc(model);

See also AIC, L_RBCM_reduced.

BIC(m::RBCM)

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

julia> solve_model!(model);

julia> BIC(model);

See also AIC, L_RBCM_reduced.

Return the reduced number of nodes in the RBCM network

L_RBCM_reduced(θ::AbstractVector, k_out::Vector, k_in::Vector, k_rec::Vector, F::Vector, nz_out::Vector, nz_in::Vector, nz_rec::Vector, n::Int=length(k_out))

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

The RBCM likelihood is (Squartini & Garlaschelli (2011), Eq. C.14, written in terms of θ = [α; β; γ])

$\mathcal{L} = -\sum_i \left[ k^{→}_i α_i + k^{←}_i β_i + k^{↔}_i γ_i \right] - \sum_{i<j} \ln\left( 1 + x_iy_j + x_jy_i + z_iz_j \right)$

with $x = e^{-α}, y = e^{-β}, z = e^{-γ}$. The dyadic normaliser is evaluated with a numerically stable four-term log-sum-exp (MaxEntropyGraphs.log1pexpsum). Channels whose constraint is zero (e.g. $k^{↔}_i = 0$) are treated as being at their analytical optimum (parameter +Inf, fitness exactly zero), so the likelihood does not depend on those coordinates.

Arguments

• θ: the maximum likelihood parameters of the model ([α; β; γ])
• k_out: the reduced non-reciprocated out-degree sequence (k→)
• k_in: the reduced non-reciprocated in-degree sequence (k←)
• k_rec: the reduced reciprocated degree sequence (k↔)
• F: the frequency of each (k→, k←, k↔) triple in the degree sequences
• nz_out: the indices of non-zero elements in k_out
• nz_in: the indices of non-zero elements in k_in
• nz_rec: the indices of non-zero elements in k_rec
• n: the number of unique node classes in the reduced model

Examples

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

julia> model = RBCM(G);

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

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

L_RBCM_reduced(m::RBCM)

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

Examples

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

julia> model = RBCM(G);

julia> solve_model!(model);

julia> L_RBCM_reduced(model);

See also L_RBCM_reduced(::Vector, ::Vector, ::Vector, ::Vector, ::Vector, ::Vector, ::Vector, ::Vector)

∇L_RBCM_reduced!(∇L::AbstractVector, θ::AbstractVector, k_out::AbstractVector, k_in::AbstractVector, k_rec::AbstractVector, F::AbstractVector, nz_out::Vector, nz_in::Vector, nz_rec::Vector, x::AbstractVector, y::AbstractVector, z::AbstractVector, n::Int)

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

For the optimisation, this function will be used to generate an anonymous function associated with a specific model. The function will update pre-allocated vectors (∇L, x, y and z) 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 non-reciprocated out-degree sequence (k→)
• k_in: the reduced non-reciprocated in-degree sequence (k←)
• k_rec: the reduced reciprocated degree sequence (k↔)
• F: the frequency of each triple in the degree sequences
• nz_out, nz_in, nz_rec: the indices of non-zero elements in the respective reduced degree sequences
• x, y, z: buffers for the exponentiated maximum likelihood parameters
• n: the number of unique node classes in the reduced model

Examples

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

julia> model = RBCM(G);

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

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

julia> ∇model_fun! = θ -> ∇L_RBCM_reduced!(∇L, θ, model.dᵣ_out, model.dᵣ_in, model.dᵣ_rec, model.f, model.dᵣ_out_nz, model.dᵣ_in_nz, model.dᵣ_rec_nz, x, y, z, model.status[:d_unique]);

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

∇L_RBCM_reduced_minus!(args...)

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

See also ∇L_RBCM_reduced!

RBCM_reduced_iter!(θ::AbstractVector, k_out::AbstractVector, k_in::AbstractVector, k_rec::AbstractVector, F::AbstractVector, nz_out::Vector, nz_in::Vector, nz_rec::Vector, x::AbstractVector, y::AbstractVector, z::AbstractVector, G::AbstractVector, n::Int)

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

Arguments

• θ: the maximum likelihood parameters of the model ([α; β; γ])
• k_out, k_in, k_rec: the reduced degree sequences (k→, k←, k↔)
• F: the frequency of each triple in the degree sequences
• nz_out, nz_in, nz_rec: the indices of non-zero elements in the respective reduced degree sequences
• x, y, z: buffers for the exponentiated maximum likelihood parameters
• G: buffer for computations
• n: the number of unique node classes in the reduced model

Examples

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

julia> model = RBCM(G);

julia> buf = zeros(length(model.θᵣ));

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

julia> RBCM_FP! = θ -> RBCM_reduced_iter!(θ, model.dᵣ_out, model.dᵣ_in, model.dᵣ_rec, model.f, model.dᵣ_out_nz, model.dᵣ_in_nz, model.dᵣ_rec_nz, x, y, z, buf, model.status[:d_unique]);

julia> RBCM_FP!(ones(size(model.θᵣ)));

log1pexpsum(a, b, c)

Numerically stable evaluation of log(1 + exp(a) + exp(b) + exp(c)) (a four-term log-sum-exp with an implicit unit term). This is the analog of softplus for the RBCM's dyadic normaliser: $\ln D_{ij} = \ln(1 + x_iy_j + x_jy_i + z_iz_j) =$ log1pexpsum(-(αᵢ+βⱼ), -(αⱼ+βᵢ), -(γᵢ+γⱼ)).

Computed by factoring out m = max(0, a, b, c), which avoids overflow for large positive arguments and degrades gracefully for -Inf arguments (channels pinned at their analytical optimum contribute an exact zero: exp(-Inf) = 0, and m ≥ 0 remains finite).

set_xᵣ!(m::RBCM)

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

set_yᵣ!(m::RBCM)

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

set_zᵣ!(m::RBCM)

Set the value of zᵣ to exp(-γᵣ) for the RBCM model m

Ĝ(m::RBCM)

Compute the expected adjacency matrix for the RBCM model m ($⟨a_{ij}⟩ = p^{→}_{ij} + p^{↔}_{ij}$).

Note: under the RBCM the entries aᵢⱼ and aⱼᵢ are correlated within a dyad, so Ĝ alone does not fully characterise the dyadic joint distribution (see MaxEntropyGraphs._dyadic_probability_matrices).

set_Ĝ!(m::RBCM)

Set the expected adjacency matrix for the RBCM model m

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

Compute the standard deviation for the elements of the adjacency matrix for the RBCM model m ($σ[a_{ij}] = \sqrt{⟨a_{ij}⟩(1-⟨a_{ij}⟩)}$).

Note: read as "sigma star". The aᵢⱼ remain Bernoulli marginally; the within-dyad correlation between aᵢⱼ and aⱼᵢ is captured separately (see MaxEntropyGraphs._cov_dyads) and is accounted for by σₓ.

set_σ!(m::RBCM)

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

_dyadic_probability_matrices(m::RBCM)

Compute the three dyadic probability matrices (P̂, R̂, Ẑ) of the RBCM model m, where P̂[i,j] = p⭢(m,i,j) (single link i→j), R̂[i,j] = p⭤(m,i,j) (reciprocated dyad) and Ẑ[i,j] = p∅(m,i,j) (absent dyad). The matrix of p⭠ probabilities is transpose(P̂). To match the conventions of the motif counting kernels (_motif_base_matrices/_motif_count), Ẑ carries a unit diagonal while P̂ and R̂ have zero diagonals.

These matrices are the correct expected counterparts of the motif base matrices P/Q/R/Z under the RBCM: because aᵢⱼ and aⱼᵢ are correlated within a dyad, they can not be obtained from the expected adjacency matrix Ĝ (e.g. ⟨aᵢⱼaⱼᵢ⟩ = p⭤ᵢⱼ ≠ ĜᵢⱼĜⱼᵢ).

_cov_dyads(m::RBCM)

Compute the within-dyad covariance matrix C of the RBCM model m, with $C_{ij} = Cov(a_{ij}, a_{ji}) = p^{↔}_{ij} - ⟨a_{ij}⟩⟨a_{ji}⟩$ for i ≠ j and zero diagonal (Squartini & Garlaschelli (2011), Eq. C.17). The matrix is symmetric; distinct dyads are independent.

precision(m::RBCM)

Determine the compute precision of the RBCM model m.

Examples

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

julia> MaxEntropyGraphs.precision(model)
Float64

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

Return the expected value of the adjacency matrix for the RBCM model m at the node pair (i,j) ($⟨a_{ij}⟩ = p^{→}_{ij} + p^{↔}_{ij}$).

❗ For performance reasons, the function does not check:

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

p⭢(m::RBCM, i::Int, j::Int)

Return the probability that the dyad (i,j) holds a single non-reciprocated link i→j under the RBCM model m ($p^{→}_{ij} = x_iy_j/D_{ij}$ with $D_{ij} = 1 + x_iy_j + x_jy_i + z_iz_j$).

❗ For performance reasons, the function does not check if the node pair is valid or if the parameters have been computed.

p⭠(m::RBCM, i::Int, j::Int)

Return the probability that the dyad (i,j) holds a single non-reciprocated link j→i under the RBCM model m ($p^{←}_{ij} = x_jy_i/D_{ij} = p^{→}_{ji}$).

❗ For performance reasons, the function does not check if the node pair is valid or if the parameters have been computed.

p⭤(m::RBCM, i::Int, j::Int)

Return the probability that the dyad (i,j) holds a reciprocated link pair (i→j and j→i) under the RBCM model m ($p^{↔}_{ij} = z_iz_j/D_{ij}$).

❗ For performance reasons, the function does not check if the node pair is valid or if the parameters have been computed.

p∅(m::RBCM, i::Int, j::Int)

Return the probability that the dyad (i,j) holds no links under the RBCM model m ($p^{∅}_{ij} = 1/D_{ij}$).

❗ For performance reasons, the function does not check if the node pair is valid or if the parameters have been computed.

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

Compute the standard deviation of metric X for the RBCM model m using the delta method of Squartini & Garlaschelli (2011) (Eqs. B.16/C.17), including the within-dyad covariance term:

$(σ^{*}[X])^2 = \sum_{i,j} \left[ σ^2[a_{ij}] \left(\frac{∂X}{∂a_{ij}}\right)^2 + Cov(a_{ij},a_{ji}) \frac{∂X}{∂a_{ij}}\frac{∂X}{∂a_{ji}} \right]_{A = ⟨A⟩}$

Under the RBCM the covariance $Cov(a_{ij},a_{ji}) = p^{↔}_{ij} - ⟨a_{ij}⟩⟨a_{ji}⟩$ is non-zero, unlike for the DBCM. This requires that both the expected values (m.Ĝ) and standard deviations (m.σ) are computed for m.