Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Conversion to logical from gpuArray is not possible. #22

Closed
alejoe91 opened this issue Mar 7, 2019 · 14 comments · Fixed by #595
Closed

Conversion to logical from gpuArray is not possible. #22

alejoe91 opened this issue Mar 7, 2019 · 14 comments · Fixed by #595

Comments

@alejoe91
Copy link

alejoe91 commented Mar 7, 2019

Hi!

First of all thanks for the amazing tool :)

I'm trying to set it up in SpikeInterface (https://github.com/SpikeInterface), a toolbox for easily running and comparing several sorters. Kilosort2 is implemented here: https://github.com/SpikeInterface/spiketoolkit/tree/kilosort2/spiketoolkit/sorters/kilosort2

It runs fine on simulated data until I get this error:

Time   0s. Determining good channels.. 
Warning: The CUDA driver must recompile the GPU libraries because your device
is more recent than the libraries. Recompiling can take several minutes. 
> In get_good_channels (line 20)
  In preprocessDataSub (line 27)
  In kilosort2_master (line 20)
  In run (line 91) 
found 3473 threshold crossings in 201.17 seconds of data 
found 0 bad channels 
Time   3s. Computing whitening matrix.. 
Getting channel whitening matrix... 
Channel-whitening filters computed. 
Time   4s. Loading raw data and applying filters... 
Time   6s. Finished preprocessing 275 batches. 
Obtained 7 PC waveforms in 0.17 seconds 
time 0.16, pre clustered 1 / 275 batches 
time 0.01, compared 1 / 275 batches 
time 3.55, Re-ordered 275 batches. 
Time   4s. Optimizing templates ...
Conversion to logical from gpuArray is not possible.

Error in median (line 153)
        if isnan(x(nCompare))        % Check last index for NaN

Error in learnAndSolve8b (line 263)
            toc, ibatch, niter, Nfilt, sum(nsp), median(mu), numel(st0), ndrop)

Error in kilosort2_master (line 26)
rez = learnAndSolve8b(rez);

I was able to run the eMouse_drift example and I'm running on Ubuntu 18.05, CUDA 9.1, MATLAB 2018b.

Have you experienced this error before?

Thanks!
Alessio

@marius10p
Copy link
Contributor

Haven't gotten this error. You could replace the median with a nanmedian, but that likely indicates some other problem elsewhere. Can you try loading the data in the GUI to make sure the channel map is configured correctly?

Also, the preprocessing seems too fast. How many channels do you have? What does the re-ordered drift matrix look like?

Finally, if your simulation doesn't have drift, it kind of beats the purpose of Kilosort2. We wrote a simulation with drift that you can find in the eMouse_drift folder and the wiki.

@alejoe91
Copy link
Author

alejoe91 commented Mar 7, 2019

Hi Marius,

thanks for the quick response. It's tetrode data, so 4 channels and 600 seconds.
The channel locations are dummy: [[0,0], [0,1], [0,2], [0,3]]. Should I increase the distance?
I will try to load it with the GUI tomorrow, maybe it's a parsing problem. Is the data format exactly the same as kilosort 1?

We will test it on real data and simulations with drifting too.
On data without drifting, is the performance expected to be the same as kilosort 1?

Thanks again!
Alessio

@marius10p
Copy link
Contributor

The distance should be fine. The data format is indeed the same. Without drift it should be similar to KS1, but it should also be better at avoiding local minima of the clustering problem and over splitting less. That said, I haven't compared it on data without drift, so please let me know how it goes.

@alejoe91
Copy link
Author

alejoe91 commented Mar 8, 2019

I tried on another batch of simulated data: 100 electrodes (10x10), 300s, 15um pitch. (So probably the error was due to the previous simulated data).
Kilosort2 works, but it only loads half of the data.
The problem is in line 69 of get_good_channels.m:

ibatch = ibatch + ceil(Nbatch/100);

In my case, Nbatch is 147 (samplesToRead=9600000, NT=65600, ntbuff=64), so the ceil(Nbatch/100) gives 2 and it skips every other batch.

What is the reason the ibatch update is not simply ibatch = ibatch + 1;?

Thanks,
Alessio

@marius10p
Copy link
Contributor

Some of the initial steps of the pipeline need just a fraction of the data, and in this case they are designed to load about 100 batches from throughout the recording. This step computes threshold crossings, and throws out channels with a very low rate (<0.1Hz by default).

@alejoe91
Copy link
Author

alejoe91 commented Mar 8, 2019

I see! It makes sense. Thanks for specifying that.

We successfully ran kilosort2 on linear probe data and it gives quite some clusters with really few spikes (<30). Should the ops.minFR allow to reject these clusters? We have 3 min recordings, so setting it to 0.2 would remove clusters with less than 0.2*180 = 36 spikes?
Is this how it is supposed to be used for?

Thanks
Alessio

@marius10p
Copy link
Contributor

Sorry for the delay. Is the totality of your recording 3 min long, or do you have multiple 3 min blocks? If the latter, I suggest concatenating and spike sorting together.

<30 is not that few for a 3 minute recording. Most neurons in cortex have low baseline firing rates. The minFR option checks periodically during optimization if a cluster is below that, and removes it. Because of that, it's not a hard threshold: sometimes you might still get some clusters with fewer than minFR spikes, sometimes such clusters might be introduced in the splitting step after optimization. If you really don't want those, discard them posthoc?

It could also be that something is wrong. How many such clusters do you get relative to good clusters?

@alejoe91
Copy link
Author

So we observed this both in real data and simulated data. Around 50% of the units found by kilosort2 were with very few spikes and most of the times they looked similar (at least in waveforms) to other templates. Maybe this can be related to #29?
On the simulated dataset kilosort1 performed better without this unit splitting.

For now we will fix it posthoc :)

@marius10p
Copy link
Contributor

Is this for the 3 minute recordings, or also for normal length recordings?

@marius10p
Copy link
Contributor

For very short recordings, the CCG and ACG calculations will be noisy, and maybe that's why it oversplits. Do you get a lot of splits at the end? Have you checked that the GUI loads the data correctly?

@alejoe91
Copy link
Author

It happened both in the 3 min recording on a 32-channel laminar probe and on a 5min recording on simulated data with 100 closely spaced electrodes.
I checked the data and they get loaded correctly on the GUI.

I will test out longer simulations (10min) and see if that reduces the problem.
I will tell you the results soon.

@marius10p
Copy link
Contributor

try 30 minutes.

@alejoe91
Copy link
Author

👍

@alejoe91
Copy link
Author

The problem is indeed that the duraiton is too short. Closing

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants