Repository: esi-neuroscience/syncopy

Issue:
Inter-trial coherence (ITC)
**Is your feature request related to a problem? Please elaborate.**

Intertrial coherence (ITC) is a measure of how consistent oscillatory phase is across an ensemble of trials.

Please provide a clear and concise description of the problem. If possible/applicable, please consider providing a minimal complete working example illustrating the problem.

https://www.mattcraddock.com/blog/2020/06/12/what-is-intertrial-coherence/

**Not really a problem per se, more like a missing feature? Please elaborate**
Please provide a description of what kind of functionality you want to see in Syncopy. If possible/applicable, please consider writing a few lines of (pseudo-)code to illustrate your thoughts.

This measure is usually the part of the connectivity analysis. I can provide the code but the need to be implemented into the connnectivity analysis

**Thank you for taking the time to help us improve Syncopy!**


Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index dfc92f07..157ac000 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,14 +1,14 @@
 # Changelog of SyNCoPy
 All notable changes to this project will be documented in this file.
 
-
-## [Unreleased]
+## [2022.12]
 
 ### NEW
-- time dependent coherence analysis #287
-- basic statistics (mean, std, var and median) for Syncopy data objects
-- timelockanalysis and new TimeLockData data type
-- PSTH method for SpikeData
+- time dependent coherence analysis
+- basic statistics (`spy.mean`, `spy.std`, `spy.var` and `spy.median`) for Syncopy data objects
+- `spy.timelockanalysis` and new `TimeLockData` data type
+- PSTH method for SpikeData - `spy.spike_psth`
+- inter trial coherence measure `spy.itc`
 - support for performing connectivity analysis from SpectralData (#364).
 - additional .info entries for Granger analysis, indicating details about the computation.
 - additional .info entries for FoooF results, e.g. Gaussian fit parameters.
diff --git a/syncopy/datatype/__init__.py b/syncopy/datatype/__init__.py
index 8e7a5140..f138e8dd 100644
--- a/syncopy/datatype/__init__.py
+++ b/syncopy/datatype/__init__.py
@@ -13,7 +13,6 @@ from .methods.padding import *
 from .methods.selectdata import *
 from .methods.show import *
 from .methods.copy import *
-from .methods.statistics import *
 
 # Populate local __all__ namespace
 __all__ = []
@@ -26,4 +25,3 @@ __all__.extend(methods.definetrial.__all__)
 __all__.extend(methods.selectdata.__all__)
 __all__.extend(methods.show.__all__)
 __all__.extend(methods.copy.__all__)
-__all__.extend(methods.statistics.__all__)
diff --git a/syncopy/datatype/continuous_data.py b/syncopy/datatype/continuous_data.py
index 0c988332..4d9b4e02 100644
--- a/syncopy/datatype/continuous_data.py
+++ b/syncopy/datatype/continuous_data.py
@@ -17,7 +17,6 @@ from collections.abc import Iterator
 # Local imports
 from .base_data import BaseData, FauxTrial, _definetrial
 from .methods.definetrial import definetrial
-from .methods import statistics
 from syncopy.shared.parsers import scalar_parser, array_parser
 from syncopy.shared.errors import SPYValueError, SPYWarning
 from syncopy.shared.tools import best_match
@@ -371,27 +370,6 @@ class ContinuousData(BaseData, ABC):
     def multipanelplot(self):
         raise NotImplementedError
 
-    # statistics
-    @statistics._attach_stat_doc(statistics.mean.__doc__)
-    def mean(self, dim, keeptrials=True, **kwargs):
-
-        return statistics.mean(self, dim, keeptrials, **kwargs)
-
-    @statistics._attach_stat_doc(statistics.std.__doc__)
-    def std(self, dim, keeptrials=True, **kwargs):
-
-        return statistics.std(self, dim, keeptrials, **kwargs)
-
-    @statistics._attach_stat_doc(statistics.var.__doc__)    
-    def var(self, dim, keeptrials=True, **kwargs):
-
-        return statistics.var(self, dim, keeptrials, **kwargs)
-
-    @statistics._attach_stat_doc(statistics.median.__doc__)        
-    def median(self, dim, keeptrials=True, **kwargs):
-
-        return statistics.median(self, dim, keeptrials, **kwargs)
-
 
 class AnalogData(ContinuousData):
     """Multi-channel, uniformly-sampled, analog (real float) data
diff --git a/syncopy/datatype/methods/show.py b/syncopy/datatype/methods/show.py
index de87a862..b02ca922 100644
--- a/syncopy/datatype/methods/show.py
+++ b/syncopy/datatype/methods/show.py
@@ -164,8 +164,8 @@ def show(data, squeeze=True, **kwargs):
     # Use an object's `_preview_trial` method fetch required indexing tuples
     idxList = []
     for trlno in data.selection.trial_ids:
-        # each dim has an entry
-        idxs = data._preview_trial(trlno).idx
+        # each dim has an entry, list only needed for mutability
+        idxs = list(data._preview_trial(trlno).idx)
 
         # time/toi is a special case, all other dims get checked
         # beforehand, e.g. foi, channel, ... but out of range toi's get mapped
@@ -177,7 +177,12 @@ def show(data, squeeze=True, **kwargs):
                 lgl = "valid `toi` selection"
                 act = sel
                 raise SPYValueError(lgl, 'show kwargs', act)
-        idxList.append(idxs)
+
+        for i, prop_idx in enumerate(idxs):            
+            if isinstance(prop_idx, list) and len(prop_idx) == 1:
+                idxs[i] = prop_idx[0]
+            
+        idxList.append(tuple(idxs))
 
     # Reset in-place subset selection
     data.selection = None
diff --git a/syncopy/plotting/sp_plotting.py b/syncopy/plotting/sp_plotting.py
index 53b0b264..c7f28050 100644
--- a/syncopy/plotting/sp_plotting.py
+++ b/syncopy/plotting/sp_plotting.py
@@ -176,9 +176,14 @@ def plot_SpectralData(data, **show_kwargs):
             data_y = data.show(**show_kwargs)
             ylabel = f'{output}'
 
+        # for itc.. needs to be improved
+        if output is None:
+            ylabel = ''
+
         # flip if required
-        if data_y.shape[1] == len(data_x):
-            data_y = data_y.T
+        if data_y.ndim > 1:
+            if data_y.shape[1] == len(data_x):
+                data_y = data_y.T
 
         fig, ax = _plotting.mk_line_figax(xlabel='frequency (Hz)',
                                           ylabel=ylabel)
diff --git a/syncopy/shared/input_processors.py b/syncopy/shared/input_processors.py
index a9a0cb73..95869076 100644
--- a/syncopy/shared/input_processors.py
+++ b/syncopy/shared/input_processors.py
@@ -285,8 +285,8 @@ def process_taper(taper,
 
         # direct mtm estimate (averaging) only valid for spectral power
         if not keeptapers and output != "pow":
-            lgl = (f"'pow'|False or '{output}'|True, set keeptapers=True"
-                   "OR `output='pow'`!")
+            lgl = (f"'pow'|False or '{output}'|True, set either keeptapers=True "
+                   "or `output='pow'`!")
             raise SPYValueError(legal=lgl, varname="output|keeptapers", actual=f"'{output}'|{keeptapers}")
 
         # --- minimal smoothing bandwidth ---
diff --git a/syncopy/statistics/__init__.py b/syncopy/statistics/__init__.py
index 61c2d568..f26a358b 100644
--- a/syncopy/statistics/__init__.py
+++ b/syncopy/statistics/__init__.py
@@ -4,8 +4,17 @@
 # 
 
 # Import __all__ routines from local modules
+from . import summary_stats
 from .spike_psth import spike_psth
 from .timelockanalysis import timelockanalysis
+from .summary_stats import (
+    mean,
+    var,
+    std,
+    median,
+    itc,
+)
 
 # Populate local __all__ namespace
 __all__ = ['spike_psth', 'timelockanalysis']
+__all__.extend(summary_stats.__all__)
diff --git a/syncopy/datatype/methods/statistics.py b/syncopy/statistics/summary_stats.py
similarity index 86%
rename from syncopy/datatype/methods/statistics.py
rename to syncopy/statistics/summary_stats.py
index a99c6233..3a95c31d 100644
--- a/syncopy/datatype/methods/statistics.py
+++ b/syncopy/statistics/summary_stats.py
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 #
-# Syncopy object simple statistics (mean, std, ...)
+# Syncopy object simple summary statistics (mean, std, ...)
 #
 
 
@@ -14,7 +14,7 @@ from syncopy.shared.errors import SPYValueError, SPYWarning
 from syncopy.shared.computational_routine import ComputationalRoutine
 from syncopy.shared.kwarg_decorators import process_io, unwrap_select, detect_parallel_client
 
-__all__ = ['mean', 'std', 'var', 'median']
+__all__ = ['mean', 'std', 'var', 'median', 'itc']
 
 
 @unwrap_select
@@ -154,7 +154,7 @@ def median(spy_data, dim, keeptrials=True, **kwargs):
     Returns
     -------
     res : Syncopy data object
-        New object with the desired dimension averaged out
+        New object with the median of the desired dimension
 
     """
 
@@ -165,6 +165,52 @@ def median(spy_data, dim, keeptrials=True, **kwargs):
                        keeptrials=keeptrials,
                        **kwargs)
 
