terminals_core/phase/

kuramoto.rs

use std::f32::consts::PI;

/// Kuramoto order parameter R·e^{iψ} = (1/N) Σ e^{iθ_j}
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct OrderParameter {
    /// Synchronization magnitude in [0, 1]
    pub r: f32,
    /// Collective phase in (-π, π]
    pub psi: f32,
}

/// Natural frequencies for oscillators: uniform or per-node.
pub enum Omega {
    Uniform(f32),
    PerNode(Vec<f32>),
}

/// Wrap θ into [0, 2π).
#[inline]
pub fn wrap_to_2pi(theta: f32) -> f32 {
    let two_pi = 2.0 * PI;
    let mut t = theta % two_pi;
    if t < 0.0 {
        t += two_pi;
    }
    t
}

/// Wrap θ into (-π, π].
#[inline]
pub fn wrap_to_pi(theta: f32) -> f32 {
    let two_pi = 2.0 * PI;
    let mut t = (theta + PI) % two_pi;
    if t < 0.0 {
        t += two_pi;
    }
    t - PI
}

/// Compute the Kuramoto order parameter from a slice of phases.
///
/// R e^{iψ} = (1/N) Σ e^{iθ_j}
pub fn order_parameter(phases: &[f32]) -> OrderParameter {
    let n = phases.len();
    if n == 0 {
        return OrderParameter { r: 0.0, psi: 0.0 };
    }
    let mut c = 0.0f32;
    let mut s = 0.0f32;
    for &theta in phases {
        c += theta.cos();
        s += theta.sin();
    }
    let inv_n = 1.0 / n as f32;
    c *= inv_n;
    s *= inv_n;
    OrderParameter {
        r: c.hypot(s),
        psi: s.atan2(c),
    }
}

