"RuntimeWarning: overflow encountered in scalar multiply" in axonarawio

[SusanL82]

Hi everyone,

I'm encountering the following error/warning in a spike-extraction script I made a while ago:
I am using the detect_peaks function on a concatenated Axona recording, which has worked fine until recently.
I've attached the script below. It's a bit long and not very elegant, but the error happens in the 'detect peaks' section, line 95-99. The main goal of the script is to detect spikes and export them so that we can use them in a manual clustering program (SpikeSort3D), I usually run it in Spyder.

For some reason, I now encounter the following:

detect peaks using locally_exclusive:  17%|#6        | 1240/7390 [00:33<02:42, 37.85it/s]
C:\Users\leemburg\AppData\Local\anaconda3\envs\SInew\lib\site-packages\neo\rawio\axonarawio.py:384: RuntimeWarning: overflow encountered in scalar multiply
  offset = i_start // 3 * (bin_dict["bytes_packet"] // 2)
detect peaks using locally_exclusive: 100%|##########| 7390/7390 [02:48<00:00, 43.92it/s] 

I have updated to Python 3.10 and installed the latest version of spikeinterface (within the latest anaconda) and the warning/error persists. What can I do about it?

import spikeinterface.extractors as se
import spikeinterface as si
import numpy as np
from probeinterface import read_prb
from spikeinterface.sortingcomponents.peak_detection import detect_peaks
from spikeinterface.preprocessing import bandpass_filter
from scipy.io import savemat
from tqdm import tqdm


InFolder = "E:/AxonaToNLXtest/"
OutFolder = "E:/AxonaToNLXtest/pos/"
Probepath ="C:/Users/leemburg/Desktop/OEphystest/"
ChanList = 'Tetlist.txt' # text file listing good and bad channels
BaseName = 'HD263-2811_'
TetList = [2,4,6] #analyse only these tetrodes (1-based, as on drive)
numsess = 3 #number of sessions to read (e.g: numsess = 4 reads bin 01 to 04)

spike_thresh = 5 # detection threshold for spike detection is spike_thresh* signal SD
spike_sign = 'pos' #detect peaks. can be: 'neg', 'pos', 'both'

# !!! need to manually edit the .set files so that all channels get read. 

##############################################################################

tetgrouping = np.array([0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6,7,7,7,7,8,8,8,8,
                        9,9,9,9,10,10,10,10,11,11,11,11,12,12,12,12,13,13,13,13,14,14,14,14,
                        15,15,15,15,15])

# read bad channel list, set bad channels to 17
tetlist = np.loadtxt(InFolder + '/' + ChanList,
                 delimiter=",")

for tetnum in range(16):
        thistet = np.where(tetgrouping==tetnum)
        thistet = np.array(thistet)
        thesewires = tetlist[tetnum,1:5]
        badwires = np.where(thesewires==0)
        badwires = np.array(badwires)
        badchans = thistet[0][badwires]
        tetgrouping[badchans]=17



# make file list
allsess = range(numsess) #count sessions zero-based
files = list()
for f in allsess:
   thissess = allsess[f]+1 #1-based session number
   InName = InFolder+BaseName+"%02d" % thissess+".bin"
   list.append(files, InName)


# make list of SI extractor recordings
recording_list = list()
for f in allsess:
   list.append(recording_list, se.AxonaRecordingExtractor(files[f]))

# merge recording objects 
MergeRec = si.concatenate_recordings(recording_list)

# add highpass filter
MergeRec_f = bandpass_filter(MergeRec, freq_min=300.0, freq_max=6000.0, margin_ms=5.0, dtype=None)

all_chan_ids = MergeRec_f.get_channel_ids()

