rust_igraph/algorithms/games/

callaway_traits.rs

//! Callaway-Hopcroft-Kleinberg-Newman-Strogatz traits game (ALGO-GN-016).
//!
//! Counterpart of `igraph_callaway_traits_game()` from
//! `references/igraph/src/games/callaway_traits.c:71`. Models a growing
//! graph with vertex types: at each step a new vertex is added, then
//! `edges_per_step` candidate edges are generated by picking **two
//! uniformly random vertices from those already added** (including the
//! new one) and accepting the edge with probability
//! `pref_matrix[type1][type2]`.
//!
//! ## Algorithm
//!
//! 1. Build a cumulative distribution from `type_dist` (or the uniform
//!    distribution when `type_dist == None`). Reject if `maxcum <= 0`.
//! 2. Pre-assign every vertex `v ∈ [0, nodes)` a categorical type by
//!    drawing a uniform sample in `[0, maxcum)` and binary-searching
//!    the cumulative table.
//! 3. For each step `i ∈ [1, nodes)` and each `_ ∈ [0, edges_per_step)`,
//!    draw `node1, node2` uniformly in `[0, i]` (inclusive — both can
//!    equal `i`, the just-added vertex). Accept the candidate edge
//!    `(node1, node2)` with probability `pref_matrix[type[node1]][type[node2]]`.
//!
//! ## Self-loops & multi-edges
//!
//! Unlike `establishment_game`, this generator can return self-loops
//! (`node1 == node2`) and parallel edges (the same `(node1, node2)`
//! pair drawn twice). The output is therefore **not simple by
//! construction**.
//!
//! ## Determinism
//!
//! Reproducible given the inputs and `seed` against the shared
//! `SplitMix64` PRNG. The stream is **not** portable
//! to upstream igraph's GLIBC RNG, so conformance assertions are
//! structural (vertex/edge counts, type-vector range, support of the
//! preference matrix) rather than bit-exact.
//!
//! ## References
//!
//! * D. S. Callaway, J. E. Hopcroft, J. M. Kleinberg, M. E. J. Newman,
//!   and S. H. Strogatz, *"Are randomly grown graphs really random?"*,
//!   Phys. Rev. E **64**, 041902 (2001).

#![allow(
    unknown_lints,
    clippy::cast_possible_truncation,
    clippy::cast_precision_loss,
    clippy::cast_sign_loss,
    clippy::float_cmp,
    clippy::too_many_arguments,
    clippy::similar_names,
    clippy::many_single_char_names,
    clippy::needless_range_loop
)]

use crate::core::rng::SplitMix64;
use crate::core::{Graph, IgraphError, IgraphResult, VertexId};

fn validate(
    nodes: u32,
    types: u32,
    type_dist: Option<&[f64]>,
    pref_matrix: &[Vec<f64>],
    directed: bool,
) -> IgraphResult<()> {
    if types < 1 {
        return Err(IgraphError::InvalidArgument(
            "The number of vertex types must be at least 1.".into(),
        ));
    }
    let k = types as usize;
    if let Some(td) = type_dist {
        if td.len() != k {
            return Err(IgraphError::InvalidArgument(format!(
                "type_dist length ({}) does not match number of types ({k})",
                td.len()
            )));
        }
        for (i, &v) in td.iter().enumerate() {
            if v.is_nan() {
                return Err(IgraphError::InvalidArgument(format!(
                    "type_dist[{i}] must not be NaN"
                )));
            }
            if v < 0.0 {
                return Err(IgraphError::InvalidArgument(format!(
                    "type_dist[{i}] = {v} must be non-negative"
                )));
            }
        }
    }
    if pref_matrix.len() != k {
        return Err(IgraphError::InvalidArgument(format!(
            "preference matrix has {} rows (expected {k})",
            pref_matrix.len()
        )));
    }
    for (i, row) in pref_matrix.iter().enumerate() {
        if row.len() != k {
            return Err(IgraphError::InvalidArgument(format!(
                "preference matrix row {i} has length {} (expected {k})",
                row.len()
            )));
        }
        for (j, &p) in row.iter().enumerate() {
            if p.is_nan() {
                return Err(IgraphError::InvalidArgument(format!(
                    "preference matrix entry [{i}][{j}] must not be NaN"
                )));
            }
            if !(0.0..=1.0).contains(&p) {
                return Err(IgraphError::InvalidArgument(format!(
                    "preference matrix entry [{i}][{j}] = {p} must lie in [0, 1]"
                )));
            }
        }
    }
    if !directed {
        for (i, row_i) in pref_matrix.iter().enumerate() {
            for (j, row_j) in pref_matrix.iter().enumerate().skip(i + 1) {
                if row_i[j] != row_j[i] {
                    return Err(IgraphError::InvalidArgument(format!(
                        "preference matrix must be symmetric for undirected graphs: \
                         pref[{i}][{j}] = {} but pref[{j}][{i}] = {}",
                        row_i[j], row_j[i],
                    )));
                }
            }
        }
    }
    let _ = nodes;
    Ok(())
}

