sc_neurocore_engine 3.15.0

// SPDX-License-Identifier: AGPL-3.0-or-later
// Commercial license available
// © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
// © Code 2020–2026 Miroslav Šotek. All rights reserved.
// ORCID: 0009-0009-3560-0851
// Contact: www.anulum.li | protoscience@anulum.li
// SC-NeuroCore — Rust N-step simulator for the Wong-Wang 2006 decision unit

//! Batch parity with `WongWangUnit.step` in
//! `src/sc_neurocore/neurons/models/wong_wang.py` (Wong & Wang 2006,
//! J. Neurosci. 26:1314–1328).
//!
//! Per step:
//!   1. `i_k = j_n * s_k - j_cross * s_(3-k) + i_0 + stim_k + sigma * xi`
//!   2. `r_k = phi(i_k)` where
//!      `phi(i) = (a*i - b) / (1 - exp(-d*(a*i - b)))`
//!      with singularity guard `|a*i - b| < 1e-6 -> 1/d`.
//!   3. `s_k += (-s_k/tau_s + (1 - s_k) * gamma * r_k) * dt`
//!   4. clamp `s_k` into `[0, 1]`.
//!
//! The Python primary draws `np.random.randn()` twice per step; the
//! Rust simulator takes `xi` pre-drawn from the Python RNG so trajectories
//! are bit-exact for matching seeds. This mirrors the ping.rs +
//! PINGCircuit pattern: Python owns the RNG, Rust owns the inner loop.

const A: f64 = 270.0;
const B: f64 = 108.0;
const D: f64 = 0.154;

#[inline]
fn phi(i_syn: f64) -> f64 {
    let x = A * i_syn - B;
    if x.abs() < 1e-6 {
        1.0 / D
    } else {
        x / (1.0 - (-D * x).exp())
    }
}

/// Simulate `n_steps` Wong-Wang iterations; write per-step state +
/// firing-rate traces. `xi` must be length `2 * n_steps` (two noise
/// samples per step, consumed in `i1, i2` order).
#[allow(clippy::too_many_arguments)]
pub fn simulate(
    mut s1: f64,
    mut s2: f64,
    tau_s: f64,
    gamma: f64,
    j_n: f64,
    j_cross: f64,
    i_0: f64,
    sigma: f64,
    dt: f64,
    stim1: &[f64],
    stim2: &[f64],
    xi: &[f64],
    s1_out: &mut [f64],
    s2_out: &mut [f64],
    r1_out: &mut [f64],
    r2_out: &mut [f64],
) -> (f64, f64) {
    let n = stim1.len();
    assert_eq!(stim2.len(), n, "stim2 length mismatch");
    assert_eq!(xi.len(), 2 * n, "xi length must be 2 * n_steps");
    assert_eq!(s1_out.len(), n, "s1_out length mismatch");
    assert_eq!(s2_out.len(), n, "s2_out length mismatch");
    assert_eq!(r1_out.len(), n, "r1_out length mismatch");
    assert_eq!(r2_out.len(), n, "r2_out length mismatch");

    for t in 0..n {
        let xi1 = xi[2 * t];
        let xi2 = xi[2 * t + 1];
        let i1 = j_n * s1 - j_cross * s2 + i_0 + stim1[t] + sigma * xi1;
        let i2 = j_n * s2 - j_cross * s1 + i_0 + stim2[t] + sigma * xi2;
        let r1 = phi(i1);
        let r2 = phi(i2);
        s1 += (-s1 / tau_s + (1.0 - s1) * gamma * r1) * dt;
        s2 += (-s2 / tau_s + (1.0 - s2) * gamma * r2) * dt;
        s1 = s1.clamp(0.0, 1.0);
        s2 = s2.clamp(0.0, 1.0);
        s1_out[t] = s1;
        s2_out[t] = s2;
        r1_out[t] = r1;
        r2_out[t] = r2;
    }
    (s1, s2)
}

#[cfg(test)]
mod tests {
    use super::*;

    fn params() -> (f64, f64, f64, f64, f64, f64, f64) {
        // tau_s, gamma, j_n, j_cross, i_0, sigma, dt
        (0.1, 0.641, 0.2609, 0.0497, 0.3255, 0.0, 0.001)
    }

