From 06807183988e10d9b2874888e156f32767822f27 Mon Sep 17 00:00:00 2001 From: tyler Date: Fri, 17 Jul 2026 00:07:56 -0700 Subject: [PATCH] unhardcoding and general cleanup --- flares.py | 1471 +++++++++++++++++------ main.py | 30 +- src/analysis/crossgroupstats.py | 274 ++++- src/analysis/intergroupstats.py | 246 ++-- src/analysis/participantfoldchannels.py | 10 +- src/shared/flaresbasewidget.py | 134 ++- src/window/viewerlauncher.py | 4 +- 7 files changed, 1614 insertions(+), 555 deletions(-) diff --git a/flares.py b/flares.py index 858854b..3ef32ee 100644 --- a/flares.py +++ b/flares.py @@ -97,6 +97,23 @@ from mne_nirs.statistics._glm_level_first import RegressionResults # type: igno from mne_connectivity.viz import plot_connectivity_circle from mne_connectivity import envelope_correlation, spectral_connectivity_epochs, spectral_connectivity_time + +import os +import json +import warnings +import numpy as np +import pandas as pd +import scipy.stats as stats +import matplotlib.pyplot as plt +import seaborn as sns +import statsmodels.formula.api as smf +from statsmodels.stats.multitest import multipletests +from statsmodels.tools.sm_exceptions import ConvergenceWarning +from mne_nirs.statistics import statsmodels_to_results +from mne_nirs.visualisation import plot_glm_group_topo +import logging + + # Needs to be set for mne os.environ["SUBJECTS_DIR"] = str(data_path()) + "/subjects" # type: ignore @@ -104,7 +121,7 @@ PRIMARY_COLORS = { "SCI only": "skyblue", # Scalp Coupling Index (Standard MNE) "SNR only": "lightgreen", # Signal-to-Noise Ratio (Your original) "PSP only": "salmon", # Power Spectral Peak (Original Noise check) - "CV only": "yellow", # Relative Noise (The CV-only check) + "Coeff_var only": "yellow", # Relative Noise (The coeff_var-only check) "Range only": "coral", # Z-Swing (The Range Outlier check) "Noise only": "plum", # High-Freq PSD (The Noise check) "Disp. only": "palegreen", # Sensor Displacement (Variance Drop) @@ -112,6 +129,7 @@ PRIMARY_COLORS = { } COMBINATION_COLOR = "gray" +NUISANCE_EXCLUDE = ("drift", "constant", "short") def get_category_color(label): """Returns the primary color if it's a single failure, otherwise gray.""" @@ -127,9 +145,9 @@ SECONDS_TO_KEEP: float OPTODE_PLACEMENT: bool SHOW_OPTODE_NAMES: bool -SHORT_CHANNEL: bool -SHORT_CHANNEL_THRESH: float -LONG_CHANNEL_THRESH: float +SHORT_CHANNELS: bool +SHORT_CHANNELS_THRESHOLD: float +LONG_CHANNELS_THRESHOLD: float HEART_RATE: bool SECONDS_TO_STRIP_HR: int @@ -150,8 +168,8 @@ PSP: bool PSP_TIME_WINDOW: int PSP_THRESHOLD: float -CV: bool -CV_THRESHOLD: int +COEFF_VAR: bool +COEFF_VAR_THRESHOLD: int MAD: bool MAD_THRESHOLD: int @@ -160,8 +178,8 @@ PSD_NOISE: bool TARGET_FREQ_DIV: int DB_LIMIT: int -CHANNEL_VAR: bool -CHANNEL_THRESH: float +SENSOR_DROPOUT: bool +SENSOR_DROPOUT_VARIANCE_THRESHOLD: float BAD_CHANNELS_HANDLING: str MAX_DIST: float @@ -244,9 +262,9 @@ REQUIRED_KEYS: dict[str, Any] = { "PSP_TIME_WINDOW": int, "PSP_THRESHOLD": float, - "SHORT_CHANNEL": bool, - "SHORT_CHANNEL_THRESH": float, - "LONG_CHANNEL_THRESH": float, + "SHORT_CHANNELS": bool, + "SHORT_CHANNELS_THRESHOLD": float, + "LONG_CHANNELS_THRESHOLD": float, "REMOVE_EVENTS": list, @@ -1171,9 +1189,9 @@ def calculate_peak_power(data: BaseRaw, l_freq: float = 0.7, h_freq: float = 1.5 return list(compress(cast(list[str], getattr(data, "ch_names")), psp < PSP_THRESHOLD)), psp1, psp2 -def mark_bads(raw, bad_sci, bad_snr, bad_psp, bad_cv, bad_range, bad_noise, bad_disp): - print(bad_sci, bad_snr, bad_psp, bad_cv, bad_range, bad_noise, bad_disp) - bads_combined = list(set(bad_snr) | set(bad_sci) | set(bad_psp) | set(bad_cv) | set(bad_range) | set(bad_noise) | set(bad_disp)) +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)) print(f"Automatically marked bad channels based on SNR and SCI: {bads_combined}") raw.info['bads'].extend(bads_combined) @@ -1183,7 +1201,7 @@ def mark_bads(raw, bad_sci, bad_snr, bad_psp, bad_cv, bad_range, bad_noise, bad (bad_sci, "SCI"), (bad_psp, "PSP"), (bad_snr, "SNR"), - (bad_cv, "CV"), + (bad_coeff_var, "coeff_var"), (bad_range, "Range"), (bad_noise, "Noise"), (bad_disp, "Disp.") @@ -1468,7 +1486,7 @@ def epochs_calculations(raw_haemo, events, event_dict): -def make_design_matrix(raw_haemo, short_chans): +def make_design_matrix(raw_haemo): events_to_remove = REMOVE_EVENTS @@ -1480,6 +1498,12 @@ def make_design_matrix(raw_haemo, short_chans): description=[ann['description'] for ann in filtered_annotations] ) + if SHORT_CHANNELS: + short_chans = get_short_channels(raw_haemo, max_dist=SHORT_CHANNELS_THRESHOLD) + raw_haemo = get_long_channels(raw_haemo, min_dist=SHORT_CHANNELS_THRESHOLD, max_dist=LONG_CHANNELS_THRESHOLD) + else: + short_chans = None + # Set the new annotations raw_haemo.set_annotations(new_annot) @@ -1550,7 +1574,7 @@ def make_design_matrix(raw_haemo, short_chans): fig, ax1 = plt.subplots(figsize=(10, 6), constrained_layout=True) _ = plot_design_matrix(design_matrix, axes=ax1) - return design_matrix, fig + return raw_haemo, design_matrix, fig @@ -2763,20 +2787,7 @@ def load_snirf(file_path: str) -> tuple[BaseRaw, Figure]: # Data science requires good statistics. This section is requiring work. Am working on figuring this out. Will likely add this to a new viewer. -import os -import json -import warnings -import numpy as np -import pandas as pd -import scipy.stats as stats -import matplotlib.pyplot as plt -import seaborn as sns -import statsmodels.formula.api as smf -from statsmodels.stats.multitest import multipletests -from statsmodels.tools.sm_exceptions import ConvergenceWarning -from mne_nirs.statistics import statsmodels_to_results -from mne_nirs.visualisation import plot_glm_group_topo -import logging + def run_roi_second_level_analysis(df_roi_all, df_cha_all=None, raw_haemo=None, @@ -3024,19 +3035,6 @@ def run_roi_second_level_analysis(df_roi_all, df_cha_all=None, raw_haemo=None, -import os -import json -import warnings -import numpy as np -import pandas as pd -import scipy.stats as stats -import matplotlib.pyplot as plt -import seaborn as sns -from statsmodels.stats.multitest import multipletests -from mne_nirs.visualisation import plot_glm_group_topo -import logging - -logger = logging.getLogger(__name__) def clean_subject_id(path_or_id): """ @@ -3322,17 +3320,471 @@ def run_cross_group_second_level_analysis(df_roi_all, file_paths_a, file_paths_b +def run_cross_group_laterality_analysis(df_roi_all_a, df_roi_all_b, roi_pairs, condition, + group_a_name="Group A", group_b_name="Group B", + target_chroma='hbo', min_subjects=3, + p_threshold=0.05, correction_method=None, + roi_contra_label=None, roi_ipsi_label=None): + """ + Compare LATERALITY between two independent groups of subjects (e.g. a + control group vs. a target group), using Welch's t-test on each + subject's within-subject laterality index rather than on raw ROI values. + + -------------------------------------------------------------------- + HOW THIS DIFFERS FROM run_cross_group_second_level_analysis + -------------------------------------------------------------------- + run_cross_group_second_level_analysis (existing): + - Compares one ROI's raw theta value between two groups directly + (Group A's Right_PFC vs Group B's Right_PFC, say). + - CLAIM IF SIGNIFICANT: this ROI's response magnitude differs between + the two populations, for this condition. + - WHAT IT DOES NOT SAY: whether that difference reflects a real, + localized neural difference or a generic between-population + difference unrelated to the specific task — e.g. different overall + vascular reactivity, arousal, skull/scalp optical properties, or + anything else that would shift a group's numbers up or down + everywhere, not just in this ROI. Two independently recruited + groups (e.g. patients vs. healthy controls) are considerably more + likely to differ in these generic ways than two subsets of one + study population, which makes this ambiguity a real risk here, not + a theoretical one. + + run_cross_group_laterality_analysis (this function): + - First computes each subject's OWN laterality index + (contralateral ROI theta - ipsilateral ROI theta, within that + subject, for one condition) — the same computation as + run_roi_paired_contrast_analysis, just not yet tested there. + - Then compares those per-subject laterality indices between the two + groups with Welch's t-test. + - CLAIM IF SIGNIFICANT: the DEGREE OF SPATIAL SPECIFICITY (how much + more one hemisphere responds than the other, within a person) + differs between the two groups — a claim about lateralization + itself, not raw magnitude. Subtracting within-subject first cancels + out whatever's common to both hemispheres for that person (general + reactivity, arousal, etc.) before ever comparing across groups, so + a significant result here is harder to explain away as a generic + between-population confound. + - WHAT IT DOES NOT SAY: anything about whether overall response + magnitude differs between groups (a group could have identical + laterality but very different raw amplitude — that's what the + existing cross-group function is for) — and it only uses subjects + who have BOTH the contra and ipsi ROI valid, so it can lose + subjects the raw-ROI comparison would have kept. + + Use both, for different questions: the existing function for "is the + raw response different between groups," this one for "is the + LATERALIZATION different between groups." + + Parameters + ---------- + df_roi_all_a, df_roi_all_b : pd.DataFrame + Individual-level ROI results (['ROI', 'Condition', 'Chroma', + 'theta', 'ID']) for Group A and Group B RESPECTIVELY. Keeping them + as separate frames (rather than one combined frame + ID lists) + avoids any risk of cross-dataset ID collisions when the two groups + come from genuinely separate studies/exports. + roi_pairs : tuple(str, str) or list of tuple(str, str) + (roi_contra, roi_ipsi) pair(s). Each pair's laterality index is + computed as theta(roi_contra) - theta(roi_ipsi), per subject. + Pass a list to test multiple hand/condition combinations in one call. + condition : str or list of str + The 'Condition' value (e.g. contrast name or event code) to use for + each pair. Single value applies to all pairs; otherwise must match + len(roi_pairs). + target_chroma : str, default 'hbo' + Chromophore to test. Never mix hbo/hbr in one laterality index. + min_subjects : int, default 3 + Minimum subjects required in EACH group (after requiring both ROIs + be present) for a pair to be tested. + p_threshold : float, default 0.05 + Significance threshold for the (optionally corrected) p-value. + correction_method : str or None, default None + Multiple comparisons correction across the pairs tested in this + call. Off by default for a single pre-specified pair; turn on + ('fdr_bh') if testing several pairs/conditions at once. + roi_contra_label, roi_ipsi_label : str or list of str, optional + Display labels for the contra/ipsi ROI in each pair. + + Returns + ------- + pd.DataFrame, one row per tested pair: + ['roi_contra', 'roi_ipsi', 'condition', 'mean_A', 'mean_B', + 'mean_diff', 't_val', 'p_val', 'p_corrected', 'significant', + 'n_A', 'n_B'] + """ + + required_cols = ['ROI', 'Condition', 'Chroma', 'theta', 'ID'] + for name, df in [('df_roi_all_a', df_roi_all_a), ('df_roi_all_b', df_roi_all_b)]: + if not all(col in df.columns for col in required_cols): + raise ValueError(f"{name} must include: {required_cols}") + + if isinstance(roi_pairs, tuple): + roi_pairs = [roi_pairs] + n_pairs = len(roi_pairs) + + if isinstance(condition, str): + conditions = [condition] * n_pairs + else: + if len(condition) != n_pairs: + raise ValueError("If passing a list of conditions, it must match len(roi_pairs).") + conditions = list(condition) + + def _expand(labels): + if labels is None: + return [None] * n_pairs + if isinstance(labels, str): + return [labels] * n_pairs + if len(labels) != n_pairs: + raise ValueError("Label list length must match len(roi_pairs).") + return list(labels) + + contra_labels = _expand(roi_contra_label) + ipsi_labels = _expand(roi_ipsi_label) + + def _laterality_per_subject(df_roi_all, roi_contra, roi_ipsi, cond): + """Collapse to one laterality index per subject, for one group.""" + df_chroma = df_roi_all[df_roi_all['Chroma'] == target_chroma] + df_cond = df_chroma[df_chroma['Condition'] == cond] + + contra_vals = df_cond[df_cond['ROI'] == roi_contra].groupby('ID', as_index=False)['theta'].mean() + ipsi_vals = df_cond[df_cond['ROI'] == roi_ipsi].groupby('ID', as_index=False)['theta'].mean() + + merged = contra_vals.merge(ipsi_vals, on='ID', suffixes=('_contra', '_ipsi')) + merged['laterality'] = merged['theta_contra'] - merged['theta_ipsi'] + return merged[['ID', 'laterality']] + + results = [] + plot_rows = [] + + for (roi_contra, roi_ipsi), cond, lbl_c, lbl_i in zip(roi_pairs, conditions, contra_labels, ipsi_labels): + lat_a = _laterality_per_subject(df_roi_all_a, roi_contra, roi_ipsi, cond) + lat_b = _laterality_per_subject(df_roi_all_b, roi_contra, roi_ipsi, cond) + + n_a, n_b = lat_a['ID'].nunique(), lat_b['ID'].nunique() + if n_a < min_subjects or n_b < min_subjects: + logger.warning( + f"Skipping pair ({roi_contra} - {roi_ipsi}) for condition '{cond}' — " + f"{group_a_name} n={n_a}, {group_b_name} n={n_b}, need at least {min_subjects} in EACH." + ) + continue + + vals_a = lat_a['laterality'].values + vals_b = lat_b['laterality'].values + + t_val, p_val = stats.ttest_ind(vals_a, vals_b, equal_var=False) + mean_a, mean_b = np.mean(vals_a), np.mean(vals_b) + mean_diff = mean_a - mean_b + + pair_label = f"{lbl_c or roi_contra} - {lbl_i or roi_ipsi}\n({cond})" + results.append({ + 'roi_contra': roi_contra, + 'roi_ipsi': roi_ipsi, + 'label': pair_label, + 'condition': cond, + 'mean_A': mean_a, + 'mean_B': mean_b, + 'mean_diff': mean_diff, + 't_val': t_val, + 'p_val': p_val, + 'n_A': n_a, + 'n_B': n_b, + }) + plot_rows.append(lat_a.assign(pair=pair_label, Group=group_a_name)) + plot_rows.append(lat_b.assign(pair=pair_label, Group=group_b_name)) + + if not results: + print("\n[ERROR] No ROI pairs met the minimum subject threshold for BOTH groups.\n") + return pd.DataFrame() + + df_group = pd.DataFrame(results) + + if correction_method is not None: + reject, p_corrected, _, _ = multipletests(df_group['p_val'].values, method=correction_method) + df_group['p_corrected'] = p_corrected + df_group['significant'] = reject + else: + df_group['p_corrected'] = df_group['p_val'] + df_group['significant'] = df_group['p_val'] <= p_threshold + + # --- Print report --- + print("\n" + "=" * 85) + print(f" CROSS-GROUP LATERALITY CONTRAST: {group_a_name.upper()} vs {group_b_name.upper()} ({target_chroma.upper()})") + print("=" * 85) + df_print = df_group.copy() + for c in ['mean_A', 'mean_B', 'mean_diff']: + df_print[c] = df_print[c].apply(lambda x: f"{x:.4f}") + df_print['t_val'] = df_print['t_val'].apply(lambda x: f"{x:.3f}") + df_print['p_corrected'] = df_print['p_corrected'].apply(lambda x: f"{x:.4f}") + print(df_print[['label', 'condition', 'mean_A', 'mean_B', 'mean_diff', + 't_val', 'p_corrected', 'significant', 'n_A', 'n_B']].to_string(index=False)) + print("=" * 85 + "\n") + + # --- Plot: grouped bar + swarm per pair, Group A vs Group B --- + sns.set_theme(style="whitegrid") + fig, ax = plt.subplots(figsize=(max(7, 3 * len(results)), 6)) + + plot_df = pd.concat(plot_rows, ignore_index=True) + pair_order = [r['label'] for r in results] + + sns.barplot( + data=plot_df, x='pair', y='laterality', hue='Group', + order=pair_order, hue_order=[group_a_name, group_b_name], + ax=ax, errorbar=('ci', 95), capsize=0.08, + palette=['#2b5c8f', '#d95f02'], edgecolor='black', linewidth=1.5, zorder=1 + ) + sns.swarmplot( + data=plot_df, x='pair', y='laterality', hue='Group', + order=pair_order, hue_order=[group_a_name, group_b_name], + ax=ax, size=6, palette=['black', 'black'], alpha=0.5, dodge=True, zorder=2, + legend=False + ) + + ax.axhline(0, color='black', linewidth=1, linestyle='--') + + global_max = plot_df['laterality'].max() + for i, row in df_group.iterrows(): + pair_points = plot_df[plot_df['pair'] == row['label']]['laterality'] + max_y = pair_points.max() if len(pair_points) > 0 else 0 + text_y = max_y + (abs(global_max) * 0.1 if global_max else 0.1) + + p_val_corr = row['p_corrected'] + if p_val_corr < p_threshold: + sig_symbol = "***" if p_val_corr < 0.001 else "**" if p_val_corr < 0.01 else "*" + ax.text(i, text_y, f"{sig_symbol}\np = {p_val_corr:.3f}", + ha='center', va='bottom', fontsize=10, fontweight='bold', color='red') + else: + ax.text(i, text_y, f"n.s.\np = {p_val_corr:.3f}", + ha='center', va='bottom', fontsize=10, color='gray') + + ax.set_ylabel(r'Laterality Index ($\Delta$ HbO, Contra $-$ Ipsi)', fontsize=12) + ax.set_xlabel('') + correction_lbl = f"({correction_method} corrected)" if correction_method else "(uncorrected)" + ax.set_title( + f"Cross-Group Laterality Comparison ({target_chroma.upper()})\n" + f"{group_a_name} vs {group_b_name}, p < {p_threshold} {correction_lbl}", + fontsize=13, fontweight='bold', pad=15 + ) + plt.tight_layout() + plt.show() + + return df_group -import numpy as np -import pandas as pd -import scipy.stats as stats -import matplotlib.pyplot as plt -import seaborn as sns -from statsmodels.stats.multitest import multipletests -import logging -logger = logging.getLogger(__name__) +def run_cross_group_contrast_analysis(df_contrasts_a, df_contrasts_b, contrast_name, roi_json_path, + group_a_name="Group A", group_b_name="Group B", + target_chroma='hbo', min_subjects=3, + p_threshold=0.05, correction_method='fdr_bh', + weighted=True): + """ + Compare a JOINT-FIT TASK CONTRAST (e.g. '2.0_vs_3.0'), aggregated to ROI + level, between two independent groups. This is the cross-group analog + of the inter-group joint-contrast method — where that method asks "does + this contrast differ from zero within one group," this asks "does the + SIZE of this contrast differ between two groups." + + -------------------------------------------------------------------- + HOW THIS DIFFERS FROM THE OTHER TWO CROSS-GROUP METHODS + -------------------------------------------------------------------- + run_cross_group_second_level_analysis: compares raw single-condition + ROI magnitude between groups — vulnerable to generic between- + population differences (vascular reactivity, arousal, etc.) that + have nothing to do with the task. + run_cross_group_laterality_analysis: compares each subject's own + contra-minus-ipsi laterality index between groups — asks whether + spatial specificity differs, says nothing about overall magnitude. + run_cross_group_contrast_analysis (this function): compares a + jointly-fit task contrast (e.g. Task A minus Task B, estimated + together within each subject's GLM) between groups — asks whether + one group differentiates between the two tasks more/less than the + other does, at this ROI. Systemic noise is cancelled at the + model-fitting stage (same GLM, both conditions) rather than left in + raw single-condition magnitude, or cancelled only by within-subject + spatial subtraction as in the laterality method. This is generally + the most statistically efficient of the three at detecting a real + between-group difference in TASK-SPECIFIC response, but — like the + inter-group version of this same idea — it does not by itself tell + you WHERE that difference is spatially localized unless you compare + the sign/pattern across multiple ROIs. + + Parameters + ---------- + df_contrasts_a, df_contrasts_b : pd.DataFrame + Combined CHANNEL-LEVEL contrast results (contrasts.csv format) for + Group A and Group B respectively. Must include: + ['ch_name', 'effect', 'stat', 'Chroma', 'contrast_name', 'ID'] + Kept as separate frames per group (not one combined frame + ID + lists) for the same reason as run_cross_group_laterality_analysis — + avoids any risk of ID-matching mismatches between two genuinely + separate dataset exports. + contrast_name : str + Which contrast to test (e.g. '2.0_vs_3.0'). Must exist in both + groups' df_contrasts for a fair comparison. + roi_json_path : str + Path to the regions.json used elsewhere in the pipeline. + target_chroma : str, default 'hbo' + min_subjects : int, default 3 + Minimum subjects required in EACH group, per ROI. + p_threshold : float, default 0.05 + correction_method : str or None, default 'fdr_bh' + Unlike the paired/laterality functions (which default to no + correction, since they test one pre-specified pair), this defaults + ON — this function screens across every ROI in regions.json, which + is an open multiple-comparisons scan, not a single planned contrast. + weighted : bool, default True + Passed through to aggregate_channel_contrasts_to_roi (inverse- + variance weighting vs. plain mean across channels within an ROI). + + Returns + ------- + pd.DataFrame, one row per ROI: + ['ROI', 'mean_A', 'mean_B', 'mean_diff', 't_val', 'p_val', + 'p_corrected', 'significant', 'n_A', 'n_B'] + """ + + required_cols = ['ch_name', 'effect', 'stat', 'Chroma', 'contrast_name', 'ID'] + for name, df in [('df_contrasts_a', df_contrasts_a), ('df_contrasts_b', df_contrasts_b)]: + if not all(col in df.columns for col in required_cols): + raise ValueError(f"{name} must include: {required_cols}") + + # Filter to the requested contrast BEFORE aggregating, so a missing + # contrast_name fails clearly here rather than silently downstream. + df_a_filt = df_contrasts_a[df_contrasts_a['contrast_name'] == contrast_name] + df_b_filt = df_contrasts_b[df_contrasts_b['contrast_name'] == contrast_name] + + if df_a_filt.empty: + print(f"[ERROR] Contrast '{contrast_name}' not found anywhere in {group_a_name}'s data.") + return pd.DataFrame() + if df_b_filt.empty: + print(f"[ERROR] Contrast '{contrast_name}' not found anywhere in {group_b_name}'s data.") + return pd.DataFrame() + + roi_a = aggregate_channel_contrasts_to_roi(df_a_filt, roi_json_path, weighted=weighted) + roi_b = aggregate_channel_contrasts_to_roi(df_b_filt, roi_json_path, weighted=weighted) + + roi_a = roi_a[roi_a['Chroma'] == target_chroma] + roi_b = roi_b[roi_b['Chroma'] == target_chroma] + + if roi_a.empty or roi_b.empty: + print(f"[ERROR] No ROI-aggregated values produced for one or both groups " + f"(check regions.json channel names against this montage).") + return pd.DataFrame() + + all_rois = sorted(set(roi_a['ROI'].unique()) | set(roi_b['ROI'].unique())) + results = [] + plot_rows = [] + + for roi in all_rois: + vals_a = roi_a[roi_a['ROI'] == roi]['theta'].values + vals_b = roi_b[roi_b['ROI'] == roi]['theta'].values + + n_a, n_b = len(vals_a), len(vals_b) + if n_a < min_subjects or n_b < min_subjects: + logger.warning( + f"Skipping ROI '{roi}' — {group_a_name} n={n_a}, {group_b_name} n={n_b}, " + f"need at least {min_subjects} in EACH." + ) + continue + + t_val, p_val = stats.ttest_ind(vals_a, vals_b, equal_var=False) + mean_a, mean_b = np.mean(vals_a), np.mean(vals_b) + mean_diff = mean_a - mean_b + + results.append({ + 'ROI': roi, 'mean_A': mean_a, 'mean_B': mean_b, 'mean_diff': mean_diff, + 't_val': t_val, 'p_val': p_val, 'n_A': n_a, 'n_B': n_b, + }) + plot_rows.append(pd.DataFrame({'theta': vals_a, 'ROI': roi, 'Group': group_a_name})) + plot_rows.append(pd.DataFrame({'theta': vals_b, 'ROI': roi, 'Group': group_b_name})) + + if not results: + print("\n[ERROR] No ROIs met the subject requirements for BOTH groups.\n") + return pd.DataFrame() + + df_group = pd.DataFrame(results) + + if correction_method is not None: + reject, p_corrected, _, _ = multipletests(df_group['p_val'].values, method=correction_method) + df_group['p_corrected'] = p_corrected + df_group['significant'] = reject + else: + df_group['p_corrected'] = df_group['p_val'] + df_group['significant'] = df_group['p_val'] <= p_threshold + + # --- Print report --- + print("\n" + "=" * 85) + print(f" CROSS-GROUP CONTRAST COMPARISON: {group_a_name.upper()} vs {group_b_name.upper()} ({target_chroma.upper()})") + print(f" Contrast: {contrast_name}") + print("=" * 85) + df_print = df_group.copy() + for c in ['mean_A', 'mean_B', 'mean_diff']: + df_print[c] = df_print[c].apply(lambda x: f"{x:.4f}") + df_print['t_val'] = df_print['t_val'].apply(lambda x: f"{x:.3f}") + df_print['p_corrected'] = df_print['p_corrected'].apply(lambda x: f"{x:.4f}") + print(df_print[['ROI', 'mean_A', 'mean_B', 'mean_diff', 't_val', + 'p_corrected', 'significant', 'n_A', 'n_B']].to_string(index=False)) + print("=" * 85 + "\n") + + # --- Plot: grouped bar + swarm per ROI, Group A vs Group B, with brackets --- + sns.set_theme(style="whitegrid") + fig, ax = plt.subplots(figsize=(max(7, 2.5 * len(results)), 6)) + + plot_df = pd.concat(plot_rows, ignore_index=True) + roi_order = [r['ROI'] for r in results] + + sns.barplot( + data=plot_df, x='ROI', y='theta', hue='Group', + order=roi_order, hue_order=[group_a_name, group_b_name], + ax=ax, errorbar=('ci', 95), capsize=0.08, + palette=['#2b5c8f', '#d95f02'], edgecolor='black', linewidth=1.5, zorder=1 + ) + sns.swarmplot( + data=plot_df, x='ROI', y='theta', hue='Group', + order=roi_order, hue_order=[group_a_name, group_b_name], + ax=ax, size=6, palette=['black', 'black'], alpha=0.5, dodge=True, zorder=2, + legend=False + ) + + global_max = plot_df['theta'].max() + global_min = plot_df['theta'].min() + y_top = global_max * 1.45 if global_max > 0 else 0.5e-6 + y_bottom = global_min * 1.15 if global_min < 0 else -0.15 * global_max + ax.set_ylim(y_bottom, y_top) + + for i, row in df_group.iterrows(): + roi_points = plot_df[plot_df['ROI'] == row['ROI']]['theta'] + max_y = roi_points.max() if len(roi_points) > 0 else 0 + x_a, x_b = i - 0.2, i + 0.2 + y_bracket = max_y + (global_max * 0.08 if global_max else 0.05) + h_tick = global_max * 0.02 if global_max else 0.01 + + p_val_corr = row['p_corrected'] + if p_val_corr < p_threshold: + sig_symbol = "***" if p_val_corr < 0.001 else "**" if p_val_corr < 0.01 else "*" + ax.plot([x_a, x_a, x_b, x_b], + [y_bracket - h_tick, y_bracket, y_bracket, y_bracket - h_tick], + color='black', lw=1.2) + ax.text(i, y_bracket + (global_max * 0.02 if global_max else 0.01), + f"{sig_symbol}\np_corr = {p_val_corr:.3f}", + ha='center', va='bottom', fontsize=9, fontweight='bold', color='red') + else: + ax.text(i, y_bracket, "n.s.", ha='center', va='bottom', fontsize=9, color='gray') + + ax.axhline(0, color='black', linewidth=1, linestyle='--') + ax.set_ylabel(r'Contrast Effect ($\Delta$ HbO)', fontsize=12) + ax.set_xlabel('Region of Interest (ROI)', fontsize=12) + correction_lbl = f"({correction_method} corrected)" if correction_method else "(uncorrected)" + ax.set_title( + f"Cross-Group Contrast Comparison: {group_a_name} vs {group_b_name}\n" + f"({target_chroma.upper()} - {contrast_name}) {correction_lbl}", + fontsize=13, fontweight='bold', pad=15 + ) + plt.tight_layout() + plt.show() + + return df_group + + def run_roi_paired_contrast_analysis(df_roi_all, roi_pairs, condition, @@ -3662,6 +4114,214 @@ def aggregate_channel_contrasts_to_roi(df_contrasts, roi_json_path, weighted=Tru + + +def collapse_fir_condition_column(df, value_col, condition_col='Condition', + group_cols=None, delay_sep='_delay_'): + """ + Collapse FIR delay-bin rows (Condition values like '2.0_delay_6') into a + single row per base condition ('2.0'), using a PLAIN (equal-weighted) + mean of `value_col` across whatever delay bins are actually present for + that condition/subject/unit. SAFE NO-OP for non-FIR data: if no + Condition value contains `delay_sep`, the input is returned unchanged — + so this can always be called unconditionally, regardless of HRF_MODEL. + + No window (start/end) parameter, deliberately: the set of delay bins + to average is read directly from whichever bins actually exist in the + data for that condition — which is itself just a reflection of + FIR_DELAYS/the design matrix that was actually built — rather than a + second, separately-configured window that could drift out of sync with + it. This also makes no assumption about response shape or timing, + appropriate when response latency is unpredictable (e.g. infant fNIRS). + + This is what lets all six group-level statistics functions work + identically whether the underlying GLM used HRF_MODEL='glover'/'spm' + (one regressor per condition already) or 'fir' (many delay-bin + regressors per condition) — the delay-bin collapsing happens once, + upstream, rather than needing separate handling inside each function. + + Parameters + ---------- + df : pd.DataFrame + Must include `condition_col` and `value_col`, plus whatever + `group_cols` you want preserved (e.g. ['ROI', 'Chroma', 'ID'] or + ['ch_name', 'Chroma', 'ID']). + value_col : str + Column to average (e.g. 'theta' for ROI data, 'effect' for + channel-level data). + condition_col : str, default 'Condition' + Column holding condition/delay-bin labels. + group_cols : list of str, optional + Columns that define one "unit" to collapse within (e.g. one + ROI/subject, or one channel/subject). Required whenever FIR rows + are present, to avoid accidentally collapsing across subjects/ROIs. + delay_sep : str, default '_delay_' + Separator used in FIR condition names, matching the naming already + used elsewhere in the pipeline ("{condition}_delay_{n}"). + + Returns + ------- + pd.DataFrame with the same columns as the input, `condition_col` + holding base condition names, and one row per (group_cols, base + condition) combination. + """ + df = df.copy() + is_fir_row = df[condition_col].astype(str).str.contains(delay_sep, regex=False) + + if not is_fir_row.any(): + # Nothing FIR-shaped here — pass through unchanged. + return df + + if group_cols is None: + raise ValueError( + "group_cols must be specified when FIR delay-bin rows are present " + "(e.g. ['ROI', 'Chroma', 'ID'] or ['ch_name', 'Chroma', 'ID']) — " + "otherwise rows from different subjects/ROIs/channels could be " + "collapsed together incorrectly." + ) + + df_fir = df[is_fir_row].copy() + df_other = df[~is_fir_row].copy() # any non-FIR rows pass through untouched + + df_fir['_base_condition'] = df_fir[condition_col].astype(str).str.split(delay_sep, n=1).str[0] + + full_group = group_cols + ['_base_condition'] + + collapsed = ( + df_fir.groupby(full_group, as_index=False)[value_col] + .mean() + .rename(columns={'_base_condition': condition_col}) + ) + + if not df_other.empty: + collapsed = pd.concat([collapsed, df_other], ignore_index=True) + + return collapsed + + + +def _channel_midpoint(ch_info): + """(x, y, z) midpoint between source and detector for one channel entry + from raw_haemo.info['chs']. Returns None if location info is missing/ + degenerate (e.g. an aux/stim channel with no real optode geometry). + MNE head-coordinate convention: x = left(-)/right(+), + y = posterior(-)/anterior(+), z = inferior(-)/superior(+).""" + loc = ch_info['loc'] + if loc is None or not np.asarray(loc).any(): + return None + 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 + channel's source-detector midpoint. + + Parameters + ---------- + axis : int + 0 = x (left/right), 1 = y (posterior/anterior). z (axis 2, superior/ + inferior) isn't offered as a fallback split — depth splits aren't a + meaningful functional distinction the way left/right or front/back + are for a 2D optode array. + names : (str, str) + Names for the (negative-side, positive-side) ROIs. + balance_threshold : float, default 0.5 + Minimum acceptable ratio of (smaller side size / larger side size). + 0.5 means the smaller side must be at least half the size of the + larger — rejects near-degenerate splits (e.g. 27 channels on one + side, 1 on the other) where "two ROIs" isn't really giving you two + usable regions, without demanding a perfect, unrealistic 50/50. + + Returns + ------- + dict of two ROIs, or None if geometry is missing/degenerate, every + channel falls on one side, or the split is too imbalanced to be useful + — signaling the caller to try the next fallback tier. + """ + neg_indices, pos_indices = [], [] + + for idx, ch in enumerate(raw_haemo.info['chs']): + mid = _channel_midpoint(ch) + if mid is None: + continue + coord = mid[axis] + if coord < 0: + neg_indices.append(idx) + elif coord > 0: + pos_indices.append(idx) + # coord == 0 (exact midline) excluded from both — genuinely + # ambiguous, not worth guessing a side for. + + if not neg_indices or not pos_indices: + return None + + balance = min(len(neg_indices), len(pos_indices)) / max(len(neg_indices), len(pos_indices)) + if balance < balance_threshold: + logger.warning( + f"Axis-{axis} split too imbalanced ({len(neg_indices)} vs " + f"{len(pos_indices)}, ratio {balance:.2f} < {balance_threshold}) — rejecting." + ) + 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 + (the most common and most interpretable axis for bilateral montages); + if that's unavailable or too imbalanced (e.g. a montage covering only + one cortical region, where every channel falls on the same side), try + a Front/Back split instead, using the exact same geometry. Returns None + if neither axis gives a usable split, signaling the caller to fall back + further to one-ROI-per-channel. + + Like the hemisphere-only version this replaces, these are coarse + geometric splits, not hand-drawn functional regions — labeled + "_Auto" so they're never mistaken for real regions.json ROIs anywhere + downstream (tables, plot titles, exported CSVs). + """ + lr = _build_axis_split_rois(raw_haemo, axis=0, names=("Left_Auto", "Right_Auto")) + if lr is not None: + logger.info("Automatic fallback ROIs: Left/Right split (axis available and balanced).") + return lr + + logger.warning("Left/Right fallback unavailable or too imbalanced — trying Front/Back split.") + fb = _build_axis_split_rois(raw_haemo, axis=1, names=("Back_Auto", "Front_Auto")) + if fb is not None: + logger.info("Automatic fallback ROIs: Front/Back split (Left/Right was not usable).") + return fb + + 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 + pair), each containing that channel's hbo AND hbr indices together — + same grouping convention as regions.json (one name -> both + chromophores), just derived automatically from whatever channels exist. + + NOT good for statistics — every "region" is a single channel, so this + provides none of ROI aggregation's noise-reduction or multiple- + comparisons benefit. Only reached if regions.json AND both automatic + geometric splits are unavailable, so processing degrades gracefully to + single-channel resolution rather than crashing, or silently averaging + unrelated regions together into one meaningless number the way a + single AllChannels ROI would. + """ + rois_formatted = {} + for ch_name in raw_haemo.ch_names: + base_name = ch_name.split()[0] # "S1_D1 hbo" -> "S1_D1" + idx = raw_haemo.ch_names.index(ch_name) + rois_formatted.setdefault(base_name, []).append(idx) + return rois_formatted + + @@ -4198,7 +4858,7 @@ def find_flatline_at_end(raw, threshold_ratio=0.05): -def detect_sensor_displacement(raw, threshold_ratio=0.05): +def detect_sensor_dropout(raw, threshold_ratio=0.05): """ Identifies channels where signal variance drops significantly. Returns flagged channel names and a summary figure. @@ -4240,7 +4900,7 @@ def detect_sensor_displacement(raw, threshold_ratio=0.05): ax.bar(range(len(ratios)), ratios, color=colors) ax.axhline(threshold_ratio, color='red', linestyle='--', label=f'Threshold ({threshold_ratio:.0%})') - ax.set_title("Sensor Displacement Check (Variance Stability)") + ax.set_title("Sensor Dropout Check (Variance Stability)") ax.set_ylabel("Variance Ratio (End / Start)") ax.set_xlabel("Channel Index") ax.set_ylim(0, max(ratios + [threshold_ratio * 2])) # Scale to see the threshold clearly @@ -4248,12 +4908,12 @@ def detect_sensor_displacement(raw, threshold_ratio=0.05): plt.close(fig_disp) - print(f"Displacement Check: Flagged {len(failed_bases)} optode pairs.") + print(f"Dropout Check: Flagged {len(failed_bases)} optode pairs.") return bad_names, fig_disp -def detect_high_freq_noise(raw, db_limit=-60, freq_div=4): +def detect_spectral_noise_spike(raw, db_limit=-60, freq_div=4): """ Identifies channels with excessive power at high frequencies (sfreq/4), usually indicating electronic interference. @@ -4294,7 +4954,8 @@ def detect_high_freq_noise(raw, db_limit=-60, freq_div=4): -def find_bad_channels_range(raw, threshold=4.0): +def find_bad_channels_by_amplitude_range(raw, threshold=4.0): + """Median absolute deviation""" picks = [ch for ch in raw.ch_names] data = raw.get_data(picks=picks) ranges = np.max(data, axis=1) - np.min(data, axis=1) @@ -4329,52 +4990,52 @@ def find_bad_channels_range(raw, threshold=4.0): -def find_bad_channels_cv(raw, cv_threshold=25.0): +def find_bad_channels_coeff_var(raw, coeff_var_threshold=25.0): """ - Identifies bad fNIRS channels using only the Coefficient of Variation (CV). + Identifies bad fNIRS channels using only the Coefficient of Variation (coeff_var). """ - print(f"\n--- Starting CV-Only Quality Check on the channels ---") + print(f"\n--- Starting coeff_var-Only Quality Check on the channels ---") picks = [ch for ch in raw.ch_names] data = raw.get_data(picks=picks) - # Calculate CV (Coefficient of Variation) + # Calculate coeff_var (Coefficient of Variation) stds = np.std(data, axis=1) means = np.mean(data, axis=1) # Using a small epsilon (1e-15) to prevent division by zero - cv_scores = (stds / (means + 1e-15)) * 100 + coeff_var_scores = (stds / (means + 1e-15)) * 100 # Find indices that exceed the threshold - bad_cv_indices = np.where(cv_scores > cv_threshold)[0] + bad_coeff_var_indices = np.where(coeff_var_scores > coeff_var_threshold)[0] # Pair-kill logic: If one wavelength (HbO or HbR) fails, flag the pair failed_bases = set() - for idx in bad_cv_indices: + for idx in bad_coeff_var_indices: base = picks[idx].split(' ')[0] failed_bases.add(base) bad_names = [ch for ch in picks if ch.split(' ')[0] in failed_bases] # Summary Prints - print(f"CV Check: Found {len(bad_cv_indices)} channels exceeding {cv_threshold}% noise threshold.") + print(f"coeff_var Check: Found {len(bad_coeff_var_indices)} channels exceeding {coeff_var_threshold}% noise threshold.") if failed_bases: print(f"Flagged {len(failed_bases)} optode pairs for removal:") for base in sorted(failed_bases): - # Find the specific CV for this base (using the first channel found for it) + # Find the specific coeff_var for this base (using the first channel found for it) ch_idx = picks.index(next(p for p in picks if p.startswith(base))) - print(f" - {base}: CV = {cv_scores[ch_idx]:.2f}%") + print(f" - {base}: coeff_var = {coeff_var_scores[ch_idx]:.2f}%") else: - print("All channels passed the CV check.") + print("All channels passed the coeff_var check.") # --- Visualization --- fig_qc, ax = plt.subplots(figsize=(10, 5), constrained_layout=True) - colors = ['coral' if c > cv_threshold else 'skyblue' for c in cv_scores] - ax.bar(range(len(cv_scores)), cv_scores, color=colors) - ax.axhline(cv_threshold, color='red', linestyle='--', label=f'Threshold ({cv_threshold}%)') + colors = ['coral' if c > coeff_var_threshold else 'skyblue' for c in coeff_var_scores] + ax.bar(range(len(coeff_var_scores)), coeff_var_scores, color=colors) + ax.axhline(coeff_var_threshold, color='red', linestyle='--', label=f'Threshold ({coeff_var_threshold}%)') ax.set_title("Coefficient of Variation (Relative Noise)") - ax.set_ylabel("CV %") + ax.set_ylabel("coeff_var %") ax.set_xlabel("Channel Index") ax.legend() @@ -4385,8 +5046,8 @@ def find_bad_channels_cv(raw, cv_threshold=25.0): def hr_calc(raw): - if SHORT_CHANNEL: - short_chans = get_short_channels(raw, max_dist=SHORT_CHANNEL_THRESH) + if SHORT_CHANNELS: + short_chans = get_short_channels(raw, max_dist=SHORT_CHANNELS_THRESHOLD) else: short_chans = None sfreq, signal_trimmed, times_trimmed = short_channel_processing_for_hr(raw, short_chans) @@ -4443,48 +5104,314 @@ def hr_calc(raw): return fig, hr1, hr2, low, high +def trim_participant_data(raw): + if hasattr(raw, 'annotations') and len(raw.annotations) > 0: + # Get time of first event + first_event_time = raw.annotations.onset[0] + trim_time = max(0, first_event_time - SECONDS_TO_KEEP) # Ensure we don't go negative + raw.crop(tmin=trim_time) + # Shift annotation onsets to match new t=0 + + ann = raw.annotations + ann_shifted = Annotations( + onset=ann.onset - trim_time, # shift to start at zero + duration=ann.duration, + description=ann.description + ) + data = raw.get_data() + info = raw.info.copy() + raw = RawArray(data, info) + raw.set_annotations(ann_shifted) + + logger.info(f"Trimmed raw data: start at {trim_time}s (5s before first event), t=0 at new start") + else: + logger.warning("No events found, skipping trim step.") + + fig_trimmed = raw.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="Trimmed Raw", show=False) + return raw, fig_trimmed + + +def remove_bad_channels(raw, bad_channels): + num_bad = len(bad_channels) + + # Check against the threshold + if num_bad > MAX_BAD_CHANNELS: + raise Exception( + f"Data Quality Error: {num_bad} channels flagged for removal, " + f"which exceeds the limit of {MAX_BAD_CHANNELS}. To avoid this, " + f"either lower your filtering parameters or increase MAX_BAD_CHANNELS." + ) + + raw.pick_types(fnirs=True, exclude='bads') + logger.info(f"Physically removed {len(bad_channels)} channels from the dataset.") + return raw + + +def make_and_run_glm(raw_haemo, df_design_matrix): + + glm_est = run_glm(raw_haemo, df_design_matrix, noise_model=NOISE_MODEL, bins=BINS, n_jobs=N_JOBS, verbose=VERBOSITY) + fir_cols = [col for col in df_design_matrix.columns if "_delay_" in col] + + if fir_cols: + # --- FIR MODEL HANDLING (Peak Delay Detection) --- + logger.info("FIR model detected. Dynamically identifying peak delays...") + + # Extract base task conditions (e.g., "Tapping_Left", "Tapping_Right") + base_conditions = list(set(col.split('_delay_')[0] for col in fir_cols)) + + 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()])) + 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) + + # Plot only the peak delays for a clean, single-column topomap per condition + fig_glm_topo = glm_est.plot_topo(conditions=peak_conditions) + + else: + # --- STANDARD HRF MODEL HANDLING --- + # Extract only task conditions (ignore drifts, constants, and short channels) + experimental_conditions = [ + col for col in df_design_matrix.columns + if not any(noise in col.lower() for noise in ['drift', 'constant', 'short']) + ] + fig_glm_topo = glm_est.plot_topo(conditions=experimental_conditions) + + plt.close(fig_glm_topo) + 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 + used for task_cols in generate_contrast_results, reused here so the + plots only ever show things worth looking at.""" + return sorted({ + v for v in values + if not any(ex in str(v).lower() for ex in exclude_list) + }) + + +def generate_channel_results(glm_est, file_path): + df_cha = glm_est.to_dataframe() + df_cha["ID"] = file_path + + df_cha = collapse_fir_condition_column( + df_cha, value_col='theta', + condition_col='Condition', group_cols=['ch_name', 'Chroma', 'ID'] + ) + + return df_cha + + +def generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path): + + rois_formatted = {} + try: + with open(JSON_LOCATION, 'r') as f: + roi_data = json.load(f) + + for region in roi_data.get("regions_of_interest", []): + roi_name = region["name"] + channels = region["channels"] + + mne_channels = [] + for ch in channels: + mne_channels.append(f"{ch} hbo") + mne_channels.append(f"{ch} hbr") + + valid_indices = [] + for ch in mne_channels: + if ch in raw_haemo.ch_names: + idx = raw_haemo.ch_names.index(ch) + valid_indices.append(idx) + + if valid_indices: + rois_formatted[roi_name] = valid_indices + else: + logger.warning(f"No channels from ROI '{roi_name}' found in raw_haemo.") + + except Exception as e: + logger.error(f"Failed to load or parse ROI JSON: {e}") + rois_formatted = {} + + # --------------------------------------------------------------------- + # Tier 2: automatic geometric split — Left/Right, or Front/Back if + # Left/Right isn't usable. Used whenever tier 1 produced nothing, whether + # from a load/parse exception OR a file that loaded fine but matched zero + # channels (the try/except alone doesn't catch that case, since no + # exception is raised). + # --------------------------------------------------------------------- + if not rois_formatted: + logger.error( + f"'{JSON_LOCATION}' produced zero valid ROIs — attempting automatic " + f"geometric fallback (Left/Right, then Front/Back)." + ) + rois_formatted = _build_geometric_fallback_rois(raw_haemo) + + # --------------------------------------------------------------------- + # Tier 3: per-channel (true last resort — only if neither geometric + # split is usable, e.g. missing/degenerate location info). + # --------------------------------------------------------------------- + if not rois_formatted: + logger.warning( + "No usable geometric fallback — falling back to one ROI per " + "channel. Statistics will run at single-channel resolution, not " + "proper ROI aggregation, until regions.json or channel geometry " + "is fixed." + ) + rois_formatted = _build_per_channel_rois(raw_haemo) + + # 3. Calculate ROI results for all conditions + conditions = df_design_matrix.columns + + # Compute output metrics by custom parsed ROIs (now passing lists of integers) + df_roi = glm_est.to_dataframe_region_of_interest(rois_formatted, conditions) + df_roi["ID"] = file_path + + + df_roi = collapse_fir_condition_column( + df_roi, value_col='theta', + condition_col='Condition', group_cols=['ROI', 'Chroma', 'ID'] + ) + + chroma = 'hbo' + value_col = 'theta' + + real_conditions = _real_conditions(df_roi['Condition'].unique(), NUISANCE_EXCLUDE) + sub = df_roi[(df_roi['Chroma'] == chroma) & (df_roi['Condition'].isin(real_conditions))] + + if sub.empty: + print(f"No ROI data for chroma '{chroma}' after excluding nuisance conditions.") + return + + subject_id = sub['ID'].iloc[0] if 'ID' in sub.columns else '' + n_conditions = sub['Condition'].nunique() + + sns.set_theme(style="whitegrid") + fig, ax = plt.subplots(figsize=(max(6, 1.5 * sub['ROI'].nunique()), 5)) + + if n_conditions > 1: + sns.barplot(data=sub, x='ROI', y=value_col, hue='Condition', ax=ax, + edgecolor='black', linewidth=1.2) + else: + sns.barplot(data=sub, x='ROI', y=value_col, ax=ax, + color='#2b5c8f', edgecolor='black', linewidth=1.2) + + ax.axhline(0, color='black', linewidth=1, linestyle='--') + ax.set_ylabel(f'{value_col} ({chroma.upper()})', fontsize=12) + ax.set_xlabel('Region of Interest (ROI)', fontsize=12) + ax.set_title(f"Individual ROI Results ({chroma.upper()})\n{subject_id}", + fontsize=13, fontweight='bold') + plt.tight_layout() + plt.close(fig) + + 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]) + basic_conts = dict( + [(column, contrast_matrix[i]) for i, column in enumerate(df_design_matrix.columns)] + ) + + if HRF_MODEL == "fir": + all_delay_cols = [col for col in df_design_matrix.columns if "_delay_" in col] + all_conditions = sorted({col.split("_delay_")[0] for col in all_delay_cols}) + if not all_conditions: + raise ValueError("No FIR regressors found in the design matrix.") + + contrast_dict = {} + for condition in all_conditions: + delay_cols = [col for col in all_delay_cols if col.startswith(f"{condition}_delay_")] + if not delay_cols: + continue + contrast_vector = np.mean([basic_conts[col] for col in delay_cols], axis=0) + contrast_dict[condition] = contrast_vector + + for cond, contrast_vector in contrast_dict.items(): + contrast = glm_est.compute_contrast(contrast_vector) + df = contrast.to_dataframe() + df["ID"] = file_path + contrast_results_dict[f"{cond}_vs_Zero"] = df + + for cond_a, cond_b in itertools.combinations(all_conditions, 2): + if cond_a not in contrast_dict or cond_b not in contrast_dict: + continue + diff_vector = contrast_dict[cond_a] - contrast_dict[cond_b] + contrast = glm_est.compute_contrast(diff_vector) + df = contrast.to_dataframe() + df["ID"] = file_path + contrast_results_dict[f"{cond_a}_vs_{cond_b}"] = df + + else: + # 0 is NOT a baseline. + # AI Explaination: + # When you regress a single condition (e.g., "Tapping_Right") against zero, the GLM asks: + # "Is the signal during Tapping_Right significantly higher than the average signal across the entire run?" + # Because of systemic physiology (the global blood pressure rise that happens during almost any active task), + # the answer is almost always "Yes, the whole head is higher than the average." + + exclude_list = ["drift", "constant", "short"] + task_cols = [c for c in df_design_matrix.columns if not any(ex in c.lower() for ex in exclude_list)] + + for cond in task_cols: + vec = np.zeros(len(df_design_matrix.columns)) + vec[list(df_design_matrix.columns).index(cond)] = 1 + contrast = glm_est.compute_contrast(vec) + df = contrast.to_dataframe() + df["ID"] = file_path + contrast_results_dict[f"{cond}_vs_Zero"] = df + + for cond_a, cond_b in itertools.combinations(task_cols, 2): + vec = np.zeros(len(df_design_matrix.columns)) + vec[list(df_design_matrix.columns).index(cond_a)] = 1 + vec[list(df_design_matrix.columns).index(cond_b)] = -1 + contrast = glm_est.compute_contrast(vec) + df = contrast.to_dataframe() + df["ID"] = file_path + contrast_results_dict[f"{cond_a}_vs_{cond_b}"] = df + + + return contrast_results_dict def process_participant(file_path, progress_callback=None): + # Step 0: Setting up fig_individual: dict[str, Figure] = {} - logger.info(f"Folding Bypass: {FOLDING_BYP}") - JSON_LOCATION = r"C:\Users\tyler\Desktop\research\flares\regions.json" + config_dict = { + k: globals()[k] + for k in __annotations__ + if k in globals() and k != "REQUIRED_KEYS" + } # Step 1: Preprocessing raw = load_snirf(file_path) fig_raw = raw.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="Loaded Raw", show=False) - fig_individual["Loaded Raw"] = fig_raw + fig_individual["Loaded Raw Data"] = fig_raw if progress_callback: progress_callback(1) logger.info("Step 1 Completed.") # Step 2: Trimming - # TODO: Clean this into a method if TRIM and not FOLDING_BYP: - if hasattr(raw, 'annotations') and len(raw.annotations) > 0: - # Get time of first event - first_event_time = raw.annotations.onset[0] - trim_time = max(0, first_event_time - SECONDS_TO_KEEP) # Ensure we don't go negative - raw.crop(tmin=trim_time) - # Shift annotation onsets to match new t=0 - - ann = raw.annotations - ann_shifted = Annotations( - onset=ann.onset - trim_time, # shift to start at zero - duration=ann.duration, - description=ann.description - ) - data = raw.get_data() - info = raw.info.copy() - raw = RawArray(data, info) - raw.set_annotations(ann_shifted) - - logger.info(f"Trimmed raw data: start at {trim_time}s (5s before first event), t=0 at new start") - else: - logger.warning("No events found, skipping trim step.") - - fig_trimmed = raw.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="Trimmed Raw", show=False) - fig_individual["Trimmed Raw"] = fig_trimmed + raw, fig_trimmed = trim_participant_data(raw) + fig_individual["Trimmed Raw Data"] = fig_trimmed if progress_callback: progress_callback(2) logger.info("Step 2 Completed.") @@ -4496,33 +5423,21 @@ def process_participant(file_path, progress_callback=None): logger.info("Step 3 Completed.") # Step 4: Short/Long Channels - if SHORT_CHANNEL and not FOLDING_BYP: - short_chans = get_short_channels(raw, max_dist=SHORT_CHANNEL_THRESH) - fig_short_chans = short_chans.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="Short Channels Only", show=False) - fig_individual["short"] = fig_short_chans - else: - short_chans = None - raw = get_long_channels(raw, min_dist=0, max_dist=LONG_CHANNEL_THRESH) # keep both short channels and all channels up to the threshold length + 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) + fig_individual["Short Channels Raw Data"] = fig_short_chans + 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 5: Heart Rate if HEART_RATE and not FOLDING_BYP: fig, hr1, hr2, low, high = hr_calc(raw) - fig_individual["PSD"] = fig - fig_individual['HeartRate_PSD'] = hr1 - fig_individual['HeartRate_Time'] = hr2 - - # --- Run it --- - # mark_bads_by_db_threshold(raw, db_limit=-60) - # dead_channels = find_flat_channels(raw) - # print(f"Dead/Flat channels removed: {dead_channels}") - # stuck_channels, movement_scores = find_truly_dead_channels(raw) - # print(f"Stuck Channels: {stuck_channels}") - # late_stage_bads = find_mid_run_flatlines(raw) - # _ = find_flatline_at_end(raw) - # print(f"Channels that died before the end: {late_stage_bads}") - + fig_individual["Power Spectral Density"] = fig + fig_individual['Heart Rate - PSD'] = hr1 + fig_individual['Heart Rate - Time'] = hr2 if progress_callback: progress_callback(5) logger.info("Step 5 Completed.") @@ -4533,8 +5448,8 @@ def process_participant(file_path, progress_callback=None): bad_sci, fig_sci_1, fig_sci_2 = calculate_scalp_coupling(raw, low, high) else: bad_sci, fig_sci_1, fig_sci_2 = calculate_scalp_coupling(raw) - fig_individual["SCI1"] = fig_sci_1 - fig_individual["SCI2"] = fig_sci_2 + fig_individual["Scalp Coupling Index Heatmap"] = fig_sci_1 + fig_individual["Scalp Coupling Index Binary Heatmap"] = fig_sci_2 if progress_callback: progress_callback(6) logger.info("Step 6 Completed.") @@ -4542,7 +5457,7 @@ def process_participant(file_path, progress_callback=None): bad_snr = [] if SNR and not FOLDING_BYP: bad_snr, fig_snr = calculate_signal_noise_ratio(raw) - fig_individual["SNR1"] = fig_snr + fig_individual["Signal To Noise Ratio"] = fig_snr if progress_callback: progress_callback(7) logger.info("Step 7 Completed.") @@ -4550,42 +5465,46 @@ def process_participant(file_path, progress_callback=None): bad_psp = [] if PSP and not FOLDING_BYP: bad_psp, fig_psp1, fig_psp2 = calculate_peak_power(raw) - fig_individual["PSP1"] = fig_psp1 - fig_individual["PSP2"] = fig_psp2 + fig_individual["Peak Spectral Power Heatmap"] = fig_psp1 + fig_individual["Peak Spectral Power Binary Heatmap"] = fig_psp2 if progress_callback: progress_callback(8) logger.info("Step 8 Completed.") - bad_cv = [] - if CV and not FOLDING_BYP: - bad_cv, fig_cv = find_bad_channels_cv(raw, cv_threshold=CV_THRESHOLD) - fig_individual['cv'] = fig_cv + # 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 if progress_callback: progress_callback(9) logger.info("Step 9 Completed.") - bad_range = [] + # Step 10: Median Absolute Deviation + bad_amplitude_range = [] if MAD and not FOLDING_BYP: - bad_range, fig_range = find_bad_channels_range(raw, threshold=MAD_THRESHOLD) - fig_individual['range'] = fig_range + bad_amplitude_range, fig_range = find_bad_channels_by_amplitude_range(raw, threshold=MAD_THRESHOLD) + fig_individual['Median Absolute Deviation'] = fig_range if progress_callback: progress_callback(10) logger.info("Step 10 Completed.") + # Step 11: Power Spectral Density Noise bad_noise = [] if PSD_NOISE and not FOLDING_BYP: - bad_noise, fig_noise = detect_high_freq_noise(raw, db_limit=DB_LIMIT, freq_div=TARGET_FREQ_DIV) - fig_individual['psd_noise'] = fig_noise + bad_noise, fig_noise = detect_spectral_noise_spike(raw, db_limit=DB_LIMIT, freq_div=TARGET_FREQ_DIV) + fig_individual['Power Spectral Density Noise'] = fig_noise if progress_callback: progress_callback(11) logger.info("Step 11 Completed.") + # Step 12: Channel Dropout bad_disp = [] - if CHANNEL_VAR and not FOLDING_BYP: - bad_disp, fig_disp = detect_sensor_displacement(raw, threshold_ratio=CHANNEL_THRESH) - fig_individual['displacement'] = fig_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 if progress_callback: progress_callback(12) logger.info("Step 12 Completed.") - # Step 9: Bad Channels Handling + # 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_cv, bad_range, bad_noise, bad_disp) + 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["fig2"] = fig_dropped fig_individual["fig3"] = fig_raw_before @@ -4595,30 +5514,18 @@ def process_participant(file_path, progress_callback=None): fig_individual["fig4"] = fig_raw_after fig_individual["Compare"] = fig_compare elif BAD_CHANNELS_HANDLING == "Remove": - num_bad = len(bad_channels) - - # Check against the threshold - if num_bad > MAX_BAD_CHANNELS: - raise Exception( - f"Data Quality Error: {num_bad} channels flagged for removal, " - f"which exceeds the limit of {MAX_BAD_CHANNELS}. To avoid this, " - f"either lower your filtering parameters or increase MAX_BAD_CHANNELS." - ) - - raw.pick_types(fnirs=True, exclude='bads') - logger.info(f"Physically removed {len(bad_channels)} channels from the dataset.") - + raw = remove_bad_channels(raw, bad_channels) if progress_callback: progress_callback(13) logger.info("Step 13 Completed.") - # Step 10: Optical Density + # 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 if progress_callback: progress_callback(14) logger.info("Step 14 Completed.") - # Step 11: Temporal Derivative Distribution Repair Filtering + # 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) @@ -4626,37 +5533,21 @@ def process_participant(file_path, progress_callback=None): if progress_callback: progress_callback(15) logger.info("Step 15 Completed.") - # Step 12: Wavelet Filtering + # Step 16: Wavelet Filtering if WAVELET and not FOLDING_BYP: raw_od, fig = calculate_and_apply_wavelet(raw_od) fig_individual["Wavelet"] = fig if progress_callback: progress_callback(16) logger.info("Step 16 Completed.") - # Step 13: Haemoglobin Concentration + # Step 17: Haemoglobin Concentration raw_haemo = beer_lambert_law(raw_od, ppf=calculate_dpf(file_path)) - - # Temporary test - if SHORT_CHANNEL and not FOLDING_BYP and short_chans is not None: - try: - logger.info("Converting raw short channels to hemoglobin concentration...") - # 1. Convert raw short channel intensity to optical density - short_od = optical_density(short_chans) - - # 2. Convert short channel optical density to hemoglobin concentration - # (using the same DPF calculation function as your main data) - short_chans = beer_lambert_law(short_od, ppf=calculate_dpf(file_path)) - - logger.info("Successfully converted short channels to HbO/HbR concentration.") - except Exception as e: - logger.error(f"Failed to convert short channels to hemoglobin: {e}") - 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["BLL"] = fig_raw_haemo_bll if progress_callback: progress_callback(17) logger.info("Step 17 Completed.") - # Step 14: Enhance Negative Correlation + # 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) @@ -4664,7 +5555,7 @@ def process_participant(file_path, progress_callback=None): if progress_callback: progress_callback(18) logger.info("Step 18 Completed.") - # Step 15: Filter + # Step 19: Filter if FILTER and not FOLDING_BYP: raw_haemo, fig_filter, fig_raw_haemo_filter = filter_the_data(raw_haemo) fig_individual["filter1"] = fig_filter @@ -4672,7 +5563,7 @@ def process_participant(file_path, progress_callback=None): if progress_callback: progress_callback(19) logger.info("Step 19 Completed.") - # Step 16: Extracting Events + # Step 20: Extracting Events if not FOLDING_BYP: events, event_dict = events_from_annotations(raw_haemo) fig_events = plot_events(events, event_id=event_dict, sfreq=raw_haemo.info["sfreq"], show=False) @@ -4680,7 +5571,7 @@ def process_participant(file_path, progress_callback=None): if progress_callback: progress_callback(20) logger.info("Step 20 Completed.") - # Step 17: Epoch Calculations + # Step 21: Epoch Calculations if not FOLDING_BYP: epochs, fig_epochs = epochs_calculations(raw_haemo, events, event_dict) for name, fig in fig_epochs: @@ -4688,30 +5579,19 @@ def process_participant(file_path, progress_callback=None): if progress_callback: progress_callback(21) logger.info("Step 21 Completed.") - # Step 18: Design Matrix - df_design_matrix, fig_design_matrix = make_design_matrix(raw_haemo, short_chans) + # Step 22: Design Matrix + raw_haemo, df_design_matrix, fig_design_matrix = make_design_matrix(raw_haemo) # Short channel is re-applied inside this method fig_individual["Design Matrix"] = fig_design_matrix if progress_callback: progress_callback(22) logger.info("Step 22 Completed.") - - # Step 19: Run GLM - glm_est = run_glm(raw_haemo, df_design_matrix, noise_model=NOISE_MODEL, bins=BINS, n_jobs=N_JOBS, verbose=VERBOSITY) - # Not used AppData\Local\Packages\PythonSoftwareFoundation.Python.3.13_qbz5n2kfra8p0\LocalCache\local-packages\Python313\site-packages\nilearn\glm\contrasts.py - # Yes used AppData\Local\Packages\PythonSoftwareFoundation.Python.3.13_qbz5n2kfra8p0\LocalCache\local-packages\Python313\site-packages\mne_nirs\utils\_io.py - - # The p-value is calculated from this t-statistic using the Student's t-distribution with appropriate degrees of freedom. - # p_value = 2 * stats.t.cdf(-abs(t_statistic), df) - # It is a two-tailed p-value. - # It says how likely it is to observe the effect you did (or something more extreme) if the true effect was zero (null hypothesis). - # A small p-value (e.g., < 0.05) suggests the effect is unlikely to be zero — it's "statistically significant." - # A large p-value means the data do not provide strong evidence that the effect is different from zero. - - + # Step 23: Run GLM + glm_est, fig_glm_topo = make_and_run_glm(raw_haemo, df_design_matrix) + fig_individual["GLM Topography"] = fig_glm_topo if progress_callback: progress_callback(23) logger.info("23") - # Step 20: Generate GLM Results + # 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: @@ -4719,210 +5599,31 @@ def process_participant(file_path, progress_callback=None): if progress_callback: progress_callback(24) logger.info("24") - # Step 21: Generate Channel Significance - # if HRF_MODEL == "fir": - # fig_significance = individual_significance(raw_haemo, glm_est) - # for name, fig in fig_significance: - # fig_individual[f"Significance {name}"] = fig + # Step 25: Generate Channel Results + df_cha = generate_channel_results(glm_est, file_path) if progress_callback: progress_callback(25) logger.info("25") - # Step 22: Generate Channel, Region of Interest, and Contrast Results - df_cha = glm_est.to_dataframe() - - # HACK: Comment out line 588 (self._renderer.show()) in _brain.py from MNE - # brain_thing = brain_3d_visualization(cha, raw_haemo) - # brain_individual.append(brain_thing) - # C++ objects made this get rendered on the fly - - - - import json - - try: - with open(JSON_LOCATION, 'r') as f: - roi_data = json.load(f) - - rois_formatted = {} - for region in roi_data.get("regions_of_interest", []): - roi_name = region["name"] - channels = region["channels"] - - # Map JSON names to MNE's expected names - mne_channels = [] - for ch in channels: - mne_channels.append(f"{ch} hbo") - mne_channels.append(f"{ch} hbr") - - # Find the INTEGER index for each valid channel in raw_haemo - valid_indices = [] - for ch in mne_channels: - if ch in raw_haemo.ch_names: - # Get the integer index of this channel name - idx = raw_haemo.ch_names.index(ch) - valid_indices.append(idx) - - if valid_indices: - rois_formatted[roi_name] = valid_indices - else: - logger.warning(f"No channels from ROI '{roi_name}' found in raw_haemo.") - - except Exception as e: - logger.error(f"Failed to load or parse ROI JSON: {e}") - # Fallback to All Channels (as indices) if the JSON fails - rois_formatted = dict(AllChannels=list(range(len(raw_haemo.ch_names)))) - - # 3. Calculate ROI results for all conditions - conditions = df_design_matrix.columns - - # Compute output metrics by custom parsed ROIs (now passing lists of integers) - df_roi = glm_est.to_dataframe_region_of_interest(rois_formatted, conditions) - df_roi["ID"] = file_path - - # rois = dict(AllChannels=range(len(raw_haemo.ch_names))) - # # Calculate ROI for all conditions - # conditions = design_matrix.columns - # # Compute output metrics by ROI - # df_roi = glm_est.to_dataframe_region_of_interest(rois, conditions) - - - # Step 18: Fold channels - # fig_fold_data, fig_fold_legend = fold_channels(raw_haemo) - # fig_individual.append(fig_fold_data) - # fig_individual.append(fig_fold_legend) - print(df_design_matrix) - - - contrast_matrix = np.eye(df_design_matrix.shape[1]) - basic_conts = dict( - [(column, contrast_matrix[i]) for i, column in enumerate(df_design_matrix.columns)] - ) - - if HRF_MODEL == "fir": - all_delay_cols = [col for col in df_design_matrix.columns if "_delay_" in col] - all_conditions = sorted({col.split("_delay_")[0] for col in all_delay_cols}) - - if not all_conditions: - raise ValueError("No FIR regressors found in the design matrix.") - - # Build contrast vectors for each condition - contrast_dict = {} - - for condition in all_conditions: - delay_cols = [ - col for col in all_delay_cols - if col.startswith(f"{condition}_delay_") and - TIME_WINDOW_START <= int(col.split("_delay_")[-1]) <= TIME_WINDOW_END - ] - - if not delay_cols: - continue # skip if no columns found (shouldn't happen?) - - # Average across all delay regressors for this condition - contrast_vector = np.sum([basic_conts[col] for col in delay_cols], axis=0) - contrast_vector /= len(delay_cols) - - contrast_dict[condition] = contrast_vector - + # Step 26: Generate Region of Interest Results + df_roi, fig_roi = generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path) + fig_individual["Region of Interest"] = fig_roi if progress_callback: progress_callback(26) logger.info("26") - # Step 23: Compute Contrast Results - contrast_results_dict = {} - - # 0 is NOT a baseline. - # AI Explaination: - # When you regress a single condition (e.g., "Tapping_Right") against zero, the GLM asks: - # "Is the signal during Tapping_Right significantly higher than the average signal across the entire run?" - # Because of systemic physiology (the global blood pressure rise that happens during almost any active task), - # the answer is almost always "Yes, the whole head is higher than the average." - - if HRF_MODEL == "fir": - - for cond, contrast_vector in contrast_dict.items(): - contrast = glm_est.compute_contrast(contrast_vector) # type: ignore - df = contrast.to_dataframe() - df["ID"] = file_path - contrast_results_dict[cond] = df - - else: - exclude_list = ["drift", "constant", "short"] - task_cols = [c for c in df_design_matrix.columns if not any(ex in c.lower() for ex in exclude_list)] - - # Dictionary to hold all our contrast results - contrast_results_dict = {} - - # 2. Loop through every task to create "Simple" contrasts (Task vs 0) - for cond in task_cols: - # Create a vector of zeros - vec = np.zeros(len(df_design_matrix.columns)) - # Set the index of our condition to 1 - vec[list(df_design_matrix.columns).index(cond)] = 1 - - contrast = glm_est.compute_contrast(vec) - df = contrast.to_dataframe() - df["ID"] = file_path - contrast_results_dict[f"{cond}_vs_Zero"] = df - - # 3. Loop through all combinations to create "Differential" contrasts (Task A vs Task B) - # This uses itertools.combinations to get every unique pair - for cond_a, cond_b in itertools.combinations(task_cols, 2): - vec = np.zeros(len(df_design_matrix.columns)) - vec[list(df_design_matrix.columns).index(cond_a)] = 1 - vec[list(df_design_matrix.columns).index(cond_b)] = -1 - - contrast = glm_est.compute_contrast(vec) - df = contrast.to_dataframe() - df["ID"] = file_path - - # Name the key nicely, e.g., "Tapping_Right_vs_Tapping_Left" - contrast_results_dict[f"{cond_a}_vs_{cond_b}"] = df - - #NOTE: Temporary - export_dir = os.path.join(os.path.dirname(file_path), "exported_contrasts") - os.makedirs(export_dir, exist_ok=True) - base = os.path.splitext(os.path.basename(file_path))[0] - - combined_contrasts = [] - for contrast_name, df in contrast_results_dict.items(): - df = df.copy() - df["contrast_name"] = contrast_name - combined_contrasts.append(df) - - if combined_contrasts: - pd.concat(combined_contrasts, ignore_index=True).to_csv( - os.path.join(export_dir, f"{base}_contrasts.csv"), index=False - ) - - df_roi.to_csv(os.path.join(export_dir, f"{base}_df_ind.csv"), index=False) - df_design_matrix.to_csv(os.path.join(export_dir, f"{base}_design_matrix.csv")) - - logger.info(f"Exported contrast/ROI/design-matrix CSVs to {export_dir}") - df_cha["ID"] = file_path - + # 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 24: Finishing Up + # Step 28: Finishing Up fig_bytes_dict = convert_fig_dict_to_png_bytes(fig_individual) - if FOLDING_BYP: epochs = None sanitize_paths_for_pickle(raw_haemo, epochs) - if progress_callback: progress_callback(28) logger.info("28") - # TODO: Tidy up - # Extract the parameters this file was ran with. No need to return age, gender, group? - config_dict = { - k: globals()[k] - for k in __annotations__ - if k in globals() and k != "REQUIRED_KEYS" - } - - print(config_dict) - + # 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 @@ -5189,7 +5890,7 @@ def get_single_subject_beta_corr(raw_hbo, event_name=None, config=None): else: print("config") if config.get("SHORT_CHANNEL_REGRESSION") == True: - short_chans = get_short_channels(raw_hbo, max_dist=config.get("SHORT_CHANNEL_THRESH")) + short_chans = get_short_channels(raw_hbo, max_dist=config.get("SHORT_CHANNELS_THRESHOLD")) design_matrix = make_first_level_design_matrix( raw=raw_hbo, diff --git a/main.py b/main.py index 76bd61e..912037d 100644 --- a/main.py +++ b/main.py @@ -104,9 +104,9 @@ SECTIONS = [ { "title": "Short/Long Channels", "params": [ - {"name": "SHORT_CHANNEL", "default": True, "type": bool, "help": "This should be set to True if the data has a short channel present in the data. For more information about short channels, please visit the Wiki."}, - {"name": "SHORT_CHANNEL_THRESH", "default": 0.015, "type": float, "depends_on": "SHORT_CHANNEL", "help": "The maximum distance the short channel can be in metres before it is no longer considered a short channel."}, - {"name": "LONG_CHANNEL_THRESH", "default": 0.045, "type": float, "help": "The maximum distance channels can be in metres. Any channel longer than this distance will be discarded."}, + {"name": "SHORT_CHANNELS", "default": True, "type": bool, "help": "This should be set to True if the data has a short channel present in the data. For more information about short channels, please visit the Wiki."}, + {"name": "SHORT_CHANNELS_THRESHOLD", "default": 0.015, "type": float, "depends_on": "SHORT_CHANNEL", "help": "The maximum distance the short channel can be in metres before it is no longer considered a short channel."}, + {"name": "LONG_CHANNELS_THRESHOLD", "default": 0.045, "type": float, "help": "The maximum distance channels can be in metres. Any channel longer than this distance will be discarded."}, ] }, { @@ -144,16 +144,16 @@ SECTIONS = [ ] }, { - "title": "Cross Validation", + "title": "Coefficient of Variation", "params": [ - {"name": "CV", "default": True, "type": bool, "help": "Identifies bad channels using the Coefficient of Variation."}, - {"name": "CV_THRESHOLD", "default": 20, "type": int, "depends_on": "CV", "help": "Noise threshold (%)."}, + {"name": "COEFF_VAR", "default": True, "type": bool, "help": "Identifies bad channels using the Coefficient of Variation."}, + {"name": "COEFF_VAR_THRESHOLD", "default": 20, "type": int, "depends_on": "COEFF_VAR", "help": "Noise threshold (%)."}, ] }, { "title": "Median Absolute Deviation", "params": [ - {"name": "MAD", "default": True, "type": bool, "help": "Identifies bad channels using Mean Absolute Deviation."}, + {"name": "MAD", "default": True, "type": bool, "help": "Identifies bad channels using Median Absolute Deviation."}, {"name": "MAD_THRESHOLD", "default": 4, "type": int, "depends_on": "MAD", "help": "Amount of deviations before the channel is flagged bad."}, ] }, @@ -166,10 +166,10 @@ SECTIONS = [ ] }, { - "title": "Channel Variance", + "title": "Sensor Dropout", "params": [ - {"name": "CHANNEL_VAR", "default": True, "type": bool, "help": "Identifies bad channels based on comparing the variance of the first 25% of the data to the last 25%."}, - {"name": "CHANNEL_THRESH", "default": 0.05, "type": float, "depends_on": "CHANNEL_VAR", "help": "If the end variance is less than this % of the start variance, the channel will be marked as bad."}, + {"name": "SENSOR_DROPOUT", "default": True, "type": bool, "help": "Identifies bad channels based on comparing the variance of the first 25% of the data to the last 25%."}, + {"name": "SENSOR_DROPOUT_VARIANCE_THRESHOLD", "default": 0.05, "type": float, "depends_on": "SENSOR_DROPOUT", "help": "If the end variance is less than this % of the start variance, the channel will be marked as bad."}, ] }, { @@ -206,6 +206,7 @@ SECTIONS = [ "title": "Haemoglobin Concentration", "params": [ # NOTE: Intentionally empty + # TODO: Manual override of PPF? ] }, { @@ -504,6 +505,7 @@ class MainApplication(QMainWindow): self.missing_events_bypass = False self.analysis_clearing_bypass = False self.folding_bypass = False + self.json_location = r"C:\Users\tyler\Desktop\research\flares\regions.json" # Initialization to ensure that saving can occur @@ -1063,7 +1065,8 @@ class MainApplication(QMainWindow): data_map["config_dict"], data_map["fig_bytes_dict"], data_map["contrast_results_dict"], - self.folding_bypass + self.folding_bypass, + self.json_location ] self.launcher_window = ViewerLauncherWidget(*args) @@ -2500,8 +2503,9 @@ def show_critical_error(error_msg): message = ( f"{APP_NAME.upper()} has encountered an unrecoverable error and needs to close.