/// Euler step of the Kuramoto model.
///
/// dθ_i/dt = ω_i + (K/degree_i) Σ_j A_ij sin(θ_j - θ_i)
///
/// - `phases`: current oscillator phases (radians)
/// - `omega`: natural frequencies (uniform scalar or per-node)
/// - `k`: global coupling strength
/// - `adjacency`: optional weighted adjacency; `None` = fully connected (weight 1)
/// - `dt`: time step
/// - `wrap`: whether to wrap result into (−π, π]
pub fn kuramoto_step(
    phases: &[f32],
    omega: &Omega,
    k: f32,
    adjacency: Option<&[&[f32]]>,
    dt: f32,
    wrap: bool,
) -> Vec<f32> {
    let n = phases.len();
    let mut next = Vec::with_capacity(n);

    for i in 0..n {
        let theta_i = phases[i];
        let w_i = match omega {
            Omega::Uniform(w) => *w,
            Omega::PerNode(ref ws) => ws[i],
        };

        let (coupling_sum, degree) = if let Some(adj) = adjacency {
            let row = adj[i];
            let mut cs = 0.0f32;
            let mut deg = 0.0f32;
            for j in 0..n {
                if j == i {
                    continue;
                }
                let aij = if j < row.len() { row[j] } else { 0.0 };
                if aij == 0.0 {
                    continue;
                }
                cs += aij * (phases[j] - theta_i).sin();
                deg += aij.abs();
            }
            (cs, deg)
        } else {
            let mut cs = 0.0f32;
            for (j, &phase_j) in phases.iter().enumerate().take(n) {
                if j == i {
                    continue;
                }
                cs += (phase_j - theta_i).sin();
            }
            let deg = (n as f32 - 1.0).max(1.0);
            (cs, deg)
        };

        let coupling = if degree > 0.0 {
            (k / degree) * coupling_sum
        } else {
            0.0
        };

        let mut val = theta_i + (w_i + coupling) * dt;
        if wrap {
            val = wrap_to_pi(val);
        }
        next.push(val);
    }

    next
}

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

    #[test]
    fn test_order_parameter_synchronized() {
        let phases = vec![0.0f32; 10];
        let op = order_parameter(&phases);
        assert!((op.r - 1.0).abs() < 1e-6, "R = {}", op.r);
    }

    #[test]
    fn test_order_parameter_uniform_spread() {
        let n = 100usize;
        let phases: Vec<f32> = (0..n).map(|i| 2.0 * PI * i as f32 / n as f32).collect();
        let op = order_parameter(&phases);
        assert!(op.r < 0.1, "R = {} (expected < 0.1)", op.r);
    }

    #[test]
    fn test_order_parameter_empty() {
        let op = order_parameter(&[]);
        assert_eq!(op.r, 0.0);
        assert_eq!(op.psi, 0.0);
    }

    #[test]
    fn test_wrap_to_pi_zero() {
        assert!((wrap_to_pi(0.0) - 0.0).abs() < 1e-6);
    }

    #[test]
    fn test_wrap_to_pi_4pi() {
        assert!(
            (wrap_to_pi(4.0 * PI) - 0.0).abs() < 1e-4,
            "got {}",
            wrap_to_pi(4.0 * PI)
        );
    }

    #[test]
    fn test_wrap_to_pi_3pi() {
        let w = wrap_to_pi(3.0 * PI);
        // 3π ≡ π (mod 2π), so wrap_to_pi should give ±π
        assert!(
            (w - (-PI)).abs() < 1e-4 || (w - PI).abs() < 1e-4,
            "wrap_to_pi(3π) = {}",
            w
        );
    }

    #[test]
    fn test_wrap_to_2pi_negative() {
        let w = wrap_to_2pi(-0.5);
        assert!(w >= 0.0);
        assert!(w < 2.0 * PI);
    }

    #[test]
    fn test_kuramoto_step_converges() {
        let phases = vec![0.0f32, 1.0, 2.0, 3.0, 4.0, 5.0];
        let r_before = order_parameter(&phases).r;
        let next = kuramoto_step(&phases, &Omega::Uniform(0.0), 5.0, None, 0.01, true);
        let r_after = order_parameter(&next).r;
        assert!(
            r_after >= r_before - 0.01,
            "r_before={} r_after={}",
            r_before,
            r_after
        );
    }

    #[test]
    fn test_kuramoto_step_preserves_length() {
        let phases = vec![0.1f32, 0.2, 0.3];
        let next = kuramoto_step(&phases, &Omega::Uniform(1.0), 1.0, None, 0.01, true);
        assert_eq!(next.len(), 3);
    }

    #[test]
    fn test_kuramoto_step_per_node_omega() {
        let phases = vec![0.0f32, 0.0, 0.0];
        let omega = Omega::PerNode(vec![1.0, 2.0, 3.0]);
        let next = kuramoto_step(&phases, &omega, 0.0, None, 0.1, false);
        assert!((next[0] - 0.1).abs() < 1e-6);
        assert!((next[1] - 0.2).abs() < 1e-6);
        assert!((next[2] - 0.3).abs() < 1e-6);
    }

    #[test]
    fn test_kuramoto_step_strong_coupling_increases_r() {
        // After many steps with strong coupling and zero natural frequencies,
        // R should approach 1.
        let mut phases = vec![0.0f32, 1.0, 2.0, 3.0, 4.0, 5.0];
        for _ in 0..1000 {
            phases = kuramoto_step(&phases, &Omega::Uniform(0.0), 5.0, None, 0.01, true);
        }
        let r = order_parameter(&phases).r;
        assert!(r > 0.99, "R = {} after convergence", r);
    }

    #[test]
    fn test_kuramoto_step_no_wrap() {
        let phases = vec![3.0f32];
        // Single oscillator, no coupling, omega = 1, dt = 1 → next = 4.0 (unwrapped)
        let next = kuramoto_step(&phases, &Omega::Uniform(1.0), 0.0, None, 1.0, false);
        assert!((next[0] - 4.0).abs() < 1e-6);
    }
}