diff --git a/.gitignore b/.gitignore index 0d85b61ef..e98bb5dce 100644 --- a/.gitignore +++ b/.gitignore @@ -133,3 +133,10 @@ docs/generated/ .vscode/ phys2bids/tests/data/* /phys2bids/tests/data/ +.idea/workspace.xml +.idea/vcs.xml +.idea/phys2bids.iml +.idea/modules.xml +.idea/misc.xml +.idea/inspectionProfiles/Project_Default.xml +.idea/inspectionProfiles/profiles_settings.xml diff --git a/.idea/.gitignore b/.idea/.gitignore new file mode 100644 index 000000000..26d33521a --- /dev/null +++ b/.idea/.gitignore @@ -0,0 +1,3 @@ +# Default ignored files +/shelf/ +/workspace.xml diff --git a/phys2bids/slice4phys.py b/phys2bids/slice4phys.py index e4e491a57..dc505d268 100644 --- a/phys2bids/slice4phys.py +++ b/phys2bids/slice4phys.py @@ -5,6 +5,8 @@ import numpy as np +from .physio_obj import BlueprintInput + LGR = logging.getLogger(__name__) @@ -121,6 +123,118 @@ def find_takes(phys_in, ntp_list, tr_list, thr=None, padding=9): return take_timestamps +def retrieve_triggers(trig_MRI): + """ + Convert trigger onsets to a list of indices & retrieve length of data + + Parameters + --------- + trig_MRI + + Returns + -------- + onsets_MRI + len_trig + """ + len_trig = len(trig_MRI) + onsets_MRI = [ + int(x) for x in np.where(trig_MRI > np.mean(trig_MRI))[0] + ] # np.where(trig_8 != 0)[0] + + return onsets_MRI, len_trig + + +def get_onsets_run(onsets_MRI, sensitivity): + """ + onsets_MRI: list + list of timepoints when MRI trigger is received + sensitivity: int + percentile of values to be considered as a 'regular' interval of two consecutive MRI triggers + """ + # get trial + values_offsets = np.array(onsets_MRI[1:]) - np.array(onsets_MRI[0:-1]) + values_offsets = np.unique(values_offsets)[1:] + check_val = np.percentile(values_offsets, sensitivity) + + idx_offset = np.where([x > check_val for x in values_offsets])[0] + idx_onset = idx_offset + 1 + + idx_offset = np.append(idx_offset, len(onsets_MRI) - 1) + idx_onset = np.insert(idx_onset, 0, 0) # get all on/ offsets + + run_onset_raw = np.array([onsets_MRI[x] for x in idx_onset]) + run_offset_raw = np.array([onsets_MRI[x] for x in idx_offset]) + + return run_onset_raw, run_offset_raw + + +def clean_onsets_run_GUI(onsets_MRI, run_onset_raw, run_offset_raw, sampling_rate, padding): + """ + Visual check of run onsets/ offsets + :param onsets_MRI: list + list of indices where the MRI triggers are received + :param run_onset_raw: numpy array + array of run onsets + :param run_offset_raw: numpy array + array of the run offsets + :param sampling_rate: int + sampling frequency of MRI trigger + :return: onsets_start: numpy array + array of run onsets + :return: onsets_end: numpy array + array of run offsets + """ + try: + import peakdet + except ImportError: + raise ImportError( + "peakdet is required for creating the data per run." "Please see install instructions." + ) + + conditions = np.zeros(max(onsets_MRI) + round(padding * sampling_rate)) + conditions[onsets_MRI] = 1 # times when MRI triggers are received + + conditions[run_onset_raw] = 1.5 + conditions[run_offset_raw] = -0.5 + + data = peakdet.Physio(conditions, fs=sampling_rate) + + data._metadata["peaks"] = run_onset_raw + data._metadata["troughs"] = run_offset_raw + data = peakdet.operations.edit_run(data) + onsets_start = data.peaks + onsets_end = data.troughs + + # data._metadata['peaks'] = np.sort(np.concatenate((run_onset_raw, run_offset_raw))) + # data._metadata['troughs'] = np.array([1]) + # onsets_start = data.peaks[::2] # odd numbered elements + # onsets_end = data.peaks[1::2] # even numbered elements + + return onsets_start, onsets_end + + +def split_signal_to_runs(phys_in, sensitivity, manual_check, padding=9): + trig_MRI = phys_in.timeseries[phys_in.trigger_idx] + sampling_rate = phys_in.freq[phys_in.trigger_idx] + + # get index of MRI triggers and length of data + onsets_MRI, len_trig = retrieve_triggers(trig_MRI) + + # calculate run onsets & offsets + run_onset_raw, run_offset_raw = get_onsets_run(onsets_MRI, sensitivity) + + # manual check of run onsets & offsets + if manual_check: + onsets_start, onsets_end = clean_onsets_run_GUI( + onsets_MRI, run_onset_raw, run_offset_raw, sampling_rate, padding + ) + else: + onsets_start = run_onset_raw + onsets_end = run_offset_raw + + return onsets_start, onsets_end + + def slice4phys(phys_in, ntp_list, tr_list, thr, padding=9): """ Slice takes for phys2bids. @@ -168,3 +282,147 @@ def slice4phys(phys_in, ntp_list, tr_list, thr, padding=9): ) return phys_in_slices + + +def get_onsets_run(trig_MRI, sensitivity=20): + """ + :param: trig_MRI: numpy array + MRI trigger time course + :param: sensitivity: int + percentile of values to be considered as a 'regular' interval of two consecutive MRI triggers + :return: onsets_MRI: list + time points corresponding to when the MRI trigger was received + :return: run_onset_raw: numpy array + array of run onsets before visual inspection + :return: run_offset_raw: numpy array + array of run offsets before visual inspection + """ + + # convert list to np.array + onsets_MRI = np.array([int(x) for x in np.where(trig_MRI > np.mean(trig_MRI))[0]]) + + # difference of timepoint n and timepoint (n-1) + difference = np.array(onsets_MRI[1:]) - np.array(onsets_MRI[0:-1]) + diff_unique = np.unique(difference)[1:] + + # get threshold of values in which the gap between two time points signify a new run + check_val = np.percentile(diff_unique, sensitivity) + + # list of gap between two time points to be kept to separate runs + list_diff_keep = diff_unique[[x > check_val for x in diff_unique]] + + # find position in onsets_MRI where the difference is in list_diff_keep + keep = [d in list_diff_keep for d in difference] + run_offset_raw = onsets_MRI[:-1] + run_offset_raw = [run_offset_raw[i] for i, flag in enumerate(keep) if flag] + run_onset_raw = onsets_MRI[1:] + run_onset_raw = [run_onset_raw[i] for i, flag in enumerate(keep) if flag] + + run_offset_raw = np.append(run_offset_raw, onsets_MRI[-1]) + run_onset_raw = np.insert(run_onset_raw, 0, onsets_MRI[0]) + + return onsets_MRI, run_onset_raw, run_offset_raw + + +def clean_onsets_run_GUI(onsets_MRI, run_onset_raw, run_offset_raw, sampling_rate, padding): + """ + Visual check of run onsets/ offsets + + :param onsets_MRI: list + list of indices where the MRI triggers are received + :param run_onset_raw: numpy array + array of run onsets + :param run_offset_raw: numpy array + array of the run offsets + :param sampling_rate: int + sampling frequency of MRI trigger + :param padding: int + number of seconds to pad the signal with + :return: onsets_start: numpy array + array of run onsets + :return: onsets_end: numpy array + array of run offsets + """ + try: + import peakdet + except ImportError: + raise ImportError( + "peakdet is required for creating the data per run." "Please see install instructions." + ) + + conditions = np.zeros(max(onsets_MRI) + round(1000 * sampling_rate)) + conditions[onsets_MRI] = 1 # times when MRI triggers are received + + conditions[run_onset_raw] = 1.5 + conditions[run_offset_raw] = -0.5 + + data = peakdet.Physio(conditions, fs=sampling_rate) + + data._metadata["peaks"] = run_onset_raw + data._metadata["troughs"] = run_offset_raw + data = peakdet.operations.edit_physio(data) + + # return start and end times at (padding) seconds before and after the real on and offsets + onsets_start = data.peaks - padding * sampling_rate + onsets_end = data.troughs + padding * sampling_rate + + return onsets_start, onsets_end + + +def slice_runs(phys_in, onsets, offsets): + """ + Slice runs according on list of onsets and offsets + + :param phys_in: BlueprintInput + Object returned by BlueprintInput class + :param onsets: numpy array + array of run onsets + :param offsets: numpy array + array of run offsets + :return run_data: dict + dictionary containing (number of runs) keys and BlueprintInput as values + """ + run_data = {} + + for i in range(len(offsets)): + signals = [phys_signal[onsets[i] : offsets[i]] for phys_signal in phys_in.timeseries] + signals_run = BlueprintInput( + signals, phys_in.freq, phys_in.ch_name, phys_in.units, phys_in.trigger_idx + ) + + run_data[i] = signals_run + + return run_data + + +def split_signal_to_runs(phys_in, manual_check=True, sensitivity=20, padding=9): + """ " + :param phys_in: BlueprintInput + Object returned by BlueprintInput class + :param manual_check: bool + whether there should be a manual check of the automated on/offsets detection + :param sensitivity: int + percentile of values to be considered as a 'regular' interval of two consecutive MRI triggers + :param padding: int + number of seconds to pad before and after each run + :return run_data: dict + dictionary containing (number of runs) keys and BlueprintInput as values + """ + trig_MRI = phys_in.timeseries[phys_in.trigger_idx] + sampling_rate = phys_in.freq[phys_in.trigger_idx] + + # calculate run onsets & offsets + onsets_MRI, run_onset_raw, run_offset_raw = get_onsets_run(trig_MRI, sensitivity) + + # manual check of run onsets & offsets + if manual_check: + onsets_start, onsets_end = clean_onsets_run_GUI( + onsets_MRI, run_onset_raw, run_offset_raw, sampling_rate, padding + ) + else: + onsets_start = run_onset_raw + onsets_end = run_offset_raw + + run_signal = slice_runs(phys_in, onsets_start, onsets_end) + + return run_signal, onsets_start, onsets_end diff --git a/phys2bids/test_slice4phys.py b/phys2bids/test_slice4phys.py new file mode 100644 index 000000000..5ab8ba495 --- /dev/null +++ b/phys2bids/test_slice4phys.py @@ -0,0 +1,40 @@ +import importlib + +import neurokit2 as nk +import numpy as np + +from phys2bids import physio_obj as po +from phys2bids import slice4phys as slice + +importlib.reload(slice) + +sub_idx = 27 +path_raw = "D:\\data\\final_experiment\\physio\\%s\\sub-%02d" % ("raw", sub_idx) +data_exp = nk.read_acqknowledge( + path_raw + "\\sub-%02d.acq" % (sub_idx), + sampling_rate=500, + resample_method="interpolation", + impute_missing=True, +) +data_exp = data_exp[0] + + +def test_split_signal_to_runs(path_data, data_exp): + # test_time = np.array([0, 1, 1, 2, 3, 5, 8, 13]) + # test_trigger = np.array([0, 1, 0, 0, 0, 0, 0, 0]) + # test_chocolate = np.array([1, 0, 0, 1, 0, 0, 1, 0]) + # test_timeseries = [test_time, test_trigger, test_chocolate] + test_timeseries = [data_exp.trig_1.to_numpy(), data_exp.trig_8.to_numpy()] + test_freq = [200, 500] + test_chn_name = ["channel_1", "trigger"] + test_units = ["s", "s"] + test_chtrig = 1 + + blueprint_in = po.BlueprintInput( + test_timeseries, test_freq, test_chn_name, test_units, test_chtrig + ) + + importlib.reload(slice) + run_onset_raw, run_offset_raw = slice.split_signal_to_runs( + blueprint_in, sensitivity=30, manual_check=True, padding=9 + ) diff --git a/phys2bids/tests/test_slice4phys.py b/phys2bids/tests/test_slice4phys.py new file mode 100644 index 000000000..deb5db4a4 --- /dev/null +++ b/phys2bids/tests/test_slice4phys.py @@ -0,0 +1,28 @@ +import numpy as np +import pandas as pd + +import phys2bids.physio_obj as po +import phys2bids.slice4phys as slice + + +def test_slice_signal_to_runs(path_data): + # read in data + data_exp = pd.read_csv(path_data, sep="\t") + + # create blueprintinput + test_chocolate = np.array(data_exp.EDA) + test_trigger = np.array(data_exp.trig_8) + test_timeseries = [test_trigger, test_chocolate] + test_freq = [500, 100] + test_chn_name = ["trigger", "chocolate"] + test_units = ["s", "sweetness"] + test_chtrig = 0 + + phys_in = po.BlueprintInput(test_timeseries, test_freq, test_chn_name, test_units, test_chtrig) + + # slice phys input into runs + run_data, onsets_start, onsets_end = slice.split_signal_to_runs( + phys_in, manual_check=True, sensitivity=20, padding=9 + ) + + assert len(run_data) == len(onsets_start) == len(onsets_end) diff --git a/setup.cfg b/setup.cfg index f9208504e..ed9f7a2a0 100644 --- a/setup.cfg +++ b/setup.cfg @@ -21,7 +21,6 @@ long_description_content_type = text/markdown; charset=UTF-8 platforms = OS Independent provides = phys2bids - [options] python_requires = >=3.6.1 install_requires = @@ -34,16 +33,25 @@ test_suite = pytest zip_safe = False packages = find: include_package_data = True - [options.extras_require] spike2 = sonpy >=1.7.5;python_version=='3.7.*,3.8.*,3.9.*' acq = bioread >=1.0.5 -mat= +mat = pymatreader >=0.0.24 +interfaces = + %(acq)s + %(mat)s + %(spike2)s +gui = + peakdet >=0.1.0 duecredit = duecredit +all = + %(interfaces)s + %(gui)s + %(duecredit)s doc = sphinx >=2.0 sphinx-argparse @@ -54,33 +62,27 @@ style = isort <6.0.0 pydocstyle codespell -interfaces = - %(acq)s - %(mat)s - %(spike2)s test = pytest >=5.3 pytest-cov coverage - %(interfaces)s + %(all)s %(style)s -all = +dev = + %(all)s %(doc)s - %(duecredit)s - %(interfaces)s %(style)s %(test)s + pre-commit [options.package_data] abagen = phys2bids/data/* phys2bids/heuristics/* phys2bids/tests/data/* - [options.entry_points] console_scripts = phys2bids=phys2bids.phys2bids:_main - [flake8] doctest = True exclude = @@ -97,7 +99,6 @@ extend-ignore = E203, E501 extend-select = B950 per-file-ignores = workflow.py:D401 - [isort] profile = black skip_gitignore = true @@ -112,24 +113,20 @@ extend_skip = phys2bids/_version.py skip_glob = docs/* - [pydocstyle] convention = numpy match = nigsp/*.py match_dir = nigsp/[^tests,^heuristics]* - [codespell] skip = venvs,.venv,versioneer.py,.git,build,./docs/_build write-changes = count = quiet-level = 3 - [tool:pytest] doctest_optionflags = NORMALIZE_WHITESPACE xfail_strict = true addopts = -rx - [versioneer] VCS = git style = pep440 @@ -137,7 +134,6 @@ versionfile_source = phys2bids/_version.py versionfile_build = phys2bids/_version.py tag_prefix = parentdir_prefix = - [coverage:run] branch = True omit =