oldies_core/
lib.rs

1//! # OldiesRules Core
2//!
3//! Shared types and utilities for legacy simulator revival.
4//!
5//! ## Supported Legacy Simulators
6//!
7//! | Simulator | Original Language | Era | Status |
8//! |-----------|------------------|-----|--------|
9//! | GENESIS | SLI + C | 1980s-2014 | Semi-abandoned |
10//! | XPPAUT | C + FORTRAN | 1990s | Hobby project |
11//! | AUTO | FORTRAN | 1980s | Legacy |
12//! | ModelDB | Various | 1996+ | Active but legacy |
13//!
14//! ## Design Philosophy
15//!
16//! 1. Preserve numerical equivalence with originals
17//! 2. Modern Rust safety and performance
18//! 3. HumanBrain integration ready
19
20use ndarray::{Array1, Array2};
21use serde::{Deserialize, Serialize};
22use thiserror::Error;
23
24/// Simulator type
25#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
26pub enum Simulator {
27    Genesis,
28    Xppaut,
29    Auto,
30    ModelDB,
31    Neuron,
32    Brian,
33}
34
35/// Common errors
36#[derive(Debug, Error)]
37pub enum OldiesError {
38    #[error("Parse error: {0}")]
39    ParseError(String),
40
41    #[error("Simulation error: {0}")]
42    SimulationError(String),
43
44    #[error("Model not found: {0}")]
45    ModelNotFound(String),
46
47    #[error("IO error: {0}")]
48    IoError(#[from] std::io::Error),
49
50    #[error("Numerical error: {0}")]
51    NumericalError(String),
52}
53
54pub type Result<T> = std::result::Result<T, OldiesError>;
55
56/// Time point
57pub type Time = f64;
58
59/// Voltage (mV)
60pub type Voltage = f64;
61
62/// Current (nA)
63pub type Current = f64;
64
65/// Conductance (mS/cm^2)
66pub type Conductance = f64;
67
68/// Concentration (mM)
69pub type Concentration = f64;
70
71/// State vector for ODE systems
72pub type StateVector = Array1<f64>;
73
74/// Time series data
75#[derive(Debug, Clone, Serialize, Deserialize)]
76pub struct TimeSeries {
77    /// Time points
78    pub time: Vec<Time>,
79    /// Values at each time point
80    pub values: Vec<f64>,
81    /// Variable name
82    pub name: String,
83    /// Units
84    pub units: Option<String>,
85}
86
87impl TimeSeries {
88    pub fn new(name: &str) -> Self {
89        Self {
90            time: Vec::new(),
91            values: Vec::new(),
92            name: name.to_string(),
93            units: None,
94        }
95    }
96
97    pub fn push(&mut self, t: Time, v: f64) {
98        self.time.push(t);
99        self.values.push(v);
100    }
101
102    pub fn len(&self) -> usize {
103        self.time.len()
104    }
105
106    pub fn is_empty(&self) -> bool {
107        self.time.is_empty()
108    }
109}
110
111/// ODE system trait (for simulators)
112pub trait OdeSystem {
113    /// System dimension
114    fn dimension(&self) -> usize;
115
116    /// Compute derivatives: dy/dt = f(t, y)
117    fn derivatives(&self, t: Time, y: &StateVector) -> StateVector;
118
119    /// Optional Jacobian for stiff systems
120    fn jacobian(&self, _t: Time, _y: &StateVector) -> Option<Array2<f64>> {
121        None
122    }
123}
124
125/// Simulation parameters
126#[derive(Debug, Clone, Serialize, Deserialize)]
127pub struct SimulationParams {
128    /// Start time
129    pub t_start: Time,
130    /// End time
131    pub t_end: Time,
132    /// Time step
133    pub dt: Time,
134    /// Output interval (for recording)
135    pub output_dt: Option<Time>,
136    /// Solver tolerance
137    pub tolerance: f64,
138}
139
140impl Default for SimulationParams {
141    fn default() -> Self {
142        Self {
143            t_start: 0.0,
144            t_end: 100.0,
145            dt: 0.01,
146            output_dt: Some(0.1),
147            tolerance: 1e-6,
148        }
149    }
150}
151
152/// Ion channel model
153#[derive(Debug, Clone, Serialize, Deserialize)]
154pub struct IonChannel {
155    /// Channel name
156    pub name: String,
157    /// Maximum conductance (mS/cm^2)
158    pub g_max: Conductance,
159    /// Reversal potential (mV)
160    pub e_rev: Voltage,
161    /// Gate variables
162    pub gates: Vec<GateVariable>,
163}
164
165/// Gate variable (e.g., m, h, n in Hodgkin-Huxley)
166#[derive(Debug, Clone, Serialize, Deserialize)]
167pub struct GateVariable {
168    /// Variable name
169    pub name: String,
170    /// Power (exponent in gating)
171    pub power: u32,
172    /// Alpha rate function parameters
173    pub alpha: RateFunction,
174    /// Beta rate function parameters
175    pub beta: RateFunction,
176}
177
178/// Rate function type
179#[derive(Debug, Clone, Serialize, Deserialize)]
180pub enum RateFunction {
181    /// Standard HH form: A*(V+B)/(exp((V+B)/C)-1)
182    HodgkinHuxley { a: f64, b: f64, c: f64 },
183    /// Exponential: A*exp((V+B)/C)
184    Exponential { a: f64, b: f64, c: f64 },
185    /// Sigmoid: A/(1+exp((V+B)/C))
186    Sigmoid { a: f64, b: f64, c: f64 },
187    /// Linear: A*(V+B)
188    Linear { a: f64, b: f64 },
189    /// Constant
190    Constant(f64),
191}
192
193impl RateFunction {
194    /// Evaluate rate at given voltage
195    pub fn eval(&self, v: Voltage) -> f64 {
196        match self {
197            Self::HodgkinHuxley { a, b, c } => {
198                let x = (v + b) / c;
199                if x.abs() < 1e-6 {
200                    // L'Hopital's rule for x -> 0
201                    a * c
202                } else {
203                    a * (v + b) / (x.exp() - 1.0)
204                }
205            }
206            Self::Exponential { a, b, c } => {
207                a * ((v + b) / c).exp()
208            }
209            Self::Sigmoid { a, b, c } => {
210                a / (1.0 + ((v + b) / c).exp())
211            }
212            Self::Linear { a, b } => {
213                a * (v + b)
214            }
215            Self::Constant(c) => *c,
216        }
217    }
218}
219
220#[cfg(test)]
221mod tests {
222    use super::*;
223
224    #[test]
225    fn test_rate_functions() {
226        let hh = RateFunction::HodgkinHuxley { a: 0.1, b: 40.0, c: 10.0 };
227        let rate = hh.eval(-65.0);
228        assert!(rate > 0.0);
229
230        let exp = RateFunction::Exponential { a: 0.1, b: 65.0, c: 80.0 };
231        let rate = exp.eval(-65.0);
232        assert!(rate > 0.0);
233    }
234
235    #[test]
236    fn test_time_series() {
237        let mut ts = TimeSeries::new("voltage");
238        ts.push(0.0, -65.0);
239        ts.push(0.1, -64.0);
240        assert_eq!(ts.len(), 2);
241    }
242}
oldies_core/lib.rs

oldies_core/
lib.rs
