From 8835464cdb0bbb480ee649902b481c7acb3d9e53 Mon Sep 17 00:00:00 2001
From: ryanhammonds <Ryan.hammonds@utexas.edu>
Date: Tue, 8 Jun 2021 17:45:11 -0700
Subject: [PATCH 1/9] add pairwise phase consistency index

---
 neurodsp/timefrequency/consistency.py | 54 +++++++++++++++++++++++++++
 1 file changed, 54 insertions(+)
 create mode 100644 neurodsp/timefrequency/consistency.py

diff --git a/neurodsp/timefrequency/consistency.py b/neurodsp/timefrequency/consistency.py
new file mode 100644
index 00000000..b31dde7b
--- /dev/null
+++ b/neurodsp/timefrequency/consistency.py
@@ -0,0 +1,54 @@
+"""Phase consistency measures."""
+
+from itertools import combinations
+
+import numpy as np
+
+###################################################################################################
+###################################################################################################
+
+def pairwise_phase_consistency(phases):
+    """Compute pairwise phase consistency.
+
+    Parameters
+    ----------
+    phases : 2d array
+        Phases from a wavelet or Fourier analysis, in radians.
+
+    Returns
+    -------
+    avg_distance : float
+        Average pairwise circular distance index.
+    distances : 2d array
+        Pairwise circular distance indices.
+
+    Notes
+    -----
+
+    - distance == -1: inverse phases (i.e. np.pi vs -np.pi)
+    - distance ==  0: pi / 2 phase difference (i.e. np.pi vs np.pi / 2)
+    - distance ==  1: equal phases (i.e. np.pi vs np.pi)
+
+    Reference
+    ---------
+    Vinck, M., van Wingerden, M., Womelsdorf, T., Fries, P., & Pennartz, C. M. A. (2010).
+    The pairwise phase consistency: A bias-free measure of rhythmic neuronal synchronization.
+    NeuroImage, 51(1), 112–122. https://doi.org/10.1016/j.neuroimage.2010.01.073
+    """
+
+    pairs = list(combinations(np.arange(len(phases)), 2))
+
+    distances = np.zeros((len(pairs), len(phases[0])))
+
+    for idx, pair in enumerate(pairs):
+
+        # Absolute angular distance
+        distances[idx] = np.abs(phases[pair[0]] - phases[pair[1]]) % np.pi
+
+        # Pairwise circular distance index (PCDI)
+        distances[idx] = (np.pi - 2 * distances[idx]) / np.pi
+
+    # Mean PCDI
+    avg_distance = (2 * np.sum(distances)) / (len(phases) * (len(phases) - 1))
+
+    return avg_distance, distances

From 08cf88b9ba7427c4d102a0b520ab99eaab243aa8 Mon Sep 17 00:00:00 2001
From: ryanhammonds <ryan.hammonds@utexas.edu>
Date: Wed, 9 Jun 2021 17:38:54 -0700
Subject: [PATCH 2/9] pairwise phase consistency fixes

---
 neurodsp/timefrequency/consistency.py | 84 +++++++++++++++++++++------
 1 file changed, 67 insertions(+), 17 deletions(-)

diff --git a/neurodsp/timefrequency/consistency.py b/neurodsp/timefrequency/consistency.py
index b31dde7b..27ef5c56 100644
--- a/neurodsp/timefrequency/consistency.py
+++ b/neurodsp/timefrequency/consistency.py
@@ -1,5 +1,7 @@
 """Phase consistency measures."""
 
+import warnings
+from importlib import import_module
 from itertools import combinations
 
 import numpy as np
@@ -7,27 +9,39 @@
 ###################################################################################################
 ###################################################################################################
 
