MEG_preprocessing/test_extract_TS.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "A module that was compiled using NumPy 1.x cannot be run in\n",
      "NumPy 2.0.0 as it may crash. To support both 1.x and 2.x\n",
      "versions of NumPy, modules must be compiled with NumPy 2.0.\n",
      "Some module may need to rebuild instead e.g. with 'pybind11>=2.12'.\n",
      "\n",
      "If you are a user of the module, the easiest solution will be to\n",
      "downgrade to 'numpy<2' or try to upgrade the affected module.\n",
      "We expect that some modules will need time to support NumPy 2.\n",
      "\n",
      "Traceback (most recent call last):  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/runpy.py\", line 197, in _run_module_as_main\n",
      "    return _run_code(code, main_globals, None,\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/runpy.py\", line 87, in _run_code\n",
      "    exec(code, run_globals)\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel_launcher.py\", line 18, in <module>\n",
      "    app.launch_new_instance()\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/traitlets/config/application.py\", line 1075, in launch_instance\n",
      "    app.start()\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/kernelapp.py\", line 739, in start\n",
      "    self.io_loop.start()\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/tornado/platform/asyncio.py\", line 205, in start\n",
      "    self.asyncio_loop.run_forever()\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/asyncio/base_events.py\", line 596, in run_forever\n",
      "    self._run_once()\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/asyncio/base_events.py\", line 1890, in _run_once\n",
      "    handle._run()\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/asyncio/events.py\", line 80, in _run\n",
      "    self._context.run(self._callback, *self._args)\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/kernelbase.py\", line 545, in dispatch_queue\n",
      "    await self.process_one()\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/kernelbase.py\", line 534, in process_one\n",
      "    await dispatch(*args)\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/kernelbase.py\", line 437, in dispatch_shell\n",
      "    await result\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/ipkernel.py\", line 362, in execute_request\n",
      "    await super().execute_request(stream, ident, parent)\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/kernelbase.py\", line 778, in execute_request\n",
      "    reply_content = await reply_content\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/ipkernel.py\", line 449, in do_execute\n",
      "    res = shell.run_cell(\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/zmqshell.py\", line 549, in run_cell\n",
      "    return super().run_cell(*args, **kwargs)\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/IPython/core/interactiveshell.py\", line 3048, in run_cell\n",
      "    result = self._run_cell(\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/IPython/core/interactiveshell.py\", line 3103, in _run_cell\n",
      "    result = runner(coro)\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/IPython/core/async_helpers.py\", line 129, in _pseudo_sync_runner\n",
      "    coro.send(None)\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/IPython/core/interactiveshell.py\", line 3308, in run_cell_async\n",
      "    has_raised = await self.run_ast_nodes(code_ast.body, cell_name,\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/IPython/core/interactiveshell.py\", line 3490, in run_ast_nodes\n",
      "    if await self.run_code(code, result, async_=asy):\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/IPython/core/interactiveshell.py\", line 3550, in run_code\n",
      "    exec(code_obj, self.user_global_ns, self.user_ns)\n",
      "  File \"/tmp/ipykernel_10162/690586036.py\", line 4, in <module>\n",
      "    import pandas as pd\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/pandas/__init__.py\", line 26, in <module>\n",
      "    from pandas.compat import (\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/pandas/compat/__init__.py\", line 27, in <module>\n",
      "    from pandas.compat.pyarrow import (\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/pandas/compat/pyarrow.py\", line 8, in <module>\n",
      "    import pyarrow as pa\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/pyarrow/__init__.py\", line 65, in <module>\n",
      "    import pyarrow.lib as _lib\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "_ARRAY_API not found",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;31mAttributeError\u001b[0m: _ARRAY_API not found"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "A module that was compiled using NumPy 1.x cannot be run in\n",
      "NumPy 2.0.0 as it may crash. To support both 1.x and 2.x\n",
      "versions of NumPy, modules must be compiled with NumPy 2.0.\n",
      "Some module may need to rebuild instead e.g. with 'pybind11>=2.12'.\n",
      "\n",
      "If you are a user of the module, the easiest solution will be to\n",
      "downgrade to 'numpy<2' or try to upgrade the affected module.\n",
      "We expect that some modules will need time to support NumPy 2.\n",
      "\n",
      "Traceback (most recent call last):  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/runpy.py\", line 197, in _run_module_as_main\n",
      "    return _run_code(code, main_globals, None,\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/runpy.py\", line 87, in _run_code\n",
      "    exec(code, run_globals)\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel_launcher.py\", line 18, in <module>\n",
      "    app.launch_new_instance()\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/traitlets/config/application.py\", line 1075, in launch_instance\n",
      "    app.start()\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/kernelapp.py\", line 739, in start\n",
      "    self.io_loop.start()\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/tornado/platform/asyncio.py\", line 205, in start\n",
      "    self.asyncio_loop.run_forever()\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/asyncio/base_events.py\", line 596, in run_forever\n",
      "    self._run_once()\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/asyncio/base_events.py\", line 1890, in _run_once\n",
      "    handle._run()\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/asyncio/events.py\", line 80, in _run\n",
      "    self._context.run(self._callback, *self._args)\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/kernelbase.py\", line 545, in dispatch_queue\n",
      "    await self.process_one()\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/kernelbase.py\", line 534, in process_one\n",
      "    await dispatch(*args)\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/kernelbase.py\", line 437, in dispatch_shell\n",
      "    await result\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/ipkernel.py\", line 362, in execute_request\n",
      "    await super().execute_request(stream, ident, parent)\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/kernelbase.py\", line 778, in execute_request\n",
      "    reply_content = await reply_content\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/ipkernel.py\", line 449, in do_execute\n",
      "    res = shell.run_cell(\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/ipykernel/zmqshell.py\", line 549, in run_cell\n",
      "    return super().run_cell(*args, **kwargs)\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/IPython/core/interactiveshell.py\", line 3048, in run_cell\n",
      "    result = self._run_cell(\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/IPython/core/interactiveshell.py\", line 3103, in _run_cell\n",
      "    result = runner(coro)\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/IPython/core/async_helpers.py\", line 129, in _pseudo_sync_runner\n",
      "    coro.send(None)\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/IPython/core/interactiveshell.py\", line 3308, in run_cell_async\n",
      "    has_raised = await self.run_ast_nodes(code_ast.body, cell_name,\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/IPython/core/interactiveshell.py\", line 3490, in run_ast_nodes\n",
      "    if await self.run_code(code, result, async_=asy):\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/IPython/core/interactiveshell.py\", line 3550, in run_code\n",
      "    exec(code_obj, self.user_global_ns, self.user_ns)\n",
      "  File \"/tmp/ipykernel_10162/690586036.py\", line 4, in <module>\n",
      "    import pandas as pd\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/pandas/__init__.py\", line 49, in <module>\n",
      "    from pandas.core.api import (\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/pandas/core/api.py\", line 9, in <module>\n",
      "    from pandas.core.dtypes.dtypes import (\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/pandas/core/dtypes/dtypes.py\", line 24, in <module>\n",
      "    from pandas._libs import (\n",
      "  File \"/headnode1/abry4213/.conda/envs/pyspi/lib/python3.9/site-packages/pyarrow/__init__.py\", line 65, in <module>\n",
      "    import pyarrow.lib as _lib\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "_ARRAY_API not found",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;31mAttributeError\u001b[0m: _ARRAY_API not found"
     ]
    }
   ],
   "source": [
    "import os.path as op\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from fpdf import FPDF\n",
    "import shutil\n",
    "\n",
    "import mne\n",
    "import matplotlib.pyplot as plt\n",
    "import mne_bids\n",
    "\n",
    "import sys\n",
    "\n",
    "import argparse\n",
    "import itertools\n",
    "import json\n",
    "import pandas as pd\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "subject_id=\"CB003\"\n",
    "visit_id=\"1\"\n",
    "bids_root=\"/headnode1/abry4213/data/Cogitate_Batch2/MEG_Data/\"\n",
    "\n",
    "# Set path to preprocessing derivatives\n",
    "prep_deriv_root = op.join(bids_root, \"derivatives\", \"preprocessing\")\n",
    "prep_figure_root =  op.join(prep_deriv_root,\n",
    "                            f\"sub-{subject_id}\",f\"ses-{visit_id}\",\"meg\",\n",
    "                            \"figures\")\n",
    "prep_report_root =  op.join(prep_deriv_root,\n",
    "                            f\"sub-{subject_id}\",f\"ses-{visit_id}\",\"meg\",\n",
    "                            \"reports\")\n",
    "prep_code_root = op.join(prep_deriv_root,\n",
    "                            f\"sub-{subject_id}\",f\"ses-{visit_id}\",\"meg\",\n",
    "                            \"codes\")\n",
    "\n",
    "# Lopp over runs\n",
    "data_path = os.path.join(bids_root,f\"sub-{subject_id}\",f\"ses-{visit_id}\",\"meg\")\n",
    "run = \"05\"\n",
    "bids_task = 'dur'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_option=\"hypothesis_driven\"\n",
    "factor = ['Category', 'Task_relevance', \"Duration\"]\n",
    "conditions = [['face', 'object', 'letter', 'false'],\n",
    "              ['Relevant non-target', 'Irrelevant'],\n",
    "              ['1000ms']]\n",
    "\n",
    "fs_deriv_root = op.join(bids_root, \"derivatives\", \"fs\")\n",
    "rois_deriv_root = op.join(bids_root, \"derivatives\", \"roilabel\")\n",
    "prep_deriv_root = op.join(bids_root, \"derivatives\", \"preprocessing\")\n",
    "\n",
    "time_series_output_path = op.join(bids_root, \"derivatives\", \"MEG_time_series\")\n",
    "if not op.exists(time_series_output_path):\n",
    "    os.makedirs(time_series_output_path, exist_ok=True)\n",
    "\n",
    "# Time series output path for this subject\n",
    "subject_time_series_output_path = op.join(time_series_output_path, f\"sub-{subject_id}\", f\"ses-{visit_id}\", \"meg\")\n",
    "if not op.exists(subject_time_series_output_path):\n",
    "    os.makedirs(subject_time_series_output_path, exist_ok=True)\n",
    "    \n",
    "# Set task\n",
    "bids_task = 'dur'\n",
    "\n",
    "# Read epoched data\n",
    "bids_path_epo = mne_bids.BIDSPath(\n",
    "    root=prep_deriv_root, \n",
    "    subject=subject_id,  \n",
    "    datatype='meg',  \n",
    "    task=bids_task,\n",
    "    session=visit_id, \n",
    "    suffix='epo',\n",
    "    extension='.fif',\n",
    "    check=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Helper function to create a dictionary of ROI labels depending on the type of region subset requested\n",
    "def compute_ROI_labels(labels_atlas, region_option, rois_deriv_root):\n",
    "    # Create dictionary to store labels and vertices\n",
    "    labels_dict = {}\n",
    "    if region_option == 'hypothesis_driven':\n",
    "        # Read GNW and IIT ROI list\n",
    "        f = open(op.join(rois_deriv_root,\n",
    "                        'hypothesis_driven_ROIs.json'))\n",
    "        hypothesis_driven_ROIs = json.load(f)\n",
    "\n",
    "        # GNWT ROIs\n",
    "        print(\"GNWT ROIs:\")\n",
    "        for lab in hypothesis_driven_ROIs['GNWT_ROIs']:\n",
    "            print(lab)\n",
    "            labels_dict[\"GNWT_\"+lab] = np.sum([l for l in labels_atlas if lab in l.name])\n",
    "\n",
    "        # IIT ROIs\n",
    "        print(\"IIT ROIs\")\n",
    "        for lab in hypothesis_driven_ROIs['IIT_ROIs']:\n",
    "            print(lab)\n",
    "            labels_dict[\"IIT_\"+lab] = np.sum([l for l in labels_atlas if lab in l.name])\n",
    "\n",
    "        # Category-selective ROIs\n",
    "        print(\"Category-selective ROIs:\")\n",
    "        for lab in hypothesis_driven_ROIs['Category_Selective_ROIs']:\n",
    "            print(lab)\n",
    "            labels_dict[\"Category_Selective_\"+lab] = np.sum([l for l in labels_atlas if lab in l.name])\n",
    "\n",
    "        # Merge all labels in a single one separatelly for GNW and IIT \n",
    "        labels_dict['GNWT_meta_ROI'] = np.sum([l for l_name, l in labels_dict.items() if 'GNWT' in l_name])\n",
    "        labels_dict['IIT_meta_ROI'] = np.sum([l for l_name, l in labels_dict.items() if 'IIT' in l_name])\n",
    "        labels_dict['Category_Selective_meta_ROI'] = np.sum([l for l_name, l in labels_dict.items() if 'Category_Selective' in l_name])\n",
    "\n",
    "        # Only keep the meta-ROIs\n",
    "        labels_dict = {k: v for k, v in labels_dict.items() if 'meta_ROI' in k}\n",
    "\n",
    "    else:\n",
    "        for label in labels_atlas: \n",
    "            label_name = label.name\n",
    "            labels_dict[label_name] =  np.sum([label])\n",
    "\n",
    "    return labels_dict\n",
    "\n",
    "# Helper function to compute covariance matrices and inverse solution \n",
    "def fit_cov_and_inverse(subject_id, visit_id, factor, conditions, bids_root, downsample=True, tmin=-0.5, tmax=1.99):\n",
    "    # Set directory paths\n",
    "    prep_deriv_root = op.join(bids_root, \"derivatives\", \"preprocessing\")\n",
    "    fwd_deriv_root = op.join(bids_root, \"derivatives\", \"forward\")\n",
    "    source_deriv_root = op.join(bids_root, \"derivatives\", \"source_dur_ERF\")\n",
    "\n",
    "    if not op.exists(source_deriv_root):\n",
    "        os.makedirs(source_deriv_root, exist_ok=True)\n",
    "\n",
    "    source_figure_root =  op.join(source_deriv_root,\n",
    "                                f\"sub-{subject_id}\",f\"ses-{visit_id}\",\"meg\",\n",
    "                                \"figures\")\n",
    "    if not op.exists(source_figure_root):\n",
    "        os.makedirs(source_figure_root)\n",
    "\n",
    "    # Set task\n",
    "    bids_task = 'dur'\n",
    "    \n",
    "    # Read epoched data\n",
    "    bids_path_epo = mne_bids.BIDSPath(\n",
    "        root=prep_deriv_root, \n",
    "        subject=subject_id,  \n",
    "        datatype='meg',  \n",
    "        task=bids_task,\n",
    "        session=visit_id, \n",
    "        suffix='epo',\n",
    "        extension='.fif',\n",
    "        check=False)\n",
    "    \n",
    "    bids_path_epo_rs = mne_bids.BIDSPath(\n",
    "        root=prep_deriv_root, \n",
    "        subject=subject_id,  \n",
    "        datatype='meg',  \n",
    "        task=bids_task,\n",
    "        session=visit_id, \n",
    "        suffix='epo_rs',\n",
    "        extension='.fif',\n",
    "        check=False)\n",
    "    \n",
    "    print(\"Loading epochs data\")\n",
    "    epochs_all = mne.read_epochs(bids_path_epo.fpath, preload=True)\n",
    "\n",
    "    if not downsample:\n",
    "        epochs_final = epochs_all\n",
    "\n",
    "    # If downsampling is requested\n",
    "    else:\n",
    "        print(\"Applying downsampling\")\n",
    "        if os.path.exists(bids_path_epo_rs.fpath):\n",
    "            epochs_rs = mne.read_epochs(bids_path_epo_rs.fpath,\n",
    "                                    preload=True)\n",
    "        else:\n",
    "            epochs_all = mne.read_epochs(bids_path_epo.fpath,\n",
    "                                    preload=True)\n",
    "            resample_epochs(epochs_all, sfreq, bids_path_epo_rs, tmin=tmin, tmax=tmax)\n",
    "            epochs_rs = mne.read_epochs(bids_path_epo_rs.fpath, preload=True)\n",
    "        epochs_final = epochs_rs\n",
    "\n",
    "    # Run baseline correction\n",
    "    print(\"Running baseline correction\")\n",
    "    b_tmin = tmin\n",
    "    b_tmax = 0.\n",
    "    baseline = (b_tmin, b_tmax)\n",
    "    epochs_final.apply_baseline(baseline=baseline)\n",
    "\n",
    "    # Compute rank\n",
    "    print(\"Computing the rank\")\n",
    "    if os.path.isfile(f\"{fwd_deriv_root}/sub-{subject_id}_ses-{visit_id}_task-{bids_task}_rank.pkl\"):\n",
    "        with open(f\"{fwd_deriv_root}/sub-{subject_id}_ses-{visit_id}_task-{bids_task}_rank.pkl\", 'rb') as f:\n",
    "            rank = pickle.load(f)\n",
    "    else: \n",
    "        rank = mne.compute_rank(epochs_final, \n",
    "                                tol=1e-6, \n",
    "                                tol_kind='relative')\n",
    "        with open(f\"{fwd_deriv_root}/sub-{subject_id}_ses-{visit_id}_task-{bids_task}_rank.pkl\", 'wb') as f:\n",
    "            pickle.dump(rank, f)\n",
    "\n",
    "    # Read forward model\n",
    "    print(\"Reading forward model\")\n",
    "    bids_path_fwd = bids_path_epo.copy().update(\n",
    "            root=fwd_deriv_root,\n",
    "            task=bids_task,\n",
    "            suffix=\"surface_fwd\",\n",
    "            extension='.fif',\n",
    "            check=False)\n",
    "    fwd = mne.read_forward_solution(bids_path_fwd.fpath)\n",
    "\n",
    "    # Compute covariance matrices\n",
    "    print(\"Computing covariance matrices\")\n",
    "    if os.path.isfile(f\"{fwd_deriv_root}/sub-{subject_id}_ses-{visit_id}_task-{bids_task}_common_cov.pkl\"): \n",
    "        with open(f\"{fwd_deriv_root}/sub-{subject_id}_ses-{visit_id}_task-{bids_task}_common_cov.pkl\", 'rb') as f:\n",
    "            common_cov = pickle.load(f)\n",
    "    else:\n",
    "        base_cov = mne.compute_covariance(epochs_final, \n",
    "                                        tmin=-0.5, \n",
    "                                        tmax=0, \n",
    "                                        n_jobs=n_jobs,\n",
    "                                        method='empirical', \n",
    "                                        rank=rank)\n",
    "\n",
    "        active_cov = mne.compute_covariance(epochs_final, \n",
    "                                        tmin=0,\n",
    "                                        tmax=None,\n",
    "                                        n_jobs=n_jobs,\n",
    "                                        method='empirical', \n",
    "                                        rank=rank)\n",
    "        common_cov = base_cov + active_cov\n",
    "\n",
    "        with open(f\"{fwd_deriv_root}/sub-{subject_id}_ses-{visit_id}_task-{bids_task}_common_cov.pkl\", 'wb') as f:\n",
    "            pickle.dump(common_cov, f)\n",
    "\n",
    "    # Make inverse operator\n",
    "    inverse_operator = mne.minimum_norm.make_inverse_operator(\n",
    "        epochs_final.info,\n",
    "        fwd, \n",
    "        common_cov,\n",
    "        loose=.2,\n",
    "        depth=.8,\n",
    "        fixed=False,\n",
    "        rank=rank,\n",
    "        use_cps=True)\n",
    "\n",
    "    # Find all combinations between variables' levels\n",
    "    if len(factor) == 1:\n",
    "        cond_combs = list(itertools.product(conditions[0]))\n",
    "    if len(factor) == 2:\n",
    "        cond_combs = list(itertools.product(conditions[0],\n",
    "                                            conditions[1]))\n",
    "    if len(factor) == 3:\n",
    "        cond_combs = list(itertools.product(conditions[0],\n",
    "                                            conditions[1],\n",
    "                                            conditions[2]))\n",
    "        \n",
    "    return epochs_final, inverse_operator, cond_combs\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading labels from parcellation...\n",
      "   read 75 labels from /headnode1/abry4213/data/Cogitate_Batch2/MEG_Data/derivatives/fs/sub-CB003/label/lh.aparc.a2009s.annot\n",
      "   read 75 labels from /headnode1/abry4213/data/Cogitate_Batch2/MEG_Data/derivatives/fs/sub-CB003/label/rh.aparc.a2009s.annot\n",
      "GNWT ROIs:\n",
      "G_and_S_cingul-Ant\n",
      "G_and_S_cingul-Mid-Ant\n",
      "G_and_S_cingul-Mid-Post\n",
      "G_front_inf-Opercular\n",
      "G_front_inf-Orbital\n",
      "G_front_inf-Triangul\n",
      "G_front_middle\n",
      "Lat_Fis-ant-Horizont\n",
      "Lat_Fis-ant-Vertical\n",
      "S_front_inf\n",
      "S_front_middle\n",
      "S_front_sup\n",
      "IIT ROIs\n",
      "G_cuneus\n",
      "G_occipital_sup\n",
      "G_oc-temp_med-Lingual\n",
      "G_oc-temp_med-Parahip\n",
      "G_temporal_inf\n",
      "Pole_occipital\n",
      "Pole_temporal\n",
      "S_calcarine\n",
      "S_intrapariet_and_P_trans\n",
      "S_oc_sup_and_transversal\n",
      "S_temporal_sup\n",
      "Category-selective ROIs:\n",
      "G_and_S_occipital_inf\n",
      "G_oc-temp_lat-fusifor\n",
      "G_occipital_middle\n",
      "S_oc_middle_and_Lunatus\n",
      "Now finding inverse operator\n",
      "Loading epochs data\n",
      "Reading /headnode1/abry4213/data/Cogitate_Batch2/MEG_Data/derivatives/preprocessing/sub-CB003/ses-1/meg/sub-CB003_ses-1_task-dur_epo.fif ...\n",
      "    Found the data of interest:\n",
      "        t =    -500.00 ...    2000.00 ms\n",
      "        0 CTF compensation matrices available\n",
      "Reading /headnode1/abry4213/data/Cogitate_Batch2/MEG_Data/derivatives/preprocessing/sub-CB003/ses-1/meg/sub-CB003_ses-1_task-dur_epo-1.fif ...\n",
      "    Found the data of interest:\n",
      "        t =    -500.00 ...    2000.00 ms\n",
      "        0 CTF compensation matrices available\n",
      "Reading /headnode1/abry4213/data/Cogitate_Batch2/MEG_Data/derivatives/preprocessing/sub-CB003/ses-1/meg/sub-CB003_ses-1_task-dur_epo-2.fif ...\n",
      "    Found the data of interest:\n",
      "        t =    -500.00 ...    2000.00 ms\n",
      "        0 CTF compensation matrices available\n",
      "Reading /headnode1/abry4213/data/Cogitate_Batch2/MEG_Data/derivatives/preprocessing/sub-CB003/ses-1/meg/sub-CB003_ses-1_task-dur_epo-3.fif ...\n",
      "    Found the data of interest:\n",
      "        t =    -500.00 ...    2000.00 ms\n",
      "        0 CTF compensation matrices available\n"
     ]
    }
   ],
   "source": [
    "# Use Desikan--Killiany atlas to compute dictionary of labels\n",
    "labels_atlas = mne.read_labels_from_annot(\n",
    "    \"sub-\"+subject_id, \n",
    "    parc='aparc.a2009s',\n",
    "    subjects_dir=fs_deriv_root)\n",
    "labels_dict = compute_ROI_labels(labels_atlas, region_option, rois_deriv_root)\n",
    "\n",
    "# Find epochs_rs, inverse_operator, cond_combs\n",
    "print(\"Now finding inverse operator\")\n",
    "epochs_final, inverse_operator, cond_combs = fit_cov_and_inverse(subject_id, visit_id, factor, conditions, bids_root, downsample=False)\n",
    "\n",
    "# extract time course in label with pca_flip mode\n",
    "src = inverse_operator['src']\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}