Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add the function to calculate noise covariance matrix #119

Merged
merged 3 commits into from
Jan 12, 2024
Merged

Add the function to calculate noise covariance matrix #119

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d32c0cd
Add the function to calculate noise covariance matrix
young-x-skyee Jan 12, 2024
7cf2e78
Merge branch 'main' into big-preprocessing-refactor-cov
neukym Jan 12, 2024
4c16649
Update invoker_estimate_noise_cov_matrix.py
neukym Jan 12, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 31 additions & 0 deletions invokers/invoker_estimate_noise_cov_matrix.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
from pathlib import Path
from colorama import Fore

from kymata.io.yaml import load_config
from kymata.io.cli import print_with_color
from kymata.preproc.data_cleansing import estimate_noise_cov


# noinspection DuplicatedCode
def main():
config = load_config(str(Path(Path(__file__).parent.parent, "kymata", "config", "dataset4.yaml")))

if config['data_location'] == "local":
data_root_dir = str(Path(Path(__file__).parent.parent, "kymata-toolbox-data", "emeg_study_data")) + "/"
elif config['data_location'] == "cbu":
data_root_dir = '/imaging/projects/cbu/kymata/data/'
elif config['data_location'] == "cbu-local":
data_root_dir = '//cbsu/data/imaging/projects/cbu/kymata/data/'
else:
raise Exception("The 'data_location' parameter in the config file must be either 'cbu' or 'local' or 'cbu-local'.")

estimate_noise_cov( data_root_dir = data_root_dir,
emeg_machine_used_to_record_data = config['emeg_machine_used_to_record_data'],
list_of_participants = config['list_of_participants'],
dataset_directory_name = config['dataset_directory_name'],
n_runs = config['number_of_runs'],
cov_method = config['cov_method']
)

if __name__ == '__main__':
main()
2 changes: 1 addition & 1 deletion invokers/invoker_run_data_cleansing.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ def main():
elif config['data_location'] == "cbu-local":
data_root_dir = '//cbsu/data/imaging/projects/cbu/kymata/data/'
else:
raise Exception("The 'data_location' parameter in the config file must be either 'cbu' or 'local'.")
raise Exception("The 'data_location' parameter in the config file must be either 'cbu' or 'local' or 'cbu-local'.")

