fix rare bug

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%
+- 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 ~50%
+- 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 ~85%
+- 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

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):
+def epochs_calculations(raw_haemo, events, event_dict, epoch_handling, max_shift, t_min, t_max, baseline, reject_epochs, reject_hbo_threshold, png_queue):
-    fig_epochs = []  # List to store figures
 
     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 = 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
+        _enqueue(f"epochs_fig_evoked_hbo_{cond}", fig_evoked_hbo, png_queue)
-        fig_epochs.append((f"fig_evoked_hbr_{cond}", fig_evoked_hbr))  # Store with a unique name
+        _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):
+    df['_effect_weight'] = df['effect'] * df['weight']
-        return np.average(g['effect'], weights=g['weight'])
+    roi_theta = df.groupby(group_cols, as_index=False).agg(
- 
+        _sum_ew=('_effect_weight', 'sum'),
-    roi_theta = (
+        _sum_w=('weight', 'sum'),
-        df.groupby(group_cols, group_keys=False)
-        .apply(lambda g: pd.Series({'theta': _weighted_mean(g)}))
-        .reset_index()
     )
+    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 ---
+    f, Pxx = welch(data, fs=sfreq, nperseg=nperseg, axis=1)  # data: (n_channels, n_samples)
-    def find_hr_from_psd(ch_data):
+    mask = (f >= hr_range[0] / 60) & (f <= hr_range[1] / 60)
-        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]
+    Pxx_masked = Pxx[:, mask]  # (n_channels, n_freq_in_range)
-        if len(Pxx_masked) < 3:
+ 
-            return np.nan
+    hr_all_channels = np.full(Pxx_masked.shape[0], np.nan)
-        peaks = [i for i in range(1, len(Pxx_masked)-1) 
+    if Pxx_masked.shape[1] >= 3:
-                if Pxx_masked[i] > Pxx_masked[i-1] and Pxx_masked[i] > Pxx_masked[i+1]]
+        # same "strictly greater than both neighbors" local-max definition as
-        if not peaks:
+        # the original per-channel loop, vectorized across all channels at once
-            return np.nan
+        interior = Pxx_masked[:, 1:-1]
-        best_idx = peaks[np.argmax([Pxx_masked[i] for i in peaks])]
+        left = Pxx_masked[:, :-2]
-        return f_masked[best_idx] * 60  # bpm
+        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
  
-    # --- Compute HR across all channels ---
-    hr_all_channels = np.array([find_hr_from_psd(data[i, :]) for i in range(len(channel_names))])
     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)
+                col_idx = 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 theta_list]))
-                avg_absolute_theta = np.mean(np.abs([ch_theta[col_idx] for ch_theta in glm_est.theta()]))
                 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)
     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
+    # --- Vectorized correlation + analytic p-values (replaces the nested
-    corr_matrix = np.zeros((n_channels, n_channels))
+    #     pearsonr loop below) ---
-    p_matrix = np.ones((n_channels, n_channels))
+    n_trials = beta_series_clean.shape[1]
+    corr_matrix = np.corrcoef(beta_series_clean)
  
-    for i in range(n_channels):
+    with np.errstate(divide='ignore', invalid='ignore'):
-        for j in range(i + 1, n_channels):
+        t_stats = corr_matrix * np.sqrt((n_trials - 2) / (1 - corr_matrix ** 2))
-            r, p = pearsonr(beta_series_clean[i, :], beta_series_clean[j, :])
+    p_matrix = 2 * t_dist.sf(np.abs(t_stats), df=n_trials - 2)
-            corr_matrix[i, j] = corr_matrix[j, i] = r
+    np.fill_diagonal(p_matrix, 1.0)  # diagonal r=1 -> nan/inf guarded explicitly
-            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)
     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]]