-def pairwise_phase_consistency(phases):
+def pairwise_phase_consistency(pha0, pha1, return_pairs=True, memory_gb=2, progress=None):
     """Compute pairwise phase consistency.
 
     Parameters
     ----------
-    phases : 2d array
-        Phases from a wavelet or Fourier analysis, in radians.
+    pha0 : 1d array
+        First phases from a wavelet analysis, in radians (i.e. lfp).
+    pha1 : 1d array
+        Second phases from a wavelet analysis, in radians (i.e. spikes).
+    return_pairs : True
+        Returns distance pairs as a 1d array if True.
+    memory_gb : float, optional, default: 2
+        Maximum size of the ``distances`` array, in gb. If the pairwise array is larger than this
+        parameter, distances will not be stored in memory to prevent OOM error. Ignored if
+        ``return_pairs`` is False.
+    progress : {None, 'tqdm', 'tqdm.notebook}
+        Displays tqdm progress bar.
 
     Returns
     -------
-    avg_distance : float
+    distance_avg : float
         Average pairwise circular distance index.
-    distances : 2d array
-        Pairwise circular distance indices.
+    distances : 2d array, optional
+        Pairwise circular distance indices. Only returned if ``return_pairs` is True. If
+        ``memory_gb`` is less than required array size, None will be returned.
 
     Notes
     -----
 
-    - distance == -1: inverse phases (i.e. np.pi vs -np.pi)
-    - distance ==  0: pi / 2 phase difference (i.e. np.pi vs np.pi / 2)
-    - distance ==  1: equal phases (i.e. np.pi vs np.pi)
+    - distance == -1: inverse phases
+    - distance ==  0: pi / 2 phase difference
+    - distance ==  1: equal phases
+
 
     Reference
     ---------
@@ -36,19 +50,55 @@ def pairwise_phase_consistency(phases):
     NeuroImage, 51(1), 112–122. https://doi.org/10.1016/j.neuroimage.2010.01.073
     """
 
-    pairs = list(combinations(np.arange(len(phases)), 2))
+    if pha0.shape != pha1.shape or pha0.ndim != 1:
+        raise ValueError("Phase arrays must be the same 1d length.")
+
+    n_combs = int((len(pha0) * (len(pha0) - 1)) / 2)
+
+    # Pairwise distance array memory limit
+    gb_per_float = 8e-9
+    limit_mem = n_combs * gb_per_float > memory_gb
+
+    if limit_mem and return_pairs:
+        warnings.warn("Memory limit is smaller than required distance array size. "
+                      "Pairwise distances will be returned as None.")
 
-    distances = np.zeros((len(pairs), len(phases[0])))
+    # Initialize variables
+    if return_pairs and not limit_mem:
+        cumulative = None
+        distances = np.zeros(n_combs)
+    else:
+        cumulative = 0
+        distances = None
 
-    for idx, pair in enumerate(pairs):
+    iterable = enumerate(combinations(np.arange(len(pha0)), 2))
+
+    # Optional progress bar
+    if progress is not None:
+        try:
+            tqdm = import_module(progress)
+            iterable = tqdm.tqdm(iterable, total=n_combs, dynamic_ncols=True,
+                                 desc='Computing Pairwise Distances')
+        except ImportError:
+            pass
+
+    # Compute distance indices
+    for idx, pair in iterable:
 
         # Absolute angular distance
-        distances[idx] = np.abs(phases[pair[0]] - phases[pair[1]]) % np.pi
+        abs_dist = np.abs(pha0[pair[0]] - pha1[pair[1]]) % np.pi
 
         # Pairwise circular distance index (PCDI)
-        distances[idx] = (np.pi - 2 * distances[idx]) / np.pi
+        distance = (np.pi - 2 * abs_dist) / np.pi
+
+        if isinstance(distances, np.ndarray):
+            distances[idx] = distance
+        else:
+            cumulative += distance
 
-    # Mean PCDI
-    avg_distance = (2 * np.sum(distances)) / (len(phases) * (len(phases) - 1))
+    distance_avg = cumulative.sum() / n_combs if distances is None else np.mean(distances)
 
-    return avg_distance, distances
+    if return_pairs:
+        return distance_avg, distances
+    else:
+        return distance_avg

From 9d6c16cb4badf447c649d6c51526500b7de7acb0 Mon Sep 17 00:00:00 2001
From: ryanhammonds <Ryan.hammonds@utexas.edu>
Date: Thu, 10 Jun 2021 13:13:07 -0700
Subject: [PATCH 3/9] move ppc to rhythm

