#![allow(
unknown_lints,
clippy::cast_possible_truncation,
clippy::cast_precision_loss,
clippy::cast_sign_loss,
clippy::float_cmp,
clippy::manual_midpoint
)]
use std::collections::HashSet;
use crate::core::rng::SplitMix64;
use crate::core::{Graph, IgraphError, IgraphResult, VertexId};
fn decode_pair(idx: u64, n: u32, directed: bool, loops: bool) -> (VertexId, VertexId) {
let n_u64 = u64::from(n);
if directed && loops {
let to = idx / n_u64;
let from = idx - to * n_u64;
debug_assert!(from < n_u64 && to < n_u64);
#[allow(clippy::cast_possible_truncation)]
((from as VertexId), (to as VertexId))
} else if directed && !loops {
let to = idx / n_u64;
let from = idx - to * n_u64;
let to = if from == to { n_u64 - 1 } else { to };
debug_assert!(from < n_u64 && to < n_u64 && from != to);
#[allow(clippy::cast_possible_truncation)]
((from as VertexId), (to as VertexId))
} else if !directed && loops {
#[allow(clippy::cast_precision_loss)]
let idx_f = idx as f64;
let to_f = ((8.0 * idx_f + 1.0).sqrt() - 1.0) / 2.0;
#[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
let mut to = to_f.trunc() as u64;
let mut from = idx - to * (to + 1) / 2;
while from > to {
to += 1;
from = idx - to * (to + 1) / 2;
}
debug_assert!(from <= to && to < n_u64);
#[allow(clippy::cast_possible_truncation)]
((from as VertexId), (to as VertexId))
} else {
#[allow(clippy::cast_precision_loss)]
let idx_f = idx as f64;
let to_f = ((8.0 * idx_f + 1.0).sqrt() + 1.0) / 2.0;
#[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
let mut to = to_f.trunc() as u64;
if to < 1 {
to = 1;
}
let mut from = idx - to * (to - 1) / 2;
while from >= to {
to += 1;
from = idx - to * (to - 1) / 2;
}
debug_assert!(from < to && to < n_u64);
#[allow(clippy::cast_possible_truncation)]
((from as VertexId), (to as VertexId))
}
}
fn max_edges(n: u32, directed: bool, loops: bool) -> u64 {
let n = u64::from(n);
match (directed, loops) {
(true, true) => n * n,
(true, false) => n * n.saturating_sub(1),
(false, true) => n * (n + 1) / 2,
(false, false) => n * n.saturating_sub(1) / 2,
}
}
fn distinct_sample(rng: &mut SplitMix64, n_pairs: u64, m: u64) -> Vec<u64> {
debug_assert!(m <= n_pairs);
let cap = usize::try_from(m).unwrap_or(usize::MAX);
let mut chosen: HashSet<u64> = HashSet::with_capacity(cap);
let mut out: Vec<u64> = Vec::with_capacity(cap);
for j in (n_pairs - m)..n_pairs {
let span = j + 1;
let span_usize = usize::try_from(span).unwrap_or(usize::MAX);
let t_usize = rng.gen_index(span_usize);
let t = t_usize as u64;
let pick = if chosen.insert(t) {
t
} else {
chosen.insert(j);
j
};
out.push(pick);
}
out.sort_unstable();
out
}
fn validate_gnp(p: f64) -> IgraphResult<()> {
if !p.is_finite() {
return Err(IgraphError::InvalidArgument(format!(
"G(n,p) probability must be finite (got {p})"
)));
}
if !(0.0..=1.0).contains(&p) {
return Err(IgraphError::InvalidArgument(format!(
"G(n,p) probability must be in [0, 1] (got {p})"
)));
}
Ok(())
}
fn validate_gnm(n: u32, m: u64, directed: bool, loops: bool) -> IgraphResult<()> {
let cap = max_edges(n, directed, loops);
if m > cap {
return Err(IgraphError::InvalidArgument(format!(
"G(n,m) requested {m} edges but the {} graph on {n} vertices admits at most {cap}",
if directed { "directed" } else { "undirected" }
)));
}
Ok(())
}
fn finalize(n: u32, directed: bool, edges: Vec<(VertexId, VertexId)>) -> IgraphResult<Graph> {
let mut g = Graph::new(n, directed)?;
g.add_edges(edges)?;
Ok(g)
}
fn complete(n: u32, directed: bool, loops: bool) -> IgraphResult<Graph> {
let cap = max_edges(n, directed, loops);
let cap_usize = usize::try_from(cap)
.map_err(|_| IgraphError::Internal("complete graph edge count exceeds usize"))?;
let mut edges = Vec::with_capacity(cap_usize);
for i in 0..n {
if directed {
for j in 0..n {
if i == j && !loops {
continue;
}
edges.push((i, j));
}
} else {
let start = if loops { i } else { i + 1 };
for j in start..n {
edges.push((i, j));
}
}
}
finalize(n, directed, edges)
}
pub fn erdos_renyi_gnp(
n: u32,
p: f64,
directed: bool,
loops: bool,
seed: u64,
) -> IgraphResult<Graph> {
validate_gnp(p)?;
let cap = max_edges(n, directed, loops);
if n == 0 || p == 0.0 || cap == 0 {
return Graph::new(n, directed);
}
if p == 1.0 {
return complete(n, directed, loops);
}
let mut rng = SplitMix64::new(seed);
#[allow(clippy::cast_precision_loss)]
let cap_f = cap as f64;
let mut last = rng.gen_geom(p);
let mut indices: Vec<u64> = Vec::new();
while last < cap_f {
#[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
let idx = last.trunc() as u64;
if idx >= cap {
break;
}
indices.push(idx);
last += rng.gen_geom(p);
last += 1.0; }
let edges: Vec<(VertexId, VertexId)> = indices
.into_iter()
.map(|idx| decode_pair(idx, n, directed, loops))
.collect();
finalize(n, directed, edges)
}
pub fn erdos_renyi_gnm(
n: u32,
m: u64,
directed: bool,
loops: bool,
seed: u64,
) -> IgraphResult<Graph> {
validate_gnm(n, m, directed, loops)?;
if m == 0 {
return Graph::new(n, directed);
}
let cap = max_edges(n, directed, loops);
if m == cap {
return complete(n, directed, loops);
}
let mut rng = SplitMix64::new(seed);
let picks = distinct_sample(&mut rng, cap, m);
let edges: Vec<(VertexId, VertexId)> = picks
.into_iter()
.map(|idx| decode_pair(idx, n, directed, loops))
.collect();
finalize(n, directed, edges)
}
#[cfg(test)]
mod tests {
use super::*;
use std::collections::HashSet;
#[test]
fn decode_undirected_no_loops_covers_strict_upper_triangle() {
let n = 5u32;
let cap = max_edges(n, false, false);
assert_eq!(cap, 10);
let mut seen: HashSet<(u32, u32)> = HashSet::new();
for idx in 0..cap {
let (u, v) = decode_pair(idx, n, false, false);
assert!(u < v && v < n);
assert!(seen.insert((u, v)), "duplicate pair at idx {idx}");
}
assert_eq!(seen.len(), cap as usize);
}
#[test]
fn decode_undirected_with_loops_covers_lower_triangle_incl_diagonal() {
let n = 4u32;
let cap = max_edges(n, false, true);
assert_eq!(cap, 10);
let mut seen: HashSet<(u32, u32)> = HashSet::new();
for idx in 0..cap {
let (u, v) = decode_pair(idx, n, false, true);
assert!(u <= v && v < n);
assert!(seen.insert((u, v)));
}
assert_eq!(seen.len(), cap as usize);
}
#[test]
fn decode_directed_no_loops_covers_off_diagonal() {
let n = 5u32;
let cap = max_edges(n, true, false);
assert_eq!(cap, 20);
let mut seen: HashSet<(u32, u32)> = HashSet::new();
for idx in 0..cap {
let (u, v) = decode_pair(idx, n, true, false);
assert!(u < n && v < n && u != v);
assert!(seen.insert((u, v)));
}
assert_eq!(seen.len(), cap as usize);
}
#[test]
fn decode_directed_with_loops_covers_full_grid() {
let n = 4u32;
let cap = max_edges(n, true, true);
assert_eq!(cap, 16);
let mut seen: HashSet<(u32, u32)> = HashSet::new();
for idx in 0..cap {
let (u, v) = decode_pair(idx, n, true, true);
assert!(u < n && v < n);
assert!(seen.insert((u, v)));
}
assert_eq!(seen.len(), cap as usize);
}
#[test]
fn gnp_p_zero_is_empty() {
let g = erdos_renyi_gnp(10, 0.0, false, false, 1).unwrap();
assert_eq!(g.vcount(), 10);
assert_eq!(g.ecount(), 0);
}
#[test]
fn gnp_p_one_is_complete() {
let g = erdos_renyi_gnp(6, 1.0, false, false, 1).unwrap();
assert_eq!(g.vcount(), 6);
assert_eq!(g.ecount(), 15); }
#[test]
fn gnp_directed_p_one_loops_is_full_grid() {
let g = erdos_renyi_gnp(4, 1.0, true, true, 1).unwrap();
assert_eq!(g.ecount(), 16);
}
#[test]
fn gnp_invalid_p_rejected() {
assert!(erdos_renyi_gnp(5, -0.1, false, false, 1).is_err());
assert!(erdos_renyi_gnp(5, 1.1, false, false, 1).is_err());
assert!(erdos_renyi_gnp(5, f64::NAN, false, false, 1).is_err());
assert!(erdos_renyi_gnp(5, f64::INFINITY, false, false, 1).is_err());
}
#[test]
fn gnp_zero_vertices_is_empty() {
let g = erdos_renyi_gnp(0, 0.5, false, false, 1).unwrap();
assert_eq!(g.vcount(), 0);
assert_eq!(g.ecount(), 0);
}
#[test]
fn gnp_deterministic_with_seed() {
let a = erdos_renyi_gnp(50, 0.3, false, false, 12345).unwrap();
let b = erdos_renyi_gnp(50, 0.3, false, false, 12345).unwrap();
assert_eq!(a.vcount(), b.vcount());
assert_eq!(a.ecount(), b.ecount());
let edges_a: Vec<_> = (0..a.ecount()).map(|e| a.edge(e as u32).unwrap()).collect();
let edges_b: Vec<_> = (0..b.ecount()).map(|e| b.edge(e as u32).unwrap()).collect();
assert_eq!(edges_a, edges_b);
}
#[test]
fn gnp_different_seeds_differ() {
let a = erdos_renyi_gnp(200, 0.05, false, false, 1).unwrap();
let b = erdos_renyi_gnp(200, 0.05, false, false, 2).unwrap();
assert_ne!(a.ecount(), b.ecount());
}
#[test]
fn gnp_expected_ecount_in_band() {
let g = erdos_renyi_gnp(100, 0.2, false, false, 31_415).unwrap();
let m = g.ecount();
assert!(m > 890 && m < 1090, "ecount = {m}");
}
#[test]
fn gnp_is_simple_no_loops() {
let g = erdos_renyi_gnp(50, 0.4, false, false, 7).unwrap();
for e in 0..g.ecount() {
let (u, v) = g.edge(e as u32).unwrap();
assert_ne!(u, v, "self-loop generated at edge {e}");
}
}
#[test]
fn gnp_no_parallel_edges_simple() {
let g = erdos_renyi_gnp(40, 0.3, false, false, 11).unwrap();
let mut seen: HashSet<(u32, u32)> = HashSet::new();
for e in 0..g.ecount() {
let pair = g.edge(e as u32).unwrap();
assert!(seen.insert(pair), "parallel edge at {e}: {pair:?}");
}
}
#[test]
fn gnm_m_zero_is_empty() {
let g = erdos_renyi_gnm(10, 0, false, false, 1).unwrap();
assert_eq!(g.ecount(), 0);
}
#[test]
fn gnm_m_max_is_complete() {
let g = erdos_renyi_gnm(6, 15, false, false, 1).unwrap();
assert_eq!(g.ecount(), 15);
}
#[test]
fn gnm_m_exceeds_capacity_rejected() {
assert!(erdos_renyi_gnm(5, 11, false, false, 1).is_err());
}
#[test]
fn gnm_exact_edge_count() {
let g = erdos_renyi_gnm(50, 100, false, false, 42).unwrap();
assert_eq!(g.ecount(), 100);
assert_eq!(g.vcount(), 50);
}
#[test]
fn gnm_deterministic_with_seed() {
let a = erdos_renyi_gnm(30, 60, false, false, 7).unwrap();
let b = erdos_renyi_gnm(30, 60, false, false, 7).unwrap();
let edges_a: Vec<_> = (0..a.ecount()).map(|e| a.edge(e as u32).unwrap()).collect();
let edges_b: Vec<_> = (0..b.ecount()).map(|e| b.edge(e as u32).unwrap()).collect();
assert_eq!(edges_a, edges_b);
}
#[test]
fn gnm_is_simple_no_loops() {
let g = erdos_renyi_gnm(20, 40, false, false, 99).unwrap();
let mut seen: HashSet<(u32, u32)> = HashSet::new();
for e in 0..g.ecount() {
let (u, v) = g.edge(e as u32).unwrap();
assert_ne!(u, v);
assert!(seen.insert((u, v)), "parallel edge: {u} {v}");
}
}
#[test]
fn gnm_directed_loops_full_grid() {
let g = erdos_renyi_gnm(3, 9, true, true, 0).unwrap();
assert_eq!(g.ecount(), 9);
let mut seen: HashSet<(u32, u32)> = HashSet::new();
for e in 0..g.ecount() {
seen.insert(g.edge(e as u32).unwrap());
}
assert_eq!(seen.len(), 9);
}
#[test]
fn gnm_directed_no_loops() {
let g = erdos_renyi_gnm(4, 12, true, false, 0).unwrap();
assert_eq!(g.ecount(), 12);
let mut seen: HashSet<(u32, u32)> = HashSet::new();
for e in 0..g.ecount() {
let (u, v) = g.edge(e as u32).unwrap();
assert_ne!(u, v);
seen.insert((u, v));
}
assert_eq!(seen.len(), 12);
}
#[cfg(all(test, feature = "proptest-harness"))]
mod prop {
use super::super::*;
use proptest::prelude::*;
proptest! {
#[test]
fn gnp_vcount_always_matches_n(
n in 0u32..50,
p in 0.0..=1.0,
directed in any::<bool>(),
loops in any::<bool>(),
seed in any::<u64>(),
) {
let g = erdos_renyi_gnp(n, p, directed, loops, seed).unwrap();
prop_assert_eq!(g.vcount(), n);
prop_assert!(g.ecount() as u64 <= max_edges(n, directed, loops));
}
#[test]
fn gnp_simple_no_self_loops_no_multi(
n in 1u32..30,
p in 0.0..=0.7,
seed in any::<u64>(),
) {
let g = erdos_renyi_gnp(n, p, false, false, seed).unwrap();
let mut seen = std::collections::HashSet::new();
for e in 0..g.ecount() {
let (u, v) = g.edge(e as u32).unwrap();
prop_assert!(u != v);
prop_assert!(seen.insert((u, v)));
}
}
#[test]
fn gnm_exact_count_and_simple(
n in 2u32..30,
seed in any::<u64>(),
) {
let cap = max_edges(n, false, false);
if cap == 0 { return Ok(()); }
let m = (seed as u64) % cap;
let g = erdos_renyi_gnm(n, m, false, false, seed).unwrap();
prop_assert_eq!(g.ecount(), m as usize);
let mut seen = std::collections::HashSet::new();
for e in 0..g.ecount() {
let (u, v) = g.edge(e as u32).unwrap();
prop_assert!(u != v);
prop_assert!(seen.insert((u, v)));
}
}
#[test]
fn gnp_determinism(
n in 1u32..40,
p in 0.05..0.95,
directed in any::<bool>(),
loops in any::<bool>(),
seed in any::<u64>(),
) {
let a = erdos_renyi_gnp(n, p, directed, loops, seed).unwrap();
let b = erdos_renyi_gnp(n, p, directed, loops, seed).unwrap();
prop_assert_eq!(a.ecount(), b.ecount());
for e in 0..a.ecount() {
prop_assert_eq!(a.edge(e as u32).unwrap(), b.edge(e as u32).unwrap());
}
}
}
}
}