Skip to main content

sc_neurocore_engine/
connectome.rs

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