Skip to content

Rust Analysis Engine -- 22 Modules, 597 Tests

Pure Rust implementations of all 22 spike_stats analysis modules, exposed to Python via PyO3 bindings. Zero NumPy or SciPy dependency on the Rust side -- all linear algebra (eigendecomposition, matrix inversion, Gaussian elimination), FFT (via rustfft), and stochastic processes (via rand_distr) are implemented natively.

Architecture

Text Only
engine/src/analysis/
  mod.rs              -- 22 pub mod declarations
  basic.rs            -- spike_times, isi, firing_rate, spike_count, bin_spike_train
  rate.rs             -- instantaneous_rate, psth, population_rate
  variability.rs      -- cv_isi, cv2, local_variation, fano_factor, entropies, Hurst
  correlation.rs      -- cross_correlation, STTC, coherence, covariance
  distance.rs         -- van_rossum, victor_purpura, ISI/SPIKE distances
  information.rs      -- mutual_information, transfer_entropy, KL MI
  causality.rs        -- pairwise/conditional/spectral Granger, PDC, DTF
  decoding.rs         -- population_vector, Bayesian, MLE, LDA, naive Bayes
  network.rs          -- functional_connectivity, unitary_events, assemblies
  surrogates.rs       -- ISI shuffle, dither, Poisson, gamma, compound Poisson
  temporal.rs         -- burst_detection, latency, onset, change points
  patterns.rs         -- spike_directionality, spike_train_order, C3
  spectral.rs         -- power_spectrum (FFT via rustfft)
  waveform.rs         -- width, amplitude, slopes, halfwidth, PT ratio
  point_process.rs    -- conditional_intensity, hazard, survivor, renewal
  statistics.rs       -- significance_bootstrap (Rust-only, takes fn pointer)
  stimulus.rs         -- STA, STC, spatial_information, place fields, tuning
  lfp.rs              -- phase_locking_value, spike_field_coherence, phase histogram
  sorting_quality.rs  -- isolation_distance, L-ratio, silhouette, d', SNR, drift
  dimensionality.rs   -- PCA, demixed PCA, factor analysis (Jacobi eigendecomp)
  gpfa.rs             -- GPFA via EM with squared-exponential GP priors
  spade.rs            -- SPADE: Apriori itemset mining + surrogate significance

Dependency Map

External crates used by analysis modules:

Crate Version Used by
rustfft 6.2 spectral.rs, lfp.rs, correlation.rs
rand_distr 0.6 surrogates.rs

All other computation (eigendecomposition, matrix inversion, complex arithmetic, Hilbert transforms, histogram binning) is self-contained.

Custom Linear Algebra

Several modules implement their own numerical primitives to avoid external LAPACK/BLAS dependencies:

  • Jacobi eigendecomposition (dimensionality.rs): symmetric eigenvalues + eigenvectors via iterative rotations, O(n^3) per sweep.
  • Gauss-Jordan inverse (sorting_quality.rs, gpfa.rs, dimensionality.rs, causality.rs): partial-pivot elimination on augmented [A|I] matrix.
  • Gaussian elimination solver (gpfa.rs, causality.rs, decoding.rs): solves Ax = b with partial pivoting.
  • Complex matrix ops (causality.rs): C64 type with inverse, determinant, multiply for spectral Granger.

PyO3 Bindings

96 functions are exposed to Python via #[pyfunction] wrappers in engine/src/lib.rs. The sole exception is significance_bootstrap in statistics.rs which takes a Fn(&[f64], &[f64]) -> f64 pointer and cannot cross the PyO3 boundary.

Functions NOT exposed via PyO3

Function Reason
significance_bootstrap Takes Fn pointer
generalized_victor_purpura Takes fn(f64) -> f64
time_rescaling_ks_test Takes fn(f64) -> f64
inhomogeneous_poisson Takes fn(f64) -> f64
surrogate_trial_shuffle Returns permutation indices only

These are callable from Rust code and tests but require Python-side wrappers with closures.

Module Reference

basic

Core spike train operations used by most other modules.

Function Signature Description
spike_times (&[i32], f64) -> Vec<f64> Extract spike times (s) from binary array
isi (&[i32], f64) -> Vec<f64> Inter-spike intervals (s)
firing_rate (&[i32], f64) -> f64 Mean firing rate (Hz)
spike_count (&[i32]) -> i64 Total spike count
bin_spike_train (&[i32], usize) -> Vec<i64> Bin into spike counts

