DBiCM

Maximum entropy model for the Directed Bipartite Configuration Model (DBiCM).

A bipartite network whose edges carry a direction has two disjoint link channels between the ⊥ and the ⊤ layer, each an N⊥ × N⊤ binary matrix:

• channel ⁺ (:to_top, ⊥ → ⊤): B⁺[i,α] = 1 iff ⊥ᵢ → ⊤α
• channel ⁻ (:to_bottom, ⊤ → ⊥): B⁻[i,α] = 1 iff ⊤α → ⊥ᵢ

Constraining all four degree sequences — the ⊥ layer's out- and in-degrees (d⊥_out, d⊥_in) and the ⊤ layer's (d⊤_out, d⊤_in) — gives the Hamiltonian

$H = Σᵢ (αᵢ k^{out}_i + βᵢ k^{in}_i) + Σ_α (γ_α h^{out}_α + δ_α h^{in}_α)$

which, collected per entry, reads Σ_{i,α} (αᵢ + δ_α) B⁺_{iα} + Σ_{i,α} (βᵢ + γ_α) B⁻_{iα}. No term couples the two channels, so the partition function factorises over all 2·N⊥·N⊤ entries and

• channel ⁺ is exactly a BiCM on the sequence pair (d⊥_out, d⊤_in), and
• channel ⁻ is exactly a BiCM on the sequence pair (d⊥_in, d⊤_out),

statistically independent of one another. The connection probabilities are

$p⁺_{iα} = \frac{x⊥^{out}_i x⊤^{in}_α}{1 + x⊥^{out}_i x⊤^{in}_α}$, $\qquad p⁻_{iα} = \frac{x⊥^{in}_i x⊤^{out}_α}{1 + x⊥^{in}_i x⊤^{out}_α}$

with x = exp(-θ). In particular B⁺_{iα} and B⁻_{iα} are independent, so the model carries no reciprocity coupling: it is the directed bipartite analogue of the DBCM, not of the RBCM. A model that constrains reciprocated bipartite degrees would be a different object.

Class reduction

Each of the four sequences reduces to unique-value classes independently. This is forced, not merely convenient: the DBCM needs joint (out, in) classes only because its sums carry the self-exclusion i ≠ j, so a node's out-equation depends on its own class. Bipartite sums range over two disjoint vertex sets and have no such term, and the cross-channel second derivatives of the log-likelihood vanish identically, so the joint pair (k^{out}_i, k^{in}_i) never enters at all. Four independent reductions are therefore never coarser than a DBCM-style joint reduction, and usually strictly finer.

Field naming

<quantity>⟨layer⟩[ᵣ]_⟨direction⟩, composing the BiCM's ⊥/⊤ with the DBCM's out/in. All four exponentiated parameter vectors share the x prefix — with four fitnesses there is no natural x/y dichotomy, and x⊥ᵣ_out * x⊤ᵣ_in reads unambiguously.

The parameter vector is laid out as two channel blocks, each in BiCM order (⊥ side, then ⊤ side): θᵣ = [α(⊥out) ; δ(⊤in) | β(⊥in) ; γ(⊤out)].

DBiCM(G::T; d⊥_out, d⊥_in, d⊤_out, d⊤_in, precision::Type{<:AbstractFloat}=Float64, kwargs...)
DBiCM(; d⊥_out::Vector, d⊥_in::Vector, d⊤_out::Vector, d⊤_in::Vector, precision=Float64, kwargs...)

Construct a DBiCM model from a directed bipartite graph G, or from the four degree sequences directly (pass G = nothing, or use the keyword-only form).

The layer membership is read off the undirected skeleton of G. This matters: Graphs.bipartite_map traverses outneighbors only, so on a directed graph it explores just the out-reachable set of its seed and leaves every other vertex at the default colour — while Graphs.is_bipartite still returns true, because an empty traversal encounters no conflict.

Arguments

• d⊥_out, d⊥_in: out- and in-degree sequences of the ⊥ layer (defaults read from G).
• d⊤_out, d⊤_in: out- and in-degree sequences of the ⊤ layer (defaults read from G).
• precision: compute precision, defaults to Float64.

