From b2562be90fde58d135763ae3dee0ad8b8c9dcf60 Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 12:58:11 -0400
Subject: [PATCH 01/16] add pycircstat and track-linearization

---
 pyproject.toml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/pyproject.toml b/pyproject.toml
index c6c3bb7..c8f6149 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -39,6 +39,8 @@ dependencies = [
     "neo>=0.13.3",
     "quantities>=0.15.0",
     "neurodsp>=2.2.1",
+    "pycircstat>=0.0.2",
+    "track-linearization>=2.3.1"
 ]
 
 [project.urls]

From 6967f7124a7d657e14d8f2d3a0ad8f408d2e8234 Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 13:19:56 -0400
Subject: [PATCH 02/16] remove pycircstats in favor of minimally usable local
 version

---
 neuro_py/process/precession_utils.py |  38 +-
 neuro_py/stats/circ_stats.py         | 548 +++++++++++++++++++++++++++
 pyproject.toml                       |   1 -
 3 files changed, 570 insertions(+), 17 deletions(-)
 create mode 100644 neuro_py/stats/circ_stats.py

diff --git a/neuro_py/process/precession_utils.py b/neuro_py/process/precession_utils.py
index d419c44..580e023 100644
--- a/neuro_py/process/precession_utils.py
+++ b/neuro_py/process/precession_utils.py
@@ -2,13 +2,13 @@
 
 import numba
 import numpy as np
-import pycircstat as pcs
 import pyfftw
 import scipy as sp
-
 from lazy_loader import attach as _attach
-from scipy.signal import find_peaks
 from scipy.ndimage import gaussian_filter1d
