From d32c0cd2c20e30a3a9956f602fd6da65f9569f3e Mon Sep 17 00:00:00 2001 From: young-x-skyee Date: Fri, 12 Jan 2024 09:43:17 +0000 Subject: [PATCH 1/2] Add the function to calculate noise covariance matrix --- invokers/invoker_estimate_noise_cov_matrix.py | 25 +++++++++ invokers/invoker_run_data_cleansing.py | 2 +- .../invoker_run_hexel_current_estimation.py | 4 +- kymata/config/dataset4.yaml | 3 + kymata/preproc/data_cleansing.py | 43 +++++++++++--- kymata/preproc/hexel_current_estimation.py | 56 +++++++++---------- 6 files changed, 93 insertions(+), 40 deletions(-) create mode 100644 invokers/invoker_estimate_noise_cov_matrix.py diff --git a/invokers/invoker_estimate_noise_cov_matrix.py b/invokers/invoker_estimate_noise_cov_matrix.py new file mode 100644 index 00000000..c1cbe322 --- /dev/null +++ b/invokers/invoker_estimate_noise_cov_matrix.py @@ -0,0 +1,25 @@ +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, config) + +if __name__ == '__main__': + main() diff --git a/invokers/invoker_run_data_cleansing.py b/invokers/invoker_run_data_cleansing.py index 5684d29a..60af6f8f 100644 --- a/invokers/invoker_run_data_cleansing.py +++ b/invokers/invoker_run_data_cleansing.py @@ -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, diff --git a/invokers/invoker_run_hexel_current_estimation.py b/invokers/invoker_run_hexel_current_estimation.py index 1c6a925b..6ecbe885 100644 --- a/invokers/invoker_run_hexel_current_estimation.py +++ b/invokers/invoker_run_hexel_current_estimation.py @@ -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) diff --git a/kymata/config/dataset4.yaml b/kymata/config/dataset4.yaml index 26bf33ae..14a7fdab 100644 --- a/kymata/config/dataset4.yaml +++ b/kymata/config/dataset4.yaml @@ -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 diff --git a/kymata/preproc/data_cleansing.py b/kymata/preproc/data_cleansing.py index ca7c2321..ea0a6a05 100644 --- a/kymata/preproc/data_cleansing.py +++ b/kymata/preproc/data_cleansing.py @@ -298,9 +298,36 @@ def _remove_ecg_eog(filt_raw, ica): ica.plot_scores(eog_scores) - - -def create_trials(dataset_directory_name: str, +def estimate_noise_cov(data_root_dir, config): + for p in config['list_of_participants']: + if config['cov_method'] == 'grand_ave': + cleaned_raws = [] + for run in range(1, config['number_of_runs'] + 1): + raw_fname = data_root_dir + config['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 + config['dataset_directory_name'] + '/intrim_preprocessing_files/3_evoked_sensor_data/covariance_grand_average/' + p + '-auto-cov-300.fif', cov) + + elif config['cov_method'] == 'empty_room': + emptyroom_raw = mne.read_raw() + emptyroom_raw = mne.preprocessing.maxwell_filter_prepare_emptyroom(emptyroom_raw) + 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) + + raw_fif_data_sss.save(saved_maxfiltered_filename, fmt='short') + 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_trials(data_root_dir: str, + dataset_directory_name: str, list_of_participants: list[str], repetitions_per_runs: int, number_of_runs: int, @@ -328,7 +355,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) @@ -395,7 +422,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))) @@ -417,11 +444,9 @@ 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 diff --git a/kymata/preproc/hexel_current_estimation.py b/kymata/preproc/hexel_current_estimation.py index 6527eb23..3d121727 100644 --- a/kymata/preproc/hexel_current_estimation.py +++ b/kymata/preproc/hexel_current_estimation.py @@ -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 @@ -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("")), @@ -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( From 4c16649959be4838a24befd5ba10dc78715e1f26 Mon Sep 17 00:00:00 2001 From: neukym Date: Fri, 12 Jan 2024 22:23:33 +0000 Subject: [PATCH 2/2] Update invoker_estimate_noise_cov_matrix.py --- invokers/invoker_estimate_noise_cov_matrix.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/invokers/invoker_estimate_noise_cov_matrix.py b/invokers/invoker_estimate_noise_cov_matrix.py index c1cbe322..bc0075d2 100644 --- a/invokers/invoker_estimate_noise_cov_matrix.py +++ b/invokers/invoker_estimate_noise_cov_matrix.py @@ -19,7 +19,13 @@ def main(): 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, config) + 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()