Rejected input

• a graph containing a vertex with no edges in either direction — its layer is not determined by the data. A vertex that only sends, or only receives, is fine.
• degree sequences that no directed bipartite graph realises: Σd⊥_out == Σd⊤_in and Σd⊤_out == Σd⊥_in are two separate link counts, one per channel, and both must hold.
• a vertex whose degree saturates the live opposite side of its channel, whose fitness diverges.

Examples

julia> using Graphs

julia> G = SimpleDiGraph(8);

julia> for (a,b) in ((1,5),(1,6),(2,6),(2,7),(3,7),(4,8),(5,2),(6,3),(7,4),(8,1)); add_edge!(G, a, b); end

julia> model = DBiCM(G)
DBiCM{Graphs.SimpleGraphs.SimpleDiGraph{Int64}, Float64} (4 + 4 vertices, 6 parameters, 0.38 compression ratio)

solve_model!(m::DBiCM; kwargs...)

Compute the maximum-likelihood parameters of the DBiCM model m.

The two channels are solved separately. They share no parameter and their cross derivatives vanish identically, so this is the exact decomposition rather than an approximation — and it keeps each solve at one gauge direction, the regime the fixed-point accelerator is measured on, at a Newton cost of n⁺³ + n⁻³ instead of (n⁺+n⁻)³.

A channel with no links at all (a purely one-directional network, which is a perfectly ordinary input) is not solved: its optimum is θ ≡ Inf, p ≡ 0 exactly, and running the fixed-point map on it would evaluate -log(0/0).

Arguments

• method: :fixedpoint (default), :BFGS, :LBFGS or :Newton.
• initial: :degrees (default), :random, :uniform or :chung_lu.
• maxiters, verbose, ftol, abstol, reltol, g_tol, AD_method, analytical_gradient: as for the other models.

Returns (m, sol), where sol is a NamedTuple (out = ..., in = ...) holding each channel's solution object, or nothing for a channel that carried no links.

Examples

julia> using Graphs

julia> G = SimpleDiGraph(8);

julia> for (a,b) in ((1,5),(1,6),(2,6),(2,7),(3,7),(4,8),(5,2),(6,3),(7,4),(8,1)); add_edge!(G, a, b); end

julia> model = DBiCM(G);

julia> solve_model!(model);

julia> model.status[:params_computed]
true

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

Compute an initial guess for the parameters of the DBiCM model m.

Supported methods: :degrees (default), :random, :uniform, :chung_lu. Unlike the BiCM's, :chung_lu does not require the model to hold a graph — each channel's link count is already in m.status.

rand(m::DBiCM; precomputed::Bool=false, rng::AbstractRNG=default_rng())

Sample a directed bipartite graph from the ensemble of the DBiCM model m, as a SimpleDiGraph in the original vertex numbering of m.G. Both channels are sampled independently.

Set precomputed = true to draw from the stored matrices (set_Ĝ! first) rather than recomputing the probabilities on the fly.

rand(m::DBiCM, n::Int; precomputed::Bool=false, rng::AbstractRNG=default_rng())

Sample n directed bipartite graphs from the ensemble of the DBiCM model m.

AIC(m::DBiCM)

Akaike information criterion of the DBiCM model m. The model has 2(N⊥ + N⊤) parameters and 2·N⊥·N⊤ observations — two independent matrices of Bernoulli entries. As for the BiCM, the gauge freedom is not discounted from the parameter count.

AICc(m::DBiCM)

Corrected Akaike information criterion of the DBiCM model m, for a small number of observations.

BIC(m::DBiCM)

Bayesian information criterion of the DBiCM model m.

Return the reduced number of parameters of the DBiCM network (both channels combined)

L_DBiCM_reduced(θ::AbstractVector, m::DBiCM)
L_DBiCM_reduced(m::DBiCM)

Log-likelihood of the reduced DBiCM model, the sum of the two channels' BiCM log-likelihoods.

∇L_DBiCM_reduced!(∇L, θ, m::DBiCM, x⊥_out, x⊤_in, x⊥_in, x⊤_out)