---
 neurodsp/rhythm/__init__.py                                | 1 +
 neurodsp/{timefrequency/consistency.py => rhythm/phase.py} | 0
 2 files changed, 1 insertion(+)
 rename neurodsp/{timefrequency/consistency.py => rhythm/phase.py} (100%)

diff --git a/neurodsp/rhythm/__init__.py b/neurodsp/rhythm/__init__.py
index 5d8872d4..342c4ada 100644
--- a/neurodsp/rhythm/__init__.py
+++ b/neurodsp/rhythm/__init__.py
@@ -2,3 +2,4 @@
 
 from .lc import compute_lagged_coherence
 from .swm import sliding_window_matching
+from .phase import pairwise_phase_consistency
diff --git a/neurodsp/timefrequency/consistency.py b/neurodsp/rhythm/phase.py
similarity index 100%
rename from neurodsp/timefrequency/consistency.py
rename to neurodsp/rhythm/phase.py

From cd592d72b56dd18864031099f422c1d563608427 Mon Sep 17 00:00:00 2001
From: ryanhammonds <Ryan.hammonds@utexas.edu>
Date: Thu, 10 Jun 2021 14:54:34 -0700
Subject: [PATCH 4/9] modulo 2pi

---
 neurodsp/rhythm/phase.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/neurodsp/rhythm/phase.py b/neurodsp/rhythm/phase.py
