# Cookbook

Short, self-contained recipes for the most common tasks. Every code block on
this page is executed as a doctest in the test suite
(`tests/test_cookbook.py`), so the recipes are guaranteed to run against the
current release.

All recipes share this setup. `time_series` has shape
`(n_time_samples, n_trials, n_signals)`; a 2-D `(n_time_samples, n_signals)`
array works too.

```python
>>> import numpy as np
>>> from spectral_connectivity import (
...     multitaper_connectivity,
...     fourier_connectivity,
...     list_measures,
... )
>>> rng = np.random.default_rng(0)
>>> time_series = rng.standard_normal((1000, 4, 3))

```

## Discover the available measures

`list_measures()` enumerates every valid `method` name, with its output
category and a one-line description. Use it instead of guessing method strings.

```python
>>> len(list_measures())
35
>>> [m.name for m in list_measures(default_only=True)][:3]
['coherence_magnitude', 'coherence_phase', 'debiased_squared_phase_lag_index']
>>> [m.name for m in list_measures(directed=True)][:2]
['pairwise_spectral_granger_prediction', 'directed_phase_lag_index']
>>> next(m for m in list_measures() if m.name == "phase_slope_index").requires_two_sided
False
>>> power = next(m for m in list_measures() if m.name == "power")
>>> power.category, power.description
('power', 'Return the one-sided power spectral density of the signal.')

```

Passing an unknown name raises a helpful error rather than an obscure
`AttributeError`:

```python
>>> multitaper_connectivity(
...     time_series, sampling_frequency=500, method="coherence"
... )
Traceback (most recent call last):
    ...
ValueError: 'coherence' is not a known connectivity measure. Did you mean: 'coherence_magnitude', 'coherence_phase', 'imaginary_coherence', 'partial_coherence', 'directed_coherence'? Call spectral_connectivity.list_measures() to see the 35 available measures.

```

## Functional connectivity: coherence

The high-level `multitaper_connectivity` runs the multitaper transform and
returns a labeled `xarray.DataArray` with `(time, frequency, source, target)`
axes.

```python
>>> coherence = multitaper_connectivity(
...     time_series,
...     sampling_frequency=500,
...     method="coherence_magnitude",
...     time_halfbandwidth_product=3,
... )
>>> type(coherence).__name__
'DataArray'
>>> coherence.dims
('time', 'frequency', 'source', 'target')
>>> coherence.name
'coherence_magnitude'

```

## Read and slice the result

Because the output is labeled, you select by name rather than by axis index.
Default signal labels are the string indices `"0"`, `"1"`, `"2"`; pass
`signal_names` to use your own.

```python
>>> pair = coherence.sel(source="0", target="1")
>>> pair.dims
('time', 'frequency')
>>> band = coherence.sel(frequency=slice(30, 50))
>>> float(band.frequency.min()) >= 30.0
True

```

## Directed connectivity: spectral Granger

Directed measures are opt-in by name. The result reads `source -> target`:
`result.sel(source="A", target="B")` is the influence **from A to B**.

```python
>>> granger = multitaper_connectivity(
...     time_series,
...     sampling_frequency=500,
...     method="pairwise_spectral_granger_prediction",
...     signal_names=["A", "B", "C"],
...     time_halfbandwidth_product=3,
... )
>>> granger.coords["source"].values.tolist()
['A', 'B', 'C']
>>> a_to_b = granger.sel(source="A", target="B")
>>> a_to_b.dims
('time', 'frequency')

```

## Compute several measures at once

Pass a list of methods to get an `xarray.Dataset` with one variable per
measure. Shared spectra are cached, so this is cheaper than separate calls.

```python
>>> result = multitaper_connectivity(
...     time_series,
...     sampling_frequency=500,
...     method=["power", "coherence_magnitude"],
...     time_halfbandwidth_product=3,
... )
>>> type(result).__name__
'Dataset'
>>> sorted(result.data_vars)
['coherence_magnitude', 'power']
>>> result["power"].dims
('time', 'frequency', 'source')

```

## Group measures: canonical coherence between areas

Group measures (`list_measures(category="group_pairwise")` and
`list_measures(category="multivariate_components")`) compare *groups* of
signals, such as all channels in one brain area against all channels in
another. Pass `group_labels`, one label per signal naming its group; a
`group_pairwise` result is indexed by `source_group` and `target_group`
instead of by signal.

```python
>>> between_areas = multitaper_connectivity(
...     time_series,
...     sampling_frequency=500,
...     method="canonical_coherence",
...     group_labels=["CA1", "CA1", "PFC"],
... )
>>> between_areas.dims
('time', 'frequency', 'source_group', 'target_group')
>>> between_areas.coords["source_group"].values.tolist()
['CA1', 'PFC']

```

## Collapse into frequency bands

Pass `frequency_bands` to average (or integrate) each measure within named
bands. The `frequency` axis is replaced by a labeled `band` axis.

```python
>>> banded = multitaper_connectivity(
...     time_series,
...     sampling_frequency=500,
...     method="coherence_magnitude",
...     frequency_bands={"theta": (4, 8), "gamma": (30, 50)},
...     time_halfbandwidth_product=3,
... )
>>> banded.dims
('time', 'band', 'source', 'target')
>>> banded.coords["band"].values.tolist()
['theta', 'gamma']

```

## Bring your own Fourier coefficients

If you already have Fourier coefficients (e.g. from a wavelet transform), skip
the multitaper step and use `fourier_connectivity`. NumPy inputs may use the
`(observation, frequency, signal)` layout shown here.

