Rotation-based metric on the Riemannian manifold of SPD matrices with applications to source data selection for brain-computer interface transfer learning.

Code ↔ Paper

3 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 3 matches

[1] § Results › Pole ratio for BCI data ↔ rotation_analysis.py, lines 1266–1342 · score 0.62 · circle represents, intra subject accuracy, south pole, pole ratio, north, class
[2] § Materials and methods › BCI data ↔ rotation_analysis.py, lines 171–245 · score 0.61 · motor imagery, covariance matrix, filtered, channels, MOABB, 35 Hz
[3] § Materials and methods › Data split for classification evaluation ↔ rotation_analysis.py, lines 436–585 · score 0.51 · target training, splits, MDM, transfer learning, classification, transformed

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 paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 1,746 lines · 59 KB · MIT · 3 matches

#!python
# venv: rotation
# Imports
# General Libraries
import os
import sys
import time
import math
import random
import datetime
import subprocess
import copy
import warnings
warnings.filterwarnings("ignore")
warnings.simplefilter(action="ignore", category=FutureWarning) # Ignore futurewarning warnings.
import shutil
from pathlib import Path
# Numpy Library
import numpy as np
# Making numpy single threaded, so I can multithread this program myself:
os.environ["OMP_NUM_THREADS"] = "1" # export OMP_NUM_THREADS=4
os.environ["OPENBLAS_NUM_THREADS"] = "1" # export OPENBLAS_NUM_THREADS=4
os.environ["MKL_NUM_THREADS"] = "1" # export MKL_NUM_THREADS=6
os.environ["VECLIB_MAXIMUM_THREADS"] = "1" # export VECLIB_MAXIMUM_THREADS=4
os.environ["NUMEXPR_NUM_THREADS"] = "1" # export NUMEXPR_NUM_THREADS=6
# Pandas Library
# import pandas as pd
# Pyrieman library
from pyriemann.classification import MDM
from pyriemann.estimation import Covariances
from pyriemann.transfer import TLCenter,TLScale,TLRotate
from pyriemann.utils.distance import distance
from pyriemann.utils.mean import mean_covariance
from pyriemann.tangentspace import TangentSpace
from pyriemann.clustering import Kmeans
# Scikit-learn library
from sklearn.metrics import accuracy_score
from sklearn.model_selection import StratifiedShuffleSplit
from sklearn.pipeline import make_pipeline
from sklearn import svm
# Matplotlib Library
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.patches import CirclePolygon
# Moabb
from moabb.datasets import BNCI2014_001,PhysionetMI, Dreyer2023,Stieger2021,Liu2024,Lee2019_MI,Cho2017
from moabb.paradigms import MotorImagery
# Multiprocessing
from multiprocessing import Pool
# Scipy
from scipy.optimize import curve_fit
from scipy.stats import ttest_ind
# Other
from PIL import Image
mpl.rcParams['font.family'] = "serif"
plt.rcParams.update({'font.size': 15})
random_state = 42
rng = np.random.default_rng(seed=random_state)
# GLOBAL PARAMETERS
NBR_THREADS = 1 # Nbr of threads for parallellisation
DATASETS = [
(BNCI2014_001(),['feet','tongue'],'BCI-IV_feet_tongue'),
(PhysionetMI(),['feet','hands'],'Physio_feet_hands'),
(Dreyer2023(),['right_hand','left_hand'],'Dreyer_left_right'),
(BNCI2014_001(),['right_hand','left_hand'],'BCI-IV_left_right'),
(PhysionetMI(),['right_hand','left_hand'],'Physio_left_right'),
(Lee2019_MI(),['right_hand','left_hand'],'Lee_left_right'),
(Cho2017(),['right_hand','left_hand'],'Cho_left_right'),
]
PROCESSING_DICT={
'channels': ['FC3', 'FC4', 'C5', 'C3', 'C1', 'C2', 'C4', 'C6', 'CP3', 'CP4'],
'fmin':7,
'fmax':35,
'tmin':1,
'tmax':2,
'baseline':None,
}
print(__doc__)
class Logger(object):
"""For saving everything printed to file as well."""
def __init__(self, filename="_console_log.txt"):
"""Initialize the class."""
self.terminal = sys.stdout
self.log = open(filename, "a")
def write(self, message):
"""Write both to stdout and to file."""
self.terminal.write(message)
self.log.write(message)
def flush(self):
"""We don't do flush."""
pass
def print_git_version():
"""We want to know what version of the source code we're dealing with."""
# Ideally, commit to git _every time_ before you run the code.
try:
git_output = subprocess.check_output("git log -n 1", shell=True, stderr=subprocess.STDOUT)
print("\n\n### git log -n 1\n%s\n\n" % git_output.decode("utf-8"))
git_output2 = subprocess.check_output("git status", shell=True, stderr=subprocess.STDOUT)
print("### git status\n%s\n\n" % git_output2.decode("utf-8"))
git_output3 = subprocess.check_output("git diff", shell=True, stderr=subprocess.STDOUT)
print("### git diff\n%s\n\n\n" % git_output3.decode("utf-8"))
except:
print("\n\n### WARNING: The code you're running is NOT under GIT version control. You can do better. Behave.\n\n")
def save_source_code(run_logdir):
print("\n\n###\n### Saving the python source code used, from the file '%s':\n###" % __file__)
shutil.copy2(__file__, f'{run_logdir}/_script.py')
def write_description(file_name='_description'):
with open(f'{run_logdir}/{file_name}.txt', 'a') as f:
# Write some text to the file
f.write("This is a description file.\n")
f.write(f'NBR_THREADS: {NBR_THREADS}\n')
for key in PROCESSING_DICT.keys():
f.write(f'{key}: {PROCESSING_DICT[key]}\n')
for (dataset,events,identifier) in DATASETS:
f.write(f'\nDataset: {dataset.code}')
f.write(f'\nClasses: {"_".join(events)}')
f.write(f'\nidentifier: {"_".join(identifier)}\n')
# -----------------------------------------------
# -----------------------------------------------
# LOADING DATA : MOABB
# -----------------------------------------------
def get_subject_dataset_MOABB(subject,dataset,events,processing_dict,run_logdir='.',data_dir='.'):
'''
subject: int
subject number
dataset: object with MOABB dataset
Dataset object.
event: list of str
Event names in the dataset to be used.
'''
print(f'\nLoading data subj: {subject}')
folder_name = f'{data_dir}/subj_{subject}'
os.makedirs(folder_name, exist_ok=True)
full_file_path_keep = os.path.join(folder_name, 'keep.npy')
full_file_path_covs = os.path.join(folder_name, 'covariance_matrices.npy')
full_file_path_str_labels = os.path.join(folder_name, 'str_labels.npy')
if os.path.exists(full_file_path_keep) and os.path.exists(full_file_path_covs) and os.path.exists(full_file_path_str_labels):
print(f'Loading stored data for subjects')
# Load covs
covs = np.load(full_file_path_covs,allow_pickle=False)
print(f'covs loaded from {full_file_path_covs}')
# Load str_labels
str_labels_filtered = np.load(full_file_path_str_labels,allow_pickle=False)
print(f'str_labels loaded from {full_file_path_str_labels}')
keep = np.load(full_file_path_keep,allow_pickle=False)
print(f'keep loaded from {full_file_path_keep}')
if not keep: # Print diagnostics to file to know why excluded.
for cla in np.unique(str_labels_filtered):
check_spd_matrices(covs[str_labels_filtered==cla], name=f'Subject: {subject} class: {cla}',run_logdir=run_logdir)
else:
keep = True # placeholder
paradigm = MotorImagery(
fmin=processing_dict['fmin'],#7,
fmax=processing_dict['fmax'],#35,
tmin=processing_dict['tmin'],#1.,
tmax=processing_dict['tmax'],#2.,
baseline=processing_dict['baseline'],#None,
channels=processing_dict['channels'])#CHANNELS) # Here is preprocessing settings to be added. Also select channels.
# https://moabb.neurotechx.com/docs/generated/moabb.paradigms.MotorImagery.html#moabb.paradigms.MotorImagery
X, str_labels, metadata = paradigm.get_data(dataset=dataset, subjects=[subject])
# Filter out data with wrong labels.
mask = np.isin(str_labels, events) # boolean mask
X_filtered = X[mask]
str_labels_filtered = str_labels[mask]
# Compute covariance matrices on scaled data
covs = Covariances().fit_transform(X_filtered)
print(f"\n\n === Diagnostics for Subject: {subject} === ")
for cla in np.unique(str_labels_filtered):
keep_cla = check_spd_matrices(covs[str_labels_filtered==cla], name=f'Subject: {subject} class: {cla}',run_logdir=run_logdir)
if not keep_cla:
keep = False
print()
# Saving data
if not os.path.exists(full_file_path_covs):
np.save(full_file_path_covs,covs,allow_pickle=False)
if not os.path.exists(full_file_path_keep):
np.save(full_file_path_keep,keep,allow_pickle=False)
if not os.path.exists(full_file_path_str_labels):
np.save(full_file_path_str_labels,str_labels_filtered,allow_pickle=False)
return covs, str_labels_filtered, keep
def check_spd_matrices(covs, name="X",run_logdir='.'):
"""
Check symmetry, SPD condition, and condition numbers for a stack of matrices.
Parameters
----------
covs : ndarray, shape (n_matrices, n, n)
Input matrices to check.
name : str, optional
Label for reporting.
"""
n_covs = covs.shape[0]
max_asym = []
min_eigs = []
conds = []
for i, A in enumerate(covs):
# enforce symmetry check
asym = np.max(np.abs(A - A.T))
max_asym.append(asym)
# eigenvalues
w = np.linalg.eigvalsh(0.5 * (A + A.T)) # symmetrize for safety
min_eigs.append(np.min(w))
# condition number (ratio of largest to smallest eigenvalue)
cond = np.max(w) / np.min(w) if np.min(w) > 0 else np.inf
conds.append(cond)
max_asym = np.array(max_asym)
min_eigs = np.array(min_eigs)
conds = np.array(conds)
mean_cov = mean_covariance(covs,metric='riemann')
print(f" --- Diagnostics for {name} --- ")
print(f"Number of matrices: {n_covs}")
print(f'ratio max/min/mean | max: {(mean_cov.max()-mean_cov.mean())/(mean_cov.max()-mean_cov.min()):.2e}, min: {(mean_cov.mean()-mean_cov.min())/(mean_cov.max()-mean_cov.min()):.2e}')
print(f"Min eigenvalue | min: {min_eigs.min():.2e}, mean: {min_eigs.mean():.2e}")
print(f"Condition number | max: {conds.max():.2e}, mean: {conds.mean():.2e}")
with open(f'{run_logdir}/_diagnostics.txt', 'a') as f:
f.write(f"\n\n --- Diagnostics for {name} --- ")
f.write(f"\nNumber of matrices: {n_covs}")
f.write(f'\nratio max/min/mean | max: {(mean_cov.max()-mean_cov.mean())/(mean_cov.max()-mean_cov.min()):.2e}, min: {(mean_cov.mean()-mean_cov.min())/(mean_cov.max()-mean_cov.min()):.2e}')
f.write(f"\nMin eigenvalue | min: {min_eigs.min():.2e}, mean: {min_eigs.mean():.2e}")
f.write(f"\nCondition number | max: {conds.max():.2e}, mean: {conds.mean():.2e}")
ratio_in_matrix = (mean_cov.max()-mean_cov.mean())/(mean_cov.max()-mean_cov.min())
# flag outliers
bad_sym = np.where(max_asym > 1e-10)[0]
bad_spd = np.where(min_eigs <= 0.01)[0]
bad_cond = np.where(conds > 1e8)[0]
keep = True
print()
if ratio_in_matrix > 0.8:
print(f"[!] bad ratio {ratio_in_matrix}.")
keep = False
if bad_sym.size:
print(f"[!] {len(bad_sym)} matrices not symmetric (idx: {bad_sym[:10]})")
if bad_spd.size:
print(f"[!] {len(bad_spd)} matrices bad eig, idx: {bad_spd[:10]})")
keep = False
if bad_cond.size:
print(f"[!] {len(bad_cond)} matrices ill-conditioned (cond > 1e8, idx: {bad_cond[:10]})")
if not (bad_sym.size or bad_spd.size or bad_cond.size or ratio_in_matrix > 0.8):
print("All matrices look fine.")
return keep
# -----------------------------------------------
# Multiprocessing helper function:
# -----------------------------------------------
def worker_function_loading_data(iteration_idx,dataset,subjects,events,processing_dict,run_logdir='.',data_dir='.'):
# This function will be passed to the Pool.
print(f"### worker_function ({subjects[iteration_idx]}/{len(subjects)}) for itr_{iteration_idx}")
covs, str_labels, keep = get_subject_dataset_MOABB(subject=subjects[iteration_idx],dataset=dataset,events=events,processing_dict=processing_dict,run_logdir=run_logdir,data_dir=data_dir)
subject_data = {
'X_cov' : covs,
'y': str_labels,
'y_enc' : np.array([f'subj_{subjects[iteration_idx]}/{str(item)}' for item in str_labels]),
'subject_nbr': subjects[iteration_idx],
'subj_id': f'subj_{subjects[iteration_idx]}'
}
return subject_data,keep
def load_data_MOABB(subjects,dataset,events,processing_dict,run_logdir='.',data_dir_identifier='.'):
all_subject_data = []
excluded_subjects = []
data_dir = f'data/raw_data/{data_dir_identifier}'
if NBR_THREADS == 1:
print('\n-------------------------------')
print( '--> Single thread started < ---')
print( '-------------------------------\n')
for subj in subjects:
covs, str_labels,keep = get_subject_dataset_MOABB(subj,dataset,events,processing_dict,run_logdir,data_dir)
if keep:
subject_data = {
'X_cov' : covs,
'y': str_labels,
'y_enc' : np.array([f'subj_{subj}/{str(item)}' for item in str_labels]),
'subject_nbr': subj,
'subj_id': f'subj_{subj}'
}
all_subject_data.append(subject_data)
else:
excluded_subjects.append(int(subj))
else:
# Multithreaded version:
print('\n---------------------------------')
print( '--> Multiprocessing started < ---')
print( '---------------------------------\n')
pool = Pool(processes=NBR_THREADS)
results = []
for iteration_idx in range(len(subjects)):
res = pool.apply_async(worker_function_loading_data, args=(iteration_idx,dataset,subjects,events,processing_dict,run_logdir,data_dir))
results.append(res)
if res.get()[1]:
all_subject_data.append(res.get()[0])
else:
excluded_subjects.append(int(subjects[iteration_idx]))
print('---> Closing pools\n')
pool.close()
pool.join()
print('---> Pools closed \n')
print("!!!!!!!!!!!!!!!!!!!!!!!!!!!")
print("!!!! EXCLUDED SUBJECTS !!!!")
print(excluded_subjects)
print("!!!!!!!!!!!!!!!!!!!!!!!!!!!")
with open(f'{run_logdir}/_diagnostics.txt', 'a') as f:
f.write(f"\n\n!!!!!!!!!!!!!!!!!!!!!!!!!!!")
f.write(f"\n!!!! EXCLUDED SUBJECTS !!!!")
f.write(f"\n{excluded_subjects}")
f.write(f"\n!!!!!!!!!!!!!!!!!!!!!!!!!!!")
return all_subject_data
# ################################
# ######## ANALYSIS #######
# ################################
def run_one_subject(data,iteration_idx,run_logdir='.',data_dir_identifier='.'):
# get data
X_covs = data[iteration_idx]['X_cov']
y = data[iteration_idx]['y']
subj_id = data[iteration_idx]['subj_id']
subject_nbr = data[iteration_idx]['subject_nbr']
print(f'\nSubj: {subject_nbr}')
find_distance_between_classes_to_pole(data,iteration_idx,run_logdir,data_dir_identifier)
find_transfer_learning_accuracies(data,iteration_idx,run_logdir,data_dir_identifier)
return
# ---- transfer learning
def transfer_learning_one_target_vs_one_source(data, target_idx, source_idx,run_logdir='.',data_dir='.',data_dir_identifier='.'):
X_covs_target = data[target_idx]['X_cov']
y_target = data[target_idx]['y']
y_enc_target = data[target_idx]['y_enc']
subj_id_target = data[target_idx]['subj_id']
subject_nbr_target = data[target_idx]['subject_nbr']
X_covs_source = data[source_idx]['X_cov']
y_source = data[source_idx]['y']
y_enc_source = data[source_idx]['y_enc']
subj_id_source = data[source_idx]['subj_id']
subject_nbr_source = data[source_idx]['subject_nbr']
full_path_accuracies = f"{data_dir}/{subj_id_source}_accuracies.npy"
if os.path.exists(full_path_accuracies):
# If it exists do nothing.
print(f'Accuracies for {subj_id_target} with {subj_id_source} already exist, skipping...')
scores = np.load(full_path_accuracies,allow_pickle=True).item()
return scores
clf = MDM(metric=dict(mean='riemann', distance='riemann'))
splitter = StratifiedShuffleSplit(n_splits=4, train_size=0.75, random_state=42) # Shuffels data.
nothing_score = np.zeros(splitter.get_n_splits())
recenter_score = np.zeros(splitter.get_n_splits())
scale_score = np.zeros(splitter.get_n_splits())
rotate_score = np.zeros(splitter.get_n_splits())
for i, (train_index,test_index) in enumerate(splitter.split(X_covs_target, y_target)):
X_target_train = copy.deepcopy(X_covs_target[train_index])
y_target_train = copy.deepcopy(y_target[train_index])
y_enc_target_train = copy.deepcopy(y_enc_target[train_index])
X_source_train = copy.deepcopy(X_covs_source)
y_source_train = copy.deepcopy(y_source)
y_enc_source_train = copy.deepcopy(y_enc_source)
X_test = copy.deepcopy(X_covs_target[test_index])
y_test = copy.deepcopy(y_target[test_index])
y_enc_test = copy.deepcopy(y_enc_target[test_index])
nbr_target_training_data = len(y_target_train)
if target_idx == source_idx: # Intra subject
clf.fit(X_target_train,y_target_train)
y_pred = clf.predict(X_test)
score = accuracy_score(y_test, y_pred)
nothing_score[i] = score
recenter_score[i] = score
scale_score[i] = score
rotate_score[i] = score
continue # Continue with next split.
# Format data in correct way for transfer learning
X_train = np.concatenate((X_source_train,X_target_train))
y_enc = np.concatenate((y_enc_source_train,y_enc_target_train))
y_train = np.concatenate((y_source_train,y_target_train))
# Before any transfer learning
clf.fit(X_train[:(-nbr_target_training_data)],y_train[:(-nbr_target_training_data)])
y_pred = clf.predict(X_test)
nothing_score[i] = accuracy_score(y_test, y_pred)
# Transfer learnig steps:
tl_recenter = TLCenter(target_domain=subj_id_target,metric='riemann')
tl_scale = TLScale(target_domain=subj_id_target, final_dispersion=1.0, centered_data=True, metric='riemann')
tl_rotate = TLRotate(target_domain=subj_id_target, metric='riemann',tol_step=1e-9, maxiter=10000)
# === RECENTER ===
# Recenter the data.
X_train = tl_recenter.fit_transform(X_train, y_enc)
X_test = tl_recenter.transform(X_test)
clf.fit(X_train[:(-nbr_target_training_data)],y_train[:(-nbr_target_training_data)])
y_pred = clf.predict(X_test)
recenter_score[i] = accuracy_score(y_test, y_pred)
# === SCALE ===
X_train = tl_scale.fit_transform(X_train, y_enc)
X_test = tl_scale.transform(X_test)
clf.fit(X_train[:(-nbr_target_training_data)],y_train[:(-nbr_target_training_data)])
y_pred = clf.predict(X_test)
scale_score[i] = accuracy_score(y_test, y_pred)
# === ROTATE ===
with warnings.catch_warnings(record=True) as w:
warnings.simplefilter("always") # capture all warnings
str_test_data = '_'.join(test_index.astype(str))
rotations_path = f"data/rotations/{data_dir_identifier}/{subj_id_target}/{subj_id_source}_iter_{i}.npy"
test_data_path = f"data/rotations/{data_dir_identifier}/{subj_id_target}/{subj_id_source}_iter_{i}_test_data.npy"
convergece_file = f"data/rotations/{data_dir_identifier}/{subj_id_target}/convergence_info.txt"
if os.path.exists(rotations_path) and os.path.exists(test_data_path):
test_str = np.load(test_data_path,allow_pickle=False)
if test_str == str_test_data:
load_rotation = True
else:
bug # You have to remove rotational data since the test/target split are wrong.
load_rotation = False
else:
load_rotation = False
if load_rotation:
Q_rot = np.load(rotations_path,allow_pickle=False)
X_to_rotate = X_train[:-nbr_target_training_data] # The first samples in the training vector are source data. Only source data is rotated
X_rotated = Q_rot @ X_to_rotate @ Q_rot.T
X_train[:(-nbr_target_training_data)] = X_rotated
print(f'Rotation loaded from {rotations_path}')
else:
X_train = tl_rotate.fit_transform(X_train, y_enc)
# Save rotations
if subj_id_target != subj_id_source:
Q_rot = tl_rotate.rotations_[subj_id_source]
os.makedirs(f'data/rotations/{data_dir_identifier}/{subj_id_target}', exist_ok=True)
np.save(f"{rotations_path}", Q_rot,allow_pickle=False)
np.save(f"{test_data_path}", str_test_data,allow_pickle=False)
# print(f'Rotation saved to {rotations_path}')
if w: # there were warnings
for warn in w:
print(f"Target: {subj_id_target} source: {subj_id_source} iter:{i}: Warning: {warn.message}")
with open(convergece_file, 'a') as f:
f.write(f"Target: {subj_id_target} source: {subj_id_source} iter:{i}: Warning: {warn.message}\n")
else:
print(f"Target: {subj_id_target} source: {subj_id_source} iter:{i}: --------------------")
with open(convergece_file, 'a') as f:
f.write(f"Target: {subj_id_target} source: {subj_id_source} iter:{i}: --------------------\n")
clf.fit(X_train[:(-nbr_target_training_data)],y_train[:(-nbr_target_training_data)])
y_pred = clf.predict(X_test)
rotate_score[i] = accuracy_score(y_test, y_pred)
scores = {
'nothing_score': nothing_score.mean(),
'recenter_score': recenter_score.mean(),
'scale_score': scale_score.mean(),
'rotate_score': rotate_score.mean(),
}
np.save(full_path_accuracies, scores,allow_pickle=True)
return scores
def find_transfer_learning_accuracies(data,iteration_idx,run_logdir='.',data_dir_identifier='.'):
X_covs_target = data[iteration_idx]['X_cov']
y_target = data[iteration_idx]['y']
y_enc_target = data[iteration_idx]['y_enc']
subj_id_target = data[iteration_idx]['subj_id']
subject_nbr_target = data[iteration_idx]['subject_nbr']
data_dir = f'data/transfer_learning/{data_dir_identifier}/{subj_id_target}'
print(f'Creating data folder: {data_dir}')
os.makedirs(f'{data_dir}', exist_ok=True)
full_path_accuracies = f"{data_dir}/all_subj_accuracies.npy"
if os.path.exists(full_path_accuracies):
# If it exists do nothing.
print(f'Accuracies for {subj_id_target} for all subj already exist, skipping...')
return
score_all_subj={}
for source_idx in range(len(data)):
subj_id_source = data[source_idx]['subj_id']
scores = transfer_learning_one_target_vs_one_source(data=data, target_idx=iteration_idx, source_idx=source_idx,run_logdir=run_logdir,data_dir=data_dir,data_dir_identifier=data_dir_identifier)
score_all_subj[subj_id_source] = scores
np.save(full_path_accuracies,score_all_subj,allow_pickle=True)
return
# ---- distance data
def find_distance_between_classes_to_pole(data,iteration_idx,run_logdir='.',data_dir_identifier='.'):
X_covs = data[iteration_idx]['X_cov']
y = data[iteration_idx]['y']
y_enc = data[iteration_idx]['y_enc']
subj_id = data[iteration_idx]['subj_id']
subject_nbr = data[iteration_idx]['subject_nbr']
data_dir = f'data/distance_pole/{data_dir_identifier}/{subj_id}'
print(f'Creating data folder: {data_dir}')
os.makedirs(f'{data_dir}', exist_ok=True)
size_matrix = X_covs[0].shape[0]
classes = np.unique(y)
full_path_distances = f"{data_dir}/distances.npy"
if os.path.exists(full_path_distances):
# If it exists do nothing.
print(f'Disntances for {subj_id} already exist, skipping...')
return
tl_recenter = TLCenter(target_domain=subj_id,metric='riemann')
tl_scale = TLScale(target_domain=subj_id, final_dispersion=1.0, centered_data=True, metric='riemann')
X = tl_recenter.fit_transform(X_covs, y_enc)
X = tl_scale.fit_transform(X, y_enc)
ref = mean_covariance(X)
X_0 = X[y==classes[0]]
X_1 = X[y==classes[1]]
X_mean_0 = mean_covariance(X_0)
X_mean_1 = mean_covariance(X_1)
dist_0 = distance(ref,X_mean_0, metric='riemann', squared=True)
dist_1 = distance(ref,X_mean_1, metric='riemann', squared=True)
pole_north_0 = np.eye(size_matrix)* np.exp(+np.sqrt(dist_0/size_matrix))
pole_south_0 = np.eye(size_matrix)* np.exp(-np.sqrt(dist_0/size_matrix))
pole_north_1 = np.eye(size_matrix)* np.exp(+np.sqrt(dist_1/size_matrix))
pole_south_1 = np.eye(size_matrix)* np.exp(-np.sqrt(dist_1/size_matrix))
# Find intra subject accuracy
clf = MDM(metric=dict(mean='riemann', distance='riemann'))
splitter = StratifiedShuffleSplit(n_splits=4, train_size=0.75, random_state=42) # Shuffels data. Test size = training size in our case.
all_score = np.zeros(splitter.get_n_splits())
for i, (train_index,test_index) in enumerate(splitter.split(X_covs, y)):
X_train = X_covs[train_index]
y_train = y[train_index]
X_test = X_covs[test_index]
y_test = y[test_index]
clf.fit(X_train,y_train)
y_pred = clf.predict(X_test)
all_score[i] = accuracy_score(y_test, y_pred)
distances = {
'subj_id': subj_id,
'Intra-subject accuracy':all_score.mean(),
'ref-mean:0': distance(ref,X_mean_0, metric='riemann', squared=False), # Same distance as for poles
'mean-southpole:0':distance(X_mean_0,pole_south_0, metric='riemann', squared=False),
'mean-northpole:0':distance(X_mean_0,pole_north_0, metric='riemann', squared=False),
'ref-mean:1': distance(ref,X_mean_1, metric='riemann', squared=False), # Same distance as for poles
'mean-southpole:1':distance(X_mean_1,pole_south_1, metric='riemann', squared=False),
'mean-northpole:1':distance(X_mean_1,pole_north_1, metric='riemann', squared=False),
}
np.save(full_path_distances,distances,allow_pickle=True)
return
# -----------------------------------------------
# Multiprocessing helper function:
# -----------------------------------------------
def worker_function(data,iteration_idx,run_logdir='.',data_dir_identifier='.'):
# This function will be passed to the Pool.
print(f"### worker_function ({iteration_idx+1}/{len(data)}) for itr_{iteration_idx}")
# return run_one_iteration(n_matrices, mean,sigma, nbr_outliers,iteration_nbr,nbr_folds=nbr_folds,run_logdir=run_logdir,plotting=plotting)
return run_one_subject(data,iteration_idx,run_logdir=run_logdir,data_dir_identifier=data_dir_identifier)
# -----------------------------------------------
# Run code single or multithreaded
# -----------------------------------------------
def run_code(data,run_logdir='.',identifier='',data_dir_identifier='.'):
if NBR_THREADS == 1:
print('\n-------------------------------')
print( '--> Single thread started < ---')
print( '-------------------------------\n')
for iteration_idx in range(len(data)):
run_one_subject(data,iteration_idx,run_logdir=run_logdir,data_dir_identifier=data_dir_identifier)
else:
# Multithreaded version:
print('\n---------------------------------')
print( '--> Multiprocessing started < ---')
print( '---------------------------------\n')
pool = Pool(processes=NBR_THREADS)
results = []
for iteration_idx in range(len(data)):
res = pool.apply_async(worker_function, args=(data,iteration_idx,run_logdir,data_dir_identifier))
results.append(res)
print('---> Closing pools\n')
pool.close()
pool.join()
print('---> Pools closed \n')
all_res_dict = {}
dist_dict_res = analyse_distance_data(data,data_dir_identifier,run_logdir)
all_res_dict.update(dist_dict_res)
acc_dict_res = analyse_accuracies_for_RPA_steps(data,data_dir_identifier,run_logdir)
all_res_dict.update(acc_dict_res)
return all_res_dict
# ------- Accuracies
def analyse_accuracies_for_RPA_steps(data, data_dir_identifier='.',run_logdir='.'):
print('Data analysis for accuracies')
accuracies_all_subj = []
distances_all_subj = []
list_of_subj = [item['subj_id'] for item in data]
data_dir_acc = f'data/transfer_learning/{data_dir_identifier}'
data_dir_dist = f'data/distance_pole/{data_dir_identifier}'
print('Loading data from iterations...')
for subj in list_of_subj:
subdir = subj
print(f'loading data from {subdir}')
# == Accuracy data
subdir_path = os.path.join(data_dir_acc, subdir)
# Check if the path is indeed a directory
if not os.path.isdir(subdir_path):
print(f'Not a folder: {subdir}')
continue # There might be a .DS_store.
file_path = os.path.join(subdir_path, 'all_subj_accuracies.npy')
# Check if the file exists
if not os.path.exists(file_path):
import pdb; pdb.set_trace()
bug # There is something wrong! The file is missing.
data_here = np.load(file_path,allow_pickle=True).item()
accuracies_all_subj.append(data_here)
# == Dist data
subdir_path = os.path.join(data_dir_dist, subdir)
# Check if the path is indeed a directory
if not os.path.isdir(subdir_path):
print(f'Not a folder: {subdir}')
continue # There might be a .DS_store.
file_path = os.path.join(subdir_path, 'distances.npy')
# Check if the file exists
if not os.path.exists(file_path):
import pdb; pdb.set_trace()
bug # There is something wrong! The file is missing.
data_here = np.load(file_path,allow_pickle=True).item()
distances_all_subj.append(data_here)
print('...done')
print('Plotting...')
all_res_dict = {}
plot_accuracies_data(accuracies_all_subj,list_of_subj,run_logdir=run_logdir,data_dir_identifier=data_dir_identifier)
p_val_dict = plot_accuracy_vs_radius_and_latitude(accuracies_all_subj,distances_all_subj,list_of_subj,run_logdir=run_logdir,data_dir_identifier=data_dir_identifier)
all_res_dict['p-val'] = p_val_dict
print('... done')
return all_res_dict
def plot_accuracy_vs_radius_and_latitude(data_acc,data_dist,list_of_subj,run_logdir='.',data_dir_identifier=''):
# --- accuracy matrix ---
nbr_subjects = len(data_acc)
# Placeholders
nothing_matrix = np.zeros((nbr_subjects,nbr_subjects))
recenter_matrix = np.zeros((nbr_subjects,nbr_subjects))
scale_matrix = np.zeros((nbr_subjects,nbr_subjects))
rotate_matrix = np.zeros((nbr_subjects,nbr_subjects))
# Create matrices
for row, row_dict in enumerate(data_acc):
for col, source_id in enumerate(list_of_subj):
nothing_matrix[row,col] = row_dict[source_id]['nothing_score']
recenter_matrix[row,col] = row_dict[source_id]['recenter_score']
scale_matrix[row,col] = row_dict[source_id]['scale_score']
rotate_matrix[row,col] = row_dict[source_id]['rotate_score']
if data_dir_identifier == 'BNCI2014-001_feet_tongue':
title = 'BCI-IV 2a (feet/tongue)'
extra_str_save = 'BCI_IV_feet_tongue'
elif data_dir_identifier == 'BNCI2014-001_left_hand_right_hand':
title = 'BCI-IV 2a (left/right)'
extra_str_save = 'BCI_IV_left_right'
elif data_dir_identifier == 'Cho2017_left_hand_right_hand':
title = 'Cho (left/right)'
extra_str_save = 'Cho'
elif data_dir_identifier == 'Dreyer2023_left_hand_right_hand':
title = 'Dreyer (left/right)'
extra_str_save = 'Dreyer'
elif data_dir_identifier == 'Lee2019-MI_left_hand_right_hand':
title = 'Lee (left/right)'
extra_str_save = 'Lee'
elif data_dir_identifier == 'PhysionetMotorImagery_feet_hands':
title = 'Physionet (feet/hands)'
extra_str_save = 'physionet_feet_hands'
elif data_dir_identifier == 'PhysionetMotorImagery_left_hand_right_hand':
title = 'Physionet (left/right)'
extra_str_save = 'physionet_left_right'
else:
title = ''
extra_str_save = ''
p_val_dict = plot_pole_ratio_to_acc_correlation_all_subj_for_paper(rotate_matrix,data_dist,list_of_subj,run_logdir=run_logdir,figname=f'lat_to_accuracy_correlation_{extra_str_save}',title=title)
return p_val_dict
def plot_pole_ratio_to_acc_correlation_all_subj_for_paper(accuracy_matrix,data_dist,list_of_subj,run_logdir='.',figname='fig',title=''):
# ----- All subject in one plot -----
nbr_subjects = len(list_of_subj)
distance_vector = np.zeros(nbr_subjects)
for i, subj_id in enumerate(list_of_subj):
distance_vector[i] = data_dist[i]['mean-southpole:0'] / (data_dist[i]['mean-northpole:0'] + data_dist[i]['mean-southpole:0'] )
distance_matrix = np.tile(distance_vector, (nbr_subjects, 1))
diff_to_target_matrix = distance_matrix - distance_vector[:,None] # Source- target
accuracy_matrix_adjusted = accuracy_matrix - np.diag(accuracy_matrix)[:,None]
mask = (np.diag(accuracy_matrix) > 0.55) & (np.diag(accuracy_matrix) < 0.95)
masked_accuracy_matrix_adjusted = accuracy_matrix_adjusted[mask][:,mask]
masked_diff_to_target_matrix = diff_to_target_matrix[mask][:,mask]
fig, axs = plt.subplots(nrows=1,ncols=1,figsize=(5,4),layout="constrained")
# Highlight areas for statistic significance.
axs.axvspan(
xmin=-0.2,
xmax=0.2,
color='lightgray',
alpha=0.5,
edgecolor=None,
)
axs.axvspan(
xmin=-0.1,
xmax=0.1,
color='darkgray',
alpha=0.5,
edgecolor=None,
)
axs.axvspan(
xmin=-0.05,
xmax=0.05,
color='dimgrey',
alpha=0.5,
edgecolor=None,
)
# Plot data
x_to_plot = masked_diff_to_target_matrix.ravel()
y_to_plot = masked_accuracy_matrix_adjusted.ravel() * 100 # * 100 to get %
noise = (rng.random(y_to_plot.shape) - 0.5)*0.05 * 100 # * 100 to get %
print(f'Included subjects in plot_pole_ratio_to_acc_correlation_one_subj() plot: {np.sum(mask)}')
print(np.arange(1,len(mask)+1)[mask])
plot_pole_ratio_to_acc_correlation_one_subj(x_to_plot, y_to_plot+noise, title='',ax=axs)
# Fit line to data.
def model(x, a,b,c):
return a*x**2 + b*x + c
popt, pcov = curve_fit(
model,
x_to_plot,
y_to_plot,
p0=[-1.0,0.0,0.0],
)
a_hat = popt[0]
b_hat = popt[1]
c_hat = popt[2]
# Plot fitted line
x_smooth = np.linspace(x_to_plot.min(), x_to_plot.max(), 300)
y_smooth = model(x_smooth, *popt)
axs.plot(x_smooth, y_smooth, color='black', label=f'{a_hat:.3f}x**2 + {b_hat:.3f}x + {c_hat:.3f}')
plt.title(title,fontsize=12)
plt.ylabel(r'$\Delta$Acc (Improvement vs Intra-subject accuracy)',fontsize=10)
plt.xlabel(r'$\Delta \rho$ (Difference in pole ratio)',fontsize=10)
save_string = f'{figname}' # Filename
plt.savefig(f'{run_logdir}/png_{save_string}.png',dpi=600) # Save the figure in the specified directory
plt.savefig(f'{run_logdir}/pdf_{save_string}.pdf') # Save the figure in the specified directory
plt.close(fig)
p_test_res = {}
# hypotestestning
mask_1 = (np.abs(x_to_plot) <= 0.2)
mask_2 = (np.abs(x_to_plot) > 0.2)
y_bin_0 = y_to_plot[mask_1]
y_bin_1 = y_to_plot[mask_2]
_, p_ttest = ttest_ind(y_bin_0, y_bin_1, equal_var=False,alternative='greater')
p_test_res['x<=0.2|0.2<x'] = p_ttest
p_test_res['Val-diff: x<=0.2|0.2<x'] = np.mean(y_bin_0) - np.mean(y_bin_1)
mask_1 = (np.abs(x_to_plot) <= 0.1)
mask_2 = (np.abs(x_to_plot) > 0.1)
y_bin_0 = y_to_plot[mask_1]
y_bin_1 = y_to_plot[mask_2]
_, p_ttest = ttest_ind(y_bin_0, y_bin_1, equal_var=False,alternative='greater')
p_test_res['x<=0.1|0.1<x'] = p_ttest
p_test_res['Val-diff: x<=0.1|0.1<x'] = np.mean(y_bin_0) - np.mean(y_bin_1)
mask_1 = (np.abs(x_to_plot) <= 0.05)
mask_2 = (np.abs(x_to_plot) > 0.05)
y_bin_0 = y_to_plot[mask_1]
y_bin_1 = y_to_plot[mask_2]
_, p_ttest = ttest_ind(y_bin_0, y_bin_1, equal_var=False,alternative='greater')
p_test_res['x<=0.05|0.05<x'] = p_ttest
p_test_res['Val-diff: x<=0.05|0.05<x'] = np.mean(y_bin_0) - np.mean(y_bin_1)
mask_1 = (np.abs(x_to_plot) <= 0.1)
mask_2 = (np.abs(x_to_plot) > 0.1) & (np.abs(x_to_plot) <= 0.2)
y_bin_0 = y_to_plot[mask_1]
y_bin_1 = y_to_plot[mask_2]
_, p_ttest = ttest_ind(y_bin_0, y_bin_1, equal_var=False,alternative='greater')
p_test_res['x<=0.1|0.1<x<=0.2'] = p_ttest
p_test_res['Val-diff: x<=0.1|0.1<x<=0.2'] = np.mean(y_bin_0) - np.mean(y_bin_1)
mask_1 = (np.abs(x_to_plot) <= 0.05)
mask_2 = (np.abs(x_to_plot) > 0.05) & (np.abs(x_to_plot) <= 0.2)
y_bin_0 = y_to_plot[mask_1]
y_bin_1 = y_to_plot[mask_2]
_, p_ttest = ttest_ind(y_bin_0, y_bin_1, equal_var=False,alternative='greater')
p_test_res['x<=0.05|0.05<x<=0.2'] = p_ttest
p_test_res['Val-diff: x<=0.05|0.05<x<=0.2'] = np.mean(y_bin_0) - np.mean(y_bin_1)
mask_1 = (np.abs(x_to_plot) <= 0.05)
mask_2 = (np.abs(x_to_plot) > 0.05) & (np.abs(x_to_plot) <= 0.1)
y_bin_0 = y_to_plot[mask_1]
y_bin_1 = y_to_plot[mask_2]
_, p_ttest = ttest_ind(y_bin_0, y_bin_1, equal_var=False,alternative='greater')
p_test_res['x<=0.05|0.05<x<=0.1'] = p_ttest
p_test_res['Val-diff: x<=0.05|0.05<x<=0.1'] = np.mean(y_bin_0) - np.mean(y_bin_1)
return p_test_res
def plot_pole_ratio_to_acc_correlation_one_subj(diff_to_target,accuracy_vector,title,ax):
print(f'Plotting {title} lattitude to accuracy correlation')
ax.scatter(diff_to_target,accuracy_vector,alpha=0.25,color='tab:red')
ax.set_title(title)
ax.axvline(0,color='black',linestyle=':')
ax.axhline(0,color='black',linestyle=':')
return
def plot_accuracies_data(data,list_of_subj,run_logdir='.',data_dir_identifier=''):
nbr_subjects = len(data)
# Placeholders
nothing_matrix = np.zeros((nbr_subjects,nbr_subjects))
recenter_matrix = np.zeros((nbr_subjects,nbr_subjects))
scale_matrix = np.zeros((nbr_subjects,nbr_subjects))
rotate_matrix = np.zeros((nbr_subjects,nbr_subjects))
# Create matrices
for row, row_dict in enumerate(data):
for col, source_id in enumerate(list_of_subj):
nothing_matrix[row,col] = row_dict[source_id]['nothing_score']
recenter_matrix[row,col] = row_dict[source_id]['recenter_score']
scale_matrix[row,col] = row_dict[source_id]['scale_score']
rotate_matrix[row,col] = row_dict[source_id]['rotate_score']
mask = (np.diag(nothing_matrix) > 0.55) & (np.diag(nothing_matrix) < 0.95)
list_of_matrices = [nothing_matrix[mask][:,mask],recenter_matrix[mask][:,mask],scale_matrix[mask][:,mask],rotate_matrix[mask][:,mask]]
sub_titles = ['No transfer learning','Recentering','Scale','Rotate (full RPA)']
plot_four_matrices_for_paper(list_of_matrices,sub_titles,run_logdir=run_logdir,filename='four_accuracies',vmin=0.5,vmax=1,title_identifier=data_dir_identifier)
return
def plot_four_matrices_for_paper(list_of_matrices,sub_titles,run_logdir='.',filename='fig',vmin=0.5,vmax=1,title_identifier=''):
figsize = (10,3.25)
fig, axs = plt.subplots(figsize=figsize,ncols=4,nrows=1,layout='constrained') #
for i, ax in enumerate(axs.flatten()):
matrix_to_plot = list_of_matrices[i]
row_sums = np.sum(matrix_to_plot, axis=1) # Sum of rows
col_sums = np.sum(matrix_to_plot, axis=0) # Sum of columns
row_sorted_indices = np.argsort(-row_sums) # Indices to sort rows in descending order
col_sorted_indices = np.argsort(-col_sums)
matrix_to_plot = matrix_to_plot[row_sorted_indices, :]
matrix_to_plot = matrix_to_plot[:, col_sorted_indices]
cax_0 = ax.matshow(matrix_to_plot, cmap='binary',vmin=vmin,vmax=vmax)
# ax.set_title(sub_titles[i],fontsize=12)
ax.set_xticks([])
ax.set_yticks([])
if i ==0:
ax.set_ylabel('Target users',fontsize=10)
ax.set_xlabel('Source users',fontsize=10)
ax.grid(False)
ax.set_title(sub_titles[i],fontsize=12)
if title_identifier == 'BNCI2014-001_feet_tongue':
title = 'BCI-IV 2a (feet/tongue)'
extra_str_save = 'BCI_IV_feet_tongue'
elif title_identifier == 'BNCI2014-001_left_hand_right_hand':
title = 'BCI-IV 2a (left/right)'
extra_str_save = 'BCI_IV_left_right'
elif title_identifier == 'Cho2017_left_hand_right_hand':
title = 'Cho (left/right)'
extra_str_save = 'Cho'
elif title_identifier == 'Dreyer2023_left_hand_right_hand':
title = 'Dreyer (left/right)'
extra_str_save = 'Dreyer'
elif title_identifier == 'Lee2019-MI_left_hand_right_hand':
title = 'Lee (left/right)'
extra_str_save = 'Lee'
elif title_identifier == 'PhysionetMotorImagery_feet_hands':
title = 'Physionet (feet/hands)'
extra_str_save = 'physionet_feet_hands'
elif title_identifier == 'PhysionetMotorImagery_left_hand_right_hand':
title = 'Physionet (left/right)'
extra_str_save = 'physionet_left_right'
else:
title = ''
extra_str_save = ''
plt.suptitle(title,fontsize=14)
# Colorbar
norm = mpl.colors.Normalize(vmin=50, vmax=100)
cmap = mpl.colormaps['binary']
sm = mpl.cm.ScalarMappable(norm=norm, cmap=cmap)
sm.set_array([]) # required but the content is irrelevant
# one colorbar for the whole figure
cbar = fig.colorbar(sm, ax=ax, shrink = 0.75)
cbar.set_label("Accuracy [%]",fontsize=10)
# optional: custom ticks
cbar.set_ticks([50, 60, 70, 80, 90, 100])
save_string = f'{filename}_{extra_str_save}' # Filename
plt.savefig(f'{run_logdir}/png_{save_string}.jpg',dpi=600) # Save the figure in the specified directory
plt.close(fig)
return
# ------ Distance data
def analyse_distance_data(data, data_dir_identifier='.',run_logdir='.'):
print('Data analysis for distances')
distances_all_subj =[]
list_of_subj = [item['subj_id'] for item in data]
data_dir = f'data/distance_pole/{data_dir_identifier}'
print('Loading data from iterations...')
# for subdir in os.listdir(data_folder_path):
for subj in list_of_subj:
subdir = subj
print(f'loading data from {subdir}')
subdir_path = os.path.join(data_dir, subdir)
# Check if the path is indeed a directory
if not os.path.isdir(subdir_path):
print(f'Not a folder: {subdir}')
continue # There might be a .DS_store.
# == Distance data
file_path = os.path.join(subdir_path, 'distances.npy')
# Check if the file exists
if not os.path.exists(file_path):
import pdb; pdb.set_trace()
bug # There is something wrong! The file is missing.
data_here = np.load(file_path,allow_pickle=True).item()
distances_all_subj.append(data_here)
print('...done')
print('Plotting...')
dict_data = plot_pole_ratio_for_paper(distances_all_subj,run_logdir=run_logdir,data_dir_identifier=data_dir_identifier)
print('... done')
return dict_data
def plot_pole_ratio_for_paper(distances_all_subj,run_logdir='.',data_dir_identifier=None):
if data_dir_identifier == 'BNCI2014-001_feet_tongue':
idx_to_highlight = [8]
title = 'BCI-IV 2a (feet/tongue)'
extra_str_save = 'BCI_IV_feet_tongue'
elif data_dir_identifier == 'BNCI2014-001_left_hand_right_hand':
idx_to_highlight = [3]
title = 'BCI-IV 2a (left/right)'
extra_str_save = 'BCI_IV_left_right'
elif data_dir_identifier == 'Cho2017_left_hand_right_hand':
idx_to_highlight = [14]
title = 'Cho (left/right)'
extra_str_save = 'Cho'
elif data_dir_identifier == 'Dreyer2023_left_hand_right_hand':
idx_to_highlight = [81]
title = 'Dreyer (left/right)'
extra_str_save = 'Dreyer'
elif data_dir_identifier == 'Lee2019-MI_left_hand_right_hand':
idx_to_highlight = [48]
title = 'Lee (left/right)'
extra_str_save = 'Lee'
elif data_dir_identifier == 'PhysionetMotorImagery_feet_hands':
idx_to_highlight = [3]
title = 'Physionet (feet/hands)'
extra_str_save = 'physionet_feet_hands'
elif data_dir_identifier == 'PhysionetMotorImagery_left_hand_right_hand':
idx_to_highlight = [76]
title = 'Physionet (left/right)'
extra_str_save = 'physionet_left_right'
else:
idx_to_highlight = [0]
title = ''
extra_str_save = ''
fig, ax = plt.subplots(nrows=1,ncols=1,figsize=(5,4),layout="constrained")
all_pole_ratios = np.zeros(len(distances_all_subj))
for i, dist_dict in enumerate(distances_all_subj):
if i in idx_to_highlight:
pole_ratio = plot_pole_ratio_one_subj_for_paper(dist_dict,ax,highlight=True)
else:
pole_ratio = plot_pole_ratio_one_subj_for_paper(dist_dict,ax,highlight=False)
all_pole_ratios[i] = pole_ratio
# Colorbar
norm = mpl.colors.Normalize(vmin=50, vmax=100)
cmap = mpl.colormaps['viridis']
sm = mpl.cm.ScalarMappable(norm=norm, cmap=cmap)
sm.set_array([]) # required but the content is irrelevant
# one colorbar for the whole figure
cbar = fig.colorbar(sm, ax=ax, shrink = 0.75)
cbar.set_label("Intra-subject accuracy [%]",fontsize=10)
# optional: custom ticks
cbar.set_ticks([50, 60, 70, 80, 90, 100])
ax.set_xlabel('Distance',fontsize=10)
ax.set_ylabel('Distance',fontsize=10)
plt.title(title,fontsize=12)
ax.set_aspect("equal")
save_string = f'distances_to_poles_for_paper_{extra_str_save}' # Filename
plt.savefig(f'{run_logdir}/png_{save_string}.png',dpi=600) # Save the figure in the specified directory
plt.savefig(f'{run_logdir}/pdf_{save_string}.pdf') # Save the figure in the specified directory
plt.close(fig)
dict_to_return = {
'mean_pole_ratio':np.mean(all_pole_ratios),
'std_pole_ratio':np.std(all_pole_ratios),
}
return dict_to_return
def plot_pole_ratio_one_subj_for_paper(dist_dict,ax,highlight=False):
print(f'Plotting: {dist_dict['subj_id']} earth/circle representation.')
# Class 0 circle
r_0 = dist_dict['ref-mean:0']
r_1 = dist_dict['ref-mean:1']
r = (r_0 + r_1) / 2
color_0 = 'tab:orange'
color_1 = 'tab:blue'
acc_score = dist_dict['Intra-subject accuracy']
norm = mpl.colors.Normalize(vmin=0.5, vmax=1.0)
cmap = mpl.colormaps['viridis']
circle_color = cmap(norm(acc_score))
# Circle
circle_poly = CirclePolygon(
xy=(0,0),
radius=r,
resolution=40,
edgecolor=circle_color,
facecolor="none",
linewidth=2,
alpha=0.5,
zorder=1,
)
ax.add_patch(circle_poly)
# Class 0 point,
ratio_away_from_northpole = dist_dict['mean-northpole:0'] / (dist_dict['mean-northpole:0'] + dist_dict['mean-southpole:0'])
angle = ratio_away_from_northpole * np.pi + np.pi/2
ax.scatter(
x = np.cos(angle) * r ,
y = np.sin(angle) * r,
color=color_0,
edgecolor='black' if highlight else color_0,
linewidth = 2 if highlight else 1,
marker='o',
s=100,
alpha=0.75 if highlight else 0.5,
zorder=5 if highlight else 4,
)
# Class 1 point
ratio_away_from_northpole = dist_dict['mean-northpole:1'] / (dist_dict['mean-northpole:1'] + dist_dict['mean-southpole:1'])
angle = ratio_away_from_northpole * np.pi - np.pi/2
ax.scatter(
x = np.cos(-angle) * r,
y = np.sin(-angle) * r,
color=color_1,
edgecolor='black' if highlight else color_1,
linewidth = 2 if highlight else 1,
marker='o',
s=100,
alpha=0.75 if highlight else 0.5,
zorder=5 if highlight else 4,
)
# north and south pole
ax.scatter(
y = [r, -r],
x = [0, 0],
color='tab:gray',
marker='x',
alpha=0.7,
zorder=3,
#s=20,
)
return ratio_away_from_northpole
def highlight_cell(row, col, ax=None, fill=False, **kwargs):
"""
Highlights a specific cell in a matrix plot by drawing a rectangle around it.
This function adds a rectangular patch to highlight a specific cell in a plot (typically a heatmap or matrix plot).
The cell is defined by its row and column indices, and the rectangle can optionally be filled with a color.
Parameters
----------
row : int
The row index of the cell to highlight.
col : int
The column index of the cell to highlight.
ax : matplotlib.axes.Axes, optional
The axes object on which to draw the rectangle. If None, the current axes will be used. Default is None.
fill : bool, optional
Whether to fill the rectangle with color. Default is False.
**kwargs : dict, optional
Additional keyword arguments to pass to `matplotlib.patches.Rectangle` (e.g., edgecolor, facecolor, linewidth).
Returns
-------
rect : matplotlib.patches.Rectangle
The rectangle patch added to the plot.
Notes
-----
This function is typically used to highlight specific cells in heatmaps or matrix plots, where each cell represents
a data point. The rectangle can be customized using the `**kwargs` to control the appearance of the highlight.
"""
y=row
x=col
rect = plt.Rectangle((x-.5, y-.5), 1,1, fill=fill, **kwargs)
ax = ax or plt.gca()
ax.add_patch(rect)
return rect
def visualize_synthetic_2x2_SPD_matrices(run_logdir='.'):
# Sample SPD(2) by sampling a,c>0 and |b|<sqrt(ac)
# matrix = [[a,b],[b,c]]
n = 2000
a = np.random.exponential(1, n)
c = np.random.exponential(1, n)
b_max = np.sqrt(a*c)
b = (2*np.random.rand(n) - 1) * b_max * 0.9 # stay inside
fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(projection="3d")
ax.scatter(1, 0, 1, s=10,color='black')
# Build a grid for the boundary surface ac - b^2 = 0 -> c = b^2 / a
a_max = a.max()
b_max = b.max()
A, C = np.meshgrid(
np.linspace(1e-6, a_max, 60),
np.linspace(1e-6, a_max, 60)
)
# Boundary: b = ±sqrt(a c)
B_pos = np.sqrt(A * C)
B_neg = -B_pos
# "Radius" in SPD metric
x = 1 # choose the distance from I
# Parameters: φ for eigenvalues, θ for rotation
n_phi = 80
n_theta = 80
phi = np.linspace(0, 2*np.pi, n_phi)
theta = np.linspace(0, np.pi, n_theta) # π is enough due to symmetry
Phi, Theta = np.meshgrid(phi, theta)
# Eigenvalues on the circle in log-space
lam1 = np.exp(x * np.cos(Phi))
lam2 = np.exp(x * np.sin(Phi))
cosT = np.cos(Theta)
sinT = np.sin(Theta)
# Map (φ, θ) -> (a, b, c)
A11 = lam1 * cosT**2 + lam2 * sinT**2
A12 = (lam1 - lam2) * cosT * sinT
A22 = lam1 * sinT**2 + lam2 * cosT**2
# Plot the "sphere" around I as a surface in (a,b,c)-space
ax.plot_wireframe(
A11, A12, A22,
rcount=15, ccount=15,
linewidth=0.4,
alpha=0.1,
edgecolor="black",
)
A_flat = np.stack(
[
A11.ravel(),
A12.ravel(),
A12.ravel(),
A22.ravel(),
],
axis=-1, # shape (N, 4)
).reshape(-1, 2, 2)
Q = random_orthogonal(2)
for i in range(10):
idx = rng.choice(len(A_flat))
A = A_flat[idx]
A_rotated = Q @ A @ Q.T
color = f"C{i % 10}"
ax.scatter([A[0,0]], [A[0,1]], [A[1,1]], s=60,marker='*',color=color)
ax.scatter([A_rotated[0,0]], [A_rotated[0,1]], [A_rotated[1,1]], s=40,color=color)
ax.plot([A[0,0], A_rotated[0,0]], [A[0,1], A_rotated[0,1]], [A[1,1], A_rotated[1,1]],color=color,linewidth=1)
lim = 1.5
ax.set_xlim([0,lim*2])
ax.set_ylim([-lim,lim])
ax.set_zlim([0,lim*2])
ax.axis('equal')
ax.set_xlabel(r"$c_{11}$",fontsize=10)
ax.set_ylabel(r"$c_{12}$",fontsize=10)
ax.set_zlabel(r"$c_{22}$",fontsize=10)
elevations = 30 - np.abs(np.linspace(-15,15,11))
aziums = np.linspace(-90,90,11)
for i,(elev,azim) in enumerate(zip(elevations,aziums)):
ax.view_init(elev=elev, azim=azim) # elevation and azimuth in degrees
save_string = f'SPD_synthetic_{i}' # Filename
plt.tight_layout()
plt.savefig(f'{run_logdir}/png_{save_string}.png',dpi=600) # Save the figure in the specified directory
plt.savefig(f'{run_logdir}/pdf_{save_string}.pdf') # Save the figure in the specified directory
plt.close(fig)
return
def random_orthogonal(n):
A = rng.normal(size=(n, n))
Q, R = np.linalg.qr(A)
# Fix sign so det(Q)=±1 stays consistent
d = np.diag(R)
Q *= np.sign(d)
return Q
def combine_images_to_figure(
paths,
out_pdf="combined.pdf",
out_png="combined.png",
ncols=None,
figsize=None,
dpi_png=600,
facecolor="white",
):
"""
Combine raster images into a single Matplotlib figure and save as PDF + PNG.
Parameters
----------
paths : list[str|Path]
Paths to input images.
out_pdf : str|Path
Output PDF path (set to None to skip).
out_png : str|Path
Output PNG path (set to None to skip).
ncols : int|None
Number of columns in the grid. Default: all images in one row.
figsize : tuple[float,float]|None
Figure size in inches. If None, choose a reasonable default.
dpi_png : int
DPI used for PNG output (controls pixel size of the saved PNG).
facecolor : str
Background color for the figure.
Returns
-------
(Path|None, Path|None)
The saved PDF and PNG paths (or None if skipped).
"""
paths = [Path(p) for p in paths]
if not paths:
raise ValueError("paths is empty")
# Load images
imgs = [Image.open(p) for p in paths]
# Layout
n = len(imgs)
if ncols is None:
ncols = n
if ncols <= 0:
raise ValueError("ncols must be >= 1")
nrows = (n + ncols - 1) // ncols
if figsize is None:
# Simple heuristic: ~4 inches per column, ~3 inches per row
figsize = (4 * ncols, 3 * nrows)
fig, axes = plt.subplots(nrows, ncols, figsize=figsize, facecolor=facecolor)
# Make axes always iterable
if nrows == 1 and ncols == 1:
axes = [axes]
else:
axes = axes.ravel()
# Plot images
for ax, im in zip(axes, imgs):
ax.imshow(im)
ax.axis("off")
# Hide unused axes
for ax in axes[len(imgs):]:
ax.axis("off")
plt.tight_layout(pad=0.5, w_pad=0.5, h_pad=0.5)
saved_pdf = None
saved_png = None
if out_pdf is not None:
out_pdf = Path(out_pdf)
fig.savefig(out_pdf, bbox_inches="tight")
saved_pdf = out_pdf
if out_png is not None:
out_png = Path(out_png)
fig.savefig(out_png, dpi=dpi_png, bbox_inches="tight")
saved_png = out_png
plt.close(fig)
return saved_pdf, saved_png
# -----------------------------------------------
# Main
# -----------------------------------------------
# execute only if run as a script
if __name__ == "__main__":
begin_time = datetime.datetime.now()
log_folder = "my_logs"
if not os.path.exists(log_folder):
os.makedirs(log_folder)
# Create run directory to store models and logs in:
root_logdir = os.path.join(os.curdir, log_folder)
run_id = "%s_on_%s" % (time.strftime("run_%Y_%m_%d-%H_%M_%S"), os.uname()[1])
run_logdir = os.path.join(root_logdir, run_id)
output_folder = run_logdir
os.mkdir(run_logdir)
sys.stderr = Logger("{}/_console_err.txt".format(run_logdir))
sys.stdout = Logger("{}/_console_log.txt".format(run_logdir))
print(
"\n### main.py was started on %s in %s with logs in %s"
% (os.uname()[1], os.getcwd(), run_logdir)
)
print(sys.argv)
print_git_version()
save_source_code(run_logdir)
now = begin_time
now_string = str(now).replace(" ", "_").replace(":","_")
write_description()
print('\nRUNNING MAIN CODE\n')
run_logdir_here = f'{run_logdir}/synthetic_data'
os.mkdir(run_logdir_here)
visualize_synthetic_2x2_SPD_matrices(run_logdir=run_logdir_here)
paths = [
f'./{run_logdir}/synthetic_data/png_SPD_synthetic_0.png',
f'./{run_logdir}/synthetic_data/png_SPD_synthetic_2.png',
f'./{run_logdir}/synthetic_data/png_SPD_synthetic_5.png'
]
combine_images_to_figure(
paths,
out_pdf=None,
out_png=f"{run_logdir}/png_synthetic_combined.jpg",
ncols=None,
figsize=(7,5),
dpi_png=600,
facecolor="white",
)
all_res = []
# -----------------------------------------------
# Run the code.
# -----------------------------------------------
for (dataset,events,identifier) in DATASETS:
print()
print()
print()
print('=======================')
print('=======================')
print('=======================')
print(f'Running {identifier}')
print('=======================')
print('=======================')
print('=======================')
print()
events.sort()
run_logdir_dataset = f'{run_logdir}/{identifier}'
data_dir_identifier = f'{dataset.code}_{"_".join(events)}'
os.mkdir(run_logdir_dataset)
all_subject_data = load_data_MOABB(subjects=dataset.subject_list,dataset=dataset,events=events,processing_dict=PROCESSING_DICT,run_logdir=run_logdir_dataset,data_dir_identifier=data_dir_identifier)
res = run_code(data=all_subject_data,run_logdir=run_logdir_dataset,identifier=identifier,data_dir_identifier=data_dir_identifier)
all_res.append(res)
print()
print()
print()
print()
print('RESULTS: ')
for i, (dataset,events,identifier) in enumerate(DATASETS):
print()
print(identifier)
res = all_res[i]
for key in res.keys():
if key == 'p-val':
print(f'p-val: z | q')
for key_2 in res[key]:
print(f'{key_2:20s}: p<0.05: {str(res[key][key_2]<0.05):6s}, p<0.01: {str(res[key][key_2]<0.01):6s}, pval={res[key][key_2]}')
else:
print(f'{key}: {res[key]:.3f} ({res[key]})')
paths = [
f'{run_logdir}/Physio_feet_hands/png_distances_to_poles_for_paper_physionet_feet_hands.png',
f'{run_logdir}/Dreyer_left_right/png_distances_to_poles_for_paper_Dreyer.png',
f'{run_logdir}/BCI-IV_feet_tongue/png_distances_to_poles_for_paper_BCI_IV_feet_tongue.png'
]
combine_images_to_figure(
paths,
out_pdf=None,
out_png=f"{run_logdir}/png_pole_ratio_combined.jpg",
ncols=None,
figsize=(7,5),
dpi_png=600,
facecolor="white",
)
paths = [
f'{run_logdir}/Physio_feet_hands/png_lat_to_accuracy_correlation_physionet_feet_hands.png',
f'{run_logdir}/Dreyer_left_right/png_lat_to_accuracy_correlation_Dreyer.png',
f'{run_logdir}/BCI-IV_feet_tongue/png_lat_to_accuracy_correlation_BCI_IV_feet_tongue.png'
]
combine_images_to_figure(
paths,
out_pdf=None,
out_png=f"{run_logdir}/png_lat_to_accuracy_correlation_combined.jpg",
ncols=None,
figsize=(7,5),
dpi_png=600,
facecolor="white",
)
plt.close('all')
print("\n...work complete!")
print("### Total run time: %s" % (datetime.datetime.now() - begin_time))

