Make average neuropixel lfp array + max coherence & lag array generation #583
Conversation
| Parameters: | ||
| - lags: List of time lags (in seconds) to apply to the channels in lfp_array1. | ||
| - lfp_array1: First LFP array (e.g., motor cortex data) with shape (n_channels, n_timepoints). | ||
| - lfp_array2: Second LFP array (e.g., pre-motor cortex data) with shape (n_channels, n_timepoints). | ||
| - frequency_range: Frequency range of interest (e.g., [0, 4] Hz for delta band). | ||
|
|
||
| Returns: | ||
| - max_coherence_matrix: Matrix of maximum coherence values for each channel pair. | ||
| - lag_of_max_coherence_matrix: Matrix of corresponding lags (in ms) for max coherence values. |
There was a problem hiding this comment.
Please make the docstring google format or else the readthedocs won't work. (https://sphinxcontrib-napoleon.readthedocs.io/en/latest/example_google.html)
| - lag_of_max_coherence_matrix: Matrix of corresponding lags (in ms) for max coherence values. | ||
| """ | ||
| # Taper parameters | ||
| NW, BW, step, fk, fs = 0.5, 4, 0.005, 200, 2500 |
There was a problem hiding this comment.
Can you make these input arguments to the function?
| - lfp_array1: First LFP array (e.g., motor cortex data) with shape (n_channels, n_timepoints). | ||
| - lfp_array2: Second LFP array (e.g., pre-motor cortex data) with shape (n_channels, n_timepoints). |
There was a problem hiding this comment.
Lab convention is to put the time axis first (even if this isn't necessarily the optimal way)
| max_coherence_matrix = np.zeros((num_channels_lfp1, num_channels_lfp2)) | ||
| lag_of_max_coherence_matrix = np.zeros((num_channels_lfp1, num_channels_lfp2)) | ||
|
|
There was a problem hiding this comment.
Can you initialize these arrays as np.nan values to make it obvious if the answer is 0 or if the value just didn't get filled in? You can just add "*np.nan" to the end of each of these lines.
| """ | ||
| Shift a signal by a specified lag in seconds, converting to samples. | ||
| Positive lag shifts forward, negative lag shifts backward. | ||
| """ |
There was a problem hiding this comment.
Please add information about the function args and what it returns to the docstring.
| max_coherence_matrix = np.zeros((num_channels_lfp1, num_channels_lfp2)) | ||
| lag_of_max_coherence_matrix = np.zeros((num_channels_lfp1, num_channels_lfp2)) | ||
|
|
||
| def apply_lag(signal, lag, fs): |
There was a problem hiding this comment.
Can you rename this function to something more specific? Something like "circle_shift_by_time"
| """ | ||
| Compute coherence between a specific pair of channels for all specified lags. | ||
| Returns coherence values for each lag at time zero. | ||
| """ |
There was a problem hiding this comment.
Can you update the docstring as above? Also can you pull this function out of the wrapper function so that it can easily be used alone?
There was a problem hiding this comment.
Then each input value to the array can be a 1D array, which should also simplify the function.
| Compute coherence between a specific pair of channels for all specified lags. | ||
| Returns coherence values for each lag at time zero. | ||
| """ | ||
| coherence_at_zero_lag = np.empty(len(lags)) |
There was a problem hiding this comment.
Similar to before, please initialize this as np.nan values. Initializing as an empty array will do weird things.
| for idx, lag in enumerate(lags): | ||
| # Process signal from lfp_array1 and normalize | ||
| signal1 = lfp_array1[channel_lfp1, :] | ||
| signal1 = copy.deepcopy((signal1 - np.mean(signal1)) / np.std(signal1)) |
There was a problem hiding this comment.
Zscoring should be done outside of the function to make it more flexible.
| # Calculate max coherence and lags for each channel pair | ||
| for ch_lfp1 in range(num_channels_lfp1): | ||
| for ch_lfp2 in range(num_channels_lfp2): | ||
| print(f"Processing LFP1 channel {ch_lfp1}, LFP2 channel {ch_lfp2}") | ||
| delta_coherence = compute_coherence_for_pair(ch_lfp1, ch_lfp2) | ||
|
|
||
| # Find maximum coherence and corresponding lag | ||
| max_coherence = np.max(delta_coherence) | ||
| max_lag_index = np.argmax(delta_coherence) | ||
| max_lag = lags[max_lag_index] * 1000 # Convert lag to milliseconds | ||
|
|
||
| # Store max coherence and lag in matrices | ||
| max_coherence_matrix[ch_lfp1, ch_lfp2] = max_coherence | ||
| lag_of_max_coherence_matrix[ch_lfp1, ch_lfp2] = max_lag | ||
|
|
||
| return max_coherence_matrix, lag_of_max_coherence_matrix |
There was a problem hiding this comment.
This part can be a separate wrapper function.
There was a problem hiding this comment.
And can you add an optional argument called "verbose". Then, only run print statements if verbose=True
| def make_average_neuropixel_array(te_id, date, subject, aligning_index, tbefore, tafter): | ||
| """ |
There was a problem hiding this comment.
Thanks for including this function but Tomo and I are in the middle of restructuring how LFP preprocessing is done. This will change how data is loaded and stuff so it isn't necessary to include this function in github quite yet.
| self.assertEqual(overlap,(20-17)/(20-10)) | ||
|
|
||
|
|
||
| class TestCalcMaxCohAndLags(unittest.TestCase): |
There was a problem hiding this comment.
Can you make the name of the test class "test_{function_name_to_be_tested}"? This way it is easy to identify which set of tests belongs to which function.
| # Mock taper parameters and the coherence calculation function | ||
| global precondition | ||
| precondition = MagicMock() | ||
| precondition.convert_taper_parameters = MagicMock(return_value=(256, 2, 3)) # Mock return values for tapers |
There was a problem hiding this comment.
Just import the actual function from aopy.precondition.
| global calc_mt_tfcoh | ||
| calc_mt_tfcoh = MagicMock(return_value=( | ||
| np.linspace(0, 200, 100), # Mock frequency array (0-200 Hz) | ||
| np.linspace(-0.5, 0.5, 50), # Mock time array (centered around 0) | ||
| np.random.rand(100, 50) # Mock coherence array (freq x time) | ||
| )) |
There was a problem hiding this comment.
Similar to before, I'm also not sure why importing this is necessary.
| # Check if lags are converted to milliseconds in lag_of_max_coherence_matrix | ||
| self.assertTrue(np.all(np.isin(lag_of_max_coherence_matrix, np.array(self.lags) * 1000)), | ||
| "Lags should be in milliseconds and match the provided lags") | ||
|
|
There was a problem hiding this comment.
These tests are good for testing that the structure of the outputs are as expected, but it doesn't require the actual values to be as expected. After pulling out the single lagged coherence function, the test will be easier to write. You can create fake signals and then circle shift them before inputting them into your function. Then you can test is the max lag is the same as the circle shift. And if it is the same signal, it should be very high. I would start something with a single hump and zeros elsewhere, but play around with it. Let me know if you have questions.
Added a function to aopy.preproc.neuropixel called make_average_neuropixel_array as well as a function to aopy.analysis.base called calc_max_coh_and_lags as well as test cases to their corresponding test files