run_first_pass_cleansing_and_maxwell_filtering(
data_root_dir = data_root_dir,
Expand Down
4 changes: 2 additions & 2 deletions invokers/invoker_run_hexel_current_estimation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,9 +17,9 @@ def main():
elif config['data_location'] == "cbu-local":
data_root_dir = '//cbsu/data/imaging/projects/cbu/kymata/data/'
else:
raise Exception("The 'data_location' parameter in the config file must be either 'cbu' or 'local'.")
raise Exception("The 'data_location' parameter in the config file must be either 'cbu' or 'local' or 'cbu-local'.")

create_current_estimation_prerequisites(data_root_dir, config=config)
# create_current_estimation_prerequisites(data_root_dir, config=config)

create_forward_model_and_inverse_solution(data_root_dir, config=config)

Expand Down
3 changes: 3 additions & 0 deletions kymata/config/dataset4.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,9 @@ supress_excessive_plots_and_prompts: True
eeg: True
meg: True

# Method to estimate noise covariance matrix
cov_method: 'grand_ave' # grand_ave | empty_room | run_start

# Creation of evoked trial pipeline
audio_delivery_latency: 26 # in milliseconds
visual_delivery_latency: 17 # in milliseconds
Expand Down
59 changes: 44 additions & 15 deletions kymata/preproc/data_cleansing.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
from numpy import zeros, nanmin
from pandas import DataFrame, Index
from pathlib import Path
from matplotlib import pyplot as plt

from kymata.io.cli import print_with_color, input_with_color
from kymata.io.yaml import load_config, modify_param_config
Expand Down Expand Up @@ -297,6 +298,43 @@ def _remove_ecg_eog(filt_raw, ica):
ica.plot_scores(eog_scores)


def estimate_noise_cov(data_root_dir: str,
emeg_machine_used_to_record_data: str,
list_of_participants: list[str],
dataset_directory_name: str,
n_runs: int,
cov_method: str,
):
for p in list_of_participants:
if cov_method == 'grand_ave':
cleaned_raws = []
for run in range(1, n_runs + 1):
raw_fname = data_root_dir + dataset_directory_name + '/intrim_preprocessing_files/2_cleaned/' + p + '_run' + str(run) + '_cleaned_raw.fif.gz'
raw = mne.io.Raw(raw_fname, preload=True)
cleaned_raws.append(raw)
raw = mne.io.concatenate_raws(raws=cleaned_raws, preload=True)
cov = mne.compute_raw_covariance(raw, tmin=300, tmax=310, return_estimators=True)
mne.write_cov(data_root_dir + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/covariance_grand_average/' + p + '-auto-cov-300.fif', cov)

elif cov_method == 'empty_room':
emptyroom_raw = mne.read_raw()
emptyroom_raw = mne.preprocessing.maxwell_filter_prepare_emptyroom(emptyroom_raw)

fine_cal_file = str(Path(Path(__file__).parent.parent.parent, 'kymata-toolbox-data', 'cbu_specific_files/SSS/sss_cal_' + emeg_machine_used_to_record_data + '.dat'))
crosstalk_file = str(Path(Path(__file__).parent.parent.parent, 'kymata-toolbox-data', 'cbu_specific_files/SSS/ct_sparse_' + emeg_machine_used_to_record_data + '.fif'))


raw_fif_data_sss = mne.preprocessing.maxwell_filter(
emptyroom_raw,
cross_talk=crosstalk_file,
calibration=fine_cal_file,
st_correlation=0.980,
st_duration=10,
verbose=True)

cov = mne.compute_raw_covariance(raw_fif_data_sss, tmin=0, tmax=10, return_estimators=True)
mne.write_cov(data_root_dir + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/covariance_grand_average/' + p + '-auto-cov-emptyroom.fif', cov)

def create_trialwise_data(dataset_directory_name: str,
list_of_participants: list[str],
repetitions_per_runs: int,
Expand Down Expand Up @@ -414,8 +452,8 @@ def create_trialwise_data(dataset_directory_name: str,
evoked = epochs.average()
evoked.save(f'{data_path}/intrim_preprocessing_files/{save_folder}/evoked_data/{p}-ave.fif', overwrite=True)


def create_trials(dataset_directory_name: str,
def create_trials(data_root_dir: str,
dataset_directory_name: str,
list_of_participants: list[str],
repetitions_per_runs: int,
number_of_runs: int,
Expand Down Expand Up @@ -443,7 +481,7 @@ def create_trials(dataset_directory_name: str,
cleaned_raws = []

for run in range(1, number_of_runs + 1):
raw_fname = '/imaging/projects/cbu/kymata/data/' + dataset_directory_name + '/intrim_preprocessing_files/2_cleaned/' + p + '_run' + str(run) + '_cleaned_raw.fif.gz'
raw_fname = data_root_dir + dataset_directory_name + '/intrim_preprocessing_files/2_cleaned/' + p + '_run' + str(run) + '_cleaned_raw.fif.gz'
raw = mne.io.Raw(raw_fname, preload=True)
cleaned_raws.append(raw)

Expand Down Expand Up @@ -510,7 +548,7 @@ def create_trials(dataset_directory_name: str,
# Log which channels are worst
dropfig = epochs.plot_drop_log(subject=p)
dropfig.savefig(
'/imaging/projects/cbu/kymata/data/' + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/logs/' + input_stream + '_drop-log_' + p + '.jpg')
data_root_dir + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/logs/' + input_stream + '_drop-log_' + p + '.jpg')

global_droplog.append('[' + input_stream + ']' + p + ':' + str(epochs.drop_log_stats(epochs.drop_log)))

Expand All @@ -532,13 +570,6 @@ def create_trials(dataset_directory_name: str,
# 'data/' + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/evoked_grand_average/' + p + '-grandave.fif',
# overwrite=True)

# save grand covs
print_with_color(f"... save grand covariance matrix", Fore.GREEN)

cov = mne.compute_raw_covariance(raw, tmin=0, tmax=10, return_estimators=True)
mne.write_cov('/imaging/projects/cbu/kymata/data/' + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/covariance_grand_average/' + p + '-auto-cov.fif', cov)


# Save global droplog
# with open('data/' + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/logs/drop-log.txt', 'a') as file:
# file.write('Average drop rate for each participant\n')
Expand All @@ -563,6 +594,7 @@ def create_trials(dataset_directory_name: str,
# 'data/' + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/evoked_data/' + input_stream + '/item' + str(
# trial) + '-ave.fif', overwrite=True)


def apply_automatic_bad_channel_detection(raw_fif_data: mne.io.Raw, machine_used: str, plot: bool = True):
"""Apply Automatic Bad Channel Detection"""
raw_check = raw_fif_data.copy()
Expand All @@ -588,10 +620,7 @@ def apply_automatic_bad_channel_detection(raw_fif_data: mne.io.Raw, machine_used

return raw_fif_data


def _plot_bad_chans(auto_scores):
import seaborn as sns
from matplotlib import pyplot as plt

# Only select the data for gradiometer channels.
ch_type = 'grad'
Expand Down Expand Up @@ -635,4 +664,4 @@ def _plot_bad_chans(auto_scores):
fig.tight_layout(rect=[0, 0.03, 1, 0.95])

# Replace the word “noisy” with “flat”, and replace
# vmin=nanmin(limits) with vmax=nanmax(limits) to print flat channels
# vmin=nanmin(limits) with vmax=nanmax(limits) to print flat channels
56 changes: 28 additions & 28 deletions kymata/preproc/hexel_current_estimation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -171,33 +171,33 @@ def create_forward_model_and_inverse_solution(data_root_dir, config: dict):
intrim_preprocessing_directory_name = Path(data_root_dir, dataset_directory_name, "intrim_preprocessing_files")
mri_structurals_directory = Path(data_root_dir, dataset_directory_name, config['mri_structurals_directory'])

# Compute forward solution
for participant in list_of_participants:
# # Compute forward solution
# for participant in list_of_participants:

fwd = mne.make_forward_solution(
# Path(Path(path.abspath("")), "data",
Path(data_root_dir,
dataset_directory_name,
'raw_emeg', participant, participant +
'_run1_raw.fif'), # note this file is only used for the sensor positions.
trans=Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","coregistration_files", participant + '-trans.fif'),
src=Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","src_files", participant + '_ico5-src.fif'),
bem=Path(mri_structurals_directory, participant, "bem", participant + '-5120-5120-5120-bem-sol.fif'),
meg=config['meg'],
eeg=config['eeg'],
mindist=5.0,
n_jobs=None,
verbose=True,
)
print(fwd)
if config['meg'] and config['eeg']:
mne.write_forward_solution(Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","forward_sol_files", participant + '-fwd.fif'), fwd)
elif config['meg']:
mne.write_forward_solution(Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","forward_sol_files", participant + '-fwd-megonly.fif'), fwd)
elif config['eeg']:
mne.write_forward_solution(Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","forward_sol_files", participant + '-fwd-eegonly.fif'), fwd)
else:
raise Exception('eeg and meg in the config file cannot be both False')
# fwd = mne.make_forward_solution(
# # Path(Path(path.abspath("")), "data",
# Path(data_root_dir,
# dataset_directory_name,
# 'raw_emeg', participant, participant +
# '_run1_raw.fif'), # note this file is only used for the sensor positions.
# trans=Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","coregistration_files", participant + '-trans.fif'),
# src=Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","src_files", participant + '_ico5-src.fif'),
# bem=Path(mri_structurals_directory, participant, "bem", participant + '-5120-5120-5120-bem-sol.fif'),
# meg=config['meg'],
# eeg=config['eeg'],
# mindist=5.0,
# n_jobs=None,
# verbose=True,
# )
# print(fwd)
# if config['meg'] and config['eeg']:
# mne.write_forward_solution(Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","forward_sol_files", participant + '-fwd.fif'), fwd)
# elif config['meg']:
# mne.write_forward_solution(Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","forward_sol_files", participant + '-fwd-megonly.fif'), fwd)
# elif config['eeg']:
# mne.write_forward_solution(Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","forward_sol_files", participant + '-fwd-eegonly.fif'), fwd)
# else:
# raise Exception('eeg and meg in the config file cannot be both False')

# Compute inverse operator

Expand All @@ -224,7 +224,7 @@ def create_forward_model_and_inverse_solution(data_root_dir, config: dict):
intrim_preprocessing_directory_name,
'3_evoked_sensor_data',
'covariance_grand_average',
participant + '-auto-cov.fif')))
participant + '-auto-cov-300.fif')))
# note this file is only used for the sensor positions.
raw = mne.io.Raw(Path(
Path(path.abspath("")),
Expand All @@ -249,7 +249,7 @@ def create_forward_model_and_inverse_solution(data_root_dir, config: dict):
intrim_preprocessing_directory_name,
'4_hexel_current_reconstruction',
'inverse-operators',
participant + '_ico5-3L-loose02-cps-nodepth.fif')),
participant + '_ico5-3L-loose02-cps-nodepth-300.fif')),
inverse_operator)
elif config['meg']:
mne.minimum_norm.write_inverse_operator(
Expand Down