modewise_fc/
lib.rs

1//! Functional connectivity: Pearson correlation and phase locking value, for a whole run or
2//! per epoch (docs/connectivity.md).
3
4mod measures;
5
6use modewise_connectivity::{FunctionalConnectivity, Measure, WindowInfo};
7use modewise_decomposition::Decomposition;
8use modewise_epochs::{EpochSource, Epochs};
9use modewise_hilbert::Padding;
10use modewise_method::MethodError;
11use modewise_signals::MultichannelSeries;
12use ndarray::{Array4, ArrayView2};
13use rayon::prelude::*;
14
15pub use measures::{correlation, phase_locking, unit_phases};
16
17/// What to compute functional connectivity on.
18#[derive(Debug, Clone, Copy)]
19pub enum Source<'a> {
20    Series(&'a MultichannelSeries),
21    /// One matrix per mode.
22    Decomposition(&'a Decomposition),
23    /// One matrix per epoch (and per mode for a decomposition source). Phases come from the
24    /// full run before cutting, so epoch edges add no Hilbert edge effects.
25    Epochs(&'a Epochs),
26}
27
28/// Functional connectivity of every component, laid out as window × component × channel × channel.
29pub fn functional_connectivity(
30    source: Source<'_>,
31    measure: Measure,
32) -> Result<FunctionalConnectivity, MethodError> {
33    let (components, labels, centres, epochs) = match source {
34        Source::Series(s) => (vec![s.data()], s.labels().to_vec(), None, None),
35        Source::Decomposition(d) => (
36            modes(d),
37            d.labels().to_vec(),
38            Some(d.centre_frequencies_hz().to_vec()),
39            None,
40        ),
41        Source::Epochs(e) => match e.source() {
42            EpochSource::Series(s) => (vec![s.data()], s.labels().to_vec(), None, Some(e)),
43            EpochSource::Decomposition(d) => (
44                modes(d),
45                d.labels().to_vec(),
46                Some(d.centre_frequencies_hz().to_vec()),
47                Some(e),
48            ),
49            _ => {
50                return Err(MethodError::invalid(
51                    "data",
52                    "functional connectivity needs epochs of a series or a decomposition",
53                ));
54            }
55        },
56    };
57    let n_samples = components.first().map_or(0, |c| c.ncols());
58    let windows: Vec<(usize, usize)> = match epochs {
59        Some(e) => e.windows().iter().map(|w| (w.start, w.end)).collect(),
60        None => vec![(0, n_samples)],
61    };
62
63    let n_channels = labels.len();
64    let mut values = Array4::zeros((windows.len(), components.len(), n_channels, n_channels));
65    let phases: Option<Vec<_>> = match measure {
66        Measure::PhaseLockingValue => Some(
67            components
68                .par_iter()
69                .map(|c| unit_phases(*c, Padding::Mirror))
70                .collect(),
71        ),
72        Measure::Pearson { .. } => None,
73    };
74    let cells: Vec<(usize, usize)> = (0..windows.len())
75        .flat_map(|w| (0..components.len()).map(move |k| (w, k)))
76        .collect();
77    let matrices: Vec<_> = cells
78        .par_iter()
79        .map(|&(w, k)| {
80            let (start, end) = windows[w];
81            let range = ndarray::s![.., start..end];
82            match (&phases, measure) {
83                (Some(p), _) => phase_locking(&p[k].0.slice(range), &p[k].1.slice(range), &p[k].2),
84                (None, Measure::Pearson { fisher_z }) => {
85                    let r = correlation(components[k].slice(range));
86                    if fisher_z { measures::fisher_z(r) } else { r }
87                }
88                (None, Measure::PhaseLockingValue) => {
89                    unreachable!("phases exist for the phase locking value")
90                }
91            }
92        })
93        .collect();
94    for (&(w, k), matrix) in cells.iter().zip(matrices) {
95        values.slice_mut(ndarray::s![w, k, .., ..]).assign(&matrix);
96    }
97
98    let mut result = FunctionalConnectivity::new(values, measure, labels)?;
99    if let Some(hz) = centres {
100        result = result.with_centre_frequencies(hz)?;
101    }
102    if let Some(e) = epochs {
103        let info = e
104            .windows()
105            .iter()
106            .map(|w| WindowInfo {
107                label: w.label.clone(),
108                event_onset_s: w.event_onset_s,
109            })
110            .collect();
111        result = result.with_windows(info, e.resolvable())?;
112    }
113    Ok(result)
114}
115
116fn modes(d: &Decomposition) -> Vec<ArrayView2<'_, f64>> {
117    (0..d.n_modes()).map(|k| d.mode(k)).collect()
118}
modewise_fc/lib.rs

modewise_fc/
lib.rs
