Skip to main content

sc_neurocore_engine/
connectome.rs

1// SPDX-License-Identifier: AGPL-3.0-or-later | Commercial license available
2// © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
3// © Code 2020–2026 Miroslav Šotek. All rights reserved.
4// ORCID: 0009-0009-3560-0851
5// Contact: www.anulum.li | protoscience@anulum.li
6// SC-NeuroCore — Biologically plausible connectivity generators
7
8//! Biologically plausible connectivity generators.
9
10use rand::{RngExt, SeedableRng};
11use rand_xoshiro::Xoshiro256PlusPlus;
12
13/// Watts-Strogatz small-world network.
14///
15/// Returns adjacency matrix as flat row-major Vec<u8> of shape [n, n].
16pub fn watts_strogatz(n: usize, k: usize, p_rewire: f64, seed: u64) -> Vec<u8> {
17    let mut rng = Xoshiro256PlusPlus::seed_from_u64(seed);
18    let mut adj = vec![0u8; n * n];
19
20    if k >= n {
21        for i in 0..n {
22            for j in 0..n {
23                if i != j {
24                    adj[i * n + j] = 1;
25                }
26            }
27        }
28        return adj;
29    }
30
31    // Ring lattice
32    for i in 0..n {
33        for d in 1..=k / 2 {
34            let j = (i + d) % n;
35            adj[i * n + j] = 1;
36            adj[j * n + i] = 1;
37        }
38    }
39
40    // Rewire
41    for i in 0..n {
42        for d in 1..=k / 2 {
43            let j = (i + d) % n;
44            if rng.random::<f64>() < p_rewire {
45                adj[i * n + j] = 0;
46                let mut new_j = i;
47                while new_j == i || adj[i * n + new_j] == 1 {
48                    new_j = (rng.random::<f64>() * n as f64) as usize % n;
49                }
50                adj[i * n + new_j] = 1;
51            }
52        }
53    }
54
55    adj
56}
57
58#[allow(clippy::needless_range_loop)]
59/// Barabási-Albert scale-free network via preferential attachment.
60///
61/// Returns adjacency matrix as flat row-major Vec<u8> of shape [n, n].
62pub fn barabasi_albert(n: usize, seed: u64) -> Vec<u8> {
63    let mut rng = Xoshiro256PlusPlus::seed_from_u64(seed);
64    let mut adj = vec![0u8; n * n];
65    let mut degree = vec![0u64; n];
66
67    if n < 2 {
68        return adj;
69    }
70
71    // Seed: node 0 ↔ node 1
72    adj[1] = 1;
73    adj[n] = 1;
74    degree[0] = 1;
75    degree[1] = 1;
76
77    for i in 2..n {
78        let total_degree: u64 = degree[..i].iter().sum();
79        if total_degree == 0 {
80            adj[i * n] = 1;
81            degree[i] += 1;
82            degree[0] += 1;
83            continue;
84        }
85
86        let r = rng.random::<f64>() * total_degree as f64;
87        let mut cumulative = 0.0;
88        let mut target = 0;
89        for j in 0..i {
90            cumulative += degree[j] as f64;
91            if cumulative > r {
92                target = j;
93                break;
94            }
95        }
96
97        adj[i * n + target] = 1;
98        degree[i] += 1;
99        degree[target] += 1;
100    }
101
102    adj
103}
104
105#[cfg(test)]
106mod tests {
107    use super::*;
108
109    #[test]
110    fn ws_no_self_loops() {
111        let adj = watts_strogatz(20, 4, 0.3, 42);
112        for i in 0..20 {
113            assert_eq!(adj[i * 20 + i], 0);
114        }
115    }
116
117    #[test]
118    fn ws_connected() {
119        let adj = watts_strogatz(10, 4, 0.0, 0);
120        for i in 0..10 {
121            let row_sum: u8 = (0..10).map(|j| adj[i * 10 + j]).sum();
122            assert!(row_sum > 0, "node {i} has no connections");
123        }
124    }
125
126    #[test]
127    fn barabasi_albert_grows_correctly() {
128        let adj = barabasi_albert(10, 99);
129        for i in 0..10 {
130            assert_eq!(adj[i * 10 + i], 0);
131        }
132        // Each node (after 0,1) adds at least one edge
133        for i in 2..10 {
134            let row_sum: u8 = (0..10).map(|j| adj[i * 10 + j]).sum();
135            assert!(row_sum > 0, "node {i} has no outgoing edges");
136        }
137    }
138
139    #[test]
140    fn deterministic() {
141        let a = watts_strogatz(10, 4, 0.3, 123);
142        let b = watts_strogatz(10, 4, 0.3, 123);
143        assert_eq!(a, b);
144    }
145}