sc_neurocore_engine/neurons/
ermentrout_kopell_pop.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 — Montbrió-Pazó-Roxin exact QIF mean-field engine
8
9//! The source prints dimensionless equations (12a-b) in `(R, v, t')`.
10//! This engine restores physical rate and time through `R = tau * r` and
11//! `t' = t / tau`, then applies atomic simultaneous-Euler batches.
12
13use std::error::Error;
14use std::fmt;
15
16/// Typed numerical and caller-contract failures for the MPR population.
17#[derive(Clone, Copy, Debug, Eq, PartialEq)]
18pub enum ErmentroutKopellPopulationError {
19    /// One state or configuration value is not finite.
20    NonFiniteConfiguration,
21    /// The initial population firing rate is negative.
22    NegativeInitialRate,
23    /// The time scale or explicit-Euler step is not positive.
24    NonPositiveTimeScale,
25    /// The Lorentzian half-width is negative.
26    NegativeHalfWidth,
27    /// One external-drive value is not finite.
28    NonFiniteInput,
29    /// A simultaneous Euler candidate contains a non-finite state.
30    NonFiniteCandidate,
31    /// A simultaneous Euler candidate has a negative firing rate.
32    NegativeCandidateRate,
33}
34
35impl fmt::Display for ErmentroutKopellPopulationError {
36    fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
37        let message = match self {
38            Self::NonFiniteConfiguration => "MPR state and parameters must be finite",
39            Self::NegativeInitialRate => "MPR firing rate must be non-negative",
40            Self::NonPositiveTimeScale => "MPR tau and dt must be positive",
41            Self::NegativeHalfWidth => "MPR Lorentzian half-width must be non-negative",
42            Self::NonFiniteInput => "MPR external input must contain only finite values",
43            Self::NonFiniteCandidate => "MPR candidate state must remain finite",
44            Self::NegativeCandidateRate => "MPR candidate firing rate became negative",
45        };
46        formatter.write_str(message)
47    }
48}
49
50impl Error for ErmentroutKopellPopulationError {}
51
52/// Legacy-named public wrapper for the Montbrió-Pazó-Roxin population model.
53#[derive(Clone, Debug)]
54pub struct ErmentroutKopellPopulation {
55    /// Population firing rate.
56    pub r: f64,
57    /// Mean membrane potential.
58    pub v: f64,
59    /// Positive membrane time scale.
60    pub tau: f64,
61    /// Non-negative Lorentzian half-width of neuronal excitability.
62    pub delta: f64,
63    /// Centre of the neuronal excitability distribution.
64    pub eta_bar: f64,
65    /// Recurrent coupling strength.
66    pub j: f64,
67    /// Positive explicit-Euler step.
68    pub dt: f64,
69}
70
71impl ErmentroutKopellPopulation {
72    /// Construct the phase-portrait parameter set used by the source paper.
73    pub fn new() -> Self {
74        Self {
75            r: 0.1,
76            v: -2.0,
77            tau: 1.0,
78            delta: 1.0,
79            eta_bar: -5.0,
80            j: 15.0,
81            dt: 0.01,
82        }
83    }
84
85    /// Construct and validate one complete numerical configuration.
86    pub fn with_parameters(
87        r: f64,
88        v: f64,
89        tau: f64,
90        delta: f64,
91        eta_bar: f64,
92        j: f64,
93        dt: f64,
94    ) -> Result<Self, ErmentroutKopellPopulationError> {
95        let population = Self {
96            r,
97            v,
98            tau,
99            delta,
100            eta_bar,
101            j,
102            dt,
103        };
104        population.validate()?;
105        Ok(population)
106    }
107
108    fn validate(&self) -> Result<(), ErmentroutKopellPopulationError> {
109        if ![
110            self.r,
111            self.v,
112            self.tau,
113            self.delta,
114            self.eta_bar,
115            self.j,
116            self.dt,
117        ]
118        .into_iter()
119        .all(f64::is_finite)
120        {
121            return Err(ErmentroutKopellPopulationError::NonFiniteConfiguration);
122        }
123        if self.r < 0.0 {
124            return Err(ErmentroutKopellPopulationError::NegativeInitialRate);
125        }
126        if self.tau <= 0.0 || self.dt <= 0.0 {
127            return Err(ErmentroutKopellPopulationError::NonPositiveTimeScale);
128        }
129        if self.delta < 0.0 {
130            return Err(ErmentroutKopellPopulationError::NegativeHalfWidth);
131        }
132        Ok(())
133    }
134
135    #[inline]
136    fn derivatives(&self, drive: f64) -> (f64, f64) {
137        let scaled_rate = std::f64::consts::PI * self.tau * self.r;
138        let dr = self.delta / (std::f64::consts::PI * self.tau * self.tau)
139            + 2.0 * self.r * self.v / self.tau;
140        let dv = (self.v * self.v + self.eta_bar + drive + self.j * self.tau * self.r
141            - scaled_rate * scaled_rate)
142            / self.tau;
143        (dr, dv)
144    }
145
146    /// Advance one simultaneous explicit-Euler step with explicit error reporting.
147    pub fn try_step(&mut self, ext_input: f64) -> Result<f64, ErmentroutKopellPopulationError> {
148        self.validate()?;
149        if !ext_input.is_finite() {
150            return Err(ErmentroutKopellPopulationError::NonFiniteInput);
151        }
152        let (dr, dv) = self.derivatives(ext_input);
153        let next_r = self.r + self.dt * dr;
154        let next_v = self.v + self.dt * dv;
155        if !next_r.is_finite() || !next_v.is_finite() {
156            return Err(ErmentroutKopellPopulationError::NonFiniteCandidate);
157        }
158        if next_r < 0.0 {
159            return Err(ErmentroutKopellPopulationError::NegativeCandidateRate);
160        }
161        self.r = next_r;
162        self.v = next_v;
163        Ok(self.r)
164    }
165
166    /// Advance one step while preserving the legacy scalar-returning Rust API.
167    ///
168    /// Invalid input is fail-closed: the state remains unchanged and the
169    /// current firing rate is returned. New callers that need diagnostics
170    /// should use [`Self::try_step`].
171    pub fn step(&mut self, ext_input: f64) -> f64 {
172        match self.try_step(ext_input) {
173            Ok(rate) => rate,
174            Err(_) => self.r,
175        }
176    }
177
178    /// Restore dynamic states while preserving the configured parameters.
179    pub fn reset(&mut self) {
180        self.r = 0.1;
181        self.v = -2.0;
182    }
183}
184
185impl Default for ErmentroutKopellPopulation {
186    fn default() -> Self {
187        Self::new()
188    }
189}
190
191/// Per-step states and final receipt from one complete drive batch.
192pub struct ErmentroutKopellPopulationTrace {
193    /// Post-update firing-rate trace.
194    pub r: Vec<f64>,
195    /// Post-update mean-voltage trace.
196    pub v: Vec<f64>,
197    /// Final ``[r, v]`` receipt, including the initial state for an empty batch.
198    pub final_state: [f64; 2],
199}
200
201/// Simulate a complete caller-owned external-drive vector atomically.
202#[expect(
203    clippy::too_many_arguments,
204    reason = "native parity surface carries the complete scientific configuration"
205)]
206pub fn simulate(
207    r: f64,
208    v: f64,
209    tau: f64,
210    delta: f64,
211    eta_bar: f64,
212    j: f64,
213    dt: f64,
214    ext_input: &[f64],
215) -> Result<ErmentroutKopellPopulationTrace, ErmentroutKopellPopulationError> {
216    let mut population =
217        ErmentroutKopellPopulation::with_parameters(r, v, tau, delta, eta_bar, j, dt)?;
218    if !ext_input.iter().all(|value| value.is_finite()) {
219        return Err(ErmentroutKopellPopulationError::NonFiniteInput);
220    }
221    let mut r_trace = Vec::with_capacity(ext_input.len());
222    let mut v_trace = Vec::with_capacity(ext_input.len());
223    for drive in ext_input {
224        population.try_step(*drive)?;
225        r_trace.push(population.r);
226        v_trace.push(population.v);
227    }
228    Ok(ErmentroutKopellPopulationTrace {
229        r: r_trace,
230        v: v_trace,
231        final_state: [population.r, population.v],
232    })
233}
234
235#[cfg(test)]
236mod tests {
237    use super::*;
238
239    #[test]
240    fn one_step_matches_equation_twelve_with_explicit_tau() {
241        let mut population =
242            ErmentroutKopellPopulation::with_parameters(0.2, -1.5, 2.0, 0.7, -3.0, 12.0, 0.005)
243                .unwrap();
244        let old_r = population.r;
245        let old_v = population.v;
246        let drive = 1.25;
247        let expected_r =
248            old_r + 0.005 * (0.7 / (std::f64::consts::PI * 4.0) + 2.0 * old_r * old_v / 2.0);
249        let expected_v = old_v
250            + 0.005
251                * (old_v * old_v + -3.0 + drive + 12.0 * 2.0 * old_r
252                    - (std::f64::consts::PI * 2.0 * old_r).powi(2))
253                / 2.0;
254        population.try_step(drive).unwrap();
255        assert_eq!(population.r, expected_r);
256        assert_eq!(population.v, expected_v);
257    }
258
259    #[test]
260    fn invalid_step_is_atomic() {
261        let mut population = ErmentroutKopellPopulation::new();
262        let before = (population.r, population.v);
263        assert_eq!(
264            population.try_step(f64::NAN),
265            Err(ErmentroutKopellPopulationError::NonFiniteInput)
266        );
267        assert_eq!((population.r, population.v), before);
268    }
269
270    #[test]
271    fn candidate_failures_are_typed_and_atomic() {
272        let mut negative =
273            ErmentroutKopellPopulation::with_parameters(1.0, -100.0, 1.0, 0.0, 0.0, 0.0, 0.1)
274                .unwrap();
275        let before = (negative.r, negative.v);
276        assert_eq!(
277            negative.try_step(0.0),
278            Err(ErmentroutKopellPopulationError::NegativeCandidateRate)
279        );
280        assert_eq!((negative.r, negative.v), before);
281
282        let mut nonfinite = ErmentroutKopellPopulation::with_parameters(
283            f64::MAX,
284            f64::MAX,
285            1.0,
286            0.0,
287            0.0,
288            0.0,
289            1.0,
290        )
291        .unwrap();
292        let before = (nonfinite.r, nonfinite.v);
293        assert_eq!(
294            nonfinite.try_step(0.0),
295            Err(ErmentroutKopellPopulationError::NonFiniteCandidate)
296        );
297        assert_eq!((nonfinite.r, nonfinite.v), before);
298    }
299
300    #[test]
301    fn batch_matches_scalar_and_empty_preserves_initial_state() {
302        let empty = simulate(0.2, -1.5, 2.0, 0.7, -3.0, 12.0, 0.005, &[]).unwrap();
303        assert!(empty.r.is_empty() && empty.v.is_empty());
304        assert_eq!(empty.final_state, [0.2, -1.5]);
305
306        let drive = [0.0, 0.25, -0.1, 1.0];
307        let batch = simulate(0.1, -2.0, 1.0, 1.0, -5.0, 15.0, 0.01, &drive).unwrap();
308        let mut scalar = ErmentroutKopellPopulation::new();
309        for value in drive {
310            scalar.try_step(value).unwrap();
311        }
312        assert_eq!(batch.final_state, [scalar.r, scalar.v]);
313    }
314}
sc_neurocore_engine/neurons/ermentrout_kopell_pop.rs

sc_neurocore_engine/neurons/
ermentrout_kopell_pop.rs