spectral

Function Signature Description
power_spectrum (&[i32], f64) -> (Vec<f64>, Vec<f64>) PSD via FFT, returns (psd, freqs_hz)

waveform

Spike waveform shape analysis. All functions take &[f64] waveform samples and dt (sample interval, default 1/30000 s).

Function Returns Reference
waveform_width f64 Trough-to-peak width (s). Bartho et al. 2004
waveform_amplitude f64 Peak-to-trough amplitude
waveform_repolarization_slope f64 Max dV/dt after trough. Bean 2007
waveform_recovery_slope f64 Min dV/dt after peak. Bean 2007
waveform_halfwidth f64 Duration at half-minimum. Bartho et al. 2004
waveform_pt_ratio f64 Post-trough peak / trough amplitude

point_process

Point process models and hazard functions.

Function Returns Reference
conditional_intensity Vec<f64> Moving-window Poisson rate (Hz). Brown et al. 2004
isi_hazard_function (Vec<f64>, Vec<f64>) h(t) = f(t)/S(t). Tuckwell 1988
isi_survivor_function (Vec<f64>, Vec<f64>) S(t) = P(ISI > t)
renewal_density (Vec<f64>, Vec<f64>) Normalised by mean rate. Cox 1962

statistics

Function Signature Notes
significance_bootstrap (F, &[f64], &[f64], usize, u64) -> (f64, f64) Rust-only. Permutation test with splitmix64 PRNG

stimulus

Spike-triggered analysis and receptive field estimation.

Function Returns Reference
spike_triggered_average Vec<f64> Mean pre-spike stimulus snippet
spike_triggered_covariance Vec<f64> (flat matrix) Covariance of pre-spike stimulus. Schwartz et al. 2006
spatial_information f64 (bits/spike) Skaggs et al. 1993
place_field_detection Vec<(f64, f64)> Contiguous high-rate bins. O'Keefe & Dostrovsky 1971
tuning_curve (Vec<f64>, Vec<f64>) Rate vs stimulus value. Dayan & Abbott 2001

lfp

Spike-LFP coupling via analytic signal (Hilbert transform with rustfft).

Function Returns Description
phase_locking_value f64 PLV =
spike_field_coherence (Vec<f64>, Vec<f64>) SFC =
spike_phase_histogram (Vec<i64>, Vec<f64>) Phase histogram in [-pi, pi]

sorting_quality

Spike sorting quality metrics. Cluster/noise inputs are row-major (n_points, n_features) flat arrays.

Function Returns Reference
isolation_distance f64 Mahalanobis at rank n_cluster. Harris et al. 2001
l_ratio f64 Normalised chi2-approximation. Schmitzer-Torbert et al. 2005
silhouette_score f64 Mean (b-a)/max(a,b). Rousseeuw 1987
d_prime f64 Projected sensitivity index. Green & Swets 1966
isi_violation_rate f64 Fraction below refractory. Hill et al. 2011
presence_ratio f64 Fraction of occupied bins. IBL 2019
amplitude_cutoff f64 Estimated missing spikes. Hill et al. 2011
snr f64 Peak / noise_std. Suner et al. 2005
nn_hit_rate f64 k-NN cluster purity. Chung et al. 2017
drift_metric f64 Amplitude drift over time. IBL 2019

dimensionality

Dimensionality reduction with built-in Jacobi eigendecomposition (no LAPACK required).

Function Returns Reference
spike_train_pca (Vec<f64>, Vec<f64>) (projected, explained_variance_ratio)
demixed_pca (Vec<f64>, Vec<f64>) Condition-dependent variance. Kobak et al. 2016
factor_analysis (Vec<f64>, Vec<f64>) (loadings, uniquenesses). Rubin & Thayer 1982

gpfa

Gaussian Process Factor Analysis. Full EM implementation with squared-exponential GP kernels, block-structured precision matrices, and approximate log-likelihood monitoring.

Function Returns Reference
gpfa GpfaResult struct Trajectories, C, d, R, tau, log-likelihoods. Yu et al. 2009
gpfa_transform Vec<f64> Project new data with learned parameters

GpfaResult fields: trajectories, c, d, r, tau, log_likelihoods, n_latents, n_bins, n_neurons.

spade

