aaron_stats.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### drift testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### one subject for testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###\n",
    "import sys\n",
    "print(sys.path)\n",
    "sys.path.append(\"C:/Users/jz421/Desktop/GlobalLocal/IEEG_Pipelines/\") #need to do this cuz otherwise ieeg isn't added to path...\n",
    "\n",
    "from ieeg.navigate import channel_outlier_marker, trial_ieeg, crop_empty_data, \\\n",
    "    outliers_to_nan\n",
    "from ieeg.io import raw_from_layout, get_data\n",
    "from ieeg.timefreq.utils import crop_pad\n",
    "from ieeg.timefreq import gamma\n",
    "from ieeg.calc.scaling import rescale\n",
    "import mne\n",
    "import os\n",
    "import numpy as np\n",
    "from ieeg.calc.reshape import make_data_same\n",
    "from ieeg.calc.stats import time_perm_cluster\n",
    "\n",
    "from misc_functions import calculate_RTs, save_channels_to_file, save_sig_chans\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def get_baseline(inst: mne.io.BaseRaw, base_times_list: list[tuple[float, float]]):\n",
    "    all_data = []\n",
    "    all_events = []\n",
    "    \n",
    "    inst = inst.copy()\n",
    "    inst.load_data()\n",
    "    ch_type = inst.get_channel_types(only_data_chs=True)[0]\n",
    "    inst.set_eeg_reference(ref_channels=\"average\", ch_type=ch_type)\n",
    "    \n",
    "    last_sample = 0\n",
    "    for base_times in base_times_list:\n",
    "        adjusted_base_times = [base_times[0] - 0.5, base_times[1] + 0.5]\n",
    "        trials = trial_ieeg(inst, \"experimentStart\", adjusted_base_times, preload=True)\n",
    "        \n",
    "        if trials is None or len(trials) == 0:\n",
    "            print(f\"Warning: No trials found for base_times {base_times}. Skipping...\")\n",
    "            continue\n",
    "        \n",
    "        all_data.append(trials.get_data())\n",
    "        all_events.extend([[last_sample + i, 0, 1] for i in range(len(trials))])\n",
    "        last_sample += len(trials)\n",
    "    \n",
    "    combined_data = np.concatenate(all_data, axis=0)\n",
    "    combined_trials = mne.EpochsArray(combined_data, inst.info, events=np.array(all_events), tmin=trials.tmin)\n",
    "    \n",
    "    HG_base = gamma.extract(combined_trials, copy=True, n_jobs=1)\n",
    "    broadband_base = gamma.extract(combined_trials, copy=True, n_jobs=1, passband=(0,1000))\n",
    "    crop_pad(HG_base, \"0.5s\")\n",
    "    crop_pad(broadband_base, \"0.5s\")\n",
    "    \n",
    "    del inst\n",
    "    return HG_base, broadband_base\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def plot_HG_and_stats(sub, task, output_name, event=None, times=(-1, 1.5),\n",
    "                      base_times=(-0.5, 0), LAB_root=None, channels=None,\n",
    "                      full_trial_base=False):\n",
    "    \"\"\"\n",
    "    Plot high gamma (HG) and statistics for a given subject and task using specified event.\n",
    "\n",
    "    Parameters:\n",
    "    - sub (str): The subject identifier.\n",
    "    - task (str): The task identifier.\n",
    "    - output_name (str): The name for the output files.\n",
    "    - event (str, optional): Event name to process. Defaults to None.\n",
    "    - times (tuple, optional): A tuple indicating the start and end times for processing. Defaults to (-1, 1.5).\n",
    "    - base_times (tuple, optional): A tuple indicating the start and end base times for processing. Defaults to (-0.5, 0).\n",
    "    - LAB_root (str, optional): The root directory for the lab. Will be determined based on OS if not provided. Defaults to None.\n",
    "    - channels (list of strings, optional): The channels to plot and get stats for. Default is all channels.\n",
    "    - full_trial_base (boolean): Whether to use the full trial as the baseline period. Default is False.\n",
    "    This function will process the provided event for a given subject and task.\n",
    "    High gamma (HG) will be computed, and statistics will be calculated and plotted.\n",
    "    The results will be saved to output files.\n",
    "    \"\"\"\n",
    "    pass\n",
    "\n",
    "\n",
    "sub = 'D0057'\n",
    "task = 'GlobalLocal'\n",
    "output_name = \"driftTest_Stimulus_experimentStartBaseNormalizedAnd0.5SecBeforeEventTrialsonXAxis\"\n",
    "events = [\"Stimulus\"]\n",
    "times = (-1,1.5)\n",
    "base_times = [(1,31),(41,71)]\n",
    "LAB_root = None\n",
    "channels = None\n",
    "full_trial_base = False\n",
    "\n",
    "if LAB_root is None:\n",
    "    HOME = os.path.expanduser(\"~\")\n",
    "    if os.name == 'nt':  # windows\n",
    "        LAB_root = os.path.join(HOME, \"Box\", \"CoganLab\")\n",
    "    else:  # mac\n",
    "        LAB_root = os.path.join(HOME, \"Library\", \"CloudStorage\", \"Box-Box\",\n",
    "                                \"CoganLab\")\n",
    "\n",
    "layout = get_data(task, root=LAB_root)\n",
    "filt = raw_from_layout(layout.derivatives['derivatives/clean'], subject=sub,\n",
    "                       extension='.edf', desc='clean', preload=False)\n",
    "save_dir = os.path.join(layout.root, 'derivatives', 'freqFilt', 'figs', sub)\n",
    "if not os.path.exists(save_dir):\n",
    "    os.makedirs(save_dir)\n",
    "\n",
    "###\n",
    "\n",
    "good = crop_empty_data(filt)\n",
    "# %%\n",
    "\n",
    "print(f\"good channels before dropping bads: {len(good.ch_names)}\")\n",
    "print(f\"filt channels before dropping bads: {len(filt.ch_names)}\")\n",
    "\n",
    "good.info['bads'] = channel_outlier_marker(good, 3, 2)\n",
    "print(\"Bad channels in 'good':\", good.info['bads'])\n",
    "\n",
    "filt.drop_channels(good.info['bads'])  # this has to come first cuz if you drop from good first, then good.info['bads'] is just empty\n",
    "good.drop_channels(good.info['bads'])\n",
    "\n",
    "print(\"Bad channels in 'good' after dropping once:\", good.info['bads'])\n",
    "\n",
    "print(f\"good channels after dropping bads: {len(good.ch_names)}\")\n",
    "print(f\"filt channels after dropping bads: {len(filt.ch_names)}\")\n",
    "\n",
    "HG_base, broadband_base = get_baseline(filt, base_times)\n",
    "good.load_data()\n",
    "\n",
    "# If channels is None, use all channels\n",
    "if channels is None:\n",
    "    channels = good.ch_names\n",
    "else:\n",
    "    # Validate the provided channels\n",
    "    invalid_channels = [ch for ch in channels if ch not in good.ch_names]\n",
    "    if invalid_channels:\n",
    "        raise ValueError(\n",
    "            f\"The following channels are not valid: {invalid_channels}\")\n",
    "\n",
    "    # Use only the specified channels\n",
    "    good.pick_channels(channels)\n",
    "\n",
    "ch_type = filt.get_channel_types(only_data_chs=True)[0]\n",
    "good.set_eeg_reference(ref_channels=\"average\", ch_type=ch_type)\n",
    "\n",
    "\n",
    "# Create a baseline EpochsTFR using the stimulus event\n",
    "\n",
    "# adjusted_base_times = [base_times[0] - 0.5, base_times[1] + 0.5]\n",
    "# trials = trial_ieeg(good, \"Stimulus\", adjusted_base_times, preload=True)\n",
    "# outliers_to_nan(trials, outliers=10)\n",
    "# HG_base_fixationCross = gamma.extract(trials, copy=False, n_jobs=1)\n",
    "# crop_pad(HG_base_fixationCross, \"0.5s\")\n",
    "\n",
    "\n",
    "all_epochs_list = []\n",
    "\n",
    "for event in events:\n",
    "# Epoching and HG extraction for each specified event. Then concatenate all trials epochs objects together (do Stimulus/c25 and Stimulus/c75 for example, and combine to get all congruent trials)\n",
    "    times_adj = [times[0] - 0.5, times[1] + 0.5]\n",
    "    trials = trial_ieeg(good, event, times_adj, preload=True,\n",
    "                        reject_by_annotation=False)\n",
    "    all_epochs_list.append(trials)\n",
    "\n",
    "# Concatenate all trials\n",
    "all_trials = mne.concatenate_epochs(all_epochs_list)\n",
    "\n",
    "outliers_to_nan(all_trials, outliers=10)\n",
    "HG_ev1 = gamma.extract(all_trials, copy=True, n_jobs=1)\n",
    "broadband_ev1 = gamma.extract(all_trials, copy=True, n_jobs=1, passband=(0,1000))\n",
    "print(\"HG_ev1 before crop_pad: \", HG_ev1.tmin, HG_ev1.tmax)\n",
    "crop_pad(HG_ev1, \"0.5s\")\n",
    "print(\"HG_ev1 after crop_pad: \", HG_ev1.tmin, HG_ev1.tmax)\n",
    "\n",
    "crop_pad(broadband_ev1, \"0.5s\")\n",
    "HG_ev1_rescaled = rescale(HG_ev1, HG_base, copy=True, mode='zscore')\n",
    "broadband_ev1_rescaled = rescale(broadband_ev1, broadband_base, copy=True, mode='zscore')\n",
    "# HG_ev1_rescaled_fixationCross = rescale(HG_ev1, HG_base_fixationCross, copy=True, mode='zscore')\n",
    "\n",
    "HG_base.decimate(2)\n",
    "HG_ev1.decimate(2)\n",
    "\n",
    "broadband_base.decimate(2)\n",
    "broadband_ev1.decimate(2)\n",
    "\n",
    "# HG_base_fixationCross.decimate(2)\n",
    "\n",
    "\n",
    "HG_ev1_avgOverTime = np.nanmean(HG_ev1.get_data(), axis=2)\n",
    "HG_ev1_rescaled_avgOverTime = np.nanmean(HG_ev1_rescaled.get_data(), axis=2)\n",
    "\n",
    "\n",
    "HG_ev1_evoke = HG_ev1.average(method=lambda x: np.nanmean(x, axis=0)) #axis=0 should be set for actually running this, the axis=2 is just for drift testing.\n",
    "HG_ev1_evoke_rescaled = HG_ev1_rescaled.average(method=lambda x: np.nanmean(x, axis=0))\n",
    "\n",
    "if event == \"Stimulus\":\n",
    "    print('plotting stimulus')\n",
    "    fig = HG_ev1_evoke_rescaled.plot(unit=False, scalings=dict(sEEG=1)) #this line is not finishing...\n",
    "    print('plotted')\n",
    "    # for ax in fig.axes:\n",
    "    #     ax.axvline(x=avg_RT, color='r', linestyle='--')\n",
    "    print('about to save')\n",
    "    fig.savefig(save_dir + '_HG_ev1_Stimulus_zscore.png')\n",
    "    print('saved')\n",
    "else:\n",
    "    print('about to plot if not stimulus')\n",
    "    fig = HG_ev1_evoke_rescaled.plot(unit=False, scalings=dict(sEEG=1))\n",
    "    print('plotted non stimulus')\n",
    "    fig.savefig(save_dir + f'_HG_ev1_{output_name}_zscore.png')\n",
    "\n",
    "###\n",
    "print(f\"Shape of HG_ev1._data: {HG_ev1._data.shape}\")\n",
    "print(f\"Shape of HG_base._data: {HG_base._data.shape}\")\n",
    "\n",
    "sig1 = HG_ev1._data\n",
    "sig2 = HG_base._data\n",
    "sig2 = make_data_same(sig2, sig1.shape)\n",
    "print(f\"Shape of sig1: {sig1.shape}\")\n",
    "print(f\"Shape of sig2: {sig2.shape}\")\n",
    "\n",
    "mat = time_perm_cluster(sig1, sig2, 0.05, n_jobs=6, ignore_adjacency=1)\n",
    "fig = plt.figure()\n",
    "plt.imshow(mat, aspect='auto')\n",
    "fig.savefig(save_dir + f'_{output_name}_stats.png', dpi=300)\n",
    "\n",
    "channels = good.ch_names\n",
    "\n",
    "#save channels with their indices \n",
    "save_channels_to_file(channels, sub, task, save_dir)\n",
    "\n",
    "# save significant channels to a json\n",
    "save_sig_chans(f'{output_name}', mat, channels, sub, save_dir)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filt[:,slice(0, 4000)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### loop over all subjects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# crop_pad(HG_ev1, \"0.5s\")\n",
    "HG_ev1.tmax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###\n",
    "import sys\n",
    "print(sys.path)\n",
    "sys.path.append(\"C:/Users/jz421/Desktop/GlobalLocal/IEEG_Pipelines/\") #need to do this cuz otherwise ieeg isn't added to path...\n",
    "\n",
    "from ieeg.navigate import channel_outlier_marker, trial_ieeg, crop_empty_data, \\\n",
    "    outliers_to_nan\n",
    "from ieeg.io import raw_from_layout, get_data\n",
    "from ieeg.timefreq.utils import crop_pad\n",
    "from ieeg.timefreq import gamma\n",
    "from ieeg.calc.scaling import rescale\n",
    "import mne\n",
    "import os\n",
    "import numpy as np\n",
    "from ieeg.calc.reshape import make_data_same\n",
    "from ieeg.calc.stats import time_perm_cluster\n",
    "\n",
    "from misc_functions import calculate_RTs, save_channels_to_file, save_sig_chans\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def get_baseline(inst: mne.io.BaseRaw, base_times: tuple[float, float]):\n",
    "    inst = inst.copy()\n",
    "    inst.load_data()\n",
    "    ch_type = inst.get_channel_types(only_data_chs=True)[0]\n",
    "    inst.set_eeg_reference(ref_channels=\"average\", ch_type=ch_type)\n",
    "\n",
    "    adjusted_base_times = [base_times[0] - 0.5, base_times[1] + 0.5]\n",
    "    trials = trial_ieeg(inst, \"experimentStart\", adjusted_base_times,\n",
    "                        preload=True)\n",
    "    # outliers_to_nan(trials, outliers=10)\n",
    "    HG_base = gamma.extract(trials, copy=False, n_jobs=1)\n",
    "    crop_pad(HG_base, \"0.5s\")\n",
    "    del inst\n",
    "    return HG_base\n",
    "\n",
    "\n",
    "def plot_HG_and_stats(sub, task, output_name, event=None, times=(-1, 1.5),\n",
    "                      base_times=(-0.5, 0), LAB_root=None, channels=None,\n",
    "                      full_trial_base=False):\n",
    "    \"\"\"\n",
    "    Plot high gamma (HG) and statistics for a given subject and task using specified event.\n",
    "\n",
    "    Parameters:\n",
    "    - sub (str): The subject identifier.\n",
    "    - task (str): The task identifier.\n",
    "    - output_name (str): The name for the output files.\n",
    "    - event (str, optional): Event name to process. Defaults to None.\n",
    "    - times (tuple, optional): A tuple indicating the start and end times for processing. Defaults to (-1, 1.5).\n",
    "    - base_times (tuple, optional): A tuple indicating the start and end base times for processing. Defaults to (-0.5, 0).\n",
    "    - LAB_root (str, optional): The root directory for the lab. Will be determined based on OS if not provided. Defaults to None.\n",
    "    - channels (list of strings, optional): The channels to plot and get stats for. Default is all channels.\n",
    "    - full_trial_base (boolean): Whether to use the full trial as the baseline period. Default is False.\n",
    "    This function will process the provided event for a given subject and task.\n",
    "    High gamma (HG) will be computed, and statistics will be calculated and plotted.\n",
    "    The results will be saved to output files.\n",
    "    \"\"\"\n",
    "    pass\n",
    "\n",
    "subjects = ['D0057', 'D0059', 'D0063', 'D0065', 'D0069', 'D0071']\n",
    "for sub in subjects:\n",
    "    sub = sub\n",
    "    output_name = \"responseIncongruent_experimentStartBase1secTo101Sec\"\n",
    "    events = [\"Response/i25\",\"Response/i75\"]\n",
    "    times = (-1, 1.5)\n",
    "    base_times = (1, 101)\n",
    "    LAB_root = None\n",
    "    channels = None\n",
    "    full_trial_base = False\n",
    "\n",
    "    if LAB_root is None:\n",
    "        HOME = os.path.expanduser(\"~\")\n",
    "        if os.name == 'nt':  # windows\n",
    "            LAB_root = os.path.join(HOME, \"Box\", \"CoganLab\")\n",
    "        else:  # mac\n",
    "            LAB_root = os.path.join(HOME, \"Library\", \"CloudStorage\", \"Box-Box\",\n",
    "                                    \"CoganLab\")\n",
    "\n",
    "    layout = get_data(task, root=LAB_root)\n",
    "    filt = raw_from_layout(layout.derivatives['derivatives/clean'], subject=sub,\n",
    "                        extension='.edf', desc='clean', preload=False)\n",
    "    save_dir = os.path.join(layout.root, 'derivatives', 'freqFilt', 'figs', sub)\n",
    "    if not os.path.exists(save_dir):\n",
    "        os.makedirs(save_dir)\n",
    "\n",
    "    ###\n",
    "\n",
    "    good = crop_empty_data(filt)\n",
    "    # %%\n",
    "\n",
    "    print(f\"good channels before dropping bads: {len(good.ch_names)}\")\n",
    "    print(f\"filt channels before dropping bads: {len(filt.ch_names)}\")\n",
    "\n",
    "    good.info['bads'] = channel_outlier_marker(good, 3, 2)\n",
    "    print(\"Bad channels in 'good':\", good.info['bads'])\n",
    "\n",
    "    filt.drop_channels(good.info['bads'])  # this has to come first cuz if you drop from good first, then good.info['bads'] is just empty\n",
    "    good.drop_channels(good.info['bads'])\n",
    "\n",
    "    print(\"Bad channels in 'good' after dropping once:\", good.info['bads'])\n",
    "\n",
    "    print(f\"good channels after dropping bads: {len(good.ch_names)}\")\n",
    "    print(f\"filt channels after dropping bads: {len(filt.ch_names)}\")\n",
    "\n",
    "    HG_base = get_baseline(filt, base_times)\n",
    "    good.load_data()\n",
    "\n",
    "    # If channels is None, use all channels\n",
    "    if channels is None:\n",
    "        channels = good.ch_names\n",
    "    else:\n",
    "        # Validate the provided channels\n",
    "        invalid_channels = [ch for ch in channels if ch not in good.ch_names]\n",
    "        if invalid_channels:\n",
    "            raise ValueError(\n",
    "                f\"The following channels are not valid: {invalid_channels}\")\n",
    "\n",
    "        # Use only the specified channels\n",
    "        good.pick_channels(channels)\n",
    "\n",
    "    ch_type = filt.get_channel_types(only_data_chs=True)[0]\n",
    "    good.set_eeg_reference(ref_channels=\"average\", ch_type=ch_type)\n",
    "\n",
    "    all_epochs_list = []\n",
    "\n",
    "    for event in events:\n",
    "    # Epoching and HG extraction for each specified event. Then concatenate all trials epochs objects together (do Stimulus/c25 and Stimulus/c75 for example, and combine to get all congruent trials)\n",
    "        times_adj = [times[0] - 0.5, times[1] + 0.5]\n",
    "        trials = trial_ieeg(good, event, times_adj, preload=True,\n",
    "                            reject_by_annotation=False)\n",
    "        all_epochs_list.append(trials)\n",
    "\n",
    "    # Concatenate all trials\n",
    "    all_trials = mne.concatenate_epochs(all_epochs_list)\n",
    "\n",
    "    outliers_to_nan(all_trials, outliers=10)\n",
    "    HG_ev1 = gamma.extract(all_trials, copy=False, n_jobs=1)\n",
    "    crop_pad(HG_ev1, \"0.5s\")\n",
    "    HG_ev1_rescaled = rescale(HG_ev1, HG_base, copy=True, mode='zscore')\n",
    "\n",
    "    HG_base.decimate(2)\n",
    "    HG_ev1.decimate(2)\n",
    "    RTs, skipped = calculate_RTs(good)\n",
    "    avg_RT = np.median(RTs)\n",
    "\n",
    "    HG_ev1_evoke = HG_ev1.average(method=lambda x: np.nanmean(x, axis=0))\n",
    "    HG_ev1_evoke_rescaled = HG_ev1_rescaled.average(method=lambda x: np.nanmean(x, axis=0))\n",
    "\n",
    "    if event == \"Stimulus\":\n",
    "        print('plotting stimulus')\n",
    "        fig = HG_ev1_evoke_rescaled.plot(unit=False, scalings=dict(sEEG=1)) #this line is not finishing...\n",
    "        print('plotted')\n",
    "        for ax in fig.axes:\n",
    "            ax.axvline(x=avg_RT, color='r', linestyle='--')\n",
    "        print('about to save')\n",
    "        fig.savefig(save_dir + '_HG_ev1_Stimulus_zscore.png')\n",
    "        print('saved')\n",
    "    else:\n",
    "        print('about to plot if not stimulus')\n",
    "        fig = HG_ev1_evoke_rescaled.plot(unit=False, scalings=dict(sEEG=1))\n",
    "        print('plotted non stimulus')\n",
    "        fig.savefig(save_dir + f'_HG_ev1_{output_name}_zscore.png')\n",
    "\n",
    "    ###\n",
    "    print(f\"Shape of HG_ev1._data: {HG_ev1._data.shape}\")\n",
    "    print(f\"Shape of HG_base._data: {HG_base._data.shape}\")\n",
    "\n",
    "    sig1 = HG_ev1._data\n",
    "    sig2 = HG_base._data\n",
    "    sig2 = make_data_same(sig2, sig1.shape)\n",
    "    print(f\"Shape of sig1: {sig1.shape}\")\n",
    "    print(f\"Shape of sig2: {sig2.shape}\")\n",
    "\n",
    "    mat = time_perm_cluster(sig1, sig2, 0.05, n_jobs=6, ignore_adjacency=1)\n",
    "    fig = plt.figure()\n",
    "    plt.imshow(mat, aspect='auto')\n",
    "    fig.savefig(save_dir + f'_{output_name}_stats.png', dpi=300)\n",
    "\n",
    "    channels = good.ch_names\n",
    "\n",
    "    #save channels with their indices \n",
    "    save_channels_to_file(channels, sub, task, save_dir)\n",
    "\n",
    "    # save significant channels to a json\n",
    "    save_sig_chans(f'{output_name}', mat, channels, sub, save_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###\n",
    "import sys\n",
    "print(sys.path)\n",
    "sys.path.append(\"C:/Users/jz421/Desktop/GlobalLocal/IEEG_Pipelines/\") #need to do this cuz otherwise ieeg isn't added to path...\n",
    "\n",
    "from ieeg.navigate import channel_outlier_marker, trial_ieeg, crop_empty_data, \\\n",
    "    outliers_to_nan\n",
    "from ieeg.io import raw_from_layout, get_data\n",
    "from ieeg.timefreq.utils import crop_pad\n",
    "from ieeg.timefreq import gamma\n",
    "from ieeg.calc.scaling import rescale\n",
    "import mne\n",
    "import os\n",
    "import numpy as np\n",
    "from ieeg.calc.reshape import make_data_same\n",
    "from ieeg.calc.stats import time_perm_cluster\n",
    "\n",
    "from misc_functions import calculate_RTs, save_channels_to_file, save_sig_chans\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def get_baseline(inst: mne.io.BaseRaw, base_times: tuple[float, float]):\n",
    "    inst = inst.copy()\n",
    "    inst.load_data()\n",
    "    ch_type = inst.get_channel_types(only_data_chs=True)[0]\n",
    "    inst.set_eeg_reference(ref_channels=\"average\", ch_type=ch_type)\n",
    "\n",
    "    adjusted_base_times = [base_times[0] - 0.5, base_times[1] + 0.5]\n",
    "    trials = trial_ieeg(inst, \"experimentStart\", adjusted_base_times,\n",
    "                        preload=True)\n",
    "    # outliers_to_nan(trials, outliers=10)\n",
    "    HG_base = gamma.extract(trials, copy=False, n_jobs=1)\n",
    "    crop_pad(HG_base, \"0.5s\")\n",
    "    del inst\n",
    "    return HG_base\n",
    "\n",
    "\n",
    "def plot_HG_and_stats(sub, task, output_name, event=None, times=(-1, 1.5),\n",
    "                      base_times=(-0.5, 0), LAB_root=None, channels=None,\n",
    "                      full_trial_base=False):\n",
    "    \"\"\"\n",
    "    Plot high gamma (HG) and statistics for a given subject and task using specified event.\n",
    "\n",
    "    Parameters:\n",
    "    - sub (str): The subject identifier.\n",
    "    - task (str): The task identifier.\n",
    "    - output_name (str): The name for the output files.\n",
    "    - event (str, optional): Event name to process. Defaults to None.\n",
    "    - times (tuple, optional): A tuple indicating the start and end times for processing. Defaults to (-1, 1.5).\n",
    "    - base_times (tuple, optional): A tuple indicating the start and end base times for processing. Defaults to (-0.5, 0).\n",
    "    - LAB_root (str, optional): The root directory for the lab. Will be determined based on OS if not provided. Defaults to None.\n",
    "    - channels (list of strings, optional): The channels to plot and get stats for. Default is all channels.\n",
    "    - full_trial_base (boolean): Whether to use the full trial as the baseline period. Default is False.\n",
    "    This function will process the provided event for a given subject and task.\n",
    "    High gamma (HG) will be computed, and statistics will be calculated and plotted.\n",
    "    The results will be saved to output files.\n",
    "    \"\"\"\n",
    "    pass\n",
    "\n",
    "subjects = ['D0057', 'D0059', 'D0063', 'D0065', 'D0069', 'D0071']\n",
    "for sub in subjects:\n",
    "    sub = sub\n",
    "    output_name = \"stimulus_experimentStartBase1secTo101Sec\"\n",
    "    events = [\"Stimulus\"]\n",
    "    times = (-1, 1.5)\n",
    "    base_times = (1, 101)\n",
    "    LAB_root = None\n",
    "    channels = None\n",
    "    full_trial_base = False\n",
    "\n",
    "    if LAB_root is None:\n",
    "        HOME = os.path.expanduser(\"~\")\n",
    "        if os.name == 'nt':  # windows\n",
    "            LAB_root = os.path.join(HOME, \"Box\", \"CoganLab\")\n",
    "        else:  # mac\n",
    "            LAB_root = os.path.join(HOME, \"Library\", \"CloudStorage\", \"Box-Box\",\n",
    "                                    \"CoganLab\")\n",
    "\n",
    "    layout = get_data(task, root=LAB_root)\n",
    "    filt = raw_from_layout(layout.derivatives['derivatives/clean'], subject=sub,\n",
    "                        extension='.edf', desc='clean', preload=False)\n",
    "    save_dir = os.path.join(layout.root, 'derivatives', 'freqFilt', 'figs', sub)\n",
    "    if not os.path.exists(save_dir):\n",
    "        os.makedirs(save_dir)\n",
    "\n",
    "    ###\n",
    "\n",
    "    good = crop_empty_data(filt)\n",
    "    # %%\n",
    "\n",
    "    print(f\"good channels before dropping bads: {len(good.ch_names)}\")\n",
    "    print(f\"filt channels before dropping bads: {len(filt.ch_names)}\")\n",
    "\n",
    "    good.info['bads'] = channel_outlier_marker(good, 3, 2)\n",
    "    print(\"Bad channels in 'good':\", good.info['bads'])\n",
    "\n",
    "    filt.drop_channels(good.info['bads'])  # this has to come first cuz if you drop from good first, then good.info['bads'] is just empty\n",
    "    good.drop_channels(good.info['bads'])\n",
    "\n",
    "    print(\"Bad channels in 'good' after dropping once:\", good.info['bads'])\n",
    "\n",
    "    print(f\"good channels after dropping bads: {len(good.ch_names)}\")\n",
    "    print(f\"filt channels after dropping bads: {len(filt.ch_names)}\")\n",
    "\n",
    "    HG_base = get_baseline(filt, base_times)\n",
    "    good.load_data()\n",
    "\n",
    "    # If channels is None, use all channels\n",
    "    if channels is None:\n",
    "        channels = good.ch_names\n",
    "    else:\n",
    "        # Validate the provided channels\n",
    "        invalid_channels = [ch for ch in channels if ch not in good.ch_names]\n",
    "        if invalid_channels:\n",
    "            raise ValueError(\n",
    "                f\"The following channels are not valid: {invalid_channels}\")\n",
    "\n",
    "        # Use only the specified channels\n",
    "        good.pick_channels(channels)\n",
    "\n",
    "    ch_type = filt.get_channel_types(only_data_chs=True)[0]\n",
    "    good.set_eeg_reference(ref_channels=\"average\", ch_type=ch_type)\n",
    "\n",
    "    all_epochs_list = []\n",
    "\n",
    "    for event in events:\n",
    "    # Epoching and HG extraction for each specified event. Then concatenate all trials epochs objects together (do Stimulus/c25 and Stimulus/c75 for example, and combine to get all congruent trials)\n",
    "        times_adj = [times[0] - 0.5, times[1] + 0.5]\n",
    "        trials = trial_ieeg(good, event, times_adj, preload=True,\n",
    "                            reject_by_annotation=False)\n",
    "        all_epochs_list.append(trials)\n",
    "\n",
    "    # Concatenate all trials\n",
    "    all_trials = mne.concatenate_epochs(all_epochs_list)\n",
    "\n",
    "    outliers_to_nan(all_trials, outliers=10)\n",
    "    HG_ev1 = gamma.extract(all_trials, copy=False, n_jobs=1)\n",
    "    crop_pad(HG_ev1, \"0.5s\")\n",
    "    HG_ev1_rescaled = rescale(HG_ev1, HG_base, copy=True, mode='zscore')\n",
    "\n",
    "    HG_base.decimate(2)\n",
    "    HG_ev1.decimate(2)\n",
    "    RTs, skipped = calculate_RTs(good)\n",
    "    avg_RT = np.median(RTs)\n",
    "\n",
    "    HG_ev1_evoke = HG_ev1.average(method=lambda x: np.nanmean(x, axis=0))\n",
    "    HG_ev1_evoke_rescaled = HG_ev1_rescaled.average(method=lambda x: np.nanmean(x, axis=0))\n",
    "\n",
    "    if event == \"Stimulus\":\n",
    "        print('plotting stimulus')\n",
    "        fig = HG_ev1_evoke_rescaled.plot(unit=False, scalings=dict(sEEG=1)) #this line is not finishing...\n",
    "        print('plotted')\n",
    "        for ax in fig.axes:\n",
    "            ax.axvline(x=avg_RT, color='r', linestyle='--')\n",
    "        print('about to save')\n",
    "        fig.savefig(save_dir + '_HG_ev1_Stimulus_zscore.png')\n",
    "        print('saved')\n",
    "    else:\n",
    "        print('about to plot if not stimulus')\n",
    "        fig = HG_ev1_evoke_rescaled.plot(unit=False, scalings=dict(sEEG=1))\n",
    "        print('plotted non stimulus')\n",
    "        fig.savefig(save_dir + f'_HG_ev1_{output_name}_zscore.png')\n",
    "\n",
    "    ###\n",
    "    print(f\"Shape of HG_ev1._data: {HG_ev1._data.shape}\")\n",
    "    print(f\"Shape of HG_base._data: {HG_base._data.shape}\")\n",
    "\n",
    "    sig1 = HG_ev1._data\n",
    "    sig2 = HG_base._data\n",
    "    sig2 = make_data_same(sig2, sig1.shape)\n",
    "    print(f\"Shape of sig1: {sig1.shape}\")\n",
    "    print(f\"Shape of sig2: {sig2.shape}\")\n",
    "\n",
    "    mat = time_perm_cluster(sig1, sig2, 0.05, n_jobs=6, ignore_adjacency=1)\n",
    "    fig = plt.figure()\n",
    "    plt.imshow(mat, aspect='auto')\n",
    "    fig.savefig(save_dir + f'_{output_name}_stats.png', dpi=300)\n",
    "\n",
    "    channels = good.ch_names\n",
    "\n",
    "    #save channels with their indices \n",
    "    save_channels_to_file(channels, sub, task, save_dir)\n",
    "\n",
    "    # save significant channels to a json\n",
    "    save_sig_chans(f'{output_name}', mat, channels, sub, save_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###\n",
    "import sys\n",
    "print(sys.path)\n",
    "sys.path.append(\"C:/Users/jz421/Desktop/GlobalLocal/IEEG_Pipelines/\") #need to do this cuz otherwise ieeg isn't added to path...\n",
    "\n",
    "from ieeg.navigate import channel_outlier_marker, trial_ieeg, crop_empty_data, \\\n",
    "    outliers_to_nan\n",
    "from ieeg.io import raw_from_layout, get_data\n",
    "from ieeg.timefreq.utils import crop_pad\n",
    "from ieeg.timefreq import gamma\n",
    "from ieeg.calc.scaling import rescale\n",
    "import mne\n",
    "import os\n",
    "import numpy as np\n",
    "from ieeg.calc.reshape import make_data_same\n",
    "from ieeg.calc.stats import time_perm_cluster\n",
    "\n",
    "from misc_functions import calculate_RTs, save_channels_to_file, save_sig_chans\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def get_baseline(inst: mne.io.BaseRaw, base_times: tuple[float, float]):\n",
    "    inst = inst.copy()\n",
    "    inst.load_data()\n",
    "    ch_type = inst.get_channel_types(only_data_chs=True)[0]\n",
    "    inst.set_eeg_reference(ref_channels=\"average\", ch_type=ch_type)\n",
    "\n",
    "    adjusted_base_times = [base_times[0] - 0.5, base_times[1] + 0.5]\n",
    "    trials = trial_ieeg(inst, \"experimentStart\", adjusted_base_times,\n",
    "                        preload=True)\n",
    "    # outliers_to_nan(trials, outliers=10)\n",
    "    HG_base = gamma.extract(trials, copy=False, n_jobs=1)\n",
    "    crop_pad(HG_base, \"0.5s\")\n",
    "    del inst\n",
    "    return HG_base\n",
    "\n",
    "\n",
    "def plot_HG_and_stats(sub, task, output_name, event=None, times=(-1, 1.5),\n",
    "                      base_times=(-0.5, 0), LAB_root=None, channels=None,\n",
    "                      full_trial_base=False):\n",
    "    \"\"\"\n",
    "    Plot high gamma (HG) and statistics for a given subject and task using specified event.\n",
    "\n",
    "    Parameters:\n",
    "    - sub (str): The subject identifier.\n",
    "    - task (str): The task identifier.\n",
    "    - output_name (str): The name for the output files.\n",
    "    - event (str, optional): Event name to process. Defaults to None.\n",
    "    - times (tuple, optional): A tuple indicating the start and end times for processing. Defaults to (-1, 1.5).\n",
    "    - base_times (tuple, optional): A tuple indicating the start and end base times for processing. Defaults to (-0.5, 0).\n",
    "    - LAB_root (str, optional): The root directory for the lab. Will be determined based on OS if not provided. Defaults to None.\n",
    "    - channels (list of strings, optional): The channels to plot and get stats for. Default is all channels.\n",
    "    - full_trial_base (boolean): Whether to use the full trial as the baseline period. Default is False.\n",
    "    This function will process the provided event for a given subject and task.\n",
    "    High gamma (HG) will be computed, and statistics will be calculated and plotted.\n",
    "    The results will be saved to output files.\n",
    "    \"\"\"\n",
    "    pass\n",
    "\n",
    "subjects = ['D0057', 'D0059', 'D0063', 'D0065', 'D0069', 'D0071']\n",
    "for sub in subjects:\n",
    "    sub = sub\n",
    "    output_name = \"response_experimentStartBase1secTo101Sec\"\n",
    "    events = [\"Response\"]\n",
    "    times = (-1, 1.5)\n",
    "    base_times = (1, 101)\n",
    "    LAB_root = None\n",
    "    channels = None\n",
    "    full_trial_base = False\n",
    "\n",
    "    if LAB_root is None:\n",
    "        HOME = os.path.expanduser(\"~\")\n",
    "        if os.name == 'nt':  # windows\n",
    "            LAB_root = os.path.join(HOME, \"Box\", \"CoganLab\")\n",
    "        else:  # mac\n",
    "            LAB_root = os.path.join(HOME, \"Library\", \"CloudStorage\", \"Box-Box\",\n",
    "                                    \"CoganLab\")\n",
    "\n",
    "    layout = get_data(task, root=LAB_root)\n",
    "    filt = raw_from_layout(layout.derivatives['derivatives/clean'], subject=sub,\n",
    "                        extension='.edf', desc='clean', preload=False)\n",
    "    save_dir = os.path.join(layout.root, 'derivatives', 'freqFilt', 'figs', sub)\n",
    "    if not os.path.exists(save_dir):\n",
    "        os.makedirs(save_dir)\n",
    "\n",
    "    ###\n",
    "\n",
    "    good = crop_empty_data(filt)\n",
    "    # %%\n",
    "\n",
    "    print(f\"good channels before dropping bads: {len(good.ch_names)}\")\n",
    "    print(f\"filt channels before dropping bads: {len(filt.ch_names)}\")\n",
    "\n",
    "    good.info['bads'] = channel_outlier_marker(good, 3, 2)\n",
    "    print(\"Bad channels in 'good':\", good.info['bads'])\n",
    "\n",
    "    filt.drop_channels(good.info['bads'])  # this has to come first cuz if you drop from good first, then good.info['bads'] is just empty\n",
    "    good.drop_channels(good.info['bads'])\n",
    "\n",
    "    print(\"Bad channels in 'good' after dropping once:\", good.info['bads'])\n",
    "\n",
    "    print(f\"good channels after dropping bads: {len(good.ch_names)}\")\n",
    "    print(f\"filt channels after dropping bads: {len(filt.ch_names)}\")\n",
    "\n",
    "    HG_base = get_baseline(filt, base_times)\n",
    "    good.load_data()\n",
    "\n",
    "    # If channels is None, use all channels\n",
    "    if channels is None:\n",
    "        channels = good.ch_names\n",
    "    else:\n",
    "        # Validate the provided channels\n",
    "        invalid_channels = [ch for ch in channels if ch not in good.ch_names]\n",
    "        if invalid_channels:\n",
    "            raise ValueError(\n",
    "                f\"The following channels are not valid: {invalid_channels}\")\n",
    "\n",
    "        # Use only the specified channels\n",
    "        good.pick_channels(channels)\n",
    "\n",
    "    ch_type = filt.get_channel_types(only_data_chs=True)[0]\n",
    "    good.set_eeg_reference(ref_channels=\"average\", ch_type=ch_type)\n",
    "\n",
    "    all_epochs_list = []\n",
    "\n",
    "    for event in events:\n",
    "    # Epoching and HG extraction for each specified event. Then concatenate all trials epochs objects together (do Stimulus/c25 and Stimulus/c75 for example, and combine to get all congruent trials)\n",
    "        times_adj = [times[0] - 0.5, times[1] + 0.5]\n",
    "        trials = trial_ieeg(good, event, times_adj, preload=True,\n",
    "                            reject_by_annotation=False)\n",
    "        all_epochs_list.append(trials)\n",
    "\n",
    "    # Concatenate all trials\n",
    "    all_trials = mne.concatenate_epochs(all_epochs_list)\n",
    "\n",
    "    outliers_to_nan(all_trials, outliers=10)\n",
    "    HG_ev1 = gamma.extract(all_trials, copy=False, n_jobs=1)\n",
    "    crop_pad(HG_ev1, \"0.5s\")\n",
    "    HG_ev1_rescaled = rescale(HG_ev1, HG_base, copy=True, mode='zscore')\n",
    "\n",
    "    HG_base.decimate(2)\n",
    "    HG_ev1.decimate(2)\n",
    "    RTs, skipped = calculate_RTs(good)\n",
    "    avg_RT = np.median(RTs)\n",
    "\n",
    "    HG_ev1_evoke = HG_ev1.average(method=lambda x: np.nanmean(x, axis=0))\n",
    "    HG_ev1_evoke_rescaled = HG_ev1_rescaled.average(method=lambda x: np.nanmean(x, axis=0))\n",
    "\n",
    "    if event == \"Stimulus\":\n",
    "        print('plotting stimulus')\n",
    "        fig = HG_ev1_evoke_rescaled.plot(unit=False, scalings=dict(sEEG=1)) #this line is not finishing...\n",
    "        print('plotted')\n",
    "        for ax in fig.axes:\n",
    "            ax.axvline(x=avg_RT, color='r', linestyle='--')\n",
    "        print('about to save')\n",
    "        fig.savefig(save_dir + '_HG_ev1_Stimulus_zscore.png')\n",
    "        print('saved')\n",
    "    else:\n",
    "        print('about to plot if not stimulus')\n",
    "        fig = HG_ev1_evoke_rescaled.plot(unit=False, scalings=dict(sEEG=1))\n",
    "        print('plotted non stimulus')\n",
    "        fig.savefig(save_dir + f'_HG_ev1_{output_name}_zscore.png')\n",
    "\n",
    "    ###\n",
    "    print(f\"Shape of HG_ev1._data: {HG_ev1._data.shape}\")\n",
    "    print(f\"Shape of HG_base._data: {HG_base._data.shape}\")\n",
    "\n",
    "    sig1 = HG_ev1._data\n",
    "    sig2 = HG_base._data\n",
    "    sig2 = make_data_same(sig2, sig1.shape)\n",
    "    print(f\"Shape of sig1: {sig1.shape}\")\n",
    "    print(f\"Shape of sig2: {sig2.shape}\")\n",
    "\n",
    "    mat = time_perm_cluster(sig1, sig2, 0.05, n_jobs=6, ignore_adjacency=1)\n",
    "    fig = plt.figure()\n",
    "    plt.imshow(mat, aspect='auto')\n",
    "    fig.savefig(save_dir + f'_{output_name}_stats.png', dpi=300)\n",
    "\n",
    "    channels = good.ch_names\n",
    "\n",
    "    #save channels with their indices \n",
    "    save_channels_to_file(channels, sub, task, save_dir)\n",
    "\n",
    "    # save significant channels to a json\n",
    "    save_sig_chans(f'{output_name}', mat, channels, sub, save_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### debug why the z-scores are negative below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "good.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HG_ev1.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ieeg",
   "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.11.3"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}