index 27ef5c56..8c1c6074 100644
--- a/neurodsp/rhythm/phase.py
+++ b/neurodsp/rhythm/phase.py
@@ -86,7 +86,7 @@ def pairwise_phase_consistency(pha0, pha1, return_pairs=True, memory_gb=2, progr
     for idx, pair in iterable:
 
         # Absolute angular distance
-        abs_dist = np.abs(pha0[pair[0]] - pha1[pair[1]]) % np.pi
+        abs_dist = abs(abs(pha0[pair[0]]) - abs(pha1[pair[1]])) % (2 * np.pi)
 
         # Pairwise circular distance index (PCDI)
         distance = (np.pi - 2 * abs_dist) / np.pi

From 5b177b3ac031fe33bdaca0d2e8b0020267c10435 Mon Sep 17 00:00:00 2001
From: ryanhammonds <Ryan.hammonds@utexas.edu>
Date: Thu, 10 Jun 2021 16:25:16 -0700
Subject: [PATCH 5/9] ppc tests added

---
 neurodsp/tests/rhythm/test_phase.py | 61 +++++++++++++++++++++++++++++
 1 file changed, 61 insertions(+)
 create mode 100644 neurodsp/tests/rhythm/test_phase.py

diff --git a/neurodsp/tests/rhythm/test_phase.py b/neurodsp/tests/rhythm/test_phase.py
new file mode 100644
index 00000000..429dddbb
--- /dev/null
+++ b/neurodsp/tests/rhythm/test_phase.py
@@ -0,0 +1,61 @@
+"""Tests for neurodsp.rhythm.phase."""
+
+from pytest import mark
+
+import numpy as np
+
+from neurodsp.tests.settings import FS, FREQ_SINE
+from neurodsp.timefrequency import phase_by_time
+from neurodsp.rhythm.phase import *
+
+###################################################################################################
+###################################################################################################
+
+@mark.parametrize('return_pairs', [True, False])
+@mark.parametrize('low_mem', [True, False])
+@mark.parametrize('phase_shift', [0, .25, .5])
+def test_pairwise_phase_consistency(tsig_sine, return_pairs, low_mem, phase_shift):
+
+    memory_gb = -np.inf if low_mem else 2
+
+    peaks = np.where(tsig_sine == 1.0)[0]
+
+    pha0 = phase_by_time(tsig_sine, FS)
+
+    # Phase shift
+    sig_shift = np.roll(tsig_sine, int((FS / FREQ_SINE) * phase_shift))
+    pha1 = phase_by_time(sig_shift, FS)
+
+    # Case where arrays are different sizes
+    try:
+        pairwise_phase_consistency(pha0[peaks], pha1[peaks][:2],
+                                   return_pairs, memory_gb, 'tqdm')
+    except ValueError:
+        pass
+
+    # Compute consistency
+    dist_avg = pairwise_phase_consistency(pha0[peaks], pha1[peaks],
+                                          return_pairs, memory_gb, 'tqdm')
+
+    # Unpack results if needed
+    if return_pairs:
+
+        dist_avg, dists = dist_avg[0], dist_avg[1]
+
+        if low_mem:
+            assert dists is None
+        else:
+            assert isinstance(dists, np.ndarray)
+            assert len(dists) == (len(peaks) * (len(peaks) - 1)) / 2
+            assert np.mean(dists) == dist_avg
+
+    # Expected consistency
+    if phase_shift == 0:
+        expected = 1
+    elif phase_shift == .25:
+        expected = 0
+    elif phase_shift == .5:
+        expected = -1
+
+    assert isinstance(dist_avg, float)
+    assert round(dist_avg) == expected

From a69c4f8727e4be179b98e858c2019b1acef7b609 Mon Sep 17 00:00:00 2001
From: ryanhammonds <Ryan.hammonds@utexas.edu>
Date: Thu, 10 Jun 2021 17:59:54 -0700
Subject: [PATCH 6/9] correct phase angles

---
 neurodsp/rhythm/phase.py | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/neurodsp/rhythm/phase.py b/neurodsp/rhythm/phase.py
index 8c1c6074..51720f60 100644
--- a/neurodsp/rhythm/phase.py
+++ b/neurodsp/rhythm/phase.py
@@ -15,9 +15,9 @@ def pairwise_phase_consistency(pha0, pha1, return_pairs=True, memory_gb=2, progr
     Parameters
     ----------
     pha0 : 1d array
-        First phases from a wavelet analysis, in radians (i.e. lfp).
+        First phases from a wavelet analysis, in radians, from -pi to pi.
     pha1 : 1d array
-        Second phases from a wavelet analysis, in radians (i.e. spikes).
+        Second phases from a wavelet analysis, in radians, from -pi to pi.
     return_pairs : True
         Returns distance pairs as a 1d array if True.
     memory_gb : float, optional, default: 2
@@ -85,8 +85,18 @@ def pairwise_phase_consistency(pha0, pha1, return_pairs=True, memory_gb=2, progr
     # Compute distance indices
     for idx, pair in iterable:
 
+        phi0= pha0[pair[0]]
+        phi1 = pha1[pair[1]]
+
+        # Convert range from (-pi, pi) to (0, 2pi)
+        phi0 = phi0 + (2*np.pi) if phi0 < 0 else phi0
+        phi1 = phi1 + (2*np.pi) if phi1 < 0 else phi1
+
         # Absolute angular distance
-        abs_dist = abs(abs(pha0[pair[0]]) - abs(pha1[pair[1]])) % (2 * np.pi)
+        abs_dist = np.abs(phi0 - phi1)
+
+        # Take smaller angle (range 0 to pi)
+        abs_dist = (2*np.pi) - abs_dist if abs_dist > np.pi else abs_dist
 
         # Pairwise circular distance index (PCDI)
         distance = (np.pi - 2 * abs_dist) / np.pi

From 598d88510d9230d56598ca0097c3de6568036662 Mon Sep 17 00:00:00 2001
From: ryanhammonds <ryan.hammonds@utexas.edu>
Date: Fri, 11 Jun 2021 16:10:07 -0700
Subject: [PATCH 7/9] optional second phase array

---
 neurodsp/rhythm/phase.py            | 47 +++++++++++++++++------------
 neurodsp/tests/rhythm/test_phase.py | 43 ++++++++++++++------------
 2 files changed, 51 insertions(+), 39 deletions(-)

diff --git a/neurodsp/rhythm/phase.py b/neurodsp/rhythm/phase.py
index 51720f60..d9016b4b 100644
--- a/neurodsp/rhythm/phase.py
+++ b/neurodsp/rhythm/phase.py
@@ -2,28 +2,24 @@
 
 import warnings
 from importlib import import_module
-from itertools import combinations
+from itertools import combinations, combinations_with_replacement
 
 import numpy as np
 
 ###################################################################################################
 ###################################################################################################
 
-def pairwise_phase_consistency(pha0, pha1, return_pairs=True, memory_gb=2, progress=None):
+def pairwise_phase_consistency(pha0, pha1=None, return_pairs=True, progress=None):
     """Compute pairwise phase consistency.
 
     Parameters
     ----------
     pha0 : 1d array
         First phases from a wavelet analysis, in radians, from -pi to pi.
-    pha1 : 1d array
+    pha1 : 1d array, optional, default: None
         Second phases from a wavelet analysis, in radians, from -pi to pi.
     return_pairs : True
         Returns distance pairs as a 1d array if True.
-    memory_gb : float, optional, default: 2
-        Maximum size of the ``distances`` array, in gb. If the pairwise array is larger than this
-        parameter, distances will not be stored in memory to prevent OOM error. Ignored if
-        ``return_pairs`` is False.
     progress : {None, 'tqdm', 'tqdm.notebook}
         Displays tqdm progress bar.
 
@@ -50,29 +46,36 @@ def pairwise_phase_consistency(pha0, pha1, return_pairs=True, memory_gb=2, progr
     NeuroImage, 51(1), 112–122. https://doi.org/10.1016/j.neuroimage.2010.01.073
     """
 
-    if pha0.shape != pha1.shape or pha0.ndim != 1:
-        raise ValueError("Phase arrays must be the same 1d length.")
+    if pha0.ndim != 1:
+        raise ValueError("Phase array must be 1-dimensional.")
 
-    n_combs = int((len(pha0) * (len(pha0) - 1)) / 2)
+    if pha1 is not None and pha0.shape != pha1.shape:
+        raise ValueError("Phase arrays must be the same length.")
 
-    # Pairwise distance array memory limit
-    gb_per_float = 8e-9
-    limit_mem = n_combs * gb_per_float > memory_gb
+    # Pairwise indices generator
+    if pha1 is None:
 
-    if limit_mem and return_pairs:
-        warnings.warn("Memory limit is smaller than required distance array size. "
-                      "Pairwise distances will be returned as None.")
+        # Number of pairwise combinations
+        n_combs = int((len(pha0) * (len(pha0) - 1)) / 2)
+
+        # Exclude self-combinations (i.e. ignore (0, 0), (1, 1)...)
+        iterable = enumerate(combinations(np.arange(len(pha0)), 2))
+
+    else:
+
+        n_combs = int((len(pha0) * (len(pha0) + 1)) / 2)
+
+        # Include self-combinations
+        iterable = enumerate(combinations_with_replacement(np.arange(len(pha0)), 2))
 
     # Initialize variables
-    if return_pairs and not limit_mem:
+    if return_pairs:
         cumulative = None
         distances = np.zeros(n_combs)
     else:
         cumulative = 0
         distances = None
 
-    iterable = enumerate(combinations(np.arange(len(pha0)), 2))
-
     # Optional progress bar
     if progress is not None:
         try:
@@ -86,7 +89,11 @@ def pairwise_phase_consistency(pha0, pha1, return_pairs=True, memory_gb=2, progr
     for idx, pair in iterable:
 
         phi0= pha0[pair[0]]
-        phi1 = pha1[pair[1]]
+
+        if pha1 is None:
+            phi1 = pha0[pair[1]]
+        else:
+            phi1 = pha1[pair[1]]
 
         # Convert range from (-pi, pi) to (0, 2pi)
         phi0 = phi0 + (2*np.pi) if phi0 < 0 else phi0
diff --git a/neurodsp/tests/rhythm/test_phase.py b/neurodsp/tests/rhythm/test_phase.py
index 429dddbb..b947f528 100644
--- a/neurodsp/tests/rhythm/test_phase.py
+++ b/neurodsp/tests/rhythm/test_phase.py
@@ -12,11 +12,8 @@
 ###################################################################################################
 
 @mark.parametrize('return_pairs', [True, False])
-@mark.parametrize('low_mem', [True, False])
 @mark.parametrize('phase_shift', [0, .25, .5])
-def test_pairwise_phase_consistency(tsig_sine, return_pairs, low_mem, phase_shift):
-
-    memory_gb = -np.inf if low_mem else 2
+def test_pairwise_phase_consistency(tsig_sine, return_pairs, phase_shift):
 
     peaks = np.where(tsig_sine == 1.0)[0]
 
@@ -26,28 +23,17 @@ def test_pairwise_phase_consistency(tsig_sine, return_pairs, low_mem, phase_shif
     sig_shift = np.roll(tsig_sine, int((FS / FREQ_SINE) * phase_shift))
     pha1 = phase_by_time(sig_shift, FS)
 
-    # Case where arrays are different sizes
-    try:
-        pairwise_phase_consistency(pha0[peaks], pha1[peaks][:2],
-                                   return_pairs, memory_gb, 'tqdm')
-    except ValueError:
-        pass
-
     # Compute consistency
-    dist_avg = pairwise_phase_consistency(pha0[peaks], pha1[peaks],
-                                          return_pairs, memory_gb, 'tqdm')
+    dist_avg = pairwise_phase_consistency(pha0[peaks], pha1[peaks], return_pairs, 'tqdm')
 
     # Unpack results if needed
     if return_pairs:
 
         dist_avg, dists = dist_avg[0], dist_avg[1]
 
-        if low_mem:
-            assert dists is None
-        else:
-            assert isinstance(dists, np.ndarray)
-            assert len(dists) == (len(peaks) * (len(peaks) - 1)) / 2
-            assert np.mean(dists) == dist_avg
+        assert isinstance(dists, np.ndarray)
+        assert len(dists) == (len(peaks) * (len(peaks) + 1)) / 2
+        assert np.mean(dists) == dist_avg
 
     # Expected consistency
     if phase_shift == 0:
@@ -59,3 +45,22 @@ def test_pairwise_phase_consistency(tsig_sine, return_pairs, low_mem, phase_shif
 
     assert isinstance(dist_avg, float)
     assert round(dist_avg) == expected
+
+    # Test self-consistency
+    dist_avg, dists = pairwise_phase_consistency(pha0[peaks], return_pairs=True)
+
+    assert dist_avg == 1
+    assert len(dists) == (len(peaks) * (len(peaks) - 1)) / 2
+
+    # Cases where arrays are invalid sizes
+    try:
+        pairwise_phase_consistency(pha0[peaks], pha1[peaks][:2], return_pairs, 'tqdm')
+        assert False
+    except ValueError:
+        pass
+
+    try:
+        pairwise_phase_consistency(np.zeros((2, 2)), return_pairs, 'tqdm')
+        assert False
+    except ValueError:
+        pass

From de24e72d37cf95e2103ec12bd4df901a04a7267d Mon Sep 17 00:00:00 2001
From: ryanhammonds <ryan.hammonds@utexas.edu>
Date: Mon, 14 Jun 2021 14:34:50 -0700
Subject: [PATCH 8/9] two array combinations fix

---
 neurodsp/rhythm/phase.py            | 20 +++++++++++---------
 neurodsp/tests/rhythm/test_phase.py |  4 ++--
 2 files changed, 13 insertions(+), 11 deletions(-)

diff --git a/neurodsp/rhythm/phase.py b/neurodsp/rhythm/phase.py
index d9016b4b..bb9c1ce0 100644
--- a/neurodsp/rhythm/phase.py
+++ b/neurodsp/rhythm/phase.py
@@ -28,8 +28,7 @@ def pairwise_phase_consistency(pha0, pha1=None, return_pairs=True, progress=None
     distance_avg : float
         Average pairwise circular distance index.
     distances : 2d array, optional
-        Pairwise circular distance indices. Only returned if ``return_pairs` is True. If
-        ``memory_gb`` is less than required array size, None will be returned.
+        Pairwise circular distance indices. Only returned if ``return_pairs` is True.
 
     Notes
     -----
@@ -63,15 +62,15 @@ def pairwise_phase_consistency(pha0, pha1=None, return_pairs=True, progress=None
 
     else:
 
-        n_combs = int((len(pha0) * (len(pha0) + 1)) / 2)
+        # Include all combinations
+        n_combs = int(len(pha0) ** 2)
 
-        # Include self-combinations
-        iterable = enumerate(combinations_with_replacement(np.arange(len(pha0)), 2))
+        iterable = enumerate((row, col) for row in range(len(pha0)) for col in range(len(pha1)))
 
     # Initialize variables
     if return_pairs:
         cumulative = None
-        distances = np.zeros(n_combs)
+        distances = np.ones((len(pha0), len(pha0)))
     else:
         cumulative = 0
         distances = None
@@ -88,7 +87,7 @@ def pairwise_phase_consistency(pha0, pha1=None, return_pairs=True, progress=None
     # Compute distance indices
     for idx, pair in iterable:
 
-        phi0= pha0[pair[0]]
+        phi0 = pha0[pair[0]]
 
         if pha1 is None:
             phi1 = pha0[pair[1]]
@@ -108,8 +107,11 @@ def pairwise_phase_consistency(pha0, pha1=None, return_pairs=True, progress=None
         # Pairwise circular distance index (PCDI)
         distance = (np.pi - 2 * abs_dist) / np.pi
 
-        if isinstance(distances, np.ndarray):
-            distances[idx] = distance
+        if isinstance(distances, np.ndarray) and pha1 is None:
+            distances[pair[0], pair[1]] = distance
+            distances[pair[1], pair[0]] = distance
+        elif isinstance(distances, np.ndarray) and pha1 is not None:
+            distances[pair[0], pair[1]] = distance
         else:
             cumulative += distance
 
diff --git a/neurodsp/tests/rhythm/test_phase.py b/neurodsp/tests/rhythm/test_phase.py
index b947f528..5fe649c4 100644
--- a/neurodsp/tests/rhythm/test_phase.py
+++ b/neurodsp/tests/rhythm/test_phase.py
@@ -32,7 +32,7 @@ def test_pairwise_phase_consistency(tsig_sine, return_pairs, phase_shift):
         dist_avg, dists = dist_avg[0], dist_avg[1]
 
         assert isinstance(dists, np.ndarray)
-        assert len(dists) == (len(peaks) * (len(peaks) + 1)) / 2
+        assert len(dists[0]) * len(dists[1]) == len(peaks) ** 2
         assert np.mean(dists) == dist_avg
 
     # Expected consistency
@@ -50,7 +50,7 @@ def test_pairwise_phase_consistency(tsig_sine, return_pairs, phase_shift):
     dist_avg, dists = pairwise_phase_consistency(pha0[peaks], return_pairs=True)
 
     assert dist_avg == 1
-    assert len(dists) == (len(peaks) * (len(peaks) - 1)) / 2
+    assert len(dists[0]) == len(dists[1]) == len(peaks)
 
     # Cases where arrays are invalid sizes
     try:

From 310a4e4d59df3e202299cafaa309b5311f05eaa5 Mon Sep 17 00:00:00 2001
From: ryanhammonds <ryan.hammonds@utexas.edu>
Date: Wed, 21 Jul 2021 15:45:02 -0700
Subject: [PATCH 9/9] fix bug when return_pairs is false

---
 neurodsp/rhythm/phase.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/neurodsp/rhythm/phase.py b/neurodsp/rhythm/phase.py
index bb9c1ce0..5d05bce5 100644
--- a/neurodsp/rhythm/phase.py
+++ b/neurodsp/rhythm/phase.py
@@ -115,7 +115,7 @@ def pairwise_phase_consistency(pha0, pha1=None, return_pairs=True, progress=None
         else:
             cumulative += distance
 
-    distance_avg = cumulative.sum() / n_combs if distances is None else np.mean(distances)
+    distance_avg = cumulative / n_combs if distances is None else np.mean(distances)
 
     if return_pairs:
         return distance_avg, distances