diff --git a/invokers/invoker_estimate_noise_cov_matrix.py b/invokers/invoker_estimate_noise_cov_matrix.py new file mode 100644 index 00000000..bc0075d2 --- /dev/null +++ b/invokers/invoker_estimate_noise_cov_matrix.py @@ -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() 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 842f9f41..ce0ac379 100755 --- 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 75654835..cc241b6b 100755 --- a/kymata/preproc/data_cleansing.py +++ b/kymata/preproc/data_cleansing.py @@ -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 @@ -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, @@ -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, @@ -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) @@ -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))) @@ -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') @@ -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() @@ -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' @@ -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 \ No newline at end of file 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(