OSCR

Probabilistic Forecasting and Information-Theoretic Analysis of Multivariate fMRI Dynamics.

Code ↔ Paper

20 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 20 matches
  1. [1] § 2. Materials and Methods › 2.4. Forecasting Models › 2.4.4. Transformer ↔ models/transformer.py, lines 19–93 · score 0.81 · transformer encoder layers, attention heads, embedded, positions, vector, stacked
  2. [2] § 2. Materials and Methods › 2.5. Training and Evaluation Protocol ↔ utils/training.py, lines 86–235 · score 0.80 · learning rate scheduling, AdamW, weight decay, optimizer, loss, training
  3. [3] § 2. Materials and Methods › 2.4. Forecasting Models › 2.4.3. Long Short-Term Memory Network ↔ models/lstm.py, lines 18–67 · score 0.80 · LSTM layers, recurrent layers, hidden state, multi layer, dropout, batch
  4. [4] § 2. Materials and Methods › 2.4. Forecasting Models › 2.4.3. Long Short-Term Memory Network ↔ models/lstm.py, lines 18–67 · score 0.74 · recurrent layer, hidden state, full forecasting, reshaped, maps, linear
  5. [5] § 2. Materials and Methods › 2.4. Forecasting Models › 2.4.4. Transformer ↔ scripts/cross_validation.py, lines 60–108 · score 0.73 · transformer encoder layers, attention heads, forecasting horizon, dropout, dimension, training
  6. [6] § 2. Materials and Methods › 2.4. Forecasting Models › 2.4.4. Transformer ↔ models/transformer.py, lines 19–93 · score 0.73 · Learned positional embeddings, ROI vector, layer, space, transformer, linear
  7. [7] § 2. Materials and Methods › 2.2. Forecasting Problem Formulation ↔ utils/measures.py, lines 741–806 · score 0.69 · Empirical residual histograms, predictive distribution, model predicted, binning, discretized, forecasting models
  8. [8] § 2. Materials and Methods › 2.1. Dataset and Region-of-Interest Definition ↔ utils/nsd_utils.py, lines 353–433 · score 0.66 · MNI space, ROI atlas, ANTs, affine, scan, warped
  9. [9] § 2. Materials and Methods › 2.1. Dataset and Region-of-Interest Definition ↔ utils/nsd_utils.py, lines 353–433 · score 0.66 · MNI space, ROI atlas, ANTs, affine, scan, warped
  10. [10] § 2. Materials and Methods › 2.5. Training and Evaluation Protocol ↔ utils/training.py, lines 430–488 · score 0.66 · LOSO folds, LOSO CV, cross validation, forecasting models, training
  11. [11] § 2. Materials and Methods › 2.6. Directed Information Post Analysis ↔ utils/measures.py, lines 741–806 · score 0.65 · empirical residual histogram, predictive distributions, model prediction, probability, horizon, ROI
  12. [12] § 2. Materials and Methods › 2.1. Dataset and Region-of-Interest Definition ↔ utils/nsd_utils.py, lines 1095–1143 · score 0.65 · LO1, VO1, LO2, PPA, V3A, VO2
  13. [13] § 2. Materials and Methods › 2.1. Dataset and Region-of-Interest Definition ↔ utils/nsd_utils.py, lines 1090–1138 · score 0.65 · LO1, VO1, LO2, PPA, V3A, VO2
  14. [14] § 2. Materials and Methods › 2.2. Forecasting Problem Formulation ↔ utils/measures.py, lines 168–172 · score 0.64 · Squared Scaled Error, step forecasting, RMSSE, Root, predictive
  15. [15] § 3. Results › Post Hoc Directed Information Analysis ↔ scripts/compute_di.py, lines 75–135 · score 0.59 · DI computation, source ROI, target ROI, likelihoods, probabilistic, training
  16. [16] § 2. Materials and Methods › 2.5. Training and Evaluation Protocol ↔ scripts/cross_validation.py, lines 111–234 · score 0.58 · cross validation, LOSO CV, fold, held, training, Model
  17. [17] § 2. Materials and Methods › 2.4. Forecasting Models › 2.4.2. Exponential Smoothing ↔ models/exponential_smoothing.py, lines 1–12 · score 0.57 · exponential smoothing model, fMRI, global, predictions, forecasting
  18. [18] § 3. Results ↔ utils/parse_data.py, lines 334–393 · score 0.56 · phase randomization, amplitude, power, spectrum, surrogate, variance
  19. [19] § 2. Materials and Methods › 2.3. Data Parsing and Sliding-Window Construction ↔ utils/parse_data.py, lines 243–331 · score 0.53 · sliding window, target windows, overlapping, Parsing
  20. [20] § 2. Materials and Methods › 2.5. Training and Evaluation Protocol ↔ utils/training.py, lines 430–488 · score 0.51 · LOSO folds, cross validated, Training, predictability, Model, forecasting

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 643 lines · 21 KB · Apache-2.0 · 3 matches

  1. import copy
  2. import inspect
  3. import pandas as pd
  4. import numpy as np
  5. from tqdm.auto import tqdm
  6. from pathlib import Path
  8. import torch
  9. import torch.nn as nn
  10. from torch.utils.data import Dataset, DataLoader
  12. from .measures import compute_eta_gauss, compute_rmse, compute_rmsse
  14. from .parse_data import split_by_subject, normalize_items, build_sliding_windows
  16. try:
  17. from models.naive_models import last_value_model_generator
  18. except ImportError:
  19. from ..models.naive_models import last_value_model_generator
  22. device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
  25. class FMRIWindowDataset(Dataset):
  26. """PyTorch dataset for fMRI windowed sequences."""
  28. def __init__(self, X, Y=None):
  29. self.X = torch.tensor(X, dtype=torch.float32)
  30. self.Y = None
  31. if Y is not None:
  32. self.Y = torch.tensor(Y, dtype=torch.float32)
  34. def __len__(self):
  35. return len(self.X)
  37. def __getitem__(self, idx):
  38. if self.Y is not None:
  39. return self.X[idx], self.Y[idx]
  40. return self.X[idx]
  43. class DeltaAwareLoss(nn.Module):
  44. """
  45. Combined loss: HuberLoss + delta (change) penalty.
  47. Combines a robust HuberLoss on absolute values with a penalty on incorrect
  48. forecast dynamics (differences between consecutive timesteps). This encourages
  49. the model to not only predict accurate values but also preserve the temporal
  50. structure and rate of change in the fMRI time series.
  52. Attributes:
  53. alpha (float): Weight for the delta loss term. Controls the trade-off between
  54. prediction accuracy and dynamics preservation.
  55. base (nn.HuberLoss): Huber loss instance used for both absolute and delta terms.
  57. Input shapes:
  58. - pred: (batch_size, horizon, n_roi) predicted fMRI connectivity values
  59. - target: (batch_size, horizon, n_roi) ground truth fMRI connectivity values
  61. Output:
  62. - scalar tensor: weighted combination of base loss and delta loss
  64. Example:
  65. >>> criterion = DeltaAwareLoss(alpha=0.3, delta=0.5)
  66. >>> pred = torch.randn(32, 5, 100) # batch=32, horizon=5, n_roi=100
  67. >>> target = torch.randn(32, 5, 100)
  68. >>> loss = criterion(pred, target) # returns scalar
  69. """
  71. def __init__(self, alpha=0.3, delta=0.5):
  72. super().__init__()
  73. self.alpha = alpha
  74. self.base = nn.HuberLoss(delta=delta)
  76. def forward(self, pred, target):
  77. base_loss = self.base(pred, target)
  79. pred_delta = pred[:, 1:, :] - pred[:, :-1, :]
  80. target_delta = target[:, 1:, :] - target[:, :-1, :]
  81. delta_loss = self.base(pred_delta, target_delta)
  83. return base_loss + self.alpha * delta_loss
  86. def train_model(
  87. model,
  88. train_loader,
  89. val_loader=None,
  90. num_epochs=30,
  91. device=None,
  92. patience=5,
  93. checkpoint_dir=None,
  94. checkpoint_prefix="forecast_model",
  95. checkpoint_every=None,
  96. save_best=True,
  97. save_last=False,
  98. verbose=True
  99. ):
  100. if device is None:
  101. device = torch.device("cpu")
  102. elif isinstance(device, str):
  103. device = torch.device(device)
  105. if type(checkpoint_dir) is str:
  106. checkpoint_dir = Path(checkpoint_dir)
  107. if checkpoint_dir is not None:
  108. checkpoint_dir.mkdir(parents=True, exist_ok=True)
  110. optimizer = torch.optim.AdamW(
  111. model.parameters(), lr=5e-4, weight_decay=1e-5
  112. )
  113. scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(
  114. optimizer, mode="min", factor=0.5, patience=patience
  115. )
  117. criterion = DeltaAwareLoss(alpha=0.3, delta=0.5)
  119. best_val_loss = float("inf")
  120. patience_counter = 0
  121. best_state = None
  123. def _save_checkpoint(path, epoch, train_loss, val_loss=None, is_best=False):
  124. if checkpoint_dir is None:
  125. return
  126. torch.save(
  127. {
  128. "epoch": epoch,
  129. "model_state_dict": model.state_dict(),
  130. "optimizer_state_dict": optimizer.state_dict(),
  131. "scheduler_state_dict": scheduler.state_dict(),
  132. "train_loss": train_loss,
  133. "val_loss": val_loss,
  134. "best_val_loss": best_val_loss,
  135. "is_best": is_best,
  136. },
  137. path,
  138. )
  140. for epoch in range(num_epochs):
  141. model.train()
  142. total_loss = 0.0
  143. loop = tqdm(train_loader, desc=f"Epoch {epoch+1}/{num_epochs}", leave=False)
  144. for X_batch, Y_batch in loop:
  145. X_batch = X_batch.to(device, non_blocking=True)
  146. Y_batch = Y_batch.to(device, non_blocking=True)
  147. optimizer.zero_grad()
  148. loss = criterion(model(X_batch), Y_batch)
  149. loss.backward()
  150. optimizer.step()
  151. total_loss += loss.item()
  152. loop.set_postfix(loss=f"{loss.item():.4f}")
  154. avg_loss = total_loss / len(train_loader)
  155. scheduler.step(avg_loss)
  157. if val_loader is not None:
  158. model.eval()
  159. val_loss = 0.0
  160. with torch.no_grad():
  161. for xb, yb in val_loader:
  162. val_loss += criterion(
  163. model(xb.to(device)), yb.to(device)
  164. ).item()
  165. val_loss /= len(val_loader)
  167. if verbose:
  168. print(f" Epoch {epoch+1:2d} | train: {avg_loss:.6f} | val: {val_loss:.6f}")
  170. if val_loss < best_val_loss:
  171. best_val_loss = val_loss
  172. best_state = {k: v.clone() for k, v in model.state_dict().items()}
  173. patience_counter = 0
  175. if save_best and checkpoint_dir is not None:
  176. _save_checkpoint(
  177. checkpoint_dir / f"{checkpoint_prefix}_best.pt",
  178. epoch=epoch + 1,
  179. train_loss=avg_loss,
  180. val_loss=val_loss,
  181. is_best=True,
  182. )
  183. else:
  184. patience_counter += 1
  185. if verbose:
  186. print(f" No improvement ({patience_counter}/{patience})")
  187. if patience_counter >= patience:
  188. if verbose:
  189. print(f" Early stopping at epoch {epoch+1}")
  190. break
  191. else:
  192. if verbose:
  193. print(f" Epoch {epoch+1:2d} | train loss: {avg_loss:.6f}")
  195. if avg_loss < best_val_loss:
  196. best_val_loss = avg_loss
  197. best_state = {k: v.detach().cpu().clone() for k, v in model.state_dict().items()}
  199. if save_best and checkpoint_dir is not None:
  200. _save_checkpoint(
  201. checkpoint_dir / f"{checkpoint_prefix}_best.pt",
  202. epoch=epoch + 1,
  203. train_loss=avg_loss,
  204. val_loss=None,
  205. is_best=True,
  206. )
  207. if (
  208. checkpoint_dir is not None
  209. and checkpoint_every is not None
  210. and checkpoint_every > 0
  211. and (epoch + 1) % checkpoint_every == 0
  212. ):
  213. _save_checkpoint(
  214. checkpoint_dir / f"{checkpoint_prefix}_epoch{epoch+1:03d}.pt",
  215. epoch=epoch + 1,
  216. train_loss=avg_loss,
  217. val_loss=val_loss if val_loader is not None else None,
  218. is_best=False,
  219. )
  221. if save_last and checkpoint_dir is not None:
  222. _save_checkpoint(
  223. checkpoint_dir / f"{checkpoint_prefix}_last.pt",
  224. epoch=epoch + 1,
  225. train_loss=avg_loss,
  226. val_loss=val_loss if val_loader is not None else None,
  227. is_best=False,
  228. )
  230. if best_state is not None:
  231. model.load_state_dict(best_state)
  232. if verbose:
  233. print(f" Best val loss: {best_val_loss:.6f} | weights restored")
  235. return model
  238. def get_predictions(model, loader, device=None):
  239. """Runs inference and returns concatenated predictions."""
  241. if device is None:
  242. device = next(model.parameters()).device
  244. model = model.to(device)
  245. model.eval()
  247. all_preds = []
  248. with torch.no_grad():
  249. for batch in loader:
  250. if isinstance(batch, (tuple, list)):
  251. xb = batch[0]
  252. else:
  253. xb = batch
  255. xb = xb.to(device, non_blocking=True)
  257. preds = model(xb)
  258. all_preds.append(preds.detach().cpu().numpy())
  259. return np.concatenate(all_preds, axis=0)
  262. def _is_torch_model(model):
  263. return isinstance(model, nn.Module)
  266. def _flatten_model_inputs(X):
  267. """Convert windowed inputs (N, M, ROI) into tabular features for sklearn models."""
  268. return X.reshape(X.shape[0], -1)
  271. def _flatten_model_targets(Y):
  272. """Convert forecasting targets into 2D multi-output targets for sklearn models."""
  273. return Y.reshape(Y.shape[0], -1)
  276. def _reshape_predictions(preds, target_shape):
  277. """Restore flattened sklearn predictions back to forecasting shape."""
  278. preds = np.asarray(preds, dtype=np.float32)
  280. if preds.shape == target_shape:
  281. return preds
  283. if preds.ndim == 1:
  284. preds = preds[:, None]
  286. if preds.ndim == 2:
  287. return preds.reshape(target_shape)
  289. raise ValueError(
  290. f"Could not reshape predictions from {preds.shape} to {target_shape}"
  291. )
  294. def _expects_windowed_input(model):
  295. return getattr(model, "expects_windowed_input", False)
  298. def _clone_model(model):
  299. """Best-effort clone for either PyTorch modules or sklearn estimators."""
  300. if _is_torch_model(model):
  301. return {k: v.detach().cpu().clone() for k, v in model.state_dict().items()}
  302. return copy.deepcopy(model)
  305. def train_forecasting_model(
  306. model,
  307. X_train,
  308. Y_train,
  309. X_val=None,
  310. Y_val=None,
  311. batch_size=512,
  312. num_epochs=30,
  313. device=None,
  314. patience=5,
  315. checkpoint_dir=None,
  316. checkpoint_prefix="forecast_model",
  317. checkpoint_every=None,
  318. save_best=True,
  319. save_last=False,
  320. verbose=True
  321. ):
  322. """
  323. Train either a PyTorch forecasting model or an sklearn-style estimator.
  324. """
  325. if _is_torch_model(model):
  326. train_loader = DataLoader(
  327. FMRIWindowDataset(X_train, Y_train),
  328. batch_size=batch_size,
  329. shuffle=True,
  330. pin_memory=True,
  331. )
  332. val_loader = None
  333. if X_val is not None and Y_val is not None and len(X_val) > 0:
  334. val_loader = DataLoader(
  335. FMRIWindowDataset(X_val, Y_val),
  336. batch_size=batch_size,
  337. shuffle=False,
  338. pin_memory=True,
  339. )
  341. return train_model(
  342. model,
  343. train_loader,
  344. val_loader=val_loader,
  345. num_epochs=num_epochs,
  346. device=device,
  347. patience=patience,
  348. checkpoint_dir=checkpoint_dir,
  349. checkpoint_prefix=checkpoint_prefix,
  350. checkpoint_every=checkpoint_every,
  351. save_best=save_best,
  352. save_last=save_last,
  353. verbose=verbose,
  354. )
  356. if hasattr(model, "fit") and hasattr(model, "predict"):
  357. if _expects_windowed_input(model):
  358. model.fit(X_train, Y_train)
  359. else:
  360. X_train_flat = _flatten_model_inputs(X_train)
  361. Y_train_flat = _flatten_model_targets(Y_train)
  362. model.fit(X_train_flat, Y_train_flat)
  363. return model
  365. raise TypeError(
  366. "Unsupported model type. Expected a torch.nn.Module or an estimator "
  367. "with fit/predict methods."
  368. )
  371. def predict_forecasting_model(model, X, batch_size=512, device=None):
  372. """Run inference for either a PyTorch forecasting model or an sklearn estimator."""
  373. if _is_torch_model(model):
  374. test_loader = DataLoader(
  375. FMRIWindowDataset(X),
  376. batch_size=batch_size,
  377. shuffle=False,
  378. pin_memory=True,
  379. )
  380. return get_predictions(model, test_loader, device)
  382. if hasattr(model, "predict"):
  383. if _expects_windowed_input(model):
  384. preds = model.predict(X)
  385. else:
  386. preds = model.predict(_flatten_model_inputs(X))
  387. return preds
  389. raise TypeError(
  390. "Unsupported model type. Expected a torch.nn.Module or an estimator "
  391. "with a predict method."
  392. )
  395. def _make_model(model_gen, n_roi=None, M=None, H=None):
  396. """
  397. Instantiate model_gen while preserving compatibility with zero-arg notebook
  398. lambdas and newer factories that accept fold dimensions.
  399. """
  400. try:
  401. signature = inspect.signature(model_gen)
  402. except (TypeError, ValueError):
  403. return model_gen()
  405. kwargs = {}
  406. for name in signature.parameters:
  407. if name in {"n_roi", "input_size", "input_dim"} and n_roi is not None:
  408. kwargs[name] = n_roi
  409. elif name in {"M", "window_size"} and M is not None:
  410. kwargs[name] = M
  411. elif name in {"H", "output_horizon", "horizon"} and H is not None:
  412. kwargs[name] = H
  414. if kwargs:
  415. return model_gen(**kwargs)
  416. return model_gen()
  419. def horizon_rmse(y_true, y_pred):
  420. """Compute and print RMSE separately for each forecast step."""
  421. horizon_scores = []
  422. print("\nHorizon-wise RMSE:")
  423. for h in range(y_true.shape[1]):
  424. r = compute_rmse(y_true[:, h, :], y_pred[:, h, :])
  425. horizon_scores.append(r)
  426. print(f" Step {h+1} RMSE: {r:.6f}")
  427. return horizon_scores
  430. def run_loso_cv(dataset_raw, model_gen, M=50, H=3, stride=1,
  431. num_epochs=20, batch_size=512, device=device,
  432. checkpoint_dir=None, checkpoint_prefix="forecast_model",
  433. checkpoint_every=None, save_best=True, save_last=False,
  434. results_path="loso_results.csv", patience=5,
  435. compute_naive_rmse=True):
  436. """
  437. Leave-One-Subject-Out Cross Validation (LOSO-CV).
  439. Supports both:
  440. - PyTorch forecasting models with the existing training loop
  441. - sklearn-style estimators exposing fit(X, y) and predict(X)
  443. Set ``compute_naive_rmse=False`` when a caller computes the last-value
  444. baseline once and joins those scores onto multiple model result tables.
  446. Returns:
  447. df - LOSO summary dataframe
  448. last_trained_model - model from the last fold
  449. last_X_test - test windows from the last fold
  450. last_Y_test - test targets from the last fold
  451. best_model - model with highest eta across all folds
  452. best_X_test - test windows of the best eta fold
  453. best_Y_test - test targets of the best eta fold
  454. """
  455. if checkpoint_dir is not None:
  456. checkpoint_dir = Path(checkpoint_dir)
  457. checkpoint_dir.mkdir(parents=True, exist_ok=True)
  459. subjects = sorted(set(d["subject"] for d in dataset_raw))
  460. n_subjects = len(subjects)
  462. print(f"\n{'='*60}")
  463. print(f"LOSO-CV | {n_subjects} subjects | M={M}, H={H}")
  464. print(f"{'='*60}")
  466. fold_results = []
  468. last_trained_model = None
  469. last_X_test = None
  470. last_Y_test = None
  472. best_eta_model = None
  473. best_eta_score = -float("inf")
  474. best_eta_subject = None
  475. best_X_test = None
  476. best_Y_test = None
  477. best_n_roi = None
  479. for fold_i, test_subj in enumerate(tqdm(subjects, desc="LOSO Folds")):
  480. print(f"\nFold {fold_i+1}/{n_subjects} | Test subject: {test_subj}")
  482. train_items, test_items = split_by_subject(
  483. dataset_raw,
  484. test_subjects=[test_subj]
  485. )
  487. train_norm = normalize_items(train_items)
  488. test_norm = normalize_items(test_items)
  490. X_tr, Y_tr = build_sliding_windows(train_norm, M, H, stride)
  491. X_te, Y_te = build_sliding_windows(test_norm, M, H, stride)
  493. if len(X_tr) == 0 or len(X_te) == 0:
  494. print("Skipping fold (no valid windows)")
  495. continue
  497. val_split = int(len(X_tr) * 0.9)
  498. X_val, Y_val = X_tr[val_split:], Y_tr[val_split:]
  499. X_tr, Y_tr = X_tr[:val_split], Y_tr[:val_split]
  501. n_roi = X_tr.shape[2]
  502. print(f"Train windows: {len(X_tr)} | Val windows: {len(X_val)} "
  503. f"| Test windows: {len(X_te)} | ROIs: {n_roi}")
  505. try:
  506. model = _make_model(model_gen, n_roi=n_roi, M=M, H=H)
  507. if _is_torch_model(model):
  508. model = model.to(device)
  509. except AttributeError as e:
  510. print(f"Model is not CUDA compatible: {e}\nContinuing with CPU...")
  511. model = _make_model(model_gen, n_roi=n_roi, M=M, H=H)
  513. fold_checkpoint_dir = None
  514. if checkpoint_dir is not None and _is_torch_model(model):
  515. fold_checkpoint_dir = checkpoint_dir
  517. print(f"Training (max {num_epochs} epochs, early stopping patience={patience})...")
  518. model = train_forecasting_model(
  519. model,
  520. X_tr,
  521. Y_tr,
  522. X_val=X_val,
  523. Y_val=Y_val,
  524. batch_size=batch_size,
  525. num_epochs=num_epochs,
  526. device=device,
  527. patience=patience,
  528. checkpoint_dir=fold_checkpoint_dir,
  529. checkpoint_prefix=f"{checkpoint_prefix}_fold{fold_i+1:02d}_{test_subj}",
  530. checkpoint_every=checkpoint_every,
  531. save_best=save_best,
  532. save_last=save_last,
  533. )
  535. all_preds = predict_forecasting_model(
  536. model,
  537. X_te,
  538. batch_size=batch_size,
  539. device=device,
  540. )
  541. all_preds = _reshape_predictions(all_preds, Y_te.shape)
  542. all_targets = Y_te
  544. model_r = compute_rmse(all_targets, all_preds)
  545. model_rmsse = compute_rmsse(all_targets, all_preds)
  546. eta = compute_eta_gauss(all_targets, all_preds)
  548. naive_r = np.nan
  549. beat_naive = pd.NA
  550. if compute_naive_rmse:
  551. naive_model = last_value_model_generator(H=H)
  552. naive_model = train_forecasting_model(
  553. naive_model,
  554. X_tr,
  555. Y_tr,
  556. verbose=False,
  557. )
  558. naive_preds = predict_forecasting_model(
  559. naive_model,
  560. X_te,
  561. batch_size=batch_size,
  562. device=device,
  563. )
  564. naive_preds = _reshape_predictions(naive_preds, Y_te.shape)
  565. naive_r = compute_rmse(all_targets, naive_preds)
  566. beat_naive = model_r < naive_r
  568. print(f"\nResults:")
  569. print(f" MODEL RMSE : {model_r:.6f}")
  570. if compute_naive_rmse:
  571. print(f" Naive RMSE : {naive_r:.6f}")
  572. print(f" eta : {eta:.4f}")
  573. if compute_naive_rmse:
  574. print(f" Beat naive : {'YES' if beat_naive else 'NO'}")
  576. hor_rmse = horizon_rmse(all_targets, all_preds)
  578. fold_results.append({
  579. "test_subject": test_subj,
  580. "Model_RMSE": round(model_r, 6),
  581. "Naive_RMSE": round(naive_r, 6),
  582. "Model_RMSSE": round(model_rmsse, 6),
  583. "eta": round(eta, 4),
  584. "beat_naive": beat_naive,
  585. })
  587. last_trained_model = model
  588. last_X_test = X_te
  589. last_Y_test = Y_te
  591. if eta > best_eta_score:
  592. best_eta_score = eta
  593. best_eta_subject = test_subj
  594. best_eta_model = _clone_model(model)
  595. best_X_test = X_te.copy()
  596. best_Y_test = Y_te.copy()
  597. best_n_roi = n_roi
  598. print(f" New best eta model saved: {test_subj} (eta={eta:.4f})")
  600. del X_tr, Y_tr, X_val, Y_val
  601. del train_norm, test_norm
  602. torch.cuda.empty_cache()
  604. print(f"\n{'='*60}")
  605. print("LOSO-CV SUMMARY")
  606. print(f"{'='*60}")
  608. df = pd.DataFrame(fold_results)
  609. if df.empty:
  610. raise ValueError(
  611. "LOSO-CV produced no valid folds. Check subject count and window "
  612. f"settings M={M}, H={H}, stride={stride}."
  613. )
  614. print(df.to_string(index=False))
  615. print(f"\nMean Model RMSE : {df['Model_RMSE'].mean():.6f}")
  616. if compute_naive_rmse:
  617. print(f"Mean Naive RMSE : {df['Naive_RMSE'].mean():.6f}")
  618. print(f"Mean eta : {df['eta'].mean():.4f}")
  619. print(f"Mean RMSSE : {df['Model_RMSSE'].mean():.6f}")
  620. if compute_naive_rmse:
  621. print(f"Folds beat naive : {df['beat_naive'].sum()} / {len(df)}")
  623. if results_path is not None:
  624. results_path = Path(results_path)
  625. results_path.parent.mkdir(parents=True, exist_ok=True)
  626. df.to_csv(results_path, index=False)
  627. print(f"\nResults saved to {results_path}")
  629. print(f"\nBest eta fold : {best_eta_subject} (eta={best_eta_score:.4f})")
  631. best_model = _make_model(model_gen, n_roi=best_n_roi, M=M, H=H)
  632. if _is_torch_model(best_model):
  633. best_model = best_model.to(device)
  634. if best_eta_model is not None:
  635. best_model.load_state_dict(best_eta_model)
  636. best_model.eval()
  637. elif best_eta_model is not None:
  638. best_model = best_eta_model
  640. return df, last_trained_model, last_X_test, last_Y_test, \
  641. best_model, best_X_test, best_Y_test