+@unwrap_select
+def itc(spec_data, **kwargs):
+    """
+    Calculates the inter trial coherence for a
+    SpectralData ``spec_data`` object, the input
+    spectrum needs to be complex.
+    The ITC of N trials is given by the length
+    of the complex mean of vectors z_i(f):
+
+        1/N \sum z_i / |z_i|
+
+    and have therefore values between 0 and 1.
+    In the literature this measure is also often
+    called the `Kuramoto order parameter`.
+
+    For time-frequency spectra the trial sizes
+    have to match exactly, and the output will
+    be also additionally time dependent.
+
+    Parameters
+    ----------
+    spec_data : :class:`~syncopy.SpectralData`
+        The input spectrum, needs at least 2 trials
+
+    Returns
+    -------
+    res  : :class:`~syncopy.SpectralData`
+        The frequency dependent order parameters,
+        the inter trial coherence
+    """
+
+    data_parser(spec_data,
+                varname='spec_data',
+                dataclass='SpectralData',
+                empty=False)
+
+    if spec_data.data.dtype != np.complex64 and spec_data.data.dtype != np.complex128:
+        lgl = "complex valued spectra, set `output='fourier` in spy.freqanalysis!"
+        act = "real valued spectral data"
+        raise SPYValueError(lgl, 'spec_data', act)
+
+    # takes care of remaining checks
+    res = _trial_statistics(spec_data, operation='itc')
+
+    return res
+
 
 def _statistics(spy_data, operation, dim, keeptrials=True, **kwargs):
 
