sc_neurocore_engine/
connectome.rs1use rand::{RngExt, SeedableRng};
11use rand_xoshiro::Xoshiro256PlusPlus;
12
13pub 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 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 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)]
59pub 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 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 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}