Repository navigation
Bridging between ByCycle and NeuroDSP functionalities #148
Description
Dear ByCycle team,
I am writing to inquire about how to bridge/transition between ByCycle and NeuroDSP. More specifically, I have an array with multiple sub-arrays corresponding to different experimental subjects that I pass through the cycle-by-cycle algorithm, which I store in the "BycycleGroup" object. As a quick background -
- Here, a "BycycleGroup" object is initialized with specific thresholds and parameters. The fit method is then called on the object, passing the signal data (sigs), the sampling frequency (fs), and a frequency range (f_alpha).
bg_dict[key] = BycycleGroup(thresholds=thresholds, center_extrema='peak')
bg_dict[key].fit(sigs, fs, f_alpha, axis=0)
- I can then plot a specific sub-array and visualize cycle detection using the code-block below (see attached image for example of output).
# Specify the key for the sub-array you want to plot
sub_array_key = 'KO_m13_week_12_ch1_filtered_2-25Hz'
# Ensure the sub_array_key exists in bg_dict
if sub_array_key not in bg_dict:
raise ValueError(f"Key {sub_array_key} not found in bg_dict")
# Get the BycycleGroup object for the specified sub-array
bg1 = bg_dict[sub_array_key]
# Number of channels in the selected BycycleGroup object
num_channels = len(bg1.df_features)
ch_names = [f"Channel {i+1}" for i in range(num_channels)] # Generate channel names for plotting
# Plot the cycle analysis for each channel
for idx in tqdm(range(num_channels), desc='Processing channels'):
bg1[idx].plot(xlim=(0, 5), figsize=(16, 3))
plt.title(ch_names[idx])
plt.show()
print(f"Processed channel {ch_names[idx]}")
- My questions are as follows and relate to the properties of the "ByCycleGroup" object :
- How do I switch to a burst visualization that includes only the black and red (overlaid) traces used in NeuroDSP (see below). Can you pass the ByCycleGroup object into something like the "plot_bursts" function in NeuroDSP (see below) and get just a red/black plot and then extract burst statistics? Or is there a way to extract the black and red traces directly from the "ByCycleGroup" object. Once I check the quality of the peak detection, I don't need to see the amp fraction, amp consistency, period consistency, and monotonicity, nor their corresponding rectangular colored highlighting.
- How do I access segments that are, or are-not bursts from the "ByCycleGroup" object and compute lagged coherence on these segments. Do I need to pass the "ByCycleGroup" object through "NeuroDSP" - using the example provided in your tutorial (see below) ? Alternatively, if I have a data array containing so much information about every cycle, (i.e. time_peak, time_trough, volt_peak, volt_trough, time_decay, time_rise, volt_decay, volt_rise, volt_amp, time_rdsym, time_ptsym), can lagged coherence (or even window matching) be applied to the array instead?
- To be clear, I already can access the "ByCycleGroup" array, which has the "is_burst" column, and from which I can take all consecutive indices that are "TRUE" - and their corresponding time values. This allows me to highlight and extract the beginning and end of bursts from the filtered timeseries (used as input to ByCycle). Afterwards, I can apply lagged coherence analysis on just the epoch containing the burst, but this seems like a sloppy work-around.
Activity
ryanhammonds commentedon Aug 12, 2024
Hi @davidglitvin, sorry for the delayed response. Lets start with some simulated data to answer your questions:
is_burstarray below isn't easy accessible and has to be created. This could added more conveniently added with a property decorator on ais_burstmethod on the base or 1d class. Or even just a utility function that computes theis_oscvariable below. PRs are welcome!Alternatively, you could use this since it's built in, but it doesn't produce the same plots as plot_bursts:
chunkvariable, this variable is defined inlagged_coherence_1freq, see neurodsp's source. I would copy these functions since there isn't an interface for user-definedchunk. Putative cycles are defined from bg[sig_ind].df_features.iloc[i_cyc_start].sample_last_trough to bg[sig_ind].df_features.iloc[i_cyc_start].sample_next_trough. These are associated with theis_burstcolumn, which should be False for this.I would then chunk, e.g. reshape, sig_non_burst above based on the length of desired lagged coherence window. You'll probably need to slice off a few samples so it's evenly divisible by the desired window length. The code above is for one non-burst segment. The length of each chunk should be constant, so the end result will be a loop of the above and a single 2d array of non-burst chunks. Then continue from this line.
davidglitvin commentedon Sep 5, 2024
Ok, thanks for the help! I had another question that is somewhat related, and probably doesn't warrant opening a new issue for this github repository. For the FOOOFGroup array containing information about the cycles that are (or are not) parts of bursts, there are a number of features that describe each cycles' shape (i.e. amp_fraction, amp_consistency, period_consistency, monotonicity, period, time_peak, time_trough, time_decay, time_rise, time_rdsym, time_ptsym, etc). I am noticing that in my data the bursts have several classes of cycle groups with common values for the attributes. Have you tried using something like ICA, or machine learning to classify discrete clusters of cycles within bursts?
ryanhammonds commentedon Sep 5, 2024
You're welcome!
I've used kmeans in the past to cluster cycles. I did this by resampling cycles to have the same number of samples per cycle. Then passed the resampled cycles into kmeans. Resampling is probably not ideal unless each cycle has a sufficient number of samples and there is little variation in the number of samples per cycle.
Using bycycle features + a clustering method sounds like a good idea and would be a nice use case! Let me know how it works out. If the data is open source or simulate-able, then it may be a nice addition to the examples page.