cortical-motor-engram-decoding / 03_hyperalignment / hyperalignment.py
hyperalignment.py
Raw
#!/usr/bin/env python
# coding: utf-8


# # Hyperalignment

# In[ ]:


"""
Whole-brain searchlight hyperalignment

This script performs searchlight-based hyperalignment using Neuroboros
and Procrustes alignment on surface data. For each subject, data from all
other subjects are used to build the common space, and the left-out subject is
aligned to this common space (leave-one-subject-out scheme).

Outputs:
- Aligned training data
- Aligned test data
- Corresponding labels (sequences)

Author: Theo Marins
"""

# ============================================================================
# Imports
# ============================================================================

import os
import pickle
import warnings
from time import time

import numpy as np
import neuroboros as nb
import nibabel as nib
from tqdm import tqdm
from sklearn.exceptions import ConvergenceWarning

from hyperalignment.searchlight import (
    searchlight_procrustes,
    searchlight_template
)

# Silence convergence warnings (expected in ridge/procrustes)
warnings.filterwarnings("ignore", category=ConvergenceWarning)

# ============================================================================
# Paths
# ============================================================================


subject_data_file = os.path.join(PKL_DIR, "darrays.pkl")
run_label_file = os.path.join(PKL_DIR, "run_label.pkl")
subject_label_file = os.path.join(PKL_DIR, "subject_label.pkl")
session_label_file = os.path.join(PKL_DIR, "session_label.pkl")
sequence_label_file = os.path.join(PKL_DIR, "sequence_label.pkl")

# ============================================================================
# Load data
# ============================================================================

with open(subject_data_file, "rb") as f:
    subject_data = pickle.load(f)

with open(run_label_file, "rb") as f:
    run_label = pickle.load(f)

with open(subject_label_file, "rb") as f:
    subject_label = pickle.load(f)

with open(session_label_file, "rb") as f:
    session_label = pickle.load(f)

with open(sequence_label_file, "rb") as f:
    sequence_label = pickle.load(f)

print("Loaded data for sessions:", subject_data.keys())

# ============================================================================
# Parameters
# ============================================================================

sessions = ["01"]
subjects = list(subject_data["01"].keys())
runs = list(subject_data["01"]["01"].keys())

radius = 20
space = "fsavg-ico32"

# ============================================================================
# Masks
# ============================================================================

mask_l = nb.mask("l", space=space, mask=True)
mask_r = nb.mask("r", space=space, mask=True)

# ============================================================================
# Output containers
# ============================================================================

train_data_dict = {}
test_data_dict = {}
train_label_dict = {}
test_label_dict = {}

# ============================================================================
# Hyperalignment helper
# ============================================================================

def hyperalign_train_data(train_data, sls, dists, radius=20):
    """
    Perform searchlight hyperalignment on training subjects.

    Parameters
    ----------
    train_data : array
        Shape (n_subjects, n_samples, n_vertices)
    sls, dists : searchlight definitions
    radius : int
        Searchlight radius

    Returns
    -------
    template : array
        Group template
    aligned_train : array
        Hyperaligned training data
    """

    print("  Building template...")
    template = searchlight_template(
        dss=train_data,
        sls=sls,
        dists=dists,
        radius=radius,
        tmpl_kind="pca",
        n_jobs=1
    )

    aligned_train = []

    for i, subj_data in enumerate(train_data):
        print(f"  Aligning training subject {i}")

        W = searchlight_procrustes(
            subj_data,
            template,
            sls,
            dists,
            radius=radius
        )

        aligned_train.append(subj_data @ W)

    aligned_train = np.stack(aligned_train, axis=0)
    return template, aligned_train

# ============================================================================
# Main loop
# ============================================================================

start_time = time()

for hemi_name, hemi_key in zip(["left", "right"], ["l", "r"]):

    print("\n" + "=" * 60)
    print(f"Processing {hemi_name.upper()} hemisphere")
    print("=" * 60)

    mask = mask_l if hemi_key == "l" else mask_r

    # Define searchlights
    sls, dists = nb.sls(
        hemi_key,
        radius,
        space=space,
        mask=True,
        return_dists=True
    )

    train_data_dict[hemi_key] = {}
    test_data_dict[hemi_key] = {}
    train_label_dict[hemi_key] = {}
    test_label_dict[hemi_key] = {}

    # Leave-one-subject-out
    for test_subject in subjects:

        print(f"\nTest subject: {test_subject}")

        train_data = []
        train_labels = []
        test_labels = []

        # ------------------------------------------------------------
        # Assemble data
        # ------------------------------------------------------------
        for ses in sessions:
            for subject in subjects:

                subj_runs = []
                subj_labels = []

                for run in runs:
                    X = subject_data[ses][subject][run][hemi_name][:, mask] 
                    y = sequence_label[ses][subject][run]

                    subj_runs.append(X) #
                    subj_labels.extend(y)

                subj_data = np.concatenate(subj_runs, axis=0)

                if subject == test_subject:
                    test_data = subj_data
                    test_labels = subj_labels
                else:
                    train_data.append(subj_data)
                    train_labels.extend(subj_labels)

        train_data = np.stack(train_data, axis=0)

        # ------------------------------------------------------------
        # Hyperalignment (train subjects)
        # ------------------------------------------------------------
        print("Running hyperalignment...")
        template, train_data_ha = hyperalign_train_data(
            train_data, sls, dists, radius
        )

        # ------------------------------------------------------------
        # Align test subject
        # ------------------------------------------------------------
        W_test = searchlight_procrustes(
            test_data,
            template,
            sls,
            dists,
            radius=radius
        )

        test_data_ha = test_data @ W_test

        # ------------------------------------------------------------
        # Store
        # ------------------------------------------------------------
        train_data_dict[hemi_key][test_subject] = train_data_ha
        test_data_dict[hemi_key][test_subject] = test_data_ha
        train_label_dict[hemi_key][test_subject] = train_labels
        test_label_dict[hemi_key][test_subject] = test_labels

print(f"\nCompleted in {time() - start_time:.2f} seconds.")

# ============================================================================
# Save results 
# ============================================================================

with open(os.path.join(PKL_DIR, "train_data.pkl"), "wb") as f:
    pickle.dump(train_data_dict, f)

with open(os.path.join(PKL_DIR, "test_data.pkl"), "wb") as f:
    pickle.dump(test_data_dict, f)

with open(os.path.join(PKL_DIR, "train_labels.pkl"), "wb") as f:
    pickle.dump(train_label_dict, f)

with open(os.path.join(PKL_DIR, "test_labels.pkl"), "wb") as f:
    pickle.dump(test_label_dict, f)

print(f"\nSaved hyperaligned data to: {PKL_DIR}")