rotation_analysis.py at commit f8522b3, under MIT · at the source

Overview

Authors: Frida Heskebeck¹, Bo Bernhardsson¹, Carolina Bergeling²

Institutions: Lund University (Sweden); Blekinge Institute of Technology (Sweden)

Journal: Frontiers in human neuroscience, volume 20, article 1824613

Dates: received 6 March 2026; accepted 30 April 2026; published online 21 May 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.3389/fnhum.2026.1824613 · PMID 42253789 · PMCID PMC13233444 · OpenAlex W7161984432

Open access: gold, a free copy (OpenAlex)

Status: code verified

Methods: Spectral & time-frequency, Statistics, Machine learning

Keywords: BCI, pole ratio, Riemann geometry, source data selection, transfer learning

Topic: EEG and Brain-Computer Interfaces (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: Knut och Alice Wallenbergs Stiftelse

Citations: not cited yet (Europe PMC); 27 references in the paper

Abstract

This paper introduces the pole ratio metric and presents a sphere-based view of symmetric positive-definite matrix rotations on the Riemannian manifold of symmetric positive-definite matrices equipped with the affine-invariant Riemannian metric. The pole ratio quantifies whether data from different users lie on this Riemannian manifold in a way that enables effective transfer learning. The sphere-based view provides insight into the rotational step of transfer learning using the Riemannian Procrustes analysis method and highlights the limitations of rotation. For effective transfer learning, selecting appropriate source data is essential for good performance. The pole ratio is shown to be an effective metric for selecting source data. The main contribution of the paper is the insight into the limitations of rotations on a Riemannian manifold; the usefulness of the pole ratio as a source selection metric is a natural extension of this insight. This paper focuses on Brain-Computer Interfaces (BCIs), but the sphere-based view of rotations of symmetric positive-definite matrix data and the pole ratio are applicable to any field that models two-class data using symmetric positive-definite matrices.

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

Repository

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

gitlab.control.lth.se/fridah/pole_ratio_public

License: MIT

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: f8522b32f7375ffefc79c31ed2e02a2e69e6b2d1, 10 March 2026

Languages: Python (1)

Size: 9 files, 1 script

Software Heritage: not archived

Found in: the text, “Materials and methods”

Holds: README, license file, environment (requirements.txt)

Not found: CITATION.cff, tests, continuous integration, documentation

Tools: Matplotlib (1 file), MOABB (1 file), NumPy (1 file), Pillow (1 file), pyRiemann (1 file), scikit-learn (1 file), SciPy (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

3 files

rotation_analysis.py, Python, 1,746 lines, 3 matches
LICENSE, License, 21 lines
README.md, Text, 42 lines

At the source: gitlab.control.lth.se/FridaH/pole_ratio_public

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:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
1 script, each with its path and the digest of its content;
3 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

Publicly available datasets were analyzed in this study. This data can be found here: https://moabb.neurotechx.com/docs/index.html.

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 2, 28 September 2026

Funding: added Knut och Alice Wallenbergs Stiftelse

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 3 authors, 5 keywords, 25 references.

Cite

This paper

Heskebeck, F., Bernhardsson, B., & Bergeling, C. (2026). Rotation-based metric on the Riemannian manifold of SPD matrices with applications to source data selection for brain-computer interface transfer learning. Frontiers in human neuroscience, 20, 1824613. https://doi.org/10.3389/fnhum.2026.1824613

BibTeX

@article{heskebeck2026rotation, author = {Heskebeck, Frida and Bernhardsson, Bo and Bergeling, Carolina}, title = {{Rotation-based metric on the Riemannian manifold of SPD matrices with applications to source data selection for brain-computer interface transfer learning}}, journal = {Frontiers in human neuroscience}, year = {2026}, month = may, volume = {20}, pages = {1824613}, publisher = {Frontiers Media SA}, issn = {1662-5161}, doi = {10.3389/fnhum.2026.1824613}, url = {https://doi.org/10.3389/fnhum.2026.1824613}, pmid = {42253789}, pmcid = {PMC13233444} }

RIS

TY - JOUR AU - Heskebeck, Frida AU - Bernhardsson, Bo AU - Bergeling, Carolina TI - Rotation-based metric on the Riemannian manifold of SPD matrices with applications to source data selection for brain-computer interface transfer learning T2 - Frontiers in human neuroscience J2 - Front Hum Neurosci PY - 2026 DA - 2026/05/21 VL - 20 SP - 1824613 SN - 1662-5161 PB - Frontiers Media SA DO - 10.3389/fnhum.2026.1824613 UR - https://doi.org/10.3389/fnhum.2026.1824613 LA - en ER -

CSL-JSON

{ "id": "10.3389/fnhum.2026.1824613", "type": "article-journal", "title": "Rotation-based metric on the Riemannian manifold of SPD matrices with applications to source data selection for brain-computer interface transfer learning", "container-title": "Frontiers in human neuroscience", "author": [ { "family": "Heskebeck", "given": "Frida" }, { "family": "Bernhardsson", "given": "Bo" }, { "family": "Bergeling", "given": "Carolina" } ], "container-title-short": "Front Hum Neurosci", "volume": "20", "page": "1824613", "DOI": "10.3389/fnhum.2026.1824613", "PMID": "42253789", "PMCID": "PMC13233444", "ISSN": "1662-5161", "publisher": "Frontiers Media SA", "URL": "https://doi.org/10.3389/fnhum.2026.1824613", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 21 ] ] } }

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

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.