/// Smallest `k ≥ 1` such that `cumdist[k] > u`. Mirrors
/// `igraph_vector_binsearch` for the cumulative-distribution case.
fn cumdist_lookup(cumdist: &[f64], u: f64) -> usize {
    let mut lo = 1usize;
    let mut hi = cumdist.len();
    while lo < hi {
        let mid = lo + (hi - lo) / 2;
        if cumdist[mid] > u {
            hi = mid;
        } else {
            lo = mid + 1;
        }
    }
    lo.min(cumdist.len() - 1).max(1)
}

/// Generates a random graph with vertex types via the Callaway et al.
/// (2001) growing model.
///
/// At each step `i ∈ [1, nodes)`, `edges_per_step` candidate edges are
/// generated by picking two vertices uniformly at random from
/// `[0, i]` (inclusive of `i` — the vertex just added) and accepting the
/// edge with probability `pref_matrix[type1][type2]`.
///
/// * `nodes` — number of vertices.
/// * `types ≥ 1` — number of vertex types.
/// * `edges_per_step` — number of candidate edges per step.
/// * `type_dist` — categorical weights over types; `None` means uniform.
///   Must be the same length as `types`, all entries non-negative,
///   total mass strictly positive.
/// * `pref_matrix` — `types × types` matrix of acceptance probabilities
///   in `[0, 1]`. When `directed = false` it must be symmetric.
/// * `directed` — generate a directed graph.
/// * `seed` — `SplitMix64` seed.
///
/// Returns the generated [`Graph`] and the per-vertex type vector.
///
/// # Self-loops & multi-edges
///
/// The output is **not** simple by construction: candidate edges allow
/// `node1 == node2` (self-loop) and the same vertex pair may be drawn
/// multiple times (multi-edge).
///
/// # Errors
///
/// Returns [`IgraphError::InvalidArgument`] when the preference matrix
/// or type distribution is malformed (wrong size, NaN, out-of-range
/// probability, asymmetric for an undirected graph) or when the type
/// distribution is non-positive everywhere.
///
/// # Example
///
/// ```
/// use rust_igraph::callaway_traits_game;
///
/// // Two types, only-cross prefs ⇒ bipartite-flavoured directed graph.
/// let pref = vec![vec![0.0, 1.0], vec![1.0, 0.0]];
/// let (g, types) =
///     callaway_traits_game(20, 2, 3, Some(&[1.0, 1.0]), &pref, true, 42).unwrap();
/// assert_eq!(g.vcount(), 20);
/// assert_eq!(types.len(), 20);
/// assert!(types.iter().all(|&t| t < 2));
/// ```
pub fn callaway_traits_game(
    nodes: u32,
    types: u32,
    edges_per_step: u32,
    type_dist: Option<&[f64]>,
    pref_matrix: &[Vec<f64>],
    directed: bool,
    seed: u64,
) -> IgraphResult<(Graph, Vec<u32>)> {
    validate(nodes, types, type_dist, pref_matrix, directed)?;

    let n_types = types as usize;
    let mut rng = SplitMix64::new(seed);

    let mut cumdist = vec![0.0f64; n_types + 1];
    if let Some(td) = type_dist {
        for i in 0..n_types {
            cumdist[i + 1] = cumdist[i] + td[i];
        }
    } else {
        for i in 0..n_types {
            cumdist[i + 1] = (i + 1) as f64;
        }
    }
    let maxcum = cumdist[n_types];
    if maxcum <= 0.0 {
        return Err(IgraphError::InvalidArgument(
            "type_dist must contain at least one positive value".into(),
        ));
    }

    let mut node_types = vec![0u32; nodes as usize];
    for v in 0..(nodes as usize) {
        let u = rng.gen_unit() * maxcum;
        let pos = cumdist_lookup(&cumdist, u);
        let t = (pos - 1).min(n_types - 1);
        node_types[v] = t as u32;
    }

    let mut edges: Vec<(VertexId, VertexId)> = Vec::new();
    if edges_per_step > 0 && nodes >= 2 {
        for i in 1..nodes {
            // Inclusive upper bound: node1, node2 ∈ [0, i].
            let span = (i as usize) + 1;
            for _ in 0..edges_per_step {
                let n1 = rng.gen_index(span) as u32;
                let n2 = rng.gen_index(span) as u32;
                let t1 = node_types[n1 as usize] as usize;
                let t2 = node_types[n2 as usize] as usize;
                let p = pref_matrix[t1][t2];
                if p > 0.0 && rng.gen_unit() < p {
                    edges.push((n1, n2));
                }
            }
        }
    }

    let mut g = Graph::new(nodes, directed)?;
    g.add_edges(edges)?;
    Ok((g, node_types))
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::collections::HashSet as Set;

    fn diag_pref(types: usize, p: f64) -> Vec<Vec<f64>> {
        (0..types)
            .map(|i| (0..types).map(|j| if i == j { p } else { 0.0 }).collect())
            .collect()
    }

    fn full_pref(types: usize, p: f64) -> Vec<Vec<f64>> {
        vec![vec![p; types]; types]
    }

    #[test]
    fn nodes_zero_returns_empty_graph() {
        let pref = full_pref(2, 0.5);
        let (g, types) = callaway_traits_game(0, 2, 5, None, &pref, false, 42).unwrap();
        assert_eq!(g.vcount(), 0);
        assert_eq!(g.ecount(), 0);
        assert_eq!(types.len(), 0);
        assert!(!g.is_directed());
    }

    #[test]
    fn nodes_one_returns_no_edges() {
        let pref = full_pref(2, 1.0);
        let (g, types) = callaway_traits_game(1, 2, 10, None, &pref, false, 7).unwrap();
        assert_eq!(g.vcount(), 1);
        assert_eq!(g.ecount(), 0);
        assert_eq!(types.len(), 1);
    }

    #[test]
    fn edges_per_step_zero_yields_no_edges() {
        let pref = full_pref(2, 1.0);
        let (g, _) = callaway_traits_game(50, 2, 0, None, &pref, false, 11).unwrap();
        assert_eq!(g.vcount(), 50);
        assert_eq!(g.ecount(), 0);
    }

    #[test]
    fn zero_pref_matrix_yields_no_edges() {
        let pref = full_pref(3, 0.0);
        let (g, types) = callaway_traits_game(40, 3, 5, None, &pref, false, 99).unwrap();
        assert_eq!(g.vcount(), 40);
        assert_eq!(g.ecount(), 0);
        assert!(types.iter().all(|&t| t < 3));
    }

    #[test]
    fn full_pref_p1_undirected_max_ecount() {
        // p=1 ⇒ every candidate edge accepted ⇒ exactly
        // (nodes - 1) * edges_per_step edges (with possible self-loops/multi).
        let pref = full_pref(2, 1.0);
        let nodes = 30u32;
        let eps = 4u32;
        let (g, _) = callaway_traits_game(nodes, 2, eps, None, &pref, false, 123).unwrap();
        assert_eq!(g.ecount() as u32, (nodes - 1) * eps);
    }

    #[test]
    fn full_pref_p1_directed_max_ecount() {
        let pref = vec![vec![1.0, 1.0], vec![1.0, 1.0]];
        let nodes = 25u32;
        let eps = 3u32;
        let (g, _) = callaway_traits_game(nodes, 2, eps, None, &pref, true, 456).unwrap();
        assert_eq!(g.ecount() as u32, (nodes - 1) * eps);
        assert!(g.is_directed());
    }

    #[test]
    fn determinism_same_seed_same_graph() {
        let pref = full_pref(3, 0.4);
        let (g1, t1) = callaway_traits_game(50, 3, 5, None, &pref, false, 0xDEAD).unwrap();
        let (g2, t2) = callaway_traits_game(50, 3, 5, None, &pref, false, 0xDEAD).unwrap();
        assert_eq!(g1.ecount(), g2.ecount());
        assert_eq!(t1, t2);
        for eid in 0..g1.ecount() as u32 {
            assert_eq!(g1.edge(eid).unwrap(), g2.edge(eid).unwrap());
        }
    }

    #[test]
    fn different_seeds_diverge() {
        let pref = full_pref(3, 0.4);
        let (g1, _) = callaway_traits_game(80, 3, 5, None, &pref, false, 1).unwrap();
        let (g2, _) = callaway_traits_game(80, 3, 5, None, &pref, false, 2).unwrap();
        // With p=0.4 over many candidates, edge counts almost certainly differ.
        let differ = g1.ecount() != g2.ecount() || {
            let mut e1: Vec<_> = (0..g1.ecount() as u32)
                .map(|e| g1.edge(e).unwrap())
                .collect();
            let mut e2: Vec<_> = (0..g2.ecount() as u32)
                .map(|e| g2.edge(e).unwrap())
                .collect();
            e1.sort_unstable();
            e2.sort_unstable();
            e1 != e2
        };
        assert!(differ);
    }

    #[test]
    fn type_dist_skewed_one_hot() {
        // Only type 0 has positive mass ⇒ every vertex is type 0.
        let pref = full_pref(3, 0.0);
        let dist = vec![1.0, 0.0, 0.0];
        let (_, types) = callaway_traits_game(40, 3, 4, Some(&dist), &pref, false, 12).unwrap();
        assert!(types.iter().all(|&t| t == 0));
    }

    #[test]
    fn type_dist_zero_everywhere_errors() {
        let pref = full_pref(2, 0.5);
        let dist = vec![0.0, 0.0];
        let err = callaway_traits_game(10, 2, 3, Some(&dist), &pref, false, 1).unwrap_err();
        assert!(matches!(err, IgraphError::InvalidArgument(_)));
    }

    #[test]
    fn types_zero_errors() {
        let pref: Vec<Vec<f64>> = Vec::new();
        let err = callaway_traits_game(10, 0, 3, None, &pref, false, 1).unwrap_err();
        assert!(matches!(err, IgraphError::InvalidArgument(_)));
    }

    #[test]
    fn pref_wrong_shape_errors() {
        let pref = vec![vec![0.5, 0.5]]; // 1 row, types=2
        let err = callaway_traits_game(10, 2, 3, None, &pref, false, 1).unwrap_err();
        assert!(matches!(err, IgraphError::InvalidArgument(_)));
    }

    #[test]
    fn pref_row_wrong_length_errors() {
        let pref = vec![vec![0.5, 0.5], vec![0.5]];
        let err = callaway_traits_game(10, 2, 3, None, &pref, false, 1).unwrap_err();
        assert!(matches!(err, IgraphError::InvalidArgument(_)));
    }

    #[test]
    fn pref_nan_errors() {
        let pref = vec![vec![0.5, f64::NAN], vec![f64::NAN, 0.5]];
        let err = callaway_traits_game(10, 2, 3, None, &pref, false, 1).unwrap_err();
        assert!(matches!(err, IgraphError::InvalidArgument(_)));
    }

    #[test]
    fn pref_out_of_range_errors() {
        let pref = vec![vec![0.5, 1.5], vec![1.5, 0.5]];
        let err = callaway_traits_game(10, 2, 3, None, &pref, false, 1).unwrap_err();
        assert!(matches!(err, IgraphError::InvalidArgument(_)));
    }

    #[test]
    fn pref_negative_errors() {
        let pref = vec![vec![0.5, -0.1], vec![-0.1, 0.5]];
        let err = callaway_traits_game(10, 2, 3, None, &pref, false, 1).unwrap_err();
        assert!(matches!(err, IgraphError::InvalidArgument(_)));
    }

    #[test]
    fn pref_asymmetric_undirected_errors() {
        let pref = vec![vec![0.5, 0.2], vec![0.7, 0.5]];
        let err = callaway_traits_game(10, 2, 3, None, &pref, false, 1).unwrap_err();
        assert!(matches!(err, IgraphError::InvalidArgument(_)));
    }

    #[test]
    fn pref_asymmetric_directed_ok() {
        let pref = vec![vec![0.5, 0.2], vec![0.7, 0.5]];
        let (g, _) = callaway_traits_game(20, 2, 3, None, &pref, true, 1).unwrap();
        assert_eq!(g.vcount(), 20);
        assert!(g.is_directed());
    }

    #[test]
    fn type_dist_nan_errors() {
        let pref = full_pref(2, 0.5);
        let dist = vec![1.0, f64::NAN];
        let err = callaway_traits_game(10, 2, 3, Some(&dist), &pref, false, 1).unwrap_err();
        assert!(matches!(err, IgraphError::InvalidArgument(_)));
    }

    #[test]
    fn type_dist_negative_errors() {
        let pref = full_pref(2, 0.5);
        let dist = vec![1.0, -0.1];
        let err = callaway_traits_game(10, 2, 3, Some(&dist), &pref, false, 1).unwrap_err();
        assert!(matches!(err, IgraphError::InvalidArgument(_)));
    }

    #[test]
    fn type_dist_wrong_length_errors() {
        let pref = full_pref(3, 0.5);
        let dist = vec![1.0, 1.0]; // length 2, types=3
        let err = callaway_traits_game(10, 3, 3, Some(&dist), &pref, false, 1).unwrap_err();
        assert!(matches!(err, IgraphError::InvalidArgument(_)));
    }

    #[test]
    fn diag_only_endpoints_share_type_when_edge_present() {
        // When pref is purely diagonal, any accepted edge connects
        // two vertices of the same type.
        let pref = diag_pref(3, 1.0);
        let (g, types) = callaway_traits_game(60, 3, 4, None, &pref, false, 7777).unwrap();
        for eid in 0..g.ecount() as u32 {
            let (s, d) = g.edge(eid).unwrap();
            assert_eq!(types[s as usize], types[d as usize]);
        }
    }

    #[test]
    fn cross_only_endpoints_differ_when_edge_present() {
        // pref is zero on the diagonal ⇒ no same-type edges
        // (excluding self-loops, which can still happen with prob 0
        // since pref[t][t] = 0 — so they cannot occur).
        let pref = vec![vec![0.0, 1.0], vec![1.0, 0.0]];
        let (g, types) = callaway_traits_game(60, 2, 4, None, &pref, false, 8888).unwrap();
        for eid in 0..g.ecount() as u32 {
            let (s, d) = g.edge(eid).unwrap();
            assert_ne!(types[s as usize], types[d as usize]);
        }
    }

    #[test]
    fn type_range_is_subset_of_types() {
        let pref = full_pref(4, 0.3);
        let (_, types) = callaway_traits_game(100, 4, 3, None, &pref, false, 314).unwrap();
        let observed: Set<u32> = types.iter().copied().collect();
        assert!(observed.iter().all(|&t| t < 4));
    }

    #[test]
    fn directed_p1_ecount_matches_formula() {
        // Directed, full-1 pref, type_dist = [1, 1, 1] uniform:
        // ecount == (nodes - 1) * edges_per_step regardless of types.
        let pref = vec![vec![1.0; 3]; 3];
        let nodes = 40u32;
        let eps = 2u32;
        let (g, _) =
            callaway_traits_game(nodes, 3, eps, Some(&[1.0, 1.0, 1.0]), &pref, true, 555).unwrap();
        assert_eq!(g.ecount() as u32, (nodes - 1) * eps);
    }
}