# select a tetrode
for tetnum in TetList:

    tet_chan_ids = all_chan_ids[np.where(tetgrouping == tetnum-1)]
    tetname = TetList[tetnum-1]

    if np.size(tet_chan_ids)>2:
        
        # 4-wire tetrode
        if np.size(tet_chan_ids) == 4:
           new_chans = [0,1,2,3]
           probename = "tet4_probe.prb"
        
        # 3-wire tetrode
        if np.size(tet_chan_ids) == 3:
            new_chans = [0,1,2]
            probename = "tet3_probe.prb"
        
    myprobe = read_prb(Probepath + "/" + probename)

    #select channels and add probe
    thistet = MergeRec_f.channel_slice(tet_chan_ids, renamed_channel_ids=new_chans)
    thistet = thistet.set_probegroup(myprobe)

    # preprocess (filter)
    thistet_f = bandpass_filter(thistet, freq_min=600, freq_max=6000)

    #detect peaks
    detectradius = 30 #tetrode channel group is 10x10um square, this radius should capture all spikes in this channelgroup
    #2 versions of function from 2 versions of spikeinterface. Yuck.
    #TetPeaks = detect_peaks(thistet_f, method='locally_exclusive', peak_sign=spike_sign, detect_threshold=spike_thresh, exclude_sweep_ms=0.1, local_radius_um=detectradius, noise_levels=None,)
    TetPeaks = detect_peaks(thistet_f, method='locally_exclusive', peak_sign=spike_sign, detect_threshold=spike_thresh, exclude_sweep_ms=0.1, radius_um=detectradius, noise_levels=None,)


    # get waveforms. must be 32 samples long, I'm setting peak at sample 8 (same as NLX)
    prepeak = 8
    postpeak = 24

    print('extracting waveforms')

    allWFs = np.zeros([32, 4, len(TetPeaks['sample_index'])], dtype='int16')
    for p in tqdm(range(len(TetPeaks['sample_index'])), desc="collecting waveforms"):
        sf = TetPeaks['sample_index'][p] - prepeak
        ef = TetPeaks['sample_index'][p] + postpeak

        thisWF = thistet_f.get_traces(segment_index=None, start_frame=sf, end_frame=ef)

        #write only complete spike waveforms (might skip last spike if too short)
        if np.size(thisWF, 0) == 32:
           allWFs[:, :, p] = thisWF

        del thisWF

   #save peaks to mat file (for later processing with Mat2NlxSpike to generate .ntt file)
    print('saving to mat file')
    outname = OutFolder+"tt"+str(tetname) + '.mat'
    savemat(outname, {"Timestamps":TetPeaks['sample_index'], "Spikes": allWFs})        

[zm711]

Hey @SusanL82 it looks like you had reported this previously #1878. Could you open the issue over on Neo. These types of overflow are operating system dependent and often are due to the fact that a numpy rather than python scalar is being used. Since this needs to be fixed on the Neo side opening an issue there will allow us to fix that in the actual package and then propagate it here. Let us know if you have questions!

Link to neo for report issue

[SusanL82]

O my gosh, you're totally right! I'm now not sure how I made the error go away last time... maybe just by running on a different PC...
I'm also not sure why I didn't find my own previous issue.
I'll post it for NEO.

[zm711]

It could also be size of the dataset. Overflow doesn't happen until it does and with datasets getting bigger and bigger maybe you're just tipping over the limits more often now.

[zm711]

@alejoe91 /@samuelgarcia , you actually want to take a look at this one. I'm not seeing how we can get a scalar overflow. The only input into axonarawio that could do this is i_start, but in this case the i_start provided would have to be from a call to get_traces inside of detect_peaks which might be a problem with the node_pipeline? Wanna confirm that?

[samuelgarcia]

Hi.
Thank you Susan for report.
And thank you Zach for the link with nodepipeline, maybe yes. I need to check. Keep this open.

[h-mayorquin]

Could you test #2854?

I think that should fix your issue.

[SusanL82]

In case it still helps, I've shared my data here: https://owncloud.cesnet.cz/index.php/s/u35wSdcL6qCU4Ie

(also: NeuralEnsemble/python-neo#1475 (comment))

[zm711]

I've tried to download so we could test, but I ended up getting an invalid zip error. Probably some corrupted communication issue between the cloud storage and my computer, but I haven't been able to test this yet. Not sure if you wanted to try to pull this down and test yourself @h-mayorquin or wait until @SusanL82 can test it?

[SusanL82]

OwnCloud has been a bit weird about one of the bin files. I thought I'd fixed it, but maybe not.
Does it work from here instead? dropbox share link (I've just uploaded the files when I post this, it'll take a little while to sync)

[zm711]

Howdy @SusanL82,

I'm a little busy for the next couple days. But I'm happy to try working on this as soon as I can!