tyler/flares
2
0
Fork 0
Code Issues 24 Pull Requests Actions Packages Projects Releases 25 Wiki Activity

unit testing and others

Browse Source
This commit is contained in:
tyler
2026-08-21 00:12:29 -07:00
parent 2f76622e52
commit 7a438b1798
7 changed files with 936 additions and 171 deletions
Download Patch File Download Diff File Expand all files Collapse all files
flares.py
+126 -149
View File
@@ -1565,19 +1565,20 @@ def make_design_matrix(
else:
short_chans = None
# Set the new annotations
raw_haemo.set_annotations(new_annot)
raw_haemo_dm = raw_haemo.copy()
if resample:
raw_haemo.resample(resample_freq, npad="auto")
raw_haemo._data = raw_haemo._data * 1e6
raw_haemo_dm.resample(resample_freq, npad="auto")
try:
short_chans.resample(resample_freq)
except:
pass
raw_haemo_dm._data = raw_haemo_dm._data * 1e6
design_matrix = make_first_level_design_matrix(
raw=raw_haemo,
raw=raw_haemo_dm,
stim_dur=stim_dur,
hrf_model=hrf_model,
drift_model=drift_model,
@@ -1635,7 +1636,7 @@ def make_design_matrix(
fig, ax1 = plt.subplots(figsize=(10, 6), constrained_layout=True)
_ = plot_design_matrix(design_matrix, axes=ax1)
return raw_haemo, design_matrix, fig
return raw_haemo, raw_haemo_dm, design_matrix, fig
@@ -5417,7 +5418,7 @@ def process_participant(file_path, file_start, progress_callback=None):
step_start = lap(step_start, timings, "Step 21")
# Step 22: Design Matrix
raw_haemo, df_design_matrix, fig_design_matrix = make_design_matrix(
raw_haemo, raw_haemo_dm, df_design_matrix, fig_design_matrix = make_design_matrix(
raw_haemo=raw_haemo,
resample=RESAMPLE,
resample_freq=RESAMPLE_FREQ,
@@ -5442,7 +5443,7 @@ def process_participant(file_path, file_start, progress_callback=None):
step_start = lap(step_start, timings, "Step 22")
# Step 23: General Linear Model
glm_est, fig_glm_topo = make_and_run_glm(raw_haemo, df_design_matrix, noise_model=NOISE_MODEL, bins=BINS, n_jobs=N_JOBS, verbosity=VERBOSITY)
glm_est, fig_glm_topo = make_and_run_glm(raw_haemo_dm, df_design_matrix, noise_model=NOISE_MODEL, bins=BINS, n_jobs=N_JOBS, verbosity=VERBOSITY)
_enqueue("GLM Topography", fig_glm_topo, png_queue)
if progress_callback: progress_callback(23)
logger.info("23")
@@ -5508,119 +5509,53 @@ def sanitize_paths_for_pickle(raw_haemo, epochs):
def functional_connectivity_spectral_epochs(
epochs: DataFrame | None,
epochs: Epochs,
n_lines: int,
vmin: float,
) -> None:
# 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"],
con_coh = spectral_connectivity_epochs(
hbo_epochs,
method="coh",
mode="multitaper",
sfreq=sfreq,
sfreq=hbo_epochs.info["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()
coh = np.squeeze(con_coh.get_data(output="dense"))
np.fill_diagonal(coh, 0)
np.fill_diagonal(plv, 0)
plot_connectivity_circle(
coh,
names,
hbo_epochs.ch_names,
title="fNIRS Functional Connectivity (HbO - Coherence)",
n_lines=n_lines,
vmin=vmin
)
def functional_connectivity_spectral_time(
epochs: DataFrame | None,
n_lines: int,
vmin: float,
) -> None:
# 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: DataFrame | None,
epochs: Epochs,
n_lines: int,
vmin: float,
) -> None:
# will crash without this load
epochs.load_data()
hbo_epochs = epochs.copy().pick(picks="hbo")
data = hbo_epochs.get_data()
hbo_epochs.filter(l_freq=0.04, h_freq=0.2, verbose=True)
env = envelope_correlation(
data,
hbo_epochs.get_data(),
orthogonalize=False,
absolute=True
)
env_data = env.get_data(output="dense")
env_corr = env_data.mean(axis=0)
env_corr = np.mean(env.get_data(output="dense"), axis=0)
env_corr = np.squeeze(env_corr)
np.fill_diagonal(env_corr, 0)
plot_connectivity_circle(
env_corr,
hbo_epochs.ch_names,
@@ -5630,107 +5565,149 @@ def functional_connectivity_envelope(
)
def functional_connectivity_spectral_time(
epochs: Epochs,
n_lines: int,
vmin: float,
) -> None:
epochs.load_data()
hbo_epochs = epochs.copy().pick(picks="hbo")
freqs = np.linspace(0.04, 0.2, 10)
n_cycles = freqs * 2
con_coh = spectral_connectivity_time(
hbo_epochs.get_data(),
freqs=freqs,
method="coh",
mode="multitaper",
sfreq=hbo_epochs.info["sfreq"],
fmin=0.04,
fmax=0.2,
n_cycles=n_cycles,
faverage=True,
verbose=True
)
coh = np.squeeze(con_coh.get_data(output="dense"))
if coh.ndim == 3:
coh = coh.mean(axis=0)
np.fill_diagonal(coh, 0)
plot_connectivity_circle(
coh,
hbo_epochs.ch_names,
title="fNIRS Functional Connectivity (HbO - Coherence, Time-Resolved)",
n_lines=n_lines,
vmin=vmin
)
def functional_connectivity_betas(
raw_hbo: BaseRaw,
n_lines: int,
vmin: float,
event_name: str | None = None,
*,
drift_model: str = "cosine",
drift_order: int = 1,
apply_gsr: bool = True,
min_effect_size: float = 0.7,
alpha: float = 0.05,
) -> 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!!!!
ann.description = np.array([
f"{desc}__trial_{i:03d}" for i, desc in enumerate(ann.description)
])
design_matrix = make_first_level_design_matrix(
raw=raw_hbo,
hrf_model='fir',
hrf_model="fir",
fir_delays=np.arange(0, 12, 1),
drift_model='cosine',
drift_order=1
drift_model=drift_model,
drift_order=drift_order,
)
# 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)
glm_results = run_glm(raw_hbo, design_matrix)
betas = np.array(glm_results.theta())
if betas.ndim == 3 and betas.shape[-1] == 1:
betas = betas.squeeze(axis=-1)
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 ("__trial_" in col) and (event_name is None or col.startswith(event_name + "__"))
})
if len(trial_tags) < 4:
raise ValueError(
f"Only {len(trial_tags)} trials found for event_name={event_name}; "
"need at least 4 to compute correlation degrees of freedom."
)
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_idx, tag in enumerate(trial_tags):
col_idx = [j for j, col in enumerate(reg_names) if col.split("_delay")[0] == tag]
beta_series[:, t_idx] = betas[:, col_idx].mean(axis=1)
# Vectorized Global Signal Regression (GSR)
if apply_gsr:
global_signal = np.mean(beta_series, axis=0)
A = np.vstack([global_signal, np.ones(len(global_signal))]).T
# Solve least squares for all channels simultaneously
params, _, _, _ = np.linalg.lstsq(A, beta_series.T, rcond=None)
beta_series_clean = (beta_series.T - A @ params).T
else:
beta_series_clean = beta_series
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)
# --- Vectorized correlation + analytic p-values (replaces the nested
# pearsonr loop below) ---
n_trials = beta_series_clean.shape[1]
corr_matrix = np.corrcoef(beta_series_clean)
with np.errstate(divide='ignore', invalid='ignore'):
t_stats = corr_matrix * np.sqrt((n_trials - 2) / (1 - corr_matrix ** 2))
# Safe t-statistic calculation avoiding division by zero on diagonal
corr_clipped = np.clip(corr_matrix, -0.999999, 0.999999)
t_stats = corr_clipped * np.sqrt((n_trials - 2) / (1 - corr_clipped ** 2))
p_matrix = 2 * t_dist.sf(np.abs(t_stats), df=n_trials - 2)
np.fill_diagonal(p_matrix, 1.0) # diagonal r=1 -> nan/inf guarded explicitly
np.fill_diagonal(p_matrix, 1.0)
triu = np.triu_indices(n_channels, k=1)
flat_p = p_matrix[triu]
reject, _ = multipletests(flat_p, method='fdr_bh', alpha=0.05)[:2]
sig_corr_matrix = np.zeros_like(corr_matrix)
reject, _, _, _ = multipletests(flat_p, method="fdr_bh", alpha=alpha)
sig_corr_matrix = np.zeros_like(corr_matrix)
for idx, is_sig in enumerate(reject):
r_val = corr_matrix[triu[0][idx], triu[1][idx]]
# Only keep the absolute strongest connections
if is_sig and abs(r_val) > 0.7:
if is_sig and abs(r_val) > min_effect_size:
sig_corr_matrix[triu[0][idx], triu[1][idx]] = r_val
sig_corr_matrix[triu[1][idx], triu[0][idx]] = r_val
# 6. Plot
gsr_tag = "GSR" if apply_gsr else "no GSR"
plot_connectivity_circle(
sig_corr_matrix,
sig_corr_matrix,
raw_hbo.ch_names,
title="Strictly Filtered Connectivity (TDDR + GSR + Z-Score)",
n_lines=None,
vmin=0.7,
title=f"Beta-Series Connectivity (FDR q<{alpha}, |r|>{min_effect_size}, {gsr_tag})",
n_lines=n_lines,
vmin=min_effect_size,
vmax=1.0,
colormap='hot' # Use 'hot' to make positive connections pop
colormap="hot",
)
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")