    #[test]
    fn phi_singularity_guard_returns_finite() {
        // a*i - b == 0 at i = b/a = 0.4
        let r = phi(B / A);
        assert!(r.is_finite());
        assert!((r - 1.0 / D).abs() < 1e-6);
    }

    #[test]
    fn phi_monotone_increasing() {
        let lo = phi(0.5);
        let hi = phi(1.0);
        assert!(hi > lo);
    }

    #[test]
    fn zero_noise_zero_stim_converges_to_fixed_point() {
        // With sigma=0, identical initial conditions, balanced dynamics:
        // s1 and s2 should stay very close and bounded in [0, 1].
        let (tau_s, gamma, j_n, j_cross, i_0, _, dt) = params();
        let n = 10_000;
        let stim = vec![0.0_f64; n];
        let xi = vec![0.0_f64; 2 * n];
        let mut s1o = vec![0.0_f64; n];
        let mut s2o = vec![0.0_f64; n];
        let mut r1o = vec![0.0_f64; n];
        let mut r2o = vec![0.0_f64; n];
        let (s1_f, s2_f) = simulate(
            0.1, 0.1, tau_s, gamma, j_n, j_cross, i_0, 0.0, dt, &stim, &stim, &xi, &mut s1o,
            &mut s2o, &mut r1o, &mut r2o,
        );
        assert!(s1_f.is_finite() && s2_f.is_finite());
        assert!((0.0..=1.0).contains(&s1_f));
        assert!((0.0..=1.0).contains(&s2_f));
        assert!(
            (s1_f - s2_f).abs() < 1e-9,
            "symmetric init must stay symmetric under zero noise"
        );
    }

    #[test]
    fn biased_stimulus_drives_winner() {
        // Strong stim1, no stim2, no noise → s1 > s2 within a few time constants.
        let (tau_s, gamma, j_n, j_cross, i_0, _, dt) = params();
        let n = 50_000;
        let stim1 = vec![0.2_f64; n];
        let stim2 = vec![0.0_f64; n];
        let xi = vec![0.0_f64; 2 * n];
        let mut s1o = vec![0.0_f64; n];
        let mut s2o = vec![0.0_f64; n];
        let mut r1o = vec![0.0_f64; n];
        let mut r2o = vec![0.0_f64; n];
        let (s1_f, s2_f) = simulate(
            0.1, 0.1, tau_s, gamma, j_n, j_cross, i_0, 0.0, dt, &stim1, &stim2, &xi, &mut s1o,
            &mut s2o, &mut r1o, &mut r2o,
        );
        assert!(s1_f > 0.5, "winner s1 should reach attractor; got {s1_f}");
        assert!(s2_f < 0.2, "loser s2 should be suppressed; got {s2_f}");
    }

    #[test]
    fn output_trace_shape_matches_input() {
        let n = 128;
        let stim = vec![0.1_f64; n];
        let xi = vec![0.0_f64; 2 * n];
        let mut s1o = vec![0.0_f64; n];
        let mut s2o = vec![0.0_f64; n];
        let mut r1o = vec![0.0_f64; n];
        let mut r2o = vec![0.0_f64; n];
        simulate(
            0.1, 0.1, 0.1, 0.641, 0.2609, 0.0497, 0.3255, 0.0, 0.001, &stim, &stim, &xi, &mut s1o,
            &mut s2o, &mut r1o, &mut r2o,
        );
        // Every step wrote; nothing still at sentinel zero (for r_out —
        // s_out might genuinely be near 0, so check r_out which is phi(I) ≥ some positive
        // floor when I > 0 on this parameter set).
        assert!(r1o.iter().all(|&r| r > 0.0));
        assert!(r2o.iter().all(|&r| r > 0.0));
    }

    #[test]
    #[should_panic(expected = "xi length must be 2 * n_steps")]
    fn mismatched_xi_length_panics() {
        let n = 10;
        let stim = vec![0.0_f64; n];
        let xi = vec![0.0_f64; n]; // wrong — should be 2*n
        let mut s1o = vec![0.0_f64; n];
        let mut s2o = vec![0.0_f64; n];
        let mut r1o = vec![0.0_f64; n];
        let mut r2o = vec![0.0_f64; n];
        simulate(
            0.1, 0.1, 0.1, 0.641, 0.2609, 0.0497, 0.3255, 0.0, 0.001, &stim, &stim, &xi, &mut s1o,
            &mut s2o, &mut r1o, &mut r2o,
        );
    }
}