From 61a9ea2f34fbc6562b7635aee41f99d8afda2b90 Mon Sep 17 00:00:00 2001 From: tyler Date: Sat, 1 Aug 2026 19:20:49 -0700 Subject: [PATCH] fix rare bug --- changelog.md | 13 +++- flares.py | 181 +++++++++++++++++++++++++-------------------------- 2 files changed, 97 insertions(+), 97 deletions(-) diff --git a/changelog.md b/changelog.md index 0d08356..5bebbdb 100644 --- a/changelog.md +++ b/changelog.md @@ -1,3 +1,10 @@ +# Verison 1.5.3 + +- Optimized calculations being performed when calculating the heart rate to speed up step 5 by up to ~35% on a per-file basis +- Optimized calculations being performed when running the General Linear Model to speed up step 5 by ~35% on a per-file basis +- Fixed an issue where participants could be skipped when processing multiple particants at one which could prevent overall processing from completing + + # Version 1.5.2 - Opening both a file or a folder now contains support for reading some metadata from the BIDS structure. Both options will still function if this metadata is not present @@ -12,9 +19,9 @@ - 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% +- Optimized some of the calculations in Scalp Coupling Index to speed up Step 6 by up to ~25% on a per-file basis +- Removed duplicate/redundant calculations in Peak Spectral Power to speed up Step 8 by up to ~50% on a per-file basis +- Changed how the figures are generated when processing to speed up Step 28 by up to ~85% on a per-file basis - 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 diff --git a/flares.py b/flares.py index 961a614..1c99138 100644 --- a/flares.py +++ b/flares.py @@ -54,7 +54,7 @@ from statsmodels.tools.sm_exceptions import ConvergenceWarning from scipy.spatial.distance import cdist from scipy.signal import welch, butter, filtfilt, periodogram # type: ignore -from scipy.stats import pearsonr, zscore, ttest_1samp, ttest_ind, sem +from scipy.stats import pearsonr, zscore, ttest_1samp, ttest_ind, sem, t as t_dist import pywt # type: ignore import neurokit2 as nk # type: ignore @@ -104,9 +104,6 @@ from src.shared.shareddata import PLATFORM_NAME, resource_path -# Needs to be set for mne -os.environ["SUBJECTS_DIR"] = str(data_path()) + "/subjects" # type: ignore - PRIMARY_COLORS = { "SCI only": "skyblue", # Scalp Coupling Index (Standard MNE) "SNR only": "lightgreen", # Signal-to-Noise Ratio (Your original) @@ -1246,7 +1243,7 @@ def mark_bads(raw, bad_sci, bad_snr, bad_psp, bad_coeff_var, bad_range, bad_noi fig_dropped.tight_layout() - raw_before = deepcopy(raw) + raw_before = raw.copy() bads_channels = [ch for ch in raw.ch_names if ch in raw.info['bads']] print(bads_channels) if bads_channels: @@ -1359,8 +1356,7 @@ def safe_create_epochs(raw, events, event_dict, tmin, tmax, baseline, max_shift, -def epochs_calculations(raw_haemo, events, event_dict, epoch_handling, max_shift, t_min, t_max, baseline, reject_epochs, reject_hbo_threshold): - fig_epochs = [] # List to store figures +def epochs_calculations(raw_haemo, events, event_dict, epoch_handling, max_shift, t_min, t_max, baseline, reject_epochs, reject_hbo_threshold, png_queue): if epoch_handling == 'shift': epochs = safe_create_epochs(raw=raw_haemo, events=events, event_dict=event_dict, tmin=t_min, tmax=t_max, baseline=baseline, max_shift=max_shift, reject_epochs=reject_epochs, reject_hbo_threshold=reject_hbo_threshold) @@ -1377,12 +1373,18 @@ def epochs_calculations(raw_haemo, events, event_dict, epoch_handling, max_shift # Plot drop log # TODO: Why show this if we never use epochs2? fig_epochs_dropped = epochs2.plot_drop_log(show=False) - fig_epochs.append(("fig_epochs_dropped", fig_epochs_dropped)) + _enqueue("epochs_fig_epochs_dropped", fig_epochs_dropped, png_queue) + + conditions = list(epochs.event_id.keys()) + evoked_cache = {} # Plot for each condition for idx, condition in enumerate(epochs.event_id.keys()): logger.info(condition) logger.info(idx) + + epo_cond = epochs[condition] + # Plot images for each condition fig_epochs_data = epochs[condition].plot_image( combine="mean", @@ -1399,22 +1401,23 @@ def epochs_calculations(raw_haemo, events, event_dict, epoch_handling, max_shift ax = fig.axes[0] original_title = ax.get_title() ax.set_title(f"{condition}: {original_title}") - fig_epochs.append((f"fig_{condition}_data_{idx}_{j}", fig)) # Store with a unique name + _enqueue(f"epochs_fig_{condition}_data_{idx}_{j}", fig, png_queue) # Evoked average figure for each condition - evoked_avg = epochs[condition].average() + evoked_avg = epo_cond.average() + evoked_cache[condition] = evoked_avg + clims = dict(hbo=[-1, 1], hbr=[1, -1]) condition_fig = evoked_avg.plot_image(clim=clims, show=False) for ax in condition_fig.axes: original_title = ax.get_title() ax.set_title(f"{original_title} - {condition}") - fig_epochs.append((f"evoked_avg_{condition}", condition_fig)) # Store with a unique name + _enqueue(f"epochs_evoked_avg_{condition}", condition_fig, png_queue) # Prepare evokeds and colors for topographic plot evokeds3 = [] colors = [] - conditions = list(epochs.event_id.keys()) cmap = plt.get_cmap("tab10", len(conditions)) for idx, cond in enumerate(conditions): @@ -1433,20 +1436,24 @@ def epochs_calculations(raw_haemo, events, event_dict, epoch_handling, max_shift lines.append(line) fig.legend(lines, conditions, loc="lower right") - fig_epochs.append(("evoked_topo", help)) # Store with a unique name + _enqueue("epochs_evoked_topo", help, png_queue) unique_annotations = set(raw_haemo.annotations.description) for cond in unique_annotations: # Evoked response for specific condition ("Activity") - evoked_stim1 = epochs[cond].average() + evoked_stim1 = evoked_cache.get(cond) + if evoked_stim1 is None: + # Not one of the epoch conditions already averaged above - fall back + # to computing it directly (matches original behavior in that case). + evoked_stim1 = epochs[cond].average() fig_evoked_hbo = evoked_stim1.copy().pick(picks='hbo').plot(time_unit='s', show=False) fig_evoked_hbr = evoked_stim1.copy().pick(picks='hbr').plot(time_unit='s', show=False) - fig_epochs.append((f"fig_evoked_hbo_{cond}", fig_evoked_hbo)) # Store with a unique name - fig_epochs.append((f"fig_evoked_hbr_{cond}", fig_evoked_hbr)) # Store with a unique name + _enqueue(f"epochs_fig_evoked_hbo_{cond}", fig_evoked_hbo, png_queue) + _enqueue(f"epochs_fig_evoked_hbr_{cond}", fig_evoked_hbr, png_queue) print("Evoked HbO peak amplitude:", evoked_stim1.copy().pick(picks='hbo').data.max()) @@ -1459,11 +1466,11 @@ def epochs_calculations(raw_haemo, events, event_dict, epoch_handling, max_shift for condition in epochs.event_id: if condition not in all_evokeds: all_evokeds[condition] = [] - all_evokeds[condition].append(epochs[condition].average()) + all_evokeds[condition].append(evoked_cache[condition]) + group_averages = {cond: evoked_cache[cond] for cond in conditions if cond in evoked_cache} group_aucs = {} - # TODO: group averages with a single person? - group_averages = {cond: grand_average(evokeds) for cond, evokeds in all_evokeds.items()} + for condition, evoked in group_averages.items(): group_aucs[condition] = {} for pick in ["hbo", "hbr"]: @@ -1513,9 +1520,9 @@ def epochs_calculations(raw_haemo, events, event_dict, epoch_handling, max_shift ax.legend(legend_labels) - fig_epochs.append((f"fig_{condition}_compare_evokeds", fig)) # Store with a unique name + _enqueue(f"epochs_fig_{condition}_compare_evokeds", fig, png_queue) - return epochs, fig_epochs + return epochs @@ -1988,7 +1995,11 @@ def plot_3d_evoked_array( ea = ea.pick(picks=picks) # type: ignore if subjects_dir is None: - subjects_dir = os.environ["SUBJECTS_DIR"] + subjects_dir = os.environ.get("SUBJECTS_DIR") + if subjects_dir is None: + subjects_dir = str(data_path()) + "/subjects" # type: ignore + os.environ["SUBJECTS_DIR"] = subjects_dir + if src is None: fname_src_fs = os.path.join( subjects_dir, "fsaverage", "bem", "fsaverage-ico-5-src.fif" @@ -3985,14 +3996,12 @@ def aggregate_channel_contrasts_to_roi( group_cols = ['ROI', 'contrast_name', 'Chroma', 'ID'] - def _weighted_mean(g): - return np.average(g['effect'], weights=g['weight']) - - roi_theta = ( - df.groupby(group_cols, group_keys=False) - .apply(lambda g: pd.Series({'theta': _weighted_mean(g)})) - .reset_index() + df['_effect_weight'] = df['effect'] * df['weight'] + roi_theta = df.groupby(group_cols, as_index=False).agg( + _sum_ew=('_effect_weight', 'sum'), + _sum_w=('weight', 'sum'), ) + roi_theta['theta'] = roi_theta['_sum_ew'] / roi_theta['_sum_w'] roi_theta = roi_theta.rename(columns={'contrast_name': 'Condition'}) return roi_theta[['ROI', 'Condition', 'Chroma', 'theta', 'ID']] @@ -4776,7 +4785,7 @@ def hr_calc(raw, seconds_to_strip_hr, l_freq, h_freq, search_min, search_max, ma hr1, hr2 = plot_heart_rate(freq_bpm_scipy, psd_scipy, freq_range_scipy, mean_hr_scipy, hr_smooth_nk, mean_hr_nk, times_trimmed, overruled, hr_window=hr_window) - fig = raw.plot_psd(show=False) + fig = raw.compute_psd().plot(show=False) raw_filtered = raw.copy().filter(0.5, 3, fir_design='firwin') sfreq = raw.info['sfreq'] data = raw_filtered.get_data() @@ -4787,23 +4796,28 @@ def hr_calc(raw, seconds_to_strip_hr, l_freq, h_freq, search_min, search_max, ma nperseg = int(sfreq / desired_bin_hz) hr_range = (search_min, search_max) - # --- Function to find strongest local peak --- - def find_hr_from_psd(ch_data): - f, Pxx = welch(ch_data, sfreq, nperseg=nperseg) - mask = (f >= hr_range[0]/60) & (f <= hr_range[1]/60) - f_masked = f[mask] - Pxx_masked = Pxx[mask] - if len(Pxx_masked) < 3: - return np.nan - peaks = [i for i in range(1, len(Pxx_masked)-1) - if Pxx_masked[i] > Pxx_masked[i-1] and Pxx_masked[i] > Pxx_masked[i+1]] - if not peaks: - return np.nan - best_idx = peaks[np.argmax([Pxx_masked[i] for i in peaks])] - return f_masked[best_idx] * 60 # bpm - - # --- Compute HR across all channels --- - hr_all_channels = np.array([find_hr_from_psd(data[i, :]) for i in range(len(channel_names))]) + f, Pxx = welch(data, fs=sfreq, nperseg=nperseg, axis=1) # data: (n_channels, n_samples) + mask = (f >= hr_range[0] / 60) & (f <= hr_range[1] / 60) + f_masked = f[mask] + Pxx_masked = Pxx[:, mask] # (n_channels, n_freq_in_range) + + hr_all_channels = np.full(Pxx_masked.shape[0], np.nan) + if Pxx_masked.shape[1] >= 3: + # same "strictly greater than both neighbors" local-max definition as + # the original per-channel loop, vectorized across all channels at once + interior = Pxx_masked[:, 1:-1] + left = Pxx_masked[:, :-2] + right = Pxx_masked[:, 2:] + is_local_peak = (interior > left) & (interior > right) + + for ch in range(Pxx_masked.shape[0]): + peak_offsets = np.where(is_local_peak[ch])[0] + if len(peak_offsets) == 0: + continue + candidate_idx = peak_offsets + 1 # shift back into Pxx_masked indexing + best_idx = candidate_idx[np.argmax(Pxx_masked[ch, candidate_idx])] + hr_all_channels[ch] = f_masked[best_idx] * 60 # bpm + hr_all_channels = hr_all_channels[~np.isnan(hr_all_channels)] hr_mode = np.round(np.median(hr_all_channels)) # Use median if some NaNs @@ -4873,20 +4887,19 @@ def make_and_run_glm(raw_haemo, df_design_matrix, noise_model, bins, n_jobs, ver # Extract base task conditions (e.g., "Tapping_Left", "Tapping_Right") base_conditions = list(set(col.split('_delay_')[0] for col in fir_cols)) + theta_list = glm_est.theta() + columns = list(df_design_matrix.columns) + peak_conditions = [] for cond in base_conditions: - # Find all delays corresponding to this specific condition cond_delays = [col for col in fir_cols if col.startswith(f"{cond}_delay_")] - - # Find the delay with the highest average absolute effect (theta) across channels + delay_impacts = {} for col in cond_delays: - col_idx = list(df_design_matrix.columns).index(col) - # glm_est.theta() returns list of theta arrays (one array per channel) - avg_absolute_theta = np.mean(np.abs([ch_theta[col_idx] for ch_theta in glm_est.theta()])) + col_idx = columns.index(col) + avg_absolute_theta = np.mean(np.abs([ch_theta[col_idx] for ch_theta in theta_list])) delay_impacts[col] = avg_absolute_theta - - # Pick the delay column with the absolute largest channel-wide effect + peak_delay_col = max(delay_impacts, key=delay_impacts.get) logger.info(f"Condition '{cond}' peak response identified at delay column: {peak_delay_col}") peak_conditions.append(peak_delay_col) @@ -5417,7 +5430,7 @@ def process_participant(file_path, file_start, progress_callback=None): # Step 21: Epoch Calculations if EPOCHS and EVENTS and not FOLDING_BYP: - epochs, fig_epochs = epochs_calculations( + epochs = epochs_calculations( raw_haemo, events, event_dict, @@ -5427,10 +5440,9 @@ def process_participant(file_path, file_start, progress_callback=None): t_max=T_MAX, baseline=(None,0), #TODO: Unhardcode this reject_epochs=REJECT_EPOCHS, - reject_hbo_threshold=dict(hbo=REJECT_HBO_THRESHOLD) + reject_hbo_threshold=dict(hbo=REJECT_HBO_THRESHOLD), + png_queue=png_queue ) - for name, fig in fig_epochs: - _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") @@ -5708,46 +5720,27 @@ def functional_connectivity_betas( # ------------------------------------------------------------------ beta_series = np.zeros((n_channels, len(trial_tags))) - for t, tag in enumerate(trial_tags): - idx = [ - i for i, col in enumerate(reg_names) - if col.startswith(f"{tag}_delay") - ] - 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)) - - # for t in range(n_trials): - # trial_indices = [i for i, col in enumerate(reg_names) if col.startswith(f"trial_{t:03d}_delay")] - # if trial_indices: - # beta_series[:, t] = np.mean(betas[:, trial_indices], axis=1).flatten() - - # Normalize each channel so they are on the same scale - # Without this, everything is connected to everything. Apparently this is a big issue in fNIRS? - beta_series = zscore(beta_series, axis=1) - - global_signal = np.mean(beta_series, axis=0) + global_signal = np.mean(beta_series, axis=0) beta_series_clean = np.zeros_like(beta_series) for i in range(n_channels): slope, _ = np.polyfit(global_signal, beta_series[i, :], 1) beta_series_clean[i, :] = beta_series[i, :] - (slope * global_signal) - - # 4. Correlation & Strict Filtering - corr_matrix = np.zeros((n_channels, n_channels)) - p_matrix = np.ones((n_channels, n_channels)) - - for i in range(n_channels): - for j in range(i + 1, n_channels): - r, p = pearsonr(beta_series_clean[i, :], beta_series_clean[j, :]) - corr_matrix[i, j] = corr_matrix[j, i] = r - p_matrix[i, j] = p_matrix[j, i] = p - - # 5. High-Bar Thresholding - reject, _ = multipletests(p_matrix[np.triu_indices(n_channels, k=1)], method='fdr_bh', alpha=0.05)[:2] - sig_corr_matrix = np.zeros_like(corr_matrix) + + # --- Vectorized correlation + analytic p-values (replaces the nested + # pearsonr loop below) --- + n_trials = beta_series_clean.shape[1] + corr_matrix = np.corrcoef(beta_series_clean) + + with np.errstate(divide='ignore', invalid='ignore'): + t_stats = corr_matrix * np.sqrt((n_trials - 2) / (1 - corr_matrix ** 2)) + p_matrix = 2 * t_dist.sf(np.abs(t_stats), df=n_trials - 2) + np.fill_diagonal(p_matrix, 1.0) # diagonal r=1 -> nan/inf guarded explicitly + triu = np.triu_indices(n_channels, k=1) + flat_p = p_matrix[triu] + + reject, _ = multipletests(flat_p, method='fdr_bh', alpha=0.05)[:2] + sig_corr_matrix = np.zeros_like(corr_matrix) for idx, is_sig in enumerate(reject): r_val = corr_matrix[triu[0][idx], triu[1][idx]]