05_extract_signals_allmasks.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import os\n",
    "import re\n",
    "import shutil\n",
    "import random\n",
    "import string\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import nilearn\n",
    "from nilearn import plotting, image\n",
    "from nilearn.input_data import NiftiMasker\n",
    "import nibabel as nib\n",
    "from nipype.interfaces import ants\n",
    "import nighres\n",
    "\n",
    "import subprocess\n",
    "import json\n",
    "import multiprocessing as mp\n",
    "from functools import partial\n",
    "import joblib\n",
    "from joblib import Parallel, delayed\n",
    "import itertools\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_rois(sub, atlas_name='MASSP', space='T1w'):\n",
    "    # THALAMUS SUBREGIONS\n",
    "    if atlas_name=='THAL':\n",
    "        if space == 'MNI152NLin2009cAsym' or space == 'mni':\n",
    "            print('')\n",
    "        else:\n",
    "            mask_dir = f'../derivatives/masks_thal_func/sub-{sub}/anat/sub-{sub}_*.nii.gz'\n",
    "            fns = sorted(glob.glob(mask_dir))\n",
    "            names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_desc-mask-(?P<label>\\S+).nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "    # ATAG ATLAS        \n",
    "    elif atlas_name == 'ATAG':\n",
    "        if space == 'MNI152NLin2009cAsym' or space == 'mni':\n",
    "            ### Rois in MNI09c-space\n",
    "            mask_dir='/home/Public/trondheim/sourcedata/masks/MNI152NLin2009cAsym_res-1p5'\n",
    "            fns = sorted(glob.glob(mask_dir + '/space-*'))\n",
    "            names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_res-1p5_label-(?P<label>[a-zA-Z0-9]+)_probseg_def-img.nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "        else:\n",
    "            mask_dir = f'../derivatives/masks_atag_func/sub-{sub}/anat/sub-{sub}_*.nii.gz'\n",
    "            fns = sorted(glob.glob(mask_dir))\n",
    "            names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_desc-mask-(?P<label>[a-zA-Z0-9]+).nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "    # MASSP ATLAS        \n",
    "    elif atlas_name == 'MASSP':\n",
    "        mask_dir = f'../derivatives/masks_massp_func/sub-{sub}/anat/sub-{sub}_*.nii.gz'\n",
    "        fns = sorted(glob.glob(mask_dir))\n",
    "        names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_desc-mask-(?P<label>\\S+).nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "    # HARVARD OXFORD ATLAS\n",
    "    elif atlas_name == 'CORT':\n",
    "        mask_dir = f'../derivatives/masks_cortex_func/sub-{sub}/anat/sub-{sub}_*.nii.gz'\n",
    "        fns = sorted(glob.glob(mask_dir))\n",
    "        names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_desc-mask-(?P<label>\\S+).nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "    #Pauli atlas\n",
    "    elif atlas_name == 'Pauli':\n",
    "        mask_dir = f'../derivatives/masks_Pauli_func/sub-{sub}/anat/sub-{sub}_*.nii.gz'\n",
    "        fns = sorted(glob.glob(mask_dir))\n",
    "        names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_desc-mask-(?P<label>\\S+).nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "    #constructed FPN masks from Brodmann areas (Pijnenburg 2022)\n",
    "    elif atlas_name == 'FPN':\n",
    "        mask_dir = f'../derivatives/masks_FPN_func/sub-{sub}/anat/sub-{sub}_*.nii.gz'\n",
    "        fns = sorted(glob.glob(mask_dir))\n",
    "        names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_desc-mask-(?P<label>\\S+).nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "    elif atlas_name == 'WM-rep':\n",
    "        mask_dir = f'../derivatives/masks_WM-rep_func/sub-{sub}/anat/sub-{sub}_*.nii.gz'\n",
    "        fns = sorted(glob.glob(mask_dir))\n",
    "        names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_desc-mask-(?P<label>\\S+).nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "    elif atlas_name == 'HCP_MMP1':\n",
    "        mask_dir = f'../derivatives/masks_HCP_MMP1_func/sub-{sub}/anat/sub-{sub}_*.nii.gz'\n",
    "        fns = sorted(glob.glob(mask_dir))\n",
    "        names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_desc-mask-(?P<label>\\S+).nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "    elif atlas_name == 'str':\n",
    "        mask_dir = f'../derivatives/masks_str_func/sub-{sub}/anat/sub-{sub}_*.nii.gz'\n",
    "        fns = sorted(glob.glob(mask_dir))\n",
    "        names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_desc-mask-(?P<label>\\S+).nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "\n",
    "    roi_dict = dict(zip(names, fns))\n",
    "    return roi_dict\n",
    "\n",
    "def load_atlas(sub, atlas_name='MASSP', space='T1w'):\n",
    "    from nilearn import image\n",
    "    \n",
    "    roi_dict = find_rois(sub, atlas_name, space)\n",
    "    if len(roi_dict) == 0:\n",
    "        warnings.warn(f'No ROIs found for sub-{sub} atlas-{atlas_name} space-{space}. Returning 0 to prevent error')\n",
    "        return 0\n",
    "    combined = image.concat_imgs(roi_dict.values())\n",
    "    \n",
    "    class AttrDict(dict):\n",
    "        def __init__(self, *args, **kwargs):\n",
    "            super(AttrDict, self).__init__(*args, **kwargs)\n",
    "            self.__dict__ = self\n",
    "            \n",
    "    roi_atlas = AttrDict({'maps': combined,\n",
    "                          'labels': roi_dict.keys()})\n",
    "    \n",
    "    return roi_atlas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Extract signals from each ROI\n",
    "## Manual coded extraction - Slow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_epi(sub, ses, task, run, use_hp=False, base_dir='../derivatives/fmriprep/fmriprep'):\n",
    "    if use_hp:\n",
    "        epi = os.path.join('../derivatives/high_passed_func', f'sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_space-T1w_desc-preproc_bold.nii.gz')\n",
    "    else:\n",
    "        epi = os.path.join(base_dir, f'sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_space-T1w_desc-preproc_bold.nii.gz')\n",
    "    return epi\n",
    "\n",
    "def _make_psc(data):\n",
    "    mean_img = image.mean_img(data)\n",
    "\n",
    "    # Replace 0s for numerical reasons\n",
    "    mean_data = mean_img.get_fdata()\n",
    "    mean_data[mean_data == 0] = 1\n",
    "    denom = image.new_img_like(mean_img, mean_data)\n",
    "\n",
    "    return image.math_img('data / denom[..., np.newaxis] * 100 - 100',\n",
    "                          data=data, denom=denom)\n",
    "\n",
    "def do_extract(to_run, atlas_name='MASSP', overwrite=False, to_psc=False, use_hp=False):\n",
    "    sub, ses, task, run = to_run\n",
    "    sub = str(sub).zfill(3)\n",
    "    print(f'Extracting from sub-{sub}/ses-{ses}/sub-{sub}_ses-{ses}_task-{task}_run-{run}', end='')\n",
    "    \n",
    "    epi_fn = get_epi(sub,ses,task,run,use_hp)\n",
    "    if not os.path.exists(epi_fn):\n",
    "        print('...doesnt exist, skipping'.format(sub,ses,task,run))\n",
    "        return None\n",
    "    \n",
    "    ## dont really need to convert to psc here\n",
    "    if to_psc:\n",
    "        epi = _make_psc(epi_fn)\n",
    "        psc_fn = '_psc'\n",
    "    else:\n",
    "        epi = nib.load(epi_fn)\n",
    "        psc_fn = ''\n",
    "    \n",
    "    # might wanna ahve the hp data handy\n",
    "    if use_hp:\n",
    "        output_fn = f'../derivatives/extracted_signals/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-{atlas_name}-signals{psc_fn}_hp.tsv'\n",
    "    else:\n",
    "        output_fn = f'../derivatives/extracted_signals/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-{atlas_name}-signals{psc_fn}.tsv'\n",
    "    \n",
    "    if os.path.exists(output_fn) and not overwrite:\n",
    "        print(f'{output_fn} already run, loading previous result...')\n",
    "        return pd.read_csv(output_fn, sep='\\t')\n",
    "    \n",
    "    #load & reshpae\n",
    "    epi_flat = epi.get_fdata().reshape((np.product(epi.shape[:3]), epi.shape[-1]))\n",
    "\n",
    "    # load atlas\n",
    "    atlas = load_atlas(sub,atlas_name=atlas_name)\n",
    "    \n",
    "    dfs = []\n",
    "    for i in np.arange(len(atlas.labels)):\n",
    "        print('.', end='')\n",
    "        label = list(atlas.labels)[i]\n",
    "        mask = image.index_img(atlas.maps, i)\n",
    "        mask_flat = mask.get_fdata().ravel()\n",
    "        print(f'There are {np.count_nonzero(mask_flat)} voxels in region {label}')\n",
    "        if mask_flat.sum() == 0: # if there are no voxels in the mask then add one voxel so code doesn't crash\n",
    "            mask_flat[-1] = 1\n",
    "#         print(mask_flat)\n",
    "#         print(label)\n",
    "        print(len(epi_flat))\n",
    "        print(len(mask_flat))\n",
    "        signal = pd.DataFrame(np.average(epi_flat, weights=mask_flat, axis=0), columns=[label])\n",
    "        signal.index.name = 'volume'\n",
    "        dfs.append(signal)\n",
    "\n",
    "    df = pd.concat(dfs, axis=1)\n",
    "    if not os.path.exists(os.path.dirname(output_fn)):\n",
    "        os.makedirs(os.path.dirname(output_fn))\n",
    "    df.to_csv(output_fn, sep='\\t')\n",
    "    print(output_fn)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find all available functional runs, extract sub/ses/task/run info\n",
    "all_runs = sorted(glob.glob('../derivatives/fmriprep/fmriprep/sub-*/ses-*/func/*space-T1w*_bold.nii.gz'))\n",
    "regex = re.compile('.*sub-(?P<sub>\\d+)_ses-(?P<ses>\\S+)_task-(?P<task>\\S+)_run-(?P<run>\\d)_space-T1w*')\n",
    "all_combs = [tuple(regex.match(x).groupdict().values()) for x in all_runs]\n",
    "all_combs[-10:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# just extract MSIT\n",
    "all_runs = sorted(glob.glob('../derivatives/fmriprep/fmriprep/sub-*/ses-sstmsit/func/*task-*space-T1w*_bold.nii.gz'))\n",
    "regex = re.compile('.*sub-(?P<sub>\\d+)_ses-(?P<ses>\\S+)_task-(?P<task>\\S+)_run-(?P<run>\\d)_space-T1w*')\n",
    "all_combs = [tuple(regex.match(x).groupdict().values()) for x in all_runs]\n",
    "# all_combs = [x for x in all_combs if not '015' in x] # bad hp data for sub 15????\n",
    "all_combs = [x for x in all_combs if '041' in x] # bad hp data for sub 26????\n",
    "# all_combs = [x for x in all_combs if x[0] in ['002','003','004','005','006','007','008','009','010','011']]\n",
    "all_combs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# just extract SST\n",
    "all_runs = sorted(glob.glob('../derivatives/fmriprep/fmriprep/sub-*/ses-sstmsit/func/*task-msit*space-T1w*_bold.nii.gz'))\n",
    "regex = re.compile('.*sub-(?P<sub>\\d+)_ses-(?P<ses>\\S+)_task-(?P<task>\\S+)_run-(?P<run>\\d)_space-T1w*')\n",
    "all_combs = [tuple(regex.match(x).groupdict().values()) for x in all_runs]\n",
    "# all_combs = [x for x in all_combs if not '015' in x] # bad hp data for sub 15????\n",
    "# all_combs = [x for x in all_combs if '041' in x] # bad hp data for sub 26????\n",
    "all_combs = [x for x in all_combs if x[0] in ['004','008','010','013','019','027']]\n",
    "all_combs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# just extract RBREVL\n",
    "all_runs = sorted(glob.glob('../derivatives/fmriprep/fmriprep/sub-*/ses-rbrevl/func/*task-rb*space-T1w*_bold.nii.gz'))\n",
    "regex = re.compile('.*sub-(?P<sub>\\d+)_ses-(?P<ses>\\S+)_task-(?P<task>\\S+)_run-(?P<run>\\d)_space-T1w*')\n",
    "all_combs = [tuple(regex.match(x).groupdict().values()) for x in all_runs]\n",
    "# all_combs = [x for x in all_combs if not '015' in x] # bad hp data for sub 15????\n",
    "all_combs = [x for x in all_combs]# if '041' in x] # bad hp data for sub 26????\n",
    "# all_combs = [x for x in all_combs if x[0] in ['002','003','004','005','006','007','008','009','010','011']]\n",
    "all_combs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# just extract MSIT\n",
    "all_runs = sorted(glob.glob('../derivatives/fmriprep/fmriprep/sub-*/ses-rlsat/func/*task-*space-T1w*_bold.nii.gz'))\n",
    "regex = re.compile('.*sub-(?P<sub>\\d+)_ses-(?P<ses>\\S+)_task-(?P<task>\\S+)_run-(?P<run>\\d)_space-T1w*')\n",
    "all_combs = [tuple(regex.match(x).groupdict().values()) for x in all_runs]\n",
    "all_combs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all_subs = np.arange(2,28)\n",
    "# all_ses = ['rlsat', 'rbrevl', 'anatomical', 'sstmsit']\n",
    "# all_tasks = ['rs', 'rlsat', 'rb', 'revl', 'sst', 'msit']\n",
    "# all_runs = [1,2,3]\n",
    "\n",
    "# all_combs = list(itertools.product(all_subs,all_ses,all_tasks,all_runs))\n",
    "# all_combs = [x for x in all_combs if (x[1]=='rlsat' and x[2]=='rlsat') or (x[1]=='rbrevl' and x[2] in ['rb', 'revl'] and x[3]<3) or (x[1]=='sstmsit' and x[2] in ['sst', 'msit'] and x[3]<3) or (x[1]=='anatomical' and x[2]=='rs' and x[3]<3)]\n",
    "# #do_extract(all_combs[0], overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_affines(sub):\n",
    "    sub = str(sub).zfill(3)\n",
    "    all_funcs = sorted(glob.glob(f'../derivatives/fmriprep/fmriprep/sub-{sub}/ses*/func/sub*_space-T1w_desc-preproc_bold.nii.gz'))\n",
    "    all_affines = [nib.load(x).affine for x in all_funcs]\n",
    "    return (np.array(all_affines)[0] == np.array(all_affines)).all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#all_atlases=['MASSP','CORT','ATAG']#,'THAL'] #['Pauli']\n",
    "# all_atlases=['HCP_MMP1']\n",
    "# all_atlases = ['CORT']\n",
    "all_atlases = ['str','MASSP','CORT','ATAG']\n",
    "\n",
    "hp_options= [True,False]\n",
    "overwrite=False\n",
    "psc=False\n",
    "\n",
    "for atlas_name in all_atlases:\n",
    "    for hp in hp_options:\n",
    "        for i, comb in enumerate(all_combs):\n",
    "            print(f'atlas-{atlas_name} hp-{hp}')\n",
    "            print(comb)\n",
    "            sub = comb[0]\n",
    "            if check_affines(sub):\n",
    "                do_extract(comb, atlas_name=atlas_name, overwrite=overwrite, to_psc=psc, use_hp=hp)\n",
    "            else:\n",
    "                print(f'Affines for sub {sub} not identical')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#all_atlases=['MASSP','CORT','ATAG']#,'THAL'] #['Pauli']\n",
    "all_atlases=['THAL']\n",
    "\n",
    "hp_options= [False]#[True,False]\n",
    "overwrite=False\n",
    "psc=False\n",
    "\n",
    "for atlas_name in all_atlases:\n",
    "    for hp in hp_options:\n",
    "        for i, comb in enumerate(all_combs):\n",
    "            print(f'atlas-{atlas_name} hp-{hp}')\n",
    "            print(comb)\n",
    "            sub = comb[0]\n",
    "            if check_affines(sub):\n",
    "                do_extract(comb, atlas_name=atlas_name, overwrite=overwrite, to_psc=psc, use_hp=hp)\n",
    "            else:\n",
    "                print(f'Affines for sub {sub} not identical')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#all_atlases=['MASSP','CORT','ATAG']#,'THAL'] #['Pauli']\n",
    "all_atlases=['WM-rep']\n",
    "\n",
    "hp_options= [False]#[True,False]\n",
    "overwrite=False\n",
    "psc=False\n",
    "\n",
    "for atlas_name in all_atlases:\n",
    "    for hp in hp_options:\n",
    "        for i, comb in enumerate(all_combs):\n",
    "            print(f'atlas-{atlas_name} hp-{hp}')\n",
    "            print(comb)\n",
    "            sub = comb[0]\n",
    "            print(sub)\n",
    "            if check_affines(sub):\n",
    "                do_extract(comb, atlas_name=atlas_name, overwrite=overwrite, to_psc=psc, use_hp=hp)\n",
    "            else:\n",
    "                print(f'Affines for sub {sub} not identical')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_funcs = sorted(glob.glob(f'../derivatives/fmriprep/fmriprep/sub-026/ses*/func/sub*_space-T1w_desc-preproc_bold.nii.gz'))\n",
    "all_funcs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "check_affines('026')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2 extract whole roi signal from thalamus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_rois(sub, atlas_name='MASSP', space='T1w'):\n",
    "    if atlas_name=='THAL':\n",
    "        if space == 'MNI152NLin2009cAsym' or space == 'mni':\n",
    "            print('')\n",
    "        else:\n",
    "            mask_dir = f'../derivatives/masks_thal_func/sub-{sub}/anat/sub-{sub}_*.nii.gz'\n",
    "            fns = sorted(glob.glob(mask_dir))\n",
    "            names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_desc-mask-(?P<label>\\S+).nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "            \n",
    "    if atlas_name == 'ATAG':\n",
    "        if space == 'MNI152NLin2009cAsym' or space == 'mni':\n",
    "            ### Rois in MNI09c-space\n",
    "            mask_dir='/home/Public/trondheim/sourcedata/masks/MNI152NLin2009cAsym_res-1p5'\n",
    "            fns = sorted(glob.glob(mask_dir + '/space-*'))\n",
    "            names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_res-1p5_label-(?P<label>[a-zA-Z0-9]+)_probseg_def-img.nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "        else:\n",
    "            mask_dir = f'../derivatives/masks_atag_func/sub-{sub}/anat/sub-{sub}_*.nii.gz'\n",
    "            fns = sorted(glob.glob(mask_dir))\n",
    "            names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_desc-mask-(?P<label>[a-zA-Z0-9]+).nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "            \n",
    "    elif atlas_name == 'MASSP':\n",
    "        mask_dir = f'../derivatives/masks_massp_func/sub-{sub}/anat/sub-{sub}_*Tha*.nii.gz' # only thalamus\n",
    "        fns = sorted(glob.glob(mask_dir))\n",
    "        names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_desc-mask-(?P<label>\\S+).nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "\n",
    "    roi_dict = dict(zip(names, fns))\n",
    "    return roi_dict\n",
    "\n",
    "def load_atlas(sub, atlas_name='MASSP', space='T1w'):\n",
    "    from nilearn import image\n",
    "    \n",
    "    roi_dict = find_rois(sub, atlas_name, space)\n",
    "    combined = image.concat_imgs(roi_dict.values())\n",
    "    \n",
    "    class AttrDict(dict):\n",
    "        def __init__(self, *args, **kwargs):\n",
    "            super(AttrDict, self).__init__(*args, **kwargs)\n",
    "            self.__dict__ = self\n",
    "            \n",
    "    roi_atlas = AttrDict({'maps': combined,\n",
    "                          'labels': roi_dict.keys()})\n",
    "    \n",
    "    return roi_atlas\n",
    "\n",
    "def get_epi(sub, ses, task, run, use_hp=False, base_dir='../derivatives/fmriprep/fmriprep'):\n",
    "    if use_hp:\n",
    "        epi = os.path.join('../derivatives/high_passed_func', f'sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_space-T1w_desc-preproc_bold.nii.gz')\n",
    "    else:\n",
    "        epi = os.path.join(base_dir, f'sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_space-T1w_desc-preproc_bold.nii.gz')\n",
    "    return epi\n",
    "\n",
    "def _make_psc(data):\n",
    "    mean_img = image.mean_img(data)\n",
    "\n",
    "    # Replace 0s for numerical reasons\n",
    "    mean_data = mean_img.get_fdata()\n",
    "    mean_data[mean_data == 0] = 1\n",
    "    denom = image.new_img_like(mean_img, mean_data)\n",
    "\n",
    "    return image.math_img('data / denom[..., np.newaxis] * 100 - 100',\n",
    "                          data=data, denom=denom)\n",
    "\n",
    "def do_extract(to_run, atlas='MASSP', overwrite=False, to_psc=False, use_hp=False):\n",
    "    sub, ses, task, run = to_run\n",
    "    sub = str(sub).zfill(3)\n",
    "    print(f'Extracting from sub-{sub}/ses-{ses}/sub-{sub}_ses-{ses}_task-{task}_run-{run}', end='')\n",
    "    \n",
    "    epi_fn = get_epi(sub,ses,task,run,use_hp)\n",
    "    if not os.path.exists(epi_fn):\n",
    "        print('...doesnt exist, skipping'.format(sub,ses,task,run))\n",
    "        return None\n",
    "\n",
    "    if atlas == 'thal':\n",
    "        toappend = '_thalamus'\n",
    "    else: \n",
    "        toappend=''\n",
    "    \n",
    "    ## dont really need to convert to psc here\n",
    "    if to_psc:\n",
    "        epi = _make_psc(epi_fn)\n",
    "        psc_fn = '_psc'\n",
    "    else:\n",
    "        epi = nib.load(epi_fn)\n",
    "        psc_fn = ''\n",
    "    \n",
    "    # might wanna ahve the hp data handy\n",
    "    if use_hp:\n",
    "        output_fn = f'../derivatives/extracted_signals_thal_voxels/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-{atlas}-signals{psc_fn}_hp.tsv'\n",
    "    else:\n",
    "        output_fn = f'../derivatives/extracted_signals_that_voxels/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-{atlas}-signals{psc_fn}.tsv'\n",
    "    \n",
    "    if os.path.exists(output_fn) and not overwrite:\n",
    "        print(f'{output_fn} already run, loading previous result...')\n",
    "        return pd.read_csv(output_fn, sep='\\t')\n",
    "    \n",
    "    #epi = nib.load(epi_fn)\n",
    "    epi_flat = epi.get_fdata().reshape((np.product(epi.shape[:3]), epi.shape[-1]))\n",
    "\n",
    "    # load atlas\n",
    "    atlas = load_atlas(sub,atlas_name=atlas)\n",
    "    \n",
    "    dfs = []\n",
    "    for i in np.arange(len(atlas.labels)):\n",
    "        print('.', end='')\n",
    "        label = list(atlas.labels)[i]\n",
    "        mask = image.index_img(atlas.maps, i)\n",
    "        mask_flat = mask.get_fdata().ravel()\n",
    "        indexes = np.where(mask_flat>0)[0] # indexes of voxel within mask. len(indexes) and np.count_nonzero(mask_flat) should be the same\n",
    "\n",
    "        print(f'There are {np.count_nonzero(mask_flat)} voxels in region {label} for sub {sub}')\n",
    "        for label_n, inds in enumerate(indexes):\n",
    "            mask_flat_voxel = mask_flat.copy()\n",
    "            mask_flat_voxel[:] = 0\n",
    "            mask_flat_voxel[inds] = 1\n",
    "            signal = pd.DataFrame(np.average(epi_flat, weights=mask_flat_voxel, axis=0), columns=[label+'_'+str(label_n)])\n",
    "            signal.index.name = 'volume'\n",
    "            dfs.append(signal)\n",
    "\n",
    "    df = pd.concat(dfs, axis=1)\n",
    "    if not os.path.exists(os.path.dirname(output_fn)):\n",
    "        os.makedirs(os.path.dirname(output_fn))\n",
    "    df.to_csv(output_fn, sep='\\t')\n",
    "    print(output_fn)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, comb in enumerate(all_combs):\n",
    "    print(comb)\n",
    "    sub = comb[0]\n",
    "    if check_affines(sub):\n",
    "        do_extract(comb, atlas='MASSP', overwrite=True, to_psc=False, use_hp=True)\n",
    "    else:\n",
    "        print(f'Affines for sub {sub} not identical')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, comb in enumerate(all_combs):\n",
    "    print(comb)\n",
    "    sub = comb[0]\n",
    "    if check_affines(sub):\n",
    "        do_extract(comb, atlas='THAL', overwrite=True, to_psc=False, use_hp=False)\n",
    "    else:\n",
    "        print(f'Affines for sub {sub} not identical')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "atlas = 'MASSP'\n",
    "to_run = [('002', 'sstmsit', 'msit', '1')]\n",
    "use_hp=True\n",
    "sub='002'\n",
    "ses='sstmsit'\n",
    "task='msit'\n",
    "run='1'\n",
    "to_psc = False\n",
    "overwrite=True\n",
    "sub = str(sub).zfill(3)\n",
    "print(f'Extracting from sub-{sub}/ses-{ses}/sub-{sub}_ses-{ses}_task-{task}_run-{run}', end='')\n",
    "\n",
    "epi_fn = get_epi(sub,ses,task,run,use_hp)\n",
    "if not os.path.exists(epi_fn):\n",
    "    print('...doesnt exist, skipping'.format(sub,ses,task,run))\n",
    "#     return None\n",
    "\n",
    "if atlas == 'thal':\n",
    "    toappend = '_thalamus'\n",
    "else: \n",
    "    toappend=''\n",
    "\n",
    "## dont really need to convert to psc here\n",
    "if to_psc:\n",
    "    epi = _make_psc(epi_fn)\n",
    "    psc_fn = '_psc'\n",
    "else:\n",
    "    epi = nib.load(epi_fn)\n",
    "    psc_fn = ''\n",
    "\n",
    "# might wanna ahve the hp data handy\n",
    "if use_hp:\n",
    "    output_fn = f'../derivatives/extracted_signals_thal_voxels/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-{atlas}-signals{psc_fn}_hp.tsv'\n",
    "else:\n",
    "    output_fn = f'../derivatives/extracted_signals_that_voxels/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-{atlas}-signals{psc_fn}.tsv'\n",
    "\n",
    "if os.path.exists(output_fn) and not overwrite:\n",
    "    print(f'{output_fn} already run, loading previous result...')\n",
    "#     return pd.read_csv(output_fn, sep='\\t')\n",
    "\n",
    "#epi = nib.load(epi_fn)\n",
    "epi_flat = epi.get_fdata().reshape((np.product(epi.shape[:3]), epi.shape[-1]))\n",
    "\n",
    "# load atlas\n",
    "atlas = load_atlas(sub,atlas_name=atlas)\n",
    "\n",
    "dfs = []\n",
    "for i in np.arange(len(atlas.labels)):\n",
    "    print('.', end='')\n",
    "    label = list(atlas.labels)[i]\n",
    "    mask = image.index_img(atlas.maps, i)\n",
    "    mask_flat = mask.get_fdata().ravel()\n",
    "    indexes = np.where(mask_flat>0)[0] # indexes of voxel within mask. len(indexes) and np.count_nonzero(mask_flat) should be the same\n",
    "    \n",
    "    print(f'There are {np.count_nonzero(mask_flat)} voxels in region {label} for sub {sub}')\n",
    "#     for j in np.arange(np.count_nonzero(mask_flat)): ### FINISH THISSS # loop over each voxel in mask .. \n",
    "    for i, kk in enumerate(indexes):\n",
    "        mask_flat_voxel = mask_flat.copy()\n",
    "        mask_flat_voxel[:] = 0\n",
    "        mask_flat_voxel[kk] = 1\n",
    "        signal = pd.DataFrame(np.average(epi_flat, weights=mask_flat_voxel, axis=0), columns=[label+'_'+str(j)])\n",
    "        signal.index.name = 'volume'\n",
    "        dfs.append(signal)\n",
    "\n",
    "# df = pd.concat(dfs, axis=1)\n",
    "# #     output_fn = f'../derivatives/extracted_signals/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-MASSP-signals.tsv'\n",
    "# if not os.path.exists(os.path.dirname(output_fn)):\n",
    "#     os.makedirs(os.path.dirname(output_fn))\n",
    "# df.to_csv(output_fn, sep='\\t')\n",
    "# print(output_fn)\n",
    "# return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Use nilearn, and high-pass filter & remove confounds along the way\n",
    "### This is much faster, but something very weird happens when multiple atlas maps overlap (eg when there's both a mask for \"M1\" as well as \"rM1\" and \"lM1\")\n",
    "the M1 case isn't very troubling, but it suggests something funny happens with overlapping maps/masks - do we have overlap? Perhaps the manually coded version is safer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## extract signals this way?\n",
    "## stolen from nideconv\n",
    "import pandas as pd\n",
    "from nilearn import input_data\n",
    "import nibabel as nb\n",
    "from nilearn._utils import check_niimg\n",
    "from nilearn import image\n",
    "import numpy as np\n",
    "\n",
    "def extract_timecourse_from_nii(atlas,\n",
    "                                nii,\n",
    "                                mask=None,\n",
    "                                confounds=None,\n",
    "                                atlas_type=None,\n",
    "                                t_r=None,\n",
    "                                low_pass=None,\n",
    "                                high_pass=1./128,\n",
    "                                to_psc=False,\n",
    "                                *args,\n",
    "                                **kwargs):\n",
    "\n",
    "\n",
    "    standardize = kwargs.pop('standardize', False)\n",
    "    detrend = kwargs.pop('detrend', False)\n",
    "\n",
    "    if atlas_type is None:\n",
    "        maps = check_niimg(atlas.maps)\n",
    "\n",
    "        if len(maps.shape) == 3:\n",
    "            atlas_type = 'labels'\n",
    "        else:\n",
    "            atlas_type = 'prob'\n",
    "\n",
    "    if atlas_type == 'labels':\n",
    "        masker = input_data.NiftiLabelsMasker(atlas.maps,\n",
    "                                              mask_img=mask,\n",
    "                                              standardize=standardize,\n",
    "                                              detrend=detrend,\n",
    "                                              t_r=t_r,\n",
    "                                              low_pass=low_pass,\n",
    "                                              high_pass=high_pass,\n",
    "                                              *args, **kwargs)\n",
    "    else:\n",
    "        masker = input_data.NiftiMapsMasker(atlas.maps,\n",
    "                                            mask_img=mask,\n",
    "                                            standardize=standardize,\n",
    "                                            detrend=detrend,\n",
    "                                            t_r=t_r,\n",
    "                                            low_pass=low_pass,\n",
    "                                            high_pass=high_pass,\n",
    "                                            *args, **kwargs)\n",
    "\n",
    "    if to_psc:\n",
    "        data = _make_psc(nii)\n",
    "    else:\n",
    "        data = nii\n",
    "\n",
    "    results = masker.fit_transform(data,\n",
    "                                   confounds=confounds)\n",
    "\n",
    "    # For weird atlases that have a label for the background\n",
    "    if len(atlas.labels) == results.shape[1] + 1:\n",
    "        atlas.labels = atlas.labels[1:]\n",
    "\n",
    "    if t_r is None:\n",
    "        t_r = 1\n",
    "    print(t_r)\n",
    "    index = pd.Index(np.arange(0,\n",
    "                               t_r*nib.load(nii).shape[-1],\n",
    "                               t_r),\n",
    "                     name='time')\n",
    "\n",
    "    columns = pd.Index(atlas.labels,\n",
    "                       name='roi')\n",
    "\n",
    "    return pd.DataFrame(results,\n",
    "                        index=index,\n",
    "                        columns=columns)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exclude_map_from_atlas(atlas, map_name):\n",
    "    \n",
    "    indx = np.where(np.array(list(atlas.labels)) == map_name)[0][0]\n",
    "\n",
    "    all_indices = np.arange(atlas.maps.shape[-1])\n",
    "    indices = [x for x in all_indices if not x == indx]\n",
    "\n",
    "    atlas.maps = nilearn.image.index_img(atlas.maps, indices)\n",
    "    atlas.labels = np.array(list(atlas.labels))[indices].tolist()\n",
    "    \n",
    "    return atlas\n",
    "\n",
    "def _make_psc(data):\n",
    "    mean_img = image.mean_img(data)\n",
    "\n",
    "    # Replace 0s for numerical reasons\n",
    "    mean_data = mean_img.get_data()\n",
    "    mean_data[mean_data == 0] = 1\n",
    "    denom = image.new_img_like(mean_img, mean_data)\n",
    "\n",
    "    return image.math_img('data / denom[..., np.newaxis] * 100 - 100',\n",
    "                          data=data, denom=denom)\n",
    "\n",
    "\n",
    "def extract_signals_nilearn(comb, include_physio=True, space='T1w', overwrite=False):\n",
    "#for sub, ses, task, run in all_combs:\n",
    "    sub,ses,task,run = comb\n",
    "    epi_fn = get_epi(sub,ses,task,run)\n",
    "    if not os.path.exists(epi_fn):\n",
    "        print('...doesnt exist, skipping'.format(sub,ses,task,run))       \n",
    "        return None\n",
    "    \n",
    "    # load confounds\n",
    "    confounds_fn = f'../derivatives/fmriprep/fmriprep/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-confounds_timeseries.tsv'\n",
    "    confounds = pd.read_csv(confounds_fn, sep='\\t')[['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z', 'dvars', 'framewise_displacement']].fillna(method='bfill')\n",
    "\n",
    "    # get retroicor\n",
    "    if include_physio:\n",
    "        retroicor_fn = f'../derivatives/retroicor/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-retroicor_regressors.tsv'\n",
    "        if not os.path.exists(retroicor_fn):\n",
    "            ## take first 20 aCompCor components\n",
    "            print(\"No retroicor found, including 20 a_comp_cor components\")\n",
    "            a_comp_cor = pd.read_csv(confounds_fn, sep='\\t')[['a_comp_cor_' + str(x).zfill(2) for x in range(20)]]\n",
    "            confounds = pd.concat([confounds, a_comp_cor], axis=1)\n",
    "        else:\n",
    "            retroicor = pd.read_csv(retroicor_fn, sep='\\t', header=None).iloc[:,:20]  ## 20 components in total\n",
    "            retroicor.columns = ['cardiac_' + str(x) for x in range(6)] + ['respiratory_' + str(x) for x in range(8)] + ['respiratoryxcardiac_' + str(x) for x in range(4)] + ['HRV', 'RVT']\n",
    "            confounds = pd.concat([confounds, retroicor], axis=1)\n",
    "\n",
    "    # get brain mask\n",
    "    brain_mask = f'../derivatives/fmriprep/fmriprep/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-1_space-{space}_desc-brain_mask.nii.gz'\n",
    "    \n",
    "    for atlas_type in ['MASSP', 'ATAG']:\n",
    "#     for atlas_type in ['ATAG']:\n",
    "        output_fn = f'../derivatives/extracted_signals_nilearn/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-{atlas_type}-signals.tsv'\n",
    "        if os.path.exists(output_fn) and not overwrite:\n",
    "            return 0\n",
    "        \n",
    "        subject_atlas = load_atlas(sub, atlas_name=atlas_type)\n",
    "        \n",
    "        if atlas_type == 'ATAG':\n",
    "            subject_atlas = exclude_map_from_atlas(subject_atlas, 'M1')\n",
    "        \n",
    "        df = extract_timecourse_from_nii(subject_atlas, epi_fn, mask=brain_mask, confounds=confounds, high_pass=1/128., t_r=1.38, to_psc=True)\n",
    "        if not os.path.exists(os.path.dirname(output_fn)):\n",
    "            os.makedirs(os.path.dirname(output_fn))\n",
    "        df.to_csv(output_fn, sep='\\t')\n",
    "\n",
    "    return 0\n",
    "#         print(output_fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include_physio = True\n",
    "space = 'T1w'\n",
    "\n",
    "all_runs = sorted(glob.glob('../derivatives/fmriprep/fmriprep/sub-*/ses-*/func/*space-T1w*_bold.nii.gz'))\n",
    "regex = re.compile('.*sub-(?P<sub>\\d+)_ses-(?P<ses>\\S+)_task-(?P<task>\\S+)_run-(?P<run>\\d)_space-T1w*')\n",
    "all_combs = [tuple(regex.match(x).groupdict().values()) for x in all_runs]\n",
    "\n",
    "\n",
    "all_combs = [x for x in all_combs if x[2] == 'msit']\n",
    "all_combs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = joblib.Parallel(n_jobs=10, verbose=True)(joblib.delayed(extract_signals_nilearn)(x, overwrite=True) for x in all_combs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Use pre-cleaned data, don't extract confounds during the process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include_physio = True\n",
    "space = 'T1w'\n",
    "\n",
    "def exclude_map_from_atlas(atlas, map_name):\n",
    "    \n",
    "    indx = np.where(np.array(list(atlas.labels)) == map_name)[0][0]\n",
    "\n",
    "    all_indices = np.arange(atlas.maps.shape[-1])\n",
    "    indices = [x for x in all_indices if not x == indx]\n",
    "\n",
    "    atlas.maps = nilearn.image.index_img(atlas.maps, indices)\n",
    "    atlas.labels = np.array(list(atlas.labels))[indices].tolist()\n",
    "    \n",
    "    return atlas\n",
    "\n",
    "def _make_psc(data):\n",
    "    mean_img = image.mean_img(data)\n",
    "\n",
    "    # Replace 0s for numerical reasons\n",
    "    mean_data = mean_img.get_data()\n",
    "    mean_data[mean_data == 0] = 1\n",
    "    denom = image.new_img_like(mean_img, mean_data)\n",
    "\n",
    "    return image.math_img('data / denom[..., np.newaxis] * 100 - 100',\n",
    "                          data=data, denom=denom)\n",
    "\n",
    "# def get_epi_fn(sub,ses,task,run, base_dir='../derivatives/fmriprep/fmriprep')\n",
    "\n",
    "def extract_signals_nilearn(comb, include_physio=True, space='T1w', overwrite=False, use_precleaned=False, use_confounds=True):\n",
    "    if use_precleaned and use_confounds:\n",
    "        raise(IOError('Cannot both use precleaned data AND clean, that''s a stupid idea!'))\n",
    "    \n",
    "    sub,ses,task,run = comb\n",
    "    output_base_dir = '../derivatives/extracted_signals_nilearn'\n",
    "    if use_precleaned:\n",
    "        base_dir = '../derivatives/cleaned_func'\n",
    "        output_base_dir += '_precleaned'\n",
    "    else:\n",
    "        base_dir = '../derivatives/fmriprep/fmriprep'\n",
    "    epi_fn = get_epi(sub,ses,task,run, base_dir=base_dir)\n",
    "    if not os.path.exists(epi_fn):\n",
    "        print('...doesnt exist, skipping'.format(sub,ses,task,run))       \n",
    "        return None\n",
    "    \n",
    "    if use_confounds:\n",
    "        output_base_dir += '_cleaned'\n",
    "        # load confounds\n",
    "        confounds_fn = f'../derivatives/fmriprep/fmriprep/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-confounds_timeseries.tsv'\n",
    "        confounds = pd.read_csv(confounds_fn, sep='\\t')[['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z', 'dvars', 'framewise_displacement']].fillna(method='bfill')\n",
    "\n",
    "        # get retroicor\n",
    "        if include_physio:\n",
    "            retroicor_fn = f'../derivatives/retroicor/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-retroicor_regressors.tsv'\n",
    "            if not os.path.exists(retroicor_fn):\n",
    "                ## take first 20 aCompCor components\n",
    "                print(\"No retroicor found, including 20 a_comp_cor components\")\n",
    "                a_comp_cor = pd.read_csv(confounds_fn, sep='\\t')[['a_comp_cor_' + str(x).zfill(2) for x in range(20)]]\n",
    "                confounds = pd.concat([confounds, a_comp_cor], axis=1)\n",
    "            else:\n",
    "                retroicor = pd.read_csv(retroicor_fn, sep='\\t', header=None).iloc[:,:20]  ## 20 components in total\n",
    "                retroicor.columns = ['cardiac_' + str(x) for x in range(6)] + ['respiratory_' + str(x) for x in range(8)] + ['respiratoryxcardiac_' + str(x) for x in range(4)] + ['HRV', 'RVT']\n",
    "                confounds = pd.concat([confounds, retroicor], axis=1)\n",
    "\n",
    "    # get brain mask\n",
    "    brain_mask = f'../derivatives/fmriprep/fmriprep/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-1_space-{space}_desc-brain_mask.nii.gz'\n",
    "    \n",
    "    for atlas_type in ['MASSP', 'ATAG']:\n",
    "        output_fn = os.path.join(output_base_dir, f'sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-{atlas_type}-signals.tsv')\n",
    "        if os.path.exists(output_fn) and not overwrite:\n",
    "            return 0\n",
    "        \n",
    "        subject_atlas = load_atlas(sub, atlas_name=atlas_type)\n",
    "        if atlas_type == 'ATAG':\n",
    "            subject_atlas = exclude_map_from_atlas(subject_atlas, 'M1')\n",
    "        \n",
    "        df = extract_timecourse_from_nii(subject_atlas, epi_fn, mask=brain_mask, high_pass=None)\n",
    "        if not os.path.exists(os.path.dirname(output_fn)):\n",
    "            os.makedirs(os.path.dirname(output_fn))\n",
    "        df.to_csv(output_fn, sep='\\t')\n",
    "\n",
    "    return 0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_runs = sorted(glob.glob('../derivatives/fmriprep/fmriprep/sub-*/ses-*/func/*space-T1w*_bold.nii.gz'))\n",
    "regex = re.compile('.*sub-(?P<sub>\\d+)_ses-(?P<ses>\\S+)_task-(?P<task>\\S+)_run-(?P<run>\\d)_space-T1w*')\n",
    "all_combs = [tuple(regex.match(x).groupdict().values()) for x in all_runs]\n",
    "\n",
    "all_combs = [x for x in all_combs if x[1] == 'rlsat']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = joblib.Parallel(n_jobs=10, verbose=True)(joblib.delayed(extract_signals_nilearn)(x, overwrite=False, use_precleaned=True, use_confounds=False) for x in all_combs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clean niftis\n",
    "1. High-pass\n",
    "2. Remove confounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include_physio=True\n",
    "def high_pass(nii, verbose=False, mask=None):\n",
    "    print('Highpass-filtering')\n",
    "    t_r = nii.header['pixdim'][4]\n",
    "    if mask is not None:\n",
    "        hp_masker = NiftiMasker(mask, high_pass=1./128, t_r=t_r)\n",
    "    else:\n",
    "        hp_masker = NiftiMasker(high_pass=1./128, t_r=t_r)\n",
    "    \n",
    "    # Generate & fit NiftiMasker\n",
    "    hp_data = hp_masker.fit_transform(nii)\n",
    "    \n",
    "    # back to brain space\n",
    "    inver = hp_masker.inverse_transform(hp_data)\n",
    "\n",
    "    # add mean of timeseries per voxel back\n",
    "    highpassed_data = inver.get_fdata() + np.mean(nii.get_fdata(), 3)[:,:,:,np.newaxis]\n",
    "    highpassed_img = nib.Nifti1Image(highpassed_data, inver.affine, header=nii.header)\n",
    "    \n",
    "    return highpassed_img\n",
    "\n",
    "\n",
    "def do_high_pass(fn, overwrite=False):\n",
    "    regex = re.compile('.*sub-(?P<sub>\\d+)_ses-(?P<ses>\\S+)_task-(?P<task>\\S+)_run-(?P<run>\\d)_space-T1w_desc-preproc_bold.*')\n",
    "    gd = regex.match(fn).groupdict()\n",
    "\n",
    "    brain_mask = nib.load(fn.replace('preproc_bold', 'brain_mask'))\n",
    "\n",
    "    # has this file been highpassed?\n",
    "    hp_save_fn = fn.replace('fmriprep/fmriprep', 'high_passed_func')\n",
    "    if os.path.exists(hp_save_fn) and not overwrite:\n",
    "        hp_data = nib.load(hp_save_fn)\n",
    "    else:\n",
    "        nii = nib.load(fn)\n",
    "        hp_data = high_pass(nii, mask=brain_mask)\n",
    "        os.makedirs(os.path.dirname(hp_save_fn), exist_ok=True)\n",
    "        hp_data.to_filename('../derivatives/high_passed_func/sub-{}/ses-{}/func/{}'.format(gd['sub'], gd['ses'], os.path.basename(fn)))\n",
    "\n",
    "        \n",
    "def do_clean(fn, overwrite=False):\n",
    "    regex = re.compile('.*sub-(?P<sub>\\d+)_ses-(?P<ses>\\S+)_task-(?P<task>\\S+)_run-(?P<run>\\d)_space-T1w_desc-preproc_bold.*')\n",
    "    gd = regex.match(fn).groupdict()\n",
    "    sub, ses, task, run = gd['sub'], gd['ses'], gd['task'], gd['run']\n",
    "    \n",
    "    # has this file been cleaned?\n",
    "    cleaned_save_fn = fn.replace('high_passed_func', 'cleaned_func')\n",
    "    if os.path.exists(cleaned_save_fn) and not overwrite:\n",
    "        print(cleaned_save_fn)\n",
    "        cleaned_data = nib.load(cleaned_save_fn)\n",
    "        print('eh')\n",
    "    else:\n",
    "        nii = nib.load(fn)\n",
    "        # load confounds\n",
    "        confounds_fn = f'../derivatives/fmriprep/fmriprep/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-confounds_timeseries.tsv'\n",
    "        confounds = pd.read_csv(confounds_fn, sep='\\t')[['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z', 'dvars', 'framewise_displacement']].fillna(method='bfill')\n",
    "\n",
    "        # get retroicor\n",
    "        if include_physio:\n",
    "            retroicor_fn = f'../derivatives/retroicor/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-retroicor_regressors.tsv'\n",
    "            if not os.path.exists(retroicor_fn):\n",
    "                ## take first 20 aCompCor components\n",
    "                print(\"No retroicor found, including 20 a_comp_cor components\")\n",
    "                a_comp_cor = pd.read_csv(confounds_fn, sep='\\t')[['a_comp_cor_' + str(x).zfill(2) for x in range(20)]]\n",
    "                confounds = pd.concat([confounds, a_comp_cor], axis=1)\n",
    "            else:\n",
    "                retroicor = pd.read_csv(retroicor_fn, sep='\\t', header=None).iloc[:,:20]  ## 20 components in total\n",
    "                retroicor.columns = ['cardiac_' + str(x) for x in range(6)] + ['respiratory_' + str(x) for x in range(8)] + ['respiratoryxcardiac_' + str(x) for x in range(4)] + ['HRV', 'RVT']\n",
    "                confounds = pd.concat([confounds, retroicor], axis=1)\n",
    "\n",
    "        cleaned_data = image.clean_img(nii, confounds=confounds, standardize=False, detrend=False)\n",
    "\n",
    "        os.makedirs(os.path.dirname(cleaned_save_fn), exist_ok=True)\n",
    "        cleaned_data.to_filename('../derivatives/cleaned_func/sub-{}/ses-{}/func/{}'.format(gd['sub'], gd['ses'], os.path.basename(fn)))\n",
    "\n",
    "all_funcs = sorted(glob.glob('../derivatives/fmriprep/fmriprep/sub-*/ses*/func/*T1w*_bold.nii.gz'))\n",
    "# all_funcs = [x for x in all_funcs if not 'sub-001' in x]\n",
    "# all_funcs\n",
    "\n",
    "all_highpassed = sorted(glob.glob('../derivatives/high_passed_func/sub*/ses*/func/*'))\n",
    "all_highpassed = [x for x in all_highpassed if 'rlsat' in x]\n",
    "\n",
    "out = joblib.Parallel(n_jobs=20, verbose=True)(joblib.delayed(do_clean)(x, overwrite=True) for x in all_highpassed)\n",
    "\n",
    "# do_clean(all_highpassed[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# calculate number of voxels in the IFG"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(scott) this was necessary to satisfy review comments, comparing the number of voxels in the IFG when using the HCP_MMP1 atlas vs the Harvard Oxford atlas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_voxel_numbers(to_run, atlas_name='MASSP', overwrite=False, to_psc=False, use_hp=False):\n",
    "    sub, ses, task, run = to_run\n",
    "    sub = str(sub).zfill(3)\n",
    "    print(f'Extracting from sub-{sub}/ses-{ses}/sub-{sub}_ses-{ses}_task-{task}_run-{run}', end='')\n",
    "    \n",
    "    epi_fn = get_epi(sub,ses,task,run,use_hp)\n",
    "    if not os.path.exists(epi_fn):\n",
    "        print('...doesnt exist, skipping'.format(sub,ses,task,run))\n",
    "        return None\n",
    "    \n",
    "    ## dont really need to convert to psc here\n",
    "    if to_psc:\n",
    "        epi = _make_psc(epi_fn)\n",
    "        psc_fn = '_psc'\n",
    "    else:\n",
    "        epi = nib.load(epi_fn)\n",
    "        psc_fn = ''\n",
    "    \n",
    "    # might wanna ahve the hp data handy\n",
    "    if use_hp:\n",
    "        output_fn = f'../derivatives/extracted_signals/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-{atlas_name}-signals{psc_fn}_hp.tsv'\n",
    "    else:\n",
    "        output_fn = f'../derivatives/extracted_signals/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-{atlas_name}-signals{psc_fn}.tsv'\n",
    "    \n",
    "    if os.path.exists(output_fn) and not overwrite:\n",
    "        print(f'{output_fn} already run, loading previous result...')\n",
    "        return pd.read_csv(output_fn, sep='\\t')\n",
    "    \n",
    "    #load & reshpae\n",
    "    epi_flat = epi.get_fdata().reshape((np.product(epi.shape[:3]), epi.shape[-1]))\n",
    "\n",
    "    # load atlas\n",
    "    atlas = load_atlas(sub,atlas_name=atlas_name)\n",
    "    \n",
    "    dfs = []\n",
    "    for i in np.arange(len(atlas.labels)):\n",
    "        print('.', end='')\n",
    "        label = list(atlas.labels)[i]\n",
    "        mask = image.index_img(atlas.maps, i)\n",
    "        mask_flat = mask.get_fdata().ravel()\n",
    "        if label in ['IFG-l', 'IFG-r','IFGhcp-l','IFGhcp-r']:\n",
    "            print(f'There are {np.count_nonzero(mask_flat)} voxels in region {label}')\n",
    "            return np.count_nonzero(mask_flat)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_atlases=['HCP_MMP1']\n",
    "# all_atlases = ['CORT']\n",
    "\n",
    "hp_options= [True]\n",
    "overwrite=True\n",
    "psc=False\n",
    "IFG_vox = []\n",
    "\n",
    "for atlas_name in all_atlases:\n",
    "    for hp in hp_options:\n",
    "        for i, comb in enumerate(all_combs):\n",
    "            print(f'atlas-{atlas_name} hp-{hp}')\n",
    "            print(comb)\n",
    "            sub = comb[0]\n",
    "            if check_affines(sub):\n",
    "                IFG_vox.append(calculate_voxel_numbers(comb, atlas_name=atlas_name, overwrite=overwrite, to_psc=psc, use_hp=hp))\n",
    "            else:\n",
    "                print(f'Affines for sub {sub} not identical')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# harvard oxford\n",
    "Harv_Ox_IFG = [1656,1544,1551,1881,1534,1622,1993,1522,2085,1484,1817,1772,2320,\n",
    " 2221,1910,1496,1549,1867,1968,1997,1485,1923,1780,1808,1626,2113,\n",
    " 1998,1825,1753,2140,1786,1925,1764,1759,1510,1836,1520]\n",
    "\n",
    "np.array(Harv_Ox_IFG).mean()\n",
    "# 1792.972972972973"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HCP_MMP1_IFG = [2311,2007,2391,2401,2344,1791,2474,1929,2576,1919,2752,2240,2777,2626,\n",
    " 2819,1849,2363,2389,2677,2836,2140,2532,2067,2325,2421,2618,2874,2648,\n",
    " 2578,2893,2536,2457,2626,2068,2080,2318,2247]\n",
    "\n",
    "np.array(HCP_MMP1_IFG).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TMP TO RUN NORMAL EXTRACTION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import os\n",
    "import re\n",
    "import shutil\n",
    "import random\n",
    "import string\n",
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import nilearn\n",
    "from nilearn import plotting, image\n",
    "from nilearn.input_data import NiftiMasker\n",
    "import nibabel as nib\n",
    "from nipype.interfaces import ants\n",
    "import nighres\n",
    "\n",
    "import subprocess\n",
    "import json\n",
    "import multiprocessing as mp\n",
    "from functools import partial\n",
    "import joblib\n",
    "from joblib import Parallel, delayed\n",
    "import itertools\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_rois(sub, atlas_name='ATAG', space='T1w'):\n",
    "    if atlas_name == 'ATAG':\n",
    "        if space == 'MNI152NLin2009cAsym' or space == 'mni':\n",
    "            ### Rois in MNI09c-space\n",
    "            mask_dir='/home/Public/trondheim/sourcedata/masks/MNI152NLin2009cAsym_res-1p5'\n",
    "            fns = sorted(glob.glob(mask_dir + '/space-*'))\n",
    "            names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_res-1p5_label-(?P<label>[a-zA-Z0-9]+)_probseg_def-img.nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "        else:\n",
    "            mask_dir = f'../derivatives/masks_atag_func/sub-{sub}/anat/sub-{sub}_*.nii.gz'\n",
    "            fns = sorted(glob.glob(mask_dir))\n",
    "            names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_desc-mask-(?P<label>[a-zA-Z0-9]+).nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "    elif atlas_name == 'MASSP':\n",
    "        mask_dir = f'../derivatives/masks_massp_func/sub-{sub}/anat/sub-{sub}_*.nii.gz'\n",
    "        fns = sorted(glob.glob(mask_dir))\n",
    "        names = [re.match('.*space-(?P<space>[a-zA-Z0-9]+)_desc-mask-(?P<label>\\S+).nii.gz', fn).groupdict()['label'] for fn in fns]\n",
    "\n",
    "    roi_dict = dict(zip(names, fns))\n",
    "    return roi_dict\n",
    "\n",
    "def load_atlas(sub, atlas_name='MASSP', space='T1w'):\n",
    "    from nilearn import image\n",
    "    \n",
    "    roi_dict = find_rois(sub, atlas_name, space)\n",
    "    if len(roi_dict) == 0:\n",
    "        warnings.warn(f'No ROIs found for sub-{sub} atlas-{atlas_name} space-{space}. Returning 0 to prevent error')\n",
    "        return 0\n",
    "    combined = image.concat_imgs(roi_dict.values())\n",
    "    \n",
    "    class AttrDict(dict):\n",
    "        def __init__(self, *args, **kwargs):\n",
    "            super(AttrDict, self).__init__(*args, **kwargs)\n",
    "            self.__dict__ = self\n",
    "            \n",
    "    roi_atlas = AttrDict({'maps': combined,\n",
    "                          'labels': roi_dict.keys()})\n",
    "    \n",
    "    return roi_atlas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_epi(sub, ses, task, run, use_hp=False, base_dir='../derivatives/fmriprep/fmriprep'):\n",
    "    if use_hp:\n",
    "        epi = os.path.join('../derivatives/high_passed_func', f'sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_space-T1w_desc-preproc_bold.nii.gz')\n",
    "    else:\n",
    "        epi = os.path.join(base_dir, f'sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_space-T1w_desc-preproc_bold.nii.gz')\n",
    "    return epi\n",
    "\n",
    "def _make_psc(data):\n",
    "    mean_img = image.mean_img(data)\n",
    "\n",
    "    # Replace 0s for numerical reasons\n",
    "    mean_data = mean_img.get_fdata()\n",
    "    mean_data[mean_data == 0] = 1\n",
    "    denom = image.new_img_like(mean_img, mean_data)\n",
    "\n",
    "    return image.math_img('data / denom[..., np.newaxis] * 100 - 100',\n",
    "                          data=data, denom=denom)\n",
    "\n",
    "def do_extract(to_run, atlas_name, overwrite=False, to_psc=False, use_hp=False):\n",
    "    sub, ses, task, run = to_run\n",
    "    sub = str(sub).zfill(3)\n",
    "    print(f'Extracting from sub-{sub}/ses-{ses}/sub-{sub}_ses-{ses}_task-{task}_run-{run}', end='')\n",
    "    \n",
    "    epi_fn = get_epi(sub,ses,task,run,use_hp)\n",
    "    if not os.path.exists(epi_fn):\n",
    "        print('...doesnt exist, skipping'.format(sub,ses,task,run))\n",
    "        return None\n",
    "    \n",
    "    # load atlas\n",
    "    atlas = load_atlas(sub, atlas_name=atlas_name)\n",
    "    if atlas == 0:\n",
    "        warnings.warn('No atlas found! skipping')\n",
    "        return None\n",
    "\n",
    "    if to_psc:\n",
    "        epi = _make_psc(epi_fn)\n",
    "        psc_fn = '_psc'\n",
    "    else:\n",
    "        epi = nib.load(epi_fn)\n",
    "        psc_fn = ''\n",
    "    \n",
    "    if use_hp:\n",
    "        output_fn = f'../derivatives/extracted_signals/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-{atasl_name}-signals{psc_fn}_hp.tsv'\n",
    "    else:\n",
    "        output_fn = f'../derivatives/extracted_signals/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-{atlas_name}-signals{psc_fn}.tsv'\n",
    "    print(output_fn)\n",
    "    if os.path.exists(output_fn) and not overwrite:\n",
    "        print(f'{output_fn} already run, loading previous result...')\n",
    "        return pd.read_csv(output_fn, sep='\\t')\n",
    "    \n",
    "    # load & reshape\n",
    "    epi_flat = epi.get_fdata().reshape((np.product(epi.shape[:3]), epi.shape[-1]))\n",
    "\n",
    "    dfs = []\n",
    "    for i in np.arange(len(atlas.labels)):\n",
    "        print('.', end='')\n",
    "        label = list(atlas.labels)[i]\n",
    "        mask = image.index_img(atlas.maps, i)\n",
    "        mask_flat = mask.get_fdata().ravel()\n",
    "        signal = pd.DataFrame(np.average(epi_flat, weights=mask_flat, axis=0), columns=[label])\n",
    "        signal.index.name = 'volume'\n",
    "        dfs.append(signal)\n",
    "\n",
    "    df = pd.concat(dfs, axis=1)\n",
    "#     output_fn = f'../derivatives/extracted_signals/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-MASSP-signals.tsv'\n",
    "    if not os.path.exists(os.path.dirname(output_fn)):\n",
    "        os.makedirs(os.path.dirname(output_fn))\n",
    "    df.to_csv(output_fn, sep='\\t')\n",
    "    print(output_fn)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find all available functional runs, extract sub/ses/task/run info\n",
    "all_runs = sorted(glob.glob('../derivatives/fmriprep/fmriprep/sub-*/ses-*/func/*space-T1w*_bold.nii.gz'))\n",
    "regex = re.compile('.*sub-(?P<sub>\\d+)_ses-(?P<ses>\\S+)_task-(?P<task>\\S+)_run-(?P<run>\\d)_space-T1w*')\n",
    "all_combs = [tuple(regex.match(x).groupdict().values()) for x in all_runs]\n",
    "# all_combs = [x for x in all_combs if x[0] in ['002','003','004','005','006','007','008','009','010','011']]\n",
    "# all_combs = [x for x in all_combs if x[0] in ['012','013','014','015','016','017','018','019','020','021','022','023','024','025','026']]\n",
    "#all_combs = [x for x in all_combs if x[0] in ['027','029','030','031','032']]\n",
    "\n",
    "# all_combs = [x for x in all_combs if x[1] == 'rlsat']\n",
    "all_combs = [x for x in all_combs if x[1] == 'sstmsit']\n",
    "# check if \n",
    "all_combs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find all available functional runs, extract sub/ses/task/run info\n",
    "all_runs = sorted(glob.glob('../derivatives/fmriprep/fmriprep/sub-*/ses-*/func/*space-T1w*_bold.nii.gz'))\n",
    "regex = re.compile('.*sub-(?P<sub>\\d+)_ses-(?P<ses>\\S+)_task-(?P<task>\\S+)_run-(?P<run>\\d)_space-T1w*')\n",
    "all_combs = [tuple(regex.match(x).groupdict().values()) for x in all_runs]\n",
    "# all_combs = [x for x in all_combs if x[0] in ['002','003','004','005','006','007','008','009','010','011']]\n",
    "# all_combs = [x for x in all_combs if x[0] in ['012','013','014','015','016','017','018','019','020','021','022','023','024','025','026']]\n",
    "#all_combs = [x for x in all_combs if x[0] in ['027','029','030','031','032']]\n",
    "\n",
    "# all_combs = [x for x in all_combs if x[1] == 'rlsat']\n",
    "# all_combs = [x for x in all_combs if x[1] == 'sstmsit']\n",
    "# check if \n",
    "all_combs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_affines(sub):\n",
    "    sub = str(sub).zfill(3)\n",
    "    all_funcs = sorted(glob.glob(f'../derivatives/fmriprep/fmriprep/sub-{sub}/ses*/func/sub*_space-T1w_desc-preproc_bold.nii.gz'))\n",
    "    all_affines = [nib.load(x).affine for x in all_funcs]\n",
    "    return (np.array(all_affines)[0] == np.array(all_affines)).all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, comb in enumerate(all_combs):\n",
    "    print(comb)\n",
    "    sub = comb[0]\n",
    "    if check_affines(sub):\n",
    "        do_extract(comb, atlas_name='ATAG',overwrite=False, to_psc=False, use_hp=False)\n",
    "    else:\n",
    "        print(f'Affines for sub {sub} not identical')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py38",
   "language": "python",
   "name": "py38"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}