# Project export: RT-Treat: Real-time neuronal analysis for epilepsy treatment This document was generated by HackStack to give an AI agent context about a hackathon project. Sections are labeled with their provenance; content marked as truncated was cut to keep this document small. ## Project metadata - Hackathon: TreeHacks 2026 - Tagline: RT-Treat is a real-time analysis software system for single neuron activity during surgery allowing for more effective, efficient, and successful treatment of epilepsy. - Devpost: https://devpost.com/software/real-time-neuronal-signal-analysis-for-epilepsy-treatment - GitHub: https://github.com/max-c-lim/RT-Treat - Video: https://www.youtube.com/embed/ysIoLNNV7QI?enablejsapi=1&hl=en_US&rel=0&start=&version=3&wmode=transparent - Team: 0 GitHub contributor(s) — ## Devpost submission (written by the team) ### Inspiration For the past year, I have been working in a hospital's neurosurgery department, focusing on studying and treating patients with intractable drug-resistant epilepsy. To determine the best course of treatment, these patients are implanted with electrodes that reach deep into their brain for multiple weeks, allowing neurologists to determine the seizure onset location by observing their neurons' activity. I have observed that for 1 in every 3 of my patients, the surgery is unsuccessful in determining the best next steps for their treatment, meaning they took the risk of having 12 holes drilled into their skull and 96 microwires implanted into their brain for no reward. During the surgery, the neurosurgeons only have access to the physical location of the electrodes, which provides little information about the epilepsy onset region. I believe that the success rate of the surgery could be much improved if they could also see the activity of individual neurons, allowing them to correct their initial location estimates and better determine where to implant. A more informed implantation area with a higher yield of neuronal activity would also reduce the amount of time the electrodes must remain in the patient's brain for observation and data collection after the surgery. ### What it does RT-Treat is a real-time analysis software that connects to the electrodes implanted in the patients' brain, receiving live streams of hundreds of neuronal signals in parallel. The software automatically analyzes hundreds of neurons' activity in parallel. This includes determining how many neurons are present, where the neurons are most active, and when they send signals to communicate with each other. Then, RT-Treat generates easy-to-understand plots that are continuously updated to provide the most pertinent information succinctly without confusing the neurosurgeons and prolonging the surgery. How I built it With C++, I use a proprietary API to connect to the electrodes and receive their signals in real time. There is one C++ process for each of 96 electrodes, and each process sends its data through TCP to a unique Python process which identifies with sub-millisecond precision when a neuron releases a signal and when. This involves a custom algorithm built from carefully chosen components of many different algorithms published in the neuroscience literature to achieve maximum performance with minimum latency. Each Python process also plots the relevant information when deemed necessary. Challenges I ran into The API to connect the electrodes does not work with Python (which is most familiar to me and easiest to write analysis pipelines on), so I had to use a custom C++ script as a intermediary between the electrodes and my analyses, resulting in 3 asynchronous processes running for each of the total 96 electrodes in parallel. Ensuring that the signal was properly communicated with millisecond latency all the way from electrode to graphed figure was extremely difficult. Additionally, creating the algorithm to decode the raw neural signals (which are just a one-dimensional voltages over time) in real time required a huge amount of planning and debugging. The algorithms and analyses commonly used in the neuroscience literature are designed to be used after all of the data has been collected and often take hours to run, meaning I had to cleverly alter them and write my own to ensure the software could provide feedback to the neurosurgeons during the experiment. It was difficult to decide exactly what should be shown to the neurosurgeons. Too much information would clutter their minds and distract them from the surgery while too little would not be much help. In the end, I decided on a simple layout with the neurons' waveforms as a middle ground. Accomplishments that I'm proud of I am very proud of RT-Treat being able to analyze neural signals and create figures showing the results in real time. It is one of the first software systems in the world (and the first in the United States) to do this. Almost all neurologists and neuroscientists use slow computationally expensive algorithms to analyze their data after the surgery or recording session has long passed. What I learned I learned how to create highly asynchronous and parallel processes that communicate with each other for millisecond response time and how to design new algorithms to analyze neural signals in real time. ### What's next I will show this demo to my research professor and head neurosurgeon to hear their thoughts. If they like it, we will go through the necessary steps (ensuring the software is robust, obtaining IRB approval, etc.) to test RT-Treat on real patients during real surgeries. If that goes well, we will create a company selling this software and real-time analysis software in general to other hospitals and companies making brain computer interfaces. While RT-Treat in its current form is for a particular use case, the general principle of analyzing neural signals in real time can be applied to many other areas of neuroscience, such as non-invasive EEG recordings. Further, all current brain-computer-interface technology uses neural signals averaged across hundreds or even thousands of neurons, resulting in low resolution decoding of user intent. On the contrary, RT-Treat deciphers the responses of individual neurons. Applying this to brain computer interfaces could result in a revolution in the field. ## README (from the GitHub repository) # RT-Treat: Real-time neuronal analysis for epilepsy treatment TreeHacks 2026 RT-Treat is a real-time analysis software system for single neuron activity during surgery allowing for more effective, efficient, and successful treatment of epilepsy. ## Inspiration For the past year, I have been working in a hospital's neurosurgery department, focusing on studying and treating patients with intractable drug-resistant epilepsy. To determine the best course of treatment, these patients are implanted with electrodes that reach deep into their brain for multiple weeks, allowing neurologists to determine the seizure onset location by observing their neurons' activity. I have observed that for 1 in every 3 of my patients, the surgery is unsuccessful in determining the best next steps for their treatment, meaning they took the risk of having 12 holes drilled into their skull and 96 microwires implanted into their brain for no reward. During the surgery, the neurosurgeons only have access to the physical location of the electrodes, which provides little information about the epilepsy onset region. I believe that the success rate of the surgery could be much improved if they could also see the activity of individual neurons, allowing them to correct their initial location estimates and better determine where to implant. A more informed implantation area with a higher yield of neuronal activity would also reduce the amount of time the electrodes must remain in the patient's brain for observation and data collection after the surgery. ## What it does RT-Treat is a real-time analysis software that connects to the electrodes implanted in the patients' brain, receiving live streams of hundreds of neuronal signals in parallel. The software automatically analyzes hundreds of neurons' activity in parallel. This includes determining how many neurons are present, where the neurons are most active, and when they send signals to communicate with each other. Then, RT-Treat generates easy-to-understand plots that are continuously updated to provide the most pertinent information succinctly without confusing the neurosurgeons and prolonging the surgery. ## How I built it With C++, I use a proprietary API to connect to the electrodes and receive their signals in real time. There is one C++ process for each of 96 electrodes, and each process sends its data through TCP to a unique Python process which identifies with sub-millisecond precision when a neuron releases a signal and when. This involves a custom algorithm built from carefully chosen components of many different algorithms published in the neuroscience literature to achieve maximum performance with minimum latency. Each Python process also plots the relevant information when deemed necessary. ## Challenges I ran into The API to connect the electrodes does not work with Python (which is most familiar to me and easiest to write analysis pipelines on), so I had to use a custom C++ script as a intermediary between the electrodes and my analyses, resulting in 3 asynchronous processes running for each of the total 96 electrodes in parallel. Ensuring that the signal was properly communicated with millisecond latency all the way from electrode to graphed figure was extremely difficult. Additionally, creating the algorithm to decode the raw neural signals (which are just a one-dimensional voltages over time) in real time required a huge amount of planning and debugging. The algorithms and analyses commonly used in the neuroscience literature are designed to be used after all of the data has been collected and often take hours to run, meaning I had to cleverly alter them and write my own to ensure the software could provide feedback to the neurosurgeons during the experiment. It was difficult to decide exactly what should be shown to the neurosurgeons. Too much information would clutter their minds and distract them from the surgery while too little would not be much help. In the end, I decided on a simple layout with the neurons' waveforms as a middle ground. ## Accomplishments that I'm proud of I am very proud of RT-Treat being able to analyze neural signals and create figures showing the results in real time. It is one of the first software systems in the world (and the first in the United States) to do this. Almost all neurologists and neuroscientists use slow computationally expensive algorithms to analyze their data after the surgery or recording session has long passed. ## What I learned I learned how to create highly asynchronous and parallel processes that communicate with each other for millisecond response time and how to design new algorithms to analyze neural signals in real time. ## What's next for RT-Treat: Real-time neuronal analysis for epilepsy treatment I will show this demo to my research professor and head neurosurgeon to hear their thoughts. If they like it, we will go through the necessary steps (ensuring the software is robust, obtaining IRB approval, etc.) to test RT-Treat on real patients during real surgeries. If that goes well, we will create a company selling this software and real-time analysis software in general to other hospitals and companies making brain computer interfaces. While RT-Treat in its current form is for a particular use case, the general principle of analyzing neural signals in real time can be applied to many other areas of neuroscience, such as non-invasive EEG recordings. Further, all current brain-computer-interface technology uses neural signals averaged across hundreds or even thousands of neurons, resulting in low resolution decoding of user intent. On the contrary, RT-Treat deciphers the responses of individual neurons. Applying this to brain computer interfaces could result in a revolution in the field. ## Detected evidence (automated analysis) Indexed codebase: 21 recognized source files, 210 KB. - C (language) — detected in the code - C++ (language) — detected in the code - Python (language) — detected in the code ## Codebase structure (from repository index) ### Files (30 of 30) ``` helpers/__init__.py helpers/comparison.py helpers/reader.py neuralynx/__init__.py neuralynx/csc_reader_wrapper.py neuralynx/CSCReader.cpp neuralynx/CSCReader.exp neuralynx/CSCReader.ilk neuralynx/CSCReader.lib neuralynx/CSCReader.obj neuralynx/CSCReader.pdb neuralynx/NetComClient.h neuralynx/NetComClient3_x64.lib neuralynx/NetComClient3.lib neuralynx/Nlx_DataTypes.h neuralynx/vc140.pdb osort/__init__.py osort/core/__init__.py osort/core/detect_signal.py osort/core/detect_spikes.py osort/core/realign_spikes.py osort/core/running_mean_std.py osort/core/sort_chunk.py osort/core/sort_spikes_online.py osort/core/test.py osort/core/utils.py osort/plotting.py osort/run.py osort/test_sorting.ipynb README.md ``` ### Dependencies No dependency index available. ### Recent commits (newest first) - Recreated repo ## Key source files (fetched from GitHub, selected and truncated for size) ### osort/__init__.py ```python """ OSort - Online Spike Sorting Package Python implementation of online spike sorting algorithm """ from .run import setup_params, run_online, run_pseudo __all__ = [ 'setup_params', 'run_online', 'run_pseudo', ] ``` ### neuralynx/csc_reader_wrapper.py ```python import os import subprocess import sys from typing import Optional def _get_default_csc_reader_path() -> str: script_dir = os.path.dirname(os.path.abspath(__file__)) return os.path.join(script_dir, "CSCReader.exe") def start_csc_reader( channel: int, port: int, neuralynx_ip: str, python_ip: str = "127.0.0.1", csc_reader_path: Optional[str] = None, stdout=None, stderr=None, ) -> Optional[subprocess.Popen]: print(f"Launching CSCReader.exe for CSC{channel} (Neuralynx at {neuralynx_ip})...") if csc_reader_path is None: csc_reader_path = _get_default_csc_reader_path() if not os.path.exists(csc_reader_path): print(f"Warning: CSCReader.exe not found at {csc_reader_path}") print("Make sure to start CSCReader.exe manually before data collection begins.") return None try: csc_process = subprocess.Popen( [csc_reader_path, str(channel), str(port), neuralynx_ip, python_ip], stdout=stdout, stderr=stderr ) print(f"CSCReader started with PID {csc_process.pid}") return csc_process except Exception as e: print(f"Failed to launch CSCReader: {e}") return None def stop_csc_reader(csc_process: Optional[subprocess.Popen], timeout_s: float = 2.0) -> None: if csc_process is None: return print("Terminating CSCReader...") csc_process.terminate() try: csc_process.wait(timeout=timeout_s) except subprocess.TimeoutExpired: csc_process.kill() ``` ### neuralynx/CSCReader.cpp ```c++ #include #include "NetComClient.h" #include "Nlx_DataTypes.h" #include #include #include #include #include #include #include #include #include #include #pragma comment(lib, "ws2_32.lib") // Global variable declarations NlxNetCom::NetComClient NetComClient; //the one and only instance of the NetComClient const std::wstring NETCOM_APP_ID(L"CSCReader.cpp"); //id string displayed in Cheetah when connected // Defaults (can be overridden by command line arguments) const int DEFAULT_CHANNEL = 177; const int DEFAULT_PORT = 150; const char* DEFAULT_PYTHON_IP = "127.0.0.1"; // IP where Python server is running const char* DEFAULT_NEURALYNX_IP = "127.0.0.1"; // IP where Neuralynx/Cheetah is running // TCP socket for sending data SOCKET tcpSocket = INVALID_SOCKET; bool socketInitialized = false; // Vector to store callback durations std::vector callbackDurations; // Flag to control program termination volatile bool keepRunning = true; // Console control handler for graceful shutdown BOOL WINAPI ConsoleHandler(DWORD signal) { if (signal == CTRL_C_EVENT || signal == CTRL_BREAK_EVENT || signal == CTRL_CLOSE_EVENT) { std::cout << "\n\nReceived termination signal. Shutting down..." << std::endl; keepRunning = false; return TRUE; } return FALSE; } // Packet structure for sending CSC data #pragma pack(push, 1) struct CSCPacket { std::uint64_t timestamp; // 8 bytes std::uint32_t numSamples; // 4 bytes std::int16_t samples[512]; // 1024 bytes (max samples) }; #pragma pack(pop) // Initialize TCP connection bool InitializeTCP(const char* host, int port) { WSADATA wsaData; // Initialize Winsock int result = WSAStartup(MAKEWORD(2, 2), &wsaData); if (result != 0) { std::cout << "WSAStartup failed: " << result << std::endl; return false; } // Create socket tcpSocket = socket(AF_INET, SOCK_STREAM, IPPROTO_TCP); if (tcpSocket == INVALID_SOCKET) { std::cout << "Socket creation failed: " << WSAGetLastError() << std::endl; WSACleanup(); return false; } // Set up server address sockaddr_in serverAddr; serverAddr.sin_family = AF_INET; serverAddr.sin_port = htons(port); inet_pton(AF_INET, host, &serverAddr.sin_addr); // Connect to server result = connect(tcpSocket, (sockaddr*)&serverAddr, sizeof(serverAddr)); if (result == SOCKET_ERROR) { std::cout << "Connect failed: " << WSAGetLastError() << std::endl; closesocket(tcpSocket); WSACleanup(); return false; } std::cout << "TCP connection established to " << host << ":" << port << std::endl; socketInitialized = true; return true; } // Send data samples over TCP // Returns true on success, false on failure bool SendDataOverTCP(const NlxDataTypes::CRRec& record) { if (!socketInitialized || tcpSocket == INVALID_SOCKET) { std::cout << "Socket not initialized!" << std::endl; return false; } // Pack data into CSCPacket CSCPacket packet; packet.timestamp = record.qwTimeStamp; packet.numSamples = record.dwNumValidSamples; // Copy samples (fill remaining with zeros if less than 512) std::memcpy(packet.samples, record.snSamples, sizeof(record.snSamples)); // Send the entire packet in one call int bytesSent = send(tcpSocket, (const char*)&packet, sizeof(CSCPacket), 0); if (bytesSent == SOCKET_ERROR) { std::cout << "Send failed: " << WSAGetLastError() << std::endl; return false; } else if (bytesSent != sizeof(CSCPacket)) { std::cout << "Warning: Partial send (" << bytesSent << " of " << sizeof(CSCPacket) << " bytes)" << std::endl; return false; } return true; } // Clean up TCP connection void CleanupTCP() { if (tcpSocket != INVALID_SOCKET) { closesocket(tcpSocket); tcpSocket = INVALID_SOCKET; } if (socketInitialized) { WSACleanup(); socketInitialized = false; } } // Callback function for CSC data void CSCCallback(void* myClassPtr, NlxDataTypes::CRRec* records, int numRecords, const wchar_t objectName[]) { auto startTime = std::chrono::high_resolution_clock::now(); // Minimal verbosity - only print every 100 callbacks to avoid filling pipes static int callbackCount = 0; callbackCount++; // if (callbackCount % 100 == 1) { // std::wcout << L"Received " << numRecords << L" records from " << objectName // << L" (callback #" << callbackCount << L")" << std::endl; // } for (int i = 0; i < numRecords; i++) { // Send entire record over TCP if (!SendDataOverTCP(records[i])) { std::cout << "\nFailed to send data over TCP. Shutting down..." << std::endl; keepRunning = false; return; // Exit callback immediately } } auto endTime = std::chrono::high_resolution_clock::now(); double duration = std::chrono::duration(endTime - startTime).count(); callbackDurations.push_back(duration); // if (callbackCount % 100 == 1) { // std::cout << "Callback duration: " << duration << " ms" << std::endl; // } } int main(int argc, char* argv[]) { // Parse command line arguments int channel = DEFAULT_CHANNEL; int port = DEFAULT_PORT; std::string pythonIP = DEFAULT_PYTHON_IP; std::string neuralynxIP = DEFAULT_NEURALYNX_IP; if (argc > 1) { channel = std::atoi(argv[1]); if (channel <= 0) { std::cout << "Invalid channel number. Using default: " << DEFAULT_CHANNEL << std::endl; channel = DEFAULT_CHANNEL; } } if (argc > 2) { port = std::atoi(argv[2]); if (port <= 0 || port > 65535) { s [truncated — 3208 more characters] ``` ### osort/plotting.py ```python from matplotlib.colors import LinearSegmentedColormap import matplotlib.pyplot as plt import numpy as np from scipy.signal import lfilter from scipy.signal.windows import gaussian def get_raster(spike_times, event_times, window): """ Generate a spike raster for given spike times and event times. Parameters: spike_times (array-like): Array of spike times (in seconds). event_times (array-like): Array of event times (in seconds). window (tuple): Time window around each event to consider (start, end) in seconds. Returns: raster (numpy.ndarray of numpy.ndarray): List where each element is an array of spike times aligned to the corresponding event. """ raster = [] for event in event_times: aligned_spikes = spike_times - event relevant_spikes = aligned_spikes[(aligned_spikes >= window[0]) & (aligned_spikes < window[1])] raster.append(relevant_spikes) return np.array(raster, dtype=object) def plot_fr_bar(spike_times, event_times, window, event_colors=None, ax=None): """ Plot a bar graph of firing rates for given spike times and event times. Parameters: spike_times (array-like): Array of spike times (in seconds). event_times (array-like): List of list of event times where each inner list is the event times for a specific identity/category/stim (in seconds). window (tuple): Time window around each event to consider (start, end) in seconds. """ # Calculate firing rates means = [] errs = [] for events in event_times: rates = [] for event in events: aligned_spikes = spike_times - event count = np.sum((aligned_spikes >= window[0]) & (aligned_spikes < window[1])) rate = count / (window[1] - window[0]) rates.append(rate) means.append(np.mean(rates)) errs.append(np.std(rates, ddof=1) / np.sqrt(len(rates))) if ax is None: fig, ax = plt.subplots(1) # Plot the bar graph ax.barh( np.arange(len(event_times)), means, color=event_colors, xerr=errs, error_kw=dict(ecolor='black', capsize=2, capthick=2, lw=2)) ax.set_xlabel('Firing Rate (Hz))') ax.set_ylim(-1, len(event_times)) def plot_spike_raster(spike_times, event_times, window, event_colors=None, ax=None): """ Plot a spike raster for given spike times and event times. Parameters: spike_times (array-like): Array of spike times (in seconds). event_times (array-like): Array of event times (in seconds). window (tuple): Time window around each event to consider (start, end) in seconds. """ # Align spike times to each event raster = get_raster(spike_times, event_times, window) if ax is None: fig, ax = plt.subplots(1) # Plot the raster ax.eventplot(raster, colors=event_colors if event_colors is not None else 'k', linelengths=2) ax.axvline(0, color='k', linestyle='--', label='Event onset', alpha=0.3) ax.set_xlim(window[0], window[1]) ax.set_xlabel('Time (s) relative to event onset') ax.set_ylim(-0.5, len(event_times)-0.5) def get_psth(spike_times, event_times, window, bin_size): """ Compute the Peri-Stimulus Time Histogram (PSTH) for given spike times and event times. Parameters: spike_times (array-like): Array of spike times (in seconds). event_times (array-like): Array of event times (in seconds). window (tuple): Time window around each event to consider (start, end) in seconds. bin_size (float): Size of each bin in seconds. Returns: bin_edges (numpy.ndarray): Edges of the bins. psth (numpy.ndarray): Counts of spikes in each bin. """ # Calculate the number of bins num_bins = round((window[1] - window[0]) / bin_size) # Initialize PSTH raster = np.zeros((len(event_times), num_bins)) # Create bin edges bin_edges = np.linspace(window[0], window[1], num_bins + 1) # Loop over each event time for e, event in enumerate(event_times): # Align spike times to the current event aligned_spikes = spike_times - event # Select spikes within the specified window relevant_spikes = aligned_spikes[(aligned_spikes >= window[0]) & (aligned_spikes < window[1])] # Bin the relevant spikes counts, _ = np.histogram(relevant_spikes, bins=bin_edges) # Accumulate counts into PSTH raster[e] = counts bin_centers = (bin_edges[1:] + bin_edges[:-1]) / 2 psth_mean = np.mean(raster, axis=0) psth_err = np.std(raster, axis=0, ddof=1) / np.sqrt(len(event_times)) return bin_centers, psth_mean, psth_err def plot_psth(spike_times, event_times, window, bin_size, color=None, ax=None, label=None): """ Plot the Peri-Stimulus Time Histogram (PSTH) for given spike times and event times. Parameters: spike_times (array-like): Array of spike times (in seconds). event_times (array-like): Array of event times (in seconds). window (tuple): Time window around each event to consider (start, end) in seconds. bin_size (float): Size of each bin in seconds. """ bin_centers, psth_mean, psth_err = get_psth(spike_times, event_times, window, bin_size) if ax is None: fig, ax = plt.subplots(1) # Plot the PSTH ax.plot(bin_centers, psth_mean, color=color if color is not None else 'k', label=label) ax.fill_between(bin_centers, psth_mean - psth_err, psth_mean + psth_err, color=color if color is not None else 'k', alpha=0.3) ax.axvline(0, color='k', linestyle='--', alpha=0.3) ax.set_xlabel('Time (s) relative to event onset') ax.set_ylabel('Firing Rate (Hz)') ax.set_xlim(window[0], window[1]) def plot_spike_pdf(waveforms, ax=None, times=None, winP=None, doploting=None, bins=None): """ Estimates and optionally plots the probability density function (PDF) of spike amplitudes over time. Args: [truncated — 14986 more characters] ``` ### osort/core/test.py ```python import tkinter as tk import matplotlib.pyplot as plt # Get the screen dimensions root = tk.Tk() screen_width = root.winfo_screenwidth() screen_height = root.winfo_screenheight() root.destroy() print(f"Screen dimensions: {screen_width}x{screen_height}") # Create figure with position at top middle of screen fig = plt.figure(figsize=(8, 6)) # Position: left corner at x=screen_width/2, top edge at y=0 mngr = fig.canvas.manager print(mngr.window.winfo_screenwidth(), mngr.window.winfo_screenheight()) mngr.window.geometry(f"+{int(screen_width/2)}+0") # Create second figure directly to the right fig2 = plt.figure(figsize=(8, 6)) # Calculate width of first figure in pixels (8 inches * dpi) fig_width_pixels = int(8 * fig.dpi) mngr2 = fig2.canvas.manager mngr2.window.geometry(f"+{int(screen_width/2) + fig_width_pixels}+0") # Create third figure directly below fig2 = plt.figure(figsize=(8, 6)) # Calculate width of first figure in pixels (8 inches * dpi) fig_height_pixels = int(6 * fig.dpi) mngr2 = fig2.canvas.manager mngr2.window.geometry(f"+{int(screen_width/2)}+{fig_height_pixels}") # mngr2.window.geometry(f"+{int(screen_width/2) + fig_width_pixels}+{fig_height_pixels}") plt.show() ``` ### osort/core/running_mean_std.py ```python """ Utility class for maintaining running mean and standard deviation Uses Welford's method for numerical stability """ import numpy as np class RunningMeanStd: """ Class to keep a mean and STD of a data stream, updates periodically (so not actually a fully running mean and std) """ def __init__(self, update_every_n_samples: int = 512000): """ Initialize RunningMeanStd Args: update_every_n_samples: Update mean/std every N samples """ self.next_mean = 0.0 self.next_m2 = 0.0 self.next_n = 0.0 self.cur_mean = np.nan self.cur_m2 = np.nan self.cur_n = np.nan self.update_every_n_samples = update_every_n_samples def is_ready(self) -> bool: """Check if mean and std are ready (have been computed at least once)""" return not np.isnan(self.cur_mean) def get_mean(self) -> float: """Get current mean""" return self.cur_mean def get_std(self) -> float: """Get current standard deviation""" return np.sqrt(self.cur_m2 / (self.cur_n - 1)) def update(self, new_data: np.ndarray) -> 'RunningMeanStd': """ Update mean and STD with new data using Welford's method std = sqrt(M2 / (N - 1)) Args: new_data: New data array to incorporate Returns: self for method chaining """ batch_n = len(new_data) batch_mean = np.mean(new_data) batch_m2 = np.sum((new_data - batch_mean) ** 2) delta = batch_mean - self.next_mean total_n = self.next_n + batch_n self.next_mean = self.next_mean + delta * (batch_n / total_n) self.next_m2 = self.next_m2 + batch_m2 + delta**2 * self.next_n * batch_n / total_n self.next_n = total_n if self.next_n >= self.update_every_n_samples: self.cur_mean = self.next_mean self.cur_m2 = self.next_m2 self.cur_n = self.next_n self.next_mean = 0.0 self.next_m2 = 0.0 self.next_n = 0.0 return self ``` ### osort/run.py ```python """ Unified spike sorting runner with online and pseudo-online modes This module combines: - Parameter setup (setup_params) - Pseudo-online spike sorting (run_pseudo) - replay recorded data - Online spike sorting (run_online) - process live Neuralynx data """ from multiprocessing import Process, Queue import os import subprocess import sys from struct import unpack from time import perf_counter, sleep, time from typing import List import matplotlib.pyplot as plt from matplotlib.axes import Axes import numpy as np from scipy.stats import t import socket sys.path.append(R'D:\rutishauser\ClosedLoop') from osort.core.sort_chunk import sort_chunk from neuralynx.csc_reader_wrapper import start_csc_reader, stop_csc_reader def setup_params(cscs: List[int] = None) -> dict: """ Setup parameters for spike sorting Args: cscs: List of CSC channel numbers to process Returns: Dictionary containing all sorting parameters """ params = {} # Sorting params params['running_mean_std_update_n'] = 5120 # 512000 params['detection_method'] = 1 # 1: power, 2: T pos, 3: T min, 4: wavelet params['peak_align_method'] = 1 # 1: find peak, 2: none, 3: power, 4: MTEO params['align_method'] = 3 # 1: pos, 2: neg, 3: mixed params['extraction_threshold'] = 5 params['correction_factor_threshold'] = 0 # minimal threshold, >0 makes threshold bigger params['running_average_length'] = 100 # Use last 100 spikes to compute running average params['merge_clusters'] = True # Whether to merge clusters online params['detection_params'] = {} params['detection_params']['kernel_size'] = 18 # for power method params['sampling_freq'] = 32000 # Hz params['before_peak'] = 24 params['after_peak'] = 59 params['prewhiten'] = 0 params['limit'] = 32556 params['seg_size'] = 512 # Number of frames per segment # Params for detecting: how many previous frames to keep for convolution kernel_size = params['detection_params']['kernel_size'] num_oob_frames = max(kernel_size - 1, params['after_peak']) # Out-of-bounds frames num_prev_frames = num_oob_frames + params['before_peak'] params['num_prev_frames'] = num_prev_frames params['cache_signal_start'] = num_prev_frames - (kernel_size - 1) + 1 # Running params # For saving sorting results params['save_root'] = R'D:\Users\maxlim\epilepsy\P111CS\MaxVarunCLScreen_12152025\sorting' params['save_name_csc_pre_body'] = '%s_preBody.mat' # When sorting in beginning and allowing merging params['save_name_csc_online_body'] = '%s_onlineBody.mat' # When sorting online (no merging) # For saving plots that can be manually looked at during experiment params['manual_root'] = R'D:\Users\maxlim\epilepsy\P111CS\MaxVarunCLScreen_12152025\manual' params['manual_selected_name'] = 'selected' params['manual_rejected_name'] = 'rejected' params['num_trials'] = 1000 # For storing responses params['ms_before_stim'] = 500 params['ms_after_stim'] = 1000 # # Additional setup # if not os.path.exists(params['save_root']): # os.makedirs(params['save_root']) # selected_path = os.path.join(params['manual_root'], params['manual_selected_name']) # if not os.path.exists(selected_path): # os.makedirs(selected_path) # params['selected_path'] = selected_path # rejected_path = os.path.join(params['manual_root'], params['manual_rejected_name']) # if not os.path.exists(rejected_path): # os.makedirs(rejected_path) # params['rejected_path'] = rejected_path # CSC setup if cscs is None: cscs = [1] # Default CSC channel to process params['cscs'] = [f'CSC{csc}' for csc in cscs] params['num_cscs'] = len(params['cscs']) params['csc_to_idx'] = {csc: idx for idx, csc in enumerate(params['cscs'])} return params def run_pseudo(trace, params=None, block_timestamps=None, block_size=512, skip_first_seconds=0, num_blocks_per_chunk=1): """ Run spike sorting in pseudo-online mode by processing data in chunks. Parameters ---------- trace : ndarray Raw signal data (1D array). params : dict Sorting parameters from setup_params(). If None, default parameters (`setup_params`) will be used. block_timestamps : ndarray, optional Timestamp of each 512-sample block in microseconds. If None, timestamps are calculated from sampling frequency. block_size : int, default=512 Number of samples per block (used for timestamp calculation). skip_first_seconds : float, default=0 Skip this many seconds at the start. num_blocks_per_chunk : int, default=1 Number of blocks to process in each iteration. 1 to simulate lowest latency possible. Offline OSort uses 1000 Returns ------- cluster_handles : dict Final cluster handles containing all spike information including: - all_spike_waveforms: All detected spike waveforms - all_spike_ids: Cluster assignments for all spikes - all_spike_timestamps: Timestamps in microseconds for all spikes """ trace = trace.astype(float) # Initialize state variables trace_idx = 0 cluster_handles = {} detect_handles = {} prev_end_trace = None if params is None: params = setup_params() # Sampling frequency in MHz for timestamp calculation samp_freq_mhz = params['sampling_freq'] / 1e6 # Optional: Pre-calculate std estimate from entire recording for perfect offline match std_estimates = [] for i in range(0, len(trace), 512000): window = trace[i:min(len(trace), i + 512000)] std_estimates.append(np.std(window)) params['sorting_std_trace'] = np.mean(std_estimates) # Process data in chunks num_frames_per_chunk = bl [truncated — 17025 more characters] ``` ### osort/core/detect_signal.py ```python """ Process signal to make it suitable for spike detection """ import numpy as np from typing import Tuple, Optional from osort.core.running_mean_std import RunningMeanStd from osort.core.utils import running_std def detect_signal(detect_handles: dict, params: dict, raw_signal: np.ndarray, prev_end_signal: Optional[np.ndarray]) -> Tuple[dict, np.ndarray, np.ndarray, np.ndarray]: """ Process signal (i.e. convolve with power kernel or CNN) to make it suitable for spike detection Args: detect_handles: Dictionary containing detection state params: Parameters dictionary raw_signal: Raw signal data prev_end_signal: Previous iteration's ending signal for continuity Returns: Tuple of (detect_handles, filtered_signal, run_std2, upperlim) """ # NOTE: IF START FILTERING SIGNAL, NEED TO HANDLE HOW prev_end_signal IS USED filtered_signal = raw_signal.copy() # Calculate the to-be-thresholded signal, depending on the method used detection_method = params['detection_method'] if detection_method == 0: # Voltage crossing upperlim = 7.2874 * 4 if prev_end_signal is not None: filtered_signal = np.concatenate([prev_end_signal, filtered_signal]) run_std2 = np.abs(filtered_signal) run_std2[:params['before_peak']] = -np.inf elif detection_method == 1: # Power method kernel_size = params['detection_params']['kernel_size'] # Calculate power if prev_end_signal is None: run_std2 = running_std(filtered_signal, kernel_size) # Pad end with last value run_std2 = np.concatenate([run_std2, np.full(kernel_size - 1, run_std2[-1])]) cache_signal = run_std2 else: filtered_signal = np.concatenate([prev_end_signal, filtered_signal]) run_std2 = running_std(filtered_signal, kernel_size) # Cache signal excludes the padding cache_signal = run_std2[params['cache_signal_start'] - 1:] # -1 for 0-indexing # Pad end with last value run_std2 = np.concatenate([run_std2, np.full(kernel_size - 1, run_std2[-1])]) # Update cache and possibly mean/std if 'running_detected_mean_std' not in detect_handles: detect_handles['running_detected_mean_std'] = RunningMeanStd(params['running_mean_std_update_n']) detect_handles['running_detected_mean_std'].update(cache_signal) # Calculate detection threshold upperlim_fixed = (detect_handles['running_detected_mean_std'].get_mean() + params['extraction_threshold'] * detect_handles['running_detected_mean_std'].get_std()) upperlim = np.ones(len(run_std2)) * upperlim_fixed else: raise ValueError(f'Unknown detection method: {detection_method}') return detect_handles, filtered_signal, run_std2, upperlim ``` ### neuralynx/NetComClient.h ```c //******************************************************************************************************************** // File Name: NetComClient.h // Copyright 1998..2014 @ Neuralynx, Inc. //******************************************************************************************************************** /** \file NetComClient.h * NetCom Client C++ API */ #pragma once #include #include //These forward declares allow us to only distribute this //header file, instead of the entire header tree. The //necessary header files are included in the cpp file class NetworkNodeClient; namespace NlxDataTypes { struct SERec; struct STRec; struct TTRec; struct CRRec; struct EventRec; struct VideoRec; } #define DllExport __declspec(dllexport) /// \summary Namespace for Neuralynx NetCom C++ interface namespace NlxNetCom { /** NetCom Client C++ API */ class NetComClient { public: /** Constructor for the NetComClient class. Creates a new instance of this class, and initializes variables. No actions are performed after this class is created. It is recommended that each application only have a single instance of a NetComClient. */ DllExport NetComClient(void); /** Destructor for the NetComClient class. Deletes the instance of this class. */ DllExport virtual ~NetComClient(void); /** Attempt to make a network connection with the Server * \param [in] serverName Either the network name (i.e. "DataAcqSysPC") or a string representing the IP address (i.e. "192.168.1.100" ) of the PC running a NetCom server. * \param [in] attemptRouterConnection (Default = true) An optional argument that tells the NetComClient to attempt to connect to the Router application prior to attempting a direct connection to the DAS. The DAS only supports a single NetComClient connection, so using the Router is advisable if multiple applications need data from the DAS. This argument defaults to True if not specified. For most applications this value should always be True. * \return Returns True if the connection to serverName was successful. Returns False on a failed connection attempt. If this function is called while currently connected to a NetCom server, it will return False and not close the current connection. */ DllExport bool ConnectToServer(const wchar_t* const serverName, bool attemptRouterConnection = true); /** Disconnects this client to a NetCom server application. * \return Returns True if the disconnection from the currently connected server was successful. Returns False on a failed disconnection attempt. On a failed disconnect, the state of the NetCom client is undetermined. It is advisable that a new NetComClient object be created. If this function is called while not connected to a NetCom server, it will return False. */ DllExport bool DisconnectFromServer(); /** Retrieves the version of the client being used. * \return The version of the NetCom client being used. */ DllExport std::wstring GetClientVersionString() { std::wstring clientVersionString(L""); wchar_t* pClientVersionString = NULL; GetClientVersionStringInternal(pClientVersionString); clientVersionString = pClientVersionString; FreeArray(pClientVersionString); return clientVersionString; } /** Opens a record stream between this client and a NetCom server. Opening a stream will cause the defined callback function corresponding to the object type of objectName to be called after receiving a record for the specified object name. * \param [in] DASObjectName The name of the object to stream data from. This name is specified in the DAS setup files. A listing of defined object names can be obtained from the GetDASObjectsAndTypes function. * \return Returns True if the the specified object name and type were found in the DAS's object list, and a stream was successfully opened. Returns False on any of the following conditions: * -# The object name specified was not found in the DAS's object list. * -# A network or other error prevented the stream from opening. * . * If this function is called while not connected to a NetCom server, it will return False and the stream will need to be reopened after a connection is established. If this function is called successfully multiple times using the same arguments, it must be closed the same number of times to halt record callbacks. Calling DisconnectFromServer automatically closes all opened streams. */ DllExport bool OpenStream(const wchar_t* const DASObjectName); /** Closes a record stream between this client and a NetCom server. Closing a stream will cause the DAS to cease sending records for the specified object. The callback function will continue to be called until all records for this object, received before calling CloseStream, have been processed. * \param [in] DASObjectName The name of the object whose stream should be closed. This name is specified in the DAS setup files. A listing of defined object names can be obtained from the GetDASObjectsAndTypes function. * \return Returns True if the the specified object name and type were found in the DAS's object list, and a stream was successfully closed. Returns false on any of the following conditions: * -# The object name specified was not found in the DAS's object list. * -# A stream for this object name and type has not yet been opened. * -# A network or other error prevented the stream from opening. * . * If this function is called while not connected to a NetCom server, it will return False. If OpenStream is called successfully multiple times using the same arguments, it must be closed the same number of times to halt record callbacks. Calling DisconnectFromServer automatically closes all opened streams. */ DllExport bool CloseStream(const wchar_t* const DASObjectName); /** Sends a command to the NetCom server. * \param [in] command An ASCII command string to send to the server. * \param [out] reply Returned messaged [truncated — 17658 more characters] ``` ### osort/core/sort_chunk.py ```python """ Main spike sorting function """ import numpy as np from typing import Tuple, Optional from osort.core.running_mean_std import RunningMeanStd from osort.core.detect_signal import detect_signal from osort.core.detect_spikes import detect_spikes from osort.core.utils import upsample_spikes from osort.core.realign_spikes import realign_spikes from osort.core.sort_spikes_online import sort_spikes_online def sort_chunk(params: dict, detect_handles: dict, cluster_handles: dict, raw_signal: np.ndarray, prev_end_signal: Optional[np.ndarray]) -> Tuple[dict, dict, np.ndarray, np.ndarray, np.ndarray]: """ Main spike sorting function Args: params: Parameters dictionary detect_handles: Detection state (empty dict {} on first call) cluster_handles: Clustering state (empty dict {} on first call) raw_signal: Raw CSC signal data prev_end_signal: Previous iteration's ending signal for continuity (None on first call) Returns: Tuple of (detect_handles, cluster_handles, spike_waveforms, spike_ids, spike_timestamps) """ spike_waveforms = np.array([]).reshape(0, params['before_peak'] + params['after_peak'] + 1) spike_ids = np.array([]) spike_timestamps = np.array([]) # Calculate STD of raw signal if 'running_raw_mean_std' not in detect_handles: detect_handles['running_raw_mean_std'] = RunningMeanStd(update_every_n_samples=params['running_mean_std_update_n']) detect_handles['running_raw_mean_std'].update(raw_signal) # Get convolved signal that will be used for detecting spikes detect_handles, filtered_signal, run_std2, upperlim = detect_signal( detect_handles, params, raw_signal, prev_end_signal ) # Only proceed if detection is ready if not detect_handles['running_detected_mean_std'].is_ready(): return detect_handles, cluster_handles, spike_waveforms, spike_ids, spike_timestamps std_trace = detect_handles['running_raw_mean_std'].get_std() spike_waveforms, spike_timestamps, detect_handles = detect_spikes( filtered_signal, std_trace, run_std2, upperlim, params, detect_handles ) if spike_waveforms.shape[0] == 0: return detect_handles, cluster_handles, spike_waveforms, spike_ids, spike_timestamps # Get current number of spikes if 'all_spike_ids' in cluster_handles: cur_nr_spikes = len(cluster_handles['all_spike_ids']) else: cur_nr_spikes = 0 # Upsample and realign spikes spike_waveforms = upsample_spikes(spike_waveforms) spike_waveforms, spike_timestamps, _ = realign_spikes( spike_waveforms, spike_timestamps, params['align_method'], std_trace ) # Handle optional sorting std_trace parameter if 'sorting_std_trace' in params: std_trace = params['sorting_std_trace'] # Sort spikes online cluster_handles = sort_spikes_online(params, cluster_handles, spike_waveforms, std_trace) # Get IDs of newly sorted spikes spike_ids = cluster_handles['all_spike_ids'][cur_nr_spikes:] return detect_handles, cluster_handles, spike_waveforms, spike_ids, spike_timestamps ``` [7 more indexed source files omitted to keep this export small. The full file list is in the Codebase structure section above.]