Gradient of the reduced DBiCM log-likelihood, computed in place. The four trailing arguments are scratch buffers for the exponentiated parameters, one per reduced sequence.

∇L_DBiCM_reduced_minus!(∇L, θ, m::DBiCM, x⊥_out, x⊤_in, x⊥_in, x⊤_out)

As ∇L_DBiCM_reduced!, returning the negative gradient (for minimisation).

DBiCM_reduced_iter!(θ, m::DBiCM, x⊥_out, x⊤_in, x⊥_in, x⊤_out, G)

One fixed-point iteration of the reduced DBiCM model, written into the buffer G, which is returned. Both channels are advanced independently.

_channel_id(channel::Symbol)

Normalise a channel selector to :plus (⊥ → ⊤) or :minus (⊤ → ⊥).

channel is absolute — it names the direction links travel — while the kind selector of the projection functions is relative to the projected layer. Keeping the two vocabularies disjoint is deliberate: under layer = :top, a layer-relative :out means "both ⊤ nodes send", which is channel ⁻, so a shared :out/:in spelling would have the same symbol select opposite matrices in the two APIs. :out/:in are therefore rejected here with a pointer rather than silently accepted.

_channel(m::DBiCM, channel::Symbol)

Everything the BiCM kernels need to evaluate one channel of m, as a NamedTuple: the θᵣ index range of that channel's block, its own/opposite reduced sequences, multiplicities, live ranges, class maps and exponentiated parameters, and its link count.

own is always the ⊥ side and opp the ⊤ side, so both channels present themselves to the BiCM kernels in the same order and no kernel ever has to know which channel it is looking at.

_check_vertex(m::DBiCM, i::Int)

Validate that m is solved and that i is a vertex of it. Called before any layer lookup, so an out-of-range vertex raises a clear ArgumentError rather than a BoundsError from inside is⊥.

_dbicm_expected_degree(m::DBiCM, i::Int, channel::Symbol, method::Symbol)

Expected number of links of graph vertex i in one channel. channel is absolute, so for a ⊥ vertex :to_top is its out-degree while for a ⊤ vertex :to_top is its in-degree.

_solve_dbicm_channel(m::DBiCM, ch::Symbol, θ_ch, method, ...)

Solve one channel of m, returning (θ, sol). The channel is an ordinary BiCM problem, so this is the BiCM solve recipe verbatim — including the Newton/Zygote second-order workaround and the "bind the model fields to locals" precaution — applied to that channel's own sequences.

set_xᵣ!(m::DBiCM)

Set the exponentiated maximum-likelihood parameters of all four reduced sequences of m.

Ĝ(m::DBiCM; channel::Symbol=:to_top)

Expected biadjacency matrix of m, N⊥ × N⊤, with rows indexed by the ⊥ layer and columns by the ⊤ layer in both cases.

channel = :to_top (default) gives ⟨B⁺⟩, the ⊥ → ⊤ links; :to_bottom gives ⟨B⁻⟩, the ⊤ → ⊥ links. channel = :both returns the two matrices stacked vertically, [⟨B⁺⟩ ; ⟨B⁻⟩] (2N⊥ × N⊤), which is the form σₓ expects for a metric that reads both channels at once.

set_Ĝ!(m::DBiCM)

Compute and store both expected biadjacency matrices of m (m.Ĝ⁺ and m.Ĝ⁻).

σˣ(m::DBiCM; channel::Symbol=:to_top)

Entry-wise standard deviations of the expected biadjacency matrix of m, σ[B] = sqrt(p(1-p)). channel is as for Ĝ, including :both.

set_σ!(m::DBiCM)

Compute and store both standard-deviation matrices of m (m.σ⁺ and m.σ⁻).

precision(m::DBiCM)

Return the compute precision of the DBiCM model m.

A(m::DBiCM, i::Int, α::Int; channel::Symbol=:to_top)

Expected value of a single biadjacency entry, with layer-local indices. channel selects the direction: :to_top for ⊥ᵢ → ⊤α, :to_bottom for ⊤α → ⊥ᵢ.