" f"We are sorry for the inconvenience. An autosave was attempted to be saved to {autosave_path}, but it may not have been saved. " - "If the file was saved, it still may not be intact, openable, or contain the correct data. Use the autosave at your discretion.

" - f"This unrecoverable error was likely due to an error with {APP_NAME.upper()} and not your data.
" + "If the file was saved, it still may not be intact, openable, or contain the correct data. Use the autosave at your own discretion.

" + f"This unrecoverable error was due to an error with {APP_NAME.upper()} and not your data.
" + f"If this crash occured inside a [BETA] branch, it is likely to eventually be fixed.
" f"Please raise an issue here and attach the error file located at {log_path2}

" f"
{error_msg}
" ) diff --git a/src/analysis/crossgroupstats.py b/src/analysis/crossgroupstats.py index 6cff112..8f9bcc3 100644 --- a/src/analysis/crossgroupstats.py +++ b/src/analysis/crossgroupstats.py @@ -9,7 +9,7 @@ License: GPL-3.0 # External library imports import pandas as pd -from flares import run_cross_group_second_level_analysis +from flares import run_cross_group_contrast_analysis, run_cross_group_laterality_analysis, run_cross_group_second_level_analysis from src.shared.flaresbasewidget import CrossGroupUIMixin, FlaresBaseWidget from src.shared.shareddata import APP_NAME @@ -17,39 +17,135 @@ from src.shared.shareddata import APP_NAME PARAMETERIZED_INDEXES = { 0: [ { - "key": "show_optodes", - "label": "Determine what is rendered above the brain. Valid values are 'sensors', 'labels', 'none', 'all'.", - "default": "all", - "type": str, - }, - { - "key": "t_or_theta", - "label": "Specify if t values or theta values should be plotted. Valid values are 't', 'theta'", - "default": "theta", - "type": str, - }, - { - "key": "show_text", - "label": "Display informative text on the top left corner about the contrast.", - "default": "True", - "type": bool, - }, - { - "key": "brain_bounds", - "label": "Graph Upper/Lower Limit", - "default": "1.0", + "key": "p_threshold", + "label": "Significance threshold P-value (e.g. 0.05)", + "default": "0.05", "type": float, }, { - "key": "is_3d", - "label": "Should we display the results in a 3D interactive window?", - "default": "True", + "key": "min_subjects", + "label": "Minimum number of participants to process", + "default": "3", + "type": int, + }, + { + "key": "correction_method", + "label": "Correction method to utilize. Valid values are 'fdr_bh', 'None'", + "default": "fdr_bh", + "type": str, + }, + { + "key": "target_chroma", + "label": "Which chroma to target. Valid values are 'hbo', 'hbr'", + "default": "hbo", + "type": str, + }, + { + "key": "roi_config", + "label": "Location of the ROI config file", + "default": r"C:\Users\tyler\Desktop\research\flares\regions.json", + "type": str, + }, + { + "key": "threshold_topo", + "label": "threshold_topo: TBD", + "default": False, "type": bool, } ], + 1: [ + { + "key": "p_threshold", + "label": "Significance threshold P-value (e.g. 0.05)", + "default": "0.05", + "type": float, + }, + { + "key": "min_subjects", + "label": "Minimum number of participants to process", + "default": "3", + "type": int, + }, + { + "key": "correction_method", + "label": "Correction method to utilize. Valid values are 'fdr_bh', 'None'", + "default": "None", + "type": str, + }, + { + "key": "target_chroma", + "label": "Which chroma to target. Valid values are 'hbo', 'hbr'", + "default": "hbo", + "type": str, + }, + { + "key": "roi_a", + "label": "ROI A (e.g. contralateral region name from regions.json)", + "default": "", + "type": str, + }, + { + "key": "roi_b", + "label": "ROI B (e.g. ipsilateral region name from regions.json)", + "default": "", + "type": str, + } + ], + 2: [ + { + "key": "p_value", + "label": "Significance threshold P-value (e.g. 0.05)", + "default": "0.05", + "type": float, + }, + { + "key": "min_subjects", + "label": "Minimum number of participants to process", + "default": "3", + "type": int, + }, + { + "key": "correction_method", + "label": "Correction method to utilize. Valid values are 'fdr_bh', 'None'", + "default": "fdr_bh", + "type": str, + }, + { + "key": "target_chroma", + "label": "Which chroma to target. Valid values are 'hbo', 'hbr'", + "default": "hbo", + "type": str, + }, + { + "key": "roi_config", + "label": "Location of the ROI config file", + "default": r"C:\Users\tyler\Desktop\research\flares\regions.json", + "type": str, + }, + { + "key": "contrast_name", + "label": "Name of the contrast to use", + "default": "", + "type": str, + }, + ], } +DESCRIPTION = """0. Raw ROI Comparison (run_cross_group_second_level_analysis) +\nCompares one ROI's raw response magnitude between two independent groups (e.g. control vs. target) for a given condition, using Welch's t-test. A significant result means the two populations differ in this ROI's response magnitude for this condition. It does not tell you whether that difference is a real, localized, task-specific effect or a generic between-population difference - different overall vascular reactivity, arousal, or skull/scalp optical properties can produce the exact same statistical signature, and two independently recruited groups (especially patients vs. healthy controls) are considerably more likely to differ this way than two subsets of one study population. +\nIf you expected a group difference and didn't find one, the most common cause is within-group heterogeneity swallowing a real between-group difference - a "target" population (e.g. a clinical group) is often more variable than a tightly-screened control group, and that added within-group variance directly weakens a between-group t-test even if the group means truly differ. Small per-group sample sizes compound this. It's also possible the true difference between your groups isn't in raw magnitude at all, but in spatial specificity or task-differentiation - which is exactly why the laterality and contrast-comparison methods exist alongside this one; a null result here doesn't rule those out. +\n\n1. Laterality Comparison (run_cross_group_laterality_analysis) +\nComputes each subject's own contralateral-minus-ipsilateral laterality index first, then compares those indices between the two groups with Welch's t-test. A significant result means the degree of spatial specificity/lateralization differs between the two populations - a claim about lateralization itself, harder to explain away as a generic population confound since person-level differences in overall reactivity largely cancel before the group comparison happens. It says nothing about overall response magnitude between groups (a group could have identical laterality but very different raw amplitude), and it only uses subjects who have both the contra and ipsi ROI valid, so it can lose subjects the raw-ROI comparison would have kept. +\nNon-significance here has two likely sources, and it's worth distinguishing them. First, the same covariance issue from the within-group paired test applies across a whole group: if contra/ipsi responses aren't well-correlated within subjects, the laterality index itself is noisier than either ROI alone, and that added noise now has to clear a between-group test on top of it - a double power cost at small N. Second, and more informative if true: the groups may genuinely have similar lateralization but differ in overall magnitude instead, in which case this test correctly returns null while method 4 (raw comparison) should be the one to look at. +\n\n2. Contrast Comparison (run_cross_group_contrast_analysis) +\nCompares a jointly-fit task contrast (e.g. Task A minus Task B, estimated together within each subject's GLM), aggregated to ROI level, between two independent groups. A significant result means one group differentiates between the two tasks more or less than the other does, at this specific ROI - with systemic noise cancelled at the model-fitting stage, the same benefit that makes the within-group version of this method the strongest of that trio. As with the within-group version, it does not by itself say where a difference is localized unless you compare sign/pattern across multiple ROIs - opposite-signed group differences across regions point to something spatially specific, same-signed differences everywhere point to a diffuse/non-specific group difference (e.g. one group simply has stronger contrast responses across the whole head). +\nIf this comes back non-significant despite an expected group difference, check first whether the underlying single-subject contrast estimates are noisy for either group - small per-group N means the joint contrast's precision depends on the same limited subject count as everything else, and a noisy input propagates all the way through the ROI aggregation. It's also possible for a real, localized sub-regional effect to get washed out by ROI averaging itself: if only part of an ROI's channels actually show the group difference while others don't, the inverse-variance-weighted average can dilute it toward null - in that case, a finer-grained ROI definition (splitting the region further) may recover the effect that a coarser ROI averaged away. Finally, FDR correction across every ROI tested reduces power exactly as it does everywhere else in this framework - a real but modest effect can fail to survive correction even when the raw p-value would have looked convincing on its own. +\n\n +\nWhy channels needed to be aggregated into ROIs: Testing every channel independently means paying a steep multiple-comparisons tax - with dozens of channels, FDR/Bonferroni correction demands very large effect sizes to call anything significant, and at small subject counts (n=5) essentially nothing survives even when a real, consistent effect exists. Collapsing channels into a handful of anatomically meaningful ROIs cuts the number of independent tests from a minimum of ~40 down to 2-8, which lets a genuinely present effect actually clear correction. It also matches the scientific question better: you have a hypothesis about regions (contralateral motor cortex, prefrontal cortex), not about individual source-detector pairs, so testing at the ROI level is testing the thing you actually believe in, using inverse-variance weighting so noisier channels contribute less to the region's combined estimate rather than diluting it equally. +\nWhy some analyses needed contrasts instead of raw values: A single condition's GLM beta is only ever measured relative to the model's implicit intercept, and that intercept absorbs whatever's happening for the rest of the recording - including systemic physiology (blood pressure, arousal, general vascular reactivity) that rises during almost any active task, not just the one you care about. Testing a raw "vs. zero" value can't tell a real, localized neural response apart from that shared full-head noise. A contrast - either a within-subject spatial subtraction (ROI A minus ROI B) or a jointly-fit task contrast (Condition A minus Condition B, estimated together in one GLM) cancels out whatever's common to both halves of the subtraction, leaving something closer to the actual differential signal. +\nWhy a minimum subject count is enforced: Every one of these tests is a t-test, and a t-test's ability to detect a real effect (its power) depends heavily on degrees of freedom - at n=5 (df=4), even a fairly large true effect can produce a middling p-value, and at n=2 (df=1) the test is barely meaningful at all regardless of the underlying data. The min_subjects floor exists to stop a channel or ROI from being silently tested (and potentially reported as significant or non-significant) on a sample too small for the resulting p-value to mean anything reliable - it's better to explicitly skip and flag an underpowered channel than to quietly produce a number that looks statistically legitimate but isn't backed by enough independent observations to trust.""" + class CrossGroupStatsWidget(CrossGroupUIMixin, FlaresBaseWidget): def __init__(self, haemo_dict, cha_dict, df_ind_dict, design_matrix_dict, contrast_results_dict, group_dict): @@ -62,7 +158,7 @@ class CrossGroupStatsWidget(CrossGroupUIMixin, FlaresBaseWidget): self.contrast_results_dict = contrast_results_dict self.group_dict = group_dict - self.setup_cross_group_ui(["0 (Compute Statistics)"]) + self.setup_cross_group_ui(["0 (Raw ROI Comparison)", "1 (Laterality Comparison)", "2 (Contrast Comparison)",], placeholder_text=DESCRIPTION) def process_request(self): @@ -94,6 +190,14 @@ class CrossGroupStatsWidget(CrossGroupUIMixin, FlaresBaseWidget): # Visualizations for idx in selected_indexes: if idx == 0: + params = param_values.get(idx, {}) + p_threshold = params.get("p_threshold", 0.05) + min_subjects = params.get("min_subjects", 3) + correction_method = params.get("correction_method", "fdr_bh") + target_chroma = params.get("target_chroma", "hbo") + roi_config = params.get("roi_config", r"C:\Users\tyler\Desktop\research\flares\regions.json") + threshold_topo = params.get("threshold_topo", False) + run_cross_group_second_level_analysis( df_roi_all=df_ind_combined, # Individual stats dataframe file_paths_a=file_paths_a, @@ -102,13 +206,119 @@ class CrossGroupStatsWidget(CrossGroupUIMixin, FlaresBaseWidget): group_b_name=self.group_b_dropdown.currentText(), df_cha_all=cha_combined, raw_haemo=p_haemo, - p_threshold=0.05, - min_subjects=3, - correction_method='fdr_bh', - target_chroma='hbo', + p_threshold=p_threshold, + min_subjects=min_subjects, + correction_method=correction_method, + target_chroma=target_chroma, selected_event=selected_event, - roi_config=r"C:\Users\tyler\Desktop\research\flares\regions.json", - threshold_topo=False # Shows the raw difference map (Unthresholded) + roi_config=roi_config, + threshold_topo=threshold_topo # Shows the raw difference map (Unthresholded) ) + elif idx == 1: + if not selected_event: + print("Laterality comparison requires a specific event/condition " + "to be selected first.") + continue + + params = param_values.get(idx, {}) + p_threshold = params.get("p_threshold", 0.05) + min_subjects = params.get("min_subjects", 3) + correction_method = params.get("correction_method", "None") + target_chroma = params.get("target_chroma", "hbo") + roi_a = params.get("roi_a", "").strip() + roi_b = params.get("roi_b", "").strip() + + if not roi_a or not roi_b: + print("Both a contralateral and ipsilateral ROI name must be specified.") + continue + + if correction_method == "None": + correction_method = None + + # Build each group's dataframe directly from the dict using + # the file-path lists as keys - no ID cleaning/matching needed. + def _build_group_df(file_paths, dict_source): + valid_dfs = [ + dict_source[fp] for fp in file_paths + if fp in dict_source and isinstance(dict_source[fp], pd.DataFrame) + and not dict_source[fp].empty + ] + return pd.concat(valid_dfs, ignore_index=True) if valid_dfs else pd.DataFrame() + + df_roi_a = _build_group_df(file_paths_a, self.df_ind_dict) + df_roi_b = _build_group_df(file_paths_b, self.df_ind_dict) + + if df_roi_a.empty or df_roi_b.empty: + print("No ROI data (df_ind) found for one or both groups.") + continue + + run_cross_group_laterality_analysis( + df_roi_all_a=df_roi_a, + df_roi_all_b=df_roi_b, + roi_pairs=(roi_a, roi_b), + condition=selected_event, + group_a_name=self.group_a_dropdown.currentText(), + group_b_name=self.group_b_dropdown.currentText(), + target_chroma=target_chroma, + min_subjects=min_subjects, + p_threshold=p_threshold, + correction_method=correction_method, + roi_contra_label=roi_a, + roi_ipsi_label=roi_b, + ) + + elif idx == 2: + params = param_values.get(idx, {}) + p_threshold = params.get("p_threshold", 0.05) + min_subjects = params.get("min_subjects", 3) + correction_method = params.get("correction_method", "fdr_bh") + target_chroma = params.get("target_chroma", "hbo") + roi_config = params.get("roi_config", r"C:\Users\tyler\Desktop\research\flares\regions.json") + contrast_name = params.get("contrast_name", "") + + if not contrast_name: + print("A contrast name must be specified (e.g. '2.0_vs_3.0').") + continue + + # Build each group's channel-level contrast dataframe + # directly from contrast_results_dict, keyed by file path - + # same dict-key approach as the laterality patch, avoids + # any ID-string matching. + def _build_group_contrast_df(file_paths, contrast_dict, name): + all_rows = [] + for fp in file_paths: + condition_dfs = contrast_dict.get(fp) + if condition_dfs is None: + print(f" [MISSING] '{fp}' not found in contrast_results.") + continue + if name in condition_dfs: + df = condition_dfs[name].copy() + df["ID"] = fp + df["contrast_name"] = name + all_rows.append(df) + else: + print(f" [MISSING CONTRAST] '{name}' not available for '{fp}'.") + return pd.concat(all_rows, ignore_index=True) if all_rows else pd.DataFrame() + + df_contrasts_a = _build_group_contrast_df(file_paths_a, self.contrast_results_dict, contrast_name) + df_contrasts_b = _build_group_contrast_df(file_paths_b, self.contrast_results_dict, contrast_name) + + if df_contrasts_a.empty or df_contrasts_b.empty: + print("No contrast data found for one or both groups.") + continue + + run_cross_group_contrast_analysis( + df_contrasts_a=df_contrasts_a, + df_contrasts_b=df_contrasts_b, + contrast_name=contrast_name, + roi_json_path=roi_config, + group_a_name=self.group_a_dropdown.currentText(), + group_b_name=self.group_b_dropdown.currentText(), + target_chroma=target_chroma, + min_subjects=min_subjects, + p_threshold=p_threshold, + correction_method=correction_method, + ) + else: print("no") \ No newline at end of file diff --git a/src/analysis/intergroupstats.py b/src/analysis/intergroupstats.py index e486b3c..af14331 100644 --- a/src/analysis/intergroupstats.py +++ b/src/analysis/intergroupstats.py @@ -17,63 +17,111 @@ from src.shared.shareddata import APP_NAME PARAMETERIZED_INDEXES = { 0: [ { - "key": "info", - "label": "Tests whether one ROI's response during one condition reliably differs from zero across subjects.\nIf significant, you can claim: This region's signal during this condition is consistently non-zero across your sample - not just noise.\nIt does NOT say: Whether that response is localized/specific to this region, or whether it reflects real neural activity versus systemic physiology (blood pressure, arousal) shared across the whole head during any active task.", - "default": "Okay.", - "type": str, - }, - { - "key": "p_value", + "key": "p_threshold", "label": "Significance threshold P-value (e.g. 0.05)", "default": "0.05", "type": float, }, + { + "key": "min_subjects", + "label": "Minimum number of participants to process", + "default": "5", + "type": int, + }, + { + "key": "correction_method", + "label": "Correction method to utilize. Valid values are 'fdr_bh', 'None'", + "default": "fdr_bh", + "type": str, + }, + { + "key": "target_chroma", + "label": "Which chroma to target. Valid values are 'hbo', 'hbr'", + "default": "hbo", + "type": str, + }, { "key": "graph_bounds", - "label": "Graph Y-Limit (Optional, e.g. 1e-5)", - "default": "0.0", # Set to 0.0 to auto-scale + "label": "Graph Upper/Lower Limit", + "default": "0.0", "type": float, } ], 1: [ { - "key": "info", - "label": "For one condition, subtracts each subject's ROI_A response from their ROI_B response, then tests whether that per-subject difference is reliably non-zero.\nIf significant, you can claim: The two regions respond differently from each other during this specific condition - a real spatial contrast, since shared systemic noise partially cancels in the subtraction.\nIt does NOT say: Anything about whether the condition itself produced meaningful activity at all (only a relative difference between two places); and its power depends on the two ROIs' noise being correlated across subjects, which isn't guaranteed.", - "default": "Okay.", + "key": "p_threshold", + "label": "Significance threshold P-value (e.g. 0.05)", + "default": "0.05", + "type": float, + }, + { + "key": "min_subjects", + "label": "Minimum number of participants to process", + "default": "5", + "type": int, + }, + { + "key": "correction_method", + "label": "Correction method to utilize. Valid values are 'fdr_bh', 'None'", + "default": "None", + "type": str, + }, + { + "key": "target_chroma", + "label": "Which chroma to target. Valid values are 'hbo', 'hbr'", + "default": "hbo", "type": str, }, { "key": "roi_a", "label": "ROI A (e.g. contralateral region name from regions.json)", - "default": "", - "type": str, + "default": [], + "type": list, }, { "key": "roi_b", "label": "ROI B (e.g. ipsilateral region name from regions.json)", - "default": "", - "type": str, - }, - { - "key": "p_value", - "label": "Significance threshold P-value (e.g. 0.05)", - "default": "0.05", - "type": float, - }, + "default": [], + "type": list, + } ], 2: [ - { - "key": "info", - "label": "Uses a contrast fit jointly within each subject's GLM (Condition A minus Condition B, estimated together), then aggregates that per-channel contrast to ROI level and tests it against zero across subjects.\nIf significant, you can claim: The two conditions produce reliably different responses at this ROI, with systemic noise largely cancelled at the model-fitting stage itself - the most statistically efficient of the three.\nIt does NOT say: Which region the difference is localized to, unless you compare the sign/pattern across multiple ROIs", - "default": "Okay.", - "type": str, - }, { "key": "p_value", "label": "Significance threshold P-value (e.g. 0.05)", "default": "0.05", "type": float, }, + { + "key": "min_subjects", + "label": "Minimum number of participants to process", + "default": "5", + "type": int, + }, + { + "key": "correction_method", + "label": "Correction method to utilize. Valid values are 'fdr_bh', 'None'", + "default": "fdr_bh", + "type": str, + }, + { + "key": "target_chroma", + "label": "Which chroma to target. Valid values are 'hbo', 'hbr'", + "default": "hbo", + "type": str, + }, + { + "key": "contrast_name", + "label": "Name of the contrast to use", + "default": [], + "type": list, + }, + { + "key": "weighted", + "label": "Use inverse-variance weighting to minimize noisy channels", + "default": True, + "type": bool, + }, { "key": "graph_bounds", "label": "Graph Upper/Lower Limit", @@ -84,9 +132,23 @@ PARAMETERIZED_INDEXES = { } +DESCRIPTION = """0. ROI vs. Zero (run_roi_second_level_analysis) +\nTests whether one ROI's response during one condition reliably differs from zero across subjects - a one-sample t-test on each subject's ROI-averaged theta. A significant result means the region's signal during this condition is consistently non-zero across your sample, not just noise. It does not tell you whether the response is localized/specific to this region, or whether it reflects real neural activity versus systemic physiology (blood pressure, arousal) shared across the whole head during any active task - a single-condition "vs. zero" test can't distinguish those two explanations on its own. +\nIf you expected significance here and didn't get it, likely causes include: the sample size is simply small relative to between-subject variability in true response magnitude or HRF shape (individual differences in timing/amplitude inflate the variance a t-test divides by); the ROI's channel composition differs slightly across subjects (missing channels get down-weighted or excluded from the inverse-variance average, diluting a real signal); FDR correction across many ROIs is suppressing a modest true effect that would clear an uncorrected threshold; or the condition itself may not reliably engage this region the way you assumed (worth checking the single-subject/individual-level results for this ROI before concluding the group effect isn't there). +\n\n1. Paired ROI Contrast (run_roi_paired_contrast_analysis) +\nFor one condition, subtracts each subject's ROI_A response from their ROI_B response, then tests whether that per-subject difference is reliably non-zero. A significant result is a genuine spatial contrast - the two regions respond differently from each other during this specific condition, with shared systemic noise partially cancelling in the subtraction. It says nothing about whether the condition produced meaningful activity at all (only a relative difference between two places), and its power depends entirely on ROI_A and ROI_B varying together across subjects - an assumption that isn't guaranteed. +\nIf this test underperforms a plain ROI-vs-zero result, which can occur, the most likely explanation is that ROI_A and ROI_B's noise isn't well-correlated across your subjects. The math is variance(A - B) = variance(A) + variance(B) - 2·covariance(A,B): subtraction only helps when the shared/systemic component is large relative to independent noise in each region. If the two regions are picking up largely independent noise sources (motion artifact affecting one side more, different channel quality, etc.), subtracting adds variance rather than removing it, and can turn a detectable single-ROI effect into an underpowered paired one. Small sample size makes this worse, since the covariance itself is poorly estimated with few subjects. +\n\n2. Joint Contrast, ROI-Aggregated (aggregate_channel_contrasts_to_roi + one-sample test) +\nUses a contrast fit jointly within each subject's GLM (Condition A minus Condition B, estimated together), then aggregates that per-channel contrast to ROI level using inverse-variance weighting, and tests it against zero across subjects. A significant result means the two conditions produce reliably different responses at this ROI, with systemic noise largely cancelled at the model-fitting stage itself - the most statistically efficient of the three within-group methods, since the correlation between conditions is handled natively rather than inferred afterward. It does not tell you where the difference is localized on its own - for that, compare the sign/pattern across multiple ROIs: opposite signs across regions indicates a real, spatially-specific effect, while the same sign everywhere suggests diffuse/systemic noise rather than localized activity (as seen when comparing a real task-vs-task contrast against a task-vs-inert-marker contrast). +\nIf this comes back non-significant despite expecting an effect, first check whether the two conditions are actually similar enough in their neural engagement of this ROI that a small or genuinely near-zero contrast is the correct answer - not every ROI should differentiate every pair of tasks, and a null result here can be the right result. Beyond that: FDR correction across every ROI in your regions file can suppress a real but modest contrast; the inverse-variance weighting can be destabilized if a few channels within the ROI have very noisy or near-zero t-statistics (their standard error estimate becomes huge or unstable); and - as always - small subject counts limit the achievable degrees of freedom regardless of how clean the underlying per-channel estimates are. +\n\n +\nWhy channels needed to be aggregated into ROIs: Testing every channel independently means paying a steep multiple-comparisons tax - with dozens of channels, FDR/Bonferroni correction demands very large effect sizes to call anything significant, and at small subject counts (n=5) essentially nothing survives even when a real, consistent effect exists. Collapsing channels into a handful of anatomically meaningful ROIs cuts the number of independent tests from a minimum of ~40 down to 2-8, which lets a genuinely present effect actually clear correction. It also matches the scientific question better: you have a hypothesis about regions (contralateral motor cortex, prefrontal cortex), not about individual source-detector pairs, so testing at the ROI level is testing the thing you actually believe in, using inverse-variance weighting so noisier channels contribute less to the region's combined estimate rather than diluting it equally. +\nWhy some analyses needed contrasts instead of raw values: A single condition's GLM beta is only ever measured relative to the model's implicit intercept, and that intercept absorbs whatever's happening for the rest of the recording - including systemic physiology (blood pressure, arousal, general vascular reactivity) that rises during almost any active task, not just the one you care about. Testing a raw "vs. zero" value can't tell a real, localized neural response apart from that shared full-head noise. A contrast - either a within-subject spatial subtraction (ROI A minus ROI B) or a jointly-fit task contrast (Condition A minus Condition B, estimated together in one GLM) cancels out whatever's common to both halves of the subtraction, leaving something closer to the actual differential signal. +\nWhy a minimum subject count is enforced: Every one of these tests is a t-test, and a t-test's ability to detect a real effect (its power) depends heavily on degrees of freedom - at n=5 (df=4), even a fairly large true effect can produce a middling p-value, and at n=2 (df=1) the test is barely meaningful at all regardless of the underlying data. The min_subjects floor exists to stop a channel or ROI from being silently tested (and potentially reported as significant or non-significant) on a sample too small for the resulting p-value to mean anything reliable - it's better to explicitly skip and flag an underpowered channel than to quietly produce a number that looks statistically legitimate but isn't backed by enough independent observations to trust.""" + class InterGroupStatsWidget(InterGroupUIMixin, FlaresBaseWidget): - def __init__(self, haemo_dict, cha, df_ind, design_matrix, contrast_results, group): + def __init__(self, haemo_dict, cha, df_ind, design_matrix, contrast_results, group, json_location): super().__init__("InterGroupStats") self.setWindowTitle(f"Inter-Group Stats Viewer - {APP_NAME.upper()}") self.haemo_dict = haemo_dict @@ -95,12 +157,12 @@ class InterGroupStatsWidget(InterGroupUIMixin, FlaresBaseWidget): self.design_matrix = design_matrix self.contrast_results = contrast_results self.group = group + self.json_location = json_location - self.setup_inter_group_ui(["0 (Significance)", "1 (More significasd)", "2 (moreeee)"]) - + self.setup_inter_group_ui(["0 (ROI vs. Zero)", "1 (Paired ROI Contrast)", "2 (Joint Contrast, ROI-Aggregated)"], placeholder_text=DESCRIPTION) def process_request(self): - request = self.get_common_request_data(PARAMETERIZED_INDEXES) + request = self.get_common_request_data(PARAMETERIZED_INDEXES, self.json_location, self.contrast_results) if request is None: return @@ -137,9 +199,15 @@ class InterGroupStatsWidget(InterGroupUIMixin, FlaresBaseWidget): for idx in selected_indexes: if idx == 0: params = param_values.get(idx, {}) - p_val = params.get("p_value", 0.05) + p_threshold = params.get("p_threshold", 0.05) + min_subjects = params.get("min_subjects", 5) + correction_method = params.get("correction_method", "fdr_bh") + target_chroma = params.get("target_chroma", "hbo") graph_bounds = params.get("graph_bounds", 0.0) + if correction_method == "None": + correction_method = None + if df_group.empty: print("No ROI data (df_ind) found for selected participants.") continue @@ -165,90 +233,77 @@ class InterGroupStatsWidget(InterGroupUIMixin, FlaresBaseWidget): else: all_cha_filtered = all_cha - # --------------------------------------------------------------------- - # run_roi_second_level_analysis - # - # Tests: is this ROI's activation reliably different from zero, for one - # condition, across subjects? (One-sample t-test per ROI.) - # - # CAUTION: "vs zero" includes systemic/global physiology shared across - # the whole head (blood pressure, arousal, etc.), not just localized - # neural response — a significant result here doesn't by itself prove - # the effect is spatially specific to this ROI. - # --------------------------------------------------------------------- run_roi_second_level_analysis( df_roi_all=df_filtered, df_cha_all=all_cha_filtered, raw_haemo=p_haemo, - p_threshold=p_val, - min_subjects=len(selected_file_paths), - correction_method='fdr_bh', - target_chroma='hbo', + p_threshold=p_threshold, + min_subjects=min_subjects, + correction_method=correction_method, + target_chroma=target_chroma, graph_bounds=graph_bounds if graph_bounds > 0.0 else None, - roi_config=r"C:\Users\tyler\Desktop\research\flares\regions.json" + roi_config=self.json_location ) elif idx == 1: + params = param_values.get(idx, {}) + p_threshold = params.get("p_threshold", 0.05) + min_subjects = params.get("min_subjects", 5) + correction_method = params.get("correction_method", "None") + target_chroma = params.get("target_chroma", "hbo") + roi_a = params.get("roi_a", "").strip() + roi_b = params.get("roi_b", "").strip() + if not selected_event: print("Paired ROI contrast requires a specific event/condition " - "to be selected — pick one from the Event dropdown first.") + "to be selected - pick one from the Event dropdown first.") continue if df_group.empty: print("No ROI data (df_ind) found for selected participants.") continue - params = param_values.get(idx, {}) - roi_a = params.get("roi_a", "").strip() - roi_b = params.get("roi_b", "").strip() - p_val = params.get("p_value", 0.05) + if correction_method == "None": + correction_method = None if not roi_a or not roi_b: print("Both ROI A and ROI B must be specified.") continue - # --------------------------------------------------------------------- - # run_roi_paired_contrast_analysis - # - # Tests: within one condition, does ROI_A's activation differ from - # ROI_B's, per subject? (Paired one-sample t-test on the per-subject - # difference, e.g. Right_PFC minus Left_PFC for a laterality check.) - # - # Only gains power over testing ROI_A and ROI_B separately if the two - # ROIs' noise is correlated across subjects (shared systemic component - # cancels in the subtraction). If they vary independently, this test - # can be WEAKER than testing either ROI alone — check per-subject - # correlation between ROI_A and ROI_B if this test underperforms. run_roi_paired_contrast_analysis( - df_roi_all=df_group, # unfiltered — function filters internally + df_roi_all=df_group, roi_pairs=(roi_a, roi_b), condition=selected_event, - target_chroma='hbo', - min_subjects=min(5, len(selected_file_paths)), - p_threshold=p_val, - correction_method=None, # single pre-specified contrast + target_chroma=target_chroma, + min_subjects=min_subjects, + p_threshold=p_threshold, + correction_method=correction_method, roi_a_label=roi_a, roi_b_label=roi_b, ) elif idx == 2: + params = param_values.get(idx, {}) + p_threshold = params.get("p_threshold", 0.05) + min_subjects = params.get("min_subjects", 5) + correction_method = params.get("correction_method", "fdr_bh") + target_chroma = params.get("target_chroma", "hbo") + contrast_name = params.get("contrast_name", "2.0_vs_3.0") + weighted = params.get("weighted", True) + graph_bounds = params.get("graph_bounds", 0.0) + if not selected_event: print("Joint contrast ROI analysis requires a specific contrast " "to be selected from the Event dropdown first.") continue - - # Build the channel-level contrast dataframe for selected - # participants + selected contrast, same pattern used in - # GroupViewerWidget.show_brain_images. - contrast_name = "15.0_vs_2.0" # <-- change this to test other contrasts - print(f"[TEMP HARDCODE] Using contrast '{contrast_name}' " - f"instead of dropdown selection ('{selected_event}') for option 2.") - - # Build the channel-level contrast dataframe for selected - # participants + selected contrast, same pattern used in - # GroupViewerWidget.show_brain_images. + + if not contrast_name: + print("Contrast name must be specified.") + continue + + all_contrasts = [] for fp in selected_file_paths: condition_dfs = self.contrast_results.get(fp) @@ -258,10 +313,6 @@ class InterGroupStatsWidget(InterGroupUIMixin, FlaresBaseWidget): if contrast_name in condition_dfs: df = condition_dfs[contrast_name].copy() df["ID"] = fp - # contrast_results dict values don't carry a - # contrast_name column themselves — that's only - # stamped on during CSV export. Add it here since - # aggregate_channel_contrasts_to_roi requires it. df["contrast_name"] = contrast_name all_contrasts.append(df) else: @@ -275,16 +326,13 @@ class InterGroupStatsWidget(InterGroupUIMixin, FlaresBaseWidget): df_contrasts = pd.concat(all_contrasts, ignore_index=True) - params = param_values.get(idx, {}) - p_val = params.get("p_value", 0.05) - graph_bounds = params.get("graph_bounds", 0.0) - try: roi_theta = aggregate_channel_contrasts_to_roi( df_contrasts, - roi_json_path=r"C:\Users\tyler\Desktop\research\flares\regions.json", - weighted=True, + roi_json_path=self.json_location, + weighted=weighted, ) + except Exception as e: print(f"Failed to aggregate contrasts to ROI: {e}") continue @@ -294,22 +342,22 @@ class InterGroupStatsWidget(InterGroupUIMixin, FlaresBaseWidget): "(check regions.json channel names against this montage).") continue + # TODO: Come back to this # df_cha_all intentionally omitted (None): the topography # section of run_roi_second_level_analysis expects # single-condition Condition values in df_cha_all, which - # doesn't semantically match a contrast name — skip it here + # doesn't semantically match a contrast name - skip it here # rather than pass mismatched data. run_roi_second_level_analysis( df_roi_all=roi_theta, df_cha_all=None, raw_haemo=p_haemo, - p_threshold=p_val, - min_subjects=min(5, len(selected_file_paths)), - correction_method='fdr_bh', - target_chroma='hbo', + p_threshold=p_threshold, + min_subjects=min_subjects, + correction_method=correction_method, + target_chroma=target_chroma, graph_bounds=graph_bounds if graph_bounds > 0.0 else None, ) - else: print(f"No method defined for index {idx}") \ No newline at end of file diff --git a/src/analysis/participantfoldchannels.py b/src/analysis/participantfoldchannels.py index 4f96d27..e768d52 100644 --- a/src/analysis/participantfoldchannels.py +++ b/src/analysis/participantfoldchannels.py @@ -69,9 +69,9 @@ def single_participant_worker(file_path, raw_data, result_queue, progress_queue) """ Runs inside its own dedicated process """ p_name = os.path.basename(file_path) try: - import flares as flares + from flares import fold_channels # Perform the heavy fold_channels logic - channel_results = flares.fold_channels(raw_data, p_name, progress_queue) + channel_results = fold_channels(raw_data, p_name, progress_queue) # Hand back results and signal completion result_queue.put({file_path: channel_results}) @@ -736,13 +736,13 @@ class ProcessOrchestrator(QObject): def run(self): try: - # 🟢 [Delay 1 Fix] Instantiate Manager completely off the main thread + # Instantiate Manager completely off the main thread manager = Manager() result_queue = manager.Queue() progress_queue = manager.Queue() active_processes = [] - # 🟢 [Delay 2 Fix] Perform heavy pickling loop safely in the background + # Perform heavy pickling loop safely in the background for file_path in self.selected_files: p = Process( target=self.worker_func, @@ -896,7 +896,7 @@ class ParticipantFoldChannelsWidget(FlaresBaseWidget): self.progress_queue = progress_queue self.active_processes = active_processes - # 🟢 Safely initialize and trigger your polling listener + # Safely initialize and trigger the polling listener self.completed_count = 0 self.result_timer = QTimer() self.result_timer.timeout.connect(self.check_parallel_results) diff --git a/src/shared/flaresbasewidget.py b/src/shared/flaresbasewidget.py index a752a7a..84e4bba 100644 --- a/src/shared/flaresbasewidget.py +++ b/src/shared/flaresbasewidget.py @@ -7,6 +7,7 @@ License: GPL-3.0 """ import os +import json from PySide6.QtWidgets import QApplication, QComboBox, QDialog, QGridLayout, QHBoxLayout, QLabel, QLineEdit, QListView, QMessageBox, QPushButton, QScrollArea, QVBoxLayout, QWidget, QFrame, QSpinBox from PySide6.QtGui import QStandardItemModel, QStandardItem, QPixmap, QIntValidator, QDoubleValidator @@ -84,28 +85,53 @@ class ParameterInputDialog(QDialog): self.params_dict = params_dict self.inputs = {} # {(idx, param_key): QLineEdit} - layout = QVBoxLayout(self) + main_layout = QVBoxLayout(self) intro_label = QLabel( "Some methods require parameters to continue:\n" "Clicking OK will simply use default values if input is left empty." ) - layout.addWidget(intro_label) + main_layout.addWidget(intro_label) + self.setMinimumWidth(400) + self.scroll = QScrollArea() + self.scroll.setWidgetResizable(True) + self.scroll.setMaximumHeight(800) + self.scroll.setVerticalScrollBarPolicy(Qt.ScrollBarPolicy.ScrollBarAsNeeded) + self.scroll.setHorizontalScrollBarPolicy(Qt.ScrollBarPolicy.ScrollBarAlwaysOff) + + self.scroll_content = QWidget() + self.scroll_layout = QVBoxLayout(self.scroll_content) + self.scroll_layout.setContentsMargins(10, 10, 10, 10) + self.scroll.setWidget(self.scroll_content) + + main_layout.addWidget(self.scroll) for idx, param_list in params_dict.items(): full_text = param_list[0].get('full_text', f"Index [{idx}]") group_label = QLabel(f"{full_text} requires parameters:") group_label.setStyleSheet("font-weight: bold; margin-top: 10px;") - layout.addWidget(group_label) + self.scroll_layout.addWidget(group_label) for param_info in param_list: label = QLabel(param_info["label"]) - layout.addWidget(label) + self.scroll_layout.addWidget(label) - line_edit = QLineEdit(self) - line_edit.setPlaceholderText(str(param_info.get("default", ""))) - layout.addWidget(line_edit) + if param_info.get("type") == list: + widget = QComboBox(self) + # Convert options to string just in case they aren't + options = [str(opt) for opt in param_info.get("options", [])] + widget.addItems(options) + + # Set default choice if it exists in the options list + default_val = str(param_info.get("default", "")) + if default_val in options: + widget.setCurrentText(default_val) + else: + widget = QLineEdit(self) + widget.setPlaceholderText(str(param_info.get("default", ""))) - self.inputs[(idx, param_info["key"])] = line_edit + self.scroll_layout.addWidget(widget) + + self.inputs[(idx, param_info["key"])] = widget # Buttons btn_layout = QHBoxLayout() @@ -113,7 +139,7 @@ class ParameterInputDialog(QDialog): cancel_btn = QPushButton("Cancel", self) btn_layout.addWidget(ok_btn) btn_layout.addWidget(cancel_btn) - layout.addLayout(btn_layout) + main_layout.addLayout(btn_layout) ok_btn.clicked.connect(self.accept) cancel_btn.clicked.connect(self.reject) @@ -131,8 +157,11 @@ class ParameterInputDialog(QDialog): Returns None if validation fails (error dialog shown). """ values = {} - for (idx, param_key), line_edit in self.inputs.items(): - text = line_edit.text().strip() + for (idx, param_key), widget in self.inputs.items(): + if isinstance(widget, QComboBox): + text = widget.currentText().strip() + else: + text = widget.text().strip() # Find param info dict param_info = None @@ -164,14 +193,15 @@ class ParameterInputDialog(QDialog): val = False else: raise ValueError(f"Invalid bool value: {text}") - elif param_type == str: + elif param_type in (str, list): val = text else: val = text # fallback except (ValueError, TypeError): + type_name = "list option" if param_type == list else param_type.__name__ self._show_error( f"Invalid input for index {idx} parameter '{param_key}': '{text}'\n" - f"Expected type: {param_type.__name__}" + f"Expected type: {type_name}" ) return None @@ -1055,7 +1085,7 @@ class FlaresBaseWidget(QWidget): class CrossGroupUIMixin: - def setup_cross_group_ui(self, index_texts): + def setup_cross_group_ui(self, index_texts, placeholder_text=""): self.group_to_paths = {} for file_path, group_name in self.group_dict.items(): @@ -1130,6 +1160,10 @@ class CrossGroupUIMixin: self.scroll_content = QWidget() self.grid_layout = QGridLayout(self.scroll_content) self.scroll_area.setWidget(self.scroll_content) + self.placeholder_label = QLabel(placeholder_text) + self.grid_layout.addWidget(self.placeholder_label, 0, 0) + self.placeholder_label.setWordWrap(True) + self.placeholder_label.setScaledContents(True) self.main_layout.addWidget(self.scroll_area) self.thumb_size = QSize(280, 180) @@ -1307,7 +1341,7 @@ class CSVUIMixin: class InterGroupUIMixin: - def setup_inter_group_ui(self, index_texts): + def setup_inter_group_ui(self, index_texts, placeholder_text=""): self.show_all_events = True self._updating_checkstates = False @@ -1366,12 +1400,16 @@ class InterGroupUIMixin: self.scroll_content = QWidget() self.grid_layout = QGridLayout(self.scroll_content) self.scroll.setWidget(self.scroll_content) + self.placeholder_label = QLabel(placeholder_text) + self.grid_layout.addWidget(self.placeholder_label, 0, 0) + self.placeholder_label.setWordWrap(True) + self.placeholder_label.setScaledContents(True) self.layout.addWidget(self.scroll) self.thumb_size = QSize(280, 180) self.showMaximized() - def get_common_request_data(self, parameterized_indexes): + def get_common_request_data(self, parameterized_indexes, json_location=None, contrast_dfs=None): selected_event = self.event_dropdown.currentText() if selected_event == "": selected_event = None @@ -1403,9 +1441,7 @@ class InterGroupUIMixin: if not selected_file_paths: print("No participants selected.") return - - # Only keep indexes 0 and 1 that need parameters - + # Inject full_text from index_texts for idx, params_list in parameterized_indexes.items(): @@ -1415,6 +1451,66 @@ class InterGroupUIMixin: indexes_needing_params = {idx: parameterized_indexes[idx] for idx in selected_indexes if idx in parameterized_indexes} + dynamic_rois = [] + + # 1. Check for the JSON file and parse ROI names + if os.path.exists(json_location): + try: + with open(json_location, 'r', encoding='utf-8') as f: + regions_data = json.load(f) + + # Extract "name" from each region under "regions_of_interest" + regions_list = regions_data.get("regions_of_interest", []) + dynamic_rois = [region["name"] for region in regions_list if "name" in region] + + except Exception as e: + # Safe log if JSON is corrupted or unreadable + print(f"Error reading ROI configurations from {json_location}: {e}") + + # Fallback to prevent UI crashes if JSON file doesn't exist or is empty + if not dynamic_rois: + dynamic_rois = ["Option 1", "Option 2"] + + dynamic_contrasts = [] + if contrast_dfs: + contrast_set = set() + for fp in selected_file_paths: + # Get the contrasts dictionary associated with this file path + file_contrasts = contrast_dfs.get(fp, {}) + for contrast_name in file_contrasts.keys(): + # If no event is selected, display all contrasts. + # If an event is selected, only keep contrasts containing the event name as a substring. + if selected_event is None or selected_event in contrast_name: + contrast_set.add(contrast_name) + + # Sort them cleanly for the UI + dynamic_contrasts = sorted(list(contrast_set)) + + # 2. Loop through the active parameters needing input and intercept 'roi_a' and 'roi_b' + for idx, params_list in indexes_needing_params.items(): + for param_info in params_list: + if param_info["key"] == "roi_a": + # Inject options list dynamically + param_info["options"] = dynamic_rois + # Default to the very first item + param_info["default"] = dynamic_rois[0] if dynamic_rois else "" + + elif param_info["key"] == "roi_b": + # Inject the same options list + param_info["options"] = dynamic_rois + # Default to the first item not taken (index 1), with safety fallbacks + if len(dynamic_rois) > 1: + param_info["default"] = dynamic_rois[1] + elif len(dynamic_rois) == 1: + param_info["default"] = dynamic_rois[0] + else: + param_info["default"] = "" + + elif param_info["key"] == "contrast_name": + param_info["options"] = dynamic_contrasts + param_info["default"] = dynamic_contrasts[0] if dynamic_contrasts else "" + + param_values = {} if indexes_needing_params: dialog = ParameterInputDialog(indexes_needing_params, parent=self) diff --git a/src/window/viewerlauncher.py b/src/window/viewerlauncher.py index 5d51b75..e3898e7 100644 --- a/src/window/viewerlauncher.py +++ b/src/window/viewerlauncher.py @@ -24,7 +24,7 @@ from src.shared.shareddata import APP_NAME class ViewerLauncherWidget(QWidget): - def __init__(self, haemo_dict, epochs_dict, cha_dict, df_ind_dict, design_matrix_dict, config_dict, fig_bytes_dict, contrast_results_dict, folding_bypass): + def __init__(self, haemo_dict, epochs_dict, cha_dict, df_ind_dict, design_matrix_dict, config_dict, fig_bytes_dict, contrast_results_dict, folding_bypass, json_location): super().__init__() self.setWindowTitle(f"Viewer Launcher - {APP_NAME.upper()}") @@ -36,7 +36,7 @@ class ViewerLauncherWidget(QWidget): ("Participant Fold Channels Viewer", ParticipantFoldChannelsWidget, [haemo_dict, cha_dict], False), ("Participant Functional Connectivity Viewer [BETA]", ParticipantFunctionalConnectivityWidget, [haemo_dict, epochs_dict], True), ("Inter-Group Functional Connectivity Viewer [BETA]", InterGroupFunctionalConnectivityWidget, [haemo_dict, group_dict, config_dict], True), - ("Inter-Group Stats Viewer", InterGroupStatsWidget, [haemo_dict, cha_dict, df_ind_dict, design_matrix_dict, contrast_results_dict, group_dict], True), + ("Inter-Group Stats Viewer", InterGroupStatsWidget, [haemo_dict, cha_dict, df_ind_dict, design_matrix_dict, contrast_results_dict, group_dict, json_location], True), ("Cross-Group Stats Viewer", CrossGroupStatsWidget, [haemo_dict, cha_dict, df_ind_dict, design_matrix_dict, contrast_results_dict, group_dict], True), ("Inter-Group Brain and Image Viewer", InterGroupBrainImageWidget, [haemo_dict, cha_dict, df_ind_dict, design_matrix_dict, contrast_results_dict, group_dict], True), ("Cross-Group Brain and Image Viewer", CrossGroupBrainImageWidget, [haemo_dict, df_ind_dict, design_matrix_dict, contrast_results_dict, group_dict], True),