#[cfg(all(test, feature = "proptest-harness"))]
mod proptest_invariants {
    use super::*;
    use proptest::prelude::*;

    proptest! {
        #![proptest_config(ProptestConfig::with_cases(32))]

        #[test]
        fn ecount_bounded_by_full_accept(
            nodes in 0u32..=80,
            types in 1u32..=5,
            eps in 0u32..=6,
            seed in any::<u64>(),
            directed in any::<bool>(),
        ) {
            let t = types as usize;
            // Use a uniform pref matrix at p=0.5; random pref values
            // would over-complicate the bounds proof.
            let pref = vec![vec![0.5; t]; t];
            let (g, types_v) = callaway_traits_game(
                nodes, types, eps, None, &pref, directed, seed,
            ).unwrap();
            prop_assert_eq!(g.vcount(), nodes);
            prop_assert_eq!(types_v.len(), nodes as usize);
            // Max possible: (n-1) * eps when nodes >= 1.
            let max_e = if nodes >= 1 { (nodes - 1) * eps } else { 0 };
            prop_assert!((g.ecount() as u32) <= max_e);
        }

        #[test]
        fn types_in_range(
            nodes in 0u32..=60,
            types in 1u32..=5,
            seed in any::<u64>(),
        ) {
            let t = types as usize;
            let pref = vec![vec![0.0; t]; t];
            let (_, types_v) = callaway_traits_game(
                nodes, types, 3, None, &pref, false, seed,
            ).unwrap();
            for &x in &types_v {
                prop_assert!(x < types);
            }
        }

        #[test]
        fn determinism(
            nodes in 0u32..=60,
            types in 1u32..=4,
            eps in 0u32..=4,
            seed in any::<u64>(),
        ) {
            let t = types as usize;
            let pref = vec![vec![0.3; t]; t];
            let (g1, t1) = callaway_traits_game(
                nodes, types, eps, None, &pref, false, seed,
            ).unwrap();
            let (g2, t2) = callaway_traits_game(
                nodes, types, eps, None, &pref, false, seed,
            ).unwrap();
            prop_assert_eq!(g1.ecount(), g2.ecount());
            prop_assert_eq!(t1, t2);
        }

        #[test]
        fn p1_full_pref_yields_exact_max_ecount(
            nodes in 1u32..=50,
            types in 1u32..=4,
            eps in 0u32..=5,
            seed in any::<u64>(),
            directed in any::<bool>(),
        ) {
            let t = types as usize;
            let pref = vec![vec![1.0; t]; t];
            let (g, _) = callaway_traits_game(
                nodes, types, eps, None, &pref, directed, seed,
            ).unwrap();
            prop_assert_eq!(g.ecount() as u32, (nodes - 1) * eps);
        }

        #[test]
        fn p0_yields_no_edges(
            nodes in 0u32..=50,
            types in 1u32..=4,
            eps in 0u32..=5,
            seed in any::<u64>(),
            directed in any::<bool>(),
        ) {
            let t = types as usize;
            let pref = vec![vec![0.0; t]; t];
            let (g, _) = callaway_traits_game(
                nodes, types, eps, None, &pref, directed, seed,
            ).unwrap();
            prop_assert_eq!(g.ecount(), 0);
        }
    }
}