p⁺(m::DBiCM, i::Int, α::Int)

Probability of the link ⊥ᵢ → ⊤α under m, with layer-local indices and no bounds checking.

p⁻(m::DBiCM, i::Int, α::Int)

Probability of the link ⊤α → ⊥ᵢ under m, with layer-local indices and no bounds checking.

σₓ(m::DBiCM, X::Function; channel::Symbol=:to_top, gradient_method::Symbol=:ReverseDiff)

Standard deviation of the metric X under the DBiCM model m, by the delta method.

Every biadjacency entry of a DBiCM is an independent Bernoulli variable — within a channel, across channels, and within a (i,α) dyad — so there are no covariance terms anywhere and

$σ²[X] = Σ (σ_{iα} ∂X/∂B_{iα})²$

With channel = :both, X is a function of the stacked 2N⊥ × N⊤ matrix [B⁺ ; B⁻] and the sum runs over both blocks; this is exact for the same reason, and is how metrics that mix the two channels (total links, reciprocity, path motifs) are propagated.

reciprocity(m::DBiCM)

Expected topological reciprocity of the DBiCM model m, as the ratio of expectations

$r = 2 Σ_{iα} p⁺_{iα} p⁻_{iα} / Σ_{iα} (p⁺_{iα} + p⁻_{iα})$

B⁺_{iα} and B⁻_{iα} are independent, so ⟨B⁺B⁻⟩ = p⁺p⁻ exactly and no approximation enters. The factor of two makes this agree with reciprocity(A::AbstractMatrix) applied to the full expected adjacency matrix [0 ⟨B⁺⟩ ; ⟨B⁻⟩ᵀ 0], since each reciprocated ⊥–⊤ dyad contributes two ordered pairs.

biadjacency_matrices(G::Graphs.SimpleDiGraph; membership, skipchecks::Bool=false)
biadjacency_matrices(m::DBiCM)

Return the pair of biadjacency matrices (B⁺, B⁻) of a directed bipartite graph, both N⊥ × N⊤: B⁺[i,α] = 1 iff ⊥ᵢ → ⊤α, and B⁻[i,α] = 1 iff ⊤α → ⊥ᵢ.

The layer membership is read off the undirected skeleton, since Graphs.bipartite_map traverses outneighbors only and mis-partitions directed input. The DBiCM method uses the model's own stored layers instead, so it works for a model built from degree sequences alone.

project(m::DBiCM; α=0.05, layer=:bottom, kind=:out, distribution=:Poisson, adjustment=BenjaminiHochberg(), multithreaded=false)

Statistically validated monopartite projection of a directed bipartite network onto one of its layers, under the DBiCM null model m.

For every observed pair the count of directed V-motifs of the requested kind is compared with its distribution under the model — exactly Poisson-binomial, or the Poisson approximation — and the resulting upper-tail p-values are corrected for multiple testing before thresholding at α.

:out and :in are symmetric, so only unordered pairs are tested and the result is a SimpleGraph; :path is not, so every ordered pair is tested and the result is a SimpleDiGraph. Note that :path therefore tests roughly twice as many hypotheses, which makes any multiple-testing correction correspondingly more conservative on the same data.

V_motifs(m::DBiCM; layer=:bottom, kind=:out)

Expected total number of directed V-motifs of the requested kind under the DBiCM model m, counted over the same index set as the matrix method (unordered pairs for :out/:in, ordered for :path).

V_motifs(m::DBiCM, i::Int, j::Int; layer=:bottom, kind=:out)

Expected number of directed V-motifs of the requested kind between nodes i and j of layer under the DBiCM model m.

V_PB_parameters(m::DBiCM, i::Int, j::Int; layer=:bottom, kind=:out)

The vector of per-shared-node success probabilities q_α whose Poisson-binomial distribution governs the directed V-motif count between i and j.

Each q_α is a product of two distinct Bernoulli entries of the model, which are independent — of each other and across α — so the Poisson-binomial law is exact, for all three kinds and including the :path diagonal i == j.

reciprocated_degree(m::DBiCM, i::Int)