+from scipy.signal import find_peaks
+
+import neuro_py.stats.circ_stats as pcs
 
 __all__ = (
     "corrcc",
@@ -26,6 +26,7 @@
 # https://github.com/seqasim/human_precession/blob/main/Precession_utils.py
 # https://doi.org/10.1016/j.cell.2021.04.017
 
+
 def corrcc(alpha1, alpha2, axis=None):
     """
     Circular correlation coefficient for two circular random variables.
@@ -57,7 +58,7 @@ def corrcc(alpha1, alpha2, axis=None):
     n = len(alpha1)
 
     # center data on circular mean
-    alpha1_centered, alpha2_centered = pcs.descriptive.center(alpha1, alpha2, axis=axis)
+    alpha1_centered, alpha2_centered = pcs.center(alpha1, alpha2, axis=axis)
 
     num = np.sum(np.sin(alpha1_centered) * np.sin(alpha2_centered), axis=axis)
     den = np.sqrt(
@@ -114,14 +115,14 @@ def corrcc_uniform(alpha1, alpha2, axis=None):
     n = len(alpha1)
 
     # center data on circular mean
-    alpha1_centered, alpha2_centered = pcs.descriptive.center(alpha1, alpha2, axis=axis)
+    alpha1_centered, alpha2_centered = pcs.center(alpha1, alpha2, axis=axis)
 
     # One of the sample means is not well defined due to uniform distribution of data
     # so take the difference of the resultant vector length for the sum and difference
     # of the alphas
-    num = pcs.descriptive.resultant_vector_length(
+    num = pcs.resultant_vector_length(
         alpha1 - alpha2
-    ) - pcs.descriptive.resultant_vector_length(alpha1 + alpha2)
+    ) - pcs.resultant_vector_length(alpha1 + alpha2)
     den = 2 * np.sqrt(
         np.sum(np.sin(alpha1_centered) ** 2, axis=axis)
         * np.sum(np.sin(alpha2_centered) ** 2, axis=axis)
@@ -206,8 +207,8 @@ def myfun1(p):
     linear_circ = np.mod(abs(sl) * lin, 2 * np.pi)
 
     # # marginal distributions:
-    p1, z1 = pcs.tests.rayleigh(circ)
-    p2, z2 = pcs.tests.rayleigh(linear_circ)
+    p1, z1 = pcs.rayleigh(circ)
+    p2, z2 = pcs.rayleigh(linear_circ)
 
     # circular-linear correlation:
     if (p1 > 0.5) | (p2 > 0.5):
@@ -297,7 +298,6 @@ def pcorrelate(t, u, bins):
 
 
 def fast_acf(counts, width, bin_width, cut_peak=True):
-
     """
     Super fast ACF function relying on numba (above).
 
@@ -339,7 +339,6 @@ def fast_acf(counts, width, bin_width, cut_peak=True):
 
 
 def acf_power(acf, norm=True):
-
     """
     Compute the power spectrum of the signal by computing the FFT of the autocorrelation.
 
@@ -384,9 +383,8 @@ def nonspatial_phase_precession(
     norm: bool = True,
     psd_lims: list = [0.65, 1.55],
     upsample: int = 4,
-    smooth_sigma=1
+    smooth_sigma=1,
 ):
-
     """
     Compute the nonspatial spike-LFP relationship modulation index.
 
@@ -425,7 +423,9 @@ def nonspatial_phase_precession(
     )
 
     frequencies = np.interp(
-        np.arange(0, len(frequencies), 1/upsample), np.arange(0, len(frequencies)), frequencies
+        np.arange(0, len(frequencies), 1 / upsample),
+        np.arange(0, len(frequencies)),
+        frequencies,
     )
 
     freqs_of_interest = np.intersect1d(
@@ -436,7 +436,7 @@ def nonspatial_phase_precession(
     psd = acf_power(acf, norm=norm)
 
     # upsample 2x psd
-    psd = np.interp(np.arange(0, len(psd), 1/upsample), np.arange(0, len(psd)), psd)
+    psd = np.interp(np.arange(0, len(psd), 1 / upsample), np.arange(0, len(psd)), psd)
     # smooth psd with gaussian filter
     psd = gaussian_filter1d(psd, smooth_sigma)
 
@@ -445,7 +445,13 @@ def nonspatial_phase_precession(
 
     # make sure there is a peak
     if ~np.any(all_peaks):
-        return np.nan, np.nan, psd[freqs_of_interest], frequencies[freqs_of_interest], acf
+        return (
+            np.nan,
+            np.nan,
+            psd[freqs_of_interest],
+            frequencies[freqs_of_interest],
+            acf,
+        )
 
     max_peak = np.max(psd[freqs_of_interest][all_peaks])
     max_idx = [all_peaks[np.argmax(psd[freqs_of_interest][all_peaks])]]
diff --git a/neuro_py/stats/circ_stats.py b/neuro_py/stats/circ_stats.py
new file mode 100644
index 0000000..fd86b4a
--- /dev/null
+++ b/neuro_py/stats/circ_stats.py
@@ -0,0 +1,548 @@
+import warnings
+from collections import namedtuple
+
+import numpy as np
+from decorator import decorator
+from scipy import stats
+
+# minimal version of pycircstat from https://github.com/circstat/pycircstat/tree/master
+
+CI = namedtuple("confidence_interval", ["lower", "upper"])
+
+
+def nd_bootstrap(data, iterations, axis=None, strip_tuple_if_one=True):
+    """
+    Bootstrap iterator for several n-dimensional data arrays.
+
+    :param data: Iterable containing the data arrays
+    :param iterations: Number of bootstrap iterations.
+    :param axis: Bootstrapping is performed along this axis.
+    """
+    shape0 = data[0].shape
+    if axis is None:
+        axis = 0
+        data = [d.ravel() for d in data]
+
+    n = len(data[0].shape)
+    K = len(data)
+    data0 = []
+
+    if axis is not None:
+        m = data[0].shape[axis]
+        to = tuple([axis]) + tuple(range(axis)) + tuple(range(axis + 1, n))
+        fro = tuple(range(1, axis + 1)) + (0,) + tuple(range(axis + 1, n))
+        for i in range(K):
+            data0.append(data[i].transpose(to))
+
+        for i in range(iterations):
+            idx = np.random.randint(m, size=(m,))
+            if len(data) == 1 and strip_tuple_if_one:
+                yield (
+                    data0[0][np.ix_(idx), ...].squeeze().transpose(fro).reshape(shape0)
+                )
+            else:
+                yield tuple(
+                    a[np.ix_(idx), ...].squeeze().transpose(fro).reshape(shape0)
+                    for a in data0
+                )
+
+
+def mod2pi(f):
+    """
+    Decorator to apply modulo 2*pi on the output of the function.
+
+    The decorated function must either return a tuple of numpy.ndarrays or a
+    numpy.ndarray itself.
+    """
+
+    def wrapper(f, *args, **kwargs):
+        ret = f(*args, **kwargs)
+
+        if isinstance(ret, tuple):
+            ret2 = []
+            for r in ret:
+                if isinstance(r, np.ndarray) or np.isscalar(r):
+                    ret2.append(r % (2 * np.pi))
+                elif isinstance(r, CI):
+                    ret2.append(CI(r.lower % (2 * np.pi), r.upper % (2 * np.pi)))
+                else:
+                    raise TypeError("Type not known!")
+            return tuple(ret2)
+        elif isinstance(ret, np.ndarray) or np.isscalar(ret):
+            return ret % (2 * np.pi)
+        else:
+            raise TypeError("Type not known!")
+
+    return decorator(wrapper, f)
+
+
+class bootstrap:
+    """
+    Decorator to implement bootstrapping. It looks for the arguments ci, axis,
+    and bootstrap_iter to determine the proper parameters for bootstrapping.
+    The argument scale determines whether the percentile is taken on a circular
+    scale or on a linear scale.
+
+    :param no_bootstrap: the number of arguments that are bootstrapped
+                        (e.g. for correlation it would be two, for median it
+                        would be one)
+    :param scale: linear or ciruclar scale (default is 'linear')
+    """
+
+    def __init__(self, no_bootstrap, scale="linear"):
+        self.no_boostrap = no_bootstrap
+        self.scale = scale
+
+    def _get_var(self, f, what, default, args, kwargs, remove=False):
+        varnames = f.__code__.co_varnames
+
+        if what in varnames:
+            what_idx = varnames.index(what)
+        else:
+            raise ValueError(
+                "Function %s does not have variable %s." % (f.__name__, what)
+            )
+
+        if len(args) >= what_idx + 1:
+            val = args[what_idx]
+            if remove:
+                args[what_idx] = default
+        else:
+            val = default
+
+        return val
+
+    def __call__(self, f):
+
+        def wrapper(f, *args, **kwargs):
+            args = list(args)
+            ci = self._get_var(f, "ci", None, args, kwargs, remove=True)
+            bootstrap_iter = self._get_var(
+                f, "bootstrap_iter", None, args, kwargs, remove=True
+            )
+            axis = self._get_var(f, "axis", None, args, kwargs)
+
+            alpha = args[: self.no_boostrap]
+            args0 = args[self.no_boostrap :]
+
+            if bootstrap_iter is None:
+                bootstrap_iter = (
+                    alpha[0].shape[axis] if axis is not None else alpha[0].size
+                )
+
+            r0 = f(*(alpha + args0), **kwargs)
+            if ci is not None:
+                r = np.asarray(
+                    [
+                        f(*(list(a) + args0), **kwargs)
+                        for a in nd_bootstrap(
+                            alpha, bootstrap_iter, axis=axis, strip_tuple_if_one=False
+                        )
+                    ]
+                )
+
+                if self.scale == "linear":
+                    ci_low, ci_high = np.percentile(
+                        r, [(1 - ci) / 2 * 100, (1 + ci) / 2 * 100], axis=0
+                    )
+                elif self.scale == "circular":
+                    ci_low, ci_high = percentile(
+                        r,
+                        [(1 - ci) / 2 * 100, (1 + ci) / 2 * 100],
+                        q0=(r0 + np.pi) % (2 * np.pi),
+                        axis=0,
+                    )
+                else:
+                    raise ValueError("Scale %s not known!" % (self.scale,))
+                return r0, CI(ci_low, ci_high)
+            else:
+                return r0
+
+        return decorator(wrapper, f)
+
+
+@mod2pi
+@bootstrap(1, "circular")
+def percentile(alpha, q, q0, axis=None, ci=None, bootstrap_iter=None):
+    """
+    Computes circular percentiles
+
+    :param alpha: array with circular samples
+    :param q: percentiles in [0,100] (single number or iterable)
+    :param q0: value of the 0 percentile
+    :param axis: percentiles will be computed along this axis.
+                 If None percentiles will be computed over the entire array
+    :param ci: if not None, confidence level is bootstrapped
+    :param bootstrap_iter: number of bootstrap iterations
+                           (number of samples if None)
+
+    :return: percentiles
+
+    """
+    if axis is None:
+        alpha = (alpha.ravel() - q0) % (2 * np.pi)
+    else:
+        if len(q0.shape) == len(alpha.shape) - 1:
+            reshaper = tuple(
+                slice(None, None) if i != axis else np.newaxis
+                for i in range(len(alpha.shape))
+            )
+            q0 = q0[reshaper]
+        elif not len(q0.shape) == len(alpha.shape):
+            raise ValueError("Dimensions of start and alpha are inconsistent!")
+
+        alpha = (alpha - q0) % (2 * np.pi)
+
+    ret = []
+    if axis is not None:
+        selector = tuple(
+            slice(None) if i != axis else 0 for i in range(len(alpha.shape))
+        )
+        q0 = q0[selector]
+
+    for qq in np.atleast_1d(q):
+        ret.append(np.percentile(alpha, qq, axis=axis) + q0)
+
+    if not hasattr(q, "__iter__"):  # if q is not some sort of list, array, etc
+        return np.asarray(ret).squeeze()
+    else:
+        return np.asarray(ret)
+
+
+def _complex_mean(alpha, w=None, axis=None, axial_correction=1):
+    if w is None:
+        w = np.ones_like(alpha)
+    alpha = np.asarray(alpha)
+
+    assert w.shape == alpha.shape, (
+        "Dimensions of data "
+        + str(alpha.shape)
+        + " and w "
+        + str(w.shape)
+        + " do not match!"
+    )
+
+    return (w * np.exp(1j * alpha * axial_correction)).sum(axis=axis) / np.sum(
+        w, axis=axis
+    )
+
+
+@bootstrap(1, "linear")
+def resultant_vector_length(
+    alpha, w=None, d=None, axis=None, axial_correction=1, ci=None, bootstrap_iter=None
+):
+    """
+    Computes mean resultant vector length for circular data.
+
+    This statistic is sometimes also called vector strength.
+
+    :param alpha: sample of angles in radians
+    :param w: number of incidences in case of binned angle data
+    :param ci: ci-confidence limits are computed via bootstrapping,
+               default None.
+    :param d: spacing of bin centers for binned data, if supplied
+              correction factor is used to correct for bias in
+              estimation of r, in radians (!)
+    :param axis: compute along this dimension, default is None
+                 (across all dimensions)
+    :param axial_correction: axial correction (2,3,4,...), default is 1
+    :param bootstrap_iter: number of bootstrap iterations
+                          (number of samples if None)
+    :return: mean resultant length
+
+    References: [Fisher1995]_, [Jammalamadaka2001]_, [Zar2009]_
+    """
+    if axis is None:
+        axis = 0
+        alpha = alpha.ravel()
+        if w is not None:
+            w = w.ravel()
+
+    cmean = _complex_mean(alpha, w=w, axis=axis, axial_correction=axial_correction)
+
+    # obtain length
+    r = np.abs(cmean)
+
+    # for data with known spacing, apply correction factor to correct for bias
+    # in the estimation of r (see Zar, p. 601, equ. 26.16)
+    if d is not None:
+        if axial_correction > 1:
+            warnings.warn("Axial correction ignored for bias correction.")
+        r *= d / 2 / np.sin(d / 2)
+    return r
+
+
+# defines synonym for resultant_vector_length
+vector_strength = resultant_vector_length
+
+
+def mean_ci_limits(alpha, ci=0.95, w=None, d=None, axis=None):
+    """
+    Computes the confidence limits on the mean for circular data.
+
+    :param alpha: sample of angles in radians
+    :param ci: ci-confidence limits are computed, default 0.95
+    :param w: number of incidences in case of binned angle data
+    :param d: spacing of bin centers for binned data, if supplied
+              correction factor is used to correct for bias in
+              estimation of r, in radians (!)
+    :param axis: compute along this dimension, default is None
+                 (across all dimensions)
+
+    :return: confidence limit width d; mean +- d yields upper/lower
+             (1-xi)% confidence limit
+
+    References: [Fisher1995]_, [Jammalamadaka2001]_, [Zar2009]_
+    """
+
+    if w is None:
+        w = np.ones_like(alpha)
+
+    assert alpha.shape == w.shape, "Dimensions of data and w do not match!"
+
+    r = np.atleast_1d(resultant_vector_length(alpha, w=w, d=d, axis=axis))
+    n = np.atleast_1d(np.sum(w, axis=axis))
+
+    R = n * r
+    c2 = stats.chi2.ppf(ci, df=1)
+
+    t = np.NaN * np.empty_like(r)
+
+    idx = (r < 0.9) & (r > np.sqrt(c2 / 2 / n))
+    t[idx] = np.sqrt(
+        (2 * n[idx] * (2 * R[idx] ** 2 - n[idx] * c2)) / (4 * n[idx] - c2)
+    )  # eq. 26.24
+
+    idx2 = r >= 0.9
+    t[idx2] = np.sqrt(
+        n[idx2] ** 2 - (n[idx2] ** 2 - R[idx2] ** 2) * np.exp(c2 / n[idx2])
+    )  # equ. 26.25
+
+    if not np.all(idx | idx2):
+        raise UserWarning(
+            """Requirements for confidence levels not met:
+                CI limits require a certain concentration of the data around the mean"""
+        )
+
+    return np.squeeze(np.arccos(t / R))
+
+
+@mod2pi
+def mean(alpha, w=None, ci=None, d=None, axis=None, axial_correction=1):
+    """
+    Compute mean direction of circular data.
+
+    :param alpha: circular data
+    :param w: 	 weightings in case of binned angle data
+    :param ci: if not None, the upper and lower 100*ci% confidence
+               interval is returned as well
+    :param d: spacing of bin centers for binned data, if supplied
+              correction factor is used to correct for bias in
+              estimation of r, in radians (!)
+    :param axis: compute along this dimension, default is None
+                 (across all dimensions)
+    :param axial_correction: axial correction (2,3,4,...), default is 1
+    :return: circular mean if ci=None, or circular mean as well as lower and
+             upper confidence interval limits
+
+    Example
+    >>> import numpy as np
+    >>> data = 2*np.pi*np.random.rand(10)
+    >>> mu, (ci_l, ci_u) = mean(data, ci=0.95)
+
+    """
+
+    cmean = _complex_mean(alpha, w=w, axis=axis, axial_correction=axial_correction)
+
+    mu = np.angle(cmean) / axial_correction
+
+    if ci is None:
+        return mu
+    else:
+        if axial_correction > 1:  # TODO: implement CI for axial correction
+            warnings.warn("Axial correction ignored for confidence intervals.")
+        t = mean_ci_limits(alpha, ci=ci, w=w, d=d, axis=axis)
+        return mu, CI(mu - t, mu + t)
+
+
+@mod2pi
+def center(*args, **kwargs):
+    """
+    Centers the data on its circular mean.
+
+    Each non-keyword argument is another data array that is centered.
+
+    :param axis: the mean is computed along this dimension (default axis=None).
+                **Must be used as a keyword argument!**
+    :return: tuple of centered data arrays
+
+    """
+
+    axis = kwargs.pop("axis", None)
+    if axis is None:
+        axis = 0
+        args = [a.ravel() for a in args]
+
+    reshaper = tuple(
+        slice(None, None) if i != axis else np.newaxis
+        for i in range(len(args[0].shape))
+    )
+    if len(args) == 1:
+        return args[0] - mean(args[0], axis=axis)
+    else:
+        return tuple(
+            [
+                a - mean(a, axis=axis)[reshaper]
+                for a in args
+                if isinstance(a, np.ndarray)
+            ]
+        )
+
+
+def get_var(f, varnames, args, kwargs):
+    fvarnames = f.__code__.co_varnames
+
+    var_idx = []
+    kwar_keys = []
+    for varname in varnames:
+        if varname in fvarnames:
+            var_pos = fvarnames.index(varname)
+        else:
+            raise ValueError(
+                "Function %s does not have variable %s." % (f.__name__, varnames)
+            )
+        if len(args) >= var_pos + 1:
+            var_idx.append(var_pos)
+        elif varname in kwargs:
+            kwar_keys.append(varname)
+        else:
+            raise ValueError("%s was not specified in  %s." % (varnames, f.__name__))
+
+    return var_idx, kwar_keys
+
+
+class swap2zeroaxis:
+    """
+    This decorator is best explained by an example::
+
+        @swap2zeroaxis(['x','y'], [0, 1])
+        def dummy(x,y,z, axis=None):
+            return np.mean(x[::2,...], axis=0), np.mean(y[::2, ...], axis=0), z
+
+    This creates a new function that
+
+    - either swaps the axes axis to zero for the arguments x and y if axis
+      is specified in dummy or ravels x and y
+    - swaps back the axes from the output arguments 0 and 1. Here it is
+      assumed that the outputs lost one dimension during the function
+      (e.g. like numpy.mean(x, axis=1) looses one axis).
+    """
+
+    def __init__(self, inputs, out_idx):
+        self.inputs = inputs
+        self.out_idx = out_idx
+
+    def __call__(self, f):
+
+        def _deco(f, *args, **kwargs):
+
+            to_swap_idx, to_swap_keys = get_var(f, self.inputs, args, kwargs)
+            args = list(args)
+
+            # extract axis parameter
+            try:
+                axis_idx, axis_kw = get_var(f, ["axis"], args, kwargs)
+                if len(axis_idx) == 0 and len(axis_kw) == 0:
+                    axis = None
+                else:
+                    if len(axis_idx) > 0:
+                        axis, args[axis_idx[0]] = args[axis_idx[0]], 0
+                    else:
+                        axis, kwargs[axis_kw[0]] = kwargs[axis_kw[0]], 0
+            except ValueError:
+                axis = None
+
+            # adjust axes or flatten
+            if axis is not None:
+                for i in to_swap_idx:
+                    if args[i] is not None:
+                        args[i] = args[i].swapaxes(0, axis)
+                for k in to_swap_keys:
+                    if kwargs[k] is not None:
+                        kwargs[k] = kwargs[k].swapaxes(0, axis)
+            else:
+                for i in to_swap_idx:
+                    if args[i] is not None:
+                        args[i] = args[i].ravel()
+                for k in to_swap_keys:
+                    if kwargs[k] is not None:
+                        kwargs[k] = kwargs[k].ravel()
+
+            # compute function
+            outputs = f(*args, **kwargs)
+
+            # swap everything back into place
+            if len(self.out_idx) > 0 and axis is not None:
+                if isinstance(outputs, tuple):
+                    outputs = list(outputs)
+                    for i in self.out_idx:
+                        outputs[i] = (
+                            outputs[i][np.newaxis, ...].swapaxes(0, axis).squeeze()
+                        )
+
+                    return tuple(outputs)
+                else:
+                    if self.out_idx != [0]:
+                        raise ValueError(
+                            "Single output argument and out_idx \
+                                         != [0] are inconsistent!"
+                        )
+                    return outputs[np.newaxis, ...].swapaxes(0, axis).squeeze()
+            else:
+                return outputs
+
+        return decorator(_deco, f)
+
+
+@swap2zeroaxis(["alpha", "w"], [0, 1])
+def rayleigh(alpha, w=None, d=None, axis=None):
+    """
+    Computes Rayleigh test for non-uniformity of circular data.
+
+    H0: the population is uniformly distributed around the circle
+    HA: the populatoin is not distributed uniformly around the circle
+
+    Assumption: the distribution has maximally one mode and the data is
+    sampled from a von Mises distribution!
+
+    :param alpha: sample of angles in radian
+    :param w:       number of incidences in case of binned angle data
+    :param d:     spacing of bin centers for binned data, if supplied
+                  correction factor is used to correct for bias in
+                  estimation of r
+    :param axis:  compute along this dimension, default is None
+                  if axis=None, array is raveled
+    :return pval: two-tailed p-value
+    :return z:    value of the z-statistic
+
+    References: [Fisher1995]_, [Jammalamadaka2001]_, [Zar2009]_
+    """
+
+    if w is None:
+        w = np.ones_like(alpha)
+
+    assert w.shape == alpha.shape, "Dimensions of alpha and w must match"
+
+    r = resultant_vector_length(alpha, w=w, d=d, axis=axis)
+    n = np.sum(w, axis=axis)
+
+    # compute Rayleigh's R (equ. 27.1)
+    R = n * r
+
+    # compute Rayleigh's z (equ. 27.2)
+    z = R**2 / n
+
+    # compute p value using approxation in Zar, p. 617
+    pval = np.exp(np.sqrt(1 + 4 * n + 4 * (n**2 - R**2)) - (1 + 2 * n))
+
+    return pval, z
diff --git a/pyproject.toml b/pyproject.toml
index c8f6149..972d1bf 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -39,7 +39,6 @@ dependencies = [
     "neo>=0.13.3",
     "quantities>=0.15.0",
     "neurodsp>=2.2.1",
-    "pycircstat>=0.0.2",
     "track-linearization>=2.3.1"
 ]
 

From 8cca2d62cd516d0d9b03698acb555b0504d4f02f Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 13:26:51 -0400
Subject: [PATCH 03/16] add pyFFTW

---
 pyproject.toml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 972d1bf..99b38c6 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -39,7 +39,8 @@ dependencies = [
     "neo>=0.13.3",
     "quantities>=0.15.0",
     "neurodsp>=2.2.1",
-    "track-linearization>=2.3.1"
+    "track-linearization>=2.3.1",
+    "pyFFTW>=0.14.0"
 ]
 
 [project.urls]

From 6ef79906c69be2b7c7ce94692561642b3caf0732 Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 13:29:39 -0400
Subject: [PATCH 04/16] Update pyproject.toml

---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 99b38c6..7940087 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -40,7 +40,7 @@ dependencies = [
     "quantities>=0.15.0",
     "neurodsp>=2.2.1",
     "track-linearization>=2.3.1",
-    "pyFFTW>=0.14.0"
+    "pyFFTW>=0.14"
 ]
 
 [project.urls]

From 15bbc99ca17c22e7ac0cd915149769cd3bced978 Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 13:31:10 -0400
Subject: [PATCH 05/16] pyFFTW ver

---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 7940087..964a429 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -40,7 +40,7 @@ dependencies = [
     "quantities>=0.15.0",
     "neurodsp>=2.2.1",
     "track-linearization>=2.3.1",
-    "pyFFTW>=0.14"
+    "pyFFTW>=0.13.1"
 ]
 
 [project.urls]

From 2c93bce272b1b3c3b769b8c544b45082be18db15 Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 13:33:57 -0400
Subject: [PATCH 06/16] add statsmodels

---
 pyproject.toml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 964a429..2ffed97 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -40,7 +40,8 @@ dependencies = [
     "quantities>=0.15.0",
     "neurodsp>=2.2.1",
     "track-linearization>=2.3.1",
-    "pyFFTW>=0.13.1"
+    "pyFFTW>=0.13.1",
+    "statsmodels>=0.14.2"
 ]
 
 [project.urls]

From 0ecc699e7e64b2799b0bda438fb65615e526afa9 Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 13:35:02 -0400
Subject: [PATCH 07/16] statsmodels ver

---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 2ffed97..dcea8df 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -41,7 +41,7 @@ dependencies = [
     "neurodsp>=2.2.1",
     "track-linearization>=2.3.1",
     "pyFFTW>=0.13.1",
-    "statsmodels>=0.14.2"
+    "statsmodels>=0.14.1"
 ]
 
 [project.urls]

From 794380486510a73fac3c412ad453581899acbbae Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 13:40:46 -0400
Subject: [PATCH 08/16] tuple to Tuple

---
 neuro_py/lfp/theta_cycles.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/neuro_py/lfp/theta_cycles.py b/neuro_py/lfp/theta_cycles.py
index 260d88a..c76a788 100644
--- a/neuro_py/lfp/theta_cycles.py
+++ b/neuro_py/lfp/theta_cycles.py
@@ -1,6 +1,6 @@
 import os
 import sys
-from typing import Union
+from typing import Union, Tuple
 
 import nelpy as nel
 import numpy as np
@@ -154,11 +154,12 @@ def save_theta_cycles(
 
 def get_theta_cycles(
     basepath: str,
-    theta_freq: tuple[int] = (6, 10),
+    theta_freq: Tuple[int] = (6, 10),
     lowpass: int = 48,
     detection_params: Union[dict, None] = None,
     ch: Union[int, None] = None,
 ):
+    # import bycycle, hidden import to avoid mandatory dependency
     from bycycle import Bycycle
 
     # load lfp as memmap

From 5ae7ecc99ca8ae3d3411b8f62be0c22cb7e604d4 Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 13:57:17 -0400
Subject: [PATCH 09/16] Update test_batch_analysis.py

---
 tests/test_batch_analysis.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/test_batch_analysis.py b/tests/test_batch_analysis.py
index 4219952..eabbbf7 100644
--- a/tests/test_batch_analysis.py
+++ b/tests/test_batch_analysis.py
@@ -92,9 +92,9 @@ def test_analysis(basepath):
 
         df = batch_analysis.load_results(save_path)
         assert df.shape[0] == 3
-        assert df["basepath"].iloc[0] == r"\test_data\test_data_1"
-        assert df["basepath"].iloc[1] == r"\test_data\test_data_2"
-        assert df["basepath"].iloc[2] == r"\test_data\test_data_3"
+        assert r"\test_data\test_data_1" in df["basepath"].values
+        assert r"\test_data\test_data_2" in df["basepath"].values
+        assert r"\test_data\test_data_3" in df["basepath"].values
 
     # test file encode/decode
     with tempfile.TemporaryDirectory() as save_path:

From bc72d01d87ec10556d0a22bc0e2424a521577d46 Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 14:03:43 -0400
Subject: [PATCH 10/16] add normpath

---
 tests/test_batch_analysis.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/tests/test_batch_analysis.py b/tests/test_batch_analysis.py
index eabbbf7..8dcd690 100644
--- a/tests/test_batch_analysis.py
+++ b/tests/test_batch_analysis.py
@@ -1,9 +1,11 @@
 import glob
 import os
 import pickle
+import tempfile
+
 import pandas as pd
+
 from neuro_py.process import batch_analysis
-import tempfile
 
 
 def test_batchanalysis():
@@ -99,7 +101,8 @@ def test_analysis(basepath):
     # test file encode/decode
     with tempfile.TemporaryDirectory() as save_path:
         file = r"C:\test_data\test_data_1"
+        file = os.path.normpath(file)
         encoded_file = batch_analysis.encode_file_path(file, save_path)
         decoded_file = batch_analysis.decode_file_path(encoded_file)
         assert decoded_file == file
-        assert encoded_file == save_path + os.sep + "C---___test_data___test_data_1.pkl"
+        assert encoded_file == os.path.normpath(save_path + os.sep + "C---___test_data___test_data_1.pkl")

From ed3de008137567574470b441febdd9e81f42ad51 Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 15:20:04 -0400
Subject: [PATCH 11/16] adjust filesep in test

---
 tests/test_batch_analysis.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/tests/test_batch_analysis.py b/tests/test_batch_analysis.py
index 8dcd690..c09f2a5 100644
--- a/tests/test_batch_analysis.py
+++ b/tests/test_batch_analysis.py
@@ -22,11 +22,10 @@ def test_analysis(basepath, a="Real", b="Python", c="Is", d="Great", e="!"):
         return results
 
     df = pd.DataFrame()
-    df["basepath"] = [
-        r"\test_data\test_data_1",
-        r"\test_data\test_data_2",
-        r"\test_data\test_data_3",
+    basepaths = [
+        os.sep + os.path.join("test_data", f"test_data_{i}") for i in range(1, 4)
     ]
+    df["basepath"] = basepaths
 
     # test serial
     with tempfile.TemporaryDirectory() as save_path:
@@ -94,15 +93,16 @@ def test_analysis(basepath):
 
         df = batch_analysis.load_results(save_path)
         assert df.shape[0] == 3
-        assert r"\test_data\test_data_1" in df["basepath"].values
-        assert r"\test_data\test_data_2" in df["basepath"].values
-        assert r"\test_data\test_data_3" in df["basepath"].values
+        for basepath in basepaths:
+            assert basepath in df["basepath"].values
 
     # test file encode/decode
     with tempfile.TemporaryDirectory() as save_path:
-        file = r"C:\test_data\test_data_1"
+        file = os.path.join("C:", os.sep, "test_data", "test_data_1")
         file = os.path.normpath(file)
         encoded_file = batch_analysis.encode_file_path(file, save_path)
         decoded_file = batch_analysis.decode_file_path(encoded_file)
         assert decoded_file == file
-        assert encoded_file == os.path.normpath(save_path + os.sep + "C---___test_data___test_data_1.pkl")
+        assert encoded_file == os.path.normpath(
+            os.path.join(save_path, "C---___test_data___test_data_1.pkl")
+        )

From 99c0ba233ade42866601baf43ca23acc6030390d Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 15:29:23 -0400
Subject: [PATCH 12/16] add drive letter

---
 tests/test_batch_analysis.py | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/tests/test_batch_analysis.py b/tests/test_batch_analysis.py
index c09f2a5..5767a2a 100644
--- a/tests/test_batch_analysis.py
+++ b/tests/test_batch_analysis.py
@@ -98,11 +98,18 @@ def test_analysis(basepath):
 
     # test file encode/decode
     with tempfile.TemporaryDirectory() as save_path:
-        file = os.path.join("C:", os.sep, "test_data", "test_data_1")
+        # extract drive letter
+        cwd = os.getcwd()
+        drive = os.path.splitdrive(cwd)[0]
+        # create a test file path
+        file = os.path.join(drive, os.sep, "test_data", "test_data_1")
         file = os.path.normpath(file)
+        # encode and decode the file path
         encoded_file = batch_analysis.encode_file_path(file, save_path)
         decoded_file = batch_analysis.decode_file_path(encoded_file)
+        # check that the decoded file path is the same as the original file path
         assert decoded_file == file
+        # check that the encoded file path is the same as the expected file path
         assert encoded_file == os.path.normpath(
-            os.path.join(save_path, "C---___test_data___test_data_1.pkl")
+            os.path.join(save_path, drive[0] + "---___test_data___test_data_1.pkl")
         )

From fd237baf515f0f99068fc1777054a75526c83297 Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 15:36:38 -0400
Subject: [PATCH 13/16] adj for linux systems

---
 tests/test_batch_analysis.py | 15 ++++++++++++---
 1 file changed, 12 insertions(+), 3 deletions(-)

diff --git a/tests/test_batch_analysis.py b/tests/test_batch_analysis.py
index 5767a2a..d9c6263 100644
--- a/tests/test_batch_analysis.py
+++ b/tests/test_batch_analysis.py
@@ -101,6 +101,12 @@ def test_analysis(basepath):
         # extract drive letter
         cwd = os.getcwd()
         drive = os.path.splitdrive(cwd)[0]
+
+        if drive:
+            drive_letter = drive[0]
+        else:
+            drive_letter = ''  # No drive letter on Linux systems
+        
         # create a test file path
         file = os.path.join(drive, os.sep, "test_data", "test_data_1")
         file = os.path.normpath(file)
@@ -110,6 +116,9 @@ def test_analysis(basepath):
         # check that the decoded file path is the same as the original file path
         assert decoded_file == file
         # check that the encoded file path is the same as the expected file path
-        assert encoded_file == os.path.normpath(
-            os.path.join(save_path, drive[0] + "---___test_data___test_data_1.pkl")
-        )
+        if drive_letter:
+            expected_encoded_file = os.path.join(save_path, drive_letter + "---___test_data___test_data_1.pkl")
+        else:
+            expected_encoded_file = os.path.join(save_path, "test_data___test_data_1.pkl")
+            
+        assert encoded_file == os.path.normpath(expected_encoded_file)
\ No newline at end of file

From 9ecad19743187a7404fb6c0c29b52bcc1500b4fe Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 16:27:52 -0400
Subject: [PATCH 14/16] Update test_batch_analysis.py

---
 tests/test_batch_analysis.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/test_batch_analysis.py b/tests/test_batch_analysis.py
index d9c6263..2fa5583 100644
--- a/tests/test_batch_analysis.py
+++ b/tests/test_batch_analysis.py
@@ -119,6 +119,6 @@ def test_analysis(basepath):
         if drive_letter:
             expected_encoded_file = os.path.join(save_path, drive_letter + "---___test_data___test_data_1.pkl")
         else:
-            expected_encoded_file = os.path.join(save_path, "test_data___test_data_1.pkl")
-            
+            expected_encoded_file = os.path.join(save_path, "___test_data___test_data_1.pkl")
+
         assert encoded_file == os.path.normpath(expected_encoded_file)
\ No newline at end of file

From ad035d0a8ce6c7e7b128784452745a30d69439e2 Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 16:46:47 -0400
Subject: [PATCH 15/16] fix: python 3.10 doesn't handle ndarray ==  []

---
 neuro_py/ensemble/assembly_reactivation.py | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/neuro_py/ensemble/assembly_reactivation.py b/neuro_py/ensemble/assembly_reactivation.py
index 81eb23c..15059bd 100644
--- a/neuro_py/ensemble/assembly_reactivation.py
+++ b/neuro_py/ensemble/assembly_reactivation.py
@@ -226,7 +226,7 @@ def get_weights(self, epoch=None):
 
         # check if st has any neurons
         if self.st.isempty:
-            self.patterns = []
+            self.patterns = None
             return
 
         if epoch is not None:
@@ -235,7 +235,7 @@ def get_weights(self, epoch=None):
             bst = self.st.bin(ds=self.weight_dt).data
 
         if (bst == 0).all():
-            self.patterns = []
+            self.patterns = None
         else:
             patterns, _, _ = assembly.runPatterns(
                 bst,
@@ -293,7 +293,7 @@ def plot(
         if not hasattr(self, "patterns"):
             return "run get_weights first"
         else:
-            if self.patterns == []:
+            if self.patterns is None:
                 return None, None
             if plot_members:
                 self.find_members()
@@ -352,9 +352,7 @@ def plot(
 
     def n_assemblies(self):
         if hasattr(self, "patterns"):
-            if self.patterns == []:
-                return 0
-            elif self.patterns is None:
+            if self.patterns is None:
                 return 0
             return self.patterns.shape[0]
 

From d97f688c5f07d28bb6e2232595b71be4b644d417 Mon Sep 17 00:00:00 2001
From: Ryan Harvey <ryanh412@gmail.com>
Date: Wed, 11 Sep 2024 16:55:04 -0400
Subject: [PATCH 16/16] linting

---
 neuro_py/ensemble/assembly.py              | 12 +++++-------
 neuro_py/ensemble/assembly_reactivation.py |  4 +---
 2 files changed, 6 insertions(+), 10 deletions(-)

diff --git a/neuro_py/ensemble/assembly.py b/neuro_py/ensemble/assembly.py
index c0115ee..fd60ad4 100644
--- a/neuro_py/ensemble/assembly.py
+++ b/neuro_py/ensemble/assembly.py
@@ -4,16 +4,14 @@
 	This implementation was written in Feb 2019.
 	Please e-mail me if you have comments, doubts, bug reports or criticism (Vítor, vtlsantos@gmail.com /  vitor.lopesdossantos@pharm.ox.ac.uk).
 """
-import warnings
-
-import numpy as np
 
+import warnings
 from typing import Tuple, Union
 
-from scipy import stats
-from sklearn.decomposition import FastICA
-from sklearn.decomposition import PCA
+import numpy as np
 from lazy_loader import attach as _attach
+from scipy import stats
+from sklearn.decomposition import PCA, FastICA
 
 __all__ = (
     "toyExample",
@@ -227,7 +225,7 @@ def runPatterns(
 
     if np.isnan(significance.nassemblies):
         return None, significance, None
-    
+
     if significance.nassemblies < 1:
         warnings.warn("no assembly detected")
 
diff --git a/neuro_py/ensemble/assembly_reactivation.py b/neuro_py/ensemble/assembly_reactivation.py
index 15059bd..d80d60d 100644
--- a/neuro_py/ensemble/assembly_reactivation.py
+++ b/neuro_py/ensemble/assembly_reactivation.py
@@ -13,9 +13,7 @@
 from neuro_py.io import loading
 from neuro_py.session.locate_epochs import compress_repeated_epochs, find_pre_task_post
 
-__all__ = (
-    "AssemblyReact",
-)
+__all__ = ("AssemblyReact",)
 __getattr__, __dir__, __all__ = _attach(f"{__name__}", submodules=__all__)
 del _attach