
\subsection{T-sne code}
For all our experiments we have used an extended version of the github available C++ t-sne code by van Maatens\\ \href{https://github.com/lvdmaaten/bhtsne/}{(https://github.com/lvdmaaten/bhtsne/)}. The bhtsne implementation in this repository allows for faster computations and makes larger data sets feasible through the use of the barnes-hut (bh) algorithm \cite{van_der_maaten_accelerating_2014}. The bh algorithm groups together distant samples into a single average sample, minimizing the number of euclidean distance computations the algorithm must perform. This step makes the calculation of the actual error minimization of the 2D landscape faster and tractable for large sample numbers (on the order of millions). Yet the calculation of the initial probabilities in the multidimensional space still requires the algorithm to calculate the euclidean distances of all sample pairs adding a significant amount to the algorithm's run time. To speed up this part we have implemented the computation of the euclidean distances between all sample pairs directly on the GPU. To do this we used the GPU euclidean distance implementation presented by Chang et al \cite{chang_compute_2008}. This allows us to run data sets of more than 10\textsuperscript{5} of spikes on a gaming desktop (i7 CPU, Titan X GPU, 64GB RAM) in just a couple of hours (see table ST1 about times for different number of samples and perplexity settings). However, storing the distances for all sample pairs on RAM soon imposes a bottleneck on the amount of samples that can be used. For N samples one requires 4*N*N (for 32 bit floats) bits of RAM, which translates to 40GB for 10\textsuperscript{5} spikes and to 1TB for 0.5*10\textsuperscript{6} samples. To overcome this bottleneck we have added the possibility to calculate the pairwise distances in groups and keep in memory only the ones that the barnes-hut algorithm will use (populating the tree structure used in the algorithm with GPU pre-calculated distances). This process requires a sorting step of the distances that increase the time required, but allows much larger sample numbers to be used.
\\For data sets that are prohibitively large, and to show that online sorting of samples is possible, we have extended the algorithm to be able to position samples on a t-sne pre-computed 2D landscape. This is done by measuring the euclidean distance of each new sample to the samples already passed through the t-sne algorithm. Then the extra samples' positions on the 2D landscape are calculated as the average of the positions of the closest 5 original samples for each new sample. For spike sorting, given that the number of spikes passed through the t-sne algorithm offers a complete representation of the spiking units in the recording, we show that the algorithm correctly places the extra samples (see Supplementary Figure 2).
\\We have also extended the python wrapper that came with the original code to accommodate the use of the extra parameters our C++ t-sne function requires and to allow the user to choose between the use of the C++ executable or the CPU only scipy implementation of t-sne. The latter allows users without access to the CUDA library, or with hardware not capable of supporting our code, to still run a full t-sne spike sorting session only on CPU and fully in python (no C++ executable is called with this option), albeit for small number of spikes and with large run times.