@@ -395,6 +441,12 @@ def _trial_statistics(in_data, operation='mean'):
         in_data.selectdata(inplace=True)
         in_data.selection._cleanup = True
 
+    nTrials = len(in_data.selection.trials)
+    if nTrials <= 1:
+        lgl = "at least 2 trials"
+        act = f"got {nTrials} trials"
+        raise SPYValueError(lgl, 'in_data', act)
+
     # we always have at least one (all-to-all) trial selection
     out_shape = in_data.selection.trials[0].shape
 
@@ -417,6 +469,9 @@ def _trial_statistics(in_data, operation='mean'):
         result = _trial_var(in_data, result)
     elif operation == 'std':
         result = np.sqrt(_trial_var(in_data, result))
+    elif operation == 'itc':
+        result = np.abs(_trial_circ_average(in_data, result))
+
     # there is no apparent clever way to achieve
     # this efficiently over multiple dimensions
     elif operation == 'median':
@@ -463,7 +518,7 @@ def _trial_average(in_data, out_arr):
 
     trials = in_data.selection.trials
     for trl in trials:
-        out_arr = np.nansum([out_arr, trl], axis=0)
+        out_arr += trl
 
     # normalize
     out_arr /= len(trials)
@@ -491,7 +546,7 @@ def _trial_var(in_data, out_arr):
     trials = in_data.selection.trials
     for trl in trials:
         # absolute value for complex numbers
-        out_arr = np.sum([out_arr, (np.abs(trl - average))**2], axis=0)
+        out_arr += np.abs(trl - average)**2
 
     # normalize
     out_arr /= len(trials)
@@ -499,8 +554,37 @@ def _trial_var(in_data, out_arr):
     return out_arr
 
 
+def _trial_circ_average(in_data, out_arr):
+    """
+    Sequential complex average, can be used for
+    inter trial coherence (order parameter) or mean phase estimation.
+    Shape checking and dealing with selections is done in `_trial_statistics`.
+
+    Parameters
+    ----------
+    in_data : Syncopy data object
+        To get a fresh trial indexer instance, pointing to the trial arrays
+    out_arr : np.ndarray
+        The empty NumPy array of correct shape to collect the results
+    """
+
+    trials = in_data.selection.trials
+    for trl in trials:
+        # add cartesian unit vectors
+        out_arr += trl / np.abs(trl)
+
+    # normalize
+    out_arr /= len(trials)
+
+    # and return complex resultant vector
+
+    return out_arr
+
+
 def _attach_stat_doc(orig_doc):
     """
+    NOT USED ATM - could be useful for other method doc strings
+
     This is a helper to attach the full doc to the statistical methods in ContinuousData.
     Including the `select` and `parallel` sections from the kwarg decorators,
     which can/should only be applied once.


Test patch:
diff --git a/syncopy/tests/local_spy.py b/syncopy/tests/local_spy.py
index 37b755d9..61e8dbbf 100644
--- a/syncopy/tests/local_spy.py
+++ b/syncopy/tests/local_spy.py
@@ -57,10 +57,10 @@ if __name__ == "__main__":
     specf2 = spy.freqanalysis(adata, tapsmofrq=2, keeptrials=False, foi=foi,
                               output="fooof_peaks", fooof_opt={'max_n_peaks': 2})
 
-    print("Start: Testing parallel computation of mtmfft")
-    spec4 = spy.freqanalysis(adata, tapsmofrq=2, keeptrials=True, foi=foi, parallel=True, output="pow")
-    print("End: Testing parallel computation of mtmfft")
+    # print("Start: Testing parallel computation of mtmfft")
+    # spec4 = spy.freqanalysis(adata, tapsmofrq=2, keeptrials=True, foi=foi, parallel=True, output="pow")
+    # print("End: Testing parallel computation of mtmfft")
 
     #spec.singlepanelplot()
     #specf.singlepanelplot()
-    #specf2.singlepanelplot()S
\ No newline at end of file
+    #specf2.singlepanelplot()S
diff --git a/syncopy/tests/test_continuousdata.py b/syncopy/tests/test_continuousdata.py
index f1902f0d..ae63bcc0 100644
--- a/syncopy/tests/test_continuousdata.py
+++ b/syncopy/tests/test_continuousdata.py
@@ -1152,228 +1152,9 @@ class TestTimeLockData:
         assert np.array_equal(avg_data3, tld2.avg)
 
 
-class TestStatistics:
-
-    # initialize rng instance
-    rng = np.random.default_rng(helpers.test_seed)
-
-    # lognormal distribution parameters
-    mu, sigma = 2, .5
-
-    nTrials = 4
-    nChannels = 3
-    nSamples = 10
-    nFreq = 10
-    nTaper = 2
-
-    ln_samples = rng.lognormal(mu, sigma, size=(nTrials, nSamples, nChannels))
-    adata = AnalogData(data=[trl for trl in ln_samples], samplerate=1)
-
-    ln_samples = rng.lognormal(mu, sigma, size=(nTrials, nSamples, nTaper, nFreq, nChannels))
-    spec_data = SpectralData(data=[trl for trl in ln_samples], samplerate=1)
-
-    ln_samples = rng.lognormal(mu, sigma, size=(nTrials, nSamples, nFreq, nChannels, nChannels))
-    crossspec_data = CrossSpectralData(data=[trl for trl in ln_samples], samplerate=1)
-
-    data_types = [adata, spec_data, crossspec_data]
-
-    def test_dim_statistics(self):
-        """
-        Tests statistics over dimensions, not trials
-        """
-
-        # check only 2nd trial
-        test_trial = 1
-
-        for spy_data in self.data_types:
-            for dim in spy_data.dimord:
-
-                # get index of dimension
-                axis = spy_data.dimord.index(dim)
-
-                # -- test average --
-
-                # data class method
-                spy_res1 = spy_data.mean(dim=dim)
-                # top-level function
-                spy_res2 = spy.mean(spy_data, dim=dim)
-                # check only one trial
-                npy_res = np.mean(spy_data.trials[test_trial], axis=axis)
-
-                self._check_trial('mean', dim, spy_res1, npy_res, trial=test_trial)
-                self._check_trial('mean', dim, spy_res2, npy_res, trial=test_trial)
-
-                # --- test variance ---
-
-                # data class method
-                spy_res1 = spy_data.var(dim=dim)
-                # top-level function
-                spy_res2 = spy.var(spy_data, dim=dim)
-                # check only one trial
-                npy_res = np.var(spy_data.trials[test_trial], axis=axis)
-
-                self._check_trial('var', dim, spy_res1, npy_res, trial=test_trial)
-                self._check_trial('var', dim, spy_res2, npy_res, trial=test_trial)
-
-                # --- test standard deviation ---
-
-                # data class method
-                spy_res1 = spy_data.std(dim=dim)
-                # top-level function
-                spy_res2 = spy.std(spy_data, dim=dim)
-                # check only one trial
-                npy_res = np.std(spy_data.trials[test_trial], axis=axis)
-
-                self._check_trial('std', dim, spy_res1, npy_res, trial=test_trial)
-                self._check_trial('std', dim, spy_res2, npy_res, trial=test_trial)
-
-                # --- test median ---
-
-                # data class method
-                spy_res1 = spy_data.median(dim=dim)
-                # top-level function
-                spy_res2 = spy.median(spy_data, dim=dim)
-                # check only one trial
-                npy_res = np.median(spy_data.trials[test_trial], axis=axis)
-
-                self._check_trial('median', dim, spy_res1, npy_res, trial=test_trial)
-                self._check_trial('median', dim, spy_res2, npy_res, trial=test_trial)
-
-    def _check_trial(self, operation, dim, spy_res, npy_res, trial=1):
-        """
-        Test that direct numpy stats give the same results
-        for a single trial
-        """
-
-        # show returns list of trials, pick only one
-        # show also squeezes out the singleton dimension
-        # which remains after the statistic got computed!
-        trial_result_spy = spy_res.show()[trial]
-        assert np.allclose(trial_result_spy, npy_res)
-
-        # check the dimension label was set to the statistical operation
-        if dim not in ['freq', 'time']:
-            assert getattr(spy_res, dim) == operation
-        # numerical dimension labels get set to 0 (axis is gone)
-        elif dim == 'time':
-            assert spy_res.time[trial] == 0
-        elif dim == 'freq':
-            assert spy_res.freq == 0
-
-    def test_trial_statistics(self):
-        """
-        Test statistics over trials
-        """
-
-        # --- test statistics trial average against CR trial average ---
-
-        spec = spy.freqanalysis(self.adata, keeptrials=True)
-        # trigger trial average after spectral estimation
-        spec1 = spec.mean(dim='trials')
-        spec1a = spy.mean(spec, dim='trials')
-
-        # trial average via CR keeptrials
-        spec2 = spy.freqanalysis(self.adata, keeptrials=False)
-
-        assert len(spec1.trials) == 1
-        assert np.allclose(spec1.data, spec2.data)
-        assert np.allclose(spec1a.data, spec2.data)
-
-        # --- test trial var and std ---
-        for spy_data in self.data_types:
-
-            spy_var = spy_data.var(dim='trials')
-
-            # reshape to get rid of trial stacking along time axis
-            # array has shape (nTrials, nSamples, ..rest-of-dims..)
-            arr = spy_data.data[()].reshape(self.nTrials, self.nSamples, *spy_data.data.shape[1:])
-            # now compute directly over the trial axis
-            npy_var = np.var(arr, axis=0)
-
-            assert len(spy_var.trials) == 1
-            assert np.allclose(npy_var, spy_var.data)
-
-            spy_std = spy_data.std(dim='trials')
-
-            # reshape to get rid of trial stacking along time axis
-            # array has shape (nTrials, nSamples, ..rest-of-dims..)
-            arr = spy_data.data[()].reshape(self.nTrials, self.nSamples, *spy_data.data.shape[1:])
-            # now compute directly over the trial axis
-            npy_std = np.std(arr, axis=0)
-
-            assert len(spy_var.trials) == 1
-            assert np.allclose(npy_std, spy_std.data)
-
-    def test_selections(self):
-
-        # got 10 samples with 1s samplerate,so time is [-1, ..., 8]
-        sdict1 = {'trials': [1, 3], 'latency': [2, 6]}
-        res = self.adata.mean(dim='channel', select=sdict1)
-        assert len(res.trials) == 2
-        assert self.adata.time[0].min() == -1
-        assert res.time[0].min() == 2
-        assert self.adata.time[0].max() == 8
-        assert res.time[0].max() == 6
-
-        # freq axis is [0, ..., 9]
-        sdict2 = {'channel': [0, 2], 'frequency': [1, 5]}
-        res = self.spec_data.var(dim='trials', select=sdict2)
-        assert np.all(res.channel == np.array(['channel1', 'channel3']))
-        assert np.all(res.freq == np.arange(1, 6))
-
-        # check at least a few times that the statistics are indeed
-        # computed correctly on the trimmed down data
-        sdict3 = {'trials': [1, 3]}
-        res = spy.mean(self.crossspec_data, dim='trials', select=sdict3)
-        # reshape to extract trial separated arrays
-        arr = self.crossspec_data.data[()].reshape(self.nTrials, self.nSamples,
-                                                   self.nFreq, self.nChannels, self.nChannels)
-        # now cut out the same 2 trials and average
-        npy_res = arr[1::2].mean(axis=0)
-        assert np.allclose(npy_res, res.data)
-
-        sdict4 = {'channel': [0, 2]}
-        res = spy.mean(self.spec_data, dim='channel', select=sdict4)
-        # now cut out the same 2 channels and average, dimord is (time, taper, freq, channel)
-        npy_res = self.spec_data.data[..., ::2].mean(axis=-1)
-        # check only 1st trial
-        assert np.allclose(npy_res[:self.nSamples], res.show(trials=0))
-
-        # one last time for the freq axis
-        sdict5 = {'frequency': [1, 4]}
-        res = spy.median(self.spec_data, dim='freq', select=sdict5)
-        # cut out same frequencies directly from the dataset array
-        npy_res = np.median(self.spec_data.data[..., 1:5, :], axis=2)
-        # check only 2nd trial
-        assert np.allclose(npy_res[self.nSamples:2 * self.nSamples], res.show(trials=1))
-
-    def test_exceptions(self):
-
-        with pytest.raises(SPYValueError) as err:
-            spy.mean(self.adata, dim='sth')
-        assert "expected one of ['time', 'channel']" in str(err.value)
-
-        # unequal trials and trial average can't work
-        with pytest.raises(SPYValueError) as err:
-            # to not screw sth up
-            adata_cpy = spy.copy(self.adata)
-            trldef = adata_cpy.trialdefinition
-            trldef[2] = [21, 25, -1]
-            adata_cpy.trialdefinition = trldef
-            spy.mean(adata_cpy, dim='trials')
-        assert "found trials of different shape" in str(err.value)
-
-    def test_stat_parallel(self, testcluster=None):
-        client = dd.Client(testcluster)
-        self.test_selections()
-        # should have no effect here
-        self.test_trial_statistics()
-
-
 if __name__ == '__main__':
 
-    T0 = TestAnalogData()
-    T1 = TestSpectralData()
-    T2 = TestTimeLockData()
-    T3 = TestStatistics()
+    T1 = TestAnalogData()
+    T2 = TestSpectralData()
+    T3 = TestTimeLockData()
     T4 = TestCrossSpectralData()
diff --git a/syncopy/tests/test_selectdata.py b/syncopy/tests/test_selectdata.py
index 506a7318..08cac004 100644
--- a/syncopy/tests/test_selectdata.py
+++ b/syncopy/tests/test_selectdata.py
@@ -505,6 +505,9 @@ class TestSelector():
                         else:
                             idx = [slice(None)] * len(dummy.dimord)
                             idx[dummy.dimord.index(prop)] = solution
+                            print(idx, solution, prop)
+                            print(np.array(dummy.data)[tuple(idx)])
+                            print(selected.data[()])
                             assert np.array_equal(np.array(dummy.data)[tuple(idx)],
                                                   selected.data)
                             assert np.array_equal(getattr(selected, prop),
diff --git a/syncopy/tests/test_statistics.py b/syncopy/tests/test_statistics.py
new file mode 100644
index 00000000..e624bd7f
--- /dev/null
+++ b/syncopy/tests/test_statistics.py
@@ -0,0 +1,309 @@
+# -*- coding: utf-8 -*-
+#
+# Test proper functionality of Syncopy's `ContinuousData` class + subclasses
+#
+
+# Builtin/3rd party package imports
+import pytest
+import numpy as np
+import dask.distributed as dd
+import matplotlib.pyplot as ppl
+
+# Local imports
+import syncopy as spy
+from syncopy.datatype import AnalogData, SpectralData, CrossSpectralData, TimeLockData
+from syncopy.shared.errors import SPYValueError, SPYTypeError
+from syncopy.shared.tools import StructDict
+from syncopy.tests import helpers
+from syncopy.tests import synth_data as sd
+
+
+class TestStatistics:
+
+    # initialize rng instance
+    rng = np.random.default_rng(helpers.test_seed)
+
+    # lognormal distribution parameters
+    mu, sigma = 2, .5
+
+    nTrials = 4
+    nChannels = 3
+    nSamples = 10
+    nFreq = 10
+    nTaper = 2
+
+    ln_samples = rng.lognormal(mu, sigma, size=(nTrials, nSamples, nChannels))
+    adata = AnalogData(data=[trl for trl in ln_samples], samplerate=1)
+
+    ln_samples = rng.lognormal(mu, sigma, size=(nTrials, nSamples, nTaper, nFreq, nChannels))
+    spec_data = SpectralData(data=[trl for trl in ln_samples], samplerate=1)
+
+    ln_samples = rng.lognormal(mu, sigma, size=(nTrials, nSamples, nFreq, nChannels, nChannels))
+    crossspec_data = CrossSpectralData(data=[trl for trl in ln_samples], samplerate=1)
+
+    data_types = [adata, spec_data, crossspec_data]
+
+    def test_dim_statistics(self):
+        """
+        Tests statistics over dimensions, not trials
+        """
+
+        # check only 2nd trial
+        test_trial = 1
+
+        for spy_data in self.data_types:
+            for dim in spy_data.dimord:
+
+                # get index of dimension
+                axis = spy_data.dimord.index(dim)
+
+                # -- test average --
+
+                # top-level function
+                spy_res1 = spy.mean(spy_data, dim=dim)
+                # check only one trial
+                npy_res = np.mean(spy_data.trials[test_trial], axis=axis)
+
+                self._check_trial('mean', dim, spy_res1, npy_res, trial=test_trial)
+
+                # --- test variance ---
+
+                # top-level function
+                spy_res1 = spy.var(spy_data, dim=dim)
+                # check only one trial
+                npy_res = np.var(spy_data.trials[test_trial], axis=axis)
+
+                self._check_trial('var', dim, spy_res1, npy_res, trial=test_trial)
+
+                # --- test standard deviation ---
+
+                # top-level function
+                spy_res1 = spy.std(spy_data, dim=dim)
+                # check only one trial
+                npy_res = np.std(spy_data.trials[test_trial], axis=axis)
+
+                self._check_trial('std', dim, spy_res1, npy_res, trial=test_trial)
+
+                # --- test median ---
+
+                # top-level function
+                spy_res1 = spy.median(spy_data, dim=dim)
+                # check only one trial
+                npy_res = np.median(spy_data.trials[test_trial], axis=axis)
+
+                self._check_trial('median', dim, spy_res1, npy_res, trial=test_trial)
+
+    def _check_trial(self, operation, dim, spy_res, npy_res, trial=1):
+        """
+        Test that direct numpy stats give the same results
+        for a single trial
+        """
+
+        # show returns list of trials, pick only one
+        # show also squeezes out the singleton dimension
+        # which remains after the statistic got computed!
+        trial_result_spy = spy_res.show()[trial]
+        assert np.allclose(trial_result_spy, npy_res)
+
+        # check the dimension label was set to the statistical operation
+        if dim not in ['freq', 'time']:
+            assert getattr(spy_res, dim) == operation
+        # numerical dimension labels get set to 0 (axis is gone)
+        elif dim == 'time':
+            assert spy_res.time[trial] == 0
+        elif dim == 'freq':
+            assert spy_res.freq == 0
+
+    def test_trial_statistics(self):
+        """
+        Test statistics over trials
+        """
+
+        # --- test statistics trial average against CR trial average ---
+
+        spec = spy.freqanalysis(self.adata, keeptrials=True)
+        # trigger trial average after spectral estimation
+        spec1 = spy.mean(spec, dim='trials')
+
+        # trial average via CR keeptrials
+        spec2 = spy.freqanalysis(self.adata, keeptrials=False)
+
+        assert len(spec1.trials) == 1
+        assert np.allclose(spec1.data, spec2.data)
+
+        # --- test trial var and std ---
+        for data in self.data_types:
+
+            spy_var = spy.var(data, dim='trials')
+
+            # reshape to get rid of trial stacking along time axis
+            # array has shape (nTrials, nSamples, ..rest-of-dims..)
+            arr = data.data[()].reshape(self.nTrials, self.nSamples, *data.data.shape[1:])
+            # now compute directly over the trial axis
+            npy_var = np.var(arr, axis=0)
+
+            assert len(spy_var.trials) == 1
+            assert np.allclose(npy_var, spy_var.data)
+
+            spy_std = spy.std(data, dim='trials')
+
+            # reshape to get rid of trial stacking along time axis
+            # array has shape (nTrials, nSamples, ..rest-of-dims..)
+            arr = data.data[()].reshape(self.nTrials, self.nSamples, *data.data.shape[1:])
+            # now compute directly over the trial axis
+            npy_std = np.std(arr, axis=0)
+
+            assert len(spy_var.trials) == 1
+            assert np.allclose(npy_std, spy_std.data)
+
+    def test_selections(self):
+
+        # got 10 samples with 1s samplerate,so time is [-1, ..., 8]
+        sdict1 = {'trials': [1, 3], 'latency': [2, 6]}
+        res = spy.mean(self.adata, dim='channel', select=sdict1)
+        assert len(res.trials) == 2
+        assert self.adata.time[0].min() == -1
+        assert res.time[0].min() == 2
+        assert self.adata.time[0].max() == 8
+        assert res.time[0].max() == 6
+
+        # freq axis is [0, ..., 9]
+        sdict2 = {'channel': [0, 2], 'frequency': [1, 5]}
+        res = spy.var(self.spec_data, dim='trials', select=sdict2)
+        assert np.all(res.channel == np.array(['channel1', 'channel3']))
+        assert np.all(res.freq == np.arange(1, 6))
+
+        # check at least a few times that the statistics are indeed
+        # computed correctly on the trimmed down data
+        sdict3 = {'trials': [1, 3]}
+        res = spy.mean(self.crossspec_data, dim='trials', select=sdict3)
+        # reshape to extract trial separated arrays
+        arr = self.crossspec_data.data[()].reshape(self.nTrials, self.nSamples,
+                                                   self.nFreq, self.nChannels, self.nChannels)
+        # now cut out the same 2 trials and average
+        npy_res = arr[1::2].mean(axis=0)
+        assert np.allclose(npy_res, res.data)
+
+        sdict4 = {'channel': [0, 2]}
+        res = spy.mean(self.spec_data, dim='channel', select=sdict4)
+        # now cut out the same 2 channels and average, dimord is (time, taper, freq, channel)
+        npy_res = self.spec_data.data[..., ::2].mean(axis=-1)
+        # check only 1st trial
+        assert np.allclose(npy_res[:self.nSamples], res.show(trials=0))
+
+        # one last time for the freq axis
+        sdict5 = {'frequency': [1, 4]}
+        res = spy.median(self.spec_data, dim='freq', select=sdict5)
+        # cut out same frequencies directly from the dataset array
+        npy_res = np.median(self.spec_data.data[..., 1:5, :], axis=2)
+        # check only 2nd trial
+        assert np.allclose(npy_res[self.nSamples:2 * self.nSamples], res.show(trials=1))
+
+    def test_exceptions(self):
+
+        with pytest.raises(SPYValueError) as err:
+            spy.mean(self.adata, dim='sth')
+        assert "expected one of ['time', 'channel']" in str(err.value)
+
+        # unequal trials and trial average can't work
+        with pytest.raises(SPYValueError) as err:
+            # to not screw sth up
+            adata_cpy = spy.copy(self.adata)
+            trldef = adata_cpy.trialdefinition
+            trldef[2] = [21, 25, -1]
+            adata_cpy.trialdefinition = trldef
+            spy.mean(adata_cpy, dim='trials')
+        assert "found trials of different shape" in str(err.value)
+
+    def test_stat_parallel(self, testcluster=None):
+        client = dd.Client(testcluster)
+        self.test_selections()
+        # should have no effect here
+        self.test_trial_statistics()
+
+    def test_itc(self, do_plot=True):
+
+        adata = sd.white_noise(100,
+                               nSamples=1000,
+                               nChannels=2,
+                               samplerate=500)
+
+        # add simple 60Hz armonic
+        adata += sd.harmonic(100,
+                             freq=60,
+                             nSamples=1000,
+                             nChannels=2,
+                             samplerate=500)
+
+        trials = []
+        # add frequency drift along trials
+        # note same initial phase
+        # so ITC will be ~1 at time 0
+        freq = 30
+        dfreq = 2 / 100  # frequency difference
+        for trl in adata.trials:
+            # start at 0s
+            dat = np.cos(2 * np.pi * freq * (adata.time[0] + 1))
+            trials.append(trl + np.c_[dat, dat])
+            freq += dfreq
+        adata = spy.AnalogData(data=trials, samplerate=500)
+
+        tf_spec = spy.freqanalysis(adata, method='mtmconvol',
+                                   t_ftimwin=0.5,
+                                   output='fourier',
+                                   foilim=[0, 100])
+
+        spec = spy.freqanalysis(adata, foilim=[0, 100], output='fourier')
+
+        # -- calculate itc --
+        itc = spy.itc(spec)
+        tf_itc = spy.itc(tf_spec)
+
+        assert isinstance(tf_itc, spy.SpectralData)
+
+        assert np.all(np.imag(itc.data[()]) == 0)
+        assert itc.data[()].max() <= 1
+        assert itc.data[()].min() >= 0
+
+        # high itc around the in phase 60Hz
+        assert np.all(itc.show(frequency=60) > 0.99)
+        # low (time averaged) itc around the drifters
+        assert np.all(itc.show(frequency=30) < 0.4)
+
+        assert np.all(np.imag(tf_itc.data[()]) == 0)
+        assert tf_itc.data[()].max() <= 1
+        assert tf_itc.data[()].min() >= 0        
+        assert np.allclose(tf_itc.time[0], tf_spec.time[0])
+
+        if do_plot:
+
+            # plot tf power spectrum
+            # after 500 samples / 1s
+            fig, ax = ppl.subplots()
+            for idx in [0, 50, 99]:
+                power500 = np.abs(tf_spec.trials[idx][500, 0, :, 0])
+                ax.plot(tf_spec.freq, power500, label=f"trial {idx}")
+            # note that the 30Hz peak wandered to ~31Hz
+            # hence the phases will decorrelate over time
+            ax.legend()
+            ax.set_xlabel('frequency (Hz)')
+
+            # plot ITCs
+            itc.singlepanelplot(channel=1)
+            tf_itc.singlepanelplot(channel=1)
+
+            # plot tf ITC time profiles
+            fig, ax = ppl.subplots()
+            ax.set_xlabel('time (s)')
+            for frq in [31, 60, 10]:
+                itc_profile = tf_itc.show(frequency=frq, channel=0)
+                ax.plot(tf_itc.time[0], itc_profile, label=f"{frq}Hz")
+            ax.legend()
+            ax.set_title("time dependent ITC")
+
+        return itc, spec
+
+
+if __name__ == '__main__':
+
+    T1 = TestStatistics()