Expected number of ⊤ nodes that node i of the ⊥ layer both points to and receives from, under the DBiCM model m (layer-local index). This is the diagonal of the :path projection, Σ_α p⁺_iα p⁻_iα, which is exact because B⁺_iα and B⁻_iα are independent.

Vn_motifs(m::DBiCM, n::Int; layer=:bottom, kind=:out, method=:exact)

Expected number of Vn-motifs (n nodes of layer sharing a partner, all in the same direction) under the DBiCM model m. See _Vn_aggregation_classes for why kind cannot be :path.

Vn_sigma(m::DBiCM, n::Int; layer=:bottom, kind=:out, method=:exact)

Standard deviation of the Vn-motif count under the DBiCM model m.

Vn_zscore(m::DBiCM, n::Int; layer=:bottom, kind=:out, method=:exact)

z-score of the observed Vn-motif count against its distribution under the DBiCM model m.

_dbicm_v_accessors(m::DBiCM, layer::Symbol, kind::Symbol)

The two per-entry probability functions whose product defines a V-motif of the requested kind, plus the number of shared nodes summed over. Both are called as f(pair_index, shared_index).

_dbicm_channel_pair(B⁺, B⁻, layer::Symbol, kind::Symbol)

The ordered pair of matrices whose product gives the projection of the requested kind, oriented so that the projection is M₁ * M₂ with rows and columns both indexed by the projected layer.

_Vn_delta_mean(cls::NTuple{5,Any}, n::Int, ::Type{N})
_Vn_delta_sigma(cls::NTuple{5,Any}, n::Int, ::Type{N})

The Saracco et al. (2015) closed forms for the mean and standard deviation of the Vn-motif count, taking the reduced-class tuple so that both the BiCM and one channel of a DBiCM can use them. The mean is only available for n ∈ {2, 3, 4}; method = :exact has no such restriction.

_Vn_aggregation_classes(m::BiCM, layer::Symbol)

Internal helper for the Vn/Λn motif family. Vn-motifs between nodes of layer aggregate over the degrees of the opposite layer: return (own_r, own_f, own_d, opp_r, opp_f) where own_r/own_f/own_d are the reduced parameters, frequencies and reduced degrees of the aggregated (opposite) layer and opp_r/opp_f those generating the connection probabilities (one per node of layer).

_Vn_aggregation_classes(m::DBiCM, layer::Symbol, kind::Symbol)

The reduced-class ingredients of the Vn machinery for one channel of a DBiCM, in the same (own_r, own_f, own_d, opp_r, opp_f) order the BiCM method returns — "own" being the layer that is aggregated over, i.e. the one opposite to layer.

Vn counts n nodes of layer sharing a partner in the same direction, so it is a single-channel quantity and kind is restricted to :out and :in. There is no :path family: a mixed-direction higher-order motif needs the joint law of a shared node's in- and out-degree, which is a two-dimensional Poisson-multinomial over the four dyad states rather than the one-dimensional convolution used here.

outdegree(m::DBiCM, i::Int; method::Symbol=:reduced)
outdegree(m::DBiCM, v::Vector{Int}=collect(1:m.status[:N]); method::Symbol=:reduced)

Expected out-degree of vertex i of the DBiCM model m (graph vertex ids, either layer). A ⊥ vertex sends through channel ⁺, a ⊤ vertex through channel ⁻.

method is :reduced (default, from the reduced parameters), :full (summing the expected entries) or :adjacency (from the stored matrices, which must have been computed with set_Ĝ!).

indegree(m::DBiCM, i::Int; method::Symbol=:reduced)
indegree(m::DBiCM, v::Vector{Int}=collect(1:m.status[:N]); method::Symbol=:reduced)

Expected in-degree of vertex i of the DBiCM model m. A ⊥ vertex receives through channel ⁻, a ⊤ vertex through channel ⁺. See outdegree for method.

degree(m::DBiCM, i::Int; method::Symbol=:reduced)
degree(m::DBiCM, v::Vector{Int}=collect(1:m.status[:N]); method::Symbol=:reduced)

Expected total degree (out plus in) of vertex i of the DBiCM model m.