6057 lines
233 KiB
Python
6057 lines
233 KiB
Python
"""
|
|
Filename: flares.py
|
|
Description: Core functionality for FLARES
|
|
|
|
Author: Tyler de Zeeuw
|
|
License: GPL-3.0
|
|
"""
|
|
|
|
# Built-in imports
|
|
import os
|
|
import sys
|
|
import platform
|
|
import threading
|
|
import logging
|
|
from io import BytesIO
|
|
from typing import Any, Optional, cast, Literal, Union
|
|
from itertools import compress
|
|
from copy import deepcopy
|
|
from multiprocessing import Queue, Pool
|
|
import os.path as op
|
|
import re
|
|
import traceback
|
|
from concurrent.futures import ProcessPoolExecutor, as_completed
|
|
from queue import Empty
|
|
import time
|
|
import multiprocessing as mp
|
|
import itertools
|
|
|
|
# External library imports
|
|
import matplotlib.pyplot as plt
|
|
import matplotlib.colors as mcolors
|
|
from matplotlib.figure import Figure
|
|
from matplotlib.axes import Axes
|
|
from matplotlib.colors import LinearSegmentedColormap
|
|
from matplotlib.lines import Line2D
|
|
|
|
import numpy as np
|
|
from numpy.typing import NDArray
|
|
from numpy import float64, floating
|
|
|
|
import pandas as pd
|
|
from pandas import DataFrame
|
|
|
|
import seaborn as sns
|
|
import h5py
|
|
|
|
from nilearn.plotting import plot_design_matrix # type: ignore
|
|
from nilearn.glm.regression import OLSModel
|
|
|
|
import statsmodels.formula.api as smf # type: ignore
|
|
from statsmodels.stats.multitest import multipletests
|
|
|
|
from scipy.spatial.distance import cdist
|
|
from scipy.signal import welch, butter, filtfilt # type: ignore
|
|
from scipy.stats import pearsonr, zscore, t
|
|
|
|
import pywt # type: ignore
|
|
import neurokit2 as nk # type: ignore
|
|
|
|
# Backend visualization needed to be defined for pyinstaller
|
|
import pyvistaqt # type: ignore
|
|
import vtkmodules.util.data_model
|
|
import vtkmodules.util.execution_model
|
|
import xlrd
|
|
|
|
# External library imports for mne
|
|
from mne import (
|
|
EvokedArray, SourceEstimate, Info, Epochs, Label, Annotations,
|
|
events_from_annotations, read_source_spaces,
|
|
stc_near_sensors, pick_types, grand_average, get_config, set_config, read_labels_from_annot
|
|
) # type: ignore
|
|
from mne.source_space import SourceSpaces
|
|
from mne.transforms import Transform # type: ignore
|
|
from mne.io import BaseRaw, RawArray, read_raw_snirf # type: ignore
|
|
from mne.preprocessing.nirs import (
|
|
beer_lambert_law, optical_density,
|
|
temporal_derivative_distribution_repair,
|
|
source_detector_distances, short_channels
|
|
) # type: ignore
|
|
from mne.viz import Brain, plot_events, plot_evoked_topo, plot_compare_evokeds
|
|
from mne.filter import filter_data # type: ignore
|
|
from mne.utils import _check_fname, _validate_type, warn
|
|
from mne.channels import make_standard_montage
|
|
from mne.datasets.sample import data_path
|
|
|
|
from mne_nirs.visualisation import plot_glm_group_topo # type: ignore
|
|
from mne_nirs.channels import get_long_channels, get_short_channels # type: ignore
|
|
from mne_nirs.experimental_design import make_first_level_design_matrix # type: ignore
|
|
from mne_nirs.statistics import run_glm, statsmodels_to_results # type: ignore
|
|
from mne_nirs.signal_enhancement import (
|
|
enhance_negative_correlation, short_channel_regression
|
|
) # type: ignore
|
|
from mne_nirs.io.fold import fold_channel_specificity # type: ignore
|
|
from mne_nirs.preprocessing import peak_power # type: ignore
|
|
from mne_nirs.statistics._glm_level_first import RegressionResults # type: ignore
|
|
|
|
from mne_connectivity.viz import plot_connectivity_circle
|
|
from mne_connectivity import envelope_correlation, spectral_connectivity_epochs, spectral_connectivity_time
|
|
|
|
|
|
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
|
|
|
|
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)
|
|
"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)
|
|
"Multiple": "orangered" # Failed 2+ categories
|
|
}
|
|
|
|
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."""
|
|
return PRIMARY_COLORS.get(label, COMBINATION_COLOR)
|
|
|
|
|
|
DOWNSAMPLE: bool
|
|
DOWNSAMPLE_FREQUENCY: int
|
|
|
|
TRIM: bool
|
|
SECONDS_TO_KEEP: float
|
|
|
|
OPTODE_PLACEMENT: bool
|
|
SHOW_OPTODE_NAMES: bool
|
|
|
|
SHORT_CHANNELS: bool
|
|
SHORT_CHANNELS_THRESHOLD: float
|
|
LONG_CHANNELS_THRESHOLD: float
|
|
|
|
HEART_RATE: bool
|
|
SECONDS_TO_STRIP_HR: int
|
|
MAX_LOW_HR: int
|
|
MAX_HIGH_HR: int
|
|
SMOOTHING_WINDOW_HR: int
|
|
HEART_RATE_WINDOW: int
|
|
|
|
SCI: bool
|
|
SCI_TIME_WINDOW: int
|
|
SCI_THRESHOLD: float
|
|
|
|
SNR: bool
|
|
# SNR_TIME_WINDOW : int #TODO: is this needed?
|
|
SNR_THRESHOLD: float
|
|
|
|
PSP: bool
|
|
PSP_TIME_WINDOW: int
|
|
PSP_THRESHOLD: float
|
|
|
|
COEFF_VAR: bool
|
|
COEFF_VAR_THRESHOLD: int
|
|
|
|
MAD: bool
|
|
MAD_THRESHOLD: int
|
|
|
|
PSD_NOISE: bool
|
|
TARGET_FREQ_DIV: int
|
|
DB_LIMIT: int
|
|
|
|
SENSOR_DROPOUT: bool
|
|
SENSOR_DROPOUT_VARIANCE_THRESHOLD: float
|
|
|
|
BAD_CHANNELS_HANDLING: str
|
|
MAX_DIST: float
|
|
MIN_NEIGHBORS: int
|
|
MAX_BAD_CHANNELS: int
|
|
|
|
TDDR: bool
|
|
|
|
WAVELET: bool
|
|
IQR: float
|
|
WAVELET_TYPE: str
|
|
WAVELET_LEVEL: int
|
|
|
|
ENHANCE_NEGATIVE_CORRELATION: bool
|
|
|
|
FILTER: bool
|
|
L_FREQ: float
|
|
H_FREQ: float
|
|
L_TRANS_BANDWIDTH: float
|
|
H_TRANS_BANDWIDTH: float
|
|
|
|
EPOCH_HANDLING: str
|
|
MAX_SHIFT: int
|
|
T_MIN: int
|
|
T_MAX: int
|
|
|
|
RESAMPLE: bool
|
|
RESAMPLE_FREQ: int
|
|
STIM_DUR: float
|
|
HRF_MODEL: str
|
|
DRIFT_MODEL: str
|
|
HIGH_PASS: float
|
|
DRIFT_ORDER: int
|
|
FIR_DELAYS: range
|
|
MIN_ONSET: int
|
|
OVERSAMPLING: int
|
|
REMOVE_EVENTS: list
|
|
SHORT_CHANNEL_REGRESSION: bool
|
|
|
|
NOISE_MODEL: str
|
|
BINS: int
|
|
N_JOBS: int
|
|
|
|
JSON_LOCATION: str
|
|
|
|
TIME_WINDOW_START: int
|
|
TIME_WINDOW_END: int
|
|
MAX_WORKERS: int
|
|
VERBOSITY: bool
|
|
|
|
AGE: int = 25 # Assume 25 if not set from the GUI. This will result in a reasonable PPF
|
|
GENDER: str = ""
|
|
GROUP: str = "Default"
|
|
|
|
FOLDING_BYP: bool = False
|
|
|
|
# These are parameters that are required for the analysis
|
|
REQUIRED_KEYS: dict[str, Any] = {
|
|
|
|
# "SECONDS_TO_STRIP": int,
|
|
"DOWNSAMPLE": bool,
|
|
"DOWNSAMPLE_FREQUENCY": int,
|
|
|
|
"TRIM": bool,
|
|
"SECONDS_TO_KEEP": float,
|
|
|
|
"OPTODE_PLACEMENT": bool,
|
|
|
|
"HEART_RATE": bool,
|
|
|
|
"SCI": bool,
|
|
"SCI_TIME_WINDOW": int,
|
|
"SCI_THRESHOLD": float,
|
|
|
|
"SNR": bool,
|
|
# SNR_TIME_WINDOW : int
|
|
"SNR_THRESHOLD": float,
|
|
|
|
"PSP": bool,
|
|
"PSP_TIME_WINDOW": int,
|
|
"PSP_THRESHOLD": float,
|
|
|
|
"SHORT_CHANNELS": bool,
|
|
"SHORT_CHANNELS_THRESHOLD": float,
|
|
"LONG_CHANNELS_THRESHOLD": float,
|
|
|
|
|
|
"REMOVE_EVENTS": list,
|
|
"TIME_WINDOW_START": int,
|
|
"TIME_WINDOW_END": int,
|
|
"L_FREQ": float,
|
|
"H_FREQ": float,
|
|
|
|
"EPOCH_HANDLING": str,
|
|
"T_MIN": int,
|
|
"T_MAX": int,
|
|
|
|
"TDDR": bool,
|
|
"WAVELET": bool,
|
|
"IQR": float,
|
|
"WAVELET_TYPE": str,
|
|
"WAVELET_LEVEL": int,
|
|
"FILTER": bool,
|
|
"DRIFT_MODEL": str,
|
|
# "REJECT_PAIRS": bool,
|
|
# "FORCE_DROP_ANNOTATIONS": list,
|
|
# "FILTER_LOW_PASS": float,
|
|
# "FILTER_HIGH_PASS": float,
|
|
# "EPOCH_PAIR_TOLERANCE_WINDOW": int,
|
|
}
|
|
|
|
|
|
# audit_log = logging.getLogger("memory_audit")
|
|
# audit_log.setLevel(logging.INFO)
|
|
# audit_log.propagate = False # This prevents it from talking to other loggers
|
|
|
|
# # 2. Add a file handler specifically for this audit logger
|
|
# if not audit_log.handlers:
|
|
# fh = logging.FileHandler('flares_memory_audit.log')
|
|
# fh.setFormatter(logging.Formatter('%(asctime)s | PID: %(process)d | %(message)s'))
|
|
# audit_log.addHandler(fh)
|
|
|
|
# def get_mem_mb():
|
|
# return psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024
|
|
|
|
|
|
|
|
class ProcessingError(Exception):
|
|
def __init__(self, message: str = "Something went wrong!"):
|
|
self.message = message
|
|
super().__init__(self.message)
|
|
|
|
|
|
# Ensure that we are working in the directory of this file
|
|
script_dir = os.path.dirname(os.path.abspath(__file__))
|
|
os.chdir(script_dir)
|
|
|
|
PLATFORM_NAME = platform.system().lower()
|
|
|
|
# Configure logging to file with timestamps and realtime flush
|
|
if PLATFORM_NAME == 'darwin':
|
|
logging.basicConfig(
|
|
filename=os.path.join(os.path.dirname(sys.executable), "../../../fnirs_analysis.log"), # Needed to get out of the bundled application
|
|
level=logging.INFO,
|
|
format='%(asctime)s - %(processName)s - %(levelname)s - %(message)s',
|
|
datefmt='%Y-%m-%d %H:%M:%S',
|
|
filemode='a'
|
|
)
|
|
|
|
else:
|
|
logging.basicConfig(
|
|
filename='fnirs_analysis.log',
|
|
level=logging.INFO,
|
|
format='%(asctime)s - %(processName)s - %(levelname)s - %(message)s',
|
|
datefmt='%Y-%m-%d %H:%M:%S',
|
|
filemode='a'
|
|
)
|
|
|
|
logger = logging.getLogger()
|
|
|
|
|
|
|
|
|
|
def set_config_me(config: dict[str, Any]) -> None:
|
|
"""
|
|
Validates and applies the given configuration dictionary.
|
|
|
|
Parameters
|
|
----------
|
|
config : dict[str, Any]
|
|
Dictionary containing configuration keys and their values.
|
|
"""
|
|
logger.info(f"[DEBUG] set_config called")
|
|
|
|
globals().update(config)
|
|
|
|
|
|
def set_metadata(file_path, metadata: dict[str, Any]) -> None:
|
|
"""
|
|
Validates and applies the given configuration dictionary.
|
|
|
|
Parameters
|
|
----------
|
|
config : dict[str, Any]
|
|
Dictionary containing configuration keys and their values.
|
|
"""
|
|
logger.info(f"[DEBUG] set_metadata called")
|
|
|
|
globals()['AGE'] = 25
|
|
globals()['GENDER'] = ""
|
|
globals()['GROUP'] = "Default"
|
|
|
|
if metadata.get(file_path) is not None:
|
|
|
|
file_metadata = metadata.get(file_path, {})
|
|
|
|
for key in ("AGE", "GENDER", "GROUP"):
|
|
val = file_metadata.get(key, None)
|
|
if val not in (None, '', [], {}, ()): # check for "empty" values
|
|
globals()[key] = val
|
|
|
|
def gui_entry(config: dict[str, Any], gui_queue: Queue, progress_queue: Queue, ack_queue: Queue) -> None:
|
|
try:
|
|
file_paths = config['SNIRF_FILES']
|
|
file_params = config['PARAMS']
|
|
file_metadata = config['METADATA']
|
|
max_workers = file_params.get("MAX_WORKERS", int(os.cpu_count()/4))
|
|
|
|
results = process_multiple_participants(
|
|
file_paths, file_params, file_metadata, progress_queue, gui_queue, max_workers
|
|
)
|
|
|
|
gui_queue.put({"type": "FINISHED_SUCCESSFULLY", "success": True})
|
|
|
|
try:
|
|
print("CHILD: Waiting for GUI acknowledgment...")
|
|
ack_queue.get(timeout=10)
|
|
except:
|
|
print("CHILD: Ack timeout, exiting anyway.")
|
|
|
|
except Exception as e:
|
|
gui_queue.put({
|
|
"type": "FINISHED_SUCCESSFULLY",
|
|
"success": False,
|
|
"error": str(e),
|
|
"traceback": traceback.format_exc()
|
|
})
|
|
|
|
finally:
|
|
pass
|
|
|
|
|
|
def process_participant_worker(file_path, file_params, file_metadata, result_queue, progress_queue):
|
|
file_name = os.path.basename(file_path)
|
|
|
|
try:
|
|
# 1. Setup
|
|
set_config_me(file_params)
|
|
set_metadata(file_path, file_metadata)
|
|
|
|
def progress_callback(step_idx):
|
|
if progress_queue:
|
|
# We use put_nowait to prevent the worker from hanging on a full queue
|
|
try:
|
|
progress_queue.put_nowait(('progress', file_path, step_idx))
|
|
except: pass
|
|
|
|
# 2. Process
|
|
result = process_participant(file_path, progress_callback=progress_callback)
|
|
|
|
# 3. Report Success
|
|
result_queue.put((file_path, result, None))
|
|
|
|
except Exception as e:
|
|
try:
|
|
result_queue.put((file_path, None, str(e)))
|
|
except: pass
|
|
|
|
finally:
|
|
try:
|
|
sys.stdout.flush()
|
|
sys.stderr.flush()
|
|
# Give the queue thread a moment to send the data before we vanish
|
|
time.sleep(0.2)
|
|
except: pass
|
|
sys.exit(0)
|
|
|
|
|
|
def process_multiple_participants(file_paths, file_params, file_metadata, progress_queue=None, gui_queue=None, max_workers=6):
|
|
#audit_log.info(f"--- SESSION START: {len(file_paths)} files ---")
|
|
|
|
pending_files = list(file_paths)
|
|
active_processes = []
|
|
results_by_file = {}
|
|
|
|
manager = mp.Manager()
|
|
result_queue = manager.Queue()
|
|
|
|
try:
|
|
while pending_files or active_processes:
|
|
for p, f_path in active_processes[:]:
|
|
if not p.is_alive():
|
|
p.join(timeout=0.1)
|
|
active_processes.remove((p, f_path))
|
|
|
|
while len(active_processes) < max_workers and pending_files:
|
|
file_path = pending_files.pop(0)
|
|
p = mp.Process(
|
|
target=process_participant_worker,
|
|
args=(file_path, file_params, file_metadata, result_queue, progress_queue)
|
|
)
|
|
p.start()
|
|
active_processes.append((p, file_path))
|
|
|
|
if progress_queue:
|
|
while not progress_queue.empty():
|
|
try:
|
|
prog_msg = progress_queue.get_nowait()
|
|
if gui_queue:
|
|
# Forward straight to GUI
|
|
gui_queue.put(prog_msg, timeout=0.1)
|
|
except: break
|
|
|
|
while not result_queue.empty():
|
|
try:
|
|
res_path, result, error = result_queue.get(timeout=0.01)
|
|
if gui_queue:
|
|
gui_queue.put({
|
|
"type": "file_done",
|
|
"file": res_path,
|
|
"success": error is None,
|
|
"result": result if error is None else None,
|
|
"error": error if error else None
|
|
}, timeout=2)
|
|
else:
|
|
results_by_file[res_path] = result
|
|
except:
|
|
break
|
|
|
|
if not pending_files and not active_processes:
|
|
print("DEBUG: Loop finished naturally.")
|
|
break
|
|
|
|
# should no longer hang because pipes are being drained
|
|
time.sleep(0.1)
|
|
|
|
except Exception as e:
|
|
print(f"MAIN LOOP ERROR: {e}")
|
|
finally:
|
|
# Cleanup
|
|
for p, _ in active_processes:
|
|
try:
|
|
if p.is_alive():
|
|
p.terminate()
|
|
p.join(timeout=0.1)
|
|
except: pass
|
|
|
|
if manager is not None:
|
|
try:
|
|
manager.shutdown()
|
|
except: pass
|
|
|
|
return results_by_file
|
|
|
|
|
|
|
|
def markbad(data, ax, ch_names: list[str]) -> None:
|
|
"""
|
|
Add a strikethrough to a plot for channels marked as bad.
|
|
|
|
Parameters
|
|
----------
|
|
data : BaseRaw
|
|
The loaded data object to process.
|
|
ax : Axes
|
|
Matplotlib Axes object where the strikethrough lines will be drawn.
|
|
ch_names : list[str]
|
|
List of channel names corresponding to the y-axis of the plot.
|
|
"""
|
|
|
|
# Iterate over all the channels
|
|
for i, ch in enumerate(ch_names):
|
|
|
|
# If it is marked as bad, place a strikethrough on the channel
|
|
if ch in data.info["bads"]:
|
|
ax.axhline(i + 0.5, ls="solid", lw=4, color="black", zorder=10) # type: ignore
|
|
|
|
|
|
|
|
def plot_timechannel_quality_metrics(data, scores, times: list[tuple[float]], color_stops: tuple[list[float], list[float]], threshold: float, title: Optional[str] = None):
|
|
|
|
"""
|
|
Generate two heatmaps visualizing channel quality metrics over time.
|
|
|
|
Parameters
|
|
----------
|
|
data : BaseRaw
|
|
The loaded data object to process.
|
|
scores : NDArray[float64]
|
|
A 2D array of quality scores for each channel over time.
|
|
times : list[tuple[float]]
|
|
List of time boundaries used to label each score column.
|
|
color_stops : tuple[list[float], list[float]]
|
|
Two lists of color values for custom colormaps.
|
|
threshold : float,
|
|
Threshold value for the color bar.
|
|
title : Optional[str], optional
|
|
Base title for the figures, (default is None).
|
|
|
|
Returns
|
|
-------
|
|
tuple[Figure, Figure]
|
|
- Figure: Heatmap of all scores across channels and time.
|
|
- Figure: Binary heatmap showing only scores above the threshold.
|
|
"""
|
|
|
|
# Get only the hbo / hbr channels once as we dont need to see the same results twice
|
|
half_ch = len(getattr(data, "ch_names")) // 2
|
|
ch_names = getattr(data, "ch_names")[:half_ch]
|
|
scores = scores[:half_ch, :]
|
|
|
|
# Extract rounded time points to use as column headers
|
|
cols = [np.round(t[0]) for t in times]
|
|
n_chans = len(ch_names)
|
|
vsize = 0.2 * n_chans
|
|
|
|
# Create the first figure
|
|
fig1, ax1 = plt.subplots(figsize=(10, vsize), layout="constrained") # type: ignore
|
|
fig1.suptitle(title + " - All Scores", fontsize=16, fontweight="bold") # type: ignore
|
|
|
|
# Create a DataFrame to structure data for the heatmap
|
|
data_to_plot = DataFrame(
|
|
data=scores,
|
|
columns=pd.Index(cols, name="Time (s)"),
|
|
index=pd.Index(ch_names, name="Channel"),
|
|
)
|
|
|
|
# Define a custom colormap using provided color stops and base colors
|
|
base_colors = ['red', 'red', 'yellow', 'green', 'green']
|
|
colors = list(zip(color_stops[0], base_colors[:len(color_stops[0])]))
|
|
cmap = mcolors.LinearSegmentedColormap.from_list('gyr', colors)
|
|
|
|
# Plot heatmap of scores
|
|
sns.heatmap( # type: ignore
|
|
data=data_to_plot,
|
|
cmap=cmap,
|
|
vmin=0,
|
|
vmax=1,
|
|
cbar_kws=dict(label="Score"),
|
|
ax=ax1,
|
|
)
|
|
|
|
# Add vertical dashed lines at each time boundary, sit the title, and place a black strikethrough through a bad channel
|
|
for x in range(1, len(times)):
|
|
ax1.axvline(x, ls="dashed", lw=0.25, dashes=(25, 15), color="gray") # type: ignore
|
|
ax1.set_title("All Scores", fontweight="bold") # type: ignore
|
|
markbad(data, ax1, ch_names)
|
|
|
|
# Calculate average score per channel and annotate to the right of the heatmap
|
|
avg_sci_subset: pd.Series[float] = data_to_plot.mean(axis=1) # type: ignore
|
|
norm = mcolors.Normalize(vmin=0, vmax=1)
|
|
text_x = data_to_plot.shape[1] + 0.5
|
|
for i, val in enumerate(avg_sci_subset):
|
|
color = cmap(norm(val))
|
|
ax1.text( # type: ignore
|
|
text_x,
|
|
i + 0.5,
|
|
f"{val:.3f}",
|
|
va='center',
|
|
ha='left',
|
|
fontsize=9,
|
|
color=color
|
|
)
|
|
ax1.set_xlim(right=text_x + 1.5)
|
|
|
|
plt.close(fig1)
|
|
|
|
# Create the second figure
|
|
fig2, ax2 = plt.subplots(figsize=(10, vsize), layout="constrained") # type: ignore
|
|
fig2.suptitle(title + " - Scores Above Threshold", fontsize=16, fontweight="bold") # type: ignore
|
|
|
|
# Create a DataFrame to structure data for the heatmap
|
|
data_to_plot = DataFrame(
|
|
data=scores > threshold,
|
|
columns=pd.Index(cols, name="Time (s)"),
|
|
index=pd.Index(ch_names, name="Channel"),
|
|
)
|
|
|
|
# Define a custom colormap using provided color stops and base colors
|
|
base_colors = ['red', 'red', 'white', 'white']
|
|
colors = list(zip(color_stops[1], base_colors[:len(color_stops[1])]))
|
|
cmap = mcolors.LinearSegmentedColormap.from_list('gyr', colors)
|
|
|
|
# Plot heatmap of scores
|
|
sns.heatmap( # type: ignore
|
|
data=data_to_plot,
|
|
vmin=0,
|
|
vmax=1,
|
|
cmap=cmap,
|
|
cbar_kws=dict(label="Score"),
|
|
ax=ax2,
|
|
)
|
|
|
|
# Add vertical dashed lines at each time boundary, sit the title, and place a black strikethrough through a bad channel
|
|
for x in range(1, len(times)):
|
|
ax2.axvline(x, ls="dashed", lw=0.25, dashes=(25, 15), color="gray") # type: ignore
|
|
ax2.set_title("Scores > Threshold", fontweight="bold") # type: ignore
|
|
markbad(data, ax2, ch_names)
|
|
|
|
plt.close(fig2)
|
|
|
|
return fig1, fig2
|
|
|
|
|
|
def scalp_coupling_index_windowed_raw(data, time_window: float = 3.0, l_freq: float = 0.7, h_freq: float = 1.5, l_trans_bandwidth: float = 0.3, h_trans_bandwidth: float = 0.3):
|
|
"""
|
|
Compute windowed scalp coupling index (SCI) across fNIRS channels.
|
|
|
|
Parameters
|
|
----------
|
|
data : BaseRaw
|
|
The loaded data object to process.
|
|
time_window : float, optional
|
|
Length of each time window in seconds (default is 3.0).
|
|
l_freq : float, optional
|
|
Low cutoff frequency for filtering in Hz (default is 0.7).
|
|
h_freq : float, optional
|
|
High cutoff frequency for filtering in Hz (default is 1.5).
|
|
l_trans_bandwidth : float, optional
|
|
Transition bandwidth for the low cutoff in Hz (default is 0.3).
|
|
h_trans_bandwidth : float, optional
|
|
Transition bandwidth for the high cutoff in Hz (default is 0.3).
|
|
|
|
Returns
|
|
-------
|
|
tuple[BaseRaw, NDArray[float64], list[tuple[float, float]]]
|
|
- BaseRaw: The original data object (unchanged). Ensures compatibility with peak_power().
|
|
- NDArray[float64]: Correlation scores for each channel and time window.
|
|
- list[tuple[float, float]]: Time intervals for each window in seconds.
|
|
"""
|
|
|
|
# Pick only fNIRS channels and sort them by channel name
|
|
picks: NDArray[np.intp] = pick_types(cast(Info, data.info), fnirs=True) # type: ignore
|
|
picks = picks[np.argsort([getattr(data, "ch_names")[pick] for pick in picks])]
|
|
|
|
# FIXME: This may happen if the heart rate calculation tries to set a value way too low
|
|
if l_freq < 0.3:
|
|
l_freq = 0.3
|
|
|
|
# Band-pass filter the selected fNIRS channels
|
|
filtered_data = filter_data(
|
|
getattr(data, "_data"),
|
|
getattr(data, "info")["sfreq"],
|
|
l_freq,
|
|
h_freq,
|
|
picks=picks,
|
|
verbose=False,
|
|
l_trans_bandwidth=l_trans_bandwidth, # type: ignore
|
|
h_trans_bandwidth=h_trans_bandwidth, # type: ignore
|
|
)
|
|
|
|
# Calculate number of samples per time window, the total number of windows, and prepare output variables
|
|
window_samples = int(np.ceil(time_window * getattr(data, "info")["sfreq"]))
|
|
n_windows = int(np.floor(len(data) / window_samples))
|
|
scores = np.zeros((len(picks), n_windows))
|
|
times: list[tuple[float, float]] = []
|
|
|
|
# Slide through the data in windows to compute scalp coupling index (SCI)
|
|
for window in range(n_windows):
|
|
start_sample = int(window * window_samples)
|
|
end_sample = start_sample + window_samples
|
|
end_sample = np.min([end_sample, len(data) - 1])
|
|
|
|
# Track time boundaries for each window
|
|
t_start = getattr(data, "times")[start_sample]
|
|
t_stop = getattr(data, "times")[end_sample]
|
|
times.append((t_start, t_stop))
|
|
|
|
# Iterate through channels in pairs (hbo, hbr). This requires them to be sorted by channel name
|
|
for ii in range(0, len(picks), 2):
|
|
c1 = filtered_data[picks[ii]][start_sample:end_sample]
|
|
c2 = filtered_data[picks[ii + 1]][start_sample:end_sample]
|
|
|
|
# Ensure the correlation data is valid
|
|
if np.std(c1) == 0 or np.std(c2) == 0 or np.any(np.isnan(c1)) or np.any(np.isnan(c2)):
|
|
c = 0
|
|
else:
|
|
c = np.corrcoef(c1, c2)[0][1]
|
|
|
|
# Assign the computed correlation to both channels in the pair
|
|
scores[ii, window] = c
|
|
scores[ii + 1, window] = c
|
|
|
|
scores = scores[np.argsort(picks)]
|
|
|
|
return data, scores, times
|
|
|
|
def calculate_scalp_coupling(data, l_freq: float = 0.7, h_freq: float = 1.5):
|
|
"""
|
|
Calculate the scalp coupling index (SCI) and identify bad channels based on a threshold.
|
|
|
|
Parameters
|
|
----------
|
|
data : BaseRaw
|
|
The loaded data object to process.
|
|
l_freq : float, optional
|
|
Low cutoff frequency for bandpass filtering in Hz (default is 0.7).
|
|
h_freq : float, optional
|
|
High cutoff frequency for bandpass filtering in Hz (default is 1.5)
|
|
|
|
Returns
|
|
-------
|
|
tuple[list[str], Figure, Figure]
|
|
- list[str]: Channel names identified as bad based on SCI threshold.
|
|
- Figure: Heatmap of all SCI scores across time and channels.
|
|
- Figure: Binary heatmap of SCI scores exceeding the threshold.
|
|
"""
|
|
|
|
print("Calculating scalp coupling index...")
|
|
|
|
# Compute the SCI
|
|
_, scores, times = scalp_coupling_index_windowed_raw(data, time_window=SCI_TIME_WINDOW, l_freq=l_freq, h_freq=h_freq)
|
|
|
|
# Identify channels that don't meet the provided threshold
|
|
print("Identifying channels that do not meet the threshold...")
|
|
sci = scores.mean(axis=1)
|
|
data.info["bads"] = list(compress(cast(list[str], getattr(data, "ch_names")), sci < SCI_THRESHOLD))
|
|
|
|
# Determine the colors based on the threshold, and create the figures
|
|
print("Creating the figures...")
|
|
color_stops = ([0.0, SCI_THRESHOLD, SCI_THRESHOLD+0.1, 0.8, 1.0], [0.0, SCI_THRESHOLD, SCI_THRESHOLD, 1.0])
|
|
fig1, fig2 = plot_timechannel_quality_metrics(data, scores, times, color_stops, SCI_THRESHOLD, "Scalp Coupling Index")
|
|
|
|
print("Successfully calculated scalp coupling index.")
|
|
|
|
return list(compress(cast(list[str], getattr(data, "ch_names")), sci < SCI_THRESHOLD)), fig1, fig2
|
|
|
|
|
|
|
|
def calculate_signal_noise_ratio(data):
|
|
"""
|
|
Calculates the signal-to-noise ratio (SNR) for each channel and identifies those below a defined threshold.
|
|
|
|
Parameters
|
|
----------
|
|
data : BaseRaw
|
|
The loaded data object to process.
|
|
|
|
Returns
|
|
-------
|
|
tuple[list[str], Figure]
|
|
- list[str]: A list of channel names that fall below the SNR threshold and are considered bad.
|
|
- Figure: A matplotlib Figure showing the channels' SNR values.
|
|
"""
|
|
|
|
print("Calculating signal to noise ratio...")
|
|
|
|
# Compute the signal-to-noise ratio values
|
|
print("Computing the signal to noise power...")
|
|
signal_band=(0.01, 0.5)
|
|
noise_band=(1.0, 10.0)
|
|
data_signal = data.copy().filter(*signal_band, verbose=False) #type: ignore
|
|
data_noise = data.copy().filter(*noise_band, verbose=False) #type: ignore
|
|
signal_power = np.mean(data_signal.get_data()**2, axis=1) #type: ignore
|
|
noise_power = np.mean(data_noise.get_data()**2, axis=1) #type: ignore
|
|
|
|
# Calculate the snr using the standard formula for dB
|
|
snr = 10 * np.log10(signal_power / (noise_power + np.finfo(float).eps))
|
|
|
|
# TODO: Understand what this does
|
|
groups: dict[str, list[str]] = {}
|
|
for ch in getattr(data, "ch_names"):
|
|
# Look for the space in the channel names and remove the characters after
|
|
# This is so we can get both oxy and deoxy to remove, as they will have the same source and detector
|
|
base = ch.rsplit(' ', 1)[0]
|
|
groups.setdefault(base, []).append(ch) # type: ignore
|
|
|
|
# If any of the channels do not meet our threshold, they will get inserted into the bad_channels set
|
|
bad_channels: set[str] = set()
|
|
for base, ch_list in groups.items():
|
|
if any(s < SNR_THRESHOLD for s, ch in zip(snr, getattr(data, "ch_names")) if ch in ch_list):
|
|
bad_channels.update(ch_list)
|
|
|
|
# Design and create the figure
|
|
print("Creating the figure...")
|
|
snr_fig, ax = plt.subplots(figsize=(12, 4), layout="constrained") # type: ignore
|
|
colors = [(0/20, 'red'), (SNR_THRESHOLD/20, 'red'), ((SNR_THRESHOLD+.5)/20, 'yellow'), ((SNR_THRESHOLD+1)/20, 'green'), (20/20, 'green')]
|
|
cmap = LinearSegmentedColormap.from_list('custom_snr_cmap', colors)
|
|
norm = mcolors.Normalize(vmin=0, vmax=20)
|
|
scatter = ax.scatter(range(len(snr)), snr, c=snr, cmap=cmap, alpha=0.8, s=100, norm=norm) # type: ignore
|
|
ax.set(xlabel="Channel Number", ylabel="Signal-to-Noise Ratio (dB)", xlim=[0, len(snr)], ylim=[0, 20])
|
|
ax.axhline(SNR_THRESHOLD, color='black', linestyle='--', alpha=0.3, linewidth=1) # type: ignore
|
|
cbar = snr_fig.colorbar(scatter, ax=ax, label="SNR Thresholds (dB)") # type: ignore
|
|
cbar.set_ticks([0, SNR_THRESHOLD, SNR_THRESHOLD+1, 20]) # type: ignore
|
|
cbar.set_ticklabels(['0', str(SNR_THRESHOLD), str(SNR_THRESHOLD+1), '20']) # type: ignore
|
|
|
|
plt.close()
|
|
|
|
print("Successfully calculated signal to noise ratio.")
|
|
|
|
return list(bad_channels), snr_fig
|
|
|
|
|
|
|
|
def build_fnirs_adjacency(raw, threshold_meters=0.03):
|
|
"""Build an adjacency dictionary for fNIRS channels using 3D distance."""
|
|
# Extract channel positions
|
|
ch_locs = []
|
|
ch_names = []
|
|
|
|
for ch in raw.info['chs']:
|
|
loc = ch['loc'][:3] # Get x, y, z coordinates
|
|
if not np.isnan(loc).any():
|
|
ch_locs.append(loc)
|
|
ch_names.append(ch['ch_name'])
|
|
|
|
ch_locs = np.array(ch_locs)
|
|
|
|
# Compute pairwise distances
|
|
dists = cdist(ch_locs, ch_locs)
|
|
|
|
# Build adjacency dictionary
|
|
adjacency = {}
|
|
for i, ch_name in enumerate(ch_names):
|
|
neighbors = [ch_names[j] for j in range(len(ch_names))
|
|
if 0 < dists[i, j] < threshold_meters]
|
|
adjacency[ch_name] = neighbors
|
|
|
|
return adjacency
|
|
|
|
|
|
|
|
def get_hbo_hbr_picks(raw):
|
|
# Pick all fNIRS channels
|
|
fnirs_picks = pick_types(raw.info, fnirs=True, exclude=[])
|
|
|
|
# Extract wavelengths from channel names (expecting something like 'S6_D4 763' or 'S6_D4 841')
|
|
wavelengths = []
|
|
for idx in fnirs_picks:
|
|
ch_name = raw.ch_names[idx]
|
|
# Extract last 3 digits from channel name using regex
|
|
match = re.search(r'(\d{3})$', ch_name)
|
|
if match:
|
|
wavelengths.append(int(match.group(1)))
|
|
else:
|
|
raise ValueError(f"Channel name '{ch_name}' does not end with 3 digits.")
|
|
|
|
wavelengths = np.array(wavelengths)
|
|
unique_wavelengths = np.unique(wavelengths)
|
|
if len(unique_wavelengths) != 2:
|
|
raise RuntimeError(f"Expected exactly 2 distinct wavelengths, found {unique_wavelengths}")
|
|
|
|
# Determine which is HbO (larger) and which is HbR (smaller)
|
|
hbr_wl = unique_wavelengths.min()
|
|
hbo_wl = unique_wavelengths.max()
|
|
|
|
print(f"HbR wavelength: {hbr_wl}, HbO wavelength: {hbo_wl}")
|
|
|
|
# Find picks corresponding to each wavelength
|
|
hbr_picks = [fnirs_picks[i] for i, wl in enumerate(wavelengths) if wl == hbr_wl]
|
|
hbo_picks = [fnirs_picks[i] for i, wl in enumerate(wavelengths) if wl == hbo_wl]
|
|
|
|
print(f"Found {len(hbr_picks)} HbR channels and {len(hbo_picks)} HbO channels.")
|
|
|
|
return hbo_picks, hbr_picks, hbo_wl, hbr_wl
|
|
|
|
|
|
def interpolate_fNIRS_bads_weighted_average(raw, max_dist=0.03, min_neighbors=2):
|
|
"""
|
|
Interpolate bad fNIRS channels using a distance-weighted average of nearby good channels.
|
|
|
|
Parameters
|
|
----------
|
|
raw : mne.io.Raw
|
|
The raw fNIRS data with bads marked in raw.info['bads'].
|
|
max_dist : float
|
|
Maximum distance (in meters) to consider for neighboring good channels.
|
|
min_neighbors : int
|
|
Minimum number of neighbors required to interpolate a bad channel.
|
|
|
|
Returns
|
|
-------
|
|
raw : mne.io.Raw
|
|
Modified raw object with bads interpolated (in-place).
|
|
"""
|
|
|
|
print("Finding fNIRS channels...")
|
|
hbo_picks, hbr_picks, hbo_wl, hbr_wl = get_hbo_hbr_picks(raw)
|
|
|
|
|
|
if len(hbo_picks) != len(hbr_picks):
|
|
raise RuntimeError("Number of HbO and HbR channels must be the same.")
|
|
|
|
# Base names without wavelength for pairing
|
|
def base_name(ch_name):
|
|
# Strip last 4 chars assuming format ' <wavelength>'
|
|
# e.g. "S6_D6 841" -> "S6_D6"
|
|
return ch_name[:-4]
|
|
|
|
hbo_names = [base_name(raw.ch_names[i]) for i in hbo_picks]
|
|
hbr_names = [base_name(raw.ch_names[i]) for i in hbr_picks]
|
|
|
|
# Sanity check: pairs must match
|
|
for i in range(len(hbo_names)):
|
|
if hbo_names[i] != hbr_names[i]:
|
|
raise RuntimeError(f"Channel pairs do not match: {hbo_names[i]} vs {hbr_names[i]}")
|
|
|
|
# Identify bad pairs if either channel in pair is bad
|
|
bad_pairs = []
|
|
good_pairs = []
|
|
for i, base in enumerate(hbo_names):
|
|
hbo_ch = raw.ch_names[hbo_picks[i]]
|
|
hbr_ch = raw.ch_names[hbr_picks[i]]
|
|
if (hbo_ch in raw.info['bads']) or (hbr_ch in raw.info['bads']):
|
|
bad_pairs.append(i)
|
|
else:
|
|
good_pairs.append(i)
|
|
|
|
print(f"Total pairs: {len(hbo_names)}")
|
|
print(f"Good pairs: {len(good_pairs)}")
|
|
print(f"Bad pairs to interpolate: {len(bad_pairs)}")
|
|
|
|
if len(bad_pairs) == 0:
|
|
print("No bad pairs found. Skipping interpolation.")
|
|
return raw, None, None
|
|
|
|
raw_before_data = raw.get_data().copy()
|
|
|
|
# Extract locations (use HbO channel loc as pair location)
|
|
locs = np.array([raw.info['chs'][hbo_picks[i]]['loc'][:3] for i in range(len(hbo_names))])
|
|
good_locs = locs[good_pairs]
|
|
bad_locs = locs[bad_pairs]
|
|
|
|
# Compute distance matrix between bad and good pairs
|
|
dist_matrix = cdist(bad_locs, good_locs)
|
|
|
|
interpolated_pairs = []
|
|
|
|
for i, bad_idx in enumerate(bad_pairs):
|
|
bad_base = hbo_names[bad_idx]
|
|
distances = dist_matrix[i]
|
|
close_idxs = np.where(distances < max_dist)[0]
|
|
|
|
print(f"\nInterpolating pair {bad_base} (index {bad_idx})")
|
|
print(f" Nearby good pairs found: {len(close_idxs)}")
|
|
|
|
if len(close_idxs) < min_neighbors:
|
|
print(f" Skipping {bad_base}: not enough neighbors (found {len(close_idxs)} < {min_neighbors})")
|
|
continue
|
|
|
|
weights = 1 / (distances[close_idxs] + 1e-6)
|
|
weights /= weights.sum()
|
|
|
|
neighbor_hbo_indices = [hbo_picks[good_pairs[idx]] for idx in close_idxs]
|
|
neighbor_hbr_indices = [hbr_picks[good_pairs[idx]] for idx in close_idxs]
|
|
|
|
neighbor_hbo_data = raw._data[neighbor_hbo_indices, :]
|
|
neighbor_hbr_data = raw._data[neighbor_hbr_indices, :]
|
|
|
|
interpolated_hbo = np.average(neighbor_hbo_data, axis=0, weights=weights)
|
|
interpolated_hbr = np.average(neighbor_hbr_data, axis=0, weights=weights)
|
|
|
|
raw._data[hbo_picks[bad_idx]] = interpolated_hbo
|
|
raw._data[hbr_picks[bad_idx]] = interpolated_hbr
|
|
|
|
interpolated_pairs.append(bad_base)
|
|
|
|
n_bad = len(bad_pairs)
|
|
n_cols = 4 # Fixed width for horizontal scaling
|
|
n_rows = int(np.ceil(n_bad / n_cols))
|
|
|
|
# Calculate height: 2.5 inches per row is usually enough for readability
|
|
fig_height = max(4, n_rows * 2.5)
|
|
fig_compare, axes = plt.subplots(n_rows, n_cols, figsize=(15, fig_height),
|
|
constrained_layout=True)
|
|
if n_bad == 1: axes = [axes] # Handle single subplot case
|
|
|
|
axes_flat = np.asarray(axes).get_data().flatten() if hasattr(axes, 'get_data') else np.asarray(axes).ravel()
|
|
for j in range(n_bad, len(axes_flat)):
|
|
if j >= n_bad:
|
|
axes_flat[j].axis('off')
|
|
|
|
times = raw.times
|
|
for i, bad_idx in enumerate(bad_pairs):
|
|
ax = axes_flat[i]
|
|
base = hbo_names[bad_idx]
|
|
|
|
# Plot "Before" (Dirty data) in light gray/red
|
|
ax.plot(times, raw_before_data[hbo_picks[bad_idx]], color='red', alpha=0.3, label='Original HbO')
|
|
|
|
# Plot "After" (Interpolated data) in solid blue/green
|
|
if base in interpolated_pairs:
|
|
ax.plot(times, raw._data[hbo_picks[bad_idx]], color='blue', label='Interpolated HbO')
|
|
status = "SUCCESS"
|
|
color = "green"
|
|
else:
|
|
status = "FAILED (Isolated)"
|
|
color = "red"
|
|
|
|
ax.set_title(f"Channel {base} | Status: {status}", color=color, fontweight='bold')
|
|
ax.set_ylabel("Amplitude")
|
|
if i == 0: ax.legend(loc='upper right')
|
|
|
|
plt.close(fig_compare)
|
|
|
|
if interpolated_pairs:
|
|
bad_ch_to_remove = []
|
|
for base_ in interpolated_pairs:
|
|
bad_ch_to_remove.append(base_ + f" {hbr_wl}") # HbR
|
|
bad_ch_to_remove.append(base_ + f" {hbo_wl}") # HbO
|
|
|
|
raw.info['bads'] = [ch for ch in raw.info['bads'] if ch not in bad_ch_to_remove]
|
|
|
|
print("\nInterpolation complete.\n")
|
|
print("Bads cleared:", raw.info['bads'])
|
|
raw.info['bads'] = []
|
|
|
|
for ch in raw.info['bads']:
|
|
print(f"Channel {ch} still marked as bad.")
|
|
|
|
fig_raw_after = raw.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="After interpolation", show=False)
|
|
|
|
return raw, fig_raw_after, fig_compare
|
|
|
|
|
|
|
|
def calculate_signal_noise_ratio(data):
|
|
"""
|
|
Calculates the signal-to-noise ratio (SNR) for each channel and identifies those below a defined threshold.
|
|
|
|
Parameters
|
|
----------
|
|
data : BaseRaw
|
|
The loaded data object to process.
|
|
|
|
Returns
|
|
-------
|
|
tuple[list[str], Figure]
|
|
- list[str]: A list of channel names that fall below the SNR threshold and are considered bad.
|
|
- Figure: A matplotlib Figure showing the channels' SNR values.
|
|
"""
|
|
|
|
print("Calculating signal to noise ratio...")
|
|
|
|
# Compute the signal-to-noise ratio values
|
|
print("Computing the signal to noise power...")
|
|
signal_band=(0.01, 0.5)
|
|
noise_band=(1.0, 10.0)
|
|
data_signal = data.copy().filter(*signal_band, verbose=False) #type: ignore
|
|
data_noise = data.copy().filter(*noise_band, verbose=False) #type: ignore
|
|
signal_power = np.mean(data_signal.get_data()**2, axis=1) #type: ignore
|
|
noise_power = np.mean(data_noise.get_data()**2, axis=1) #type: ignore
|
|
|
|
# Calculate the snr using the standard formula for dB
|
|
snr = 10 * np.log10(signal_power / (noise_power + np.finfo(float).eps))
|
|
|
|
# TODO: Understand what this does
|
|
groups: dict[str, list[str]] = {}
|
|
for ch in getattr(data, "ch_names"):
|
|
# Look for the space in the channel names and remove the characters after
|
|
# This is so we can get both oxy and deoxy to remove, as they will have the same source and detector
|
|
base = ch.rsplit(' ', 1)[0]
|
|
groups.setdefault(base, []).append(ch) # type: ignore
|
|
|
|
# If any of the channels do not meet our threshold, they will get inserted into the bad_channels set
|
|
bad_channels: set[str] = set()
|
|
for base, ch_list in groups.items():
|
|
if any(s < SNR_THRESHOLD for s, ch in zip(snr, getattr(data, "ch_names")) if ch in ch_list):
|
|
bad_channels.update(ch_list)
|
|
|
|
# Design and create the figure
|
|
print("Creating the figure...")
|
|
snr_fig, ax = plt.subplots(figsize=(12, 4), layout="constrained") # type: ignore
|
|
colors = [(0/20, 'red'), (SNR_THRESHOLD/20, 'red'), ((SNR_THRESHOLD+.5)/20, 'yellow'), ((SNR_THRESHOLD+1)/20, 'green'), (20/20, 'green')]
|
|
cmap = LinearSegmentedColormap.from_list('custom_snr_cmap', colors)
|
|
norm = mcolors.Normalize(vmin=0, vmax=20)
|
|
scatter = ax.scatter(range(len(snr)), snr, c=snr, cmap=cmap, alpha=0.8, s=100, norm=norm) # type: ignore
|
|
ax.set(xlabel="Channel Number", ylabel="Signal-to-Noise Ratio (dB)", xlim=[0, len(snr)], ylim=[0, 20])
|
|
ax.axhline(SNR_THRESHOLD, color='black', linestyle='--', alpha=0.3, linewidth=1) # type: ignore
|
|
cbar = snr_fig.colorbar(scatter, ax=ax, label="SNR Thresholds (dB)") # type: ignore
|
|
cbar.set_ticks([0, SNR_THRESHOLD, SNR_THRESHOLD+1, 20]) # type: ignore
|
|
cbar.set_ticklabels(['0', str(SNR_THRESHOLD), str(SNR_THRESHOLD+1), '20']) # type: ignore
|
|
|
|
plt.close()
|
|
|
|
print("Successfully calculated signal to noise ratio.")
|
|
|
|
return list(bad_channels), snr_fig
|
|
|
|
|
|
|
|
|
|
def calculate_peak_power(data: BaseRaw, l_freq: float = 0.7, h_freq: float = 1.5) -> tuple[list[str], Figure, Figure]:
|
|
"""
|
|
Calculate peak spectral power (PSP) for fNIRS channels and identify bad channels.
|
|
|
|
Parameters
|
|
----------
|
|
data : BaseRaw
|
|
The loaded data object to process.
|
|
l_freq : float, optional
|
|
Low cutoff frequency for filtering in Hz (default is 0.7)
|
|
h_freq : float, optional
|
|
High cutoff frequency for filtering in Hz (default is 1.5)
|
|
|
|
Returns
|
|
-------
|
|
tuple[list[str], Figure, Figure]
|
|
- list[str]: Names of channels below the PSP threshold.
|
|
- Figure: Heatmap of all PSP scores.
|
|
- Figure: Heatmap of scores above the PSP threshold.
|
|
"""
|
|
|
|
|
|
# Compute the PSP
|
|
_, scores, times = cast(tuple[NDArray[float64], NDArray[float64], list[tuple[float]]], peak_power(data, time_window=PSP_TIME_WINDOW, threshold=PSP_THRESHOLD, l_freq=l_freq, h_freq=h_freq))
|
|
|
|
# Identify channels that don't meet the provided threshold
|
|
psp = scores.mean(axis=1)
|
|
data.info["bads"] = list(compress(cast(list[str], getattr(data, "ch_names")), psp < PSP_THRESHOLD))
|
|
|
|
# Determine the colors based on the threshold, and create the figures
|
|
color_stops = ([0.0, PSP_THRESHOLD, PSP_THRESHOLD+0.1, PSP_THRESHOLD+0.2, 1.0], [0.0, PSP_THRESHOLD, PSP_THRESHOLD, 1.0])
|
|
psp1, psp2 = plot_timechannel_quality_metrics(data, scores, times, color_stops, PSP_THRESHOLD, "Peak Spectral Power")
|
|
|
|
|
|
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_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)
|
|
|
|
# Organize channels into categories
|
|
sets = [
|
|
(bad_sci, "SCI"),
|
|
(bad_psp, "PSP"),
|
|
(bad_snr, "SNR"),
|
|
(bad_coeff_var, "coeff_var"),
|
|
(bad_range, "Range"),
|
|
(bad_noise, "Noise"),
|
|
(bad_disp, "Disp.")
|
|
]
|
|
|
|
# Graph what channels were dropped and why they were dropped
|
|
channel_categories: dict[str, str] = {}
|
|
|
|
for ch in bads_combined:
|
|
present_in = [name for s, name in sets if ch in s]
|
|
# Create a label for the category
|
|
if len(present_in) == 1:
|
|
label = f"{present_in[0]} only"
|
|
else:
|
|
label = " + ".join(sorted(present_in))
|
|
channel_categories[ch] = label
|
|
|
|
# Sort channels alphabetically within categories for nicer visualization
|
|
categories = sorted(set(channel_categories.values()))
|
|
channel_names: list[str] = []
|
|
category_labels: list[str] = []
|
|
for cat in categories:
|
|
chs_in_cat = sorted([ch for ch, c in channel_categories.items() if c == cat])
|
|
channel_names.extend(chs_in_cat)
|
|
category_labels.extend([cat] * len(chs_in_cat))
|
|
|
|
colors = {cat: get_category_color(cat) for cat in categories}
|
|
# Create the figure
|
|
fig_dropped, ax = plt.subplots(figsize=(10, max(3, len(channel_names) * 0.3))) # type: ignore
|
|
y_pos = range(len(channel_names))
|
|
ax.barh(y_pos, [1]*len(channel_names), color=[colors[cat] for cat in category_labels]) # type: ignore
|
|
ax.set_yticks(y_pos) # type: ignore
|
|
ax.set_yticklabels(channel_names) # type: ignore
|
|
ax.set_xlabel("Marked as Bad") # type: ignore
|
|
ax.set_title(f"Bad Channels by Method for") # type: ignore
|
|
ax.set_xlim(0, 1)
|
|
ax.set_xticks([]) # type: ignore
|
|
ax.grid(axis='x', linestyle='--', alpha=0.7) # type: ignore
|
|
|
|
# Add a legend denoting why the channels were bad
|
|
for label, color in colors.items():
|
|
ax.bar(0, 0, color=color, label=label) # type: ignore
|
|
ax.legend() # type: ignore
|
|
|
|
fig_dropped.tight_layout()
|
|
|
|
|
|
raw_before = deepcopy(raw)
|
|
bads_channels = [ch for ch in raw.ch_names if ch in raw.info['bads']]
|
|
print(bads_channels)
|
|
if bads_channels:
|
|
fig_raw_before = raw_before.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], picks=bads_channels, title="What they were BEFORE", show=False)
|
|
else:
|
|
fig_dropped = None
|
|
fig_raw_before = None
|
|
|
|
return raw, fig_dropped, fig_raw_before, bads_channels
|
|
|
|
|
|
|
|
|
|
def filter_the_data(raw_haemo):
|
|
# --- STEP 5: Filtering (0.01-0.2 Hz bandpass) ---
|
|
fig_filter = raw_haemo.compute_psd(fmax=3).plot(
|
|
average=True, color="r", show=False, amplitude=True
|
|
)
|
|
|
|
if L_FREQ == 0 and H_FREQ != 0:
|
|
raw_haemo = raw_haemo.filter(l_freq=None, h_freq=H_FREQ, h_trans_bandwidth=H_TRANS_BANDWIDTH)
|
|
elif L_FREQ != 0 and H_FREQ == 0:
|
|
raw_haemo = raw_haemo.filter(l_freq=L_FREQ, h_freq=None, l_trans_bandwidth=L_TRANS_BANDWIDTH)
|
|
elif L_FREQ != 0 and H_FREQ != 0:
|
|
raw_haemo = raw_haemo.filter(l_freq=L_FREQ, h_freq=H_FREQ, l_trans_bandwidth=L_TRANS_BANDWIDTH, h_trans_bandwidth=H_TRANS_BANDWIDTH)
|
|
else:
|
|
print("No filter")
|
|
#raw_haemo = raw_haemo.filter(l_freq=None, h_freq=0.4, h_trans_bandwidth=0.2)
|
|
#raw_haemo = raw_haemo.filter(l_freq=None, h_freq=0.7, h_trans_bandwidth=0.2)
|
|
#raw_haemo = raw_haemo.filter(0.005, 0.7, h_trans_bandwidth=0.02, l_trans_bandwidth=0.002)
|
|
|
|
raw_haemo.compute_psd(fmax=3).plot(
|
|
average=True, axes=fig_filter.axes, color="g", amplitude=True, show=False
|
|
)
|
|
|
|
fig_raw_haemo_filter = raw_haemo.plot(duration=raw_haemo.times[-1], n_channels=raw_haemo.info['nchan'], title="Filtered HbO and HbR", show=False)
|
|
|
|
return raw_haemo, fig_filter, fig_raw_haemo_filter
|
|
|
|
|
|
|
|
def safe_create_epochs(raw, events, event_dict, tmin, tmax, baseline):
|
|
"""
|
|
Attempts to create epochs, shifting event times slightly if
|
|
sample collisions are detected.
|
|
"""
|
|
shift_increment = 1.0 / raw.info['sfreq'] # The duration of exactly one sample
|
|
#TODO: User expose this
|
|
reject_criteria = dict(hbo=80e-7)
|
|
for attempt in range(MAX_SHIFT): # Limit attempts to avoid infinite loops
|
|
try:
|
|
epochs = Epochs(
|
|
raw, events, event_id=event_dict,
|
|
tmin=tmin, tmax=tmax, baseline=baseline,
|
|
reject=reject_criteria,
|
|
preload=True, verbose=False
|
|
)
|
|
return epochs
|
|
except RuntimeError as e:
|
|
if "Event time samples were not unique" in str(e):
|
|
# Find duplicates in the events array (column 0 is the sample index)
|
|
vals, counts = np.unique(events[:, 0], return_counts=True)
|
|
duplicates = vals[counts > 1]
|
|
|
|
# Shift the second occurrence of every duplicate by 1 sample
|
|
for dup in duplicates:
|
|
idx = np.where(events[:, 0] == dup)[0][1:] # Get all but the first
|
|
events[idx, 0] += 1
|
|
|
|
print(f"Collision detected. Nudging events by {shift_increment:.4f}s and retrying...")
|
|
continue
|
|
else:
|
|
raise e # Raise if it's a different Runtime Error
|
|
|
|
raise RuntimeError("Could not resolve event collisions after 10 attempts.")
|
|
|
|
|
|
|
|
def epochs_calculations(raw_haemo, events, event_dict):
|
|
fig_epochs = [] # List to store figures
|
|
|
|
if EPOCH_HANDLING == 'shift':
|
|
epochs = safe_create_epochs(raw=raw_haemo, events=events, event_dict=event_dict, tmin=T_MIN, tmax=T_MAX, baseline=(None, 0))
|
|
else:
|
|
epochs = Epochs(raw_haemo, events, event_id=event_dict, tmin=T_MIN, tmax=T_MAX, baseline=(None, 0))
|
|
|
|
# Make a copy of the epochs and drop bad ones
|
|
epochs2 = epochs.copy()
|
|
epochs2.drop_bad()
|
|
|
|
# Plot drop log
|
|
# TODO: Why show this if we never use epochs2?
|
|
fig_epochs_dropped = epochs2.plot_drop_log(show=False)
|
|
fig_epochs.append(("fig_epochs_dropped", fig_epochs_dropped))
|
|
|
|
# Plot for each condition
|
|
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 4))
|
|
for idx, condition in enumerate(epochs.event_id.keys()):
|
|
logger.info(condition)
|
|
logger.info(idx)
|
|
# Plot images for each condition
|
|
fig_epochs_data = epochs[condition].plot_image(
|
|
combine="mean",
|
|
vmin=-1,
|
|
vmax=1,
|
|
ts_args=dict(ylim=dict(hbo=[-1, 1], hbr=[-1, 1])),
|
|
show=False
|
|
)
|
|
for j, fig in enumerate(fig_epochs_data):
|
|
logger.info("------------------------------------------")
|
|
logger.info(j)
|
|
logger.info(fig)
|
|
|
|
ax = fig.axes[0]
|
|
original_title = ax.get_title()
|
|
ax.set_title(f"{condition}: {original_title}")
|
|
fig_epochs.append((f"fig_{condition}_data_{idx}_{j}", fig)) # Store with a unique name
|
|
|
|
# Evoked average figure for each condition
|
|
evoked_avg = epochs[condition].average()
|
|
clims = dict(hbo=[-1, 1], hbr=[1, -1])
|
|
condition_fig = evoked_avg.plot_image(clim=clims, show=False)
|
|
|
|
for ax in condition_fig.axes:
|
|
original_title = ax.get_title()
|
|
ax.set_title(f"{original_title} - {condition}")
|
|
fig_epochs.append((f"evoked_avg_{condition}", condition_fig)) # Store with a unique name
|
|
|
|
# Prepare evokeds and colors for topographic plot
|
|
evokeds3 = []
|
|
colors = []
|
|
conditions = list(epochs.event_id.keys())
|
|
cmap = plt.get_cmap("tab10", len(conditions))
|
|
|
|
for idx, cond in enumerate(conditions):
|
|
evoked = epochs[cond].average(picks="hbo")
|
|
evokeds3.append(evoked)
|
|
colors.append(cmap(idx))
|
|
|
|
# Create the topographic plot
|
|
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 4))
|
|
help = plot_evoked_topo(evokeds3, color=colors, axes=axes, legend=False, show=False)
|
|
|
|
# Build custom legend
|
|
lines = []
|
|
for color in colors:
|
|
line = plt.Line2D([0], [0], color=color, lw=2)
|
|
lines.append(line)
|
|
|
|
fig.legend(lines, conditions, loc="lower right")
|
|
fig_epochs.append(("evoked_topo", help)) # Store with a unique name
|
|
|
|
unique_annotations = set(raw_haemo.annotations.description)
|
|
|
|
for cond in unique_annotations:
|
|
|
|
# Evoked response for specific condition ("Activity")
|
|
evoked_stim1 = epochs[cond].average()
|
|
|
|
fig_evoked_hbo = evoked_stim1.copy().pick(picks='hbo').plot(time_unit='s', show=False)
|
|
fig_evoked_hbr = evoked_stim1.copy().pick(picks='hbr').plot(time_unit='s', show=False)
|
|
|
|
fig_epochs.append((f"fig_evoked_hbo_{cond}", fig_evoked_hbo)) # Store with a unique name
|
|
fig_epochs.append((f"fig_evoked_hbr_{cond}", fig_evoked_hbr)) # Store with a unique name
|
|
|
|
print("Evoked HbO peak amplitude:", evoked_stim1.copy().pick(picks='hbo').data.max())
|
|
|
|
evokeds = {}
|
|
for condition in epochs2.event_id:
|
|
evokeds[condition] = epochs2[condition].average()
|
|
print(f"Condition '{condition}': {len(epochs2[condition])} epochs averaged.")
|
|
|
|
all_evokeds = {}
|
|
for condition in epochs.event_id:
|
|
if condition not in all_evokeds:
|
|
all_evokeds[condition] = []
|
|
all_evokeds[condition].append(epochs[condition].average())
|
|
|
|
group_aucs = {}
|
|
# TODO: group averages with a single person?
|
|
group_averages = {cond: grand_average(evokeds) for cond, evokeds in all_evokeds.items()}
|
|
for condition, evoked in group_averages.items():
|
|
group_aucs[condition] = {}
|
|
for pick in ["hbo", "hbr"]:
|
|
picks_idx = [i for i, ch in enumerate(evoked.ch_names) if pick in ch]
|
|
if not picks_idx:
|
|
continue
|
|
|
|
data = evoked.data[picks_idx, :].mean(axis=0)
|
|
t_start, t_end = 0, 15 #TODO: Is this in seconds? or is it 1hz input that makes it 15s?
|
|
times_mask = (evoked.times >= t_start) & (evoked.times <= t_end)
|
|
data_segment = data[times_mask]
|
|
times_segment = evoked.times[times_mask]
|
|
|
|
auc = np.trapezoid(data_segment, times_segment)
|
|
group_aucs[condition][pick] = auc
|
|
|
|
# Final evoked comparison plot for each condition
|
|
for condition in conditions:
|
|
if condition not in evokeds:
|
|
continue
|
|
evoked = evokeds[condition]
|
|
|
|
fig, ax = plt.subplots(figsize=(6, 5))
|
|
legend_labels = ["Oxyhaemoglobin"]
|
|
|
|
for pick, color in zip(["hbo", "hbr"], ["r", "b"]):
|
|
plot_compare_evokeds(
|
|
evoked,
|
|
combine="mean",
|
|
picks=pick,
|
|
axes=ax,
|
|
show=False,
|
|
colors=[color],
|
|
legend=False,
|
|
title=f"Participant: nCondition: {condition}",
|
|
ylim=dict(hbo=[-0.5, 1], hbr=[-0.5, 1]),
|
|
show_sensors=False,
|
|
)
|
|
auc_value = group_aucs.get(condition, {}).get(pick, None)
|
|
if auc_value is not None:
|
|
label = f"{pick.upper()} AUC: {auc_value * 1e6:.4f} µM·s"
|
|
else:
|
|
label = f"{pick.upper()} AUC: N/A"
|
|
legend_labels.append(label)
|
|
if len(legend_labels) == 2:
|
|
legend_labels.append("Deoxyhaemoglobin")
|
|
|
|
ax.legend(legend_labels)
|
|
|
|
fig_epochs.append((f"fig_{condition}_compare_evokeds", fig)) # Store with a unique name
|
|
|
|
return epochs, fig_epochs
|
|
|
|
|
|
|
|
def make_design_matrix(raw_haemo):
|
|
|
|
events_to_remove = REMOVE_EVENTS
|
|
|
|
filtered_annotations = [ann for ann in raw_haemo.annotations if ann['description'] not in events_to_remove]
|
|
|
|
new_annot = Annotations(
|
|
onset=[ann['onset'] for ann in filtered_annotations],
|
|
duration=[ann['duration'] for ann in filtered_annotations],
|
|
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)
|
|
|
|
if RESAMPLE:
|
|
raw_haemo.resample(RESAMPLE_FREQ, npad="auto")
|
|
raw_haemo._data = raw_haemo._data * 1e6
|
|
try:
|
|
short_chans.resample(RESAMPLE_FREQ)
|
|
except:
|
|
pass
|
|
|
|
design_matrix = make_first_level_design_matrix(
|
|
raw=raw_haemo,
|
|
stim_dur=STIM_DUR,
|
|
hrf_model=HRF_MODEL,
|
|
drift_model=DRIFT_MODEL,
|
|
high_pass=HIGH_PASS,
|
|
drift_order=DRIFT_ORDER,
|
|
fir_delays=FIR_DELAYS,
|
|
min_onset=MIN_ONSET,
|
|
oversampling=OVERSAMPLING
|
|
)
|
|
|
|
# 3) Average and Append Short Channels
|
|
if SHORT_CHANNEL_REGRESSION and not FOLDING_BYP:
|
|
if short_chans is not None and len(short_chans.ch_names) > 0:
|
|
ch_types = short_chans.get_channel_types()
|
|
|
|
# Scenario A: Short channels are already converted to Hemoglobin (hbo/hbr)
|
|
if "hbo" in ch_types or "hbr" in ch_types:
|
|
hbo_chs = [ch for ch, t in zip(short_chans.ch_names, ch_types) if t == "hbo"]
|
|
hbr_chs = [ch for ch, t in zip(short_chans.ch_names, ch_types) if t == "hbr"]
|
|
|
|
if hbo_chs:
|
|
hbo_data = short_chans.copy().pick(hbo_chs).get_data()
|
|
design_matrix["ShortHbO"] = np.mean(hbo_data, axis=0)
|
|
if hbr_chs:
|
|
hbr_data = short_chans.copy().pick(hbr_chs).get_data()
|
|
design_matrix["ShortHbR"] = np.mean(hbr_data, axis=0)
|
|
print(f"Successfully added averaged ShortHbO ({len(hbo_chs)} chs) and ShortHbR ({len(hbr_chs)} chs) to the matrix.")
|
|
|
|
# Scenario B: Short channels are raw wavelengths (760nm, 850nm, etc.)
|
|
else:
|
|
wavelength_groups = {}
|
|
for ch_name in short_chans.ch_names:
|
|
# Look for the wavelength number (digits) at the end of the channel name
|
|
match = re.search(r'(\d+)$', ch_name)
|
|
if match:
|
|
wl = match.group(1)
|
|
wavelength_groups.setdefault(wl, []).append(ch_name)
|
|
|
|
if wavelength_groups:
|
|
for wl, chs in wavelength_groups.items():
|
|
wl_data = short_chans.copy().pick(chs).get_data()
|
|
col_name = f"Short_{wl}"
|
|
design_matrix[col_name] = np.mean(wl_data, axis=0)
|
|
print(f"Successfully added averaged short channels by wavelength: {list(wavelength_groups.keys())}")
|
|
else:
|
|
# Emergency fallback: if names have no digits, average all of them together
|
|
design_matrix["Short_Avg"] = np.mean(short_chans.get_data(), axis=0)
|
|
print("Could not detect wavelengths. Averaged all short channels into 'Short_Avg'.")
|
|
else:
|
|
print("Warning: SHORT_CHANNEL_REGRESSION is True, but no short channels were found.")
|
|
|
|
print(design_matrix.head())
|
|
print(design_matrix.columns)
|
|
|
|
fig, ax1 = plt.subplots(figsize=(10, 6), constrained_layout=True)
|
|
_ = plot_design_matrix(design_matrix, axes=ax1)
|
|
|
|
return raw_haemo, design_matrix, fig
|
|
|
|
|
|
|
|
def generate_montage_locations():
|
|
"""Get standard MNI montage locations in dataframe.
|
|
|
|
Data is returned in the same format as the eeg_positions library.
|
|
"""
|
|
# standard_1020 and standard_1005 are in MNI (fsaverage) space already,
|
|
# but we need to undo the scaling that head_scale will do
|
|
montage = make_standard_montage(
|
|
"standard_1005", head_size=0.09700884729534559
|
|
)
|
|
for d in montage.dig:
|
|
d["coord_frame"] = 2003
|
|
montage.dig[:] = montage.dig[3:]
|
|
montage.add_mni_fiducials() # now in fsaverage space
|
|
coords = pd.DataFrame.from_dict(montage.get_positions()["ch_pos"]).T
|
|
coords["label"] = coords.index
|
|
coords = coords.rename(columns={0: "x", 1: "y", 2: "z"})
|
|
|
|
return coords.reset_index(drop=True)
|
|
|
|
|
|
def _find_closest_standard_location(position, reference, *, out="label"):
|
|
"""Return closest montage label to coordinates.
|
|
|
|
Parameters
|
|
----------
|
|
position : array, shape (3,)
|
|
Coordinates.
|
|
reference : dataframe
|
|
As generated by _generate_montage_locations.
|
|
trans_pos : str
|
|
Apply a transformation to positions to specified frame.
|
|
Use None for no transformation.
|
|
"""
|
|
|
|
p0 = np.array(position)
|
|
p0.shape = (-1, 3)
|
|
# head_mri_t, _ = _get_trans("fsaverage", "head", "mri")
|
|
# p0 = apply_trans(head_mri_t, p0)
|
|
dists = cdist(p0, np.asarray(reference[["x", "y", "z"]], float))
|
|
|
|
if out == "label":
|
|
min_idx = np.argmin(dists)
|
|
return reference["label"][min_idx]
|
|
else:
|
|
assert out == "dists"
|
|
return dists
|
|
|
|
|
|
|
|
def _source_detector_fold_table(raw, cidx, reference, fold_tbl, interpolate):
|
|
src = raw.info["chs"][cidx]["loc"][3:6]
|
|
det = raw.info["chs"][cidx]["loc"][6:9]
|
|
|
|
ref_lab = list(reference["label"])
|
|
dists = _find_closest_standard_location([src, det], reference, out="dists")
|
|
src_min, det_min = np.argmin(dists, axis=1)
|
|
src_name, det_name = ref_lab[src_min], ref_lab[det_min]
|
|
|
|
tbl = fold_tbl.query("Source == @src_name and Detector == @det_name")
|
|
dist = np.linalg.norm(dists[[0, 1], [src_min, det_min]])
|
|
# Try reversing source and detector
|
|
if len(tbl) == 0:
|
|
tbl = fold_tbl.query("Source == @det_name and Detector == @src_name")
|
|
if len(tbl) == 0 and interpolate:
|
|
# Try something hopefully not too terrible: pick the one with the
|
|
# smallest net distance
|
|
good = np.isin(fold_tbl["Source"], reference["label"]) & np.isin(
|
|
fold_tbl["Detector"], reference["label"]
|
|
)
|
|
assert good.any()
|
|
tbl = fold_tbl[good]
|
|
assert len(tbl)
|
|
src_idx = [ref_lab.index(src) for src in tbl["Source"]]
|
|
det_idx = [ref_lab.index(det) for det in tbl["Detector"]]
|
|
# Original
|
|
tot_dist = np.linalg.norm([dists[0, src_idx], dists[1, det_idx]], axis=0)
|
|
assert tot_dist.shape == (len(tbl),)
|
|
idx = np.argmin(tot_dist)
|
|
dist_1 = tot_dist[idx]
|
|
src_1, det_1 = ref_lab[src_idx[idx]], ref_lab[det_idx[idx]]
|
|
# And the reverse
|
|
tot_dist = np.linalg.norm([dists[0, det_idx], dists[1, src_idx]], axis=0)
|
|
idx = np.argmin(tot_dist)
|
|
dist_2 = tot_dist[idx]
|
|
src_2, det_2 = ref_lab[det_idx[idx]], ref_lab[src_idx[idx]]
|
|
if dist_1 < dist_2:
|
|
new_dist, src_use, det_use = dist_1, src_1, det_1
|
|
else:
|
|
new_dist, src_use, det_use = dist_2, det_2, src_2
|
|
|
|
|
|
tbl = fold_tbl.query("Source == @src_use and Detector == @det_use")
|
|
tbl = tbl.copy()
|
|
tbl["BestSource"] = src_name
|
|
tbl["BestDetector"] = det_name
|
|
tbl["BestMatchDistance"] = dist
|
|
tbl["MatchDistance"] = new_dist
|
|
assert len(tbl)
|
|
else:
|
|
tbl = tbl.copy()
|
|
tbl["BestSource"] = src_name
|
|
tbl["BestDetector"] = det_name
|
|
tbl["BestMatchDistance"] = dist
|
|
tbl["MatchDistance"] = dist
|
|
|
|
tbl = tbl.copy() # don't get warnings about setting values later
|
|
return tbl
|
|
|
|
|
|
|
|
def _read_fold_xls(fname, atlas="Juelich"):
|
|
"""Read fOLD toolbox xls file.
|
|
|
|
The values are then manipulated in to a tidy dataframe.
|
|
|
|
Note the xls files are not included as no license is provided.
|
|
|
|
Parameters
|
|
----------
|
|
fname : str
|
|
Path to xls file.
|
|
atlas : str
|
|
Requested atlas.
|
|
"""
|
|
page_reference = {"AAL2": 2, "AICHA": 5, "Brodmann": 8, "Juelich": 11, "Loni": 14}
|
|
|
|
tbl = pd.read_excel(fname, sheet_name=page_reference[atlas])
|
|
|
|
# Remove the spacing between rows
|
|
empty_rows = np.where(np.isnan(tbl["Specificity"]))[0]
|
|
tbl = tbl.drop(empty_rows).reset_index(drop=True)
|
|
|
|
# Empty values in the table mean its the same as above
|
|
for row_idx in range(1, tbl.shape[0]):
|
|
for col_idx, col in enumerate(tbl.columns):
|
|
if not isinstance(tbl[col][row_idx], str):
|
|
if np.isnan(tbl[col][row_idx]):
|
|
tbl.iloc[row_idx, col_idx] = tbl.iloc[row_idx - 1, col_idx]
|
|
|
|
tbl["Specificity"] = tbl["Specificity"] * 100
|
|
tbl["brainSens"] = tbl["brainSens"] * 100
|
|
return tbl
|
|
|
|
|
|
def _check_load_fold(fold_files, atlas):
|
|
# _validate_type(fold_files, (list, "path-like", None), "fold_files")
|
|
if fold_files is None:
|
|
fold_files = get_config("MNE_NIRS_FOLD_PATH")
|
|
if fold_files is None:
|
|
raise ValueError(
|
|
"MNE_NIRS_FOLD_PATH not set, either set it using "
|
|
"mne.set_config or pass fold_files as str or list"
|
|
)
|
|
if not isinstance(fold_files, list): # path-like
|
|
fold_files = _check_fname(
|
|
fold_files,
|
|
overwrite="read",
|
|
must_exist=True,
|
|
name="fold_files",
|
|
need_dir=True,
|
|
)
|
|
fold_files = [op.join(fold_files, f"10-{x}.xls") for x in (5, 10)]
|
|
|
|
fold_tbl = pd.DataFrame()
|
|
for fi, fname in enumerate(fold_files):
|
|
fname = _check_fname(
|
|
fname, overwrite="read", must_exist=True, name=f"fold_files[{fi}]"
|
|
)
|
|
fold_tbl = pd.concat(
|
|
[fold_tbl, _read_fold_xls(fname, atlas=atlas)], ignore_index=True
|
|
)
|
|
return fold_tbl
|
|
|
|
|
|
|
|
def fold_channel_specificity_normal(raw, fold_files=None, atlas="Juelich", interpolate=False):
|
|
"""Return the landmarks and specificity a channel is sensitive to.
|
|
|
|
Parameters
|
|
|
|
""" # noqa: E501
|
|
_validate_type(raw, BaseRaw, "raw")
|
|
|
|
reference_locations = generate_montage_locations()
|
|
|
|
fold_tbl = _check_load_fold(fold_files, atlas)
|
|
|
|
chan_spec = list()
|
|
for cidx in range(len(raw.ch_names)):
|
|
tbl = _source_detector_fold_table(
|
|
raw, cidx, reference_locations, fold_tbl, interpolate
|
|
)
|
|
chan_spec.append(tbl.reset_index(drop=True))
|
|
|
|
return chan_spec
|
|
|
|
|
|
|
|
def resource_path(relative_path):
|
|
"""
|
|
Get absolute path to resource regardless of running directly or packaged using PyInstaller
|
|
"""
|
|
|
|
if hasattr(sys, '_MEIPASS'):
|
|
# PyInstaller bundle path
|
|
base_path = sys._MEIPASS
|
|
else:
|
|
base_path = os.path.abspath(".")
|
|
|
|
return os.path.join(base_path, relative_path)
|
|
|
|
|
|
|
|
# def fold_channels(raw: BaseRaw) -> None:
|
|
|
|
# # Locate the fOLD excel files
|
|
# if getattr(sys, 'frozen', False):
|
|
# set_config('MNE_NIRS_FOLD_PATH', resource_path("./mne_data/fOLD/fOLD-public-master/Supplementary")) # type: ignore
|
|
# else:
|
|
# path = os.path.expanduser("~") + "/mne_data/fOLD/fOLD-public-master/Supplementary"
|
|
# set_config('MNE_NIRS_FOLD_PATH', resource_path(path)) # type: ignore
|
|
|
|
# output = None
|
|
|
|
# # List to store the results
|
|
# landmark_specificity_data: list[dict[str, Any]] = []
|
|
|
|
# # Filter the data to only what we want
|
|
# hbo_channel_names = cast(list[str], getattr(raw.copy().pick(picks='hbo'), "ch_names")) # type: ignore
|
|
|
|
# # Format the output to make it slightly easier to read
|
|
|
|
# if True:
|
|
# num_channels = len(hbo_channel_names)
|
|
# rows, cols = 4, 7 # 6 rows and 4 columns of pie charts
|
|
# fig, axes = plt.subplots(rows, cols, figsize=(16, 10), constrained_layout=True)
|
|
# axes = axes.flatten() # Flatten the axes array for easier indexing
|
|
|
|
# # If more pie charts than subplots, create extra subplots
|
|
# if num_channels > rows * cols:
|
|
# fig, axes = plt.subplots((num_channels // cols) + 1, cols, figsize=(16, 10), constrained_layout=True)
|
|
# axes = axes.flatten()
|
|
|
|
# # Create a list for consistent color mapping
|
|
# landmarks = [
|
|
# "1 - Primary Somatosensory Cortex",
|
|
# "2 - Primary Somatosensory Cortex",
|
|
# "3 - Primary Somatosensory Cortex",
|
|
# "4 - Primary Motor Cortex",
|
|
# "5 - Somatosensory Association Cortex",
|
|
# "6 - Pre-Motor and Supplementary Motor Cortex",
|
|
# "7 - Somatosensory Association Cortex",
|
|
# "8 - Includes Frontal eye fields",
|
|
# "9 - Dorsolateral prefrontal cortex",
|
|
# "10 - Frontopolar area",
|
|
# "11 - Orbitofrontal area",
|
|
# "17 - Primary Visual Cortex (V1)",
|
|
# "18 - Visual Association Cortex (V2)",
|
|
# "19 - V3",
|
|
# "20 - Inferior Temporal gyrus",
|
|
# "21 - Middle Temporal gyrus",
|
|
# "22 - Superior Temporal Gyrus",
|
|
# "23 - Ventral Posterior cingulate cortex",
|
|
# "24 - Ventral Anterior cingulate cortex",
|
|
# "25 - Subgenual cortex",
|
|
# "32 - Dorsal anterior cingulate cortex",
|
|
# "37 - Fusiform gyrus",
|
|
# "38 - Temporopolar area",
|
|
# "39 - Angular gyrus, part of Wernicke's area",
|
|
# "40 - Supramarginal gyrus part of Wernicke's area",
|
|
# "41 - Primary and Auditory Association Cortex",
|
|
# "42 - Primary and Auditory Association Cortex",
|
|
# "43 - Subcentral area",
|
|
# "44 - pars opercularis, part of Broca's area",
|
|
# "45 - pars triangularis Broca's area",
|
|
# "46 - Dorsolateral prefrontal cortex",
|
|
# "47 - Inferior prefrontal gyrus",
|
|
# "48 - Retrosubicular area",
|
|
# "Brain_Outside",
|
|
# ]
|
|
|
|
# cmap1 = plt.get_cmap('tab20') # First 20 colors
|
|
# cmap2 = plt.get_cmap('tab20b') # Next 20 colors
|
|
|
|
# # Combine the colors from both colormaps
|
|
# colors = [cmap1(i) for i in range(20)] + [cmap2(i) for i in range(20)] # Total 40 colors
|
|
|
|
# landmarks.sort(key=lambda x: (int(x.split(" - ")[0]) if x.split(" - ")[0].isdigit() else float('inf')))
|
|
|
|
# landmark_color_map = {landmark: colors[i % len(colors)] for i, landmark in enumerate(landmarks)}
|
|
|
|
# # Iterate over each channel
|
|
# print(len(hbo_channel_names))
|
|
|
|
# for idx, channel_name in enumerate(hbo_channel_names):
|
|
|
|
# print(idx, channel_name)
|
|
# # Run the fOLD on the selected channel
|
|
# channel_data = raw.copy().pick(picks=channel_name) # type: ignore
|
|
|
|
# output = cast(list[DataFrame], fold_channel_specificity_normal(channel_data, interpolate=True, atlas='Brodmann'))
|
|
|
|
# # Process each DataFrame that fold_channel_specificity returns
|
|
# for df_data in output:
|
|
|
|
# # Extract the relevant columns
|
|
# useful_data = df_data[['Landmark', 'Specificity']]
|
|
|
|
# # Store the results
|
|
# landmark_specificity_data.append({
|
|
# 'Channel': channel_name,
|
|
# 'Data': useful_data,
|
|
# })
|
|
|
|
|
|
# # Plot the results
|
|
# # TODO: Fix this
|
|
# if True:
|
|
# unique_landmarks = sorted(useful_data['Landmark'].unique())
|
|
# color_list = [landmark_color_map[landmark] for landmark in useful_data['Landmark']]
|
|
|
|
# # Plot specificity for each channel
|
|
# ax = axes[idx]
|
|
|
|
# labels = [f'{landmark.split(" - ")[0]}' if landmark != 'Brain_Outside' else 'B' for landmark in useful_data['Landmark']]
|
|
|
|
# wedges, texts, autotexts = ax.pie(
|
|
# useful_data['Specificity'],
|
|
# autopct='%1.1f%%',
|
|
# startangle=90,
|
|
# labels=labels,
|
|
# labeldistance=1.05,
|
|
# colors=color_list)
|
|
|
|
# ax.set_title(f'{channel_name}')
|
|
# ax.axis('equal')
|
|
|
|
# landmark_specificity_data = []
|
|
|
|
|
|
# # TODO: Fix this
|
|
# if True:
|
|
# handles = [
|
|
# plt.Line2D([0], [0], marker='o', color='w', label=landmark, markersize=10,
|
|
# markerfacecolor=landmark_color_map[landmark])
|
|
# for landmark in landmarks
|
|
# ]
|
|
# n_landmarks = len(landmarks)
|
|
|
|
# # Calculate the figure size based on number of rows and columns
|
|
# fig_width = 5
|
|
# fig_height = n_landmarks / 4
|
|
|
|
# # Create a new figure window for the legend
|
|
# legend_fig = plt.figure(figsize=(fig_width, fig_height))
|
|
# legend_axes = legend_fig.add_subplot(111)
|
|
# legend_axes.axis('off') # Turn off axis for the legend window
|
|
# legend_axes.legend(handles=handles, loc='center', fontsize=10, title="Landmarks")
|
|
|
|
# for ax in axes[len(hbo_channel_names):]:
|
|
# ax.axis('off')
|
|
|
|
# #plt.show()
|
|
# fig_dict = {"main": fig, "legend": legend_fig}
|
|
# return convert_fig_dict_to_png_bytes(fig_dict)
|
|
|
|
|
|
|
|
def fold_channels(raw: BaseRaw, p_name: str, progress_queue=None) -> dict[str, list[dict[str, Any]]]:
|
|
"""Runs in background process.
|
|
Does only heavy math/lookup. Returns data instead of a static image.
|
|
"""
|
|
if getattr(sys, 'frozen', False):
|
|
set_config('MNE_NIRS_FOLD_PATH', resource_path("./mne_data/fOLD/fOLD-public-master/Supplementary"))
|
|
else:
|
|
path = os.path.expanduser("~") + "/mne_data/fOLD/fOLD-public-master/Supplementary"
|
|
set_config('MNE_NIRS_FOLD_PATH', resource_path(path))
|
|
|
|
hbo_channel_names = cast(list[str], getattr(raw.copy().pick(picks='hbo'), "ch_names"))
|
|
|
|
# Store clean, picklable data lists instead of complex DataFrames
|
|
channel_results = {}
|
|
|
|
step_idx = 0
|
|
|
|
for channel_name in hbo_channel_names:
|
|
channel_data = raw.copy().pick(picks=channel_name)
|
|
output = cast(list[DataFrame], fold_channel_specificity_normal(channel_data, interpolate=True, atlas='Brodmann'))
|
|
|
|
channel_results[channel_name] = []
|
|
for df_data in output:
|
|
# Extract just raw primitive types so they transfer over process channels flawlessly
|
|
for _, row in df_data.iterrows():
|
|
channel_results[channel_name].append({
|
|
'Landmark': str(row['Landmark']),
|
|
'Specificity': float(row['Specificity'])
|
|
})
|
|
step_idx += 1
|
|
if progress_queue is not None:
|
|
progress_queue.put((p_name, step_idx))
|
|
|
|
# Return raw data dictionary to the result_queue
|
|
return channel_results
|
|
|
|
|
|
|
|
|
|
def plot_glm_results(file_path, raw_haemo, glm_est, design_matrix):
|
|
|
|
fig_glms = [] # List to store figures
|
|
|
|
dm = design_matrix.copy()
|
|
logger.info(design_matrix.shape)
|
|
logger.info(design_matrix.columns)
|
|
logger.info(design_matrix.head())
|
|
|
|
rois = dict(AllChannels=range(len(raw_haemo.ch_names)))
|
|
conditions = design_matrix.columns
|
|
df_individual = glm_est.to_dataframe_region_of_interest(rois, conditions)
|
|
|
|
df_individual["ID"] = file_path
|
|
# df_individual["theta"] = [t * 1.0e6 for t in df_individual["theta"]]
|
|
|
|
first_onset_for_cond = {}
|
|
for onset, desc in zip(raw_haemo.annotations.onset, raw_haemo.annotations.description):
|
|
if desc not in first_onset_for_cond:
|
|
first_onset_for_cond[desc] = onset
|
|
|
|
# Get unique condition names from annotations (descriptions)
|
|
unique_annotations = set(raw_haemo.annotations.description)
|
|
|
|
for cond in unique_annotations:
|
|
logger.info(cond)
|
|
df_individual_filtered = df_individual.copy()
|
|
|
|
# Filter for the condition of interest and FIR delays
|
|
df_individual_filtered["isCondition"] = [cond in n for n in df_individual_filtered["Condition"]]
|
|
df_individual_filtered["isDelay"] = ["delay" in n for n in df_individual_filtered["Condition"]]
|
|
df_individual_filtered = df_individual_filtered.query("isDelay and isCondition")
|
|
|
|
# Remove other conditions from design matrix
|
|
dm_condition_cols = [col for col in dm.columns if cond in col]
|
|
dm_cond = dm[dm_condition_cols]
|
|
|
|
# Add a numeric delay column
|
|
def extract_delay_number(condition_str):
|
|
# Extracts the number at the end of a string like 'Activity_delay_5'
|
|
return int(condition_str.split("_")[-1])
|
|
|
|
df_individual_filtered["DelayNum"] = df_individual_filtered["Condition"].apply(extract_delay_number)
|
|
|
|
# Now separate and sort using numeric delay
|
|
df_hbo = df_individual_filtered[df_individual_filtered["Chroma"] == "hbo"].sort_values("DelayNum")
|
|
df_hbr = df_individual_filtered[df_individual_filtered["Chroma"] == "hbr"].sort_values("DelayNum")
|
|
|
|
vals_hbo = df_hbo["theta"].values
|
|
vals_hbr = df_hbr["theta"].values
|
|
|
|
# Create the plot
|
|
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(19, 10))
|
|
|
|
# Scale design matrix components using numpy arrays instead of pandas operations
|
|
dm_cond_values = dm_cond.values
|
|
dm_cond_scaled_hbo = dm_cond_values * vals_hbo.reshape(1, -1)
|
|
dm_cond_scaled_hbr = dm_cond_values * vals_hbr.reshape(1, -1)
|
|
|
|
# Create time axis relative to stimulus onset
|
|
time = dm_cond.index - np.ceil(first_onset_for_cond.get(cond, 0))
|
|
|
|
# Plot
|
|
axes[0].plot(time, dm_cond_values)
|
|
axes[1].plot(time, dm_cond_scaled_hbo)
|
|
axes[2].plot(time, np.sum(dm_cond_scaled_hbo, axis=1), 'r')
|
|
axes[2].plot(time, np.sum(dm_cond_scaled_hbr, axis=1), 'b')
|
|
|
|
# Format plots
|
|
for ax in range(3):
|
|
axes[ax].set_xlim(-5, 25)
|
|
axes[ax].set_xlabel("Time (s)")
|
|
axes[0].set_ylim(-0.2, 1.2)
|
|
axes[1].set_ylim(-0.5, 1)
|
|
axes[2].set_ylim(-0.5, 1)
|
|
axes[0].set_title(f"FIR Model (Unscaled)")
|
|
axes[1].set_title(f"FIR Components (Scaled by {cond} GLM Estimates)")
|
|
axes[2].set_title(f"Evoked Response ({cond})")
|
|
axes[0].set_ylabel("FIR Model")
|
|
axes[1].set_ylabel("Oxyhaemoglobin (ΔμMol)")
|
|
axes[2].set_ylabel("Haemoglobin (ΔμMol)")
|
|
axes[2].legend(["Oxyhaemoglobin", "Deoxyhaemoglobin"])
|
|
|
|
|
|
print(f"Number of FIR bins: {len(vals_hbo)}")
|
|
print(f"Mean theta (HbO): {np.mean(vals_hbo):.4f}")
|
|
print(f"Sum of theta (HbO): {np.sum(vals_hbo):.4f}")
|
|
print(f"Mean theta (HbR): {np.mean(vals_hbr):.4f}")
|
|
print(f"Sum of theta (HbR): {np.sum(vals_hbr):.4f}")
|
|
|
|
fig_glms.append((f"Condition {cond}", fig))
|
|
|
|
return fig_glms
|
|
|
|
|
|
def plot_3d_evoked_array(
|
|
inst: Union[BaseRaw, EvokedArray, Info],
|
|
statsmodel_df: DataFrame,
|
|
picks: Optional[Union[str, list[str]]] = "hbo",
|
|
value: str = "Coef.",
|
|
background: str = "w",
|
|
figure: Optional[object] = None,
|
|
clim: Union[str, dict[str, Union[str, list[float]]]] = "auto",
|
|
mode: str = "weighted",
|
|
colormap: str = "RdBu_r",
|
|
surface: str = "pial",
|
|
hemi: str = "both",
|
|
size: int = 800,
|
|
view: Optional[Union[str, dict[str, float]]] = None,
|
|
colorbar: bool = True,
|
|
distance: float = 0.03,
|
|
subjects_dir: Optional[str] = None,
|
|
src: Optional[SourceSpaces] = None,
|
|
verbose: bool = False,
|
|
) -> Brain:
|
|
'''Ported from MNE'''
|
|
|
|
info: Info = cast(Info, deepcopy(inst if isinstance(inst, Info) else inst.info)) # type: ignore
|
|
if not (getattr(info, "ch_names") == list(statsmodel_df["ch_name"].values)): # type: ignore
|
|
raise RuntimeError(
|
|
'MNE data structure does not match dataframe '
|
|
f'results.\nMNE = {getattr(info, "ch_names")}.\n'
|
|
f'GLM = {list(statsmodel_df["ch_name"].values)}' # type: ignore
|
|
)
|
|
|
|
ea = EvokedArray(np.tile(statsmodel_df[value].values.T, (1, 1)).T, info.copy()) # type: ignore
|
|
|
|
# TODO: mimic behaviour of other MNE-NIRS glm plotting options
|
|
if picks is not None:
|
|
ea = ea.pick(picks=picks) # type: ignore
|
|
|
|
if subjects_dir is None:
|
|
subjects_dir = os.environ["SUBJECTS_DIR"]
|
|
if src is None:
|
|
fname_src_fs = os.path.join(
|
|
subjects_dir, "fsaverage", "bem", "fsaverage-ico-5-src.fif"
|
|
)
|
|
src = read_source_spaces(fname_src_fs, verbose=verbose)
|
|
|
|
picks = getattr(ea, "info")["ch_names"]
|
|
|
|
# Set coord frame
|
|
for idx in range(len(getattr(ea, "ch_names"))):
|
|
getattr(ea, "info")["chs"][idx]["coord_frame"] = 4
|
|
|
|
# Generate source estimate
|
|
kwargs = dict(
|
|
evoked=ea,
|
|
subject="fsaverage",
|
|
trans=Transform('head', 'mri', np.eye(4)),
|
|
distance=distance,
|
|
mode=mode,
|
|
surface=surface,
|
|
subjects_dir=subjects_dir,
|
|
src=src,
|
|
project=True,
|
|
)
|
|
|
|
stc = stc_near_sensors(picks=picks, **kwargs, verbose=verbose) # type: ignore
|
|
|
|
assert isinstance(stc, SourceEstimate)
|
|
|
|
# Produce brain plot
|
|
brain: Brain = stc.plot( # type: ignore
|
|
src=src,
|
|
subjects_dir=subjects_dir,
|
|
hemi=hemi,
|
|
surface=surface,
|
|
initial_time=0,
|
|
clim=clim, # type: ignore
|
|
size=size,
|
|
colormap=colormap,
|
|
figure=figure,
|
|
background=background,
|
|
colorbar=colorbar,
|
|
verbose=verbose,
|
|
)
|
|
if view is not None:
|
|
brain.show_view(view) # type: ignore
|
|
|
|
return brain
|
|
|
|
|
|
def aggregate_fnirs_group_geometry(raw_list):
|
|
"""
|
|
Averages fNIRS geometry across participants in two tiers:
|
|
1. Average by Channel Pairing (S_D).
|
|
2. Average by Individual Optode (S, D) across all averaged pairings.
|
|
Returns a unified MNE Raw object with exactly one dot per optode.
|
|
"""
|
|
import mne
|
|
import numpy as np
|
|
|
|
channel_locs = {}
|
|
all_ch_names = []
|
|
|
|
for raw in raw_list:
|
|
if raw is None: continue
|
|
raw_hbo = raw.copy().pick(picks="hbo")
|
|
|
|
for i, ch_name in enumerate(raw_hbo.ch_names):
|
|
if ch_name not in channel_locs:
|
|
channel_locs[ch_name] = []
|
|
all_ch_names.append(ch_name)
|
|
|
|
channel_locs[ch_name].append(raw_hbo.info['chs'][i]['loc'])
|
|
|
|
avg_pairings = {name: np.nanmean(locs, axis=0) for name, locs in channel_locs.items()}
|
|
|
|
optode_collections = {'sources': {}, 'detectors': {}}
|
|
|
|
for ch_name, loc in avg_pairings.items():
|
|
parts = ch_name.split()[0].split('_')
|
|
s_name, d_name = parts[0], parts[1]
|
|
|
|
optode_collections['sources'].setdefault(s_name, []).append(loc[3:6])
|
|
optode_collections['detectors'].setdefault(d_name, []).append(loc[6:9])
|
|
|
|
final_sources = {s: np.nanmean(coords, axis=0) for s, coords in optode_collections['sources'].items()}
|
|
final_detectors = {d: np.nanmean(coords, axis=0) for d, coords in optode_collections['detectors'].items()}
|
|
|
|
ref_raw = raw_list[0].copy().pick(picks="hbo")
|
|
template_lookup = {ch['ch_name']: ch for ch in ref_raw.info['chs']}
|
|
final_chs = []
|
|
|
|
for ch_name in all_ch_names:
|
|
unified_loc = avg_pairings[ch_name].copy()
|
|
parts = ch_name.split()[0].split('_')
|
|
s_name, d_name = parts[0], parts[1]
|
|
|
|
unified_loc[3:6] = final_sources[s_name]
|
|
unified_loc[6:9] = final_detectors[d_name]
|
|
unified_loc[0:3] = (final_sources[s_name] + final_detectors[d_name]) / 2.0
|
|
|
|
# Create the new channel object
|
|
new_ch = template_lookup.get(ch_name, ref_raw.info['chs'][0]).copy()
|
|
new_ch['ch_name'] = ch_name
|
|
new_ch['loc'] = unified_loc
|
|
final_chs.append(new_ch)
|
|
|
|
# Create the final MNE Info
|
|
fake_info = mne.create_info(ch_names=all_ch_names, sfreq=ref_raw.info['sfreq'], ch_types='hbo')
|
|
with fake_info._unlock():
|
|
fake_info['chs'] = final_chs
|
|
|
|
return mne.io.RawArray(np.zeros((len(all_ch_names), 1)), fake_info)
|
|
|
|
|
|
|
|
def brain_3d_visualization(raw_haemo, df_cha, selected_event, t_or_theta: Literal['t', 'theta'] = 'theta', show_optodes: Literal['sensors', 'labels', 'none', 'all'] = 'all', show_text: bool = True, brain_bounds: float = 1.0) -> None:
|
|
|
|
clim = dict(kind="value", pos_lims=(0, brain_bounds/2, brain_bounds))
|
|
|
|
# Get all activity conditions
|
|
for cond in [f'{selected_event}']:
|
|
|
|
ch_summary = df_cha.query(f"Condition.str.startswith('{cond}_delay_') and Chroma == 'hbo'", engine='python') # type: ignore
|
|
|
|
print(ch_summary)
|
|
|
|
if ch_summary.empty:
|
|
#not fir model
|
|
print("No data found for this condition.")
|
|
ch_summary = df_cha.query(f"Condition in [@cond] and Chroma == 'hbo'", engine='python')
|
|
|
|
# Use ordinary least squares (OLS) if only one participant
|
|
# TODO: Fix.
|
|
if True:
|
|
|
|
# t values
|
|
if t_or_theta == 't':
|
|
ch_model = smf.ols("t ~ -1 + ch_name", ch_summary).fit() # type: ignore
|
|
|
|
# theta values
|
|
elif t_or_theta == 'theta':
|
|
ch_model = smf.ols("theta ~ -1 + ch_name", ch_summary).fit() # type: ignore
|
|
|
|
print("OLS model is being used as there is only one participant!")
|
|
|
|
# Convert model results
|
|
model_df = cast(DataFrame, statsmodels_to_results(ch_model, order=ch_summary["ch_name"].unique())) # type: ignore
|
|
|
|
valid_channels = ch_summary["ch_name"].unique().tolist() # type: ignore
|
|
raw_for_plot = raw_haemo.copy().pick(picks=valid_channels) # type: ignore
|
|
|
|
print(f"DEBUG: Model DF rows: {len(model_df)}")
|
|
print(f"DEBUG: Raw channels: {len(raw_for_plot.ch_names)}")
|
|
|
|
brain = plot_3d_evoked_array(raw_for_plot.pick(picks="hbo"), model_df, view="dorsal", distance=0.02, colorbar=True, clim=clim, mode="weighted", size=(800, 700)) # type: ignore
|
|
|
|
if show_optodes == 'all' or show_optodes == 'sensors':
|
|
brain.add_sensors(raw_for_plot.pick(picks="hbo").info, trans=Transform('head', 'mri', np.eye(4)), fnirs=["channels", "pairs", "sources", "detectors"], verbose=False) # type: ignore
|
|
|
|
if True:
|
|
display_text = ('Folder: ' + '\nGroup: ' + '\nCondition: '+ cond + '\nShort Channel Regression: '
|
|
+ '\nLooking at: ' + t_or_theta + ' values') + '\nBrain Distance: '
|
|
|
|
# Apply the text onto the brain
|
|
if show_text:
|
|
brain.add_text(0.12, 0.64, display_text, "title", font_size=11, color="k") # type: ignore
|
|
|
|
return brain
|
|
|
|
|
|
|
|
|
|
def brain_landmarks_3d(raw_haemo: BaseRaw, show_optodes: Literal['sensors', 'labels', 'none', 'all'] = 'all', show_brodmann: bool = True) -> None:
|
|
|
|
brain = Brain("fsaverage", background="white", size=(800, 700)) # type: ignore
|
|
|
|
distances = source_detector_distances(raw_haemo.info)
|
|
|
|
# Add optode text labels manually
|
|
if show_optodes == 'all' or show_optodes == 'sensors':
|
|
brain.add_sensors(getattr(raw_haemo, "info"), trans=Transform('head', 'mri', np.eye(4)), fnirs=["channels", "pairs", "sources", "detectors"], verbose=False) # type: ignore
|
|
|
|
if show_optodes == 'all' or show_optodes == 'labels':
|
|
labeled_srcs = set()
|
|
labeled_dets = set()
|
|
label_counts = {}
|
|
|
|
for idx, ch in enumerate(raw_haemo.info['chs']):
|
|
ch_name = ch['ch_name']
|
|
|
|
if not ch_name.endswith('hbo'):
|
|
continue
|
|
|
|
loc = ch['loc']
|
|
logger.info(f"Channel: {ch_name}")
|
|
logger.info(f"loc length: {len(loc)}")
|
|
logger.info("loc contents:")
|
|
for i, val in enumerate(loc):
|
|
logger.info(f" loc[{i}]: {val}")
|
|
logger.info("-" * 30)
|
|
if not ch_name or not ch['loc'].any():
|
|
continue
|
|
|
|
parts = ch_name.split()[0]
|
|
src_str, det_str = parts.split('_')
|
|
|
|
src_num = int(src_str[1:])
|
|
det_num = int(det_str[1:])
|
|
|
|
if src_num not in labeled_srcs:
|
|
src_xyz = ch['loc'][3:6] * 1000
|
|
brain._renderer.text3d(src_xyz[0], src_xyz[1], src_xyz[2], src_str,
|
|
color='red', scale=0.002)
|
|
labeled_srcs.add(src_num)
|
|
|
|
if det_num not in labeled_dets:
|
|
det_xyz = ch['loc'][6:9] * 1000
|
|
brain._renderer.text3d(det_xyz[0], det_xyz[1], det_xyz[2], det_str,
|
|
color='blue', scale=0.002)
|
|
labeled_dets.add(det_num)
|
|
|
|
# Get the source-detector distance for this channel (in meters)
|
|
dist_m = distances[idx]
|
|
dist_mm = dist_m * 1000
|
|
|
|
label_text = f"{dist_mm:.1f} mm"
|
|
label_counts[label_text] = label_counts.get(label_text, 0) + 1
|
|
if label_counts[label_text] > 1:
|
|
label_text += f" ({label_counts[label_text]})"
|
|
|
|
# Label at channel midpoint
|
|
mid_xyz = loc[0:3] * 1000
|
|
|
|
logger.info(f"Channel: {ch_name} | Midpoint (mm): x={mid_xyz[0]:.2f}, y={mid_xyz[1]:.2f}, z={mid_xyz[2]:.2f} | Distance: {dist_mm:.1f} mm")
|
|
|
|
brain._renderer.text3d(
|
|
mid_xyz[0], mid_xyz[1], mid_xyz[2],
|
|
label_text,
|
|
color='gray',
|
|
scale=0.002
|
|
)
|
|
|
|
|
|
if show_brodmann:# Add Brodmann labels
|
|
labels = cast(list[Label], read_labels_from_annot("fsaverage", "PALS_B12_Brodmann", "lh", verbose=False)) # type: ignore
|
|
|
|
label_colors = {
|
|
"Brodmann.1-lh": "red",
|
|
"Brodmann.2-lh": "red",
|
|
"Brodmann.3-lh": "red",
|
|
"Brodmann.4-lh": "orange",
|
|
"Brodmann.5-lh": "green",
|
|
"Brodmann.6-lh": "yellow",
|
|
"Brodmann.7-lh": "green",
|
|
"Brodmann.17-lh": "blue",
|
|
"Brodmann.18-lh": "blue",
|
|
"Brodmann.19-lh": "blue",
|
|
"Brodmann.39-lh": "pink",
|
|
"Brodmann.40-lh": "purple",
|
|
"Brodmann.42-lh": "white",
|
|
"Brodmann.44-lh": "white",
|
|
"Brodmann.48-lh": "white",
|
|
|
|
}
|
|
|
|
for label in labels:
|
|
name = getattr(label, "name", None)
|
|
if not isinstance(name, str):
|
|
continue
|
|
if name in label_colors:
|
|
brain.add_label(label, borders=False, color=label_colors[name]) # type: ignore
|
|
|
|
|
|
return brain
|
|
|
|
|
|
def verify_channel_positions(data: BaseRaw) -> None:
|
|
"""
|
|
Visualizes the sensor/channel positions of the raw data for verification.
|
|
|
|
Parameters
|
|
----------
|
|
data : BaseRaw
|
|
The loaded data object to process.
|
|
"""
|
|
|
|
def convert_fig_dict_to_png_bytes(fig_dict: dict[str, Figure]) -> dict[str, bytes]:
|
|
png_dict = {}
|
|
for label, fig in fig_dict.items():
|
|
buf = BytesIO()
|
|
fig.savefig(buf, format="png", bbox_inches="tight")
|
|
buf.seek(0)
|
|
png_dict[label] = buf.read()
|
|
plt.close(fig)
|
|
return png_dict
|
|
|
|
|
|
|
|
def brain_3d_contrast(con_model_df: DataFrame, con_model_df_filtered: BaseRaw, common_channels: list[str], first_name: str, second_name: str, t_or_theta: Literal['t', 'theta'] = 'theta', show_optodes: Literal['sensors', 'labels', 'none', 'all'] = 'all', show_text: bool = True, brain_bounds: float = 1.0) -> None:
|
|
# Filter DataFrame to only common channels, and sort by raw order
|
|
con_model = con_model_df
|
|
|
|
con_model["ch_name"] = pd.Categorical(
|
|
con_model["ch_name"], categories=common_channels, ordered=True
|
|
)
|
|
con_model = con_model.sort_values("ch_name").reset_index(drop=True) # type: ignore
|
|
|
|
|
|
clim=dict(kind="value", pos_lims=(0, brain_bounds/2, brain_bounds))
|
|
|
|
# Plot brain figure
|
|
brain = plot_3d_evoked_array(con_model_df_filtered.copy().pick(picks="hbo"), con_model, view="dorsal", distance=0.02, colorbar=True, mode="weighted", clim=clim, size=(800, 700), verbose=False) # type: ignore
|
|
|
|
if show_optodes == 'all' or show_optodes == 'sensors':
|
|
brain.add_sensors(getattr(con_model_df_filtered, "info"), trans=Transform('head', 'mri', np.eye(4)), fnirs=["channels", "pairs", "sources", "detectors"], verbose=False) # type: ignore
|
|
|
|
display_text = ('Contrast: ' + first_name + ' - ' + second_name + '\nLooking at: ' + t_or_theta + ' values')
|
|
|
|
# Apply the text onto the brain
|
|
if show_text:
|
|
brain.add_text(0.12, 0.70, display_text, "title", font_size=11, color="k") # type: ignore
|
|
|
|
|
|
|
|
def plot_2d_3d_contrasts_between_groups(
|
|
contrast_df_a: pd.DataFrame,
|
|
contrast_df_b: pd.DataFrame,
|
|
raw_haemo: BaseRaw,
|
|
group_a_name: str,
|
|
group_b_name: str,
|
|
is_3d: bool = True,
|
|
t_or_theta: Literal['t', 'theta'] = 'theta',
|
|
show_optodes: Literal['sensors', 'labels', 'none', 'all'] = 'all',
|
|
show_text: bool = True,
|
|
brain_bounds: float = 1.0,
|
|
) -> None:
|
|
|
|
|
|
logger.info("-----")
|
|
contrast_df_a = contrast_df_a.copy()
|
|
contrast_df_a["group"] = group_a_name
|
|
contrast_df_b = contrast_df_b.copy()
|
|
contrast_df_b["group"] = group_b_name
|
|
logger.info("-----")
|
|
|
|
df_combined = pd.concat([contrast_df_a, contrast_df_b], ignore_index=True)
|
|
|
|
con_summary = df_combined.query("Chroma == 'hbo'").copy()
|
|
logger.info("-----")
|
|
|
|
valid_channels = (pd.crosstab(con_summary["group"], con_summary["ch_name"]) > 1).all()
|
|
valid_channels = valid_channels[valid_channels].index.tolist()
|
|
con_summary = con_summary[con_summary["ch_name"].isin(valid_channels)]
|
|
logger.info("-----")
|
|
|
|
model_formula = "effect ~ -1 + group:ch_name:Chroma"
|
|
con_model = smf.mixedlm(model_formula, con_summary, groups=con_summary["ID"]).fit(method="nm")
|
|
logger.info("-----")
|
|
|
|
if t_or_theta == "t":
|
|
group1_vals = con_model.tvalues.filter(like=f"group[{group_a_name}]")
|
|
group2_vals = con_model.tvalues.filter(like=f"group[{group_b_name}]")
|
|
else:
|
|
group1_vals = con_model.params.filter(like=f"group[{group_a_name}]")
|
|
group2_vals = con_model.params.filter(like=f"group[{group_b_name}]")
|
|
logger.info("-----")
|
|
|
|
group1_channels = [name.split(":")[1].split("[")[1].split("]")[0] for name in group1_vals.index]
|
|
group2_channels = [name.split(":")[1].split("[")[1].split("]")[0] for name in group2_vals.index]
|
|
|
|
df_group1 = DataFrame({"Coef.": group1_vals.values}, index=group1_channels)
|
|
df_group2 = DataFrame({"Coef.": group2_vals.values}, index=group2_channels)
|
|
|
|
df_contrast = df_group1.join(df_group2, how="inner", lsuffix=f"_{group_a_name}", rsuffix=f"_{group_b_name}")
|
|
logger.info("-----")
|
|
|
|
# A - B
|
|
df_contrast["Coef."] = df_contrast[f"Coef._{group_a_name}"] - df_contrast[f"Coef._{group_b_name}"]
|
|
con_model_df_1_2 = DataFrame({
|
|
"ch_name": df_contrast.index,
|
|
"Coef.": df_contrast["Coef."],
|
|
"Chroma": "hbo"
|
|
})
|
|
logger.info("-----")
|
|
|
|
mne_ch_names = raw_haemo.copy().pick(picks="hbo").ch_names
|
|
glm_ch_names = con_model_df_1_2["ch_name"].tolist()
|
|
common_channels = [ch for ch in mne_ch_names if ch in glm_ch_names]
|
|
|
|
con_model_df_filtered = raw_haemo.copy().pick(picks=common_channels)
|
|
con_model_df_1_2 = con_model_df_1_2.set_index("ch_name").loc[common_channels].reset_index()
|
|
logger.info("-----")
|
|
|
|
if is_3d:
|
|
brain_3d_contrast(
|
|
con_model_df_1_2,
|
|
con_model_df_filtered,
|
|
common_channels,
|
|
group_a_name,
|
|
group_b_name,
|
|
t_or_theta,
|
|
show_optodes,
|
|
show_text,
|
|
brain_bounds
|
|
)
|
|
else:
|
|
plot_glm_group_topo(con_model_df_filtered.copy().pick(picks="hbo"), con_model_df_1_2, names=True, res=128, vlim=(-brain_bounds, brain_bounds)) # type: ignore
|
|
|
|
# TODO: The title currently goes on the colorbar. Low priority
|
|
plt.title(f"Contrast: {group_a_name} vs {group_b_name}") # type: ignore
|
|
plt.show() # type: ignore
|
|
|
|
# plt.title(f"Contrast: {group_a_name} vs {group_b_name}")
|
|
# plt.show()
|
|
|
|
# B - A
|
|
df_contrast["Coef."] = df_contrast[f"Coef._{group_b_name}"] - df_contrast[f"Coef._{group_a_name}"]
|
|
con_model_df_2_1 = DataFrame({
|
|
"ch_name": df_contrast.index,
|
|
"Coef.": df_contrast["Coef."],
|
|
"Chroma": "hbo"
|
|
})
|
|
|
|
glm_ch_names = con_model_df_2_1["ch_name"].tolist()
|
|
common_channels = [ch for ch in mne_ch_names if ch in glm_ch_names]
|
|
|
|
con_model_df_filtered = raw_haemo.copy().pick(picks=common_channels)
|
|
con_model_df_2_1 = con_model_df_2_1.set_index("ch_name").loc[common_channels].reset_index()
|
|
|
|
if is_3d:
|
|
brain_3d_contrast(
|
|
con_model_df_2_1,
|
|
con_model_df_filtered,
|
|
common_channels,
|
|
group_b_name,
|
|
group_a_name,
|
|
t_or_theta,
|
|
show_optodes,
|
|
show_text,
|
|
brain_bounds
|
|
)
|
|
else:
|
|
plot_glm_group_topo(con_model_df_filtered.copy().pick(picks="hbo"), con_model_df_2_1, names=True, res=128, vlim=(-brain_bounds, brain_bounds)) # type: ignore
|
|
|
|
# TODO: The title currently goes on the colorbar. Low priority
|
|
plt.title(f"Contrast: {group_b_name} vs {group_a_name}") # type: ignore
|
|
plt.show() # type: ignore
|
|
|
|
|
|
|
|
|
|
def plot_fir_model_results(df, raw_haemo, dm, selected_event, l_bound, u_bound):
|
|
|
|
|
|
df["isActivity"] = [f"{selected_event}" in n for n in df["Condition"]]
|
|
df["isDelay"] = ["delay" in n for n in df["Condition"]]
|
|
df = df.query("isDelay in [True]")
|
|
df = df.query("isActivity in [True]")
|
|
# Make a new column that stores the condition name for tidier model below
|
|
df.loc[:, "TidyCond"] = ""
|
|
df.loc[df["isActivity"] == True, "TidyCond"] = f"{selected_event}" # noqa: E712
|
|
# Finally, extract the FIR delay in to its own column in data frame
|
|
df.loc[:, "delay"] = [n.split("_")[-1] for n in df.Condition]
|
|
|
|
# To simplify this example we will only look at the activity
|
|
# condition so we now remove the other conditions from the
|
|
# design matrix and GLM results
|
|
dm_cols_activity = np.where([f"{selected_event}" in c for c in dm.columns])[0]
|
|
dm = dm[[dm.columns[i] for i in dm_cols_activity]]
|
|
|
|
try:
|
|
lme = smf.mixedlm("theta ~ -1 + delay:TidyCond:Chroma", df, groups=df["ID"]).fit()
|
|
except:
|
|
lme = smf.ols("theta ~ -1 + delay:TidyCond:Chroma", df, groups=df["ID"]).fit() # type: ignore
|
|
|
|
df_sum = statsmodels_to_results(lme)
|
|
df_sum["delay"] = [int(n) for n in df_sum["delay"]]
|
|
df_sum = df_sum.sort_values("delay")
|
|
|
|
# Print the result for the oxyhaemoglobin data in the target condition
|
|
df_sum.query(f"TidyCond in ['{selected_event}']").query("Chroma in ['hbo']")
|
|
|
|
|
|
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(19, 10))
|
|
|
|
print("dm columns:", dm.columns.tolist())
|
|
|
|
# Extract design matrix columns that correspond to the condition of interest
|
|
dm_cond_idxs = np.where([f"{selected_event}" in n for n in dm.columns])[0]
|
|
dm_cond = dm[[dm.columns[i] for i in dm_cond_idxs]]
|
|
|
|
# Extract the corresponding estimates from the lme dataframe for hbo
|
|
df_hbo = df_sum.query(f"TidyCond in ['{selected_event}']").query("Chroma in ['hbo']")
|
|
vals_hbo = [float(v) for v in df_hbo["Coef."]]
|
|
|
|
# print("--------------------------------------")
|
|
# print(f"dm_cond shape: {dm_cond.shape}")
|
|
# print(f"dm_cond columns: {dm_cond.columns.tolist()}")
|
|
# print(f"vals_hbo length: {len(vals_hbo)}")
|
|
# print(f"vals_hbo sample: {vals_hbo[:5]}")
|
|
# print(f"vals_hbo type: {type(vals_hbo)}")
|
|
# print(f"vals_hbo element type: {type(vals_hbo[0]) if len(vals_hbo) > 0 else 'N/A'}")
|
|
|
|
dm_cond_scaled_hbo = dm_cond * vals_hbo
|
|
|
|
# Extract the corresponding estimates from the lme dataframe for hbr
|
|
df_hbr = df_sum.query(f"TidyCond in ['{selected_event}']").query("Chroma in ['hbr']")
|
|
vals_hbr = [float(v) for v in df_hbr["Coef."]]
|
|
dm_cond_scaled_hbr = dm_cond * vals_hbr
|
|
|
|
first_onset = None
|
|
for desc, onset in zip(raw_haemo.annotations.description, raw_haemo.annotations.onset):
|
|
if selected_event in desc:
|
|
first_onset = onset
|
|
break
|
|
|
|
if first_onset is None:
|
|
raise ValueError(f"Selected event '{selected_event}' not found in annotations.")
|
|
|
|
# Align index values (time axis) to the first occurrence of selected_event
|
|
index_values = dm_cond_scaled_hbo.index - np.ceil(first_onset)
|
|
index_values = np.asarray(index_values)
|
|
|
|
# Plot the result
|
|
axes[0].plot(index_values, np.asarray(dm_cond))
|
|
axes[1].plot(index_values, np.asarray(dm_cond_scaled_hbo))
|
|
axes[2].plot(index_values, np.sum(dm_cond_scaled_hbo, axis=1), "r")
|
|
axes[2].plot(index_values, np.sum(dm_cond_scaled_hbr, axis=1), "b")
|
|
|
|
valid_mask = (index_values >= 0) & (index_values <= 15)
|
|
hbo_sum_window = np.sum(dm_cond_scaled_hbo.loc[valid_mask, :], axis=1)
|
|
peak_idx_in_window = np.argmax(hbo_sum_window)
|
|
peak_idx = np.where(valid_mask)[0][peak_idx_in_window]
|
|
peak_time = float(round(index_values[peak_idx], 2)) # type: ignore
|
|
|
|
axes[2].axvline(x=peak_time, color='k', linestyle='--', linewidth=1.5, label='Peak time') # type: ignore
|
|
|
|
# Format the plot
|
|
for ax in range(3):
|
|
axes[ax].set_xlim(-5, 25)
|
|
axes[ax].set_xlabel("Time (s)")
|
|
axes[0].set_ylim(-0.1, 1.1)
|
|
axes[1].set_ylim(l_bound, u_bound)
|
|
axes[2].set_ylim(l_bound, u_bound)
|
|
axes[0].set_title("FIR Model (Unscaled by GLM estimates)")
|
|
axes[1].set_title(f"FIR Components (Scaled by {selected_event} GLM Estimates)")
|
|
axes[2].set_title(f"Evoked Response {selected_event}")
|
|
axes[0].set_ylabel("FIR Model")
|
|
axes[1].set_ylabel("Oyxhaemoglobin (ΔμMol)")
|
|
axes[2].set_ylabel("Haemoglobin (ΔμMol)")
|
|
axes[2].legend(["Oyxhaemoglobin", "Deoyxhaemoglobin"])
|
|
|
|
# We can also extract the 95% confidence intervals of the estimates too
|
|
l95_hbo = [float(v) for v in df_hbo["[0.025"]] # type: ignore
|
|
u95_hbo = [float(v) for v in df_hbo["0.975]"]] # type: ignore
|
|
dm_cond_scaled_hbo_l95 = dm_cond * l95_hbo
|
|
dm_cond_scaled_hbo_u95 = dm_cond * u95_hbo
|
|
l95_hbr = [float(v) for v in df_hbr["[0.025"]] # type: ignore
|
|
u95_hbr = [float(v) for v in df_hbr["0.975]"]] # type: ignore
|
|
dm_cond_scaled_hbr_l95 = dm_cond * l95_hbr
|
|
dm_cond_scaled_hbr_u95 = dm_cond * u95_hbr
|
|
|
|
axes2: Axes
|
|
fig2, axes2 = plt.subplots(nrows=1, ncols=1, figsize=(7, 7)) # type: ignore
|
|
|
|
# Plot the result
|
|
axes2.plot(index_values, np.sum(dm_cond_scaled_hbo, axis=1), "r") # type: ignore
|
|
axes2.plot(index_values, np.sum(dm_cond_scaled_hbr, axis=1), "b") # type: ignore
|
|
axes2.axvline(x=peak_time, color='k', linestyle='--', linewidth=1.5, label='Peak time') # type: ignore
|
|
|
|
axes2.fill_between( # type: ignore
|
|
index_values,
|
|
np.asarray(np.sum(dm_cond_scaled_hbo_l95, axis=1)),
|
|
np.asarray(np.sum(dm_cond_scaled_hbo_u95, axis=1)),
|
|
facecolor="red",
|
|
alpha=0.25,
|
|
)
|
|
axes2.fill_between( # type: ignore
|
|
index_values,
|
|
np.asarray(np.sum(dm_cond_scaled_hbr_l95, axis=1)),
|
|
np.asarray(np.sum(dm_cond_scaled_hbr_u95, axis=1)),
|
|
facecolor="blue",
|
|
alpha=0.25,
|
|
)
|
|
|
|
# Format the plot
|
|
axes2.set_xlim(-5, 20)
|
|
axes2.set_ylim(l_bound, u_bound)
|
|
axes2.set_title(f"Evoked Response with 95% confidence intervals for )") # type: ignore
|
|
axes2.set_ylabel("Haemoglobin (ΔμMol)") # type: ignore
|
|
axes2.legend(["Oyxhaemoglobin", "Deoyxhaemoglobin", f"Peak {peak_time}s"]) # type: ignore
|
|
axes2.set_xlabel("Time (s)") # type: ignore
|
|
|
|
fig2.tight_layout()
|
|
|
|
fig.show()
|
|
fig2.show()
|
|
|
|
|
|
|
|
def load_snirf(file_path: str) -> tuple[BaseRaw, Figure]:
|
|
"""
|
|
Loads a snirf file, optionally drops channels, downsamples, and creates a figure showing the results.
|
|
|
|
Parameters
|
|
----------
|
|
file_path : str
|
|
Path of the snirf file to load.
|
|
ID : str
|
|
File name of the the snirf file that was loaded.
|
|
drop_prefixes : list[str]
|
|
List of channel name prefixes to drop from the data.
|
|
|
|
Returns
|
|
-------
|
|
tuple[BaseRaw, Figure]
|
|
- BaseRaw: The processed data object.
|
|
- Figure: The corresponding Matplotlib figure.
|
|
"""
|
|
|
|
# Read the snirf file
|
|
raw = read_raw_snirf(file_path, preload=True, verbose=VERBOSITY) # type: ignore
|
|
raw.load_data(verbose=VERBOSITY) # type: ignore
|
|
|
|
# TODO: Why was this commented again?
|
|
# Maybe this should be a bypass parameter?
|
|
|
|
# # Strip the specified amount of seconds from the start of the file
|
|
# total_duration = getattr(raw, "times")[-1]
|
|
# if total_duration > SECONDS_TO_STRIP:
|
|
# raw.crop(tmin=SECONDS_TO_STRIP, tmax=total_duration, verbose=VERBOSITY) # type: ignore
|
|
# logger.info(f"Stripped first {SECONDS_TO_STRIP} second(s) of data.")
|
|
# else:
|
|
# logger.info(f"Data length ({total_duration:.2f}s) less than strip duration; no cropping applied.")
|
|
|
|
# If the user forcibly dropped channels, remove them now before any processing occurs
|
|
# logger.info("Checking if there are channels to forcibly drop...")
|
|
# if drop_prefixes:
|
|
# logger.info("Force dropped channels was specified.")
|
|
# channels_to_drop = [ch for ch in cast(list[str], getattr(raw, "ch_names")) if any(ch.startswith(prefix) for prefix in drop_prefixes)]
|
|
# raw.drop_channels(channels_to_drop, "raise") # type: ignore
|
|
# logger.info("Force dropped channels:", channels_to_drop)
|
|
|
|
# If the user wants to downsample, do it right away
|
|
logger.info("Checking if we should downsample...")
|
|
if DOWNSAMPLE:
|
|
logger.info("Downsample was specified.")
|
|
sfreq_old = getattr(raw, "info")["sfreq"]
|
|
raw.resample(DOWNSAMPLE_FREQUENCY, verbose=VERBOSITY) # type: ignore
|
|
sfreq_new = getattr(raw, "info")["sfreq"]
|
|
logger.info(f"Finished downsampling. Old frequency: {sfreq_old}. New frequency: {sfreq_new}.")
|
|
|
|
logger.info("Successfully loaded the snirf file.")
|
|
|
|
return raw
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
# Data science requires good statistics. This section is requiring work. Am working on figuring this out. Will likely add this to a new viewer.
|
|
|
|
|
|
|
|
|
|
|
|
def run_roi_second_level_analysis(df_roi_all, df_cha_all=None, raw_haemo=None,
|
|
p_threshold=0.05, min_subjects=5,
|
|
correction_method='fdr_bh', target_chroma='hbo',
|
|
graph_bounds=None, roi_config=None,
|
|
threshold_topo=False): # Added parameter
|
|
"""
|
|
Perform group-level ROI analysis, prints stats to console, plots the ROI bar chart,
|
|
and dynamically plots isolated channel-level group topography maps based on a JSON config.
|
|
"""
|
|
# 1. Validation checks
|
|
required_cols = ['ROI', 'Condition', 'Chroma', 'theta', 'ID']
|
|
if not all(col in df_roi_all.columns for col in required_cols):
|
|
raise ValueError(f"Input ROI DataFrame must include: {required_cols}")
|
|
|
|
# 2. Filter ROI data for the targeted chromophore
|
|
df_chroma = df_roi_all[df_roi_all['Chroma'] == target_chroma].copy()
|
|
df_chroma = df_chroma.dropna(subset=['theta'])
|
|
|
|
# 3. Perform 1-sample t-test against zero for each ROI
|
|
rois = df_chroma['ROI'].unique()
|
|
group_results = []
|
|
|
|
for roi in rois:
|
|
roi_data = df_chroma[df_chroma['ROI'] == roi]
|
|
sub_data = roi_data.groupby('ID', as_index=False)['theta'].mean()
|
|
|
|
n_subs = sub_data['ID'].nunique()
|
|
if n_subs < min_subjects:
|
|
continue
|
|
|
|
Y = sub_data['theta'].values
|
|
t_val, p_val = stats.ttest_1samp(Y, 0)
|
|
mean_beta = np.mean(Y)
|
|
std_err = stats.sem(Y)
|
|
|
|
group_results.append({
|
|
'ROI': roi,
|
|
't_val': t_val,
|
|
'p_val': p_val,
|
|
'mean_beta': mean_beta,
|
|
'std_err': std_err,
|
|
'n_subjects': n_subs
|
|
})
|
|
|
|
if not group_results:
|
|
print("\n[ERROR] No ROIs met the minimum subject threshold.\n")
|
|
return pd.DataFrame()
|
|
|
|
df_group = pd.DataFrame(group_results)
|
|
|
|
# 4. Multiple comparisons correction
|
|
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 results table to terminal
|
|
print("\n" + "="*65)
|
|
print(f" GROUP-LEVEL ROI STATISTICAL RESULTS ({target_chroma.upper()})")
|
|
print("="*65)
|
|
df_print = df_group.copy()
|
|
df_print['mean_beta'] = df_print['mean_beta'].apply(lambda x: f"{x:.4f}")
|
|
df_print['t_val'] = df_print['t_val'].apply(lambda x: f"{x:.3f}")
|
|
df_print['p_val'] = df_print['p_val'].apply(lambda x: f"{x:.4f}")
|
|
df_print['p_corrected'] = df_print['p_corrected'].apply(lambda x: f"{x:.4f}")
|
|
print(df_print[['ROI', 'mean_beta', 't_val', 'p_val', 'p_corrected', 'significant']].to_string(index=False))
|
|
print("="*65 + "\n")
|
|
|
|
# 5. Plotting ROI Bar Chart
|
|
sns.set_theme(style="whitegrid")
|
|
fig, ax = plt.subplots(figsize=(8, 6))
|
|
|
|
df_sub_avg = df_chroma.groupby(['ROI', 'ID'], as_index=False)['theta'].mean()
|
|
|
|
sns.barplot(
|
|
data=df_sub_avg, x='ROI', y='theta',
|
|
ax=ax, errorbar=('ci', 95), capsize=0.1,
|
|
color='lightgray', edgecolor='black', linewidth=1.5, zorder=1
|
|
)
|
|
sns.swarmplot(
|
|
data=df_sub_avg, x='ROI', y='theta',
|
|
ax=ax, color='darkblue', size=8, alpha=0.7, zorder=2
|
|
)
|
|
|
|
global_max = df_sub_avg['theta'].max()
|
|
global_min = df_sub_avg['theta'].min()
|
|
y_top = global_max * 1.35 if global_max > 0 else 0.5e-6
|
|
y_bottom = global_min * 1.1 if global_min < 0 else -0.1 * global_max
|
|
ax.set_ylim(y_bottom, y_top)
|
|
|
|
if graph_bounds is not None and graph_bounds > 0.0:
|
|
if graph_bounds < 0.5:
|
|
ax.set_ylim(-graph_bounds, graph_bounds)
|
|
|
|
for idx, row in df_group.iterrows():
|
|
roi_name = row['ROI']
|
|
p_val_corr = row['p_corrected']
|
|
|
|
roi_points = df_sub_avg[df_sub_avg['ROI'] == roi_name]['theta']
|
|
max_y = roi_points.max() if len(roi_points) > 0 else 0
|
|
text_y = max_y + (global_max * 0.03)
|
|
|
|
if p_val_corr < 0.001:
|
|
sig_symbol = "***"
|
|
elif p_val_corr < 0.01:
|
|
sig_symbol = "**"
|
|
elif p_val_corr < p_threshold:
|
|
sig_symbol = "*"
|
|
else:
|
|
sig_symbol = "n.s."
|
|
|
|
sig_text = f"{sig_symbol}\np_corr = {p_val_corr:.3f}"
|
|
x_pos = list(rois).index(roi_name)
|
|
ax.text(
|
|
x_pos, text_y, sig_text,
|
|
ha='center', va='bottom', fontsize=11,
|
|
fontweight='bold', color='red' if p_val_corr < p_threshold else 'gray'
|
|
)
|
|
|
|
ax.axhline(0, color='black', linewidth=1, linestyle='--')
|
|
ax.set_ylabel(r'Hemodynamic Response ($\Delta$ HbO / $\mu$mol/L)' if global_max > 1e-3 else r'Hemodynamic Response ($\Delta$ HbO / mol/L)', 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"Group-Level ROI Activation ({target_chroma.upper()})\nSignificance threshold: p < {p_threshold} {correction_lbl}",
|
|
fontsize=13, fontweight='bold', pad=15
|
|
)
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
# === 6. Segmented Topography Plotting (No Hardcoded Regions) ===
|
|
if df_cha_all is not None and raw_haemo is not None:
|
|
print(f"--> Fitting group-level channel LME for {target_chroma.upper()} topography...")
|
|
try:
|
|
val_col = 'effect' if 'effect' in df_cha_all.columns else 'theta'
|
|
ch_col = 'ch_name' if 'ch_name' in df_cha_all.columns else 'channel'
|
|
|
|
con_summary = df_cha_all[df_cha_all['Chroma'] == target_chroma].copy()
|
|
raw_picked = raw_haemo.copy().pick(picks=target_chroma)
|
|
|
|
# Fit channel LME (suppress ConvergenceWarning locally)
|
|
model_formula = f"{val_col} ~ -1 + {ch_col}:Chroma"
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("ignore", ConvergenceWarning)
|
|
con_model = smf.mixedlm(model_formula, con_summary, groups=con_summary["ID"]).fit(method="nm")
|
|
|
|
# Map statsmodels output to MNE result format
|
|
con_model_df = statsmodels_to_results(con_model, order=raw_picked.ch_names)
|
|
|
|
# --- DYNAMIC ROI PARSING ---
|
|
roi_mapping = {}
|
|
if roi_config is not None:
|
|
raw_json = None
|
|
if isinstance(roi_config, str) and os.path.exists(roi_config):
|
|
with open(roi_config, 'r') as f:
|
|
raw_json = json.load(f)
|
|
elif isinstance(roi_config, dict):
|
|
raw_json = roi_config
|
|
|
|
if raw_json:
|
|
if "regions_of_interest" in raw_json:
|
|
for roi_item in raw_json["regions_of_interest"]:
|
|
roi_name = roi_item.get("name")
|
|
channels = roi_item.get("channels", [])
|
|
if roi_name and channels:
|
|
roi_mapping[roi_name] = channels
|
|
else:
|
|
roi_mapping = raw_json
|
|
|
|
if roi_mapping:
|
|
ch_to_roi = {}
|
|
for roi_name, channels in roi_mapping.items():
|
|
for ch in channels:
|
|
ch_to_roi[ch] = roi_name
|
|
ch_to_roi[ch.split()[0]] = roi_name
|
|
|
|
con_summary['ROI'] = con_summary[ch_col].apply(
|
|
lambda x: ch_to_roi.get(x, ch_to_roi.get(x.split()[0], None) if isinstance(x, str) else None)
|
|
)
|
|
|
|
unique_rois = []
|
|
if 'ROI' in con_summary.columns:
|
|
unique_rois = [r for r in con_summary['ROI'].dropna().unique() if r != ""]
|
|
|
|
if not unique_rois:
|
|
print("--> Warning: No ROI mappings detected. Plotting as a unified grid.")
|
|
unique_rois = ['All_Channels']
|
|
con_summary['ROI'] = 'All_Channels'
|
|
|
|
# Calculate shared symmetric color limits
|
|
vlim = (None, None)
|
|
if 'Coef.' in con_model_df.columns:
|
|
clean_vals = con_model_df['Coef.'].dropna().values
|
|
if len(clean_vals) > 0:
|
|
max_abs = np.max(np.abs(clean_vals))
|
|
if max_abs > 1e-9:
|
|
vlim = (-max_abs, max_abs)
|
|
|
|
fig_topo, ax_topo = plt.subplots(figsize=(6, 6))
|
|
|
|
# Dynamic loop: Plot each region independently
|
|
for i, roi_name in enumerate(unique_rois):
|
|
roi_ch_names = con_summary[con_summary['ROI'] == roi_name][ch_col].unique().tolist()
|
|
roi_ch_names = [ch for ch in roi_ch_names if ch in raw_picked.ch_names]
|
|
|
|
if not roi_ch_names:
|
|
continue
|
|
|
|
raw_roi = raw_picked.copy().pick(picks=roi_ch_names)
|
|
show_colorbar = (i == len(unique_rois) - 1)
|
|
|
|
# === FIX 2: Filter stats dataframe first to prevent "Reducing GLM results..." warnings ===
|
|
roi_con_model_df = con_model_df[con_model_df['ch_name'].isin(roi_ch_names)].copy()
|
|
|
|
plot_glm_group_topo(
|
|
raw_roi,
|
|
roi_con_model_df,
|
|
colorbar=show_colorbar,
|
|
threshold=threshold_topo, # Now uses the parameter!
|
|
axes=ax_topo,
|
|
cmap='RdBu_r',
|
|
vlim=vlim
|
|
)
|
|
|
|
threshold_text = "p < 0.05 Masked" if threshold_topo else "Unthresholded"
|
|
ax_topo.set_title(
|
|
f"Group-Level {target_chroma.upper()} Activation Map\n(Regions Isolated Dynamically, {threshold_text})",
|
|
fontsize=11, fontweight='bold', pad=10
|
|
)
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
except Exception as e:
|
|
logger.error(f"Could not generate topography plot: {e}")
|
|
|
|
return df_group
|
|
|
|
|
|
|
|
|
|
def clean_subject_id(path_or_id):
|
|
"""
|
|
Cleans file paths and ID strings to get a standardized subject identifier.
|
|
E.g., 'C:/path/Sub-01_haemo.snirf' -> 'Sub-01'
|
|
"""
|
|
if not isinstance(path_or_id, str):
|
|
return str(path_or_id)
|
|
base = os.path.basename(path_or_id)
|
|
for ext in ['.snirf', '.nirs', '.fif', '.csv', '.pkl', '_haemo']:
|
|
if base.endswith(ext):
|
|
base = base[:-len(ext)]
|
|
if base.endswith('_haemo'):
|
|
base = base[:-6]
|
|
return base
|
|
|
|
|
|
|
|
def run_cross_group_second_level_analysis(df_roi_all, file_paths_a, file_paths_b,
|
|
group_a_name="Group A", group_b_name="Group B",
|
|
df_cha_all=None, raw_haemo=None,
|
|
p_threshold=0.05, min_subjects=3,
|
|
correction_method='fdr_bh', target_chroma='hbo',
|
|
selected_event=None, graph_bounds=None,
|
|
roi_config=None, threshold_topo=False):
|
|
"""
|
|
Perform cross-group independent statistical analyses (Group A vs Group B),
|
|
renders a grouped bar chart with significance brackets, and plots a group-contrast topography map.
|
|
"""
|
|
# 1. Align IDs and filter dataset to selected Event & Chromophore
|
|
clean_a = set(file_paths_a)
|
|
clean_b = set(file_paths_b)
|
|
|
|
df_roi_all = df_roi_all.copy()
|
|
# df_roi_all['clean_ID'] = df_roi_all['ID'].apply(clean_subject_id)
|
|
df_roi_all['clean_ID'] = df_roi_all['ID']
|
|
|
|
# Filter for active experimental conditions
|
|
df_filtered = df_roi_all[
|
|
(df_roi_all['Chroma'] == target_chroma) &
|
|
(df_roi_all['Condition'] == selected_event)
|
|
].copy()
|
|
|
|
|
|
df_a = df_filtered[df_filtered['clean_ID'].isin(clean_a)].copy()
|
|
df_b = df_filtered[df_filtered['clean_ID'].isin(clean_b)].copy()
|
|
|
|
print(f"DEBUG: Filtering for Event: {selected_event}")
|
|
print(f"DEBUG: Unique IDs in df_filtered: {df_filtered['clean_ID'].unique()}")
|
|
print(f"DEBUG: Clean IDs from Group A: {clean_a}")
|
|
print(f"DEBUG: Clean IDs from Group B: {clean_b}")
|
|
print(f"DEBUG: Rows in df_a: {len(df_a)}, Rows in df_b: {len(df_b)}")
|
|
|
|
|
|
if df_a.empty or df_b.empty:
|
|
print("[ERROR] Missing data for one or both cohorts. Check file selection/IDs.")
|
|
return pd.DataFrame()
|
|
|
|
# 2. ROI-Level Welch's T-Test (Independent Two-Sample)
|
|
rois = df_filtered['ROI'].dropna().unique()
|
|
group_results = []
|
|
|
|
for roi in rois:
|
|
vals_a = df_a[df_a['ROI'] == roi].groupby('clean_ID')['theta'].mean().values
|
|
vals_b = df_b[df_b['ROI'] == roi].groupby('clean_ID')['theta'].mean().values
|
|
|
|
n_a, n_b = len(vals_a), len(vals_b)
|
|
if n_a < min_subjects or n_b < min_subjects:
|
|
continue
|
|
|
|
# Welch's t-test (assumes unequal variances)
|
|
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)
|
|
diff_val = mean_a - mean_b
|
|
|
|
group_results.append({
|
|
'ROI': roi,
|
|
'mean_A': mean_a,
|
|
'mean_B': mean_b,
|
|
'mean_diff': diff_val,
|
|
't_val': t_val,
|
|
'p_val': p_val,
|
|
'n_A': n_a,
|
|
'n_B': n_b
|
|
})
|
|
|
|
if not group_results:
|
|
print("\n[ERROR] No ROIs met the subject requirements for BOTH groups.\n")
|
|
return pd.DataFrame()
|
|
|
|
df_group = pd.DataFrame(group_results)
|
|
|
|
# Apply FDR correction
|
|
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 clean terminal report
|
|
print("\n" + "="*85)
|
|
print(f" CROSS-GROUP ROI CONTRAST: {group_a_name.upper()} vs {group_b_name.upper()} ({target_chroma.upper()})")
|
|
print(f" Event Condition: {selected_event}")
|
|
print("="*85)
|
|
df_print = df_group.copy()
|
|
df_print['mean_A'] = df_print['mean_A'].apply(lambda x: f"{x:.4f}")
|
|
df_print['mean_B'] = df_print['mean_B'].apply(lambda x: f"{x:.4f}")
|
|
df_print['mean_diff'] = df_print['mean_diff'].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']].to_string(index=False))
|
|
print("="*85 + "\n")
|
|
|
|
# 3. Double Grouped Bar Plot (Side-by-Side)
|
|
sns.set_theme(style="whitegrid")
|
|
fig, ax = plt.subplots(figsize=(10, 6))
|
|
|
|
# Construct unified dataframe for seaborn grouped layouts
|
|
df_a_tidy = df_a.groupby(['ROI', 'clean_ID'], as_index=False)['theta'].mean()
|
|
df_a_tidy['Group'] = group_a_name
|
|
df_b_tidy = df_b.groupby(['ROI', 'clean_ID'], as_index=False)['theta'].mean()
|
|
df_b_tidy['Group'] = group_b_name
|
|
combined_df = pd.concat([df_a_tidy, df_b_tidy], ignore_index=True)
|
|
|
|
# Plot Bars
|
|
sns.barplot(
|
|
data=combined_df, x='ROI', y='theta', hue='Group',
|
|
hue_order=[group_a_name, group_b_name], order=rois,
|
|
ax=ax, errorbar=('ci', 95), capsize=0.08,
|
|
palette=['#2b5c8f', '#d95f02'], edgecolor='black', linewidth=1.5, zorder=1
|
|
)
|
|
# Plot Individual Dots (Dodged over the specific bar widths)
|
|
sns.swarmplot(
|
|
data=combined_df, x='ROI', y='theta', hue='Group',
|
|
hue_order=[group_a_name, group_b_name], order=rois,
|
|
ax=ax, size=6, color='black', alpha=0.5, dodge=True, zorder=2,
|
|
legend=False
|
|
)
|
|
|
|
global_max = combined_df['theta'].max()
|
|
global_min = combined_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)
|
|
|
|
if graph_bounds is not None and graph_bounds > 0.0 and graph_bounds < 0.5:
|
|
ax.set_ylim(-graph_bounds, graph_bounds)
|
|
|
|
# Draw professional brackets over paired bars
|
|
for idx, row in df_group.iterrows():
|
|
roi_name = row['ROI']
|
|
p_val_corr = row['p_corrected']
|
|
|
|
roi_points = combined_df[combined_df['ROI'] == roi_name]['theta']
|
|
max_y = roi_points.max() if len(roi_points) > 0 else 0
|
|
|
|
x_a = idx - 0.2 # Approximate left bar X offset
|
|
x_b = idx + 0.2 # Approximate right bar X offset
|
|
y_bracket = max_y + (global_max * 0.08)
|
|
h_tick = global_max * 0.02
|
|
|
|
if p_val_corr < p_threshold:
|
|
sig_symbol = "***" if p_val_corr < 0.001 else "**" if p_val_corr < 0.01 else "*"
|
|
# Draw standard bracket line
|
|
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(
|
|
idx, y_bracket + (global_max * 0.02), f"{sig_symbol}\np_corr = {p_val_corr:.3f}",
|
|
ha='center', va='bottom', fontsize=9, fontweight='bold', color='red'
|
|
)
|
|
else:
|
|
ax.text(
|
|
idx, y_bracket, "n.s.",
|
|
ha='center', va='bottom', fontsize=9, color='gray'
|
|
)
|
|
|
|
ax.axhline(0, color='black', linewidth=1, linestyle='--')
|
|
ax.set_ylabel(r'Hemodynamic Response ($\Delta$ HbO)', fontsize=12)
|
|
ax.set_xlabel('Region of Interest (ROI)', fontsize=12)
|
|
ax.set_title(f"Cross-Group Comparison: {group_a_name} vs {group_b_name}\n({target_chroma.upper()} - {selected_event})", fontsize=13, fontweight='bold', pad=15)
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
# 4. Channel-by-Channel Group-Contrast Topography Map (Zero Hardcoding)
|
|
if df_cha_all is not None and raw_haemo is not None:
|
|
print(f"--> Computing group-level channel contrasts for topography...")
|
|
try:
|
|
val_col = 'effect' if 'effect' in df_cha_all.columns else 'theta'
|
|
ch_col = 'ch_name' if 'ch_name' in df_cha_all.columns else 'channel'
|
|
|
|
# Match channel levels and clean IDs
|
|
con_summary = df_cha_all[
|
|
(df_cha_all['Chroma'] == target_chroma) &
|
|
(df_cha_all['Condition'] == selected_event)
|
|
].copy()
|
|
con_summary['clean_ID'] = con_summary['ID'].apply(clean_subject_id)
|
|
|
|
raw_picked = raw_haemo.copy().pick(picks=target_chroma)
|
|
|
|
# --- Perform manual Channel-by-Channel Two-Sample t-tests ---
|
|
contrast_data = []
|
|
for ch in raw_picked.ch_names:
|
|
ch_a = con_summary[(con_summary['clean_ID'].isin(clean_a)) & (con_summary[ch_col] == ch)]
|
|
ch_b = con_summary[(con_summary['clean_ID'].isin(clean_b)) & (con_summary[ch_col] == ch)]
|
|
|
|
vals_a = ch_a[val_col].dropna().values
|
|
vals_b = ch_b[val_col].dropna().values
|
|
|
|
if len(vals_a) >= min_subjects and len(vals_b) >= min_subjects:
|
|
t_stat, p_val = stats.ttest_ind(vals_a, vals_b, equal_var=False)
|
|
mean_diff = np.mean(vals_a) - np.mean(vals_b)
|
|
else:
|
|
t_stat, p_val, mean_diff = 0.0, 1.0, 0.0
|
|
|
|
contrast_data.append({
|
|
'ch_name': ch,
|
|
'Coef.': mean_diff, # Represents Mean A - Mean B
|
|
't': t_stat,
|
|
'P>|t|': p_val # For threshold masking
|
|
})
|
|
|
|
con_model_df = pd.DataFrame(contrast_data)
|
|
|
|
# --- DYNAMIC ROI PARSING ---
|
|
roi_mapping = {}
|
|
if roi_config is not None and os.path.exists(roi_config):
|
|
with open(roi_config, 'r') as f:
|
|
raw_json = json.load(f)
|
|
if "regions_of_interest" in raw_json:
|
|
for roi_item in raw_json["regions_of_interest"]:
|
|
roi_mapping[roi_item.get("name")] = roi_item.get("channels", [])
|
|
|
|
if roi_mapping:
|
|
ch_to_roi = {}
|
|
for roi_name, channels in roi_mapping.items():
|
|
for ch in channels:
|
|
ch_to_roi[ch] = roi_name
|
|
ch_to_roi[ch.split()[0]] = roi_name
|
|
con_summary['ROI'] = con_summary[ch_col].apply(lambda x: ch_to_roi.get(x, ch_to_roi.get(x.split()[0], None) if isinstance(x, str) else None))
|
|
|
|
unique_rois = [r for r in con_summary['ROI'].dropna().unique() if r != ""] if 'ROI' in con_summary.columns else ['All_Channels']
|
|
|
|
# Shared symmetric limits for the color bar
|
|
max_abs = np.max(np.abs(con_model_df['Coef.'].dropna().values)) if len(con_model_df['Coef.']) > 0 else 1.0
|
|
vlim = (-max_abs, max_abs) if max_abs > 1e-9 else (None, None)
|
|
|
|
fig_topo, ax_topo = plt.subplots(figsize=(6, 6))
|
|
|
|
# Isolated dynamic plotting loop to prevent spatial bleeding
|
|
for i, roi_name in enumerate(unique_rois):
|
|
roi_ch_names = con_summary[con_summary['ROI'] == roi_name][ch_col].unique().tolist() if 'ROI' in con_summary.columns else raw_picked.ch_names
|
|
roi_ch_names = [ch for ch in roi_ch_names if ch in raw_picked.ch_names]
|
|
|
|
if not roi_ch_names:
|
|
continue
|
|
|
|
raw_roi = raw_picked.copy().pick(picks=roi_ch_names)
|
|
show_colorbar = (i == len(unique_rois) - 1)
|
|
|
|
# Filter contrast DF to current ROI channels
|
|
roi_con_df = con_model_df[con_model_df['ch_name'].isin(roi_ch_names)].copy()
|
|
|
|
plot_glm_group_topo(
|
|
raw_roi,
|
|
roi_con_df,
|
|
colorbar=show_colorbar,
|
|
threshold=threshold_topo,
|
|
axes=ax_topo,
|
|
cmap='RdBu_r',
|
|
vlim=vlim
|
|
)
|
|
|
|
threshold_text = "p < 0.05 Masked" if threshold_topo else "Unthresholded Contrast"
|
|
ax_topo.set_title(f"Group Contrast: {group_a_name} - {group_b_name}\n({target_chroma.upper()} - {threshold_text})", fontsize=11, fontweight='bold', pad=10)
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
except Exception as e:
|
|
logger.error(f"Could not generate group-contrast topography plot: {e}", exc_info=True)
|
|
|
|
return df_group
|
|
|
|
|
|
|
|
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
|
|
|
|
|
|
|
|
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,
|
|
target_chroma='hbo', min_subjects=5,
|
|
p_threshold=0.05, correction_method=None,
|
|
roi_a_label=None, roi_b_label=None):
|
|
"""
|
|
Paired within-subject ROI contrast (e.g. contralateral minus ipsilateral
|
|
motor ROI), as a companion to run_roi_second_level_analysis rather than a
|
|
replacement for it. Where run_roi_second_level_analysis tests each ROI's
|
|
theta against zero independently (still contaminated by systemic/global
|
|
physiology shared across the whole head), this function computes, per
|
|
subject, (ROI_A theta - ROI_B theta) for a single condition and tests
|
|
THAT difference against zero. Any systemic component that's roughly equal
|
|
in both ROIs cancels out in the subtraction itself, rather than being
|
|
inferred afterwards by comparing two separate p-values.
|
|
|
|
This is the more powerful, more directly interpretable test whenever you
|
|
already have a specific hypothesis about which two ROIs should differ
|
|
(e.g. laterality) — use run_roi_second_level_analysis for open-ended
|
|
per-ROI screening, and this function for a pre-specified paired
|
|
comparison you want to report as a single confirmatory statistic.
|
|
|
|
Parameters
|
|
----------
|
|
df_roi_all : pd.DataFrame
|
|
Combined individual-level ROI results across subjects.
|
|
Must include: ['ROI', 'Condition', 'Chroma', 'theta', 'ID']
|
|
roi_pairs : tuple(str, str) or list of tuple(str, str)
|
|
One (roi_a, roi_b) pair, or several. Each pair is tested
|
|
independently as (roi_a - roi_b). Passing several pairs lets you
|
|
e.g. test left-hand-tap laterality and right-hand-tap laterality
|
|
(different `condition` values) in one call/figure.
|
|
condition : str or list of str
|
|
The 'Condition' value to filter to for the paired test. If
|
|
`roi_pairs` has multiple pairs and you want a different condition
|
|
per pair, pass a list of the same length as `roi_pairs`; otherwise
|
|
a single value is used for every pair.
|
|
target_chroma : str, default 'hbo'
|
|
Chromophore to test. HbO and HbR should never be tested together.
|
|
min_subjects : int, default 5
|
|
Minimum number of subjects with BOTH ROI_A and ROI_B present (after
|
|
dropping NaNs) required to run the test. Below this, the pair is
|
|
skipped with a warning rather than silently reported.
|
|
p_threshold : float, default 0.05
|
|
Significance threshold applied to the (optionally corrected) p-value.
|
|
correction_method : str or None, default None
|
|
Multiple comparisons correction across the pairs tested in this call
|
|
(statsmodels.stats.multitest.multipletests method name, e.g.
|
|
'fdr_bh'). Left off by default since a single pre-specified paired
|
|
contrast typically doesn't need correction — turn it on if you're
|
|
testing several pairs in the same call and want to control for that.
|
|
roi_a_label, roi_b_label : str or list of str, optional
|
|
Display labels for each pair's ROI_A/ROI_B (defaults to the raw ROI
|
|
names). If testing multiple pairs, pass lists matching `roi_pairs`.
|
|
|
|
Returns
|
|
-------
|
|
pd.DataFrame with one row per tested pair:
|
|
['roi_a', 'roi_b', 'condition', 't_val', 'p_val', 'p_corrected',
|
|
'significant', 'mean_diff', 'n_subjects']
|
|
"""
|
|
|
|
required_cols = ['ROI', 'Condition', 'Chroma', 'theta', 'ID']
|
|
if not all(col in df_roi_all.columns for col in required_cols):
|
|
raise ValueError(f"Input ROI DataFrame must include: {required_cols}")
|
|
|
|
# Normalize inputs to lists so single-pair and multi-pair calls share code.
|
|
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(labels, default_from):
|
|
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)
|
|
|
|
roi_a_labels = _expand_labels(roi_a_label, roi_pairs)
|
|
roi_b_labels = _expand_labels(roi_b_label, roi_pairs)
|
|
|
|
df_chroma = df_roi_all[df_roi_all['Chroma'] == target_chroma].copy()
|
|
df_chroma = df_chroma.dropna(subset=['theta'])
|
|
|
|
results = []
|
|
diff_data_for_plot = [] # keep per-subject diffs around for plotting
|
|
|
|
for (roi_a, roi_b), cond, lbl_a, lbl_b in zip(roi_pairs, conditions, roi_a_labels, roi_b_labels):
|
|
df_cond = df_chroma[df_chroma['Condition'] == cond]
|
|
|
|
a_vals = df_cond[df_cond['ROI'] == roi_a].groupby('ID', as_index=False)['theta'].mean()
|
|
b_vals = df_cond[df_cond['ROI'] == roi_b].groupby('ID', as_index=False)['theta'].mean()
|
|
|
|
# Inner join on ID: only subjects with BOTH ROIs present for this
|
|
# condition contribute to the paired test.
|
|
merged = a_vals.merge(b_vals, on='ID', suffixes=('_a', '_b'))
|
|
merged['diff'] = merged['theta_a'] - merged['theta_b']
|
|
|
|
n_subs = merged['ID'].nunique()
|
|
if n_subs < min_subjects:
|
|
logger.warning(
|
|
f"Skipping pair ({roi_a} - {roi_b}) for condition '{cond}' — "
|
|
f"only {n_subs} subject(s) have both ROIs, need at least {min_subjects}."
|
|
)
|
|
continue
|
|
|
|
Y = merged['diff'].values
|
|
t_val, p_val = stats.ttest_1samp(Y, 0)
|
|
mean_diff = np.mean(Y)
|
|
|
|
results.append({
|
|
'roi_a': roi_a,
|
|
'roi_b': roi_b,
|
|
'label_a': lbl_a or roi_a,
|
|
'label_b': lbl_b or roi_b,
|
|
'condition': cond,
|
|
't_val': t_val,
|
|
'p_val': p_val,
|
|
'mean_diff': mean_diff,
|
|
'n_subjects': n_subs,
|
|
})
|
|
diff_data_for_plot.append(merged.assign(pair=f"{lbl_a or roi_a} - {lbl_b or roi_b}\n({cond})"))
|
|
|
|
if not results:
|
|
print("\n[ERROR] No ROI pairs met the minimum subject threshold.\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" + "=" * 70)
|
|
print(f" PAIRED ROI CONTRAST RESULTS ({target_chroma.upper()})")
|
|
print("=" * 70)
|
|
df_print = df_group.copy()
|
|
df_print['mean_diff'] = df_print['mean_diff'].apply(lambda x: f"{x:.4f}")
|
|
df_print['t_val'] = df_print['t_val'].apply(lambda x: f"{x:.3f}")
|
|
df_print['p_val'] = df_print['p_val'].apply(lambda x: f"{x:.4f}")
|
|
df_print['p_corrected'] = df_print['p_corrected'].apply(lambda x: f"{x:.4f}")
|
|
print(df_print[['label_a', 'label_b', 'condition', 'mean_diff', 't_val',
|
|
'p_val', 'p_corrected', 'significant', 'n_subjects']].to_string(index=False))
|
|
print("=" * 70 + "\n")
|
|
|
|
# --- Plot: one bar per pair, individual subject differences overlaid ---
|
|
sns.set_theme(style="whitegrid")
|
|
fig, ax = plt.subplots(figsize=(max(6, 2.2 * len(results)), 6))
|
|
|
|
plot_df = pd.concat(diff_data_for_plot, ignore_index=True)
|
|
|
|
sns.barplot(
|
|
data=plot_df, x='pair', y='diff', ax=ax,
|
|
errorbar=('ci', 95), capsize=0.1,
|
|
color='lightgray', edgecolor='black', linewidth=1.5, zorder=1
|
|
)
|
|
sns.swarmplot(
|
|
data=plot_df, x='pair', y='diff', ax=ax,
|
|
color='darkblue', size=8, alpha=0.7, zorder=2
|
|
)
|
|
|
|
ax.axhline(0, color='black', linewidth=1, linestyle='--')
|
|
|
|
global_max = plot_df['diff'].max()
|
|
for i, row in df_group.iterrows():
|
|
pair_label = f"{row['label_a']} - {row['label_b']}\n({row['condition']})"
|
|
pair_points = plot_df[plot_df['pair'] == pair_label]['diff']
|
|
max_y = pair_points.max() if len(pair_points) > 0 else 0
|
|
text_y = max_y + (abs(global_max) * 0.08 if global_max else 0.1)
|
|
|
|
p_val_corr = row['p_corrected']
|
|
if p_val_corr < 0.001:
|
|
sig_symbol = "***"
|
|
elif p_val_corr < 0.01:
|
|
sig_symbol = "**"
|
|
elif p_val_corr < p_threshold:
|
|
sig_symbol = "*"
|
|
else:
|
|
sig_symbol = "n.s."
|
|
|
|
ax.text(
|
|
i, text_y, f"{sig_symbol}\np = {p_val_corr:.3f}",
|
|
ha='center', va='bottom', fontsize=11, fontweight='bold',
|
|
color='red' if p_val_corr < p_threshold else 'gray'
|
|
)
|
|
|
|
ax.set_ylabel(r'Paired ROI Difference ($\Delta$ HbO, A $-$ B)', fontsize=12)
|
|
ax.set_xlabel('')
|
|
correction_lbl = f"({correction_method} corrected)" if correction_method else "(uncorrected — single pre-specified contrast)"
|
|
ax.set_title(
|
|
f"Paired ROI Contrast ({target_chroma.upper()})\n"
|
|
f"Significance threshold: p < {p_threshold} {correction_lbl}",
|
|
fontsize=13, fontweight='bold', pad=15
|
|
)
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
return df_group
|
|
|
|
|
|
|
|
def aggregate_channel_contrasts_to_roi(df_contrasts, roi_json_path, weighted=True):
|
|
"""
|
|
Combine already-computed per-channel CONTRAST results (e.g. your
|
|
'2.0_vs_3.0' rows from contrasts.csv / contrast_results) into
|
|
per-subject, per-ROI values — so a joint-fit task contrast can be tested
|
|
at the ROI level using the same one-sample machinery as
|
|
run_roi_second_level_analysis / run_roi_paired_contrast_analysis.
|
|
|
|
This exists because mne_nirs.statistics.RegressionResults has a built-in
|
|
.to_dataframe_region_of_interest() that does inverse-variance-weighted
|
|
channel combination, but the ContrastResults object returned by
|
|
glm_est.compute_contrast() does NOT have that method. This function
|
|
replicates the same weighting logic (weight each channel by the inverse
|
|
of its GLM fit's variance) manually, on the already-exported contrast
|
|
dataframe, rather than requiring you to go back and recompute anything
|
|
from raw GLM objects.
|
|
|
|
Parameters
|
|
----------
|
|
df_contrasts : pd.DataFrame
|
|
Combined per-channel contrast results across subjects/contrasts
|
|
(i.e. your contrasts.csv format). Must include:
|
|
['ch_name', 'effect', 'stat', 'Chroma', 'contrast_name', 'ID']
|
|
`stat` must be the t-statistic (ContrastType == 't'), since standard
|
|
error is recovered as effect / stat.
|
|
roi_json_path : str
|
|
Path to the same regions.json used elsewhere in the pipeline, with
|
|
the structure: {"regions_of_interest": [{"name": ..., "channels": [...]}]}
|
|
`channels` entries should be bare source-detector names (e.g. "S1_D1"),
|
|
matching the convention already used for the GLM-level ROI loading.
|
|
weighted : bool, default True
|
|
If True, combine channels within an ROI using inverse-variance
|
|
weighting (weight = 1 / se^2), matching MNE-NIRS's own default
|
|
behavior for to_dataframe_region_of_interest. If False, channels are
|
|
weighted equally (a plain mean).
|
|
|
|
Returns
|
|
-------
|
|
pd.DataFrame with columns ['ROI', 'Condition', 'Chroma', 'theta', 'ID'],
|
|
directly usable as `df_roi_all` in run_roi_second_level_analysis or
|
|
run_roi_paired_contrast_analysis. 'Condition' holds the contrast name
|
|
(e.g. '2.0_vs_3.0'), and 'theta' holds the ROI-combined contrast effect.
|
|
"""
|
|
|
|
required_cols = ['ch_name', 'effect', 'stat', 'Chroma', 'contrast_name', 'ID']
|
|
print(df_contrasts.columns)
|
|
if not all(col in df_contrasts.columns for col in required_cols):
|
|
raise ValueError(f"Input contrast DataFrame must include: {required_cols}")
|
|
|
|
# --- Load ROI definitions and build a channel-base -> ROI lookup ---
|
|
# Channel base names (e.g. "S1_D1") map to both hbo/hbr rows via the
|
|
# ch_name column ("S1_D1 hbo" / "S1_D1 hbr"), so we key on the base name.
|
|
with open(roi_json_path, 'r') as f:
|
|
roi_data = json.load(f)
|
|
|
|
ch_base_to_roi = {}
|
|
for region in roi_data.get("regions_of_interest", []):
|
|
roi_name = region["name"]
|
|
for ch_base in region["channels"]:
|
|
if ch_base in ch_base_to_roi:
|
|
logger.warning(
|
|
f"Channel '{ch_base}' assigned to multiple ROIs "
|
|
f"('{ch_base_to_roi[ch_base]}' and '{roi_name}') — "
|
|
f"using '{roi_name}' (last one wins)."
|
|
)
|
|
ch_base_to_roi[ch_base] = roi_name
|
|
|
|
df = df_contrasts.copy()
|
|
df['ch_base'] = df['ch_name'].str.split().str[0] # "S1_D1 hbo" -> "S1_D1"
|
|
df['ROI'] = df['ch_base'].map(ch_base_to_roi)
|
|
|
|
n_unassigned = df['ROI'].isna().sum()
|
|
if n_unassigned:
|
|
logger.warning(
|
|
f"{n_unassigned} channel-rows did not match any ROI in "
|
|
f"'{roi_json_path}' and will be excluded."
|
|
)
|
|
df = df.dropna(subset=['ROI'])
|
|
|
|
# Recover standard error from the t-statistic: t = effect / se -> se = effect / t
|
|
with np.errstate(divide='ignore', invalid='ignore'):
|
|
df['se'] = df['effect'] / df['stat']
|
|
# A zero or near-zero t-stat gives an undefined/huge se; drop those rows
|
|
# from the weighting rather than let them explode the ROI average.
|
|
bad_se = ~np.isfinite(df['se']) | (df['se'] == 0)
|
|
if bad_se.any():
|
|
logger.warning(f"Dropping {bad_se.sum()} channel-rows with non-finite "
|
|
f"standard error (t-stat ~ 0) from ROI aggregation.")
|
|
df = df[~bad_se]
|
|
|
|
if weighted:
|
|
df['weight'] = 1.0 / (df['se'] ** 2)
|
|
else:
|
|
df['weight'] = 1.0
|
|
|
|
group_cols = ['ROI', 'contrast_name', 'Chroma', 'ID']
|
|
|
|
def _weighted_mean(g):
|
|
return np.average(g['effect'], weights=g['weight'])
|
|
|
|
roi_theta = (
|
|
df.groupby(group_cols, group_keys=False)
|
|
.apply(lambda g: pd.Series({'theta': _weighted_mean(g)}))
|
|
.reset_index()
|
|
)
|
|
|
|
roi_theta = roi_theta.rename(columns={'contrast_name': 'Condition'})
|
|
return roi_theta[['ROI', 'Condition', 'Chroma', 'theta', 'ID']]
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def collapse_fir_condition_column(df, value_col, condition_col='Condition',
|
|
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
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def calculate_dpf(file_path):
|
|
# order is hbo / hbr
|
|
with h5py.File(file_path, 'r') as f:
|
|
wavelengths = f['/nirs/probe/wavelengths'][:]
|
|
logger.info(f"Wavelengths (nm): {wavelengths}")
|
|
wavelengths = sorted(wavelengths, reverse=True)
|
|
age = float(AGE)
|
|
logger.info(f"Their age was {AGE}")
|
|
# where the hell did I get these from again?
|
|
a = 223.3
|
|
b = 0.05624
|
|
c = 0.8493
|
|
d = -5.723e-7
|
|
e = 0.001245
|
|
f = -0.9025
|
|
dpf = []
|
|
for w in wavelengths:
|
|
logger.info(w)
|
|
dpf.append(a + b * (age**c) + d* (w**3) + e * (w**2) + f*w)
|
|
logger.info(dpf)
|
|
return dpf
|
|
|
|
|
|
|
|
def iqr_threshold(coeffs: NDArray[float64], k: float = 1.5) -> floating[Any]:
|
|
|
|
"""
|
|
Calculate the interquartile range (IQR) threshold scaled by a factor, k.
|
|
|
|
Parameters
|
|
----------
|
|
coeffs : NDArray[float64]
|
|
Array of coefficients to compute the IQR from.
|
|
k : float, optional
|
|
Scaling factor for the IQR (default is 1.5).
|
|
|
|
Returns
|
|
-------
|
|
floating[Any]
|
|
The scaled IQR threshold value.
|
|
"""
|
|
|
|
# Calculate the IQR
|
|
q1 = np.percentile(coeffs, 25)
|
|
q3 = np.percentile(coeffs, 75)
|
|
iqr = q3 - q1
|
|
|
|
return k * iqr
|
|
|
|
|
|
|
|
def wavelet_iqr_denoise(signal: NDArray[float64], wavelet: str = 'db4', level: int = 3) -> NDArray[float64]:
|
|
"""
|
|
Denoises a signal using wavelet decomposition and IQR-based thresholding on detail coefficients.
|
|
|
|
Parameters
|
|
----------
|
|
signal : NDArray[float64]
|
|
The input signal array to denoise.
|
|
wavelet : str, optional
|
|
The type of wavelet to use for decomposition (default is 'db4').
|
|
level : int, optional
|
|
Decomposition level for wavelet transform (default is 3).
|
|
|
|
Returns
|
|
-------
|
|
NDArray[float64]
|
|
The denoised signal array, with the same length as the input.
|
|
"""
|
|
|
|
# Decompose the signal using wavelet transform and initialize a list with approximation coefficients
|
|
coeffs: list[NDArray[float64]] = pywt.wavedec(signal, wavelet, level=level) # type: ignore
|
|
cA = coeffs[0]
|
|
denoised_coeffs = [cA]
|
|
|
|
# Threshold detail coefficients to reduce noise
|
|
for cD in coeffs[1:]:
|
|
threshold = iqr_threshold(cD, IQR)
|
|
cD_thresh = np.sign(cD) * np.maximum(np.abs(cD) - threshold, 0.0) # np.where((cD < lower) | (cD > upper), 0, cD)
|
|
cD_thresh = cD_thresh.astype(float64)
|
|
denoised_coeffs.append(cD_thresh)
|
|
|
|
# Reconstruct the denoised signal
|
|
denoised_signal = cast(NDArray[float64], pywt.waverec(denoised_coeffs, wavelet)) # type: ignore
|
|
return denoised_signal[:len(signal)]
|
|
|
|
|
|
|
|
def calculate_and_apply_wavelet(data: BaseRaw) -> tuple[BaseRaw, Figure]:
|
|
"""
|
|
Applies a wavelet IQR denoising filter to the data and generates a plot.
|
|
|
|
Parameters
|
|
----------
|
|
data : BaseRaw
|
|
The loaded data object to process.
|
|
ID : str
|
|
File name of the the snirf file that was loaded.
|
|
|
|
Returns
|
|
-------
|
|
tuple[BaseRaw, Figure]
|
|
- BaseRaw: The processed data object.
|
|
- Figure: The corresponding Matplotlib figure.
|
|
"""
|
|
|
|
logger.info("Applying the wavelet filter...")
|
|
|
|
# Denoise the data
|
|
logger.info("Denoising the data...")
|
|
loaded_data: NDArray[float64] = data.get_data(verbose=VERBOSITY) # type: ignore
|
|
denoised_data = np.zeros_like(loaded_data)
|
|
|
|
logger.info("Calculating the IQR, decomposing the signal, and thresholding the coefficients...")
|
|
for ch in range(loaded_data.shape[0]):
|
|
denoised_data[ch, :] = wavelet_iqr_denoise(loaded_data[ch, :], wavelet=WAVELET_TYPE, level=WAVELET_LEVEL)
|
|
|
|
# Reconstruct the data with the annotations
|
|
logger.info("Reconstructing the data with annotations...")
|
|
raw_with_tddr_and_wavelet = RawArray(denoised_data, cast(Info, data.info), verbose=VERBOSITY)
|
|
raw_with_tddr_and_wavelet.set_annotations(data.annotations.copy(), verbose=VERBOSITY) # type: ignore
|
|
|
|
# Create a figure for the results
|
|
logger.info("Creating the figure...")
|
|
fig = cast(Figure, raw_with_tddr_and_wavelet.plot(show=False, n_channels=len(getattr(data, "ch_names")), duration=data.times[-1]).figure) # type: ignore
|
|
fig.suptitle(f"Wavelet for ", fontsize=16) # type: ignore
|
|
fig.subplots_adjust(top=0.92)
|
|
plt.close(fig)
|
|
|
|
logger.info("Successfully applied the wavelet filter.")
|
|
|
|
return raw_with_tddr_and_wavelet, fig
|
|
|
|
|
|
|
|
def short_channel_processing_for_hr(data: BaseRaw, short_chans: BaseRaw | None) -> tuple[float, NDArray[float64], NDArray[float64]]:
|
|
"""
|
|
Extract and trim short-channel fNIRS signal for heart rate analysis.
|
|
|
|
Parameters
|
|
----------
|
|
data : BaseRaw
|
|
The loaded data object to process.
|
|
short_chans : BaseRaw | None
|
|
Data object with only short separation channels, or None if unavailable.
|
|
|
|
Returns
|
|
-------
|
|
tuple[float, NDArray[float64], NDArray[float64]]
|
|
- float: Sampling frequency of the signal.
|
|
- NDArray[float64]: Trimmed short-channel signal.
|
|
- NDArray[float64]: Corresponding time values.
|
|
"""
|
|
|
|
# Find the short channel (or best candidate) and extract signal data and sampling frequency
|
|
logger.info("Extracting the signal and calculating the sampling frequency...")
|
|
|
|
# If a short channel exists, use it for our signal. Otherwise just take the first channel in the data
|
|
# TODO: Find a better way around this
|
|
if short_chans is not None:
|
|
signal = cast(NDArray[float64], short_chans.get_data(picks=[0], verbose=VERBOSITY))[0] # type: ignore
|
|
else:
|
|
signal = cast(NDArray[float64], data.get_data(picks=[0], verbose=VERBOSITY))[0] # type: ignore
|
|
|
|
# Calculate the sampling frequency
|
|
sfreq = cast(int, data.info['sfreq'])
|
|
|
|
# Trim start and end of the signal to remove edge artifacts
|
|
logger.info(f"Removing {SECONDS_TO_STRIP_HR} seconds from the beginning and end of the file...")
|
|
strip_samples = int(sfreq * SECONDS_TO_STRIP_HR)
|
|
signal_trimmed = signal[strip_samples:-strip_samples]
|
|
times_trimmed = cast(NDArray[float64], getattr(data, "times"))[strip_samples:-strip_samples]
|
|
|
|
return sfreq, signal_trimmed, times_trimmed
|
|
|
|
|
|
|
|
def calculate_heart_rate_neurokit(sfreq: float, signal_trimmed: NDArray[float64]) -> tuple[NDArray[float64], float]:
|
|
"""
|
|
Calculate and smooth heart rate from a trimmed signal using NeuroKit.
|
|
|
|
Parameters
|
|
----------
|
|
sfreq : float
|
|
Sampling frequency of the signal.
|
|
signal_trimmed : NDArray[float64]
|
|
Preprocessed and trimmed fNIRS signal.
|
|
|
|
Returns
|
|
-------
|
|
tuple[NDArray[float64], float]
|
|
- NDArray[float64]: Smoothed heart rate time series (BPM).
|
|
- float: Mean heart rate.
|
|
"""
|
|
|
|
logger.info("Calculating heart rate using NeuroKit...")
|
|
|
|
# Filter signal to isolate heart rate frequencies and detect peaks
|
|
logger.info("Filtering the signal and detecting peaks...")
|
|
signal_filtered = cast(NDArray[float64], nk.signal_filter(signal_trimmed, sampling_rate=sfreq, lowcut=0.8, highcut=2.5)) # type: ignore
|
|
peaks_dict = cast(dict[str, Any], nk.signal_findpeaks(signal_filtered)) # type: ignore
|
|
peaks = peaks_dict['Peaks']
|
|
hr = cast(NDArray[float64], nk.signal_rate(peaks, sampling_rate=sfreq, desired_length=len(signal_trimmed))) # type: ignore
|
|
hr_clean = np.clip(hr, MAX_LOW_HR, MAX_HIGH_HR)
|
|
|
|
# Smooth heart rate time series by replacing spikes with local rolling mean and calculate the mean
|
|
logger.info("Smoothing the signal and calculating the mean...")
|
|
hr_series = pd.Series(hr_clean)
|
|
local_median = hr_series.rolling(window=SMOOTHING_WINDOW_HR, center=True, min_periods=1).median()
|
|
spikes = hr_series > (local_median + 10)
|
|
smoothed_values = hr_series.copy()
|
|
smoothed_spikes = hr_series.rolling(window=SMOOTHING_WINDOW_HR, center=True, min_periods=1).mean()
|
|
smoothed_values[spikes] = smoothed_spikes[spikes]
|
|
hr_smooth_nk = cast(NDArray[float64], smoothed_values.to_numpy()) # type: ignore
|
|
mean_hr_nk = hr_smooth_nk.mean()
|
|
|
|
logger.info("Original HR min/max: %f, %f", hr_clean.min(), hr_clean.max())
|
|
logger.info("Smoothed HR min/max:%f, %f", hr_smooth_nk.min(), hr_smooth_nk.max())
|
|
logger.info(f"Estimated mean HR nk: {mean_hr_nk:.1f} BPM")
|
|
|
|
logger.info("Successfully calculated heart rate using NeuroKit.")
|
|
|
|
return hr_smooth_nk, mean_hr_nk
|
|
|
|
|
|
|
|
def calculate_heart_rate_scipy(sfreq: float, signal_trimmed: NDArray[float64]) -> tuple[NDArray[floating[Any]], NDArray[float64], np.ndarray[Any, np.dtype[np.bool_]], float]:
|
|
"""
|
|
Estimate heart rate using spectral analysis on a high-pass filtered signal.
|
|
|
|
Parameters
|
|
----------
|
|
sfreq : float
|
|
Sampling frequency of the input signal.
|
|
signal_trimmed : NDArray[float64]
|
|
Trimmed fNIRS signal to analyze.
|
|
|
|
Returns
|
|
-------
|
|
tuple[NDArray[floating[Any]], NDArray[float64], np.ndarray[Any, np.dtype[np.bool_]], float]
|
|
- NDArray[floating[Any]]: Frequencies converted to beats per minute (BPM).
|
|
- NDArray[float64]: Power spectral density (PSD) of the signal.
|
|
- np.ndarray[Any, np.dtype[np.bool_]]: Boolean mask indicating frequencies within heart rate range (30-300 BPM).
|
|
- float: Estimated mean heart rate in BPM corresponding to the PSD peak within the range.
|
|
"""
|
|
|
|
logger.info("Calculating heart rate using SciPy...")
|
|
|
|
# Apply a high-pass Butterworth filter to remove slow trends below 0.5 Hz from the trimmed signal (actual data)
|
|
logger.info("Applying a butterworth filter...")
|
|
b, a = cast(tuple[NDArray[float64], NDArray[float64]], butter(2, 0.5 / (sfreq / 2), btype='high'))
|
|
signal_hp = cast(NDArray[float64],filtfilt(b, a, signal_trimmed))
|
|
|
|
# Calculate the Power Spectral Density (PSD) of the filtered signal using Welch's method
|
|
logger.info("Calculating the PSD...")
|
|
nperseg = min(len(signal_hp), 4096)
|
|
frequencies_scipy, psd_scipy = cast(tuple[NDArray[float64], NDArray[float64]], welch(signal_hp, fs=sfreq, nperseg=nperseg, noverlap=nperseg//2))
|
|
|
|
# Convert frequency values to beats per minute (BPM) and set a heart rate range (30-300 BPM)
|
|
logger.info("Converting to BPM...")
|
|
freq_bpm_scipy = frequencies_scipy * 60
|
|
freq_range_scipy = (freq_bpm_scipy > 30) & (freq_bpm_scipy < 300)
|
|
|
|
# Identify the peak frequency within the heart rate range and estimate the mean heart rate in BPM
|
|
logger.info("Finding the mean...")
|
|
peak_index = np.argmax(psd_scipy[freq_range_scipy])
|
|
mean_hr_scipy = freq_bpm_scipy[freq_range_scipy][peak_index]
|
|
|
|
logger.info("Successfully calculated heart rate using SciPy.")
|
|
|
|
return freq_bpm_scipy, psd_scipy, freq_range_scipy, mean_hr_scipy
|
|
|
|
|
|
def plot_heart_rate(
|
|
freq_bpm_scipy: NDArray[floating[Any]],
|
|
psd_scipy: NDArray[float64],
|
|
freq_range_scipy: np.ndarray[Any, np.dtype[np.bool_]],
|
|
mean_hr_scipy: float,
|
|
hr_smooth_nk: NDArray[floating[Any]],
|
|
mean_hr_nk: float,
|
|
times_trimmed: NDArray[floating[Any]],
|
|
overruled: bool
|
|
) -> tuple[Figure, Figure]:
|
|
"""
|
|
Generate plots comparing heart rate estimates from SciPy PSD and NeuroKit2.
|
|
|
|
Parameters
|
|
----------
|
|
freq_bpm_scipy : NDArray[floating[Any]]
|
|
Frequencies in beats per minute from SciPy PSD analysis.
|
|
psd_scipy : NDArray[float64]
|
|
Power spectral density values corresponding to freq_bpm_scipy.
|
|
freq_range_scipy : np.ndarray[Any, np.dtype[np.bool_]]
|
|
Boolean mask indicating the heart rate frequency range used in PSD.
|
|
mean_hr_scipy : float
|
|
Mean heart rate estimated from SciPy PSD peak.
|
|
hr_smooth_nk : NDArray[floating[Any]]
|
|
Smoothed instantaneous heart rate from NeuroKit2.
|
|
mean_hr_nk : float
|
|
Mean heart rate estimated from NeuroKit2 data.
|
|
times_trimmed : NDArray[floating[Any]]
|
|
Time points corresponding to hr_smooth_nk values.
|
|
overruled: bool
|
|
True if the heart rate from NeuroKit2 is overriding the results from the PSD.
|
|
|
|
Returns
|
|
-------
|
|
tuple[Figure, Figure]
|
|
- Figure showing the PSD and SciPy heart rate estimate.
|
|
- Figure showing the time series comparison of heart rates.
|
|
"""
|
|
|
|
# Create the first plot for the PSD. Add a yellow range to show what we will be filtering to.
|
|
logger.info("Creating the figure...")
|
|
fig1, ax1 = plt.subplots(figsize=(10, 5)) # type: ignore
|
|
ax1.set_xlim(30, 300)
|
|
ax1.plot(freq_bpm_scipy[freq_range_scipy], psd_scipy[freq_range_scipy]) # type: ignore
|
|
ax1.axvline(x=mean_hr_scipy, color='red', linestyle='--', label=f'Mean HR: {mean_hr_scipy:.1f} BPM') # type: ignore
|
|
ax1.axvspan(min(mean_hr_nk - HEART_RATE_WINDOW, mean_hr_scipy - HEART_RATE_WINDOW), max(mean_hr_nk + HEART_RATE_WINDOW, mean_hr_scipy + HEART_RATE_WINDOW), color='yellow', alpha=0.3, label=f'HR Range ±{HEART_RATE_WINDOW} BPM') # type: ignore
|
|
ax1.set_xlabel('Heart Rate (BPM)') # type: ignore
|
|
ax1.set_ylabel('Power Spectral Density') # type: ignore
|
|
ax1.set_title('PSD of fNIRS signal - Peak indicates Heart Rate') # type: ignore
|
|
ax1.grid(True) # type: ignore
|
|
|
|
# Was the value we reported here correct for the data on the graph or was it overruled?
|
|
if overruled:
|
|
note = (
|
|
'\n'
|
|
'Note: Calculation was bad!\n'
|
|
'Data has been set to match\n'
|
|
'the value from NeuroKit2.'
|
|
)
|
|
phantom = Line2D([0], [0], color='none', label=note)
|
|
handles, _ = ax1.get_legend_handles_labels()
|
|
ax1.legend(handles=handles + [phantom]) # type: ignore
|
|
|
|
else:
|
|
ax1.legend() # type: ignore
|
|
plt.close(fig1)
|
|
|
|
# Create the second plot showing the rolling heart rate, as well as the two averages that were calculated
|
|
logger.info("Creating the figure...")
|
|
fig2, ax2 = plt.subplots(figsize=(14, 6)) # type: ignore
|
|
ax2.plot(times_trimmed, hr_smooth_nk, label='Instantaneous HR (NeuroKit2)', color='blue', alpha=0.7) # type: ignore
|
|
ax2.axhline(mean_hr_nk, color='red', linestyle='--', label=f'Mean HR NeuroKit2: {mean_hr_nk:.1f} BPM') # type: ignore
|
|
ax2.axhline(mean_hr_scipy, color='orange', linestyle=':', label=f'SciPy Welch PSD (HP filtered): {mean_hr_scipy:.1f} BPM') # type: ignore
|
|
ax2.set_xlabel('Time (seconds)') # type: ignore
|
|
ax2.set_ylabel('Heart Rate (BPM)') # type: ignore
|
|
ax2.set_title('Heart Rate Estimates Comparison') # type: ignore
|
|
ax2.legend() # type: ignore
|
|
ax2.grid(True) # type: ignore
|
|
fig2.tight_layout()
|
|
plt.close(fig2)
|
|
|
|
return fig1, fig2
|
|
|
|
|
|
# import numpy as np
|
|
|
|
# def mark_bads_by_db_threshold(raw, db_limit=-60):
|
|
# """
|
|
# Converts a dB threshold to absolute power and marks channels
|
|
# exceeding it at 6.25 Hz (1/2 Nyquist).
|
|
# """
|
|
# # 1. Convert your "Magic Number" from dB to absolute units
|
|
# abs_threshold = 10 ** (db_limit / 10)
|
|
|
|
# # 2. Compute PSD specifically around the frequency of interest
|
|
# sfreq = raw.info['sfreq']
|
|
# target_freq = (sfreq / 2) - 1 # 6.25 Hz
|
|
|
|
# # We use a small window around 6.25 Hz to get a stable average
|
|
# spectrum = raw.compute_psd(fmin=target_freq - 0.5, fmax=target_freq + 0.5)
|
|
# psd_data, freqs = spectrum.get_data(return_freqs=True)
|
|
|
|
# # 3. Get the power at that frequency (averaging across the small window)
|
|
# # psd_data shape: (n_channels, n_freqs)
|
|
# power_at_6Hz = np.mean(psd_data, axis=1)
|
|
|
|
# # 4. Identify the "loud" channels
|
|
# bad_indices = np.where(power_at_6Hz > abs_threshold)[0]
|
|
# new_bads = [raw.ch_names[i] for i in bad_indices]
|
|
|
|
# # 5. Update raw.info['bads'] without duplicates
|
|
# # raw.info['bads'] = list(set(raw.info['bads'] + new_bads))
|
|
|
|
# print(f"Threshold: {db_limit} dB -> {abs_threshold:.2e} V^2/Hz")
|
|
# print(f"Newly potentially bad channels: {new_bads}")
|
|
|
|
# return new_bads
|
|
|
|
|
|
|
|
# def find_flat_channels(raw, threshold=1e-15):
|
|
# """
|
|
# Identifies channels that are essentially flat lines (zero variance).
|
|
# """
|
|
# data = raw.get_data()
|
|
# # Calculate standard deviation of each channel
|
|
# std_devs = np.std(data, axis=1)
|
|
|
|
# # Identify channels with almost zero movement
|
|
# flat_idx = np.where(std_devs < threshold)[0]
|
|
# flat_names = [raw.ch_names[i] for i in flat_idx]
|
|
|
|
# return flat_names
|
|
|
|
|
|
# def find_truly_dead_channels(raw, threshold=1e-20):
|
|
# """
|
|
# Looks for channels where the signal literally doesn't move
|
|
# between samples (successive differences are zero).
|
|
# """
|
|
# data = raw.get_data()
|
|
# # Calculate the absolute difference between every sample (t and t+1)
|
|
# # diffs shape: (n_channels, n_times - 1)
|
|
# diffs = np.abs(np.diff(data, axis=1))
|
|
|
|
# # Calculate the mean 'step' size for each channel
|
|
# mean_diffs = np.mean(diffs, axis=1)
|
|
|
|
# # Identify channels where the signal is 'stuck'
|
|
# stuck_idx = np.where(mean_diffs < threshold)[0]
|
|
# stuck_names = [raw.ch_names[i] for i in stuck_idx]
|
|
|
|
# return stuck_names, mean_diffs
|
|
|
|
|
|
# def find_mid_run_flatlines(raw, window_size=10.0, threshold=1e-12):
|
|
# """
|
|
# Checks for channels that go flat partway through the recording.
|
|
# window_size: Duration in seconds to check for flatness.
|
|
# """
|
|
# sfreq = raw.info['sfreq']
|
|
# data = raw.get_data()
|
|
# n_samples_win = int(window_size * sfreq)
|
|
|
|
# bad_channels = []
|
|
|
|
# for i, ch_name in enumerate(raw.ch_names):
|
|
# ch_data = data[i]
|
|
# # Create sliding windows (non-overlapping for speed)
|
|
# windows = np.array_split(ch_data, len(ch_data) // n_samples_win)
|
|
|
|
# # Calculate variance for each window
|
|
# win_vars = [np.var(w) for w in windows]
|
|
|
|
# # If the LAST window (or any significant portion at the end) is flat
|
|
# if win_vars[-1] < threshold:
|
|
# bad_channels.append(ch_name)
|
|
|
|
# return bad_channels
|
|
|
|
|
|
def find_flatline_at_end(raw, threshold_ratio=0.05):
|
|
"""
|
|
Compares the variance of the first 25% of the data to the last 25%.
|
|
If the end variance is less than 5% of the start variance, mark it bad.
|
|
"""
|
|
data = raw.get_data()
|
|
n_samples = data.shape[1]
|
|
quarter = n_samples // 4
|
|
|
|
bad_channels = []
|
|
|
|
for i, ch_name in enumerate(raw.ch_names):
|
|
start_var = np.var(data[i, :quarter])
|
|
end_var = np.var(data[i, -quarter:])
|
|
|
|
# Avoid division by zero for truly dead channels
|
|
if start_var == 0:
|
|
bad_channels.append(ch_name)
|
|
continue
|
|
|
|
ratio = end_var / start_var
|
|
|
|
if ratio < threshold_ratio:
|
|
bad_channels.append(ch_name)
|
|
print(f"Flagged {ch_name}: Variance dropped to {ratio:.2%} of original.")
|
|
|
|
return bad_channels
|
|
|
|
|
|
# def plot_with_death_lines(raw, picks, window_sec=5.0, var_threshold=0.05):
|
|
# sfreq = raw.info['sfreq']
|
|
# data, times = raw.get_data(picks=picks, return_times=True)
|
|
|
|
# fig, ax = plt.subplots(figsize=(12, 6))
|
|
|
|
# for i, ch_name in enumerate(picks):
|
|
# ch_data = data[i]
|
|
# ax.plot(times, ch_data, label=ch_name, alpha=0.8)
|
|
|
|
# # --- Find the 'Death Point' ---
|
|
# # Calculate rolling variance in 5-second chunks
|
|
# win_samples = int(window_sec * sfreq)
|
|
# initial_var = np.var(ch_data[:win_samples])
|
|
|
|
# # Check from the end backwards to find where it "died"
|
|
# death_time = None
|
|
# n_samples = len(ch_data)
|
|
# for start_idx in range(n_samples - win_samples - 1, 0, -win_samples):
|
|
# current_var = np.var(ch_data[start_idx : start_idx + win_samples])
|
|
# if current_var > (initial_var * var_threshold):
|
|
# # The point right after this is where it stays dead
|
|
# death_time = times[start_idx + win_samples]
|
|
# break
|
|
|
|
# if death_time and death_time < (times[-1] - 10): # Only plot if it died early
|
|
# ax.axvline(x=death_time, color='red', linestyle='--', alpha=0.5)
|
|
# ax.text(death_time, ax.get_ylim()[1], 'Signal Loss',
|
|
# color='red', rotation=90, verticalalignment='top')
|
|
|
|
# ax.set_title("HbO/HbR with Automated Signal Loss Detection")
|
|
# ax.set_ylabel("Concentration (Δ μmol)")
|
|
# ax.set_xlabel("Time (s)")
|
|
# ax.legend(loc='lower right')
|
|
# plt.grid(True, alpha=0.3)
|
|
# plt.show(block=True)
|
|
|
|
|
|
# import numpy as np
|
|
|
|
|
|
|
|
def detect_sensor_dropout(raw, threshold_ratio=0.05):
|
|
"""
|
|
Identifies channels where signal variance drops significantly.
|
|
Returns flagged channel names and a summary figure.
|
|
"""
|
|
ch_names = raw.ch_names
|
|
data = raw.get_data()
|
|
n_samples = data.shape[1]
|
|
quarter = n_samples // 4
|
|
|
|
ratios = []
|
|
dead_idx = []
|
|
|
|
for i in range(len(ch_names)):
|
|
start_var = np.var(data[i, :quarter])
|
|
end_var = np.var(data[i, -quarter:])
|
|
|
|
# Handle zero variance (dead from the start)
|
|
if start_var == 0:
|
|
ratio = 0.0
|
|
else:
|
|
ratio = end_var / start_var
|
|
|
|
ratios.append(ratio)
|
|
|
|
if ratio < threshold_ratio:
|
|
dead_idx.append(i)
|
|
print(f"Flagged {ch_names[i]}: Variance dropped to {ratio:.2%} of original.")
|
|
|
|
# Pair-kill logic
|
|
failed_bases = {ch_names[i].split(' ')[0] for i in dead_idx}
|
|
bad_names = [ch for ch in ch_names if ch.split(' ')[0] in failed_bases]
|
|
|
|
# --- Visualization ---
|
|
fig_disp, ax = plt.subplots(figsize=(10, 5), constrained_layout=True)
|
|
|
|
# Color logic: Coral for channels below the threshold
|
|
colors = ['coral' if r < threshold_ratio else 'skyblue' for r in ratios]
|
|
|
|
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 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
|
|
ax.legend()
|
|
|
|
plt.close(fig_disp)
|
|
|
|
print(f"Dropout Check: Flagged {len(failed_bases)} optode pairs.")
|
|
return bad_names, fig_disp
|
|
|
|
|
|
|
|
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.
|
|
"""
|
|
ch_names = raw.ch_names
|
|
sfreq = raw.info['sfreq']
|
|
target_freq = sfreq / freq_div
|
|
|
|
# Compute PSD
|
|
spectrum = raw.compute_psd(fmin=0.1, fmax=sfreq/2)
|
|
psd_data, freqs = spectrum.get_data(return_freqs=True)
|
|
|
|
# Find power near the target frequency
|
|
f_idx = np.where((freqs >= target_freq - 0.2) & (freqs <= target_freq + 0.2))[0]
|
|
power_at_target = np.mean(psd_data[:, f_idx], axis=1)
|
|
|
|
abs_threshold = 10 ** (db_limit / 10)
|
|
noisy_idx = np.where(power_at_target > abs_threshold)[0]
|
|
|
|
# Pair-kill logic
|
|
failed_bases = {ch_names[i].split(' ')[0] for i in noisy_idx}
|
|
bad_names = [ch for ch in ch_names if ch.split(' ')[0] in failed_bases]
|
|
|
|
# --- Visualization ---
|
|
fig, ax = plt.subplots(figsize=(8, 4))
|
|
ax.plot(freqs, 10 * np.log10(psd_data.T), color='gray', alpha=0.2)
|
|
if len(noisy_idx) > 0:
|
|
ax.plot(freqs, 10 * np.log10(psd_data[noisy_idx].T), color='plum', label='Noisy Pairs')
|
|
|
|
ax.axhline(db_limit, color='red', linestyle='--', label='Threshold')
|
|
ax.set_title(f"PSD Noise Analysis (Target: {target_freq}Hz)")
|
|
ax.set_ylabel("Power (dB)")
|
|
ax.legend()
|
|
plt.close(fig)
|
|
|
|
print(f"Noise Check: Flagged {len(failed_bases)} optode pairs.")
|
|
return bad_names, fig
|
|
|
|
|
|
|
|
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)
|
|
|
|
# Calculate Z-Scores
|
|
median_range = np.median(ranges)
|
|
mad = np.median(np.abs(ranges - median_range))
|
|
z_scores = 0.6745 * (ranges - median_range) / (mad if mad > 0 else 1e-15)
|
|
|
|
# Identify failed bases
|
|
failed_indices = np.where(np.abs(z_scores) > threshold)[0]
|
|
failed_bases = {picks[i].split(' ')[0] for i in failed_indices}
|
|
|
|
# Flag entire pairs
|
|
bad_names = [ch for ch in picks if ch.split(' ')[0] in failed_bases]
|
|
|
|
# --- Visualization ---
|
|
fig_swing, ax = plt.subplots(figsize=(8, 4))
|
|
# We color bars by the specific Z-score of that individual channel
|
|
colors = ['coral' if np.abs(z) > threshold else 'skyblue' for z in z_scores]
|
|
|
|
ax.bar(range(len(z_scores)), z_scores, color=colors)
|
|
ax.axhline(threshold, color='red', linestyle='--', label='Outlier Threshold')
|
|
ax.axhline(-threshold, color='red', linestyle='--')
|
|
ax.set_title("Physiological Swing Analysis (Z-Scores)")
|
|
ax.set_ylabel("Standardized Deviation")
|
|
ax.set_xlabel("Channel Index")
|
|
ax.legend()
|
|
plt.close(fig_swing)
|
|
|
|
return bad_names, fig_swing
|
|
|
|
|
|
|
|
def find_bad_channels_coeff_var(raw, coeff_var_threshold=25.0):
|
|
"""
|
|
Identifies bad fNIRS channels using only the Coefficient of Variation (coeff_var).
|
|
"""
|
|
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 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
|
|
coeff_var_scores = (stds / (means + 1e-15)) * 100
|
|
|
|
# Find indices that exceed the threshold
|
|
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_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"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 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}: coeff_var = {coeff_var_scores[ch_idx]:.2f}%")
|
|
else:
|
|
print("All channels passed the coeff_var check.")
|
|
|
|
# --- Visualization ---
|
|
fig_qc, ax = plt.subplots(figsize=(10, 5), constrained_layout=True)
|
|
|
|
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("coeff_var %")
|
|
ax.set_xlabel("Channel Index")
|
|
ax.legend()
|
|
|
|
plt.close(fig_qc)
|
|
|
|
return bad_names, fig_qc
|
|
|
|
|
|
|
|
def hr_calc(raw):
|
|
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)
|
|
hr_smooth_nk, mean_hr_nk = calculate_heart_rate_neurokit(sfreq, signal_trimmed)
|
|
freq_bpm_scipy, psd_scipy, freq_range_scipy, mean_hr_scipy = calculate_heart_rate_scipy(sfreq, signal_trimmed)
|
|
|
|
# HACK: This sucks but looking at the graphs I trust neurokit2 more
|
|
overruled = False
|
|
if mean_hr_scipy < mean_hr_nk - 15:
|
|
mean_hr_scipy = mean_hr_nk
|
|
overruled = True
|
|
if mean_hr_scipy > mean_hr_nk + 15:
|
|
mean_hr_scipy = mean_hr_nk
|
|
overruled = True
|
|
|
|
hr1, hr2 = plot_heart_rate(freq_bpm_scipy, psd_scipy, freq_range_scipy, mean_hr_scipy, hr_smooth_nk, mean_hr_nk, times_trimmed, overruled)
|
|
|
|
fig = raw.plot_psd(show=False)
|
|
raw_filtered = raw.copy().filter(0.5, 3, fir_design='firwin')
|
|
sfreq = raw.info['sfreq']
|
|
data = raw_filtered.get_data()
|
|
channel_names = raw.ch_names
|
|
|
|
# --- Parameters for PSD ---
|
|
desired_bin_hz = 0.1
|
|
nperseg = int(sfreq / desired_bin_hz)
|
|
hr_range = (30, 180) # TODO: Should this not use the user defined values?
|
|
|
|
# --- Function to find strongest local peak ---
|
|
def find_hr_from_psd(ch_data):
|
|
f, Pxx = welch(ch_data, sfreq, nperseg=nperseg)
|
|
mask = (f >= hr_range[0]/60) & (f <= hr_range[1]/60)
|
|
f_masked = f[mask]
|
|
Pxx_masked = Pxx[mask]
|
|
if len(Pxx_masked) < 3:
|
|
return np.nan
|
|
peaks = [i for i in range(1, len(Pxx_masked)-1)
|
|
if Pxx_masked[i] > Pxx_masked[i-1] and Pxx_masked[i] > Pxx_masked[i+1]]
|
|
if not peaks:
|
|
return np.nan
|
|
best_idx = peaks[np.argmax([Pxx_masked[i] for i in peaks])]
|
|
return f_masked[best_idx] * 60 # bpm
|
|
|
|
# --- Compute HR across all channels ---
|
|
hr_all_channels = np.array([find_hr_from_psd(data[i, :]) for i in range(len(channel_names))])
|
|
hr_all_channels = hr_all_channels[~np.isnan(hr_all_channels)]
|
|
hr_mode = np.round(np.median(hr_all_channels)) # Use median if some NaNs
|
|
|
|
print(f"Estimated Heart Rate: {hr_mode} bpm")
|
|
|
|
hr_freq = hr_mode / 60 # Hz
|
|
low = hr_freq - 0.3
|
|
high = hr_freq + 0.3
|
|
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] = {}
|
|
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 Data"] = fig_raw
|
|
if progress_callback: progress_callback(1)
|
|
logger.info("Step 1 Completed.")
|
|
|
|
# Step 2: Trimming
|
|
if TRIM and not FOLDING_BYP:
|
|
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.")
|
|
|
|
# Step 3: Verify Optode Placement
|
|
if OPTODE_PLACEMENT:
|
|
fig_optodes = raw.plot_sensors(show_names=SHOW_OPTODE_NAMES, to_sphere=True, show=False) # type: ignore
|
|
fig_individual["Plot Sensors"] = fig_optodes
|
|
if progress_callback: progress_callback(3)
|
|
logger.info("Step 3 Completed.")
|
|
|
|
# Step 4: Short/Long Channels
|
|
if SHORT_CHANNELS and not FOLDING_BYP:
|
|
#NOTE: Have to split again later but since needed for heart rate, this will stay at step 4. Will split later again.
|
|
_short_chans = get_short_channels(raw, max_dist=SHORT_CHANNELS_THRESHOLD) # JUST FOR PLOTTING THEM SEPERATELY
|
|
fig_short_chans = _short_chans.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="Short Channels Only", show=False)
|
|
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["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.")
|
|
|
|
# Step 6: Scalp Coupling Index
|
|
bad_sci = []
|
|
if SCI and not FOLDING_BYP:
|
|
if HEART_RATE:
|
|
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["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.")
|
|
|
|
# Step 7: Signal to Noise Ratio
|
|
bad_snr = []
|
|
if SNR and not FOLDING_BYP:
|
|
bad_snr, fig_snr = calculate_signal_noise_ratio(raw)
|
|
fig_individual["Signal To Noise Ratio"] = fig_snr
|
|
if progress_callback: progress_callback(7)
|
|
logger.info("Step 7 Completed.")
|
|
|
|
# Step 8: Peak Spectral Power
|
|
bad_psp = []
|
|
if PSP and not FOLDING_BYP:
|
|
bad_psp, fig_psp1, fig_psp2 = calculate_peak_power(raw)
|
|
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.")
|
|
|
|
# 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.")
|
|
|
|
# Step 10: Median Absolute Deviation
|
|
bad_amplitude_range = []
|
|
if MAD and not FOLDING_BYP:
|
|
bad_amplitude_range, fig_range = find_bad_channels_by_amplitude_range(raw, threshold=MAD_THRESHOLD)
|
|
fig_individual['Median Absolute Deviation'] = fig_range
|
|
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_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 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 13: Bad Channels Handling
|
|
if BAD_CHANNELS_HANDLING != "None" and not FOLDING_BYP:
|
|
raw, fig_dropped, fig_raw_before, bad_channels = mark_bads(raw, bad_sci, bad_snr, bad_psp, bad_coeff_var, bad_amplitude_range, bad_noise, bad_disp)
|
|
if fig_dropped and fig_raw_before is not None:
|
|
fig_individual["fig2"] = fig_dropped
|
|
fig_individual["fig3"] = fig_raw_before
|
|
if bad_channels:
|
|
if BAD_CHANNELS_HANDLING == "Interpolate":
|
|
raw, fig_raw_after, fig_compare = interpolate_fNIRS_bads_weighted_average(raw, max_dist=MAX_DIST, min_neighbors=MIN_NEIGHBORS)
|
|
fig_individual["fig4"] = fig_raw_after
|
|
fig_individual["Compare"] = fig_compare
|
|
elif BAD_CHANNELS_HANDLING == "Remove":
|
|
raw = remove_bad_channels(raw, bad_channels)
|
|
if progress_callback: progress_callback(13)
|
|
logger.info("Step 13 Completed.")
|
|
|
|
# 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 15: Temporal Derivative Distribution Repair Filtering
|
|
if TDDR and not FOLDING_BYP:
|
|
raw_od = temporal_derivative_distribution_repair(raw_od)
|
|
fig_raw_od_tddr = raw_od.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="After TDDR (Motion Correction)", show=False)
|
|
fig_individual["TDDR"] = fig_raw_od_tddr
|
|
if progress_callback: progress_callback(15)
|
|
logger.info("Step 15 Completed.")
|
|
|
|
# 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 17: Haemoglobin Concentration
|
|
raw_haemo = beer_lambert_law(raw_od, ppf=calculate_dpf(file_path))
|
|
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 18: Enhance Negative Correlation
|
|
if ENHANCE_NEGATIVE_CORRELATION and not FOLDING_BYP:
|
|
raw_haemo = enhance_negative_correlation(raw_haemo)
|
|
fig_raw_haemo_enc = raw_haemo.plot(duration=raw_haemo.times[-1], n_channels=raw_haemo.info['nchan'], title="Enhance Negative Correlation", show=False)
|
|
fig_individual["ENC"] = fig_raw_haemo_enc
|
|
if progress_callback: progress_callback(18)
|
|
logger.info("Step 18 Completed.")
|
|
|
|
# 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
|
|
fig_individual["filter2"] = fig_raw_haemo_filter
|
|
if progress_callback: progress_callback(19)
|
|
logger.info("Step 19 Completed.")
|
|
|
|
# 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)
|
|
fig_individual["events"] = fig_events
|
|
if progress_callback: progress_callback(20)
|
|
logger.info("Step 20 Completed.")
|
|
|
|
# Step 21: Epoch Calculations
|
|
if not FOLDING_BYP:
|
|
epochs, fig_epochs = epochs_calculations(raw_haemo, events, event_dict)
|
|
for name, fig in fig_epochs:
|
|
fig_individual[f"epochs_{name}"] = fig
|
|
if progress_callback: progress_callback(21)
|
|
logger.info("Step 21 Completed.")
|
|
|
|
# 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 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 24: Generate GLM Results
|
|
if "derivative" not in HRF_MODEL.lower():
|
|
fig_glm_result = plot_glm_results(file_path, raw_haemo, glm_est, df_design_matrix)
|
|
for name, fig in fig_glm_result:
|
|
fig_individual[f"GLM {name}"] = fig
|
|
if progress_callback: progress_callback(24)
|
|
logger.info("24")
|
|
|
|
# Step 25: Generate Channel Results
|
|
df_cha = generate_channel_results(glm_est, file_path)
|
|
if progress_callback: progress_callback(25)
|
|
logger.info("25")
|
|
|
|
# Step 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 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 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")
|
|
|
|
# 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
|
|
|
|
|
|
|
|
|
|
def sanitize_paths_for_pickle(raw_haemo, epochs):
|
|
# Fix raw_haemo._filenames
|
|
if hasattr(raw_haemo, '_filenames'):
|
|
raw_haemo._filenames = [str(p) for p in raw_haemo._filenames]
|
|
|
|
# Fix epochs._raw._filenames
|
|
if hasattr(epochs, '_raw') and hasattr(epochs._raw, '_filenames'):
|
|
epochs._raw._filenames = [str(p) for p in epochs._raw._filenames]
|
|
|
|
|
|
def functional_connectivity_spectral_epochs(epochs, n_lines, vmin):
|
|
|
|
# will crash without this load
|
|
epochs.load_data()
|
|
hbo_epochs = epochs.copy().pick(picks="hbo")
|
|
data = hbo_epochs.get_data()
|
|
names = hbo_epochs.ch_names
|
|
sfreq = hbo_epochs.info["sfreq"]
|
|
con = spectral_connectivity_epochs(
|
|
data,
|
|
method=["coh", "plv"],
|
|
mode="multitaper",
|
|
sfreq=sfreq,
|
|
fmin=0.04,
|
|
fmax=0.2,
|
|
faverage=True,
|
|
verbose=True
|
|
)
|
|
|
|
con_coh, con_plv = con
|
|
|
|
coh = con_coh.get_data(output="dense").squeeze()
|
|
plv = con_plv.get_data(output="dense").squeeze()
|
|
|
|
np.fill_diagonal(coh, 0)
|
|
np.fill_diagonal(plv, 0)
|
|
|
|
plot_connectivity_circle(
|
|
coh,
|
|
names,
|
|
title="fNIRS Functional Connectivity (HbO - Coherence)",
|
|
n_lines=n_lines,
|
|
vmin=vmin
|
|
)
|
|
|
|
|
|
|
|
|
|
|
|
def functional_connectivity_spectral_time(epochs, n_lines, vmin):
|
|
|
|
# will crash without this load
|
|
epochs.load_data()
|
|
hbo_epochs = epochs.copy().pick(picks="hbo")
|
|
data = hbo_epochs.get_data()
|
|
names = hbo_epochs.ch_names
|
|
sfreq = hbo_epochs.info["sfreq"]
|
|
|
|
freqs = np.linspace(0.04, 0.2, 10)
|
|
n_cycles = freqs * 2
|
|
|
|
con = spectral_connectivity_time(
|
|
data,
|
|
freqs=freqs,
|
|
method=["coh", "plv"],
|
|
mode="multitaper",
|
|
sfreq=sfreq,
|
|
fmin=0.04,
|
|
fmax=0.2,
|
|
n_cycles=n_cycles,
|
|
faverage=True,
|
|
verbose=True
|
|
)
|
|
|
|
con_coh, con_plv = con
|
|
|
|
coh = con_coh.get_data(output="dense").squeeze()
|
|
plv = con_plv.get_data(output="dense").squeeze()
|
|
|
|
np.fill_diagonal(coh, 0)
|
|
np.fill_diagonal(plv, 0)
|
|
|
|
plot_connectivity_circle(
|
|
coh,
|
|
names,
|
|
title="fNIRS Functional Connectivity (HbO - Coherence)",
|
|
n_lines=n_lines,
|
|
vmin=vmin
|
|
)
|
|
|
|
|
|
|
|
|
|
def functional_connectivity_envelope(epochs, n_lines, vmin):
|
|
# will crash without this load
|
|
|
|
epochs.load_data()
|
|
hbo_epochs = epochs.copy().pick(picks="hbo")
|
|
data = hbo_epochs.get_data()
|
|
|
|
|
|
env = envelope_correlation(
|
|
data,
|
|
orthogonalize=False,
|
|
absolute=True
|
|
)
|
|
env_data = env.get_data(output="dense")
|
|
|
|
env_corr = env_data.mean(axis=0)
|
|
|
|
env_corr = np.squeeze(env_corr)
|
|
|
|
np.fill_diagonal(env_corr, 0)
|
|
|
|
plot_connectivity_circle(
|
|
env_corr,
|
|
hbo_epochs.ch_names,
|
|
title="fNIRS HbO Envelope Correlation (Task Connectivity)",
|
|
n_lines=n_lines,
|
|
vmin=vmin
|
|
)
|
|
|
|
|
|
def functional_connectivity_betas(raw_hbo, n_lines, vmin, event_name=None):
|
|
|
|
raw_hbo = raw_hbo.copy().pick(picks="hbo")
|
|
onsets = raw_hbo.annotations.onset
|
|
|
|
# CRITICAL: Update the Raw object's annotations so the GLM sees unique events
|
|
ann = raw_hbo.annotations
|
|
new_desc = []
|
|
|
|
for i, desc in enumerate(ann.description):
|
|
new_desc.append(f"{desc}__trial_{i:03d}")
|
|
|
|
ann.description = np.array(new_desc)
|
|
|
|
|
|
# shoudl use user defiuned!!!!
|
|
design_matrix = make_first_level_design_matrix(
|
|
raw=raw_hbo,
|
|
hrf_model='fir',
|
|
fir_delays=np.arange(0, 12, 1),
|
|
drift_model='cosine',
|
|
drift_order=1
|
|
)
|
|
|
|
|
|
# 3. Run GLM & Extract Betas
|
|
glm_results = run_glm(raw_hbo, design_matrix)
|
|
betas = np.array(glm_results.theta())
|
|
reg_names = list(design_matrix.columns)
|
|
|
|
|
|
|
|
|
|
|
|
n_channels = betas.shape[0]
|
|
|
|
# ------------------------------------------------------------------
|
|
# 5. Find unique trial tags (optionally filtered by event)
|
|
# ------------------------------------------------------------------
|
|
trial_tags = sorted({
|
|
col.split("_delay")[0]
|
|
for col in reg_names
|
|
if (
|
|
("__trial_" in col)
|
|
and (event_name is None or col.startswith(event_name + "__"))
|
|
)
|
|
})
|
|
|
|
if len(trial_tags) == 0:
|
|
raise ValueError(f"No trials found for event_name={event_name}")
|
|
|
|
# ------------------------------------------------------------------
|
|
# 6. Build beta series (average across FIR delays per trial)
|
|
# ------------------------------------------------------------------
|
|
beta_series = np.zeros((n_channels, len(trial_tags)))
|
|
|
|
for t, tag in enumerate(trial_tags):
|
|
idx = [
|
|
i for i, col in enumerate(reg_names)
|
|
if col.startswith(f"{tag}_delay")
|
|
]
|
|
beta_series[:, t] = np.mean(betas[:, idx], axis=1).flatten()
|
|
|
|
|
|
|
|
# n_channels, n_trials = betas.shape[0], len(onsets)
|
|
# beta_series = np.zeros((n_channels, n_trials))
|
|
|
|
# for t in range(n_trials):
|
|
# trial_indices = [i for i, col in enumerate(reg_names) if col.startswith(f"trial_{t:03d}_delay")]
|
|
# if trial_indices:
|
|
# beta_series[:, t] = np.mean(betas[:, trial_indices], axis=1).flatten()
|
|
|
|
# Normalize each channel so they are on the same scale
|
|
# Without this, everything is connected to everything. Apparently this is a big issue in fNIRS?
|
|
beta_series = zscore(beta_series, axis=1)
|
|
|
|
global_signal = np.mean(beta_series, axis=0)
|
|
beta_series_clean = np.zeros_like(beta_series)
|
|
for i in range(n_channels):
|
|
slope, _ = np.polyfit(global_signal, beta_series[i, :], 1)
|
|
beta_series_clean[i, :] = beta_series[i, :] - (slope * global_signal)
|
|
|
|
# 4. Correlation & Strict Filtering
|
|
corr_matrix = np.zeros((n_channels, n_channels))
|
|
p_matrix = np.ones((n_channels, n_channels))
|
|
|
|
for i in range(n_channels):
|
|
for j in range(i + 1, n_channels):
|
|
r, p = pearsonr(beta_series_clean[i, :], beta_series_clean[j, :])
|
|
corr_matrix[i, j] = corr_matrix[j, i] = r
|
|
p_matrix[i, j] = p_matrix[j, i] = p
|
|
|
|
# 5. High-Bar Thresholding
|
|
reject, _ = multipletests(p_matrix[np.triu_indices(n_channels, k=1)], method='fdr_bh', alpha=0.05)[:2]
|
|
sig_corr_matrix = np.zeros_like(corr_matrix)
|
|
triu = np.triu_indices(n_channels, k=1)
|
|
|
|
for idx, is_sig in enumerate(reject):
|
|
r_val = corr_matrix[triu[0][idx], triu[1][idx]]
|
|
# Only keep the absolute strongest connections
|
|
if is_sig and abs(r_val) > 0.7:
|
|
sig_corr_matrix[triu[0][idx], triu[1][idx]] = r_val
|
|
sig_corr_matrix[triu[1][idx], triu[0][idx]] = r_val
|
|
|
|
# 6. Plot
|
|
plot_connectivity_circle(
|
|
sig_corr_matrix,
|
|
raw_hbo.ch_names,
|
|
title="Strictly Filtered Connectivity (TDDR + GSR + Z-Score)",
|
|
n_lines=None,
|
|
vmin=0.7,
|
|
vmax=1.0,
|
|
colormap='hot' # Use 'hot' to make positive connections pop
|
|
)
|
|
|
|
|
|
|
|
|
|
def get_single_subject_beta_corr(raw_hbo, event_name=None, config=None):
|
|
"""Processes one participant and returns their correlation matrix."""
|
|
raw_hbo = raw_hbo.copy().pick(picks="hbo")
|
|
ann = raw_hbo.annotations
|
|
|
|
# Rename for trial-level GLM
|
|
new_desc = [f"{desc}__trial_{i:03d}" for i, desc in enumerate(ann.description)]
|
|
ann.description = np.array(new_desc)
|
|
|
|
if config == None:
|
|
print("no config")
|
|
design_matrix = make_first_level_design_matrix(
|
|
raw=raw_hbo, hrf_model='fir',
|
|
fir_delays=np.arange(0, 12, 1),
|
|
drift_model='cosine', drift_order=1
|
|
)
|
|
else:
|
|
print("config")
|
|
if config.get("SHORT_CHANNEL_REGRESSION") == True:
|
|
short_chans = get_short_channels(raw_hbo, max_dist=config.get("SHORT_CHANNELS_THRESHOLD"))
|
|
|
|
design_matrix = make_first_level_design_matrix(
|
|
raw=raw_hbo,
|
|
stim_dur=config.get("STIM_DUR"),
|
|
hrf_model=config.get("HRF_MODEL"),
|
|
drift_model=config.get("DRIFT_MODEL"),
|
|
high_pass=config.get("HIGH_PASS"),
|
|
drift_order=config.get("DRIFT_ORDER"),
|
|
fir_delays=config.get("FIR_DELAYS"),
|
|
add_regs=short_chans.get_data().T,
|
|
add_reg_names=short_chans.ch_names,
|
|
min_onset=config.get("MIN_ONSET"),
|
|
oversampling=config.get("OVERSAMPLING")
|
|
)
|
|
print("yep")
|
|
else:
|
|
design_matrix = make_first_level_design_matrix(
|
|
raw=raw_hbo,
|
|
stim_dur=config.get("STIM_DUR"),
|
|
hrf_model=config.get("HRF_MODEL"),
|
|
drift_model=config.get("DRIFT_MODEL"),
|
|
high_pass=config.get("HIGH_PASS"),
|
|
drift_order=config.get("DRIFT_ORDER"),
|
|
fir_delays=config.get("FIR_DELAYS"),
|
|
min_onset=config.get("MIN_ONSET"),
|
|
oversampling=config.get("OVERSAMPLING")
|
|
)
|
|
|
|
|
|
glm_results = run_glm(raw_hbo, design_matrix)
|
|
betas = np.array(glm_results.theta())
|
|
reg_names = list(design_matrix.columns)
|
|
n_channels = betas.shape[0]
|
|
|
|
# Filter trials by event name
|
|
trial_tags = sorted({
|
|
col.split("_delay")[0] for col in reg_names
|
|
if "__trial_" in col and (event_name is None or col.startswith(event_name + "__"))
|
|
})
|
|
|
|
if not trial_tags:
|
|
return None, None
|
|
|
|
# Build Beta Series
|
|
beta_series = np.zeros((n_channels, len(trial_tags)))
|
|
for t, tag in enumerate(trial_tags):
|
|
idx = [i for i, col in enumerate(reg_names) if col.startswith(f"{tag}_delay")]
|
|
beta_series[:, t] = np.mean(betas[:, idx], axis=1).flatten()
|
|
#beta_series[:, t] = np.max(betas[:, idx], axis=1).flatten() #TODO: Figure out which one to use
|
|
|
|
# Z-score and GSR (Global Signal Regression)
|
|
beta_series = zscore(beta_series, axis=1)
|
|
global_signal = np.mean(beta_series, axis=0)
|
|
for i in range(n_channels):
|
|
slope, _ = np.polyfit(global_signal, beta_series[i, :], 1)
|
|
beta_series[i, :] -= (slope * global_signal)
|
|
|
|
# Correlation Matrix
|
|
corr_matrix = np.corrcoef(beta_series)
|
|
return corr_matrix, raw_hbo.ch_names
|
|
|
|
|
|
def run_group_functional_connectivity(haemo_dict, config_dict, selected_paths, event_name, n_lines, vmin):
|
|
"""Aggregates multiple participants and triggers the plot."""
|
|
all_z_matrices = []
|
|
common_names = None
|
|
|
|
for path in selected_paths:
|
|
raw = haemo_dict.get(path)
|
|
config = config_dict.get(path)
|
|
if raw is None: continue
|
|
print(config)
|
|
|
|
corr, names = get_single_subject_beta_corr(raw, event_name, config)
|
|
|
|
if corr is not None:
|
|
# Fisher Z-transform for averaging
|
|
z_mat = np.arctanh(np.clip(corr, -0.99, 0.99))
|
|
all_z_matrices.append(z_mat)
|
|
common_names = names
|
|
|
|
from scipy.stats import ttest_1samp
|
|
# 1. Convert list to 3D array: (Participants, Channels, Channels)
|
|
group_z_data = np.array(all_z_matrices)
|
|
|
|
print("1")
|
|
# 2. Perform a T-Test across the participant dimension (axis 0)
|
|
# We test if the mean Z-score is different from 0
|
|
# C:\Users\tyler\Desktop\research\.venv\Lib\site-packages\scipy\stats\_axis_nan_policy.py:611: RuntimeWarning: Precision loss occurred in moment calculation due to catastrophic cancellation. This occurs when the data are nearly identical. Results may be unreliable.
|
|
# res = hypotest_fun_out(*samples, axis=axis, **kwds)
|
|
|
|
print("--- Variance Check ---")
|
|
|
|
# ADD THIS LINE: Define n_channels based on the data shape
|
|
# group_z_data.shape is (n_participants, n_channels, n_channels)
|
|
n_channels = group_z_data.shape[1]
|
|
|
|
variance_matrix = np.var(group_z_data, axis=0)
|
|
|
|
# Find where variance is exactly 0 (or very close to it)
|
|
zero_var_indices = np.where(variance_matrix < 1e-15)
|
|
coords = list(zip(zero_var_indices[0], zero_var_indices[1]))
|
|
|
|
diag_count = 0
|
|
non_diag_pairs = []
|
|
|
|
for r, c in coords:
|
|
if r == c:
|
|
diag_count += 1
|
|
else:
|
|
non_diag_pairs.append((r, c))
|
|
|
|
print(f"Total pairs with zero variance: {len(coords)}")
|
|
print(f"Identical diagonals: {diag_count}/{n_channels}")
|
|
|
|
if non_diag_pairs:
|
|
print(f"Warning: {len(non_diag_pairs)} non-diagonal pairs have zero variance!")
|
|
for r, c in non_diag_pairs[:10]: # Print first 10
|
|
print(f" - Pair: Channel {r} & Channel {c}")
|
|
else:
|
|
print("Clean! Zero variance only exists on the diagonals.")
|
|
print("----------------------")
|
|
|
|
t_stats, p_values = ttest_1samp(group_z_data, popmean=0, axis=0)
|
|
print("2")
|
|
|
|
# 3. Multiple Comparisons Correction (FDR)
|
|
# We only care about the upper triangle (unique connections)
|
|
n_channels = p_values.shape[0]
|
|
triu_indices = np.triu_indices(n_channels, k=1)
|
|
flat_p = p_values[triu_indices]
|
|
|
|
reject, corrected_p = multipletests(flat_p, method='fdr_bh', alpha=0.05)[:2]
|
|
|
|
# 4. Create the final "Significant" Matrix
|
|
avg_r = np.tanh(np.mean(group_z_data, axis=0))
|
|
sig_avg_r = np.zeros_like(avg_r)
|
|
|
|
# Only keep connections that are Significant AND above your VMIN (r-threshold)
|
|
for idx, is_sig in enumerate(reject):
|
|
row, col = triu_indices[0][idx], triu_indices[1][idx]
|
|
r_val = avg_r[row, col]
|
|
|
|
if is_sig and abs(r_val) >= vmin:
|
|
sig_avg_r[row, col] = sig_avg_r[col, row] = r_val
|
|
|
|
# 5. Plot the significant results
|
|
|
|
|
|
# if not all_z_matrices:
|
|
# return
|
|
|
|
# # Average and convert back to R
|
|
# avg_z = np.mean(all_z_matrices, axis=0)
|
|
# avg_r = np.tanh(avg_z)
|
|
|
|
# # Thresholding
|
|
# avg_r[np.abs(avg_r) < vmin] = 0
|
|
|
|
plot_connectivity_circle(
|
|
sig_avg_r, common_names, n_lines=n_lines,
|
|
title=f"Group Connectivity: {event_name if event_name else 'All Events'}",
|
|
vmin=vmin, vmax=1.0, colormap='hot'
|
|
) |