From defd4ec1c8387ff6ea2df2416e675018a7f70159 Mon Sep 17 00:00:00 2001 From: tyler Date: Fri, 31 Jul 2026 01:17:29 -0700 Subject: [PATCH] 1.5.2 time --- changelog.md | 4 + changelog_major.md | 4 + flares.py | 1022 ++++++++++++++++++-------------------------- 3 files changed, 426 insertions(+), 604 deletions(-) diff --git a/changelog.md b/changelog.md index b184e97..0d08356 100644 --- a/changelog.md +++ b/changelog.md @@ -12,6 +12,10 @@ - New parameters have been added to the right side of the screen! This allows for more flexibility and customizability when processing - A new Preference Menu option has been added: Show Advanced Parameters. This keeps some of the parameters hidden when not checked. Since this is a preference, it will be saved when reopening the application - Advanced parameters should only be changed if you know what you are doing, and will have a yellow warning symbol next to them to avoid potential confusion on what parameters are advanced +- Optimized some of the calculations in Scalp Coupling Index to speed up Step 6 by ~25% +- Removed duplicate/redundant calculations in Peak Spectral Power to speed up Step 8 by ~50% +- Changed how the figures are generated when processing to speed up Step 28 by ~85% +- Removed unused methods inside the processing file to slightly speed up application load time - Fixed an issue with the build script not properly updating the version string causing the application to falsely think that an update was always available - Fixed an issue where parameters that were dependent on SHORT_CHANNELS were not properly being updated - Fixed an issue where files that failed processing were not having their progress bar turn red at the step that failed diff --git a/changelog_major.md b/changelog_major.md index a72f492..844566c 100644 --- a/changelog_major.md +++ b/changelog_major.md @@ -12,6 +12,10 @@ - New parameters have been added to the right side of the screen! This allows for more flexibility and customizability when processing - A new Preference Menu option has been added: Show Advanced Parameters. This keeps some of the parameters hidden when not checked. Since this is a preference, it will be saved when reopening the application - Advanced parameters should only be changed if you know what you are doing, and will have a yellow warning symbol next to them to avoid potential confusion on what parameters are advanced +- Optimized some of the calculations in Scalp Coupling Index to speed up Step 6 by ~25% +- Removed duplicate/redundant calculations in Peak Spectral Power to speed up Step 8 by ~50% +- Changed how the figures are generated when processing to speed up Step 28 by ~85% +- Removed unused methods inside the processing file to slightly speed up application load time - Fixed an issue with the build script not properly updating the version string causing the application to falsely think that an update was always available - Fixed an issue where parameters that were dependent on SHORT_CHANNELS were not properly being updated - Fixed an issue where files that failed processing were not having their progress bar turn red at the step that failed diff --git a/flares.py b/flares.py index 405a3e7..961a614 100644 --- a/flares.py +++ b/flares.py @@ -14,28 +14,26 @@ import sys import json import time import logging -import platform import warnings import threading import traceback import itertools import os.path as op from io import BytesIO -from queue import Empty from pathlib import Path from copy import deepcopy import multiprocessing as mp from itertools import compress -from multiprocessing import Queue +from queue import Empty, Queue from typing import Any, Optional, Sequence, cast, Literal, Union # External library imports import matplotlib.pyplot as plt import matplotlib.colors as mcolors -from matplotlib.figure import Figure from matplotlib.axes import Axes -from matplotlib.colors import LinearSegmentedColormap from matplotlib.lines import Line2D +from matplotlib.figure import Figure +from matplotlib.colors import LinearSegmentedColormap import numpy as np from numpy.typing import NDArray @@ -44,8 +42,8 @@ from numpy import float64, floating import pandas as pd from pandas import DataFrame -import seaborn as sns import h5py +import seaborn as sns from nilearn.plotting import plot_design_matrix # type: ignore from nilearn.glm.regression import OLSModel @@ -55,7 +53,7 @@ from statsmodels.stats.multitest import multipletests from statsmodels.tools.sm_exceptions import ConvergenceWarning from scipy.spatial.distance import cdist -from scipy.signal import welch, butter, filtfilt # type: ignore +from scipy.signal import welch, butter, filtfilt, periodogram # type: ignore from scipy.stats import pearsonr, zscore, ttest_1samp, ttest_ind, sem import pywt # type: ignore @@ -96,10 +94,13 @@ from mne_nirs.signal_enhancement import ( ) # type: ignore from mne_nirs.io.fold import fold_channel_specificity # type: ignore from mne_nirs.preprocessing import peak_power # type: ignore +from mne.preprocessing.nirs.nirs import _validate_nirs_info from mne_nirs.statistics._glm_level_first import RegressionResults # type: ignore from mne_connectivity.viz import plot_connectivity_circle -from mne_connectivity import envelope_correlation, spectral_connectivity_epochs, spectral_connectivity_time +from mne_connectivity import envelope_correlation, spectral_connectivity_epochs, spectral_connectivity_time + +from src.shared.shareddata import PLATFORM_NAME, resource_path @@ -124,7 +125,6 @@ def get_category_color(label): """Returns the primary color if it's a single failure, otherwise gray.""" return PRIMARY_COLORS.get(label, COMBINATION_COLOR) - DOWNSAMPLE: bool DOWNSAMPLE_FREQUENCY: int @@ -263,87 +263,11 @@ GROUP: str = "Default" FOLDING_BYP: bool = False -# These are parameters that are required for the analysis -REQUIRED_KEYS: dict[str, Any] = { - - # "SECONDS_TO_STRIP": int, - "DOWNSAMPLE": bool, - "DOWNSAMPLE_FREQUENCY": int, - - "TRIM": bool, - "SECONDS_TO_KEEP": float, - - "OPTODE_PLACEMENT": bool, - - "HEART_RATE": bool, - - "SCI": bool, - "SCI_TIME_WINDOW": int, - "SCI_THRESHOLD": float, - - "SNR": bool, - # SNR_TIME_WINDOW : int - "SNR_THRESHOLD": float, - - "PSP": bool, - "PSP_TIME_WINDOW": int, - "PSP_THRESHOLD": float, - - "SHORT_CHANNELS": bool, - "SHORT_CHANNELS_THRESHOLD": float, - "LONG_CHANNELS_THRESHOLD": float, - - - "REMOVE_EVENTS": list, - "L_FREQ": float, - "H_FREQ": float, - - "EPOCH_HANDLING": str, - "T_MIN": int, - "T_MAX": int, - - "TDDR": bool, - "WAVELET": bool, - "IQR": float, - "WAVELET_TYPE": str, - "WAVELET_LEVEL": int, - "FILTER": bool, - "DRIFT_MODEL": str, - # "REJECT_PAIRS": bool, - # "FORCE_DROP_ANNOTATIONS": list, - # "FILTER_LOW_PASS": float, - # "FILTER_HIGH_PASS": float, - # "EPOCH_PAIR_TOLERANCE_WINDOW": int, -} - - -# audit_log = logging.getLogger("memory_audit") -# audit_log.setLevel(logging.INFO) -# audit_log.propagate = False # This prevents it from talking to other loggers - -# # 2. Add a file handler specifically for this audit logger -# if not audit_log.handlers: -# fh = logging.FileHandler('flares_memory_audit.log') -# fh.setFormatter(logging.Formatter('%(asctime)s | PID: %(process)d | %(message)s')) -# audit_log.addHandler(fh) - -# def get_mem_mb(): -# return psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 - - - -class ProcessingError(Exception): - def __init__(self, message: str = "Something went wrong!"): - self.message = message - super().__init__(self.message) - # Ensure that we are working in the directory of this file script_dir = os.path.dirname(os.path.abspath(__file__)) os.chdir(script_dir) -PLATFORM_NAME = platform.system().lower() - # Configure logging to file with timestamps and realtime flush if PLATFORM_NAME == 'darwin': logging.basicConfig( @@ -367,7 +291,6 @@ logger = logging.getLogger() - def set_config_me(config: dict[str, Any]) -> None: """ Validates and applies the given configuration dictionary. @@ -381,6 +304,7 @@ def set_config_me(config: dict[str, Any]) -> None: globals().update(config) + def set_metadata(file_path, metadata: dict[str, Any]) -> None: """ @@ -406,7 +330,9 @@ def set_metadata(file_path, metadata: dict[str, Any]) -> None: if val not in (None, '', [], {}, ()): # check for "empty" values globals()[key] = val -def gui_entry(config: dict[str, Any], gui_queue: Queue, progress_queue: Queue, ack_queue: Queue) -> None: + + +def gui_entry(config: dict[str, Any], gui_queue: mp.Queue, progress_queue: mp.Queue, ack_queue: mp.Queue) -> None: start_time = time.time() try: file_paths = config['SNIRF_FILES'] @@ -472,7 +398,7 @@ def process_participant_worker(file_path, file_params, file_metadata, result_que except Exception: pass - result = process_participant(file_path, progress_callback=progress_callback) + result = process_participant(file_path, file_start, progress_callback=progress_callback) duration = time.time() - file_start result_queue.put((file_path, result, None, duration, 1.0)) @@ -491,6 +417,7 @@ def process_participant_worker(file_path, file_params, file_metadata, result_que pass + def process_multiple_participants(file_paths, file_params, file_metadata, progress_queue=None, gui_queue=None, max_workers=6): ctx = mp.get_context("spawn") @@ -781,6 +708,7 @@ def plot_timechannel_quality_metrics(data, scores, times: list[tuple[float]], co return fig1, fig2 + def scalp_coupling_index_windowed_raw(data, time_window: float = 3.0, l_freq: float = 0.7, h_freq: float = 1.5, l_trans_bandwidth: float = 0.3, h_trans_bandwidth: float = 0.3): """ Compute windowed scalp coupling index (SCI) across fNIRS channels. @@ -809,8 +737,14 @@ def scalp_coupling_index_windowed_raw(data, time_window: float = 3.0, l_freq: fl """ # Pick only fNIRS channels and sort them by channel name - picks: NDArray[np.intp] = pick_types(cast(Info, data.info), fnirs=True) # type: ignore - picks = picks[np.argsort([getattr(data, "ch_names")[pick] for pick in picks])] + sfreq: float = data.info["sfreq"] + ch_names = data.ch_names + times_arr = data.times + + # 1. Pick and sort fNIRS channel indices + picks: NDArray[np.intp] = pick_types(cast(Info, data.info), fnirs=True) + sort_idx = np.argsort([ch_names[p] for p in picks]) + picks = picks[sort_idx] # FIXME: This may happen if the heart rate calculation tries to set a value way too low if l_freq < 0.3: @@ -829,41 +763,51 @@ def scalp_coupling_index_windowed_raw(data, time_window: float = 3.0, l_freq: fl ) # Calculate number of samples per time window, the total number of windows, and prepare output variables - window_samples = int(np.ceil(time_window * getattr(data, "info")["sfreq"])) - n_windows = int(np.floor(len(data) / window_samples)) + window_samples = int(np.ceil(time_window * sfreq)) + n_windows = int(np.floor(filtered_data.shape[1] / window_samples)) + total_samples = n_windows * window_samples + + starts = np.arange(n_windows) * window_samples + stops = np.minimum(starts + window_samples, len(times_arr) - 1) + times = [(times_arr[s], times_arr[e]) for s, e in zip(starts, stops)] + + # 5. Vectorized Correlation Calculation + # Truncate to exact window boundary: shape (n_channels, n_windows, window_samples) + truncated = filtered_data[picks, :total_samples].reshape(len(picks), n_windows, window_samples) + + # Pair channels: c1 (even rows), c2 (odd rows) -> shape (n_pairs, n_windows, window_samples) + c1 = truncated[0::2] + c2 = truncated[1::2] + + # Zero-mean center along the window sample axis + c1_zero = c1 - np.mean(c1, axis=-1, keepdims=True) + c2_zero = c2 - np.mean(c2, axis=-1, keepdims=True) + + # Standard deviations along window sample axis + std1 = np.std(c1, axis=-1) + std2 = np.std(c2, axis=-1) + + # Covariance along window sample axis + cov = np.mean(c1_zero * c2_zero, axis=-1) + denom = std1 * std2 + + # Vectorized Pearson r calculation with NaN/Zero-std handling + with np.errstate(divide="ignore", invalid="ignore"): + corrs = np.where((denom == 0) | np.isnan(denom), 0.0, cov / denom) + + # Assign pair correlations back to output score matrix scores = np.zeros((len(picks), n_windows)) - times: list[tuple[float, float]] = [] + scores[0::2, :] = corrs + scores[1::2, :] = corrs - # Slide through the data in windows to compute scalp coupling index (SCI) - for window in range(n_windows): - start_sample = int(window * window_samples) - end_sample = start_sample + window_samples - end_sample = np.min([end_sample, len(data) - 1]) - - # Track time boundaries for each window - t_start = getattr(data, "times")[start_sample] - t_stop = getattr(data, "times")[end_sample] - times.append((t_start, t_stop)) - - # Iterate through channels in pairs (hbo, hbr). This requires them to be sorted by channel name - for ii in range(0, len(picks), 2): - c1 = filtered_data[picks[ii]][start_sample:end_sample] - c2 = filtered_data[picks[ii + 1]][start_sample:end_sample] - - # Ensure the correlation data is valid - if np.std(c1) == 0 or np.std(c2) == 0 or np.any(np.isnan(c1)) or np.any(np.isnan(c2)): - c = 0 - else: - c = np.corrcoef(c1, c2)[0][1] - - # Assign the computed correlation to both channels in the pair - scores[ii, window] = c - scores[ii + 1, window] = c - - scores = scores[np.argsort(picks)] + # Revert scores to the original pick ordering if needed + inv_sort = np.argsort(sort_idx) + scores = scores[inv_sort] return data, scores, times + + def calculate_scalp_coupling(data, l_freq: float = 0.7, h_freq: float = 1.5, time_window: int = 3, threshold: float = 0.6): """ Calculate the scalp coupling index (SCI) and identify bad channels based on a threshold. @@ -969,6 +913,7 @@ def get_hbo_hbr_picks(raw): return hbo_picks, hbr_picks, hbo_wl, hbr_wl + def interpolate_fNIRS_bads_weighted_average(raw, max_dist=0.03, min_neighbors=2, short_channels_threshold=0.015): """ Interpolate bad fNIRS channels using a distance-weighted average of nearby good channels. @@ -1157,7 +1102,7 @@ def calculate_signal_noise_ratio(data): print("Calculating signal to noise ratio...") # Compute the signal-to-noise ratio values - print("Computing the signal to noise power...") + print("Computing the signal to noise ratio...") signal_band=(0.01, 0.5) noise_band=(1.0, 10.0) data_signal = data.copy().filter(*signal_band, verbose=False) #type: ignore @@ -1203,7 +1148,6 @@ def calculate_signal_noise_ratio(data): - def calculate_peak_power(data: BaseRaw, time_window: int = 3, threshold: float = 0.1, l_freq: float = 0.7, h_freq: float = 1.5) -> tuple[list[str], Figure, Figure]: """ Calculate peak spectral power (PSP) for fNIRS channels and identify bad channels. @@ -1227,7 +1171,7 @@ def calculate_peak_power(data: BaseRaw, time_window: int = 3, threshold: float = # Compute the PSP - _, scores, times = cast(tuple[NDArray[float64], NDArray[float64], list[tuple[float]]], peak_power(data, time_window=time_window, threshold=threshold, l_freq=l_freq, h_freq=h_freq)) + _, scores, times = cast(tuple[NDArray[float64], NDArray[float64], list[tuple[float]]], peak_power_fast(data, time_window=time_window, threshold=threshold, l_freq=l_freq, h_freq=h_freq)) # Identify channels that don't meet the provided threshold psp = scores.mean(axis=1) @@ -1241,6 +1185,7 @@ def calculate_peak_power(data: BaseRaw, time_window: int = 3, threshold: float = return list(compress(cast(list[str], getattr(data, "ch_names")), psp < threshold)), psp1, psp2 + def mark_bads(raw, bad_sci, bad_snr, bad_psp, bad_coeff_var, bad_range, bad_noise, bad_disp): print(bad_sci, bad_snr, bad_psp, bad_coeff_var, bad_range, bad_noise, bad_disp) bads_combined = list(set(bad_snr) | set(bad_sci) | set(bad_psp) | set(bad_coeff_var) | set(bad_range) | set(bad_noise) | set(bad_disp)) @@ -1314,7 +1259,6 @@ def mark_bads(raw, bad_sci, bad_snr, bad_psp, bad_coeff_var, bad_range, bad_noi - def filter_the_data( raw_haemo, filter_algorithm, @@ -1322,21 +1266,21 @@ def filter_the_data( h_freq, l_trans_bandwidth, h_trans_bandwidth, - iir_type, - iir_order, + # iir_type, + # iir_order, filter_length, filter_phase, fir_window, fir_design, - iir_output, - passband_ripple, - stopband_attenuation, + # iir_output, + # passband_ripple, + # stopband_attenuation, filter_pad, - skip_by_annotation, + # skip_by_annotation, filter_n_jobs, verbosity ): - # --- STEP 5: Filtering (0.01-0.2 Hz bandpass) --- + fig_filter = raw_haemo.compute_psd(fmax=3).plot( average=True, color="r", show=False, amplitude=True ) @@ -1360,9 +1304,6 @@ def filter_the_data( raw_haemo = raw_haemo.filter(l_freq=l_freq, h_freq=h_freq, filter_length=filter_length, l_trans_bandwidth=l_trans_bandwidth, h_trans_bandwidth=h_trans_bandwidth, n_jobs=filter_n_jobs, method=filter_algorithm, pad=filter_pad, verbose=verbosity) else: print("No filter") - #raw_haemo = raw_haemo.filter(l_freq=None, h_freq=0.4, h_trans_bandwidth=0.2) - #raw_haemo = raw_haemo.filter(l_freq=None, h_freq=0.7, h_trans_bandwidth=0.2) - #raw_haemo = raw_haemo.filter(0.005, 0.7, h_trans_bandwidth=0.02, l_trans_bandwidth=0.002) raw_haemo.compute_psd(fmax=3).plot( average=True, axes=fig_filter.axes, color="g", amplitude=True, show=False @@ -1439,7 +1380,6 @@ def epochs_calculations(raw_haemo, events, event_dict, epoch_handling, max_shift fig_epochs.append(("fig_epochs_dropped", fig_epochs_dropped)) # Plot for each condition - fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 4)) for idx, condition in enumerate(epochs.event_id.keys()): logger.info(condition) logger.info(idx) @@ -1580,24 +1520,24 @@ def epochs_calculations(raw_haemo, events, event_dict, epoch_handling, max_shift def make_design_matrix( - raw_haemo, - resample, - resample_freq, - stim_dur, - hrf_model, - drift_model, - high_pass, - drift_order, - fir_delays, - min_onset, - oversampling, - short_channel_regression, - short_channels, - long_channels, - short_channels_threshold, - long_channels_threshold, - folding_bypass - ): + raw_haemo, + resample, + resample_freq, + stim_dur, + hrf_model, + drift_model, + high_pass, + drift_order, + fir_delays, + min_onset, + oversampling, + short_channel_regression, + short_channels, + long_channels, + short_channels_threshold, + long_channels_threshold, + folding_bypass +): # events_to_remove = REMOVE_EVENTS events_to_remove = "" @@ -1613,7 +1553,7 @@ def make_design_matrix( if short_channels: short_chans = get_short_channels(raw_haemo, max_dist=short_channels_threshold) if long_channels: - raw_haemo = get_long_channels(raw_haemo, min_dist=long_channels_threshold, max_dist=long_channels_threshold) + raw_haemo = get_long_channels(raw_haemo, min_dist=short_channels_threshold, max_dist=long_channels_threshold) else: short_chans = None @@ -1713,6 +1653,7 @@ def generate_montage_locations(): return coords.reset_index(drop=True) + def _find_closest_standard_location(position, reference, *, out="label"): """Return closest montage label to coordinates. @@ -1837,6 +1778,7 @@ def _read_fold_xls(fname, atlas="Juelich"): return tbl + def _check_load_fold(fold_files, atlas): # _validate_type(fold_files, (list, "path-like", None), "fold_files") if fold_files is None: @@ -1866,176 +1808,6 @@ def _check_load_fold(fold_files, atlas): ) return fold_tbl - - -def resource_path(relative_path): - """ - Get absolute path to resource regardless of running directly or packaged using PyInstaller - """ - - if hasattr(sys, '_MEIPASS'): - # PyInstaller bundle path - base_path = sys._MEIPASS - else: - base_path = os.path.abspath(".") - - return os.path.join(base_path, relative_path) - - - -# def fold_channels(raw: BaseRaw) -> None: - -# # Locate the fOLD excel files -# if getattr(sys, 'frozen', False): -# set_config('MNE_NIRS_FOLD_PATH', resource_path("./mne_data/fOLD/fOLD-public-master/Supplementary")) # type: ignore -# else: -# path = os.path.expanduser("~") + "/mne_data/fOLD/fOLD-public-master/Supplementary" -# set_config('MNE_NIRS_FOLD_PATH', resource_path(path)) # type: ignore - -# output = None - -# # List to store the results -# landmark_specificity_data: list[dict[str, Any]] = [] - -# # Filter the data to only what we want -# hbo_channel_names = cast(list[str], getattr(raw.copy().pick(picks='hbo'), "ch_names")) # type: ignore - -# # Format the output to make it slightly easier to read - -# if True: -# num_channels = len(hbo_channel_names) -# rows, cols = 4, 7 # 6 rows and 4 columns of pie charts -# fig, axes = plt.subplots(rows, cols, figsize=(16, 10), constrained_layout=True) -# axes = axes.flatten() # Flatten the axes array for easier indexing - -# # If more pie charts than subplots, create extra subplots -# if num_channels > rows * cols: -# fig, axes = plt.subplots((num_channels // cols) + 1, cols, figsize=(16, 10), constrained_layout=True) -# axes = axes.flatten() - -# # Create a list for consistent color mapping -# landmarks = [ -# "1 - Primary Somatosensory Cortex", -# "2 - Primary Somatosensory Cortex", -# "3 - Primary Somatosensory Cortex", -# "4 - Primary Motor Cortex", -# "5 - Somatosensory Association Cortex", -# "6 - Pre-Motor and Supplementary Motor Cortex", -# "7 - Somatosensory Association Cortex", -# "8 - Includes Frontal eye fields", -# "9 - Dorsolateral prefrontal cortex", -# "10 - Frontopolar area", -# "11 - Orbitofrontal area", -# "17 - Primary Visual Cortex (V1)", -# "18 - Visual Association Cortex (V2)", -# "19 - V3", -# "20 - Inferior Temporal gyrus", -# "21 - Middle Temporal gyrus", -# "22 - Superior Temporal Gyrus", -# "23 - Ventral Posterior cingulate cortex", -# "24 - Ventral Anterior cingulate cortex", -# "25 - Subgenual cortex", -# "32 - Dorsal anterior cingulate cortex", -# "37 - Fusiform gyrus", -# "38 - Temporopolar area", -# "39 - Angular gyrus, part of Wernicke's area", -# "40 - Supramarginal gyrus part of Wernicke's area", -# "41 - Primary and Auditory Association Cortex", -# "42 - Primary and Auditory Association Cortex", -# "43 - Subcentral area", -# "44 - pars opercularis, part of Broca's area", -# "45 - pars triangularis Broca's area", -# "46 - Dorsolateral prefrontal cortex", -# "47 - Inferior prefrontal gyrus", -# "48 - Retrosubicular area", -# "Brain_Outside", -# ] - -# cmap1 = plt.get_cmap('tab20') # First 20 colors -# cmap2 = plt.get_cmap('tab20b') # Next 20 colors - -# # Combine the colors from both colormaps -# colors = [cmap1(i) for i in range(20)] + [cmap2(i) for i in range(20)] # Total 40 colors - -# landmarks.sort(key=lambda x: (int(x.split(" - ")[0]) if x.split(" - ")[0].isdigit() else float('inf'))) - -# landmark_color_map = {landmark: colors[i % len(colors)] for i, landmark in enumerate(landmarks)} - -# # Iterate over each channel -# print(len(hbo_channel_names)) - -# for idx, channel_name in enumerate(hbo_channel_names): - -# print(idx, channel_name) -# # Run the fOLD on the selected channel -# channel_data = raw.copy().pick(picks=channel_name) # type: ignore - -# output = cast(list[DataFrame], fold_channel_specificity_normal(channel_data, interpolate=True, atlas='Brodmann')) - -# # Process each DataFrame that fold_channel_specificity returns -# for df_data in output: - -# # Extract the relevant columns -# useful_data = df_data[['Landmark', 'Specificity']] - -# # Store the results -# landmark_specificity_data.append({ -# 'Channel': channel_name, -# 'Data': useful_data, -# }) - - -# # Plot the results -# # TODO: Fix this -# if True: -# unique_landmarks = sorted(useful_data['Landmark'].unique()) -# color_list = [landmark_color_map[landmark] for landmark in useful_data['Landmark']] - -# # Plot specificity for each channel -# ax = axes[idx] - -# labels = [f'{landmark.split(" - ")[0]}' if landmark != 'Brain_Outside' else 'B' for landmark in useful_data['Landmark']] - -# wedges, texts, autotexts = ax.pie( -# useful_data['Specificity'], -# autopct='%1.1f%%', -# startangle=90, -# labels=labels, -# labeldistance=1.05, -# colors=color_list) - -# ax.set_title(f'{channel_name}') -# ax.axis('equal') - -# landmark_specificity_data = [] - - -# # TODO: Fix this -# if True: -# handles = [ -# plt.Line2D([0], [0], marker='o', color='w', label=landmark, markersize=10, -# markerfacecolor=landmark_color_map[landmark]) -# for landmark in landmarks -# ] -# n_landmarks = len(landmarks) - -# # Calculate the figure size based on number of rows and columns -# fig_width = 5 -# fig_height = n_landmarks / 4 - -# # Create a new figure window for the legend -# legend_fig = plt.figure(figsize=(fig_width, fig_height)) -# legend_axes = legend_fig.add_subplot(111) -# legend_axes.axis('off') # Turn off axis for the legend window -# legend_axes.legend(handles=handles, loc='center', fontsize=10, title="Landmarks") - -# for ax in axes[len(hbo_channel_names):]: -# ax.axis('off') - -# #plt.show() -# fig_dict = {"main": fig, "legend": legend_fig} -# return convert_fig_dict_to_png_bytes(fig_dict) - def fold_channels(raw: BaseRaw, p_name: str, atlas: str='Brodmann', progress_queue=None) -> dict[str, list[dict[str, Any]]]: @@ -2083,7 +1855,6 @@ def fold_channels(raw: BaseRaw, p_name: str, atlas: str='Brodmann', progress_que - def plot_glm_results(file_path, raw_haemo, glm_est, design_matrix): fig_glms = [] # List to store figures @@ -2179,6 +1950,7 @@ def plot_glm_results(file_path, raw_haemo, glm_est, design_matrix): return fig_glms + def plot_3d_evoked_array( inst: Union[BaseRaw, EvokedArray, Info], statsmodel_df: DataFrame, @@ -2267,6 +2039,7 @@ def plot_3d_evoked_array( return brain + def aggregate_fnirs_group_geometry(raw_list: Sequence[BaseRaw | None]) -> BaseRaw: """ Averages fNIRS geometry across participants in two tiers: @@ -2395,7 +2168,6 @@ def brain_3d_visualization( - def brain_landmarks_3d(raw_haemo: BaseRaw, show_optodes: Literal['sensors', 'labels', 'none', 'all'] = 'all', show_brodmann: bool = True) -> None: brain = Brain("fsaverage", background="white", size=(800, 700)) # type: ignore @@ -2470,6 +2242,7 @@ def brain_landmarks_3d(raw_haemo: BaseRaw, show_optodes: Literal['sensors', 'lab if show_brodmann:# Add Brodmann labels labels = cast(list[Label], read_labels_from_annot("fsaverage", "PALS_B12_Brodmann", "lh", verbose=False)) # type: ignore + #TODO: This has been hardcoded here for the entire applications lifecycle. About time to user expose? label_colors = { "Brodmann.1-lh": "red", "Brodmann.2-lh": "red", @@ -2500,17 +2273,6 @@ def brain_landmarks_3d(raw_haemo: BaseRaw, show_optodes: Literal['sensors', 'lab return brain -def convert_fig_dict_to_png_bytes(fig_dict: dict[str, Figure]) -> dict[str, bytes]: - png_dict = {} - for label, fig in fig_dict.items(): - buf = BytesIO() - fig.savefig(buf, format="png", bbox_inches="tight") - buf.seek(0) - png_dict[label] = buf.read() - plt.close(fig) - return png_dict - - def brain_3d_contrast(con_model_df: DataFrame, con_model_df_filtered: BaseRaw, common_channels: list[str], first_name: str, second_name: str, t_or_theta: Literal['t', 'theta'] = 'theta', show_optodes: Literal['sensors', 'labels', 'none', 'all'] = 'all', show_text: bool = True, brain_bounds: float = 1.0) -> None: # Filter DataFrame to only common channels, and sort by raw order @@ -2551,7 +2313,6 @@ def plot_2d_3d_contrasts_between_groups( brain_bounds: float = 1.0, ) -> None: - logger.info("-----") contrast_df_a = contrast_df_a.copy() contrast_df_a["group"] = group_a_name @@ -2664,7 +2425,6 @@ def plot_2d_3d_contrasts_between_groups( - def plot_fir_model_results( df: DataFrame, raw_haemo: BaseRaw | None, @@ -2674,7 +2434,6 @@ def plot_fir_model_results( u_bound: float, ) -> None: - df["isActivity"] = [f"{selected_event}" in n for n in df["Condition"]] df["isDelay"] = ["delay" in n for n in df["Condition"]] df = df.query("isDelay in [True]") @@ -2848,14 +2607,6 @@ def load_snirf(file_path: str, downsample_frequency: int, verbosity: bool) -> tu # TODO: Why was this commented again? # Maybe this should be a bypass parameter? - # # Strip the specified amount of seconds from the start of the file - # total_duration = getattr(raw, "times")[-1] - # if total_duration > SECONDS_TO_STRIP: - # raw.crop(tmin=SECONDS_TO_STRIP, tmax=total_duration, verbose=VERBOSITY) # type: ignore - # logger.info(f"Stripped first {SECONDS_TO_STRIP} second(s) of data.") - # else: - # logger.info(f"Data length ({total_duration:.2f}s) less than strip duration; no cropping applied.") - # If the user forcibly dropped channels, remove them now before any processing occurs # logger.info("Checking if there are channels to forcibly drop...") # if drop_prefixes: @@ -2879,18 +2630,6 @@ def load_snirf(file_path: str, downsample_frequency: int, verbosity: bool) -> tu - - - - - - -# Data science requires good statistics. This section is requiring work. Am working on figuring this out. Will likely add this to a new viewer. - - - - - def run_roi_second_level_analysis( df_roi_all: DataFrame, df_cha_all: DataFrame | None = None, @@ -3144,7 +2883,6 @@ def run_roi_second_level_analysis( - def clean_subject_id(path_or_id): """ Cleans file paths and ID strings to get a standardized subject identifier. @@ -3925,7 +3663,6 @@ def run_cross_group_contrast_analysis( - def run_roi_paired_contrast_analysis( df_roi_all: DataFrame, roi_pairs: Sequence[tuple[str, str]] | list[list[str]], @@ -4260,10 +3997,6 @@ def aggregate_channel_contrasts_to_roi( roi_theta = roi_theta.rename(columns={'contrast_name': 'Condition'}) return roi_theta[['ROI', 'Condition', 'Chroma', 'theta', 'ID']] - - - - def collapse_fir_condition_column(df, value_col, condition_col='Condition', @@ -4362,8 +4095,9 @@ def _channel_midpoint(ch_info): src = loc[3:6] det = loc[6:9] return tuple((s + d) / 2.0 for s, d in zip(src, det)) - - + + + def _build_axis_split_rois(raw_haemo, axis, names, balance_threshold=0.5): """ Split channels into two ROIs by the sign of one coordinate axis of each @@ -4417,8 +4151,9 @@ def _build_axis_split_rois(raw_haemo, axis, names, balance_threshold=0.5): return None return {names[0]: neg_indices, names[1]: pos_indices} - - + + + def _build_geometric_fallback_rois(raw_haemo): """ Generalized, zero-configuration fallback: try a Left/Right split first @@ -4447,8 +4182,9 @@ def _build_geometric_fallback_rois(raw_haemo): logger.warning("Neither Left/Right nor Front/Back split is usable for this montage.") return None - - + + + def _build_per_channel_rois(raw_haemo): """ Last-resort failsafe: one ROI per physical channel (source-detector @@ -4473,16 +4209,6 @@ def _build_per_channel_rois(raw_haemo): - - - - - - - - - - def calculate_dpf(file_path): # order is hbo / hbr with h5py.File(file_path, 'r') as f: @@ -4754,7 +4480,8 @@ def calculate_heart_rate_scipy(sfreq: float, signal_trimmed: NDArray[float64], s logger.info("Successfully calculated heart rate using SciPy.") return freq_bpm_scipy, psd_scipy, freq_range_scipy, mean_hr_scipy - + + def plot_heart_rate( freq_bpm_scipy: NDArray[floating[Any]], @@ -4841,169 +4568,6 @@ def plot_heart_rate( return fig1, fig2 -# def mark_bads_by_db_threshold(raw, db_limit=-60): -# """ -# Converts a dB threshold to absolute power and marks channels -# exceeding it at 6.25 Hz (1/2 Nyquist). -# """ -# # 1. Convert your "Magic Number" from dB to absolute units -# abs_threshold = 10 ** (db_limit / 10) - -# # 2. Compute PSD specifically around the frequency of interest -# sfreq = raw.info['sfreq'] -# target_freq = (sfreq / 2) - 1 # 6.25 Hz - -# # We use a small window around 6.25 Hz to get a stable average -# spectrum = raw.compute_psd(fmin=target_freq - 0.5, fmax=target_freq + 0.5) -# psd_data, freqs = spectrum.get_data(return_freqs=True) - -# # 3. Get the power at that frequency (averaging across the small window) -# # psd_data shape: (n_channels, n_freqs) -# power_at_6Hz = np.mean(psd_data, axis=1) - -# # 4. Identify the "loud" channels -# bad_indices = np.where(power_at_6Hz > abs_threshold)[0] -# new_bads = [raw.ch_names[i] for i in bad_indices] - -# # 5. Update raw.info['bads'] without duplicates -# # raw.info['bads'] = list(set(raw.info['bads'] + new_bads)) - -# print(f"Threshold: {db_limit} dB -> {abs_threshold:.2e} V^2/Hz") -# print(f"Newly potentially bad channels: {new_bads}") - -# return new_bads - - - -# def find_flat_channels(raw, threshold=1e-15): -# """ -# Identifies channels that are essentially flat lines (zero variance). -# """ -# data = raw.get_data() -# # Calculate standard deviation of each channel -# std_devs = np.std(data, axis=1) - -# # Identify channels with almost zero movement -# flat_idx = np.where(std_devs < threshold)[0] -# flat_names = [raw.ch_names[i] for i in flat_idx] - -# return flat_names - - -# def find_truly_dead_channels(raw, threshold=1e-20): -# """ -# Looks for channels where the signal literally doesn't move -# between samples (successive differences are zero). -# """ -# data = raw.get_data() -# # Calculate the absolute difference between every sample (t and t+1) -# # diffs shape: (n_channels, n_times - 1) -# diffs = np.abs(np.diff(data, axis=1)) - -# # Calculate the mean 'step' size for each channel -# mean_diffs = np.mean(diffs, axis=1) - -# # Identify channels where the signal is 'stuck' -# stuck_idx = np.where(mean_diffs < threshold)[0] -# stuck_names = [raw.ch_names[i] for i in stuck_idx] - -# return stuck_names, mean_diffs - - -# def find_mid_run_flatlines(raw, window_size=10.0, threshold=1e-12): -# """ -# Checks for channels that go flat partway through the recording. -# window_size: Duration in seconds to check for flatness. -# """ -# sfreq = raw.info['sfreq'] -# data = raw.get_data() -# n_samples_win = int(window_size * sfreq) - -# bad_channels = [] - -# for i, ch_name in enumerate(raw.ch_names): -# ch_data = data[i] -# # Create sliding windows (non-overlapping for speed) -# windows = np.array_split(ch_data, len(ch_data) // n_samples_win) - -# # Calculate variance for each window -# win_vars = [np.var(w) for w in windows] - -# # If the LAST window (or any significant portion at the end) is flat -# if win_vars[-1] < threshold: -# bad_channels.append(ch_name) - -# return bad_channels - - -def find_flatline_at_end(raw, threshold_ratio=0.05): - """ - Compares the variance of the first 25% of the data to the last 25%. - If the end variance is less than 5% of the start variance, mark it bad. - """ - data = raw.get_data() - n_samples = data.shape[1] - quarter = n_samples // 4 - - bad_channels = [] - - for i, ch_name in enumerate(raw.ch_names): - start_var = np.var(data[i, :quarter]) - end_var = np.var(data[i, -quarter:]) - - # Avoid division by zero for truly dead channels - if start_var == 0: - bad_channels.append(ch_name) - continue - - ratio = end_var / start_var - - if ratio < threshold_ratio: - bad_channels.append(ch_name) - print(f"Flagged {ch_name}: Variance dropped to {ratio:.2%} of original.") - - return bad_channels - - -# def plot_with_death_lines(raw, picks, window_sec=5.0, var_threshold=0.05): -# sfreq = raw.info['sfreq'] -# data, times = raw.get_data(picks=picks, return_times=True) - -# fig, ax = plt.subplots(figsize=(12, 6)) - -# for i, ch_name in enumerate(picks): -# ch_data = data[i] -# ax.plot(times, ch_data, label=ch_name, alpha=0.8) - -# # --- Find the 'Death Point' --- -# # Calculate rolling variance in 5-second chunks -# win_samples = int(window_sec * sfreq) -# initial_var = np.var(ch_data[:win_samples]) - -# # Check from the end backwards to find where it "died" -# death_time = None -# n_samples = len(ch_data) -# for start_idx in range(n_samples - win_samples - 1, 0, -win_samples): -# current_var = np.var(ch_data[start_idx : start_idx + win_samples]) -# if current_var > (initial_var * var_threshold): -# # The point right after this is where it stays dead -# death_time = times[start_idx + win_samples] -# break - -# if death_time and death_time < (times[-1] - 10): # Only plot if it died early -# ax.axvline(x=death_time, color='red', linestyle='--', alpha=0.5) -# ax.text(death_time, ax.get_ylim()[1], 'Signal Loss', -# color='red', rotation=90, verticalalignment='top') - -# ax.set_title("HbO/HbR with Automated Signal Loss Detection") -# ax.set_ylabel("Concentration (Δ μmol)") -# ax.set_xlabel("Time (s)") -# ax.legend(loc='lower right') -# plt.grid(True, alpha=0.3) -# plt.show(block=True) - - - def detect_sensor_dropout(raw, threshold_ratio=0.05): """ @@ -5251,6 +4815,7 @@ def hr_calc(raw, seconds_to_strip_hr, l_freq, h_freq, search_min, search_max, ma return fig, hr1, hr2, low, high + def trim_participant_data(raw, seconds_to_keep: float): if hasattr(raw, 'annotations') and len(raw.annotations) > 0: # Get time of first event @@ -5278,6 +4843,7 @@ def trim_participant_data(raw, seconds_to_keep: float): return raw, fig_trimmed + def remove_bad_channels(raw, bad_channels, max_bad_channels: int): num_bad = len(bad_channels) @@ -5294,6 +4860,7 @@ def remove_bad_channels(raw, bad_channels, max_bad_channels: int): return raw + def make_and_run_glm(raw_haemo, df_design_matrix, noise_model, bins, n_jobs, verbosity): glm_est = run_glm(raw_haemo, df_design_matrix, noise_model=noise_model, bins=bins, n_jobs=n_jobs, verbose=verbosity) @@ -5340,6 +4907,7 @@ def make_and_run_glm(raw_haemo, df_design_matrix, noise_model, bins, n_jobs, ver return glm_est, fig_glm_topo + def _real_conditions(values, exclude_list=NUISANCE_EXCLUDE): """Filter out drift/constant/short-style nuisance regressor names, keeping only actual task conditions — same filtering logic already @@ -5351,6 +4919,7 @@ def _real_conditions(values, exclude_list=NUISANCE_EXCLUDE): }) + def generate_channel_results(glm_est, file_path): df_cha = glm_est.to_dataframe() df_cha["ID"] = file_path @@ -5363,6 +4932,7 @@ def generate_channel_results(glm_est, file_path): return df_cha + def generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path, json_location): rois_formatted = {} @@ -5468,8 +5038,8 @@ def generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path, json_l return df_roi, fig + def generate_contrast_results(df_design_matrix, glm_est, file_path): - contrast_results_dict = {} contrast_matrix = np.eye(df_design_matrix.shape[1]) @@ -5538,6 +5108,7 @@ def generate_contrast_results(df_design_matrix, glm_est, file_path): return contrast_results_dict + def haemoglobin_concentration(raw_od, file_path, override_ppf=False, ppf_lower_wavelength=6.0, ppf_upper_wavelength=6.0): if override_ppf: raw_haemo = beer_lambert_law(raw_od, ppf=(ppf_lower_wavelength, ppf_upper_wavelength)) @@ -5546,48 +5117,114 @@ def haemoglobin_concentration(raw_od, file_path, override_ppf=False, ppf_lower_w return raw_haemo -def process_participant(file_path, progress_callback=None): - """Goal is to have all of the constants be located here, and nowhere else in the code to remove any ambiguity.""" - # Step 0: Setting up - fig_individual: dict[str, Figure] = {} +def _png_worker(png_queue: Queue, fig_bytes_dict: dict, dpi: int = 100): + """ + Runs on a single dedicated background thread for the lifetime of one + process_participant() call. Consumes (label, fig) pairs and renders them + to PNG bytes. + + Safety: figures are `plt.close()`-d by the producer (main thread) the + instant they're created, BEFORE being queued. That means this thread + never touches pyplot's global figure registry (Gcf) - it only calls + fig.savefig(), which operates on the Figure/Canvas object directly. + Since this is the only thread that ever calls savefig/draw, there's no + concurrent-draw hazard, even while the main thread keeps creating and + closing new figures elsewhere in the pipeline. + """ + while True: + item = png_queue.get() + if item is None: # sentinel - no more figures coming + png_queue.task_done() + break + label, fig = item + try: + buf = BytesIO() + fig.savefig(buf, format="png", dpi=dpi, pil_kwargs={"compress_level": 1}) + fig_bytes_dict[label] = buf.getvalue() + finally: + fig.clf() # drop axes/artists now that we're done, frees memory sooner + png_queue.task_done() + + + +def initial_setup(): + timings = {} + step_start = time.perf_counter() config_dict = { k: globals()[k] for k in __annotations__ - if k in globals() and k != "REQUIRED_KEYS" + if k in globals() } + fig_bytes_dict: dict[str, bytes] = {} + png_queue: Queue = Queue() + png_thread = threading.Thread( + target=_png_worker, args=(png_queue, fig_bytes_dict), daemon=True + ) + png_thread.start() + + return fig_bytes_dict, config_dict, png_queue, timings, step_start + + + +def _enqueue(label, fig, png_queue): + if fig is None: + return + plt.close(fig) # detach from pyplot's global registry - main thread only + png_queue.put((label, fig)) + + + +def lap(start, timings, name): + now = time.perf_counter() + timings[name] = now - start + return now + + + +def process_participant(file_path, file_start, progress_callback=None): + + print(f"File was started with {time.time() - file_start:2f} seconds elapsed.") + # Step 0: Setting up + fig_bytes_dict, config_dict, png_queue, timings, step_start = initial_setup() + step_start = lap(step_start, timings, "Step 0") + # Step 1: Preprocessing raw = load_snirf(file_path=file_path, downsample_frequency=DOWNSAMPLE_FREQUENCY, verbosity=VERBOSITY) fig_raw = raw.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="Loaded Raw", show=False) - fig_individual["Loaded Raw Data"] = fig_raw + _enqueue("Loaded Raw Data", fig_raw, png_queue) if progress_callback: progress_callback(1) logger.info("Step 1 Completed.") - + step_start = lap(step_start, timings, "Step 1") + # Step 2: Trimming if TRIM and not FOLDING_BYP: raw, fig_trimmed = trim_participant_data(raw, seconds_to_keep=SECONDS_TO_KEEP) - fig_individual["Trimmed Raw Data"] = fig_trimmed + _enqueue("Trimmed Raw Data", fig_trimmed, png_queue) if progress_callback: progress_callback(2) logger.info("Step 2 Completed.") + step_start = lap(step_start, timings, "Step 2") # Step 3: Verify Optode Placement if OPTODE_PLACEMENT: fig_optodes = raw.plot_sensors(show_names=SHOW_OPTODE_NAMES, to_sphere=True, show=False, verbose=VERBOSITY) # type: ignore - fig_individual["Plot Sensors"] = fig_optodes + _enqueue("Plot Sensors", fig_optodes, png_queue) if progress_callback: progress_callback(3) logger.info("Step 3 Completed.") + step_start = lap(step_start, timings, "Step 3") # Step 4: Short/Long Channels if SHORT_CHANNELS and not FOLDING_BYP: #NOTE: Have to split again later but since needed for heart rate, this will stay at step 4. Will split later again. _short_chans = get_short_channels(raw, max_dist=SHORT_CHANNELS_THRESHOLD) # JUST FOR PLOTTING THEM SEPERATELY fig_short_chans = _short_chans.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="Short Channels Only", show=False, verbose=VERBOSITY) - fig_individual["Short Channels Raw Data"] = fig_short_chans + _enqueue("Short Channels Raw Data", fig_short_chans, png_queue) if LONG_CHANNELS: raw = get_long_channels(raw, min_dist=0, max_dist=LONG_CHANNELS_THRESHOLD) if progress_callback: progress_callback(4) logger.info("Step 4 Completed.") + step_start = lap(step_start, timings, "Step 4") # Step 5: Heart Rate if HEART_RATE and not FOLDING_BYP: @@ -5606,12 +5243,13 @@ def process_participant(file_path, progress_callback=None): short_channels_threshold=SHORT_CHANNELS_THRESHOLD, verbosity=VERBOSITY ) - fig_individual["Power Spectral Density"] = fig - fig_individual['Heart Rate - PSD'] = hr1 - fig_individual['Heart Rate - Time'] = hr2 + _enqueue("Power Spectral Density", fig, png_queue) + _enqueue('Heart Rate - PSD', hr1, png_queue) + _enqueue('Heart Rate - Time', hr2, png_queue) if progress_callback: progress_callback(5) logger.info("Step 5 Completed.") - + step_start = lap(step_start, timings, "Step 5") + # Step 6: Scalp Coupling Index bad_sci = [] if SCI and not FOLDING_BYP: @@ -5619,112 +5257,125 @@ def process_participant(file_path, progress_callback=None): bad_sci, fig_sci_1, fig_sci_2 = calculate_scalp_coupling(raw, l_freq=low, h_freq=high, time_window=SCI_TIME_WINDOW, threshold=SCI_THRESHOLD) else: bad_sci, fig_sci_1, fig_sci_2 = calculate_scalp_coupling(raw, l_freq=SCI_LOW_FREQ, h_freq=SCI_HIGH_FREQ, time_window=SCI_TIME_WINDOW, threshold=SCI_THRESHOLD) - fig_individual["Scalp Coupling Index Heatmap"] = fig_sci_1 - fig_individual["Scalp Coupling Index Binary Heatmap"] = fig_sci_2 + _enqueue("Scalp Coupling Index Heatmap", fig_sci_1, png_queue) + _enqueue("Scalp Coupling Index Binary Heatmap", fig_sci_2, png_queue) if progress_callback: progress_callback(6) logger.info("Step 6 Completed.") + step_start = lap(step_start, timings, "Step 6") # Step 7: Signal to Noise Ratio bad_snr = [] if SNR and not FOLDING_BYP: bad_snr, fig_snr = calculate_signal_noise_ratio(raw) - fig_individual["Signal To Noise Ratio"] = fig_snr + _enqueue("Signal To Noise Ratio", fig_snr, png_queue) if progress_callback: progress_callback(7) logger.info("Step 7 Completed.") + step_start = lap(step_start, timings, "Step 7") # Step 8: Peak Spectral Power bad_psp = [] if PSP and not FOLDING_BYP: bad_psp, fig_psp1, fig_psp2 = calculate_peak_power(raw, time_window=PSP_TIME_WINDOW, threshold=PSP_THRESHOLD, l_freq=PSP_LOW_FREQ, h_freq=PSP_HIGH_FREQ) - fig_individual["Peak Spectral Power Heatmap"] = fig_psp1 - fig_individual["Peak Spectral Power Binary Heatmap"] = fig_psp2 + _enqueue("Peak Spectral Power Heatmap", fig_psp1, png_queue) + _enqueue("Peak Spectral Power Binary Heatmap", fig_psp2, png_queue) if progress_callback: progress_callback(8) logger.info("Step 8 Completed.") + step_start = lap(step_start, timings, "Step 8") # Step 9: Coefficient of Variation bad_coeff_var = [] if COEFF_VAR and not FOLDING_BYP: bad_coeff_var, fig_coeff_var = find_bad_channels_coeff_var(raw, coeff_var_threshold=COEFF_VAR_THRESHOLD) - fig_individual['Coefficient of Variation'] = fig_coeff_var + _enqueue('Coefficient of Variation', fig_coeff_var, png_queue) if progress_callback: progress_callback(9) logger.info("Step 9 Completed.") + step_start = lap(step_start, timings, "Step 9") # Step 10: Median Absolute Deviation bad_amplitude_range = [] if MAD and not FOLDING_BYP: bad_amplitude_range, fig_range = find_bad_channels_by_amplitude_range(raw, threshold=MAD_THRESHOLD) - fig_individual['Median Absolute Deviation'] = fig_range + _enqueue('Median Absolute Deviation', fig_range, png_queue) if progress_callback: progress_callback(10) logger.info("Step 10 Completed.") + step_start = lap(step_start, timings, "Step 10") # Step 11: Power Spectral Density Noise bad_noise = [] if PSD_NOISE and not FOLDING_BYP: bad_noise, fig_noise = detect_spectral_noise_spike(raw, db_limit=DB_LIMIT, freq_div=TARGET_FREQ_DIV, min_freq=PSD_MIN_FREQ, target_bandwith=PSD_TARGET_BANDWIDTH) - fig_individual['Power Spectral Density Noise'] = fig_noise + _enqueue('Power Spectral Density Noise', fig_noise, png_queue) if progress_callback: progress_callback(11) logger.info("Step 11 Completed.") + step_start = lap(step_start, timings, "Step 11") # Step 12: Channel Dropout bad_disp = [] if SENSOR_DROPOUT and not FOLDING_BYP: bad_disp, fig_disp = detect_sensor_dropout(raw, threshold_ratio=SENSOR_DROPOUT_VARIANCE_THRESHOLD) - fig_individual['Sensor Dropout'] = fig_disp + _enqueue('Sensor Dropout', fig_disp, png_queue) if progress_callback: progress_callback(12) logger.info("Step 12 Completed.") + step_start = lap(step_start, timings, "Step 12") # Step 13: Bad Channels Handling if BAD_CHANNELS_HANDLING != "None" and not FOLDING_BYP: raw, fig_dropped, fig_raw_before, bad_channels = mark_bads(raw, bad_sci, bad_snr, bad_psp, bad_coeff_var, bad_amplitude_range, bad_noise, bad_disp) if fig_dropped and fig_raw_before is not None: - fig_individual["Bad Channels by Method"] = fig_dropped - fig_individual["Bad Channels Data"] = fig_raw_before + _enqueue("Bad Channels by Method", fig_dropped, png_queue) + _enqueue("Bad Channels Data", fig_raw_before, png_queue) if bad_channels: if BAD_CHANNELS_HANDLING == "Interpolate": raw, fig_raw_after, fig_compare = interpolate_fNIRS_bads_weighted_average(raw, max_dist=MAX_DIST, min_neighbors=MIN_NEIGHBORS, short_channels_threshold=SHORT_CHANNELS_THRESHOLD) - fig_individual["Data after Interpolating Bad Channels"] = fig_raw_after - fig_individual["Bad Channels Interpolation Results"] = fig_compare + _enqueue("Data after Interpolating Bad Channels", fig_raw_after, png_queue) + _enqueue("Bad Channels Interpolation Results", fig_compare, png_queue) elif BAD_CHANNELS_HANDLING == "Remove": raw = remove_bad_channels(raw, bad_channels, max_bad_channels=MAX_BAD_CHANNELS) if progress_callback: progress_callback(13) logger.info("Step 13 Completed.") + step_start = lap(step_start, timings, "Step 13") # Step 14: Optical Density raw_od = optical_density(raw) fig_raw_od = raw_od.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="Optical Density", show=False) - fig_individual["Optical Density"] = fig_raw_od + _enqueue("Optical Density", fig_raw_od, png_queue) if progress_callback: progress_callback(14) logger.info("Step 14 Completed.") + step_start = lap(step_start, timings, "Step 14") # Step 15: Temporal Derivative Distribution Repair Filtering if TDDR and not FOLDING_BYP: raw_od = temporal_derivative_distribution_repair(raw_od) fig_raw_od_tddr = raw_od.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="After TDDR (Motion Correction)", show=False) - fig_individual["Temporal Derivative Distribution Repair"] = fig_raw_od_tddr + _enqueue("Temporal Derivative Distribution Repair", fig_raw_od_tddr, png_queue) if progress_callback: progress_callback(15) logger.info("Step 15 Completed.") + step_start = lap(step_start, timings, "Step 15") # Step 16: Wavelet Filtering if WAVELET and not FOLDING_BYP: raw_od, fig = calculate_and_apply_wavelet(data=raw_od, wavelet_type=WAVELET_TYPE, wavelet_level=WAVELET_LEVEL, iqr=IQR, verbosity=VERBOSITY) - fig_individual["Wavelet"] = fig + _enqueue("Wavelet", fig, png_queue) if progress_callback: progress_callback(16) logger.info("Step 16 Completed.") + step_start = lap(step_start, timings, "Step 16") # Step 17: Haemoglobin Concentration raw_haemo = haemoglobin_concentration(raw_od, file_path, OVERRIDE_PPF, PPF_LOWER_WAVELENGTH, PPF_UPPER_WAVELENGTH) fig_raw_haemo_bll = raw_haemo.plot(duration=raw_haemo.times[-1], n_channels=raw_haemo.info['nchan'], title="HbO and HbR Signals", show=False) - fig_individual["Modified Beer Lambert Law"] = fig_raw_haemo_bll + _enqueue("Modified Beer Lambert Law", fig_raw_haemo_bll, png_queue) if progress_callback: progress_callback(17) logger.info("Step 17 Completed.") + step_start = lap(step_start, timings, "Step 17") # Step 18: Enhance Negative Correlation if ENHANCE_NEGATIVE_CORRELATION and not FOLDING_BYP: raw_haemo = enhance_negative_correlation(raw_haemo) fig_raw_haemo_enc = raw_haemo.plot(duration=raw_haemo.times[-1], n_channels=raw_haemo.info['nchan'], title="Enhance Negative Correlation", show=False) - fig_individual["Enhance Negative Correlation"] = fig_raw_haemo_enc + _enqueue("Enhance Negative Correlation", fig_raw_haemo_enc, png_queue) if progress_callback: progress_callback(18) logger.info("Step 18 Completed.") + step_start = lap(step_start, timings, "Step 18") # Step 19: Filter if FILTER and not FOLDING_BYP: @@ -5735,32 +5386,34 @@ def process_participant(file_path, progress_callback=None): h_freq=H_FREQ, l_trans_bandwidth=L_TRANS_BANDWIDTH, h_trans_bandwidth=H_TRANS_BANDWIDTH, - iir_type=IIR_TYPE, - iir_order=IIR_ORDER, + # iir_type=IIR_TYPE, + # iir_order=IIR_ORDER, filter_length=FILTER_LENGTH, filter_phase=FILTER_PHASE, fir_window=FIR_WINDOW, fir_design=FIR_DESIGN, - iir_output=IIR_OUTPUT, - passband_ripple=PASSBAND_RIPPLE, - stopband_attenuation=STOPBAND_ATTENUATION, + # iir_output=IIR_OUTPUT, + # passband_ripple=PASSBAND_RIPPLE, + # stopband_attenuation=STOPBAND_ATTENUATION, filter_pad=FILTER_PAD, - skip_by_annotation=SKIP_BY_ANNOTATION, + # skip_by_annotation=SKIP_BY_ANNOTATION, filter_n_jobs=FILTER_N_JOBS, verbosity=VERBOSITY ) - fig_individual["Filter_1"] = fig_filter - fig_individual["Filter_2"] = fig_raw_haemo_filter + _enqueue("Filter_1", fig_filter, png_queue) + _enqueue("Filter_2", fig_raw_haemo_filter, png_queue) if progress_callback: progress_callback(19) logger.info("Step 19 Completed.") + step_start = lap(step_start, timings, "Step 19") # Step 20: Extracting Events if EVENTS and not FOLDING_BYP: events, event_dict = events_from_annotations(raw_haemo, event_id=EVENT_ID, regexp=EVENT_REGEX, verbose=VERBOSITY) #TODO: Implement the chunk duration fig_events = plot_events(events, event_id=event_dict, sfreq=raw_haemo.info["sfreq"], show=False) - fig_individual["Events"] = fig_events + _enqueue("Events", fig_events, png_queue) if progress_callback: progress_callback(20) logger.info("Step 20 Completed.") + step_start = lap(step_start, timings, "Step 20") # Step 21: Epoch Calculations if EPOCHS and EVENTS and not FOLDING_BYP: @@ -5777,9 +5430,10 @@ def process_participant(file_path, progress_callback=None): reject_hbo_threshold=dict(hbo=REJECT_HBO_THRESHOLD) ) for name, fig in fig_epochs: - fig_individual[f"epochs_{name}"] = fig + _enqueue(f"epochs_{name}", fig, png_queue) if progress_callback: progress_callback(21) logger.info("Step 21 Completed.") + step_start = lap(step_start, timings, "Step 21") # Step 22: Design Matrix raw_haemo, df_design_matrix, fig_design_matrix = make_design_matrix( @@ -5801,51 +5455,63 @@ def process_participant(file_path, progress_callback=None): long_channels_threshold=LONG_CHANNELS_THRESHOLD, folding_bypass=FOLDING_BYP ) - fig_individual["Design Matrix"] = fig_design_matrix + _enqueue("Design Matrix", fig_design_matrix, png_queue) if progress_callback: progress_callback(22) logger.info("Step 22 Completed.") + step_start = lap(step_start, timings, "Step 22") # Step 23: General Linear Model glm_est, fig_glm_topo = make_and_run_glm(raw_haemo, df_design_matrix, noise_model=NOISE_MODEL, bins=BINS, n_jobs=N_JOBS, verbosity=VERBOSITY) - fig_individual["GLM Topography"] = fig_glm_topo + _enqueue("GLM Topography", fig_glm_topo, png_queue) if progress_callback: progress_callback(23) logger.info("23") + step_start = lap(step_start, timings, "Step 23") # Step 24: Generate GLM Results if "derivative" not in HRF_MODEL.lower(): fig_glm_result = plot_glm_results(file_path, raw_haemo, glm_est, df_design_matrix) for name, fig in fig_glm_result: - fig_individual[f"GLM {name}"] = fig + _enqueue(f"GLM {name}", fig, png_queue) if progress_callback: progress_callback(24) logger.info("24") + step_start = lap(step_start, timings, "Step 24") # Step 25: Generate Channel Results df_cha = generate_channel_results(glm_est, file_path) if progress_callback: progress_callback(25) logger.info("25") + step_start = lap(step_start, timings, "Step 25") # Step 26: Generate Region of Interest Results df_roi, fig_roi = generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path, json_location=JSON_LOCATION) - fig_individual["Region of Interest"] = fig_roi + _enqueue("Region of Interest", fig_roi, png_queue) if progress_callback: progress_callback(26) logger.info("26") + step_start = lap(step_start, timings, "Step 26") # Step 27: Generate Contrast Results contrast_results_dict = generate_contrast_results(df_design_matrix, glm_est, file_path) if progress_callback: progress_callback(27) logger.info("27") - + step_start = lap(step_start, timings, "Step 27") + # Step 28: Finishing Up - fig_bytes_dict = convert_fig_dict_to_png_bytes(fig_individual) + png_queue.put(None) # sentinel + png_queue.join() # blocks only on remaining unfinished work if FOLDING_BYP: epochs = None sanitize_paths_for_pickle(raw_haemo, epochs) if progress_callback: progress_callback(28) logger.info("28") + step_start = lap(step_start, timings, "Step 28") # Step 28.5: Return the results - return raw_haemo, epochs, df_cha, df_roi, df_design_matrix, config_dict, fig_bytes_dict, contrast_results_dict, True + logger.info("Step timings:") + for name, elapsed in timings.items(): + logger.info(f" {name:<25} {elapsed:7.3f}s") + logger.info(f"Total processing time: {sum(timings.values()):.3f}s") + return raw_haemo, epochs, df_cha, df_roi, df_design_matrix, config_dict, fig_bytes_dict, contrast_results_dict, True @@ -5859,6 +5525,7 @@ def sanitize_paths_for_pickle(raw_haemo, epochs): epochs._raw._filenames = [str(p) for p in epochs._raw._filenames] + def functional_connectivity_spectral_epochs( epochs: DataFrame | None, n_lines: int, @@ -5900,8 +5567,6 @@ def functional_connectivity_spectral_epochs( - - def functional_connectivity_spectral_time( epochs: DataFrame | None, n_lines: int, @@ -5949,7 +5614,6 @@ def functional_connectivity_spectral_time( - def functional_connectivity_envelope( epochs: DataFrame | None, n_lines: int, @@ -5985,6 +5649,7 @@ def functional_connectivity_envelope( ) + def functional_connectivity_betas( raw_hbo: BaseRaw, n_lines: int, @@ -6021,9 +5686,6 @@ def functional_connectivity_betas( reg_names = list(design_matrix.columns) - - - n_channels = betas.shape[0] # ------------------------------------------------------------------ @@ -6054,7 +5716,6 @@ def functional_connectivity_betas( beta_series[:, t] = np.mean(betas[:, idx], axis=1).flatten() - # n_channels, n_trials = betas.shape[0], len(onsets) # beta_series = np.zeros((n_channels, n_trials)) @@ -6108,7 +5769,6 @@ def functional_connectivity_betas( - def get_single_subject_beta_corr(raw_hbo, event_name=None, config=None): """Processes one participant and returns their correlation matrix.""" raw_hbo = raw_hbo.copy().pick(picks="hbo") @@ -6191,6 +5851,7 @@ def get_single_subject_beta_corr(raw_hbo, event_name=None, config=None): return corr_matrix, raw_hbo.ch_names + def run_group_functional_connectivity( haemo_dict: dict[str | Path, BaseRaw], config_dict: dict[str, Any], @@ -6302,6 +5963,7 @@ def run_group_functional_connectivity( ) + def sparks_csv_export( haemo_obj: BaseRaw, save_path: str, @@ -6324,3 +5986,155 @@ def sparks_csv_export( df.insert(0, "annotation", ann_col) df.insert(0, "time", times) df.to_csv(save_path, index=False) + + + +def peak_power_fast( + raw, + time_window=10, + threshold=0.1, + l_freq=0.7, + h_freq=1.5, + l_trans_bandwidth=0.3, + h_trans_bandwidth=0.3, + verbose=False, +): + """ + Numerically-equivalent drop-in for mne_nirs' peak_power(). + + Same algorithm, same math, same return shape/order. The only change: + the periodogram() call is batched across ALL (channel-pair, window) + combinations in a single vectorized call instead of one call per + combination. np.correlate() is still called per (pair, window) - it's + cheap relative to periodogram's per-call overhead - but if profiling + shows it's now the bottleneck, that's the next thing to vectorize + (via FFT-based batch cross-correlation), flagged separately since + it's a bigger, riskier change to get bit-exact. + + Validate before trusting in production: + _, scores_orig, times_orig = peak_power(raw, ...) + _, scores_fast, times_fast = peak_power_fast(raw, ...) + assert np.allclose(scores_orig, scores_fast) + """ + raw = raw.copy().load_data() + _validate_type(raw, BaseRaw, "raw") + + picks = _validate_nirs_info(raw.info) + sfreq = raw.info["sfreq"] + + filtered_data = filter_data( + raw._data, + sfreq, + l_freq, + h_freq, + picks=picks, + verbose=verbose, + l_trans_bandwidth=l_trans_bandwidth, + h_trans_bandwidth=h_trans_bandwidth, + ) + + window_samples = int(np.ceil(time_window * sfreq)) + n_windows = int(np.floor(len(raw) / window_samples)) + n_pairs = len(picks) // 2 + + scores = np.zeros((len(picks), n_windows)) + times = [] + + if n_windows == 0: + scores = scores[np.argsort(picks)] + return raw, scores, times + + # All windows except the last are guaranteed to be exactly window_samples + # long (the last window can be shorter due to the min() clamp below) - + # batch-process those; handle the last window with the original scalar path. + full_windows = n_windows - 1 + corr_len = 2 * window_samples - 1 + + if full_windows > 0: + c1_stack = np.empty((n_pairs, full_windows, window_samples)) + c2_stack = np.empty((n_pairs, full_windows, window_samples)) + + for window in range(full_windows): + start = window * window_samples + end = start + window_samples + for pi, ii in enumerate(range(0, len(picks), 2)): + c1_stack[pi, window] = filtered_data[picks[ii]][start:end] + c2_stack[pi, window] = filtered_data[picks[ii + 1]][start:end] + + std1 = c1_stack.std(axis=-1, keepdims=True) + std1[std1 == 0] = 1 + std2 = c2_stack.std(axis=-1, keepdims=True) + std2[std2 == 0] = 1 + c1_stack = c1_stack / std1 + c2_stack = c2_stack / std2 + + corr_stack = np.empty((n_pairs, full_windows, corr_len)) + for pi in range(n_pairs): + for window in range(full_windows): + corr_stack[pi, window] = ( + np.correlate(c1_stack[pi, window], c2_stack[pi, window], "full") + / window_samples + ) + + # single vectorized call replaces n_pairs * full_windows separate calls + _, pxx = periodogram(corr_stack, fs=sfreq, window="hamming", axis=-1) + window_scores = pxx.max(axis=-1) # shape (n_pairs, full_windows) + + scores[0::2, :full_windows] = window_scores + scores[1::2, :full_windows] = window_scores + + for window in range(full_windows): + start = window * window_samples + end = start + window_samples + times.append((raw.times[start], raw.times[min(end, len(raw) - 1)])) + if threshold is not None: + for pi in np.where(window_scores[:, window] < threshold)[0]: + ii = pi * 2 + raw.annotations.append( + raw.times[start], + time_window, + "BAD_PeakPower", + ch_names=[raw.ch_names[ii : ii + 2]], + ) + + # last (possibly truncated) window - original scalar path, unchanged + window = n_windows - 1 + start_sample = window * window_samples + end_sample = int(np.min([start_sample + window_samples, len(raw) - 1])) + t_start, t_stop = raw.times[start_sample], raw.times[end_sample] + times.append((t_start, t_stop)) + + for ii in range(0, len(picks), 2): + c1 = filtered_data[picks[ii]][start_sample:end_sample] + c2 = filtered_data[picks[ii + 1]][start_sample:end_sample] + c1 = c1 / (np.std(c1) or 1) + c2 = c2 / (np.std(c2) or 1) + c = np.correlate(c1, c2, "full") / window_samples + f, pxx = periodogram(c, fs=sfreq, window="hamming") + scores[ii, window] = max(pxx) + scores[ii + 1, window] = max(pxx) + if (threshold is not None) and (max(pxx) < threshold): + raw.annotations.append( + t_start, time_window, "BAD_PeakPower", + ch_names=[raw.ch_names[ii : ii + 2]], + ) + + scores = scores[np.argsort(picks)] + return raw, scores, times + + +if __name__ == "__main__": + print("This file has no functionality when not used in tandem with the FLARES application.") + + # audit_log = logging.getLogger("memory_audit") + # audit_log.setLevel(logging.INFO) + # audit_log.propagate = False # This prevents it from talking to other loggers + + # # 2. Add a file handler specifically for this audit logger + # if not audit_log.handlers: + # fh = logging.FileHandler('flares_memory_audit.log') + # fh.setFormatter(logging.Formatter('%(asctime)s | PID: %(process)d | %(message)s')) + # audit_log.addHandler(fh) + + # def get_mem_mb(): + # return psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 \ No newline at end of file