From a4dc222ff602998fa7ceb1387e089321776e3a69 Mon Sep 17 00:00:00 2001 From: Vinay Raghavan <42253618+vinaysraghavan@users.noreply.github.com> Date: Fri, 10 Jan 2025 15:29:01 -0500 Subject: [PATCH 01/10] Update rereference.py Change the make_contact_rereference_arr function to handle input channelnames in any order (not strictly ascending) and to properly handle local referencing on ECoG grid --- naplib/preprocessing/rereference.py | 92 +++++++++++++++++++++++------ 1 file changed, 74 insertions(+), 18 deletions(-) diff --git a/naplib/preprocessing/rereference.py b/naplib/preprocessing/rereference.py index 97b5293d..431081d3 100644 --- a/naplib/preprocessing/rereference.py +++ b/naplib/preprocessing/rereference.py @@ -128,13 +128,12 @@ def make_contact_rereference_arr(channelnames, extent=None): be alphanumeric, with any numbers only being on the right. 2) The numeric portion specifies a different electrode number, while the character portion in the left of the channelname specifies the contact name. E.g. ['RT1','RT2','RT3','Ls1','Ls2'] indicates two contacts, the first with 3 electrodes - and the second with 2 electrodes. 3) Electrodes from the same contact must be contiguous. + and the second with 2 electrodes. extent : int, optional, default=None If provided, then only contacts from the same group which are within ``extent`` electrodes away - from each other (inclusive) are still grouped together. Only used if ``method='contact'``. For - example, if ``extent=1``, only the nearest electrode on either side of a given electrode on the - same contact is still grouped with it. For example, extent=1 produces the traditional local - average reference scheme. + from each other (inclusive) are still grouped together. For example, if ``extent=1``, only the + nearest electrode on either side of a given electrode on the same contact is still grouped with it. + This ``extent=1`` produces the traditional local average reference scheme. Returns ------- @@ -145,18 +144,75 @@ def make_contact_rereference_arr(channelnames, extent=None): -------- rereference """ - contact_arrays = pd.Series([x.rstrip('0123456789') for x in channelnames]) - connections = np.zeros((len(contact_arrays),) * 2, dtype=float) - for _, inds in contact_arrays.groupby(contact_arrays): - for i in inds.index: - connections[i, inds.index] = 1.0 + def _find_adjacent_numbers(a, b, number, extent): + ''' + Used to determine electrodes for local averaging ECoG grid" + ''' + # Validate if the number is within the valid range + if number < 1 or number > a * b: + raise ValueError("The number is outside the range of the grid.") - # remove longer than extent if desired - if extent is not None: - if extent < 1: - raise ValueError(f'Invalid extent. Must be no less than 1 but got extent={extent}') - connections *= np.tri(*connections.shape, k=extent) - connections *= np.fliplr(np.flipud(np.tri(*connections.shape, k=extent))) - connections = connections - + # Calculate the row and column of the given number + row = (number - 1) // b + col = (number - 1) % b + + # Find all adjacent numbers within the extent + adjacent_numbers = [] + for dr in range(-extent, extent + 1): # Rows within the extent + for dc in range(-extent, extent + 1): # Columns within the extent + if dr == 0 and dc == 0: + continue # Skip the number itself + new_row, new_col = row + dr, col + dc + if 0 <= new_row < a and 0 <= new_col < b: + adjacent_num = new_row * b + new_col + 1 + adjacent_numbers.append(adjacent_num) + + return np.array(adjacent_numbers, dtype=int) + connections = np.zeros((len(channelnames),) * 2, dtype=float) + channelnames = np.array(channelnames) + contact_arrays = np.array([x.rstrip('0123456789') for x in channelnames]) + contacts, ch_per_contact = np.unique([x.rstrip('0123456789') for x in channelnames], return_counts=True) + if extent is None: + # Common average referencing per electrode array (ECoG grid or sEEG shank) + # CAR will end up subtracting parts of channel ch from itself + for contact, num_ch in zip(contacts, ch_per_contact): + for ch in range(1,num_ch+1): + curr = np.where(channelnames==f'{contact}{ch}')[0] + inds = np.where(contact_arrays==contact)[0] + connections[curr,inds] = 1 + elif extent < 1: + raise ValueError(f'Invalid extent. Must be no less than 1 but got extent={extent}') + else: + # Local average referencing within each electrode array + # LAR will NOT subtract parts of channel ch from itself + for contact, num_ch in zip(contacts, ch_per_contact): + for ch in range(1,num_ch+1): + # Local referencing for ECoG grids + if 'grid' in contact.lower(): + side = np.sqrt(num_ch) + half_side = np.sqrt(num_ch/2) + if np.isclose(side, int(side)): + nrows, ncols = side, side + elif np.isclose(half_side, int(half_side)): + nrows, ncols = half_side, half_side*2 + else: + raise Exception('Cannot automatically determine grid layout') + adjacent = _find_adjacent_numbers(nrows, ncols, ch, extent) + curr = np.where(channelnames==f'{contact}{ch}')[0] + inds = [] + print(ch, adjacent) + for adj in adjacent: + inds.append(np.where(channelnames==f'{contact}{adj}')[0]) + + # Local referencing for sEEG shanks and strips + else: + curr = np.where(channelnames==f'{contact}{ch}')[0] + inds = [] + for cc in range(ch-extent, ch+extent+1): + if cc != ch: + inds.append(np.where(channelnames==f'{contact}{cc}')[0]) + + inds = np.concatenate(inds) + connections[curr,inds] = 1 + return connections From 8be3445d1d4c6917f68867088237bd7e2cd688d3 Mon Sep 17 00:00:00 2001 From: Vinay Raghavan <42253618+vinaysraghavan@users.noreply.github.com> Date: Fri, 10 Jan 2025 15:59:00 -0500 Subject: [PATCH 02/10] Update rereference.py Added input `grid_sizes` to allow custom input of grid sizes. After talking with Max Nentwich, I realized that, while square and 1 x 2 rectangle grids are relatively standard/common, grids are occasionally cut, so custom shapes as input is needed. --- naplib/preprocessing/rereference.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/naplib/preprocessing/rereference.py b/naplib/preprocessing/rereference.py index 431081d3..0e76f13c 100644 --- a/naplib/preprocessing/rereference.py +++ b/naplib/preprocessing/rereference.py @@ -116,7 +116,7 @@ def _rereference(data_arr, method='avg', return_ref=False): return data_rereferenced -def make_contact_rereference_arr(channelnames, extent=None): +def make_contact_rereference_arr(channelnames, extent=None, grid_sizes={}): """ Create grid which defines re-referencing scheme based on electrodes being on the same contact as each other. @@ -134,6 +134,10 @@ def make_contact_rereference_arr(channelnames, extent=None): from each other (inclusive) are still grouped together. For example, if ``extent=1``, only the nearest electrode on either side of a given electrode on the same contact is still grouped with it. This ``extent=1`` produces the traditional local average reference scheme. + grid_sizes : dict, optional, default={} + If provided, contains {'contact_name': (nrow, ncol)} values for any known ECoG grid sizes. + E.g. {'GridA': (8, 16)} indicates that electrodes on contact 'GridA' are arranged in an 8 x 16 grid, + which is needed to determine adjacent electrodes for local average referencing with ``extent >= 1``. Returns ------- @@ -191,12 +195,17 @@ def _find_adjacent_numbers(a, b, number, extent): if 'grid' in contact.lower(): side = np.sqrt(num_ch) half_side = np.sqrt(num_ch/2) - if np.isclose(side, int(side)): + # Check grid_sizes dict + if contact in grid_sizes: + nrows, ncols = grid_sizes[contact] + # Assume a square + elif np.isclose(side, int(side)): nrows, ncols = side, side + # Assume a 1 x 2 rectangle elif np.isclose(half_side, int(half_side)): nrows, ncols = half_side, half_side*2 else: - raise Exception('Cannot automatically determine grid layout') + raise Exception(f'Cannot determine {contact} layout. Please include layout in `grid_sizes`') adjacent = _find_adjacent_numbers(nrows, ncols, ch, extent) curr = np.where(channelnames==f'{contact}{ch}')[0] inds = [] From ca3e27031a1f629581227488a63eaf78e3adac0b Mon Sep 17 00:00:00 2001 From: Vinay Raghavan <42253618+vinaysraghavan@users.noreply.github.com> Date: Fri, 10 Jan 2025 16:20:16 -0500 Subject: [PATCH 03/10] Update test_rereference.py Update test of make_contact_rereference_arr to improve coverage --- tests/preprocessing/test_rereference.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/preprocessing/test_rereference.py b/tests/preprocessing/test_rereference.py index d9fa02a5..ce5068e4 100644 --- a/tests/preprocessing/test_rereference.py +++ b/tests/preprocessing/test_rereference.py @@ -5,8 +5,8 @@ def test_create_contact_rereference_arr(): expected = np.array([[1,1,0,0],[1,1,0,0],[0,0,1,1],[0,0,1,1]]) - g = ['LT1','LT2','RT1','RT2'] - arr = make_contact_rereference_arr(g) + g = ['LT1','LT2','RT1','RT2'] + [f'GridA{n}' for n in range(1,5)] + arr = make_contact_rereference_arr(g, extent=1) assert np.allclose(expected, arr) def test_rereference_avg(): From d0879f438fafb2c5994de6c857632ef161764fd9 Mon Sep 17 00:00:00 2001 From: Vinay Raghavan <42253618+vinaysraghavan@users.noreply.github.com> Date: Fri, 10 Jan 2025 16:22:33 -0500 Subject: [PATCH 04/10] Update test_rereference.py Change expected test outcome --- tests/preprocessing/test_rereference.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/tests/preprocessing/test_rereference.py b/tests/preprocessing/test_rereference.py index ce5068e4..77cdfcce 100644 --- a/tests/preprocessing/test_rereference.py +++ b/tests/preprocessing/test_rereference.py @@ -4,7 +4,15 @@ def test_create_contact_rereference_arr(): - expected = np.array([[1,1,0,0],[1,1,0,0],[0,0,1,1],[0,0,1,1]]) + expected = np.array([[1,1,0,0,0,0,0,0], + [1,1,0,0,0,0,0,0], + [0,0,1,1,0,0,0,0], + [0,0,1,1,0,0,0,0], + [0,0,0,0,0,1,1,1], + [0,0,0,0,1,0,1,1], + [0,0,0,0,1,1,0,1], + [0,0,0,0,1,1,1,0], + ]) g = ['LT1','LT2','RT1','RT2'] + [f'GridA{n}' for n in range(1,5)] arr = make_contact_rereference_arr(g, extent=1) assert np.allclose(expected, arr) From d43bb2b27d76d1ecb2d519b262c54d06741b67bf Mon Sep 17 00:00:00 2001 From: Vinay Raghavan <42253618+vinaysraghavan@users.noreply.github.com> Date: Fri, 10 Jan 2025 16:32:31 -0500 Subject: [PATCH 05/10] Updated rerefererence array test to improve coverage --- tests/preprocessing/test_rereference.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/preprocessing/test_rereference.py b/tests/preprocessing/test_rereference.py index 77cdfcce..c24f07ff 100644 --- a/tests/preprocessing/test_rereference.py +++ b/tests/preprocessing/test_rereference.py @@ -4,10 +4,10 @@ def test_create_contact_rereference_arr(): - expected = np.array([[1,1,0,0,0,0,0,0], - [1,1,0,0,0,0,0,0], - [0,0,1,1,0,0,0,0], - [0,0,1,1,0,0,0,0], + expected = np.array([[0,1,0,0,0,0,0,0], + [1,0,0,0,0,0,0,0], + [0,0,0,1,0,0,0,0], + [0,0,1,0,0,0,0,0], [0,0,0,0,0,1,1,1], [0,0,0,0,1,0,1,1], [0,0,0,0,1,1,0,1], From a5f4d5eb4243ecbbfccbae047dfac943ba88235d Mon Sep 17 00:00:00 2001 From: Vinay Raghavan <42253618+vinaysraghavan@users.noreply.github.com> Date: Fri, 10 Jan 2025 16:50:42 -0500 Subject: [PATCH 06/10] Improve test coverage for make contact rereference array --- tests/preprocessing/test_rereference.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/tests/preprocessing/test_rereference.py b/tests/preprocessing/test_rereference.py index c24f07ff..a5c1cdea 100644 --- a/tests/preprocessing/test_rereference.py +++ b/tests/preprocessing/test_rereference.py @@ -4,6 +4,11 @@ def test_create_contact_rereference_arr(): + + g = ['LT1','LT2','RT1','RT2'] + arr = make_contact_rereference_arr(g) + assert np.allclose(expected, arr) + expected = np.array([[0,1,0,0,0,0,0,0], [1,0,0,0,0,0,0,0], [0,0,0,1,0,0,0,0], @@ -13,9 +18,20 @@ def test_create_contact_rereference_arr(): [0,0,0,0,1,1,0,1], [0,0,0,0,1,1,1,0], ]) + expected1 = np.array([[1,1,0,0,0,0,0,0], + [1,1,0,0,0,0,0,0], + [0,0,1,1,0,0,0,0], + [0,0,1,1,0,0,0,0], + [0,0,0,0,1,1,1,1], + [0,0,0,0,1,1,1,1], + [0,0,0,0,1,1,1,1], + [0,0,0,0,1,1,1,1], + ]) g = ['LT1','LT2','RT1','RT2'] + [f'GridA{n}' for n in range(1,5)] arr = make_contact_rereference_arr(g, extent=1) + arr1 = make_contact_rereference_arr(g) assert np.allclose(expected, arr) + assert np.allclose(expected1, arr1) def test_rereference_avg(): arr = np.array([[1,1,0,0],[1,1,0,0],[0,0,1,1],[0,0,1,1]]) From d3b60da54b861ba85401d2d58221f8f7e33ab750 Mon Sep 17 00:00:00 2001 From: Vinay Raghavan <42253618+vinaysraghavan@users.noreply.github.com> Date: Fri, 10 Jan 2025 16:53:27 -0500 Subject: [PATCH 07/10] Improve test coverage further --- tests/preprocessing/test_rereference.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/preprocessing/test_rereference.py b/tests/preprocessing/test_rereference.py index a5c1cdea..1c5dedce 100644 --- a/tests/preprocessing/test_rereference.py +++ b/tests/preprocessing/test_rereference.py @@ -27,8 +27,8 @@ def test_create_contact_rereference_arr(): [0,0,0,0,1,1,1,1], [0,0,0,0,1,1,1,1], ]) - g = ['LT1','LT2','RT1','RT2'] + [f'GridA{n}' for n in range(1,5)] - arr = make_contact_rereference_arr(g, extent=1) + g = ['LT1','LT2','GridA1','GridA2'] + [f'GridB{n}' for n in range(1,5)] + arr = make_contact_rereference_arr(g, extent=1, grid_sizes={'GridA':(1,2)}) arr1 = make_contact_rereference_arr(g) assert np.allclose(expected, arr) assert np.allclose(expected1, arr1) From 6897058b1c70241f348b0faeb44d777fb2308d4e Mon Sep 17 00:00:00 2001 From: Vinay Raghavan <42253618+vinaysraghavan@users.noreply.github.com> Date: Fri, 10 Jan 2025 16:54:49 -0500 Subject: [PATCH 08/10] Remove extra lines --- tests/preprocessing/test_rereference.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/tests/preprocessing/test_rereference.py b/tests/preprocessing/test_rereference.py index 1c5dedce..39c4bfa3 100644 --- a/tests/preprocessing/test_rereference.py +++ b/tests/preprocessing/test_rereference.py @@ -4,11 +4,6 @@ def test_create_contact_rereference_arr(): - - g = ['LT1','LT2','RT1','RT2'] - arr = make_contact_rereference_arr(g) - assert np.allclose(expected, arr) - expected = np.array([[0,1,0,0,0,0,0,0], [1,0,0,0,0,0,0,0], [0,0,0,1,0,0,0,0], From fae2e2e17e2339fed45b4a43632f05ab38b71308 Mon Sep 17 00:00:00 2001 From: Vinay Raghavan <42253618+vinaysraghavan@users.noreply.github.com> Date: Mon, 13 Jan 2025 12:09:11 -0500 Subject: [PATCH 09/10] Update rereference.py Remove debugging print statement --- naplib/preprocessing/rereference.py | 1 - 1 file changed, 1 deletion(-) diff --git a/naplib/preprocessing/rereference.py b/naplib/preprocessing/rereference.py index 0e76f13c..8aac19f3 100644 --- a/naplib/preprocessing/rereference.py +++ b/naplib/preprocessing/rereference.py @@ -209,7 +209,6 @@ def _find_adjacent_numbers(a, b, number, extent): adjacent = _find_adjacent_numbers(nrows, ncols, ch, extent) curr = np.where(channelnames==f'{contact}{ch}')[0] inds = [] - print(ch, adjacent) for adj in adjacent: inds.append(np.where(channelnames==f'{contact}{adj}')[0]) From cb550542e1899764e57870734f3ac027fab1eef2 Mon Sep 17 00:00:00 2001 From: Vinay Raghavan <42253618+vinaysraghavan@users.noreply.github.com> Date: Wed, 15 Jan 2025 10:41:30 -0500 Subject: [PATCH 10/10] Update rereference.py Removed the assumption that channels 1-N are always provided for a given contact. Now, any channel numbers on a contact can be provided, and local average referencing will be robust to this. This is useful in cases where depth electrodes are inserted through the brain, e.g. some electrodes are in the CSF beyond the brain, some electrodes are in the brain, and some are in the skull/scalp, since only electrodes in the brain should be used for re-referencing. It now also notifies the user if channels cannot properly be re-referenced because `extent` is too small due to a lack of spatially local sites. --- naplib/preprocessing/rereference.py | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/naplib/preprocessing/rereference.py b/naplib/preprocessing/rereference.py index 8aac19f3..6ebf8350 100644 --- a/naplib/preprocessing/rereference.py +++ b/naplib/preprocessing/rereference.py @@ -134,6 +134,7 @@ def make_contact_rereference_arr(channelnames, extent=None, grid_sizes={}): from each other (inclusive) are still grouped together. For example, if ``extent=1``, only the nearest electrode on either side of a given electrode on the same contact is still grouped with it. This ``extent=1`` produces the traditional local average reference scheme. + The default ``extent=None`` produces the traditional common average reference scheme. grid_sizes : dict, optional, default={} If provided, contains {'contact_name': (nrow, ncol)} values for any known ECoG grid sizes. E.g. {'GridA': (8, 16)} indicates that electrodes on contact 'GridA' are arranged in an 8 x 16 grid, @@ -172,15 +173,21 @@ def _find_adjacent_numbers(a, b, number, extent): adjacent_numbers.append(adjacent_num) return np.array(adjacent_numbers, dtype=int) + connections = np.zeros((len(channelnames),) * 2, dtype=float) channelnames = np.array(channelnames) contact_arrays = np.array([x.rstrip('0123456789') for x in channelnames]) - contacts, ch_per_contact = np.unique([x.rstrip('0123456789') for x in channelnames], return_counts=True) + contacts = np.unique(contact_arrays) + # Determine the channel numbers on each contact + ch_per_contact = {contact:[int(x.replace(contact,'')) for x in channelnames + if x.rstrip('0123456789')==contact] + for contact in contacts} + if extent is None: # Common average referencing per electrode array (ECoG grid or sEEG shank) # CAR will end up subtracting parts of channel ch from itself - for contact, num_ch in zip(contacts, ch_per_contact): - for ch in range(1,num_ch+1): + for contact in contacts: + for ch in ch_per_contact[contact]: curr = np.where(channelnames==f'{contact}{ch}')[0] inds = np.where(contact_arrays==contact)[0] connections[curr,inds] = 1 @@ -189,10 +196,11 @@ def _find_adjacent_numbers(a, b, number, extent): else: # Local average referencing within each electrode array # LAR will NOT subtract parts of channel ch from itself - for contact, num_ch in zip(contacts, ch_per_contact): - for ch in range(1,num_ch+1): + for contact in contacts: + for ch in ch_per_contact[contact]: # Local referencing for ECoG grids if 'grid' in contact.lower(): + num_ch = len(ch_per_contact[contact]) side = np.sqrt(num_ch) half_side = np.sqrt(num_ch/2) # Check grid_sizes dict @@ -221,6 +229,9 @@ def _find_adjacent_numbers(a, b, number, extent): inds.append(np.where(channelnames==f'{contact}{cc}')[0]) inds = np.concatenate(inds) - connections[curr,inds] = 1 + if len(inds) < 1: + print(f'{contact}{cc} has no re-references.') + else: + connections[curr,inds] = 1 return connections