Spike Pattern Detection and Evaluation. Apriori-style frequent itemset mining extended to spatiotemporal patterns with lag search and surrogate significance testing.

Function Returns Reference
spade_detect Vec<SpadePattern> Significant patterns with p-values. Torre et al. 2013

SpadePattern fields: neurons, lags, count, p_value.

Test Coverage

597 tests across all modules. Every public function has at least:

  • Positive case (typical input producing expected output)
  • Edge case (empty input, single element, degenerate dimensions)
  • Numerical accuracy check (comparison against known values or bounds)

Run tests:

Bash
export PATH="$HOME/.rustup/toolchains/stable-x86_64-unknown-linux-gnu/bin:$PATH"
cd engine && cargo test --lib

Run a single module's tests:

Bash
cargo test --lib analysis::spectral
cargo test --lib analysis::gpfa

Benchmark Results

Measured with Criterion on i5-11600K @ 3.90 GHz, DDR4-3200. Values are median latency. Pure CPU benchmarks.

Last measured: 2026-05-02 (cargo bench --bench analysis_bench)

basic

Function 100 10K 100K
spike_times 86 ns 5.4 µs 85.0 µs
isi 106 ns 5.6 µs 86.3 µs
firing_rate 13 ns 1.2 µs 15.3 µs
bin_spike_train 36 ns 2.5 µs 29.0 µs

rate

Function 100 10K 100K
instantaneous_rate 4.3 µs 352 µs 3.5 ms

variability

Function 100 10K 100K
cv_isi 108 ns 6.5 µs 94.7 µs
fano_factor 76 ns 4.8 µs 50.7 µs
sample_entropy 32.6 µs 3.23 ms O(n²) — not benchmarked

correlation

Function 100 10K
cross_correlation 1.46 µs 213 µs
event_synchronization 142 ns 63.1 µs

distance

Function 100 5K
van_rossum 1.07 µs 314 µs
victor_purpura 162 ns 274 µs
isi_distance 165 ns 5.95 µs

information

Function 100 10K
mutual_information 714 ns 48.7 µs
transfer_entropy 1.03 µs 60.8 µs

causality

Function 100 5K
pairwise_granger 102 ns 81.7 µs

decoding

Function Latency
population_vector_decode 14.9 µs
bayesian_decode 11.0 µs

network

Function Latency
functional_connectivity (10n) 814 µs

surrogates

Function 1K 100K
isi_shuffle 865 ns 103 µs
homogeneous_poisson 2.55 µs 242 µs

temporal

Function 1K 100K
burst_detection 708 ns 91.4 µs
change_point_detection 296 ns 31.8 µs

patterns

Function Latency
spike_directionality (5K) 3.85 µs
cubic_higher_order (5K) 409 µs

spectral

Function 256 10K 100K
power_spectrum 3.81 µs 177 µs 4.91 ms

waveform (50 samples)

Function Latency
waveform_width 64 ns
waveform_amplitude 38 ns
waveform_repolarization_slope 72 ns
waveform_halfwidth 168 ns
waveform_pt_ratio 58 ns

point_process

Function 1K 100K
conditional_intensity 22.9 µs 2.33 ms
isi_hazard 1.07 µs 103 µs

statistics

Function Latency
significance_bootstrap (200 surr) 661 µs

stimulus

Function 1K 50K
STA 1.18 µs 67.8 µs
spatial_information 3.89 µs 205 µs

lfp

Function 500 10K
phase_locking_value 26.2 µs 540 µs
spike_field_coherence 11.9 µs 237 µs

sorting_quality

Function Latency
isolation_distance (50 pts, 4D) 6.82 µs
isolation_distance (200 pts, 4D) 27.0 µs
silhouette_score (50 pts, 4D) 35.8 µs
silhouette_score (200 pts, 4D) 536 µs
isi_violation_rate (5K) 2.93 µs

dimensionality

Function Latency
PCA (10n, 2000t) 43.4 µs
factor_analysis (10n, 2000t) 408 µs

gpfa

Function Latency
GPFA (4n, 500t, 5 EM iter) 4.84 ms

spade

Function Latency
SPADE detect (3n, 500t, 50 surr) 748 µs

Performance Notes

All functions are single-threaded pure CPU. Benchmarks run via cargo bench --bench analysis_bench. Parallelisation is left to the Python caller via concurrent.futures or joblib.

Full benchmark results with scaling analysis: see Analysis Benchmarks.