```python
>>> coefficients = rng.standard_normal((20, 16, 2)) + 1j * rng.standard_normal(
...     (20, 16, 2)
... )
>>> frequencies = np.linspace(0, 250, 16)
>>> byo = fourier_connectivity(
...     coefficients, frequencies=frequencies, method="coherence_magnitude"
... )
>>> byo.dims
('time', 'frequency', 'source', 'target')
>>> byo.sizes["frequency"]
16

```

## Plug in your own transform

`Connectivity.from_transform` accepts any object that satisfies the
`SpectralTransform` protocol: an `fft()` method returning coefficients shaped
`(n_time_windows, n_trials, n_tapers, n_fft_samples, n_signals)`, plus
`frequencies` and `time`. No subclassing is needed. `help(SpectralTransform)`
lists the optional attributes such as `is_one_sided` (absent ones describe a
two-sided, unweighted, independent spectrum) and the scaling that makes
`power()` a density. `fft()` must return fresh, unshared storage on each call.
Neither the transform nor its caller may subsequently mutate it through any
alias, because `Connectivity` keeps it without copying. Read-only flags are
applied where the backend supports them. The example below is not scaled to a
density, so its `power()` is in arbitrary units; normalized measures such as
coherence are unaffected.

```python
>>> from spectral_connectivity import Connectivity, SpectralTransform
>>> class HannTransform:
...     """One Hann-windowed FFT per trial, non-negative frequencies only."""
...
...     is_one_sided = True  # optional; omit for a two-sided FFT-order spectrum
...
...     def __init__(self, time_series, sampling_frequency):
...         self.time_series = time_series  # (n_time_samples, n_trials, n_signals)
...         n_time_samples = time_series.shape[0]
...         self.frequencies = np.fft.rfftfreq(n_time_samples, d=1 / sampling_frequency)
...         self.time = np.array([n_time_samples / 2 / sampling_frequency])
...
...     def fft(self):
...         window = np.hanning(self.time_series.shape[0])[:, np.newaxis, np.newaxis]
...         coefficients = np.fft.rfft(window * self.time_series, axis=0)
...         # (frequency, trial, signal) -> (time, trial, taper, frequency, signal)
...         return coefficients.transpose(1, 0, 2)[np.newaxis, :, np.newaxis]
>>> transform = HannTransform(time_series, sampling_frequency=500)
>>> isinstance(transform, SpectralTransform)
True
>>> Connectivity.from_transform(transform).coherence_magnitude().shape
(1, 501, 3, 3)

```

## Pass a labeled DataArray

`time_series` may be an `xarray.DataArray`. Dimension names, not positions,
define the roles, so a `time` coordinate in seconds supplies the sampling
frequency and the signal labels become the `source`/`target` coordinates:

```python
>>> import xarray as xr
>>> da = xr.DataArray(
...     rng.standard_normal((1000, 4, 3)),
...     dims=("time", "trial", "signal"),
...     coords={"time": np.arange(1000) / 500, "signal": ["CA1", "CA3", "PFC"]},
... )
>>> coherence = multitaper_connectivity(da, method="coherence_magnitude")
>>> coherence.attrs["mt_sampling_frequency"]
500.0
>>> coherence.source.values.tolist()
['CA1', 'CA3', 'PFC']

```

For other dimension names, say which dimension plays which role:

```python
>>> custom = da.rename(time="t", trial="epoch", signal="channel")
>>> coherence = multitaper_connectivity(
...     custom,
...     method="coherence_magnitude",
...     time_dim="t",
...     trial_dim="epoch",
...     signal_dim="channel",
... )
>>> coherence.attrs["mt_sampling_frequency"], coherence.target.values.tolist()
(500.0, ['CA1', 'CA3', 'PFC'])

```

Common cases:

- Common time, trial and signal dimension names are recognized and transposed
  automatically; pass `time_dim`, `trial_dim` and `signal_dim` for others.
- A numeric `time` coordinate is elapsed seconds and supplies
  `sampling_frequency`. A numeric `sample` coordinate is sample numbers, which
  have no time scale, so pass `sampling_frequency` with it. Either one labels
  the output window centers, and a `sampling_frequency` you pass is checked
  against it.
- A 1-D index on the signal dimension, including its label type, becomes the
  `source`/`target` coordinates unless you pass `signal_names`.

If you hit an error or a warning:

- Ambiguous dimension names raise instead of falling back to axis position:
  name the roles with `time_dim`, `trial_dim` and `signal_dim`. When a single
  unrecognized dimension is left for the one remaining role, it is assigned by
  elimination and a warning names the mapping.
- Inferring the rate needs enough coordinate precision: pass
  `sampling_frequency` for low-precision or large-offset time coordinates.
- Signal labels must be unique, non-missing, NetCDF-compatible scalars
  (strings, real numbers, datetimes or timedeltas); integer labels must fit the
  signed 32-bit range for portable NetCDF3 files.
- Datetime, timedelta and object-valued time coordinates are not supported
  yet; convert them to elapsed seconds as below. (Datetime and timedelta
  *signal labels* are fine.)
- A dask-backed DataArray is rejected; call `.compute()` (or `.load()`) first.

```python
>>> stamped = da.assign_coords(
...     time=np.datetime64("2024-01-01T00:00", "ns") + np.arange(1000) * np.timedelta64(2, "ms")
... )
>>> elapsed = stamped.assign_coords(
...     time=(stamped.time - stamped.time[0]) / np.timedelta64(1, "s")
... )
>>> multitaper_connectivity(elapsed, method="coherence_magnitude").attrs["mt_sampling_frequency"]
500.0

```

## Where to go next

- Value ranges for every measure: `docs/CONNECTIVITY_METRIC_RANGES.md`.
- The full lower-level API: the `Connectivity` class (each `method` above is a
  method on it).
- End-to-end tutorials: the notebooks under `examples/`.