training.py at commit a064b79, under Apache-2.0 · at the source

Overview

Authors: Arda Bayer1, Zhiyao Zhang1, Ahmet Emre Ipek2, Rose Khavari3, Behnaam Aazhang1
ORCID iDs: Arda Bayer, Zhiyao Zhang, Ahmet Emre Ipek, Rose Khavari, Behnaam Aazhang
  1. Department of Electrical & Computer Engineering, Rice University, Houston, TX 77005, USA; (Z.Z.); (B.A.)
  2. Department of Electrical & Electronics Engineering, Özyeğin University, 34794 Istanbul, Türkiye
  3. Department of Urology, Houston Methodist, Houston, TX 77030, USA
Institutions: Rice University (United States); Özyeğin University (Türkiye); Houston Methodist (United States)
Journal: Entropy (Basel, Switzerland), volume 28, issue 7, article 738
Dates: received 22 May 2026; accepted 22 June 2026; published online 1 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.3390/e28070738 · PMID 42511348 · PMCID PMC13409728 · OpenAlex W7166854284
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: fMRI (modality), computational (subfield)
Methods: Statistics, Machine learning, Preprocessing, Connectivity, fMRI & imaging
Keywords: functional magnetic resonance imaging, BOLD signal, probabilistic forecasting, information theory, entropy, directed information, brain dynamics, stochastic processes, transformer models, recurrent neural networks
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Rice University (Provost's TMC Collaborator Seed Fund)
Citations: not cited yet (Europe PMC); 32 references in the paper

Abstract

Functional magnetic resonance imaging (fMRI) signals exhibit complex temporal structure arising from multivariate neural dynamics, physiological variability, and measurement uncertainty. In this work, we formulate region-of-interest-level fMRI analysis as a probabilistic multi-step forecasting problem and investigate the predictability of blood-oxygen-level-dependent (BOLD) activity from an information-theoretic perspective. Using the Natural Scenes Dataset, we model multiregional BOLD activity as a stochastic process with finite memory and train multiple forecasting architectures, including linear regression, exponential smoothing, recurrent neural networks, and transformer-based models, to predict future BOLD samples from preceding temporal observations. Forecasting performance is analyzed together with entropy-based quantities, including marginal entropy, conditional entropy, and normalized predictive information measures estimated directly from model-derived predictive distributions without imposing restrictive Gaussian assumptions on the underlying BOLD dynamics. The transformer model achieved significant improvement over a naive persistence baseline (p=0.001) while yielding a high predictive information fraction (η=75.49%). Post hoc directed information analysis revealed that short-horizon prediction was dominated primarily by autoregressive, within-ROI, temporal structure. Overall, the proposed framework demonstrates how probabilistic forecasting and information-theoretic analysis can be integrated to characterize the predictability, uncertainty structure, and directional organization of large-scale fMRI dynamics and may support future downstream neuroengineering and neural-state inference applications.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 20 matches between paragraphs and lines of code.

ab126/fmri_forecasting

License: Apache-2.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: a064b79113462d62a1c3808f3e3008fe8cf8df31, 15 September 2026
Languages: Python (22), Jupyter (6)
Size: 31 files, 28 scripts
Software Heritage: not archived
Found in: “Data Availability Statement”
Holds: README, license file, tests, 5 notebooks
Not found: CITATION.cff, environment file, continuous integration, documentation
Tools: NumPy (20 files), PyTorch (13 files), Matplotlib (6 files), scikit-learn (4 files), NiBabel (3 files), pandas (3 files), seaborn (3 files), Nilearn (2 files), ANTs (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
30 files
At the source: github.com/ab126/fmri_forecasting

Zenodo 20341604

License: Apache-2.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Data Availability Statement”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (19 files), PyTorch (13 files), Matplotlib (6 files), scikit-learn (4 files), pandas (3 files), seaborn (3 files), NiBabel (2 files), Nilearn (2 files), ANTs (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
29 files
At the source: doi.org/10.5281/zenodo.20341604

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 55 scripts, each with its path and the digest of its content;
  • 20 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability Statement

The Natural Scenes Dataset (NSD) analyzed in this study is publicly available from the original NSD release [17]. The code used for forecasting, entropy estimation, and directed information analysis is publicly available at https://github.com/ab126/fmri_forecasting (accessed on 21 June 2026). An archived release of the code associated with this manuscript is available through Zenodo: https://doi.org/10.5281/zenodo.20341604.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 10 keywords, 1 funder, 21 references.

Cite

This paper

Bayer, A., Zhang, Z., Ipek, A. E., Khavari, R., & Aazhang, B. (2026). Probabilistic Forecasting and Information-Theoretic Analysis of Multivariate fMRI Dynamics. Entropy (Basel, Switzerland), 28(7), 738. https://doi.org/10.3390/e28070738

BibTeX

@article{bayer2026probabilistic,
author = {Bayer, Arda and Zhang, Zhiyao and Ipek, Ahmet Emre and Khavari, Rose and Aazhang, Behnaam},
title = {{Probabilistic Forecasting and Information-Theoretic Analysis of Multivariate fMRI Dynamics}},
journal = {Entropy (Basel, Switzerland)},
year = {2026},
month = jul,
volume = {28},
number = {7},
pages = {738},
publisher = {Multidisciplinary Digital Publishing Institute (MDPI)},
issn = {1099-4300},
doi = {10.3390/e28070738},
url = {https://doi.org/10.3390/e28070738},
pmid = {42511348},
pmcid = {PMC13409728}
}

RIS

TY - JOUR
AU - Bayer, Arda
AU - Zhang, Zhiyao
AU - Ipek, Ahmet Emre
AU - Khavari, Rose
AU - Aazhang, Behnaam
TI - Probabilistic Forecasting and Information-Theoretic Analysis of Multivariate fMRI Dynamics
T2 - Entropy (Basel, Switzerland)
J2 - Entropy (Basel)
PY - 2026
DA - 2026/07/01
VL - 28
IS - 7
SP - 738
SN - 1099-4300
PB - Multidisciplinary Digital Publishing Institute (MDPI)
DO - 10.3390/e28070738
UR - https://doi.org/10.3390/e28070738
LA - en
ER -

CSL-JSON

{
"id": "10.3390/e28070738",
"type": "article-journal",
"title": "Probabilistic Forecasting and Information-Theoretic Analysis of Multivariate fMRI Dynamics",
"container-title": "Entropy (Basel, Switzerland)",
"author": [
{
"family": "Bayer",
"given": "Arda"
},
{
"family": "Zhang",
"given": "Zhiyao"
},
{
"family": "Ipek",
"given": "Ahmet Emre"
},
{
"family": "Khavari",
"given": "Rose"
},
{
"family": "Aazhang",
"given": "Behnaam"
}
],
"container-title-short": "Entropy (Basel)",
"volume": "28",
"issue": "7",
"page": "738",
"DOI": "10.3390/e28070738",
"PMID": "42511348",
"PMCID": "PMC13409728",
"ISSN": "1099-4300",
"publisher": "Multidisciplinary Digital Publishing Institute (MDPI)",
"URL": "https://doi.org/10.3390/e28070738",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
1
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.7554/elife.107933 [code]
Modality-agnostic decoding of vision and language from fMRI.
Journal: eLife
In common: Nilearn, statsmodels, NiBabel, 6 other tools, fMRI, 3 references
[2] doi:10.1038/s41597-026-07248-6 [code]
A large-scale fMRI dataset for vision-language semantic association.
Journal: Scientific data
In common: ANTs, Nilearn, NiBabel, 5 other tools, fMRI, 2 references
[3] doi:10.1162/netn.a.547 [code]
An evaluation of the efficacy of single-echo and multi-echo fMRI denoising strategies.
Journal: Network neuroscience (Cambridge, Mass.)
In common: ANTs, Nilearn, statsmodels, 6 other tools, fMRI, 1 reference
[4] doi:10.1016/j.neuroimage.2026.122171 [code]
A conserved node degree-based backbone and flexible hub organization of brain connectome during naturalistic movie watching.
Journal: NeuroImage
In common: ANTs, Nilearn, NiBabel, 4 other tools, fMRI, 2 references
[5] doi:10.1162/imag.a.1286 [code]
Behavioral imitation with artificial neural networks leads to personalized models of brain dynamics during videogame play.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Nilearn, statsmodels, NiBabel, 6 other tools, fMRI, 1 reference
[6] doi:10.1162/imag.a.1256 [code]
Gamer in the scanner: Event-related analysis of fMRI activity during retro videogame play guided by automated annotations of game content.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Nilearn, statsmodels, NiBabel, 6 other tools, fMRI, 1 reference
[7] doi:10.1038/s41467-026-76452-0 [code]
Music evokes shared neural representations of imagined narratives across sensory modalities.
Journal: Nature communications
In common: ANTs, Nilearn, statsmodels, 7 other tools
[8] doi:10.1371/journal.pbio.3003856 [code]
Aging and metabolism contribute separately to brain-body health.
Journal: PLoS biology
In common: ANTs, Nilearn, statsmodels, 7 other tools
[9] doi:10.1038/s41467-026-71151-2 [code]
Common and distinct neural correlates of social interaction processing and theory of mind in narratives.
Journal: Nature communications
In common: ANTs, Nilearn, statsmodels, 6 other tools, 1 reference
[10] doi:10.1162/imag.a.1252 [code]
Does the brain's E:I balance really shape long-range temporal correlations? Lessons learned from 3T MRI.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: ANTs, NiBabel, seaborn, 4 other tools, fMRI, 2 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.