Forebrain_Manuscript_Figures_Resubmission.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "87a6c935",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pathlib\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image\n",
    "from matplotlib.axes import Axes\n",
    "from matplotlib.ticker import MultipleLocator\n",
    "import scipy.stats\n",
    "from scipy.fftpack import fft  # fourier transform\n",
    "from powltools.analysis.recording import group_by_param\n",
    "from mytools import (\n",
    "    plot_bracket,\n",
    "    anova_tukey,\n",
    "    condition_batch,\n",
    "    subplot_indicator,\n",
    "    figure_add_axes_inch,\n",
    "    figure_add_axes_group_inch,\n",
    ")\n",
    "from xcorr_tools import binary_spiketrain, make_psth_bins, get_peak\n",
    "from xcorr_tools import wrap_to_pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0195da18",
   "metadata": {},
   "outputs": [],
   "source": [
    "c_forebrain = [\"orange\", \"orange\", \"orange\", \"orange\", \"orange\", \"orange\"]\n",
    "c_ot = [\"b\", \"b\", \"b\", \"b\", \"b\", \"b\"]\n",
    "\n",
    "c_forebrain_am55 = [\n",
    "    \"mediumseagreen\",\n",
    "    \"mediumseagreen\",\n",
    "    \"mediumseagreen\",\n",
    "    \"mediumseagreen\",\n",
    "    \"mediumseagreen\",\n",
    "    \"mediumseagreen\",\n",
    "]\n",
    "c_forebrain_am75 = [\n",
    "    \"cornflowerblue\",\n",
    "    \"cornflowerblue\",\n",
    "    \"cornflowerblue\",\n",
    "    \"cornflowerblue\",\n",
    "    \"cornflowerblue\",\n",
    "    \"cornflowerblue\",\n",
    "]\n",
    "\n",
    "# OUTDIR = pathlib.Path(\"./forebrain_figures_output\").absolute()\n",
    "OUTDIR = pathlib.Path(\"./forebrain_output_resubmission\").absolute()\n",
    "\n",
    "OUTDIR.mkdir(exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d473f90e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# matplotlib settings\n",
    "mpl.rcParams[\"mathtext.default\"] = \"regular\"\n",
    "mpl.rcParams[\"font.family\"] = \"sans-serif\"\n",
    "mpl.rcParams[\"font.sans-serif\"] = \"Arial\"\n",
    "\n",
    "# Reproducible SVG files (including xml ids):\n",
    "mpl.rcParams[\"svg.hashsalt\"] = \"123\"\n",
    "mpl.rcParams[\"savefig.dpi\"] = \"300\"\n",
    "\n",
    "# Font sizes:\n",
    "mpl.rcParams[\"font.size\"] = 12\n",
    "mpl.rcParams[\"axes.labelsize\"] = 12  # default: 'medium' == 10\n",
    "mpl.rcParams[\"xtick.labelsize\"] = 10  # default: 'medium' == 10\n",
    "mpl.rcParams[\"ytick.labelsize\"] = 10  # default: 'medium' == 10\n",
    "mpl.rcParams[\"legend.fontsize\"] = 10\n",
    "mpl.rcParams[\"xtick.major.size\"] = 2\n",
    "mpl.rcParams[\"ytick.major.size\"] = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6300f839",
   "metadata": {},
   "source": [
    "# Load Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "384cb56c",
   "metadata": {},
   "outputs": [],
   "source": [
    "FEATHER_DIR = pathlib.Path(\"./forebrain_intermediate_results\")\n",
    "FEATHER_DIR.absolute().mkdir(exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8320be82",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Flat Noise Competition\n",
    "\n",
    "single_srf = pd.read_feather(FEATHER_DIR / \"single_srf_dualregion.feather\").set_index(\n",
    "    [\"date\", \"owl\", \"channel\"]\n",
    ")\n",
    "\n",
    "\n",
    "single_rlf = pd.read_feather(FEATHER_DIR / \"single_rlf_dualregion.feather\").set_index(\n",
    "    [\"date\", \"owl\", \"channel\"]\n",
    ")\n",
    "\n",
    "single_rlf_out = pd.read_feather(\n",
    "    FEATHER_DIR / \"single_rlf_out_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel\"])\n",
    "\n",
    "\n",
    "single_ccg = pd.read_feather(FEATHER_DIR / \"single_ccg_dualregion.feather\").set_index(\n",
    "    [\"date\", \"owl\", \"channel1\", \"channel2\"]\n",
    ")\n",
    "\n",
    "single_ccg_out = pd.read_feather(\n",
    "    FEATHER_DIR / \"single_ccg_dualregion_out.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel1\", \"channel2\"])\n",
    "\n",
    "\n",
    "single_gamma_power = pd.read_feather(\n",
    "    FEATHER_DIR / \"single_gamma_power_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel\"])\n",
    "\n",
    "\n",
    "within_area_sfc = pd.read_feather(\n",
    "    FEATHER_DIR / \"twostim_sfc_within_area_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel\"])\n",
    "\n",
    "\n",
    "twostim_rlf = pd.read_feather(FEATHER_DIR / \"twostim_rlf_dualregion.feather\").set_index(\n",
    "    [\"date\", \"owl\", \"channel\"]\n",
    ")\n",
    "\n",
    "twostim_rlf_switch = pd.read_feather(\n",
    "    FEATHER_DIR / \"twostim_rlf_switch_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel\"])\n",
    "\n",
    "\n",
    "twostim_ccg = pd.read_feather(FEATHER_DIR / \"twostim_ccg_dualregion.feather\").set_index(\n",
    "    [\"date\", \"owl\", \"channel1\", \"channel2\"]\n",
    ")\n",
    "\n",
    "\n",
    "twostim_ccg_switch = pd.read_feather(\n",
    "    FEATHER_DIR / \"twostim_ccg_switch_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel1\", \"channel2\"])\n",
    "\n",
    "\n",
    "cross_region_sfc = pd.read_feather(\n",
    "    FEATHER_DIR / \"sfc_cross_region_df_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel1\", \"channel2\"])\n",
    "\n",
    "\n",
    "twostim_gamma_power = pd.read_feather(\n",
    "    FEATHER_DIR / \"twostim_gamma_power_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel\"])\n",
    "\n",
    "\n",
    "## AM Competition\n",
    "am_single_rlf = pd.read_feather(\n",
    "    FEATHER_DIR / \"am_single_rlf_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel\"])\n",
    "\n",
    "am_single_ccg = pd.read_feather(\n",
    "    FEATHER_DIR / \"am_single_ccg_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel1\", \"channel2\"])\n",
    "am_single_stim_phaselocking = pd.read_feather(\n",
    "    FEATHER_DIR / \"am_single_stim_phaselocking_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel\"])\n",
    "\n",
    "\n",
    "cross_region_gamma = pd.read_feather(\n",
    "    FEATHER_DIR / \"cross_region_gamma_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel1\", \"channel2\"])\n",
    "\n",
    "am_single_gamma_power = pd.read_feather(\n",
    "    FEATHER_DIR / \"am_single_gamma_power_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel\"])\n",
    "\n",
    "am_twostim_rlf = pd.read_feather(\n",
    "    FEATHER_DIR / \"am_twostim_rlf_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel\"])\n",
    "am_twostim_ccg = pd.read_feather(\n",
    "    FEATHER_DIR / \"am_twostim_ccg_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel1\", \"channel2\"])\n",
    "am_twostim_stim_phaselocking = pd.read_feather(\n",
    "    FEATHER_DIR / \"am_twostim_stim_phaselocking_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel\"])\n",
    "am_twostim_gamma_power = pd.read_feather(\n",
    "    FEATHER_DIR / \"am_twostim_gamma_power_dualregion.feather\"\n",
    ").set_index([\"date\", \"owl\", \"channel\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e87a963c",
   "metadata": {},
   "outputs": [],
   "source": [
    "## example data\n",
    "example_srf = pd.read_feather(FEATHER_DIR / \"example_srf_20240315_256.feather\")\n",
    "example_spiketrains_forebrain_flat = pd.read_feather(\n",
    "    FEATHER_DIR / \"example_spiketrains_forebrain_flat_20240315_256.feather\"\n",
    ")\n",
    "\n",
    "example_spiketrains_am55 = pd.read_feather(\n",
    "    FEATHER_DIR / \"example_spiketrains_forebrain_driver55_20240330_40.feather\"\n",
    ")\n",
    "example_spiketrains_am75 = pd.read_feather(\n",
    "    FEATHER_DIR / \"example_spiketrains_forebrain_driver75_20240330_40.feather\"\n",
    ")\n",
    "\n",
    "\n",
    "lowcorr_srf_ot = pd.read_feather(\n",
    "    FEATHER_DIR / \"example_srf_lowcorr_OT_20240420_54.feather\"\n",
    ")\n",
    "lowcorr_srf_fb = pd.read_feather(\n",
    "    FEATHER_DIR / \"example_srf_lowcorr_Forebrain_20240420_54.feather\"\n",
    ")\n",
    "\n",
    "midlowcorr_srf_ot = pd.read_feather(\n",
    "    FEATHER_DIR / \"example_srf_midlowcorr_OT_20240423_54.feather\"\n",
    ")\n",
    "midlowcorr_srf_fb = pd.read_feather(\n",
    "    FEATHER_DIR / \"example_srf_midlowcorr_Forebrain_20240423_54.feather\"\n",
    ")\n",
    "\n",
    "midhighcorr_srf_ot = pd.read_feather(\n",
    "    FEATHER_DIR / \"example_srf_midhighcorr_OT_20240420_54.feather\"\n",
    ")\n",
    "midhighcorr_srf_fb = pd.read_feather(\n",
    "    FEATHER_DIR / \"example_srf_midhighcorr_Forebrain_20240420_54.feather\"\n",
    ")\n",
    "\n",
    "highcorr_srf_ot = pd.read_feather(\n",
    "    FEATHER_DIR / \"example_srf_highcorr_OT_20240420_54.feather\"\n",
    ")\n",
    "highcorr_srf_fb = pd.read_feather(\n",
    "    FEATHER_DIR / \"example_srf_highcorr_Forebrain_20240420_54.feather\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aba1f655",
   "metadata": {},
   "source": [
    "# Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "17536bdc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def confidence_intervale(sample_data):\n",
    "    # Sample statistics\n",
    "    sample_mean = np.mean(sample_data)\n",
    "    sample_std = np.std(sample_data, ddof=1)\n",
    "    n = len(sample_data)\n",
    "    confidence_level = 0.95\n",
    "\n",
    "    # Calculate standard error (population standard deviation unknown)\n",
    "    standard_error = sample_std / np.sqrt(n)\n",
    "\n",
    "    # Degrees of freedom for t-distribution\n",
    "    df = n - 1\n",
    "\n",
    "    # Calculate confidence interval using t.interval()\n",
    "    confidence_interval = scipy.stats.t.interval(\n",
    "        confidence_level, df, loc=sample_mean, scale=standard_error\n",
    "    )\n",
    "\n",
    "    print(f\"Confidence Interval: {confidence_interval}\")\n",
    "    return confidence_interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5dedc9de",
   "metadata": {},
   "outputs": [],
   "source": [
    "def phase_properties(phases):\n",
    "    ## convert to complex numbers\n",
    "    pi_bin = np.linspace(-np.pi, np.pi, 9)\n",
    "    #    pi_bin = np.arange(-np.pi, np.pi + np.pi / 6, np.pi / 6)\n",
    "    sfc = np.histogram(phases, pi_bin)\n",
    "    bins = sfc[1]\n",
    "\n",
    "    ##calculate mean angle of frequency grouped data\n",
    "    y_angle = sum(sfc[0] * np.sin(sfc[1][:-1])) / sum(sfc[0])\n",
    "    x_angle = sum(sfc[0] * np.cos(sfc[1][:-1])) / sum(sfc[0])\n",
    "    # r_val = np.sqrt((y_angle**2) + (x_angle**2))\n",
    "    mean_ang = np.arctan2((y_angle), (x_angle))\n",
    "\n",
    "    r_val = np.sqrt((y_angle**2) + (x_angle**2))\n",
    "    d = np.diff(bins)[0]\n",
    "    corr_fac = d / 2 / np.sin(d / 2)\n",
    "    r_val = corr_fac * r_val\n",
    "\n",
    "    total_n = len(phases)\n",
    "    rayleigh_r = total_n * r_val\n",
    "    pval = np.exp(\n",
    "        np.sqrt(1 + 4 * total_n + 4 * (total_n**2 - rayleigh_r**2)) - (1 + 2 * total_n)\n",
    "    )\n",
    "    return {\"vector_strength\": r_val, \"mean_angle\": mean_ang, \"pval\": pval}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a573089a",
   "metadata": {},
   "source": [
    "# Plotting + Analysis Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7ef28e6",
   "metadata": {},
   "source": [
    "## 2.1 example spatial receptive field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9f798bfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def figure_spatial_receptive_field(\n",
    "    srf_data: pd.DataFrame,\n",
    "    # rec: Recording,\n",
    "    # channel_number: int,\n",
    "    ax=None,\n",
    "    cax=None,\n",
    "):\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots()\n",
    "    else:\n",
    "        fig = ax.figure\n",
    "\n",
    "    ax.xaxis.set_major_locator(MultipleLocator(30))\n",
    "    ax.set_xlabel(\"Azimuth [deg]\")\n",
    "    ax.yaxis.set_major_locator(MultipleLocator(30))\n",
    "    ax.set_ylabel(\"Elevation [deg]\")\n",
    "\n",
    "    azimuths = srf_data[\"azimuth\"].values\n",
    "    elevations = srf_data[\"elevation\"].values\n",
    "    response_rates = srf_data[\"response_mean\"].values\n",
    "\n",
    "    cmap = mpl.colormaps[\"Greys\"]\n",
    "    srf = ax.scatter(\n",
    "        x=azimuths,\n",
    "        y=elevations,\n",
    "        c=response_rates,\n",
    "        cmap=cmap,\n",
    "        edgecolors=\"0.7\",\n",
    "        vmin=0.0,\n",
    "        linewidths=0.5,\n",
    "    )\n",
    "    cb = fig.colorbar(srf, cax=cax)\n",
    "    cb.set_label(\"Response [spks/stim]\")\n",
    "    return"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1e6db4a",
   "metadata": {},
   "source": [
    "## Azimuth tuning heatmap + population average"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f60e2d87",
   "metadata": {},
   "outputs": [],
   "source": [
    "def azimuth_heatplot_alt(single_srf_df, axs=None, caxs=None):\n",
    "\n",
    "    if axs is None:\n",
    "        fig, axs = plt.subplots(1, 2)\n",
    "\n",
    "    srf_df = single_srf_df.copy()\n",
    "    rlf_df = single_rlf_df.copy()\n",
    "\n",
    "    itd_range = np.linspace(-100, 100, 21)\n",
    "    itd_response_hm = np.empty([len(np.unique(srf_df.index)), len(itd_range)])\n",
    "\n",
    "    for krow, ind in enumerate(np.unique(srf_df.index)):\n",
    "        unit_df = srf_df.loc[ind]\n",
    "        means = unit_df.groupby(\"azimuth\")[\"norm_resp\"].max()\n",
    "\n",
    "        azimuths = means.index\n",
    "        spkresponse = means.values\n",
    "\n",
    "        if unit_df.iloc[0][\"hemisphere\"] == \"right\":\n",
    "            itd_response_hm[krow, :] = spkresponse[::-1]\n",
    "        else:\n",
    "            itd_response_hm[krow, :] = spkresponse\n",
    "    print(f\"nRT units {itd_response_hm.shape[0]}\")\n",
    "    cmap_reversed = matplotlib.cm.get_cmap(\"hot\")\n",
    "    plot_xticks = np.linspace(-100, 100, 5)\n",
    "    sc = axs[0].imshow(\n",
    "        itd_response_hm,\n",
    "        extent=[-100, 100, 0, itd_response_hm.shape[0]],\n",
    "        aspect=\"auto\",\n",
    "        cmap=cmap_reversed,\n",
    "        vmin=0,\n",
    "        vmax=1.0,\n",
    "    )\n",
    "    axs[0].set_xticks(\n",
    "        [int(val) for val in plot_xticks], [int(val) for val in plot_xticks]\n",
    "    )\n",
    "    axs[0].set_ylabel(\"Units\")\n",
    "    axs[0].set_xlabel(\"Azimuth\")\n",
    "    cb = plt.colorbar(sc, cax=caxs)\n",
    "    cb.set_label(\"Norm. Response\")\n",
    "\n",
    "    mean_azimuth_curve = np.mean(itd_response_hm, axis=0)\n",
    "    sem = scipy.stats.sem(itd_response_hm, axis=0)\n",
    "\n",
    "    axs[1].errorbar(itd_range, mean_azimuth_curve, sem, color=\"k\")\n",
    "    axs[1].set_ylim(bottom=0.30, top=0.75)\n",
    "    axs[1].set_xlabel(\"Azimuth\")\n",
    "    axs[1].set_ylabel(\"Norm. Response\")\n",
    "    axs[1].spines[\"top\"].set_visible(False)\n",
    "    axs[1].spines[\"right\"].set_visible(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "139fe51b",
   "metadata": {},
   "source": [
    "## Elevation tuning heatmap and population average"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "17ec5ba5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def elevation_heatplot_alt(single_srf_df, axs=None, caxs=None):\n",
    "\n",
    "    if axs is None:\n",
    "        fig, axs = plt.subplots(1, 2)\n",
    "\n",
    "    srf_df = single_srf_df.copy()\n",
    "\n",
    "    ele_range = np.linspace(-80, 80, 17)\n",
    "    ele_response_hm = np.empty([len(np.unique(srf_df.index)), len(ele_range)])\n",
    "\n",
    "    for krow, ind in enumerate(np.unique(srf_df.index)):\n",
    "        unit_df = srf_df.loc[ind]\n",
    "        means = unit_df.groupby(\"elevation\")[\"norm_resp\"].max()\n",
    "\n",
    "        elevations = means.index\n",
    "        spkresponse = means.values\n",
    "\n",
    "        ele_response_hm[krow, :] = spkresponse\n",
    "\n",
    "    cmap_reversed = matplotlib.cm.get_cmap(\"hot\")\n",
    "\n",
    "    plot_xticks = np.linspace(-80, 80, 5)\n",
    "    sc = axs[0].imshow(\n",
    "        ele_response_hm,\n",
    "        extent=[-80, 80, 0, ele_response_hm.shape[0]],\n",
    "        aspect=\"auto\",\n",
    "        cmap=cmap_reversed,\n",
    "        vmin=0,\n",
    "        vmax=1.0,\n",
    "    )\n",
    "    axs[0].set_xticks(\n",
    "        [int(val) for val in plot_xticks], [int(val) for val in plot_xticks]\n",
    "    )\n",
    "    axs[0].set_ylabel(\"Units\")\n",
    "    axs[0].set_xlabel(\"Elevation\")\n",
    "    cb = plt.colorbar(sc, cax=caxs)\n",
    "    cb.set_label(\"Norm. Response\")\n",
    "\n",
    "    mean_ele_curve = np.mean(ele_response_hm, axis=0)\n",
    "    sem = scipy.stats.sem(ele_response_hm, axis=0)\n",
    "\n",
    "    axs[1].errorbar(ele_range, mean_ele_curve, sem, color=\"k\")\n",
    "    axs[1].set_xlabel(\"Elevation\")\n",
    "    axs[1].set_ylabel(\"Norm. Response\")\n",
    "    axs[1].set_ylim(bottom=0.30, top=0.75)\n",
    "    axs[1].spines[\"top\"].set_visible(False)\n",
    "    axs[1].spines[\"right\"].set_visible(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0852933",
   "metadata": {},
   "source": [
    "## competition boxplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "eb58628d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def figure_competition_boxplot(\n",
    "    twostim_rlf: pd.DataFrame,\n",
    "    single_rlf: pd.DataFrame,\n",
    "    ax=None,\n",
    "    colors=c_forebrain,\n",
    "    brackets: dict[tuple[int, int] : float] | None = None,\n",
    "    anova_align=\"right\",\n",
    "):\n",
    "    \"\"\"Boxplot of Change in Spike Rate as function of Relative Level\"\"\"\n",
    "\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots()\n",
    "\n",
    "    merged_rlf = twostim_rlf.join(single_rlf, how=\"inner\", rsuffix=\"_single\")\n",
    "    merged_rlf = merged_rlf.loc[\n",
    "        (merged_rlf[\"fixedintensity\"] == merged_rlf[\"intensity\"])\n",
    "    ]\n",
    "    merged_rlf = merged_rlf.loc[(merged_rlf[\"fixedazi\"] == merged_rlf[\"azimuth\"])]\n",
    "    merged_rlf = merged_rlf.loc[(merged_rlf[\"fixedele\"] == merged_rlf[\"elevation\"])]\n",
    "\n",
    "    print(\"Number of units:\", np.unique(merged_rlf.index).size)\n",
    "    print(\n",
    "        \"Number of sessions:\",\n",
    "        len(sorted(set([idx[:2] for idx in merged_rlf.index.values]))),\n",
    "    )\n",
    "    merged_rlf.set_index(\"relative_level\", inplace=True)\n",
    "    # display(merged_rlf_ot)\n",
    "\n",
    "    merged_rlf[\"change_in_response\"] = (\n",
    "        merged_rlf[\"resp\"] - merged_rlf[\"resp_single\"]\n",
    "    ) / merged_rlf[\"resp_single\"]\n",
    "    relative_levels = np.unique(merged_rlf.index)\n",
    "    groupdata = [\n",
    "        merged_rlf.loc[relative_level, \"change_in_response\"].values\n",
    "        for relative_level in relative_levels\n",
    "    ]\n",
    "    for k, lev in enumerate(relative_levels):\n",
    "        print(\n",
    "            f\"{lev}: {np.median(groupdata[k])} , sem: {scipy.stats.sem(groupdata[k])}\"\n",
    "        )\n",
    "\n",
    "    bp = ax.boxplot(\n",
    "        groupdata,\n",
    "        patch_artist=True,\n",
    "        notch=False,\n",
    "        showfliers=False,\n",
    "        whis=(5, 95),\n",
    "    )\n",
    "\n",
    "    ax.set_xticklabels(relative_levels.astype(int))\n",
    "\n",
    "    stats = anova_tukey(\n",
    "        merged_rlf, val_col=\"change_in_response\", group_col=\"relative_level\"\n",
    "    )\n",
    "\n",
    "    for patch, color in zip(bp[\"boxes\"], colors):\n",
    "        patch.set_facecolor(color)\n",
    "\n",
    "    for median in bp[\"medians\"]:\n",
    "        median.set(color=\"white\", linewidth=2)\n",
    "\n",
    "    if True:  # show_individual_data\n",
    "        width = 0.7\n",
    "        for k, groupvalues in enumerate(groupdata):\n",
    "            pos = k + 1\n",
    "            y = np.sort(groupvalues)\n",
    "            vals_x = np.linspace(\n",
    "                pos - width / 3,\n",
    "                pos + width / 3,\n",
    "                y.size,\n",
    "                endpoint=True,\n",
    "            )\n",
    "            ax.plot(\n",
    "                vals_x,\n",
    "                y,\n",
    "                color=\"k\",\n",
    "                ls=\"None\",\n",
    "                marker=\".\",\n",
    "                markersize=2,\n",
    "                markeredgewidth=0.0,\n",
    "                zorder=5,\n",
    "            )\n",
    "\n",
    "    if True:  # show_brackets\n",
    "        brackets_shrink = 1.0\n",
    "        if brackets is None:\n",
    "            brackets = {(k, k + 1): 0.8 for k in range(len(stats[\"tukey\"]) - 1)}\n",
    "            brackets_shrink = 0.8\n",
    "        for (idx_left, idx_right), y in brackets.items():\n",
    "            plot_bracket(\n",
    "                ax,\n",
    "                left=idx_left + 1,\n",
    "                right=idx_right + 1,\n",
    "                text=f\"{stats['tukey'].iloc[idx_left, idx_right]:.3f}\"[1:],\n",
    "                y=y,\n",
    "                shrink=brackets_shrink,\n",
    "            )\n",
    "\n",
    "    ax.spines[\"top\"].set_visible(False)\n",
    "    ax.spines[\"right\"].set_visible(False)\n",
    "\n",
    "    if anova_align == \"right\":\n",
    "        ax.text(\n",
    "            1,\n",
    "            1,\n",
    "            f\"p = {stats['anova_p']:.3g}\",\n",
    "            horizontalalignment=\"right\",\n",
    "            verticalalignment=\"top\",\n",
    "            fontsize=10,\n",
    "            transform=ax.transAxes,\n",
    "        )\n",
    "    else:\n",
    "        ax.text(\n",
    "            0,\n",
    "            1,\n",
    "            f\"p = {stats['anova_p']:.3g}\",\n",
    "            horizontalalignment=\"left\",\n",
    "            verticalalignment=\"top\",\n",
    "            fontsize=10,\n",
    "            transform=ax.transAxes\n",
    "            + mpl.transforms.ScaledTranslation(\n",
    "                +2 / 72,\n",
    "                0 / 72,\n",
    "                ax.figure.dpi_scale_trans,\n",
    "            ),\n",
    "        )\n",
    "\n",
    "    ax.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, pos: f\"{x:.0%}\"))\n",
    "\n",
    "    ax.axhline(0, color=\"k\", ls=\"--\", zorder=-1)\n",
    "    ax.set_xlabel(\"Relative Level [dB]\")\n",
    "    ax.set_ylabel(\"Change in\\nSpike Rate\")\n",
    "    ax.set_ylim(-1, 1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cd4dbc0",
   "metadata": {},
   "source": [
    "## Raster Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a2359178",
   "metadata": {},
   "outputs": [],
   "source": [
    "def figure_coincident_rasterplot(\n",
    "    spiketrain_data: pd.DataFrame,\n",
    "    # filename: str,\n",
    "    # chan1: int,\n",
    "    # chan2: int,\n",
    "    relative_levels=[-15, +10],\n",
    "    axs=None,\n",
    "):\n",
    "\n",
    "    if axs is None:\n",
    "\n",
    "        fig, ax = plt.subplots(1, len(relative_levels))\n",
    "    axs = axs.flatten()\n",
    "\n",
    "    spiketrains_unit1 = spiketrain_data[\"spiketrain_unit1\"].values\n",
    "    spiketrains_unit2 = spiketrain_data[\"spiketrain_unit2\"].values\n",
    "    trial_relative_levels = spiketrain_data[\"relative_level\"]\n",
    "\n",
    "    stim_spikes1 = group_by_param(\n",
    "        np.array(\n",
    "            spiketrains_unit1,\n",
    "            dtype=\"object\",\n",
    "        ),\n",
    "        trial_relative_levels,\n",
    "    )\n",
    "    stim_spikes2 = group_by_param(\n",
    "        np.array(\n",
    "            spiketrains_unit2,\n",
    "            dtype=\"object\",\n",
    "        ),\n",
    "        trial_relative_levels,\n",
    "    )\n",
    "\n",
    "    psth_bins = make_psth_bins(\n",
    "        stimdelay=1.0,\n",
    "        stimduration=1.0,\n",
    "        binsize=0.001,\n",
    "        offset=0.000,\n",
    "    )\n",
    "\n",
    "    for k, relative_level in enumerate(relative_levels):\n",
    "        ax = axs[k]\n",
    "        spiketrains1 = stim_spikes1[relative_level][:10]\n",
    "        spiketrains2 = stim_spikes2[relative_level][:10]\n",
    "        coincidences = [\n",
    "            np.where(\n",
    "                binary_spiketrain(st1, bins=psth_bins)\n",
    "                & binary_spiketrain(st2, bins=psth_bins)\n",
    "            )[0]\n",
    "            * 0.001\n",
    "            + psth_bins[0]\n",
    "            for st1, st2 in zip(spiketrains1, spiketrains2)\n",
    "        ]\n",
    "        ax.eventplot(\n",
    "            spiketrains1,\n",
    "            lineoffsets=4 * np.arange(spiketrains1.size),\n",
    "            linelengths=1,\n",
    "            color=\"#00dd33\",\n",
    "            alpha=0.3,\n",
    "        )\n",
    "        ax.eventplot(\n",
    "            spiketrains2,\n",
    "            lineoffsets=2 + 4 * np.arange(spiketrains2.size),\n",
    "            linelengths=1,\n",
    "            color=\"#ffaa00\",\n",
    "            alpha=0.3,\n",
    "        )\n",
    "        ax.eventplot(\n",
    "            coincidences,\n",
    "            lineoffsets=1 + 4 * np.arange(len(coincidences)),\n",
    "            linelengths=1.4,\n",
    "            color=\"k\",\n",
    "        )\n",
    "        ax.axvspan(1.0, 1.050, color=\"#ddf\", zorder=-50)\n",
    "        ax.set_xlim(left=1.0, right=2.0)\n",
    "        ax.set_xticks([1.0, 1.5, 2.0], [\"0\", \"0.5\", \"1\"], fontsize=8)\n",
    "        ax.set_ylim(bottom=-0.5, top=2 + 4 * (spiketrains2.size - 1) + 0.5)\n",
    "\n",
    "    for ax in axs[1:]:\n",
    "        ax.set_yticks([])\n",
    "\n",
    "    axs[0].set_xlabel(\"Time [s]\", fontsize=8)\n",
    "    ytick_trial = 1\n",
    "    axs[0].set_yticks(\n",
    "        [4 * ytick_trial, 2 + 4 * ytick_trial, 1 + 4 * (ytick_trial + 2)],\n",
    "        [\"\", \"\", \"\"],\n",
    "    )\n",
    "    axs[0].tick_params(axis=\"y\", length=4)\n",
    "    axs[0].text(\n",
    "        0,\n",
    "        4 * ytick_trial,\n",
    "        \"Unit A\",\n",
    "        fontsize=10,\n",
    "        ha=\"right\",\n",
    "        va=\"top\",\n",
    "        color=\"g\",\n",
    "        transform=mpl.transforms.blended_transform_factory(\n",
    "            axs[0].transAxes, axs[0].transData\n",
    "        )\n",
    "        + mpl.transforms.ScaledTranslation(\n",
    "            -4 / 72,\n",
    "            0,\n",
    "            axs[0].figure.dpi_scale_trans,\n",
    "        ),\n",
    "    )\n",
    "\n",
    "    axs[0].text(\n",
    "        0,\n",
    "        2 + 4 * ytick_trial,\n",
    "        \"Unit B\",\n",
    "        fontsize=10,\n",
    "        ha=\"right\",\n",
    "        va=\"bottom\",\n",
    "        color=\"orange\",\n",
    "        transform=mpl.transforms.blended_transform_factory(\n",
    "            axs[0].transAxes, axs[0].transData\n",
    "        )\n",
    "        + mpl.transforms.ScaledTranslation(\n",
    "            -4 / 72,\n",
    "            0,\n",
    "            axs[0].figure.dpi_scale_trans,\n",
    "        ),\n",
    "    )\n",
    "\n",
    "    axs[0].text(\n",
    "        0,\n",
    "        1 + 4 * (ytick_trial + 2),\n",
    "        \"Coinc.\",\n",
    "        fontsize=10,\n",
    "        ha=\"right\",\n",
    "        va=\"bottom\",\n",
    "        color=\"k\",\n",
    "        transform=mpl.transforms.blended_transform_factory(\n",
    "            axs[0].transAxes, axs[0].transData\n",
    "        )\n",
    "        + mpl.transforms.ScaledTranslation(\n",
    "            -4 / 72,\n",
    "            0,\n",
    "            axs[0].figure.dpi_scale_trans,\n",
    "        ),\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef274bbf",
   "metadata": {},
   "source": [
    "## CCGs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b023a76c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def figure_xcorr_ccg(\n",
    "    twostim_ccg: pd.DataFrame, single_ccg: pd.DataFrame, axs=None, colors=c_forebrain\n",
    "):\n",
    "    \"\"\"Boxplot of Change in Spike Rate as function of Relative Level\"\"\"\n",
    "\n",
    "    if axs is None:\n",
    "        fig, axs = plt.subplots(1, 6)\n",
    "    axs = axs.flatten()\n",
    "\n",
    "    merged_ccg = twostim_ccg.join(single_ccg, how=\"inner\", rsuffix=\"_single\")\n",
    "    merged_ccg = merged_ccg.loc[\n",
    "        (merged_ccg[\"fixedintensity\"] == merged_ccg[\"intensity\"])\n",
    "    ]\n",
    "    merged_ccg = merged_ccg.loc[(merged_ccg[\"fixedazi\"] == merged_ccg[\"azimuth\"])]\n",
    "    merged_ccg = merged_ccg.loc[(merged_ccg[\"fixedele\"] == merged_ccg[\"elevation\"])]\n",
    "\n",
    "    print(\"Number of unit pairs:\", np.unique(merged_ccg.index).size)\n",
    "    print(\n",
    "        \"Number of sessions:\",\n",
    "        len(sorted(set([idx[:2] for idx in merged_ccg.index.values]))),\n",
    "    )\n",
    "\n",
    "    # Before reindexing by relative_level:\n",
    "    single_mean_peak = np.max(\n",
    "        np.mean(\n",
    "            np.vstack(merged_ccg[\"ccg_single\"].groupby(merged_ccg.index.names).first()),\n",
    "            axis=0,\n",
    "        )\n",
    "    )\n",
    "    single_std_peak = np.std(\n",
    "        np.max(\n",
    "            np.vstack(merged_ccg[\"ccg_single\"].groupby(merged_ccg.index.names).first()),\n",
    "            axis=0,\n",
    "        )\n",
    "    )\n",
    "    print(f\"mean: {single_mean_peak = :.4g} std: {single_std_peak}\")\n",
    "\n",
    "    merged_ccg.set_index(\"relative_level\", inplace=True)\n",
    "\n",
    "    psth_len = (merged_ccg.iloc[0][\"ccg_single\"].size + 1) / 2\n",
    "    lags = scipy.signal.correlation_lags(psth_len, psth_len) / 1000\n",
    "    lags_mask = np.abs(lags) <= 0.05\n",
    "\n",
    "    relative_levels = np.unique(merged_ccg.index)\n",
    "\n",
    "    for k, relative_level in enumerate(relative_levels):\n",
    "        mean_ccg = np.mean(np.vstack(merged_ccg.loc[relative_level, \"ccg\"]), axis=0)\n",
    "        n_ccg = merged_ccg.loc[relative_level, \"xcorr_peak_single\"].count()\n",
    "\n",
    "        axs[k].axhline(single_mean_peak, lw=1, color=colors[0])\n",
    "        print(\n",
    "            f\"{relative_level:5} | {np.mean(merged_ccg.loc[relative_level, 'xcorr_peak_single']):.4g} {n_ccg = }\"\n",
    "            f\"| {np.mean(merged_ccg.loc[relative_level, 'xcorr_peak']):.4g} {n_ccg = }\"\n",
    "            f\"| peak curve: {np.max(mean_ccg)}\"\n",
    "        )\n",
    "        axs[k].plot(lags[lags_mask], mean_ccg[lags_mask], color=\"k\", ls=\"-\", lw=1)\n",
    "        axs[k].text(\n",
    "            0.5,\n",
    "            1,\n",
    "            f\"{n_ccg}\",\n",
    "            ha=\"center\",\n",
    "            va=\"top\",\n",
    "            fontsize=8,\n",
    "            transform=axs[k].transAxes,\n",
    "        )\n",
    "        axs[k].spines[\"top\"].set_visible(False)\n",
    "        axs[k].spines[\"right\"].set_visible(False)\n",
    "\n",
    "        axs[k].set_ylim(bottom=0, top=0.000075)\n",
    "        axs[k].set_xlim(left=-0.060, right=0.060)\n",
    "\n",
    "    axs[0].set_ylabel(\"Coinc./spk\", labelpad=12)\n",
    "\n",
    "    def formatter_y(x, pos):\n",
    "        if x > 0:\n",
    "            return f\"${{{x*1e5:.0f}}}$\"\n",
    "            # return f\"$10^{{{np.log10(x):.0f}}}$\"\n",
    "        else:\n",
    "            return \"0\"\n",
    "\n",
    "    axs[0].set_yticks([0, 5e-5])\n",
    "    axs[0].yaxis.set_major_formatter(mpl.ticker.FuncFormatter(formatter_y))\n",
    "    axs[0].set_yticks(np.array([1, 2, 3, 4, 6, 7]) * 1e-5, minor=True)\n",
    "    axs[0].text(\n",
    "        0,\n",
    "        1,\n",
    "        \"$10^{-5}$\",\n",
    "        ha=\"right\",\n",
    "        va=\"center\",\n",
    "        fontsize=8,\n",
    "        transform=axs[0].transAxes\n",
    "        + mpl.transforms.ScaledTranslation(\n",
    "            -2 / 72,\n",
    "            0,\n",
    "            axs[0].figure.dpi_scale_trans,\n",
    "        ),\n",
    "    )\n",
    "\n",
    "    for ax in axs[1:]:\n",
    "        ax.set_yticks([])\n",
    "        ax.set_yticks([], minor=True)\n",
    "        ax.spines[\"left\"].set_visible(False)\n",
    "\n",
    "    axs[0].set_xticks([-0.05, 0, +0.05], [\"-50\", \"0\", \"+50\"], fontsize=8)\n",
    "    axs[0].set_xlabel(\"Lags [ms]\", fontsize=8, labelpad=0)\n",
    "    for ax in axs[1:]:\n",
    "        ax.set_xticks([-0.05, 0, +0.05], [])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b657c060",
   "metadata": {},
   "source": [
    "## CCG and power spectra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "fe42f26d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def power_spectrum(signal, samplerate):\n",
    "    n = len(signal)\n",
    "    AudioFreq = fft(signal)\n",
    "    AudioFreq = AudioFreq[0 : int(np.ceil((n + 1) / 2.0))]  # Half of the spectrum\n",
    "    MagFreq = np.abs(AudioFreq)  # Magnitude\n",
    "    MagFreq = MagFreq / float(n)\n",
    "    # power spectrum\n",
    "    MagFreq = MagFreq**2\n",
    "    if n % 2 > 0:  # ffte odd\n",
    "        MagFreq[1 : len(MagFreq)] = MagFreq[1 : len(MagFreq)] * 2\n",
    "    else:  # fft even\n",
    "        MagFreq[1 : len(MagFreq) - 1] = MagFreq[1 : len(MagFreq) - 1] * 2\n",
    "\n",
    "    freqAxis = np.arange(0, int(np.ceil((n + 1) / 2.0)), 1.0) * (samplerate / n)\n",
    "    return freqAxis, MagFreq\n",
    "\n",
    "\n",
    "def figure_xcorr_pspectra(\n",
    "    twostim_ccg: pd.DataFrame, single_ccg: pd.DataFrame, axs=None, colors=c_forebrain\n",
    "):\n",
    "    \"\"\"Boxplot of Change in Spike Rate as function of Relative Level\"\"\"\n",
    "\n",
    "    if axs is None:\n",
    "        fig, axs = plt.subplots(1, 6)\n",
    "    axs = axs.flatten()\n",
    "\n",
    "    merged_ccg = twostim_ccg.join(single_ccg, how=\"inner\", rsuffix=\"_single\")\n",
    "    merged_ccg = merged_ccg.loc[\n",
    "        (merged_ccg[\"fixedintensity\"] == merged_ccg[\"intensity\"])\n",
    "    ]\n",
    "    merged_ccg = merged_ccg.loc[(merged_ccg[\"fixedazi\"] == merged_ccg[\"azimuth\"])]\n",
    "    merged_ccg = merged_ccg.loc[(merged_ccg[\"fixedele\"] == merged_ccg[\"elevation\"])]\n",
    "\n",
    "    print(\"Number of unit pairs:\", np.unique(merged_ccg.index).size)\n",
    "    print(\n",
    "        \"Number of sessions:\",\n",
    "        len(sorted(set([idx[:2] for idx in merged_ccg.index.values]))),\n",
    "    )\n",
    "\n",
    "    # Before reindexing by relative_level:\n",
    "\n",
    "    merged_ccg.set_index(\"relative_level\", inplace=True)\n",
    "\n",
    "    psth_len = (merged_ccg.iloc[0][\"ccg_single\"].size + 1) / 2\n",
    "    lags = scipy.signal.correlation_lags(psth_len, psth_len) / 1000\n",
    "    lags_mask = np.abs(lags) <= 0.05\n",
    "\n",
    "    relative_levels = np.unique(merged_ccg.index)\n",
    "\n",
    "    for k, relative_level in enumerate(relative_levels):\n",
    "        mean_ccg = np.mean(np.vstack(merged_ccg.loc[relative_level, \"ccg\"]), axis=0)\n",
    "        n_ccg = merged_ccg.loc[relative_level, \"xcorr_peak_single\"].count()\n",
    "\n",
    "        avg_psd = []\n",
    "        for trace in np.vstack(merged_ccg.loc[relative_level, \"ccg\"]):\n",
    "            freq, psd = power_spectrum(\n",
    "                trace[(lags > -0.200) & (lags < 0.200)]\n",
    "                - np.mean(trace[(lags > -0.200) & (lags < 0.200)]),\n",
    "                1000,\n",
    "            )\n",
    "            avg_psd.append(psd)\n",
    "\n",
    "        norm_psd = [trace - np.mean(trace) for trace in avg_psd]\n",
    "        avg_psd = np.mean(norm_psd, axis=0)\n",
    "\n",
    "        ind_correct = np.searchsorted(freq, 50)\n",
    "        ind_max = np.argmax(avg_psd[(freq > 50) & (freq <= 80)])\n",
    "        freq_max = freq[ind_max + ind_correct]\n",
    "\n",
    "        axs[k].plot(freq, avg_psd, color=\"k\", ls=\"-\", lw=1)\n",
    "\n",
    "        if relative_level >= 0:\n",
    "            axs[k].plot(\n",
    "                freq_max, avg_psd[ind_max] + 2.5e-13, \"v\", color=colors[0], markersize=4\n",
    "            )\n",
    "\n",
    "        axs[k].spines[\"top\"].set_visible(False)\n",
    "        axs[k].spines[\"right\"].set_visible(False)\n",
    "\n",
    "        axs[k].set_ylim(bottom=0, top=7.5e-13)\n",
    "        axs[k].set_xlim(left=20, right=100)\n",
    "\n",
    "    axs[0].set_ylabel(\"Power\", labelpad=12)\n",
    "\n",
    "    def formatter_y(x, pos):\n",
    "        if x > 0:\n",
    "            return f\"${{{x*1e13:.0f}}}$\"\n",
    "            # return f\"$10^{{{np.log10(x):.0f}}}$\"\n",
    "        else:\n",
    "            return \"0\"\n",
    "\n",
    "    axs[0].set_yticks([0, 6.5e-13])\n",
    "    axs[0].yaxis.set_major_formatter(mpl.ticker.FuncFormatter(formatter_y))\n",
    "    axs[0].set_yticks(np.array([0, 3.5]) * 1e-13, minor=True)\n",
    "    axs[0].text(\n",
    "        0,\n",
    "        1,\n",
    "        \"$10^{-13}$\",\n",
    "        ha=\"right\",\n",
    "        va=\"center\",\n",
    "        fontsize=8,\n",
    "        transform=axs[0].transAxes\n",
    "        + mpl.transforms.ScaledTranslation(\n",
    "            -2 / 72,\n",
    "            0,\n",
    "            axs[0].figure.dpi_scale_trans,\n",
    "        ),\n",
    "    )\n",
    "\n",
    "    for ax in axs[1:]:\n",
    "        ax.set_yticks([])\n",
    "        ax.set_yticks([], minor=True)\n",
    "        ax.spines[\"left\"].set_visible(False)\n",
    "\n",
    "    axs[0].set_xticks([50, 75], [\"50\", \"75\"], fontsize=8, minor=False)\n",
    "    axs[0].set_xticks([50, 75], [], fontsize=8, minor=True)\n",
    "\n",
    "    axs[0].set_xlabel(\"Freq. [Hz]\", fontsize=8, labelpad=0)\n",
    "    for ax in axs[1:]:\n",
    "        ax.set_xticks([50, 75], [\"50\", \"75\"], fontsize=8, minor=False)\n",
    "\n",
    "        ax.set_xticks([50, 75], [\"50\", \"75\"], minor=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9391fb76",
   "metadata": {},
   "source": [
    "## CSI boxplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "524197a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def figure_xcorr_boxplot(\n",
    "    twostim_rlf: pd.DataFrame,\n",
    "    single_rlf: pd.DataFrame,\n",
    "    ax=None,\n",
    "    colors=c_forebrain,\n",
    "    brackets: dict[tuple[int, int] : float] | None = None,\n",
    "    anova_align=\"right\",\n",
    "):\n",
    "    \"\"\"Boxplot of Change in Spike Rate as function of Relative Level\"\"\"\n",
    "\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots()\n",
    "\n",
    "    merged_rlf = twostim_rlf.join(single_rlf, how=\"inner\", rsuffix=\"_single\")\n",
    "    merged_rlf = merged_rlf.loc[\n",
    "        (merged_rlf[\"fixedintensity\"] == merged_rlf[\"intensity\"])\n",
    "    ]\n",
    "    merged_rlf = merged_rlf.loc[(merged_rlf[\"fixedazi\"] == merged_rlf[\"azimuth\"])]\n",
    "    merged_rlf = merged_rlf.loc[(merged_rlf[\"fixedele\"] == merged_rlf[\"elevation\"])]\n",
    "\n",
    "    print(\"Number of units:\", np.unique(merged_rlf.index).size)\n",
    "    print(\n",
    "        \"Number of sessions:\",\n",
    "        len(sorted(set([idx[:2] for idx in merged_rlf.index.values]))),\n",
    "    )\n",
    "    merged_rlf.set_index(\"relative_level\", inplace=True)\n",
    "    # display(merged_rlf_ot)\n",
    "\n",
    "    merged_rlf[\"csi\"] = (merged_rlf[\"xcorr_peak\"] - merged_rlf[\"xcorr_peak_single\"]) / (\n",
    "        merged_rlf[\"xcorr_peak\"] + merged_rlf[\"xcorr_peak_single\"]\n",
    "    )\n",
    "    relative_levels = np.unique(merged_rlf.index)\n",
    "    groupdata = [\n",
    "        merged_rlf.loc[relative_level, \"csi\"].values\n",
    "        for relative_level in relative_levels\n",
    "    ]\n",
    "\n",
    "    bp = ax.boxplot(\n",
    "        groupdata,\n",
    "        patch_artist=True,\n",
    "        notch=False,\n",
    "        showfliers=False,\n",
    "        whis=(5, 95),\n",
    "    )\n",
    "\n",
    "    ax.set_xticklabels(relative_levels.astype(int))\n",
    "\n",
    "    stats = anova_tukey(merged_rlf, val_col=\"csi\", group_col=\"relative_level\")\n",
    "\n",
    "    for patch, color in zip(bp[\"boxes\"], colors):\n",
    "        patch.set_facecolor(color)\n",
    "\n",
    "    for median in bp[\"medians\"]:\n",
    "        median.set(color=\"white\", linewidth=2)\n",
    "\n",
    "    if True:  # show_individual_data\n",
    "        width = 0.7\n",
    "        for k, groupvalues in enumerate(groupdata):\n",
    "            pos = k + 1\n",
    "            y = np.sort(groupvalues)\n",
    "            vals_x = np.linspace(\n",
    "                pos - width / 3,\n",
    "                pos + width / 3,\n",
    "                y.size,\n",
    "                endpoint=True,\n",
    "            )\n",
    "            ax.plot(\n",
    "                vals_x,\n",
    "                y,\n",
    "                color=\"k\",\n",
    "                ls=\"None\",\n",
    "                marker=\".\",\n",
    "                markersize=2,\n",
    "                markeredgewidth=0.0,\n",
    "                zorder=5,\n",
    "            )\n",
    "\n",
    "    if True:  # show_brackets\n",
    "        brackets_shrink = 1.0\n",
    "        if brackets is None:\n",
    "            brackets = {(k, k + 1): 0.8 for k in range(len(stats[\"tukey\"]) - 1)}\n",
    "            brackets_shrink = 0.8\n",
    "        for (idx_left, idx_right), y in brackets.items():\n",
    "            plot_bracket(\n",
    "                ax,\n",
    "                left=idx_left + 1,\n",
    "                right=idx_right + 1,\n",
    "                text=f\"{stats['tukey'].iloc[idx_left, idx_right]:.3f}\"[1:],\n",
    "                y=y,\n",
    "                shrink=brackets_shrink,\n",
    "            )\n",
    "\n",
    "    ax.spines[\"top\"].set_visible(False)\n",
    "    ax.spines[\"right\"].set_visible(False)\n",
    "\n",
    "    if anova_align == \"right\":\n",
    "        ax.text(\n",
    "            1,\n",
    "            1,\n",
    "            f\"p = {stats['anova_p']:.3g}\",\n",
    "            horizontalalignment=\"right\",\n",
    "            verticalalignment=\"top\",\n",
    "            fontsize=10,\n",
    "            transform=ax.transAxes,\n",
    "        )\n",
    "    else:\n",
    "        ax.text(\n",
    "            0,\n",
    "            1,\n",
    "            f\"p = {stats['anova_p']:.3g}\",\n",
    "            horizontalalignment=\"left\",\n",
    "            verticalalignment=\"top\",\n",
    "            fontsize=10,\n",
    "            transform=ax.transAxes\n",
    "            + mpl.transforms.ScaledTranslation(\n",
    "                +2 / 72,\n",
    "                0 / 72,\n",
    "                ax.figure.dpi_scale_trans,\n",
    "            ),\n",
    "        )\n",
    "\n",
    "    # ax.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, pos: f\"{x:.0%}\"))\n",
    "\n",
    "    ax.axhline(0, color=\"k\", ls=\"--\", zorder=-1)\n",
    "    ax.set_xlabel(\"Relative Level [dB]\")\n",
    "    ax.set_ylabel(\"CSI\")\n",
    "    ax.set_ylim(-1, 1.45)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4e957b4",
   "metadata": {},
   "source": [
    "## histogram of tuning selectivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "1968ed60",
   "metadata": {},
   "outputs": [],
   "source": [
    "def histogram_selectivity_resubmit(\n",
    "    twostim_rlf: pd.DataFrame,\n",
    "    single_rlf: pd.DataFrame,\n",
    "    single_rlf_out: pd.DataFrame,\n",
    "    ax=None,\n",
    "    colors=c_forebrain,\n",
    "):\n",
    "    \"\"\"Boxplot of Change in Spike Rate as function of Relative Level\"\"\"\n",
    "\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots(1, 1)\n",
    "\n",
    "    ## merge single rlf data and calculate selectivity index\n",
    "    single_in = single_rlf.copy()\n",
    "    single_out = single_rlf_out.copy()\n",
    "\n",
    "    frontal_driver = single_in[abs(single_in[\"azimuth\"]) <= 45]\n",
    "    peripheral_driver = single_in[abs(single_in[\"azimuth\"]) > 45]\n",
    "\n",
    "    frontal_competitor = single_out[abs(single_out[\"azimuth\"]) <= 45]\n",
    "    peripheral_competitor = single_out[abs(single_out[\"azimuth\"]) > 45]\n",
    "\n",
    "    merged_rlf1 = frontal_driver.join(\n",
    "        peripheral_competitor, how=\"inner\", rsuffix=\"_out\"\n",
    "    )\n",
    "    merged_rlf1 = merged_rlf1.loc[\n",
    "        merged_rlf1[\"intensity\"] == merged_rlf1[\"intensity_out\"]\n",
    "    ]\n",
    "    merged_rlf1[\"selectivity_index\"] = (\n",
    "        merged_rlf1[\"resp\"] - merged_rlf1[\"resp_out\"]\n",
    "    ) / (merged_rlf1[\"resp\"] + merged_rlf1[\"resp_out\"])\n",
    "\n",
    "    merged_rlf2 = frontal_competitor.join(\n",
    "        peripheral_driver, how=\"inner\", rsuffix=\"_out\"\n",
    "    )\n",
    "    merged_rlf2 = merged_rlf2.loc[\n",
    "        merged_rlf2[\"intensity\"] == merged_rlf2[\"intensity_out\"]\n",
    "    ]\n",
    "    merged_rlf2[\"selectivity_index\"] = (\n",
    "        merged_rlf2[\"resp\"] - merged_rlf2[\"resp_out\"]\n",
    "    ) / (merged_rlf2[\"resp\"] + merged_rlf2[\"resp_out\"])\n",
    "\n",
    "    tuning_index_1 = merged_rlf1.filter(\n",
    "        [\n",
    "            \"azimuth\",\n",
    "            \"elevation\",\n",
    "            \"azimuth_out\",\n",
    "            \"elevation_out\",\n",
    "            \"intensity\",\n",
    "            \"resp\",\n",
    "            \"selectivity_index\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "    tuning_index_2 = merged_rlf2.filter(\n",
    "        [\n",
    "            \"azimuth\",\n",
    "            \"elevation\",\n",
    "            \"azimuth_out\",\n",
    "            \"elevation_out\",\n",
    "            \"intensity\",\n",
    "            \"resp\",\n",
    "            \"selectivity_index\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "    tuning_index = pd.concat([tuning_index_1, tuning_index_2])\n",
    "\n",
    "    merged_1 = twostim_rlf.join(tuning_index_1, how=\"inner\", rsuffix=\"_single\")\n",
    "    merged_1 = merged_1.loc[(merged_1[\"fixedintensity\"] == merged_1[\"intensity\"])]\n",
    "    merged_1 = merged_1.loc[(merged_1[\"fixedazi\"] == merged_1[\"azimuth\"])]\n",
    "    merged_1 = merged_1.loc[(merged_1[\"fixedele\"] == merged_1[\"elevation\"])]\n",
    "\n",
    "    merged_1 = merged_1.loc[(merged_1[\"varyingazi\"] == merged_1[\"azimuth_out\"])]\n",
    "    merged_1 = merged_1.loc[(merged_1[\"varyingele\"] == merged_1[\"elevation_out\"])]\n",
    "\n",
    "    merged_2 = twostim_rlf.join(tuning_index_2, how=\"inner\", rsuffix=\"_single\")\n",
    "    merged_2 = merged_2.loc[(merged_2[\"fixedintensity\"] == merged_2[\"intensity\"])]\n",
    "\n",
    "    merged_2 = merged_2.loc[(merged_2[\"fixedazi\"] == merged_2[\"azimuth_out\"])]\n",
    "    merged_2 = merged_2.loc[(merged_2[\"fixedele\"] == merged_2[\"elevation_out\"])]\n",
    "\n",
    "    merged_2 = merged_2.loc[(merged_2[\"varyingazi\"] == merged_2[\"azimuth\"])]\n",
    "    merged_2 = merged_2.loc[(merged_2[\"varyingele\"] == merged_2[\"elevation\"])]\n",
    "\n",
    "    merged_rlf = pd.concat([merged_1, merged_2])\n",
    "\n",
    "    merged_rlf = merged_rlf[merged_rlf[\"relative_level\"] == 0]\n",
    "\n",
    "    quantiles = np.quantile(merged_rlf[\"selectivity_index\"], [0.25, 0.5, 0.75])\n",
    "    # print(quantiles[1])\n",
    "\n",
    "    print(\"Number of units:\", np.unique(merged_rlf.index).size)\n",
    "    print(\n",
    "        \"Number of frontal driver units:\",\n",
    "        np.unique(\n",
    "            merged_rlf[merged_rlf[\"selectivity_index\"] >= quantiles[2]].index\n",
    "        ).size,\n",
    "    )\n",
    "    print(\n",
    "        \"Number of peripheral driver units:\",\n",
    "        np.unique(\n",
    "            merged_rlf[merged_rlf[\"selectivity_index\"] <= quantiles[0]].index\n",
    "        ).size,\n",
    "    )\n",
    "    print(\n",
    "        \"Number of intermediate units:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (merged_rlf[\"selectivity_index\"] > quantiles[0])\n",
    "                & (merged_rlf[\"selectivity_index\"] < quantiles[2])\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "\n",
    "    print(\n",
    "        \"Number of sessions:\",\n",
    "        len(sorted(set([idx[:2] for idx in merged_rlf.index.values]))),\n",
    "    )\n",
    "\n",
    "    tuning_bins = np.linspace(-1, 1, 100)\n",
    "    tuning_index_plot = ax.hist(\n",
    "        merged_rlf[\"selectivity_index\"], tuning_bins, color=colors[0]\n",
    "    )\n",
    "    mean_SI = np.mean(merged_rlf[\"selectivity_index\"])\n",
    "    std_SI = np.std(merged_rlf[\"selectivity_index\"])\n",
    "\n",
    "    print(f\"MEAN:{mean_SI:.2f}\")\n",
    "    print(f\"STD: {std_SI:.2f}\")\n",
    "\n",
    "    ax.axvline(0, color=\"k\", ls=\":\")\n",
    "\n",
    "    ax.spines[\"top\"].set_visible(False)\n",
    "    ax.spines[\"right\"].set_visible(False)\n",
    "    ax.legend(loc=\"upper left\", frameon=False)\n",
    "\n",
    "    ax.set_xlabel(\"Tuning Selectivity Index\")\n",
    "    ax.set_ylabel(\"# Units\")\n",
    "    ax.set_ylim(0, 20)\n",
    "\n",
    "    ax.text(\n",
    "        0.8,\n",
    "        1,\n",
    "        \"Frontal Driver Selective\",\n",
    "        horizontalalignment=\"right\",\n",
    "        verticalalignment=\"top\",\n",
    "        fontsize=10,\n",
    "        transform=ax.transAxes,\n",
    "    )\n",
    "    ax.text(\n",
    "        0.45,\n",
    "        1,\n",
    "        \"Lateral Driver Selective\",\n",
    "        horizontalalignment=\"right\",\n",
    "        verticalalignment=\"top\",\n",
    "        fontsize=10,\n",
    "        transform=ax.transAxes,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3deecbc",
   "metadata": {},
   "source": [
    "## competition by selectivity index line plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "7a6c84a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def figure_rlf_selectivity_resubmit(\n",
    "    twostim_rlf: pd.DataFrame,\n",
    "    single_rlf: pd.DataFrame,\n",
    "    single_rlf_out: pd.DataFrame,\n",
    "    ax=None,\n",
    "    colors=c_forebrain,\n",
    "    brackets: dict[tuple[int, int] : float] | None = None,\n",
    "    anova_align=\"right\",\n",
    "):\n",
    "    \"\"\"Boxplot of Change in Spike Rate as function of Relative Level\"\"\"\n",
    "\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots()\n",
    "\n",
    "    ## merge single rlf data and calculate selectivity index\n",
    "    single_in = single_rlf.copy()\n",
    "    single_out = single_rlf_out.copy()\n",
    "\n",
    "    frontal_driver = single_in[abs(single_in[\"azimuth\"]) <= 45]\n",
    "    peripheral_driver = single_in[abs(single_in[\"azimuth\"]) > 45]\n",
    "\n",
    "    frontal_competitor = single_out[abs(single_out[\"azimuth\"]) <= 45]\n",
    "    peripheral_competitor = single_out[abs(single_out[\"azimuth\"]) > 45]\n",
    "\n",
    "    merged_rlf1 = frontal_driver.join(\n",
    "        peripheral_competitor, how=\"inner\", rsuffix=\"_out\"\n",
    "    )\n",
    "    merged_rlf1 = merged_rlf1.loc[\n",
    "        merged_rlf1[\"intensity\"] == merged_rlf1[\"intensity_out\"]\n",
    "    ]\n",
    "    merged_rlf1[\"selectivity_index\"] = (\n",
    "        merged_rlf1[\"resp\"] - merged_rlf1[\"resp_out\"]\n",
    "    ) / (merged_rlf1[\"resp\"] + merged_rlf1[\"resp_out\"])\n",
    "\n",
    "    merged_rlf2 = frontal_competitor.join(\n",
    "        peripheral_driver, how=\"inner\", rsuffix=\"_out\"\n",
    "    )\n",
    "    merged_rlf2 = merged_rlf2.loc[\n",
    "        merged_rlf2[\"intensity\"] == merged_rlf2[\"intensity_out\"]\n",
    "    ]\n",
    "    merged_rlf2[\"selectivity_index\"] = (\n",
    "        merged_rlf2[\"resp\"] - merged_rlf2[\"resp_out\"]\n",
    "    ) / (merged_rlf2[\"resp\"] + merged_rlf2[\"resp_out\"])\n",
    "\n",
    "    tuning_index_1 = merged_rlf1.filter(\n",
    "        [\n",
    "            \"azimuth\",\n",
    "            \"elevation\",\n",
    "            \"azimuth_out\",\n",
    "            \"elevation_out\",\n",
    "            \"intensity\",\n",
    "            \"resp\",\n",
    "            \"selectivity_index\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "    tuning_index_2 = merged_rlf2.filter(\n",
    "        [\n",
    "            \"azimuth\",\n",
    "            \"elevation\",\n",
    "            \"azimuth_out\",\n",
    "            \"elevation_out\",\n",
    "            \"intensity\",\n",
    "            \"resp_out\",\n",
    "            \"selectivity_index\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "    tuning_index_2 = tuning_index_2.rename(columns={\"resp_out\": \"resp\"})\n",
    "\n",
    "    merged_1 = twostim_rlf.join(tuning_index_1, how=\"inner\", rsuffix=\"_single\")\n",
    "    merged_1 = merged_1.loc[(merged_1[\"fixedintensity\"] == merged_1[\"intensity\"])]\n",
    "    merged_1 = merged_1.loc[(merged_1[\"fixedazi\"] == merged_1[\"azimuth\"])]\n",
    "    merged_1 = merged_1.loc[(merged_1[\"fixedele\"] == merged_1[\"elevation\"])]\n",
    "\n",
    "    merged_1 = merged_1.loc[(merged_1[\"varyingazi\"] == merged_1[\"azimuth_out\"])]\n",
    "    merged_1 = merged_1.loc[(merged_1[\"varyingele\"] == merged_1[\"elevation_out\"])]\n",
    "\n",
    "    merged_2 = twostim_rlf.join(tuning_index_2, how=\"inner\", rsuffix=\"_single\")\n",
    "    merged_2 = merged_2.loc[(merged_2[\"fixedintensity\"] == merged_2[\"intensity\"])]\n",
    "\n",
    "    merged_2 = merged_2.loc[(merged_2[\"fixedazi\"] == merged_2[\"azimuth_out\"])]\n",
    "    merged_2 = merged_2.loc[(merged_2[\"fixedele\"] == merged_2[\"elevation_out\"])]\n",
    "    merged_2 = merged_2.loc[(merged_2[\"varyingazi\"] == merged_2[\"azimuth\"])]\n",
    "    merged_2 = merged_2.loc[(merged_2[\"varyingele\"] == merged_2[\"elevation\"])]\n",
    "    merged_rlf = pd.concat([merged_1, merged_2])\n",
    "\n",
    "    quantiles = np.quantile(merged_rlf[\"selectivity_index\"], [0.25, 0.5, 0.75])\n",
    "\n",
    "    print(\n",
    "        \"Number of frontal driver units + frontal driver:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (abs(merged_rlf[\"fixedazi\"]) <= 45)\n",
    "                & (abs(merged_rlf[\"varyingazi\"]) > 45)\n",
    "                & (merged_rlf[\"selectivity_index\"] > 0)\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "    print(\n",
    "        \"Number of frontal driver units + lateral driver:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (abs(merged_rlf[\"fixedazi\"]) > 45)\n",
    "                & (abs(merged_rlf[\"varyingazi\"]) <= 45)\n",
    "                & (merged_rlf[\"selectivity_index\"] > 0)\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "\n",
    "    ##lateral driver selective units\n",
    "    print(\n",
    "        \"Number of lateral driver units + frontal driver:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (abs(merged_rlf[\"fixedazi\"]) <= 45)\n",
    "                & (abs(merged_rlf[\"varyingazi\"]) > 45)\n",
    "                & (merged_rlf[\"selectivity_index\"] <= 0)\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "    print(\n",
    "        \"Number of lateral driver units + lateral driver:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (abs(merged_rlf[\"fixedazi\"]) > 45)\n",
    "                & (abs(merged_rlf[\"varyingazi\"]) <= 45)\n",
    "                & (merged_rlf[\"selectivity_index\"] <= 0)\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "\n",
    "    merged_rlf[\"change_in_response\"] = (\n",
    "        merged_rlf[\"resp\"] - merged_rlf[\"resp_single\"]\n",
    "    ) / merged_rlf[\"resp_single\"]\n",
    "\n",
    "    merged_rlf.set_index(\"relative_level\", inplace=True)\n",
    "\n",
    "    ##sort data based on frontal driver + lateral competitor\n",
    "    frontal_driver = merged_rlf[\n",
    "        (abs(merged_rlf[\"fixedazi\"]) < 45)\n",
    "        & (abs(merged_rlf[\"varyingazi\"]) >= 45)\n",
    "        & (merged_rlf[\"selectivity_index\"] > 0)\n",
    "    ]\n",
    "\n",
    "    peripheral_driver = merged_rlf[\n",
    "        (abs(merged_rlf[\"fixedazi\"]) > 45)\n",
    "        & (abs(merged_rlf[\"varyingazi\"]) <= 45)\n",
    "        & (merged_rlf[\"selectivity_index\"] > 0)\n",
    "    ]\n",
    "\n",
    "    stats_frontal_driver = anova_tukey(\n",
    "        frontal_driver, val_col=\"change_in_response\", group_col=\"relative_level\"\n",
    "    )\n",
    "    stats_peripheral_driver = anova_tukey(\n",
    "        peripheral_driver, val_col=\"change_in_response\", group_col=\"relative_level\"\n",
    "    )\n",
    "\n",
    "    frontal_driver_means = frontal_driver.groupby(\"relative_level\")[\n",
    "        \"change_in_response\"\n",
    "    ].mean()\n",
    "    peripheral_driver_means = peripheral_driver.groupby(\"relative_level\")[\n",
    "        \"change_in_response\"\n",
    "    ].mean()\n",
    "\n",
    "    frontal_driver_sems = frontal_driver.groupby(\"relative_level\")[\n",
    "        \"change_in_response\"\n",
    "    ].sem()\n",
    "    peripheral_driver_sems = peripheral_driver.groupby(\"relative_level\")[\n",
    "        \"change_in_response\"\n",
    "    ].sem()\n",
    "\n",
    "    relative_levels = np.unique(frontal_driver_means.index)\n",
    "\n",
    "    frontal_driver_plot = ax.errorbar(\n",
    "        relative_levels,\n",
    "        frontal_driver_means,\n",
    "        frontal_driver_sems,\n",
    "        color=colors[0],\n",
    "        label=f\"Frontal Driver, p = {stats_frontal_driver['anova_p']:.3g}\",\n",
    "    )\n",
    "    peripheral_driver_plot = ax.errorbar(\n",
    "        relative_levels,\n",
    "        peripheral_driver_means,\n",
    "        peripheral_driver_sems,\n",
    "        color=\"b\",\n",
    "        label=f\"Lateral Driver, p = {stats_peripheral_driver['anova_p']:.3g}\",\n",
    "    )\n",
    "\n",
    "    ax.set_xticks(relative_levels.astype(int), relative_levels.astype(int))\n",
    "\n",
    "    ax.spines[\"top\"].set_visible(False)\n",
    "    ax.spines[\"right\"].set_visible(False)\n",
    "    ax.legend(loc=\"lower left\", frameon=False, fontsize=9)\n",
    "\n",
    "    ax.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, pos: f\"{x:.0%}\"))\n",
    "\n",
    "    ax.axhline(0, color=\"k\", ls=\"--\", zorder=-1)\n",
    "    ax.set_xlabel(\"Relative Level [dB]\")\n",
    "    ax.set_ylabel(\"Change in\\nSpike Rate\")\n",
    "    ax.set_ylim(-0.5, 0.95)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3604ea7a",
   "metadata": {},
   "source": [
    "## Single sound rate level function by selectivity index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "6dbd506b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def figure_singlesound_rlf_selectivity_resubmit(\n",
    "    single_rlf: pd.DataFrame,\n",
    "    single_rlf_out: pd.DataFrame,\n",
    "    ax=None,\n",
    "    colors=c_forebrain,\n",
    "    brackets: dict[tuple[int, int] : float] | None = None,\n",
    "    anova_align=\"right\",\n",
    "):\n",
    "    \"\"\"Boxplot of Change in Spike Rate as function of Relative Level\"\"\"\n",
    "\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots(1, 1, sharey=True)\n",
    "\n",
    "    INTENSITY_0_EQUALS_DB_SPL = +63.0\n",
    "\n",
    "    single_in = single_rlf.copy()\n",
    "    single_out = single_rlf_out.copy()\n",
    "\n",
    "    frontal_driver = single_in[abs(single_in[\"azimuth\"]) <= 45]\n",
    "    peripheral_driver = single_in[abs(single_in[\"azimuth\"]) > 45]\n",
    "\n",
    "    frontal_competitor = single_out[abs(single_out[\"azimuth\"]) <= 45]\n",
    "    peripheral_competitor = single_out[abs(single_out[\"azimuth\"]) > 45]\n",
    "\n",
    "    merged_rlf1 = frontal_driver.join(\n",
    "        peripheral_competitor, how=\"inner\", rsuffix=\"_out\"\n",
    "    )\n",
    "    merged_rlf1 = merged_rlf1.loc[\n",
    "        merged_rlf1[\"intensity\"] == merged_rlf1[\"intensity_out\"]\n",
    "    ]\n",
    "    merged_rlf1[\"selectivity_index\"] = (\n",
    "        merged_rlf1[(merged_rlf1[\"intensity\"] == -20)][\"resp\"]\n",
    "        - merged_rlf1[(merged_rlf1[\"intensity\"] == -20)][\"resp_out\"]\n",
    "    ) / (\n",
    "        merged_rlf1[(merged_rlf1[\"intensity\"] == -20)][\"resp\"]\n",
    "        + merged_rlf1[(merged_rlf1[\"intensity\"] == -20)][\"resp_out\"]\n",
    "    )\n",
    "\n",
    "    merged_rlf2 = frontal_competitor.join(\n",
    "        peripheral_driver, how=\"inner\", rsuffix=\"_out\"\n",
    "    )\n",
    "    merged_rlf2 = merged_rlf2.loc[\n",
    "        merged_rlf2[\"intensity\"] == merged_rlf2[\"intensity_out\"]\n",
    "    ]\n",
    "    merged_rlf2[\"selectivity_index\"] = (\n",
    "        merged_rlf2[(merged_rlf2[\"intensity\"] == -20)][\"resp\"]\n",
    "        - merged_rlf2[(merged_rlf2[\"intensity\"] == -20)][\"resp_out\"]\n",
    "    ) / (\n",
    "        merged_rlf2[(merged_rlf2[\"intensity\"] == -20)][\"resp\"]\n",
    "        + merged_rlf2[(merged_rlf2[\"intensity\"] == -20)][\"resp_out\"]\n",
    "    )\n",
    "\n",
    "    tuning_index_1 = merged_rlf1.filter(\n",
    "        [\n",
    "            \"azimuth\",\n",
    "            \"elevation\",\n",
    "            \"azimuth_out\",\n",
    "            \"elevation_out\",\n",
    "            \"intensity\",\n",
    "            \"norm_resp\",\n",
    "            \"norm_resp_out\",\n",
    "            \"selectivity_index\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "    tuning_index_2 = merged_rlf2.filter(\n",
    "        [\n",
    "            \"azimuth\",\n",
    "            \"elevation\",\n",
    "            \"azimuth_out\",\n",
    "            \"elevation_out\",\n",
    "            \"intensity\",\n",
    "            \"norm_resp\",\n",
    "            \"norm_resp_out\",\n",
    "            \"selectivity_index\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "    merged_rlf = pd.concat([tuning_index_1, tuning_index_2])\n",
    "\n",
    "    quantiles = np.quantile(merged_rlf[\"selectivity_index\"], [0.25, 0.5, 0.75])\n",
    "\n",
    "    competitor_selective = merged_rlf[merged_rlf[\"selectivity_index\"] <= 0]\n",
    "\n",
    "    driver_selective = merged_rlf[(merged_rlf[\"selectivity_index\"] > 0)]\n",
    "\n",
    "    print(\"Number of units:\", np.unique(merged_rlf.index).size)\n",
    "    print(\n",
    "        \"Number of driver units:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (merged_rlf[\"selectivity_index\"] > 0) & (merged_rlf[\"intensity\"] == -20)\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "    print(\n",
    "        \"Number of competitor units:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (merged_rlf[\"selectivity_index\"] <= 0)\n",
    "                & (merged_rlf[\"intensity\"] == -20)\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "\n",
    "    print(\n",
    "        \"Number of sessions:\",\n",
    "        len(sorted(set([idx[:2] for idx in merged_rlf.index.values]))),\n",
    "    )\n",
    "    merged_rlf = merged_rlf.set_index(\"intensity\")\n",
    "\n",
    "    intensities = np.unique(merged_rlf.index) + INTENSITY_0_EQUALS_DB_SPL\n",
    "\n",
    "    DRIVERALONE_driverselective_means = driver_selective.groupby(\"intensity\")[\n",
    "        \"norm_resp\"\n",
    "    ].mean()\n",
    "    DRIVERALONE_driverselective_sems = driver_selective.groupby(\"intensity\")[\n",
    "        \"norm_resp\"\n",
    "    ].sem()\n",
    "\n",
    "    COMPETITORONE_driverselective_means = driver_selective.groupby(\"intensity\")[\n",
    "        \"norm_resp_out\"\n",
    "    ].mean()\n",
    "    COMPETITORALONE_driverselective_sems = driver_selective.groupby(\"intensity\")[\n",
    "        \"norm_resp_out\"\n",
    "    ].sem()\n",
    "\n",
    "    display(len(np.unique(driver_selective.index)))\n",
    "\n",
    "    print(np.asarray(DRIVERALONE_driverselective_means.values))\n",
    "    print(COMPETITORONE_driverselective_means.values)\n",
    "\n",
    "    ax.set_title(\"Frontal Driver Selective\")\n",
    "    ax.errorbar(\n",
    "        intensities,\n",
    "        DRIVERALONE_driverselective_means.values,\n",
    "        DRIVERALONE_driverselective_sems.values,\n",
    "        color=colors[0],\n",
    "        label=\"Frontal stimulus alone\",\n",
    "    )\n",
    "    ax.errorbar(\n",
    "        intensities,\n",
    "        COMPETITORONE_driverselective_means.values,\n",
    "        COMPETITORALONE_driverselective_sems.values,\n",
    "        color=\"b\",\n",
    "        label=\"Lateral stimulus alone\",\n",
    "    )\n",
    "\n",
    "    ax.axvspan(\n",
    "        -35 + INTENSITY_0_EQUALS_DB_SPL,\n",
    "        -10 + INTENSITY_0_EQUALS_DB_SPL,\n",
    "        color=\"pink\",\n",
    "        alpha=0.5,\n",
    "        zorder=-1,\n",
    "    )\n",
    "\n",
    "    ax.set_xlabel(\"Average binaural level [dB]\")\n",
    "    ax.set_ylabel(\"Norm Resp.\")\n",
    "    ax.legend(loc=\"lower right\", frameon=False)\n",
    "\n",
    "    ax.set_ylim(bottom=0, top=0.80)\n",
    "    ax.spines[\"top\"].set_visible(False)\n",
    "    ax.spines[\"right\"].set_visible(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "218ad7cb",
   "metadata": {},
   "source": [
    "## Extreme selectivity competition "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "09d5713e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def extreme_selectivity_competition_resubmit(\n",
    "    twostim_rlf: pd.DataFrame,\n",
    "    single_rlf: pd.DataFrame,\n",
    "    single_rlf_out: pd.DataFrame,\n",
    "    axs=None,\n",
    "    colors=c_forebrain,\n",
    "):\n",
    "\n",
    "    if axs is None:\n",
    "        fig, axs = plt.subplots(1, 2)\n",
    "\n",
    "    ## merge single rlf data and calculate selectivity index\n",
    "    single_in = single_rlf.copy()\n",
    "    single_out = single_rlf_out.copy()\n",
    "\n",
    "    frontal_driver = single_in[abs(single_in[\"azimuth\"]) <= 45]\n",
    "    peripheral_driver = single_in[abs(single_in[\"azimuth\"]) > 45]\n",
    "\n",
    "    frontal_competitor = single_out[abs(single_out[\"azimuth\"]) <= 45]\n",
    "    peripheral_competitor = single_out[abs(single_out[\"azimuth\"]) > 45]\n",
    "\n",
    "    merged_rlf1 = frontal_driver.join(\n",
    "        peripheral_competitor, how=\"inner\", rsuffix=\"_out\"\n",
    "    )\n",
    "    merged_rlf1 = merged_rlf1.loc[\n",
    "        merged_rlf1[\"intensity\"] == merged_rlf1[\"intensity_out\"]\n",
    "    ]\n",
    "    merged_rlf1[\"selectivity_index\"] = (\n",
    "        merged_rlf1[\"resp\"] - merged_rlf1[\"resp_out\"]\n",
    "    ) / (merged_rlf1[\"resp\"] + merged_rlf1[\"resp_out\"])\n",
    "\n",
    "    merged_rlf2 = frontal_competitor.join(\n",
    "        peripheral_driver, how=\"inner\", rsuffix=\"_out\"\n",
    "    )\n",
    "    merged_rlf2 = merged_rlf2.loc[\n",
    "        merged_rlf2[\"intensity\"] == merged_rlf2[\"intensity_out\"]\n",
    "    ]\n",
    "    merged_rlf2[\"selectivity_index\"] = (\n",
    "        merged_rlf2[\"resp\"] - merged_rlf2[\"resp_out\"]\n",
    "    ) / (merged_rlf2[\"resp\"] + merged_rlf2[\"resp_out\"])\n",
    "\n",
    "    tuning_index_1 = merged_rlf1.filter(\n",
    "        [\n",
    "            \"azimuth\",\n",
    "            \"elevation\",\n",
    "            \"azimuth_out\",\n",
    "            \"elevation_out\",\n",
    "            \"intensity\",\n",
    "            \"resp\",\n",
    "            \"selectivity_index\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "    tuning_index_2 = merged_rlf2.filter(\n",
    "        [\n",
    "            \"azimuth\",\n",
    "            \"elevation\",\n",
    "            \"azimuth_out\",\n",
    "            \"elevation_out\",\n",
    "            \"intensity\",\n",
    "            \"resp_out\",\n",
    "            \"selectivity_index\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "    tuning_index_2 = tuning_index_2.rename(columns={\"resp_out\": \"resp\"})\n",
    "\n",
    "    merged_1 = twostim_rlf.join(tuning_index_1, how=\"inner\", rsuffix=\"_single\")\n",
    "    merged_1 = merged_1.loc[(merged_1[\"fixedintensity\"] == merged_1[\"intensity\"])]\n",
    "    merged_1 = merged_1.loc[(merged_1[\"fixedazi\"] == merged_1[\"azimuth\"])]\n",
    "    merged_1 = merged_1.loc[(merged_1[\"fixedele\"] == merged_1[\"elevation\"])]\n",
    "\n",
    "    merged_1 = merged_1.loc[(merged_1[\"varyingazi\"] == merged_1[\"azimuth_out\"])]\n",
    "    merged_1 = merged_1.loc[(merged_1[\"varyingele\"] == merged_1[\"elevation_out\"])]\n",
    "\n",
    "    merged_2 = twostim_rlf.join(tuning_index_2, how=\"inner\", rsuffix=\"_single\")\n",
    "    merged_2 = merged_2.loc[(merged_2[\"fixedintensity\"] == merged_2[\"intensity\"])]\n",
    "\n",
    "    merged_2 = merged_2.loc[(merged_2[\"fixedazi\"] == merged_2[\"azimuth_out\"])]\n",
    "    merged_2 = merged_2.loc[(merged_2[\"fixedele\"] == merged_2[\"elevation_out\"])]\n",
    "    merged_2 = merged_2.loc[(merged_2[\"varyingazi\"] == merged_2[\"azimuth\"])]\n",
    "    merged_2 = merged_2.loc[(merged_2[\"varyingele\"] == merged_2[\"elevation\"])]\n",
    "    merged_rlf = pd.concat([merged_1, merged_2])\n",
    "\n",
    "    quantiles = np.quantile(merged_rlf[\"selectivity_index\"], [0.25, 0.5, 0.75])\n",
    "\n",
    "    print(\n",
    "        \"Number of driver units:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (abs(merged_rlf[\"fixedazi\"]) <= 45)\n",
    "                & (merged_rlf[\"selectivity_index\"] > 0.30)\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "    print(\n",
    "        \"Number of competitor units:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (abs(merged_rlf[\"fixedazi\"]) <= 45)\n",
    "                & (merged_rlf[\"selectivity_index\"] <= -0.50)\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "\n",
    "    merged_rlf[\"change_in_response\"] = (\n",
    "        merged_rlf[\"resp\"] - merged_rlf[\"resp_single\"]\n",
    "    ) / merged_rlf[\"resp_single\"]\n",
    "\n",
    "    merged_rlf.set_index(\"relative_level\", inplace=True)\n",
    "\n",
    "    ## sort data based on selectivity index\n",
    "    driver_selective = merged_rlf[\n",
    "        (abs(merged_rlf[\"fixedazi\"]) <= 45) & (merged_rlf[\"selectivity_index\"] > 0.30)\n",
    "    ]\n",
    "    competitor_selective = merged_rlf[\n",
    "        (abs(merged_rlf[\"fixedazi\"]) <= 45) & (merged_rlf[\"selectivity_index\"] < -0.50)\n",
    "    ]\n",
    "\n",
    "    driver_means = driver_selective.groupby(\"relative_level\")[\n",
    "        \"change_in_response\"\n",
    "    ].mean()\n",
    "    competitor_means = competitor_selective.groupby(\"relative_level\")[\n",
    "        \"change_in_response\"\n",
    "    ].mean()\n",
    "\n",
    "    driver_sems = driver_selective.groupby(\"relative_level\")[\"change_in_response\"].sem()\n",
    "    competitor_sems = competitor_selective.groupby(\"relative_level\")[\n",
    "        \"change_in_response\"\n",
    "    ].sem()\n",
    "\n",
    "    relative_levels = np.unique(driver_means.index)\n",
    "\n",
    "    # fig, ax = plt.subplots(1,2)\n",
    "    driver_plot = axs[0].errorbar(\n",
    "        driver_means.index.values,\n",
    "        driver_means,\n",
    "        driver_sems,\n",
    "        color=\"orange\",\n",
    "        label=\"Driver Selective\",\n",
    "    )\n",
    "    competitor_plot = axs[1].errorbar(\n",
    "        competitor_means.index.values,\n",
    "        competitor_means,\n",
    "        competitor_sems,\n",
    "        color=\"b\",\n",
    "        label=\"Competitor Selective\",\n",
    "    )\n",
    "\n",
    "    axs[0].yaxis.set_major_formatter(\n",
    "        mpl.ticker.FuncFormatter(lambda x, pos: f\"{x:.0%}\")\n",
    "    )\n",
    "    axs[1].yaxis.set_major_formatter(\n",
    "        mpl.ticker.FuncFormatter(lambda x, pos: f\"{x:.0%}\")\n",
    "    )\n",
    "\n",
    "    axs[0].axhline(0, ls=\":\", color=\"k\")\n",
    "    axs[0].set_xticks(relative_levels.astype(int), relative_levels.astype(int))\n",
    "    axs[0].set_xlabel(\"Relative Level [dB]\")\n",
    "    axs[0].set_ylabel(\"Change in\\nSpike Rate\")\n",
    "    axs[0].set_ylim(-0.60, 0.1)\n",
    "\n",
    "    axs[1].axhline(0, ls=\":\", color=\"k\")\n",
    "    axs[1].set_xticks(relative_levels.astype(int), relative_levels.astype(int))\n",
    "    axs[1].set_xlabel(\"Relative Level [dB]\")\n",
    "    axs[1].set_ylabel(\"Change in\\nSpike Rate\")\n",
    "\n",
    "    axs[1].set_ylim(-0.10, 2.00)\n",
    "\n",
    "    axs[0].text(\n",
    "        0.02,\n",
    "        0.20,\n",
    "        f\"Frontal Driver+Lateral competitor\",\n",
    "        ha=\"left\",\n",
    "        va=\"top\",\n",
    "        fontsize=10,\n",
    "        transform=axs[0].transAxes,\n",
    "    )\n",
    "\n",
    "    axs[1].text(\n",
    "        0.02,\n",
    "        0.20,\n",
    "        f\"Frontal Driver+Lateral competitor\",\n",
    "        ha=\"left\",\n",
    "        va=\"top\",\n",
    "        fontsize=10,\n",
    "        transform=axs[1].transAxes,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2562d01d",
   "metadata": {},
   "source": [
    "## Modeling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "72c6c65a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(RL, b0, b1):\n",
    "    return 1 / (1 + np.exp(-(b0 + b1 * RL)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "48e3977b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def logit(p):\n",
    "    return np.log(p / (1 - p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "41506b3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_sigmoid_parameters(switch_type):\n",
    "\n",
    "    # Specify the relative level values where 90% and 10% of the max response are reached.\n",
    "    if switch_type:\n",
    "        relative_level_90, relative_level_10 = -10, -5\n",
    "    else:\n",
    "        relative_level_90, relative_level_10 = -15, 15\n",
    "\n",
    "    # Use the constraints that the sigmoid function takes on 90% and 10% at the points above to find its parameters.\n",
    "    p_90 = 0.9\n",
    "    p_10 = 0.1\n",
    "\n",
    "    A = np.array([[1, relative_level_90], [1, relative_level_10]])\n",
    "    b = np.array([[logit(p_90)], [logit(p_10)]])\n",
    "\n",
    "    x = np.linalg.solve(A, b)\n",
    "\n",
    "    b0 = x[0][0]\n",
    "    b1 = x[1][0]\n",
    "\n",
    "    return b0, b1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "be385710",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rlf_modeling_data(\n",
    "    single_rlf: pd.DataFrame,\n",
    "    single_rlf_out: pd.DataFrame,\n",
    "):\n",
    "\n",
    "    INTENSITY_0_EQUALS_DB_SPL = +63.0\n",
    "\n",
    "    single_in = single_rlf.copy()\n",
    "    single_out = single_rlf_out.copy()\n",
    "\n",
    "    frontal_driver = single_in[abs(single_in[\"azimuth\"]) <= 45]\n",
    "    peripheral_driver = single_in[abs(single_in[\"azimuth\"]) > 45]\n",
    "\n",
    "    frontal_competitor = single_out[abs(single_out[\"azimuth\"]) <= 45]\n",
    "    peripheral_competitor = single_out[abs(single_out[\"azimuth\"]) > 45]\n",
    "\n",
    "    merged_rlf1 = frontal_driver.join(\n",
    "        peripheral_competitor, how=\"inner\", rsuffix=\"_out\"\n",
    "    )\n",
    "    merged_rlf1 = merged_rlf1.loc[\n",
    "        merged_rlf1[\"intensity\"] == merged_rlf1[\"intensity_out\"]\n",
    "    ]\n",
    "    merged_rlf1[\"selectivity_index\"] = (\n",
    "        merged_rlf1[(merged_rlf1[\"intensity\"] == -20)][\"resp\"]\n",
    "        - merged_rlf1[(merged_rlf1[\"intensity\"] == -20)][\"resp_out\"]\n",
    "    ) / (\n",
    "        merged_rlf1[(merged_rlf1[\"intensity\"] == -20)][\"resp\"]\n",
    "        + merged_rlf1[(merged_rlf1[\"intensity\"] == -20)][\"resp_out\"]\n",
    "    )\n",
    "\n",
    "    merged_rlf2 = frontal_competitor.join(\n",
    "        peripheral_driver, how=\"inner\", rsuffix=\"_out\"\n",
    "    )\n",
    "    merged_rlf2 = merged_rlf2.loc[\n",
    "        merged_rlf2[\"intensity\"] == merged_rlf2[\"intensity_out\"]\n",
    "    ]\n",
    "    merged_rlf2[\"selectivity_index\"] = (\n",
    "        merged_rlf2[(merged_rlf2[\"intensity\"] == -20)][\"resp\"]\n",
    "        - merged_rlf2[(merged_rlf2[\"intensity\"] == -20)][\"resp_out\"]\n",
    "    ) / (\n",
    "        merged_rlf2[(merged_rlf2[\"intensity\"] == -20)][\"resp\"]\n",
    "        + merged_rlf2[(merged_rlf2[\"intensity\"] == -20)][\"resp_out\"]\n",
    "    )\n",
    "\n",
    "    tuning_index_1 = merged_rlf1.filter(\n",
    "        [\n",
    "            \"azimuth\",\n",
    "            \"elevation\",\n",
    "            \"azimuth_out\",\n",
    "            \"elevation_out\",\n",
    "            \"intensity\",\n",
    "            \"norm_resp\",\n",
    "            \"norm_resp_out\",\n",
    "            \"selectivity_index\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "    tuning_index_2 = merged_rlf2.filter(\n",
    "        [\n",
    "            \"azimuth\",\n",
    "            \"elevation\",\n",
    "            \"azimuth_out\",\n",
    "            \"elevation_out\",\n",
    "            \"intensity\",\n",
    "            \"norm_resp\",\n",
    "            \"norm_resp_out\",\n",
    "            \"selectivity_index\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "    merged_rlf = pd.concat([tuning_index_1, tuning_index_2])\n",
    "\n",
    "    quantiles = np.quantile(merged_rlf[\"selectivity_index\"], [0.25, 0.5, 0.75])\n",
    "\n",
    "    competitor_selective = merged_rlf[merged_rlf[\"selectivity_index\"] <= 0]\n",
    "\n",
    "    driver_selective = merged_rlf[(merged_rlf[\"selectivity_index\"] > 0)]\n",
    "\n",
    "    print(\"Number of units:\", np.unique(merged_rlf.index).size)\n",
    "    print(\n",
    "        \"Number of driver units:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (merged_rlf[\"selectivity_index\"] > 0) & (merged_rlf[\"intensity\"] == -20)\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "    print(\n",
    "        \"Number of competitor units:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (merged_rlf[\"selectivity_index\"] <= 0)\n",
    "                & (merged_rlf[\"intensity\"] == -20)\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "\n",
    "    print(\n",
    "        \"Number of sessions:\",\n",
    "        len(sorted(set([idx[:2] for idx in merged_rlf.index.values]))),\n",
    "    )\n",
    "    merged_rlf = merged_rlf.set_index(\"intensity\")\n",
    "\n",
    "    intensities = np.unique(merged_rlf.index) + INTENSITY_0_EQUALS_DB_SPL\n",
    "\n",
    "    DRIVERALONE_driverselective_means = driver_selective.groupby(\"intensity\")[\n",
    "        \"norm_resp\"\n",
    "    ].mean()\n",
    "    DRIVERALONE_driverselective_sems = driver_selective.groupby(\"intensity\")[\n",
    "        \"norm_resp\"\n",
    "    ].sem()\n",
    "\n",
    "    COMPETITORONE_driverselective_means = driver_selective.groupby(\"intensity\")[\n",
    "        \"norm_resp_out\"\n",
    "    ].mean()\n",
    "    COMPETITORALONE_driverselective_sems = driver_selective.groupby(\"intensity\")[\n",
    "        \"norm_resp_out\"\n",
    "    ].sem()\n",
    "\n",
    "    return (\n",
    "        intensities,\n",
    "        DRIVERALONE_driverselective_means.values,\n",
    "        COMPETITORONE_driverselective_means.values,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "1f5b3aef",
   "metadata": {},
   "outputs": [],
   "source": [
    "def competition_modeling_data(\n",
    "    twostim_rlf: pd.DataFrame,\n",
    "    single_rlf: pd.DataFrame,\n",
    "    single_rlf_out: pd.DataFrame,\n",
    "):\n",
    "\n",
    "    ## merge single rlf data and calculate selectivity index\n",
    "    single_in = single_rlf.copy()\n",
    "    single_out = single_rlf_out.copy()\n",
    "\n",
    "    frontal_driver = single_in[abs(single_in[\"azimuth\"]) <= 45]\n",
    "    peripheral_driver = single_in[abs(single_in[\"azimuth\"]) > 45]\n",
    "\n",
    "    frontal_competitor = single_out[abs(single_out[\"azimuth\"]) <= 45]\n",
    "    peripheral_competitor = single_out[abs(single_out[\"azimuth\"]) > 45]\n",
    "\n",
    "    merged_rlf1 = frontal_driver.join(\n",
    "        peripheral_competitor, how=\"inner\", rsuffix=\"_out\"\n",
    "    )\n",
    "    merged_rlf1 = merged_rlf1.loc[\n",
    "        merged_rlf1[\"intensity\"] == merged_rlf1[\"intensity_out\"]\n",
    "    ]\n",
    "    merged_rlf1[\"selectivity_index\"] = (\n",
    "        merged_rlf1[\"resp\"] - merged_rlf1[\"resp_out\"]\n",
    "    ) / (merged_rlf1[\"resp\"] + merged_rlf1[\"resp_out\"])\n",
    "\n",
    "    merged_rlf2 = frontal_competitor.join(\n",
    "        peripheral_driver, how=\"inner\", rsuffix=\"_out\"\n",
    "    )\n",
    "    merged_rlf2 = merged_rlf2.loc[\n",
    "        merged_rlf2[\"intensity\"] == merged_rlf2[\"intensity_out\"]\n",
    "    ]\n",
    "    merged_rlf2[\"selectivity_index\"] = (\n",
    "        merged_rlf2[\"resp\"] - merged_rlf2[\"resp_out\"]\n",
    "    ) / (merged_rlf2[\"resp\"] + merged_rlf2[\"resp_out\"])\n",
    "\n",
    "    tuning_index_1 = merged_rlf1.filter(\n",
    "        [\n",
    "            \"azimuth\",\n",
    "            \"elevation\",\n",
    "            \"azimuth_out\",\n",
    "            \"elevation_out\",\n",
    "            \"intensity\",\n",
    "            \"resp\",\n",
    "            \"selectivity_index\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "    tuning_index_2 = merged_rlf2.filter(\n",
    "        [\n",
    "            \"azimuth\",\n",
    "            \"elevation\",\n",
    "            \"azimuth_out\",\n",
    "            \"elevation_out\",\n",
    "            \"intensity\",\n",
    "            \"resp_out\",\n",
    "            \"selectivity_index\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "    tuning_index_2 = tuning_index_2.rename(columns={\"resp_out\": \"resp\"})\n",
    "\n",
    "    merged_1 = twostim_rlf.join(tuning_index_1, how=\"inner\", rsuffix=\"_single\")\n",
    "    merged_1 = merged_1.loc[(merged_1[\"fixedintensity\"] == merged_1[\"intensity\"])]\n",
    "    merged_1 = merged_1.loc[(merged_1[\"fixedazi\"] == merged_1[\"azimuth\"])]\n",
    "    merged_1 = merged_1.loc[(merged_1[\"fixedele\"] == merged_1[\"elevation\"])]\n",
    "\n",
    "    merged_1 = merged_1.loc[(merged_1[\"varyingazi\"] == merged_1[\"azimuth_out\"])]\n",
    "    merged_1 = merged_1.loc[(merged_1[\"varyingele\"] == merged_1[\"elevation_out\"])]\n",
    "\n",
    "    merged_2 = twostim_rlf.join(tuning_index_2, how=\"inner\", rsuffix=\"_single\")\n",
    "    merged_2 = merged_2.loc[(merged_2[\"fixedintensity\"] == merged_2[\"intensity\"])]\n",
    "\n",
    "    merged_2 = merged_2.loc[(merged_2[\"fixedazi\"] == merged_2[\"azimuth_out\"])]\n",
    "    merged_2 = merged_2.loc[(merged_2[\"fixedele\"] == merged_2[\"elevation_out\"])]\n",
    "    merged_2 = merged_2.loc[(merged_2[\"varyingazi\"] == merged_2[\"azimuth\"])]\n",
    "    merged_2 = merged_2.loc[(merged_2[\"varyingele\"] == merged_2[\"elevation\"])]\n",
    "    merged_rlf = pd.concat([merged_1, merged_2])\n",
    "\n",
    "    quantiles = np.quantile(merged_rlf[\"selectivity_index\"], [0.25, 0.5, 0.75])\n",
    "\n",
    "    print(\n",
    "        \"Number of frontal driver units + frontal driver:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (abs(merged_rlf[\"fixedazi\"]) <= 45)\n",
    "                & (abs(merged_rlf[\"varyingazi\"]) > 45)\n",
    "                & (merged_rlf[\"selectivity_index\"] > 0)\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "    print(\n",
    "        \"Number of frontal driver units + lateral driver:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (abs(merged_rlf[\"fixedazi\"]) > 45)\n",
    "                & (abs(merged_rlf[\"varyingazi\"]) <= 45)\n",
    "                & (merged_rlf[\"selectivity_index\"] > 0)\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "\n",
    "    ##lateral driver selective units\n",
    "    print(\n",
    "        \"Number of lateral driver units + frontal driver:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (abs(merged_rlf[\"fixedazi\"]) <= 45)\n",
    "                & (abs(merged_rlf[\"varyingazi\"]) > 45)\n",
    "                & (merged_rlf[\"selectivity_index\"] <= 0)\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "    print(\n",
    "        \"Number of lateral driver units + lateral driver:\",\n",
    "        np.unique(\n",
    "            merged_rlf[\n",
    "                (abs(merged_rlf[\"fixedazi\"]) > 45)\n",
    "                & (abs(merged_rlf[\"varyingazi\"]) <= 45)\n",
    "                & (merged_rlf[\"selectivity_index\"] <= 0)\n",
    "            ].index\n",
    "        ).size,\n",
    "    )\n",
    "\n",
    "    merged_rlf[\"change_in_response\"] = (\n",
    "        merged_rlf[\"resp\"] - merged_rlf[\"resp_single\"]\n",
    "    ) / merged_rlf[\"resp_single\"]\n",
    "\n",
    "    merged_rlf.set_index(\"relative_level\", inplace=True)\n",
    "\n",
    "    ##sort data based on frontal driver + lateral competitor\n",
    "    frontal_driver = merged_rlf[\n",
    "        (abs(merged_rlf[\"fixedazi\"]) < 45)\n",
    "        & (abs(merged_rlf[\"varyingazi\"]) >= 45)\n",
    "        & (merged_rlf[\"selectivity_index\"] > 0)\n",
    "    ]\n",
    "\n",
    "    peripheral_driver = merged_rlf[\n",
    "        (abs(merged_rlf[\"fixedazi\"]) > 45)\n",
    "        & (abs(merged_rlf[\"varyingazi\"]) <= 45)\n",
    "        & (merged_rlf[\"selectivity_index\"] > 0)\n",
    "    ]\n",
    "\n",
    "    frontal_driver_means = frontal_driver.groupby(\"relative_level\")[\n",
    "        \"change_in_response\"\n",
    "    ].mean()\n",
    "    peripheral_driver_means = peripheral_driver.groupby(\"relative_level\")[\n",
    "        \"change_in_response\"\n",
    "    ].mean()\n",
    "\n",
    "    frontal_driver_sems = frontal_driver.groupby(\"relative_level\")[\n",
    "        \"change_in_response\"\n",
    "    ].sem()\n",
    "    peripheral_driver_sems = peripheral_driver.groupby(\"relative_level\")[\n",
    "        \"change_in_response\"\n",
    "    ].sem()\n",
    "\n",
    "    relative_levels = np.unique(frontal_driver_means.index)\n",
    "\n",
    "    return relative_levels, frontal_driver_means.values, peripheral_driver_means.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "472ceb92",
   "metadata": {},
   "outputs": [],
   "source": [
    "def model(relative_resp, axs=None):\n",
    "    driver_level = 43\n",
    "    competitor_level = driver_level + relative_level\n",
    "\n",
    "    r_d = rate_driver_4b[sound_level == driver_level]\n",
    "    r_c = rate_competitor_4b[np.in1d(sound_level, competitor_level)]\n",
    "    relative_response = relative_resp\n",
    "    switch_type = True\n",
    "    b0, b1 = get_sigmoid_parameters(switch_type)\n",
    "\n",
    "    gain_d = 0.7 * sigmoid(relative_level, b0, b1) + 0.3\n",
    "\n",
    "    gain_c = r_d * (1 + relative_response - gain_d) / r_c\n",
    "\n",
    "    axs[0].plot(relative_level, gain_c, color=\"b\")\n",
    "    axs[0].plot(relative_level, gain_d, color=\"orange\")\n",
    "    axs[0].axhline(1.0, ls=\":\", color=\"k\")\n",
    "    axs[0].set_ylim(-0.2, 2)\n",
    "\n",
    "    axs[0].set_xlabel(\"Relative level [dB]\")\n",
    "    axs[0].set_ylabel(\"Gain on input\")\n",
    "    axs[0].legend(\n",
    "        [\"Competitor position\", \"Driver position\"], frameon=False, loc=\"upper left\"\n",
    "    )\n",
    "    axs[0].set_xticks(\n",
    "        [int(lev) for lev in relative_level], [int(lev) for lev in relative_level]\n",
    "    )\n",
    "\n",
    "    r = r_d * gain_d + r_c * gain_c\n",
    "\n",
    "    percent_change_d = (r - r_d) / r_d  # * 100\n",
    "\n",
    "    axs[1].plot(relative_level, percent_change_d, color=\"k\")\n",
    "    axs[1].axhline(0, ls=\":\", color=\"k\")\n",
    "\n",
    "    axs[1].set_ylim(-1, 1)\n",
    "    axs[1].set_xticks(\n",
    "        [int(lev) for lev in relative_level], [int(lev) for lev in relative_level]\n",
    "    )\n",
    "    axs[1].set_xlabel(\"Relative level [dB]\")\n",
    "    axs[1].set_ylabel(\"Change in\\nSpike Rate\")\n",
    "    axs[1].yaxis.set_major_formatter(\n",
    "        mpl.ticker.FuncFormatter(lambda x, pos: f\"{x:.0%}\")\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6fb3051",
   "metadata": {},
   "source": [
    "## SRF correlation distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "ac85a14a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def srf_correlation_distribution(single_srf_df, twostim_ccg_df, axs=None):\n",
    "    if axs is None:\n",
    "        fig, axs = plt.subplots(1, 1)\n",
    "    srf_df = single_srf_df.copy()\n",
    "    twostim_ccg_df = twostim_ccg_df.copy()\n",
    "\n",
    "    twostim_ccg_df = twostim_ccg_df[twostim_ccg_df[\"corr_type\"] == \"cross_region\"]\n",
    "\n",
    "    srf_corr_df = []\n",
    "    for index in np.unique(twostim_ccg_df.index):\n",
    "        try:\n",
    "            date = index[0]\n",
    "            owl = index[1]\n",
    "            channel_ot = index[-2]\n",
    "            channel_fore = index[-1]\n",
    "\n",
    "            channel_ot_srf = srf_df.loc[date, owl, channel_ot]\n",
    "            channel_fore_srf = srf_df.loc[date, owl, channel_fore]\n",
    "\n",
    "            srf_corr = np.corrcoef(\n",
    "                channel_ot_srf[\"norm_resp\"], channel_fore_srf[\"norm_resp\"]\n",
    "            )[0][1]\n",
    "\n",
    "            tmp = {\n",
    "                \"date\": date,\n",
    "                \"owl\": owl,\n",
    "                \"channel1\": channel_ot,\n",
    "                \"channel2\": channel_fore,\n",
    "                \"srf_corr\": srf_corr,\n",
    "            }\n",
    "            srf_corr_df.append(tmp)\n",
    "        except:\n",
    "            pass\n",
    "    srf_corr_df = pd.DataFrame(srf_corr_df)\n",
    "    srf_corr_df = srf_corr_df.set_index([\"date\", \"owl\", \"channel1\", \"channel2\"])\n",
    "    srf_corr_df = srf_corr_df.filter(\n",
    "        [\n",
    "            \"srf_corr\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "    merge = twostim_ccg_df.join(srf_corr_df, how=\"inner\", rsuffix=\"_srf\").set_index(\n",
    "        \"relative_level\"\n",
    "    )\n",
    "    ##for just one level\n",
    "    merge = merge.loc[0, \"srf_corr\"]\n",
    "\n",
    "    mean_srf_corr = np.mean(merge)\n",
    "    std_srf_corr = np.std(merge)\n",
    "    confidence_interval = confidence_intervale(merge)\n",
    "\n",
    "    print(f\"Mean {mean_srf_corr:.2f}\")\n",
    "    # print(f'CI {confidence_interval:.2f}')\n",
    "    print(f\"STD {std_srf_corr:.2f}\")\n",
    "\n",
    "    quantiles = np.quantile(merge, [0.25, 0.50, 0.75])\n",
    "    colors = [\"slategrey\", \"royalblue\", \"blue\", \"navy\"]\n",
    "\n",
    "    bins = np.arange(-1.0, 1.0, 0.01)\n",
    "    # plt.figure()\n",
    "    axs.hist(\n",
    "        merge[(merge <= quantiles[0])], bins, facecolor=colors[0], edgecolor=\"None\"\n",
    "    )\n",
    "    axs.hist(\n",
    "        merge[(merge > quantiles[0]) & ((merge <= quantiles[1]))],\n",
    "        bins,\n",
    "        facecolor=colors[1],\n",
    "        edgecolor=\"None\",\n",
    "    )\n",
    "    axs.hist(\n",
    "        merge[(merge > quantiles[1]) & ((merge <= quantiles[2]))],\n",
    "        bins,\n",
    "        facecolor=colors[2],\n",
    "        edgecolor=\"None\",\n",
    "    )\n",
    "    axs.hist(merge[(merge > quantiles[2])], bins, facecolor=colors[3], edgecolor=\"None\")\n",
    "\n",
    "    axs.set_ylabel(\"Number of unit pairs\")\n",
    "    axs.set_xlabel(\"Spatial Tuning Correlation between OT and nRt\")\n",
    "    axs.spines[\"right\"].set_visible(False)\n",
    "    axs.spines[\"top\"].set_visible(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "243670ae",
   "metadata": {},
   "source": [
    "## SRF correlation by frontal OT and peripheral OT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "69567c9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def srf_correlation_frontalperiperpheral(single_srf_df, twostim_ccg_df, axs=None):\n",
    "    if axs is None:\n",
    "        fig, axs = plt.subplots(1, 1)\n",
    "    srf_df = single_srf_df.copy()\n",
    "    twostim_ccg_df = twostim_ccg_df.copy()\n",
    "\n",
    "    twostim_ccg_df = twostim_ccg_df[twostim_ccg_df[\"corr_type\"] == \"cross_region\"]\n",
    "\n",
    "    srf_corr_df = []\n",
    "    for index in np.unique(twostim_ccg_df.index):\n",
    "        try:\n",
    "            date = index[0]\n",
    "            owl = index[1]\n",
    "            channel_ot = index[-2]\n",
    "            channel_fore = index[-1]\n",
    "\n",
    "            channel_ot_srf = srf_df.loc[date, owl, channel_ot]\n",
    "            channel_fore_srf = srf_df.loc[date, owl, channel_fore]\n",
    "\n",
    "            srf_corr = np.corrcoef(\n",
    "                channel_ot_srf[\"norm_resp\"], channel_fore_srf[\"norm_resp\"]\n",
    "            )[0][1]\n",
    "\n",
    "            tmp = {\n",
    "                \"date\": date,\n",
    "                \"owl\": owl,\n",
    "                \"channel1\": channel_ot,\n",
    "                \"channel2\": channel_fore,\n",
    "                \"srf_corr\": srf_corr,\n",
    "            }\n",
    "            srf_corr_df.append(tmp)\n",
    "        except:\n",
    "            pass\n",
    "    srf_corr_df = pd.DataFrame(srf_corr_df)\n",
    "    zz = srf_corr_df.copy()\n",
    "    srf_corr_df = srf_corr_df.set_index([\"date\", \"owl\", \"channel1\", \"channel2\"])\n",
    "\n",
    "    unique_ot = zz.set_index([\"date\", \"owl\", \"channel1\"])\n",
    "    unique_forebrain = zz.set_index([\"date\", \"owl\", \"channel2\"])\n",
    "    print(\"Number of OT units:\", np.unique(unique_ot.index).size)\n",
    "    print(\"Number of Forebrain units:\", np.unique(unique_forebrain.index).size)\n",
    "\n",
    "    srf_corr_df = srf_corr_df.filter(\n",
    "        [\n",
    "            \"srf_corr\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "    merge = twostim_ccg_df.join(srf_corr_df, how=\"inner\", rsuffix=\"_srf\")\n",
    "    merge = merge[merge[\"relative_level\"] == 0]\n",
    "\n",
    "    print(\n",
    "        \"Number of sessions:\",\n",
    "        len(sorted(set([idx[:2] for idx in merge.index.values]))),\n",
    "    )\n",
    "\n",
    "    bins = np.arange(-1.0, 1.0, 0.01)\n",
    "    frontal_ot_units = merge[(abs(merge[\"fixedazi\"]) <= 45)][\"srf_corr\"]\n",
    "    peripheral_ot_units = merge[(abs(merge[\"fixedazi\"]) > 45)][\"srf_corr\"]\n",
    "\n",
    "    CI_frontal = confidence_intervale(frontal_ot_units)\n",
    "    CI_peripheral = confidence_intervale(peripheral_ot_units)\n",
    "\n",
    "    tval, p = scipy.stats.ttest_ind(frontal_ot_units, peripheral_ot_units)\n",
    "\n",
    "    print(tval, p)\n",
    "    print(f\"frontal stats: {np.mean(frontal_ot_units)} +/- {CI_frontal}\")\n",
    "    print(f\"peripheral stats: {np.mean(peripheral_ot_units)} +/- {CI_peripheral}\")\n",
    "\n",
    "    axs.hist(\n",
    "        peripheral_ot_units,\n",
    "        bins,\n",
    "        edgecolor=\"r\",\n",
    "        facecolor=\"r\",\n",
    "        alpha=0.5,\n",
    "        density=True,\n",
    "        label=\"lateral\",\n",
    "    )\n",
    "    axs.hist(\n",
    "        frontal_ot_units,\n",
    "        bins,\n",
    "        edgecolor=\"b\",\n",
    "        facecolor=\"b\",\n",
    "        alpha=0.5,\n",
    "        density=True,\n",
    "        label=\"frontal\",\n",
    "    )\n",
    "    axs.plot(np.mean(frontal_ot_units), 8, \"v\", color=\"b\")\n",
    "    axs.plot(np.mean(peripheral_ot_units), 8, \"v\", color=\"r\")\n",
    "    axs.legend(frameon=False)\n",
    "    axs.set_xlabel(\"Spatial tuning correlation\\nbetween OT and Forebrain\")\n",
    "    axs.set_ylabel(\"Density\")\n",
    "    axs.spines[\"right\"].set_visible(False)\n",
    "    axs.spines[\"top\"].set_visible(False)\n",
    "    axs.text(\n",
    "        0,\n",
    "        1,\n",
    "        f\"p = {p:.3g}\",\n",
    "        horizontalalignment=\"left\",\n",
    "        verticalalignment=\"top\",\n",
    "        fontsize=10,\n",
    "        transform=axs.transAxes\n",
    "        + mpl.transforms.ScaledTranslation(\n",
    "            +2 / 72,\n",
    "            0 / 72,\n",
    "            axs.figure.dpi_scale_trans,\n",
    "        ),\n",
    "    )\n",
    "    axs.set_ylim(bottom=0, top=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "155884ec",
   "metadata": {},
   "source": [
    "## VS competitior as a function of VS driver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "e2ac70e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def figure_stimphase_driver_competitor(\n",
    "    df: pd.DataFrame,\n",
    "    axs=None,\n",
    "    colors=[c_forebrain_am55, c_forebrain_am75],\n",
    "    axs_pol=None,\n",
    "):\n",
    "    if axs is None:\n",
    "        fig, axs = plt.subplots(1, df[\"relative_level\"].unique().size)\n",
    "    axs = axs.flatten()\n",
    "\n",
    "    df = df.set_index(\"relative_level\")\n",
    "\n",
    "    for k, level in enumerate(np.unique(df.index)):\n",
    "        ax = axs[k]\n",
    "\n",
    "        level_rayleigh_p = df.loc[level, \"fixedstim_plv_p\"]\n",
    "\n",
    "        level_rayleigh_p_varying = df.loc[level, \"varyingstim_plv_p\"]\n",
    "        level_plv_driver = df.loc[level, \"fixedstim_plv\"]\n",
    "        level_plv_competitor = df.loc[level, \"varyingstim_plv\"]\n",
    "\n",
    "        # Filter phases:\n",
    "        mask = level_rayleigh_p < 0.05\n",
    "        #        mask_comp = level_rayleigh_p_varying < 0.05\n",
    "\n",
    "        ax.plot(\n",
    "            level_plv_driver[~mask],\n",
    "            level_plv_competitor[~mask],\n",
    "            ls=\"none\",\n",
    "            marker=\".\",\n",
    "            color=\".5\",\n",
    "            mec=\"none\",\n",
    "        )\n",
    "\n",
    "        ax.plot(\n",
    "            level_plv_driver[mask],\n",
    "            level_plv_competitor[mask],\n",
    "            ls=\"none\",\n",
    "            marker=\".\",\n",
    "            color=colors[0][0],\n",
    "            mec=\"none\",\n",
    "            zorder=-1,\n",
    "        )\n",
    "\n",
    "        ax.plot([0, 0.5], [0, 0.5], color=\"k\", lw=0.8, ls=\"--\")\n",
    "        ax.text(\n",
    "            0,\n",
    "            1,\n",
    "            f\"{np.mean(mask):.0%}\",\n",
    "            ha=\"left\",\n",
    "            va=\"top\",\n",
    "            fontsize=8,\n",
    "            color=colors[0][0],\n",
    "            transform=ax.transAxes\n",
    "            + mpl.transforms.ScaledTranslation(\n",
    "                +1 / 72,\n",
    "                -1 / 72,\n",
    "                ax.figure.dpi_scale_trans,\n",
    "            ),\n",
    "        )\n",
    "\n",
    "        ax.set_title(f\"{level:.0f}\", fontsize=10)\n",
    "\n",
    "        ax.set_xlim(left=0, right=0.5)\n",
    "        ax.set_xticks([0, 0.5], [])\n",
    "        ax.set_xticks([0.25], minor=True)\n",
    "\n",
    "        ax.set_ylim(bottom=0, top=0.5)\n",
    "        ax.set_yticks([0, 0.5], [])\n",
    "        ax.set_yticks([0.25], minor=True)\n",
    "        ax.spines[\"right\"].set_visible(False)\n",
    "        ax.spines[\"top\"].set_visible(False)\n",
    "        if k == 0:\n",
    "            ax.set_yticklabels([\"0\", \"0.5\"])\n",
    "            ax.set_xticklabels([\"0\", \"0.5\"])\n",
    "        if k == 3:\n",
    "            ax.set_xlabel(\"VS Driver\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17dfc79c",
   "metadata": {},
   "source": [
    "## AM VS change with competition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "b1864327",
   "metadata": {},
   "outputs": [],
   "source": [
    "def figure_stimphase_switch(\n",
    "    df: pd.DataFrame,\n",
    "    ax=None,\n",
    "    colors=[c_forebrain_am55, c_forebrain_am75],\n",
    "    axs_pol=None,\n",
    "):\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots(1, 1)\n",
    "\n",
    "    driver_vs = []\n",
    "    competitor_vs = []\n",
    "    for index in df.index:\n",
    "        sub = df.loc[index]\n",
    "        driver_plv = list(sub.fixedstim_plv.values)\n",
    "        competitor_plv = list(sub.varyingstim_plv.values)\n",
    "        driver_vs.append(driver_plv)\n",
    "        if len(driver_plv) > 6:\n",
    "            display(sub)\n",
    "        competitor_vs.append(competitor_plv)\n",
    "\n",
    "    print([len(x) for x in driver_vs])\n",
    "\n",
    "    ax.errorbar(\n",
    "        np.unique(df[\"relative_level\"]),\n",
    "        np.mean(driver_vs, axis=0),\n",
    "        scipy.stats.sem(driver_vs),\n",
    "        color=colors[0][0],\n",
    "        label=\"Driver\",\n",
    "    )\n",
    "    ax.errorbar(\n",
    "        np.unique(df[\"relative_level\"]),\n",
    "        np.mean(competitor_vs, axis=0),\n",
    "        scipy.stats.sem(competitor_vs),\n",
    "        color=colors[1][0],\n",
    "        label=\"Competitor\",\n",
    "    )\n",
    "    ax.set_xticks(\n",
    "        np.unique(df[\"relative_level\"]),\n",
    "        [int(val) for val in np.unique(df[\"relative_level\"])],\n",
    "    )\n",
    "    ax.set_xlabel(\"Relative Level [dB]\")\n",
    "    ax.set_ylabel(\"Vector Strength\")\n",
    "    ax.spines[\"right\"].set_visible(False)\n",
    "    ax.spines[\"top\"].set_visible(False)\n",
    "\n",
    "    fixfreq = df.iloc[0][\"fixed_modulation_frequency\"]\n",
    "    varfreq = df.iloc[0][\"varying_modulation_frequency\"]\n",
    "\n",
    "    mylabels = [f\"{fixfreq}Hz driver\", f\"{varfreq}Hz competitor\"]\n",
    "\n",
    "    ax.legend(labels=mylabels, frameon=False, loc=\"lower left\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e4ebe39",
   "metadata": {},
   "source": [
    "## Dual region xcorr by srf correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "d7f4ae1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fig_dual_region_xcorr_grid(single_srf_df, twostim_ccg_df, axs=None, caxs=None):\n",
    "    if axs is None:\n",
    "        fig, axs = plt.subplots(4, 6)\n",
    "    # axs = axs.flatten()\n",
    "\n",
    "    srf_df = single_srf_df.copy()\n",
    "    twostim_ccg_df = twostim_ccg_df.copy()\n",
    "\n",
    "    twostim_ccg_df = twostim_ccg_df[twostim_ccg_df[\"corr_type\"] == \"cross_region\"]\n",
    "\n",
    "    srf_corr_df = []\n",
    "    for index in np.unique(twostim_ccg_df.index):\n",
    "        try:\n",
    "            date = index[0]\n",
    "            owl = index[1]\n",
    "            channel_ot = index[-2]\n",
    "            channel_fore = index[-1]\n",
    "\n",
    "            channel_ot_srf = srf_df.loc[date, owl, channel_ot]\n",
    "            channel_fore_srf = srf_df.loc[date, owl, channel_fore]\n",
    "\n",
    "            srf_corr = np.corrcoef(\n",
    "                channel_ot_srf[\"norm_resp\"], channel_fore_srf[\"norm_resp\"]\n",
    "            )[0][1]\n",
    "\n",
    "            tmp = {\n",
    "                \"date\": date,\n",
    "                \"owl\": owl,\n",
    "                \"channel1\": channel_ot,\n",
    "                \"channel2\": channel_fore,\n",
    "                \"srf_corr\": srf_corr,\n",
    "            }\n",
    "            srf_corr_df.append(tmp)\n",
    "        except:\n",
    "            pass\n",
    "    srf_corr_df = pd.DataFrame(srf_corr_df)\n",
    "    srf_corr_df = srf_corr_df.set_index([\"date\", \"owl\", \"channel1\", \"channel2\"])\n",
    "    srf_corr_df = srf_corr_df.filter(\n",
    "        [\n",
    "            \"srf_corr\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "    merge = twostim_ccg_df.join(srf_corr_df, how=\"inner\", rsuffix=\"_srf\")\n",
    "\n",
    "    quantiles = np.quantile(merge[\"srf_corr\"], [0.25, 0.50, 0.75])\n",
    "\n",
    "    low_corr = merge[(merge[\"srf_corr\"] <= quantiles[0])].set_index(\"relative_level\")\n",
    "    lower_mid_corr = merge[\n",
    "        (merge[\"srf_corr\"] > quantiles[0]) & ((merge[\"srf_corr\"] <= quantiles[1]))\n",
    "    ].set_index(\"relative_level\")\n",
    "    higher_mid_corr = merge[\n",
    "        (merge[\"srf_corr\"] > quantiles[1]) & ((merge[\"srf_corr\"] <= quantiles[2]))\n",
    "    ].set_index(\"relative_level\")\n",
    "    high_corr = merge[(merge[\"srf_corr\"] > quantiles[2])].set_index(\"relative_level\")\n",
    "\n",
    "    merge = merge.set_index(\"relative_level\")\n",
    "\n",
    "    psth_len = (merge.iloc[0][\"ccg\"].size + 1) / 2\n",
    "    lags = scipy.signal.correlation_lags(psth_len, psth_len) / 1000\n",
    "    lags_mask = np.abs(lags) <= 0.95\n",
    "    lags_mask_std = np.abs(lags) >= 0.825\n",
    "\n",
    "    relative_levels = np.unique(merge.index)\n",
    "\n",
    "    groupdata = [low_corr, lower_mid_corr, higher_mid_corr, high_corr]\n",
    "    label = [\"lowcorr\", \"lowermid\", \"highermid\", \"high\"]\n",
    "    colors = [\"slategrey\", \"royalblue\", \"blue\", \"navy\"]\n",
    "    display(high_corr.head(5))\n",
    "    for ndata, data in enumerate(groupdata):\n",
    "        for k, relative_level in enumerate(relative_levels):\n",
    "            mean_ccg = np.mean(np.vstack(data.loc[relative_level, \"ccg\"]), axis=0)[::-1]\n",
    "\n",
    "            mean_n_ccg = data.loc[relative_level, \"ccg\"].count()\n",
    "\n",
    "            if relative_level == -15:\n",
    "                peak_properties = get_peak(mean_ccg, lags)\n",
    "                peak_time = round(peak_properties[\"peak_lag\"], 3)\n",
    "                print(label[ndata], relative_level, peak_time)\n",
    "\n",
    "            axs[ndata, k].plot(\n",
    "                lags[lags_mask], mean_ccg[lags_mask], color=colors[ndata], ls=\"-\", lw=1\n",
    "            )\n",
    "\n",
    "            axs[ndata, k].text(\n",
    "                0.5,\n",
    "                0.95,\n",
    "                f\"{mean_n_ccg}\",\n",
    "                ha=\"center\",\n",
    "                va=\"top\",\n",
    "                fontsize=8,\n",
    "                transform=axs[ndata, k].transAxes,\n",
    "            )\n",
    "            axs[ndata, k].spines[\"top\"].set_visible(False)\n",
    "            axs[ndata, k].spines[\"right\"].set_visible(False)\n",
    "            # axs[k].axvline(0, ls = ':', color = 'r')\n",
    "            axs[ndata, k].set_ylim(bottom=-1.5e-6, top=6e-6)\n",
    "            axs[ndata, k].set_xlim(left=-0.950, right=0.950)\n",
    "            axs[ndata, 0].set_ylabel(\"Coinc./spk\", labelpad=12)\n",
    "\n",
    "            if ndata == 0:\n",
    "                axs[ndata, k].set_title(relative_level)\n",
    "                axs[ndata, 0].text(\n",
    "                    0.80,\n",
    "                    0.80,\n",
    "                    f\"OT to\\nnRt\",\n",
    "                    ha=\"center\",\n",
    "                    va=\"top\",\n",
    "                    fontsize=8,\n",
    "                    transform=axs[ndata, k].transAxes,\n",
    "                )\n",
    "                axs[ndata, 0].text(\n",
    "                    0.25,\n",
    "                    0.80,\n",
    "                    f\"nRt to\\nOT\",\n",
    "                    ha=\"center\",\n",
    "                    va=\"top\",\n",
    "                    fontsize=8,\n",
    "                    transform=axs[ndata, k].transAxes,\n",
    "                )\n",
    "\n",
    "                axs[ndata, 0].axvline(0, ls=\":\", color=\"k\", lw=0.4, zorder=-1)\n",
    "\n",
    "        for ax in axs[ndata, 1:]:\n",
    "            ax.set_yticks([])\n",
    "            ax.set_yticks([], minor=True)\n",
    "            ax.spines[\"left\"].set_visible(False)\n",
    "\n",
    "        if ndata == 0:\n",
    "            axs[0, 0].set_xticks([-0.95, 0, +0.95], [\"-950\", \"0\", \"+950\"], fontsize=8)\n",
    "            axs[0, 0].set_xlabel(\"Lags [ms]\", fontsize=8, labelpad=0)\n",
    "            axs[0, 0].set_ylabel(\"Coinc./spk\", labelpad=12)\n",
    "            for ax in axs[ndata, 1:]:\n",
    "                ax.set_xticks([-0.95, 0, 0.95], [])\n",
    "\n",
    "        elif ndata != 0:\n",
    "            for ax in axs[ndata, :]:\n",
    "                ax.set_xticks([-0.95, 0, 0.95], [])\n",
    "\n",
    "        def formatter_y(x, pos):\n",
    "            if x > 0:\n",
    "                return f\"${{{x*1e6:.0f}}}$\"\n",
    "                # return f\"$10^{{{np.log10(x):.0f}}}$\"\n",
    "            else:\n",
    "                return \"0\"\n",
    "\n",
    "        axs[ndata, 0].yaxis.set_major_formatter(mpl.ticker.FuncFormatter(formatter_y))\n",
    "        axs[ndata, 0].set_yticks(np.array([0, 2.5, 5]) * 1e-6, minor=True)\n",
    "        axs[ndata, 0].text(\n",
    "            0,\n",
    "            1,\n",
    "            \"$10^{-6}$\",\n",
    "            ha=\"right\",\n",
    "            va=\"center\",\n",
    "            fontsize=8,\n",
    "            transform=axs[ndata, 0].transAxes\n",
    "            + mpl.transforms.ScaledTranslation(\n",
    "                -2 / 72,\n",
    "                0,\n",
    "                axs[ndata, 0].figure.dpi_scale_trans,\n",
    "            ),\n",
    "        )\n",
    "    cmap = mpl.colors.ListedColormap(colors[::-1]).with_extremes(\n",
    "        under=colors[0], over=colors[-1]\n",
    "    )\n",
    "\n",
    "    qrt = [0, 0.25, 0.5, 0.75, 1.0]\n",
    "    qrt = [1.0, 0.75, 0.5, 0.25, 0]\n",
    "    norm = mpl.colors.BoundaryNorm(qrt, cmap.N)\n",
    "    plt.colorbar(\n",
    "        mpl.cm.ScalarMappable(cmap=cmap, norm=norm),\n",
    "        cax=caxs,\n",
    "        spacing=\"proportional\",\n",
    "        label=\"Spatial Tuning Correlation Quartile\",\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00b6503b",
   "metadata": {},
   "source": [
    "## Dual Region CSI by relative level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "6df26a02",
   "metadata": {},
   "outputs": [],
   "source": [
    "def csi_boxplot_dualregion(\n",
    "    single_srf_df,\n",
    "    twostim_ccg_df,\n",
    "    single_ccg_df,\n",
    "    ax=None,\n",
    "    brackets: dict[tuple[int, int] : float] | None = None,\n",
    "    anova_align=\"right\",\n",
    "):\n",
    "    srf_df = single_srf_df.copy()\n",
    "    twostim_ccg_df = twostim_ccg_df.copy()\n",
    "    single_ccg_df = single_ccg_df.copy()\n",
    "\n",
    "    if ax == None:\n",
    "        fig, ax = plt.subplots(1, len(np.unique(twostim_ccg_df[\"relative_level\"])))\n",
    "    # axs.flatten()\n",
    "\n",
    "    twostim_ccg_df = twostim_ccg_df[twostim_ccg_df[\"corr_type\"] == \"cross_region\"]\n",
    "    single_ccg_df = single_ccg_df[single_ccg_df[\"corr_type\"] == \"cross_region\"]\n",
    "\n",
    "    srf_corr_df = []\n",
    "    for index in np.unique(twostim_ccg_df.index):\n",
    "        try:\n",
    "            date = index[0]\n",
    "            owl = index[1]\n",
    "            channel_ot = index[-2]\n",
    "            channel_fore = index[-1]\n",
    "\n",
    "            channel_ot_srf = srf_df.loc[date, owl, channel_ot]\n",
    "            channel_fore_srf = srf_df.loc[date, owl, channel_fore]\n",
    "\n",
    "            srf_corr = np.corrcoef(\n",
    "                channel_ot_srf[\"norm_resp\"], channel_fore_srf[\"norm_resp\"]\n",
    "            )[0][1]\n",
    "\n",
    "            tmp = {\n",
    "                \"date\": date,\n",
    "                \"owl\": owl,\n",
    "                \"channel1\": channel_ot,\n",
    "                \"channel2\": channel_fore,\n",
    "                \"srf_corr\": srf_corr,\n",
    "            }\n",
    "            srf_corr_df.append(tmp)\n",
    "        except:\n",
    "            pass\n",
    "    srf_corr_df = pd.DataFrame(srf_corr_df)\n",
    "    srf_corr_df = srf_corr_df.set_index([\"date\", \"owl\", \"channel1\", \"channel2\"])\n",
    "    srf_corr_df = srf_corr_df.filter(\n",
    "        [\n",
    "            \"srf_corr\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "    merge = twostim_ccg_df.join(srf_corr_df, how=\"inner\", rsuffix=\"_srf\")\n",
    "\n",
    "    merge = merge.join(\n",
    "        single_ccg_df,\n",
    "        on=[\"date\", \"owl\", \"channel1\", \"channel2\"],\n",
    "        how=\"inner\",\n",
    "        rsuffix=\"_single\",\n",
    "    )\n",
    "\n",
    "    merge = merge.loc[\n",
    "        (merge[\"fixedintensity\"] == merge[\"intensity\"])\n",
    "        & (merge[\"fixedazi\"] == merge[\"azimuth\"])\n",
    "        & (merge[\"fixedele\"] == merge[\"elevation\"])\n",
    "    ]\n",
    "    print(\"Number of units:\", np.unique(merge.index).size)\n",
    "    print(\n",
    "        \"Number of sessions:\",\n",
    "        len(sorted(set([idx[:2] for idx in merge.index.values]))),\n",
    "    )\n",
    "    merge[\"csi\"] = (merge[\"xcorr_peak\"] - merge[\"xcorr_peak_single\"]) / (\n",
    "        merge[\"xcorr_peak\"] + merge[\"xcorr_peak_single\"]\n",
    "    )\n",
    "\n",
    "    quantiles = np.quantile(merge[\"srf_corr\"], [0.25, 0.50, 0.75])\n",
    "\n",
    "    low_corr = merge[(merge[\"srf_corr\"] <= quantiles[0])].set_index(\"relative_level\")\n",
    "    lower_mid_corr = merge[\n",
    "        (merge[\"srf_corr\"] > quantiles[0]) & ((merge[\"srf_corr\"] <= quantiles[1]))\n",
    "    ].set_index(\"relative_level\")\n",
    "    higher_mid_corr = merge[\n",
    "        (merge[\"srf_corr\"] > quantiles[1]) & ((merge[\"srf_corr\"] <= quantiles[2]))\n",
    "    ].set_index(\"relative_level\")\n",
    "    high_corr = merge[(merge[\"srf_corr\"] > quantiles[2])].set_index(\"relative_level\")\n",
    "\n",
    "    merge = merge.set_index(\"relative_level\")\n",
    "\n",
    "    relative_levels = np.unique(merge.index)\n",
    "\n",
    "    groupdata = [\n",
    "        merge.loc[relative_level, \"csi\"].values for relative_level in relative_levels\n",
    "    ]\n",
    "\n",
    "    bp = ax.boxplot(\n",
    "        groupdata,\n",
    "        patch_artist=True,\n",
    "        notch=False,\n",
    "        showfliers=False,\n",
    "        whis=(5, 95),\n",
    "    )\n",
    "\n",
    "    ax.set_xticklabels(relative_levels.astype(int))\n",
    "\n",
    "    stats = anova_tukey(merge, val_col=\"csi\", group_col=\"relative_level\")\n",
    "    colors = [\"orange\"] * len(relative_levels)\n",
    "    for patch, color in zip(bp[\"boxes\"], colors):\n",
    "        patch.set_facecolor(color)\n",
    "\n",
    "    for median in bp[\"medians\"]:\n",
    "        median.set(color=\"white\", linewidth=2)\n",
    "\n",
    "    if True:  # show_individual_data\n",
    "        width = 0.7\n",
    "        for k, groupvalues in enumerate(groupdata):\n",
    "            pos = k + 1\n",
    "            y = np.sort(groupvalues)\n",
    "            vals_x = np.linspace(\n",
    "                pos - width / 3,\n",
    "                pos + width / 3,\n",
    "                y.size,\n",
    "                endpoint=True,\n",
    "            )\n",
    "            ax.plot(\n",
    "                vals_x,\n",
    "                y,\n",
    "                color=\"k\",\n",
    "                ls=\"None\",\n",
    "                marker=\".\",\n",
    "                markersize=2,\n",
    "                markeredgewidth=0.0,\n",
    "                zorder=5,\n",
    "            )\n",
    "\n",
    "    if True:  # show_brackets\n",
    "        brackets_shrink = 0.85\n",
    "        if brackets is None:\n",
    "            brackets = {(k, k + 1): 0.8 for k in range(len(stats[\"tukey\"]) - 1)}\n",
    "            brackets_shrink = 0.8\n",
    "        for (idx_left, idx_right), y in brackets.items():\n",
    "            plot_bracket(\n",
    "                ax,\n",
    "                left=idx_left + 1,\n",
    "                right=idx_right + 1,\n",
    "                text=f\"{stats['tukey'].iloc[idx_left, idx_right]:.3f}\"[1:],\n",
    "                y=y,\n",
    "                shrink=brackets_shrink,\n",
    "            )\n",
    "\n",
    "    ax.spines[\"top\"].set_visible(False)\n",
    "    ax.spines[\"right\"].set_visible(False)\n",
    "\n",
    "    if anova_align == \"right\":\n",
    "        ax.text(\n",
    "            1,\n",
    "            1,\n",
    "            f\"p = {stats['anova_p']:.3g}\",\n",
    "            horizontalalignment=\"right\",\n",
    "            verticalalignment=\"top\",\n",
    "            fontsize=10,\n",
    "            transform=ax.transAxes,\n",
    "        )\n",
    "    else:\n",
    "        ax.text(\n",
    "            0,\n",
    "            1,\n",
    "            f\"p = {stats['anova_p']:.3g}\",\n",
    "            horizontalalignment=\"left\",\n",
    "            verticalalignment=\"top\",\n",
    "            fontsize=10,\n",
    "            transform=ax.transAxes\n",
    "            + mpl.transforms.ScaledTranslation(\n",
    "                +2 / 72,\n",
    "                0 / 72,\n",
    "                ax.figure.dpi_scale_trans,\n",
    "            ),\n",
    "        )\n",
    "\n",
    "    # ax.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, pos: f\"{x:.0%}\"))\n",
    "\n",
    "    ax.axhline(0, color=\"k\", ls=\"--\", zorder=-1)\n",
    "    ax.set_xlabel(\"Relative Level [dB]\")\n",
    "    ax.set_ylabel(\"CSI\")\n",
    "    ax.set_ylim(-1, 1.45)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8aa07d21",
   "metadata": {},
   "source": [
    "## Dual Region CSI by SRF Corr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "a3e78645",
   "metadata": {},
   "outputs": [],
   "source": [
    "def csi_boxplot_dualregion_srfcorr(\n",
    "    single_srf_df,\n",
    "    twostim_ccg_df,\n",
    "    single_ccg_df,\n",
    "    ax=None,\n",
    "    brackets: dict[tuple[int, int] : float] | None = None,\n",
    "    anova_align=\"right\",\n",
    "):\n",
    "    srf_df = single_srf_df.copy()\n",
    "    twostim_ccg_df = twostim_ccg_df.copy()\n",
    "    single_ccg_df = single_ccg_df.copy()\n",
    "\n",
    "    if ax == None:\n",
    "        fig, ax = plt.subplots(1, len(np.unique(twostim_ccg_df[\"relative_level\"])))\n",
    "\n",
    "    twostim_ccg_df = twostim_ccg_df[twostim_ccg_df[\"corr_type\"] == \"cross_region\"]\n",
    "    single_ccg_df = single_ccg_df[single_ccg_df[\"corr_type\"] == \"cross_region\"]\n",
    "\n",
    "    srf_corr_df = []\n",
    "    for index in np.unique(twostim_ccg_df.index):\n",
    "        try:\n",
    "            date = index[0]\n",
    "            owl = index[1]\n",
    "            channel_ot = index[-2]\n",
    "            channel_fore = index[-1]\n",
    "\n",
    "            channel_ot_srf = srf_df.loc[date, owl, channel_ot]\n",
    "            channel_fore_srf = srf_df.loc[date, owl, channel_fore]\n",
    "\n",
    "            srf_corr = np.corrcoef(\n",
    "                channel_ot_srf[\"norm_resp\"], channel_fore_srf[\"norm_resp\"]\n",
    "            )[0][1]\n",
    "\n",
    "            tmp = {\n",
    "                \"date\": date,\n",
    "                \"owl\": owl,\n",
    "                \"channel1\": channel_ot,\n",
    "                \"channel2\": channel_fore,\n",
    "                \"srf_corr\": srf_corr,\n",
    "            }\n",
    "            srf_corr_df.append(tmp)\n",
    "        except:\n",
    "            pass\n",
    "    srf_corr_df = pd.DataFrame(srf_corr_df)\n",
    "    srf_corr_df = srf_corr_df.set_index([\"date\", \"owl\", \"channel1\", \"channel2\"])\n",
    "    srf_corr_df = srf_corr_df.filter(\n",
    "        [\n",
    "            \"srf_corr\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "    merge = twostim_ccg_df.join(srf_corr_df, how=\"inner\", rsuffix=\"_srf\")\n",
    "\n",
    "    merge = merge.join(\n",
    "        single_ccg_df,\n",
    "        on=[\"date\", \"owl\", \"channel1\", \"channel2\"],\n",
    "        how=\"inner\",\n",
    "        rsuffix=\"_single\",\n",
    "    )\n",
    "\n",
    "    merge = merge.loc[\n",
    "        (merge[\"fixedintensity\"] == merge[\"intensity\"])\n",
    "        & (merge[\"fixedazi\"] == merge[\"azimuth\"])\n",
    "        & (merge[\"fixedele\"] == merge[\"elevation\"])\n",
    "    ]\n",
    "    print(\"Number of units:\", np.unique(merge.index).size)\n",
    "    print(\n",
    "        \"Number of sessions:\",\n",
    "        len(sorted(set([idx[:2] for idx in merge.index.values]))),\n",
    "    )\n",
    "    merge[\"csi\"] = (merge[\"xcorr_peak\"] - merge[\"xcorr_peak_single\"]) / (\n",
    "        merge[\"xcorr_peak\"] + merge[\"xcorr_peak_single\"]\n",
    "    )\n",
    "\n",
    "    quantiles = np.quantile(merge[\"srf_corr\"], [0.25, 0.50, 0.75])\n",
    "\n",
    "    low_corr = merge[(merge[\"srf_corr\"] <= quantiles[0])].set_index(\"relative_level\")\n",
    "    lower_mid_corr = merge[\n",
    "        (merge[\"srf_corr\"] > quantiles[0]) & ((merge[\"srf_corr\"] <= quantiles[1]))\n",
    "    ].set_index(\"relative_level\")\n",
    "    higher_mid_corr = merge[\n",
    "        (merge[\"srf_corr\"] > quantiles[1]) & ((merge[\"srf_corr\"] <= quantiles[2]))\n",
    "    ].set_index(\"relative_level\")\n",
    "    high_corr = merge[(merge[\"srf_corr\"] > quantiles[2])].set_index(\"relative_level\")\n",
    "\n",
    "    srf_corr_groups = [low_corr, lower_mid_corr, higher_mid_corr, high_corr]\n",
    "    srf_corr_labels = [\"<25 %\", \"25-50th %\", \"50-75th %\", \">75th %\"]\n",
    "\n",
    "    relative_levels = np.unique(high_corr.index)\n",
    "\n",
    "    groupdata = {group: [] for group in srf_corr_labels}\n",
    "    stats_df = []\n",
    "    for ndata, data in enumerate(srf_corr_groups):\n",
    "        plot_values = data.loc[[-15, -10], \"csi\"].values\n",
    "        groupdata[srf_corr_labels[ndata]] = plot_values\n",
    "        for val in plot_values:\n",
    "            tmp = {\"csi\": val, \"tuning_index\": srf_corr_labels[ndata]}\n",
    "            stats_df.append(tmp)\n",
    "    stats_df = pd.DataFrame(stats_df).set_index(\"tuning_index\")\n",
    "\n",
    "    plot_data = [values for values in groupdata.values()]\n",
    "\n",
    "    bp = ax.boxplot(\n",
    "        plot_data,\n",
    "        patch_artist=True,\n",
    "        notch=False,\n",
    "        showfliers=False,\n",
    "        whis=(5, 95),\n",
    "    )\n",
    "\n",
    "    ax.set_xticklabels(srf_corr_labels, fontsize=9)\n",
    "    # stats = scipy.stats.f_oneway(*plot_data)\n",
    "    stats = anova_tukey(stats_df, val_col=\"csi\", group_col=\"tuning_index\")\n",
    "    colors = [\"slategrey\", \"royalblue\", \"blue\", \"navy\"]\n",
    "    for patch, color in zip(bp[\"boxes\"], colors):\n",
    "        patch.set_facecolor(color)\n",
    "\n",
    "    for median in bp[\"medians\"]:\n",
    "        median.set(color=\"white\", linewidth=2)\n",
    "\n",
    "    if True:  # show_individual_data\n",
    "        width = 0.7\n",
    "        for k, groupvalues in enumerate(plot_data):\n",
    "            pos = k + 1\n",
    "            y = np.sort(groupvalues)\n",
    "            vals_x = np.linspace(\n",
    "                pos - width / 3,\n",
    "                pos + width / 3,\n",
    "                y.size,\n",
    "                endpoint=True,\n",
    "            )\n",
    "            ax.plot(\n",
    "                vals_x,\n",
    "                y,\n",
    "                color=\"k\",\n",
    "                ls=\"None\",\n",
    "                marker=\".\",\n",
    "                markersize=2,\n",
    "                markeredgewidth=0.0,\n",
    "                zorder=5,\n",
    "            )\n",
    "\n",
    "    if True:  # show_brackets\n",
    "        brackets_shrink = 1.0\n",
    "        for (idx_left, idx_right), y in brackets.items():\n",
    "            plot_bracket(\n",
    "                ax,\n",
    "                left=idx_left + 1,\n",
    "                right=idx_right + 1,\n",
    "                text=f\"{stats['tukey'].iloc[idx_left, idx_right]:.3f}\"[1:],\n",
    "                y=y,\n",
    "                shrink=brackets_shrink,\n",
    "            )\n",
    "\n",
    "    ax.spines[\"top\"].set_visible(False)\n",
    "    ax.spines[\"right\"].set_visible(False)\n",
    "\n",
    "    if anova_align == \"right\":\n",
    "        ax.text(\n",
    "            1,\n",
    "            1,\n",
    "            f\"p = {stats['anova_p']:.3g}\",\n",
    "            horizontalalignment=\"right\",\n",
    "            verticalalignment=\"top\",\n",
    "            fontsize=10,\n",
    "            transform=ax.transAxes,\n",
    "        )\n",
    "    else:\n",
    "        ax.text(\n",
    "            0,\n",
    "            1,\n",
    "            f\"p = {stats['anova_p']:.3g}\",\n",
    "            horizontalalignment=\"left\",\n",
    "            verticalalignment=\"top\",\n",
    "            fontsize=10,\n",
    "            transform=ax.transAxes\n",
    "            + mpl.transforms.ScaledTranslation(\n",
    "                +2 / 72,\n",
    "                0 / 72,\n",
    "                ax.figure.dpi_scale_trans,\n",
    "            ),\n",
    "        )\n",
    "\n",
    "    # ax.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, pos: f\"{x:.0%}\"))\n",
    "\n",
    "    ax.axhline(0, color=\"k\", ls=\"--\", zorder=-1)\n",
    "    ax.set_xlabel(\"Spatial Tuning Quartile\")\n",
    "    ax.set_ylabel(\"CSI\")\n",
    "    ax.set_ylim(-1, 1.45)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dd19bd8",
   "metadata": {},
   "source": [
    "## Gamma heat map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "c416bd48",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gamma_heatmap_scatter_resubmit(\n",
    "    single_srf_df, twostim_ccg_df, twostim_gamma_power_df, axs=None, caxs=None\n",
    "):\n",
    "    srf_df = single_srf_df.copy()\n",
    "    twostim_ccg_df = twostim_ccg_df.copy()\n",
    "    gamma_df = twostim_gamma_power_df.copy()\n",
    "    if axs == None:\n",
    "        fig, axs = plt.subplots(1, 3)\n",
    "    # axs.flatten()\n",
    "\n",
    "    twostim_ccg_df = twostim_ccg_df[twostim_ccg_df[\"corr_type\"] == \"cross_region\"]\n",
    "    srf_corr_df = []\n",
    "\n",
    "    for index in np.unique(twostim_ccg_df.index):\n",
    "        try:\n",
    "            date = index[0]\n",
    "            owl = index[1]\n",
    "            channel_ot = index[-2]\n",
    "            channel_fore = index[-1]\n",
    "            channel_ot_srf = srf_df.loc[date, owl, channel_ot]\n",
    "            channel_fore_srf = srf_df.loc[date, owl, channel_fore]\n",
    "            srf_corr = np.corrcoef(\n",
    "                channel_ot_srf[\"norm_resp\"], channel_fore_srf[\"norm_resp\"]\n",
    "            )[0][1]\n",
    "\n",
    "            tmp = {\n",
    "                \"date\": date,\n",
    "                \"owl\": owl,\n",
    "                \"channel1\": channel_ot,\n",
    "                \"channel2\": channel_fore,\n",
    "                \"srf_corr\": srf_corr,\n",
    "            }\n",
    "            srf_corr_df.append(tmp)\n",
    "\n",
    "        except:\n",
    "            pass\n",
    "\n",
    "    srf_corr_df = pd.DataFrame(srf_corr_df)\n",
    "    srf_corr_df = srf_corr_df.set_index([\"date\", \"owl\", \"channel1\", \"channel2\"])\n",
    "    srf_corr_df = srf_corr_df.filter(\n",
    "        [\n",
    "            \"srf_corr\",\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "    merge = twostim_ccg_df.join(srf_corr_df, how=\"inner\", rsuffix=\"_srf\")\n",
    "\n",
    "    merge = merge.join(\n",
    "        gamma_df,\n",
    "        on=[\"date\", \"owl\", \"channel1\"],\n",
    "        how=\"inner\",\n",
    "        rsuffix=\"1\",\n",
    "    )\n",
    "\n",
    "    merge = merge.join(\n",
    "        gamma_df,\n",
    "        on=[\"date\", \"owl\", \"channel2\"],\n",
    "        how=\"inner\",\n",
    "        rsuffix=\"2\",\n",
    "    )\n",
    "\n",
    "    merge = merge.loc[\n",
    "        (merge[\"fixedintensity\"] == merge[\"fixedintensity1\"])\n",
    "        & (merge[\"fixedintensity\"] == merge[\"fixedintensity2\"])\n",
    "        & (merge[\"fixedintensity1\"] == merge[\"fixedintensity2\"])\n",
    "        & (merge[\"varyingintensity\"] == merge[\"varyingintensity1\"])\n",
    "        & (merge[\"varyingintensity\"] == merge[\"varyingintensity2\"])\n",
    "        & (merge[\"varyingintensity1\"] == merge[\"varyingintensity2\"])\n",
    "        & (merge[\"fixedazi\"] == merge[\"fixedazi1\"])\n",
    "        & (merge[\"fixedazi\"] == merge[\"fixedazi2\"])\n",
    "        & (merge[\"fixedazi1\"] == merge[\"fixedazi2\"])\n",
    "        & (merge[\"fixedele\"] == merge[\"fixedele1\"])\n",
    "        & (merge[\"fixedele\"] == merge[\"fixedele2\"])\n",
    "        & (merge[\"fixedele1\"] == merge[\"fixedele2\"])\n",
    "    ]\n",
    "\n",
    "    merge[\"mean_gamma_power\"] = (merge[\"gammapower\"] + merge[\"gammapower2\"]) / 2\n",
    "    merge[\"mean_gamma_plv\"] = (merge[\"gamma_plv\"] + merge[\"gamma_plv2\"]) / 2\n",
    "    merge[\"diff_gamma_angle\"] = np.abs(\n",
    "        merge[\"gamma_plv_angle\"] - merge[\"gamma_plv_angle2\"]\n",
    "    )\n",
    "\n",
    "    quantiles = np.quantile(merge[\"srf_corr\"], [0.25, 0.50, 0.75])\n",
    "    low_corr = merge[(merge[\"srf_corr\"] <= quantiles[0])].set_index(\"relative_level\")\n",
    "    lower_mid_corr = merge[\n",
    "        (merge[\"srf_corr\"] > quantiles[0]) & ((merge[\"srf_corr\"] <= quantiles[1]))\n",
    "    ].set_index(\"relative_level\")\n",
    "\n",
    "    higher_mid_corr = merge[\n",
    "        (merge[\"srf_corr\"] > quantiles[1]) & ((merge[\"srf_corr\"] <= quantiles[2]))\n",
    "    ].set_index(\"relative_level\")\n",
    "\n",
    "    high_corr = merge[(merge[\"srf_corr\"] > quantiles[2])].set_index(\"relative_level\")\n",
    "\n",
    "    merge = merge.set_index(\"relative_level\")\n",
    "\n",
    "    psth_len = (merge.iloc[0][\"ccg\"].size + 1) / 2\n",
    "\n",
    "    lags = scipy.signal.correlation_lags(psth_len, psth_len) / 1000\n",
    "    lags_mask = np.abs(lags) <= 0.95\n",
    "    lags_mask_std = np.abs(lags) >= 0.825\n",
    "    relative_levels = np.unique(merge.index)\n",
    "    groupdata = [low_corr, lower_mid_corr, higher_mid_corr, high_corr]\n",
    "\n",
    "    gamma_xcorr_df = {\n",
    "        \"mean_gamma_power\": [],\n",
    "        \"gamma_power_ot\": [],\n",
    "        \"gamma_power_forebrain\": [],\n",
    "        \"xregion_corr\": [],\n",
    "        \"relative_level\": [],\n",
    "    }\n",
    "\n",
    "    gamma_heatmap = np.empty([len(groupdata), len(relative_levels)])\n",
    "    gamma_heatmap_ot = np.empty([len(groupdata), len(relative_levels)])\n",
    "    gamma_heatmap_fb = np.empty([len(groupdata), len(relative_levels)])\n",
    "    xregion_synchrony = np.empty([len(groupdata), len(relative_levels)])\n",
    "\n",
    "    for ndata, data in enumerate(groupdata):\n",
    "        for k, relative_level in enumerate(relative_levels):\n",
    "            mean_ccg = np.mean(np.vstack(data.loc[relative_level, \"ccg\"]), axis=0)[::-1]\n",
    "            mean_condition_gamma = np.mean(\n",
    "                np.vstack(data.loc[relative_level, \"mean_gamma_power\"])\n",
    "            )\n",
    "            mean_gamma_ot = np.mean(np.vstack(data.loc[relative_level, \"gammapower\"]))\n",
    "            mean_gamma_fb = np.mean(np.vstack(data.loc[relative_level, \"gammapower2\"]))\n",
    "            mean_srf_corr = np.mean(np.vstack(data.loc[relative_level, \"srf_corr\"]))\n",
    "            peak_properties = get_peak(mean_ccg, lags)\n",
    "            gamma_heatmap[ndata, k] = mean_condition_gamma\n",
    "            gamma_heatmap_ot[ndata, k] = mean_gamma_ot\n",
    "            gamma_heatmap_fb[ndata, k] = mean_gamma_fb\n",
    "            xregion_synchrony[ndata, k] = peak_properties[\"peak_corr\"]\n",
    "            gamma_xcorr_df[\"mean_gamma_power\"].append(mean_condition_gamma)\n",
    "            gamma_xcorr_df[\"gamma_power_ot\"].append(mean_gamma_ot)\n",
    "            gamma_xcorr_df[\"gamma_power_forebrain\"].append(mean_gamma_fb)\n",
    "            gamma_xcorr_df[\"xregion_corr\"].append(peak_properties[\"peak_corr\"])\n",
    "            gamma_xcorr_df[\"relative_level\"].append(relative_level)\n",
    "\n",
    "    ot_corr_val = np.corrcoef(\n",
    "        gamma_xcorr_df[\"gamma_power_ot\"], gamma_xcorr_df[\"xregion_corr\"]\n",
    "    )[0][1]\n",
    "\n",
    "    nrt_corr_val = np.corrcoef(\n",
    "        gamma_xcorr_df[\"gamma_power_forebrain\"], gamma_xcorr_df[\"xregion_corr\"]\n",
    "    )[0][1]\n",
    "\n",
    "    xregion_corr_val = np.corrcoef(\n",
    "        gamma_xcorr_df[\"mean_gamma_power\"], gamma_xcorr_df[\"xregion_corr\"]\n",
    "    )[0][1]\n",
    "\n",
    "    sc3 = axs[3].scatter(\n",
    "        gamma_xcorr_df[\"gamma_power_ot\"],\n",
    "        gamma_xcorr_df[\"xregion_corr\"],\n",
    "        label=f\"R = {ot_corr_val:.2f}\",\n",
    "        color=\"k\",\n",
    "    )\n",
    "\n",
    "    axs[3].set_title(\"Within OT\")\n",
    "    axs[3].text(\n",
    "        0.05,\n",
    "        0.92,\n",
    "        f\"R = {ot_corr_val:.2f}\",\n",
    "        ha=\"left\",\n",
    "        va=\"top\",\n",
    "        fontsize=12,\n",
    "        transform=axs[3].transAxes,\n",
    "    )\n",
    "    axs[3].set_xlabel(\"Mean OT gamma power\")\n",
    "    axs[3].set_ylabel(\"X-Region Synchrony\")\n",
    "    axs[3].spines[\"top\"].set_visible(False)\n",
    "    axs[3].spines[\"right\"].set_visible(False)\n",
    "\n",
    "    sc4 = axs[4].scatter(\n",
    "        gamma_xcorr_df[\"gamma_power_forebrain\"],\n",
    "        gamma_xcorr_df[\"xregion_corr\"],\n",
    "        label=f\"R = {nrt_corr_val:.2f}\",\n",
    "        color=\"k\",\n",
    "    )\n",
    "\n",
    "    axs[4].set_title(\"Within nRt\")\n",
    "    axs[4].text(\n",
    "        0.05,\n",
    "        0.92,\n",
    "        f\"R = {nrt_corr_val:.2f}\",\n",
    "        ha=\"left\",\n",
    "        va=\"top\",\n",
    "        fontsize=12,\n",
    "        transform=axs[4].transAxes,\n",
    "    )\n",
    "    axs[4].set_xlabel(\"Mean nRt gamma power\")\n",
    "    axs[4].set_ylabel(\"X-Region Synchrony\")\n",
    "    axs[4].spines[\"top\"].set_visible(False)\n",
    "    axs[4].spines[\"right\"].set_visible(False)\n",
    "\n",
    "    sc5 = axs[5].scatter(\n",
    "        gamma_xcorr_df[\"mean_gamma_power\"],\n",
    "        gamma_xcorr_df[\"xregion_corr\"],\n",
    "        label=f\"R = {xregion_corr_val:.2f}\",\n",
    "        color=\"k\",\n",
    "    )\n",
    "    axs[5].set_title(\"Cross Region\")\n",
    "    axs[5].text(\n",
    "        0.05,\n",
    "        0.92,\n",
    "        f\"R = {xregion_corr_val:.2f}\",\n",
    "        ha=\"left\",\n",
    "        va=\"top\",\n",
    "        fontsize=12,\n",
    "        transform=axs[5].transAxes,\n",
    "    )\n",
    "    axs[5].set_xlabel(\"Mean X-region gamma power\")\n",
    "    axs[5].set_ylabel(\"X-Region Synchrony\")\n",
    "    axs[5].spines[\"top\"].set_visible(False)\n",
    "    axs[5].spines[\"right\"].set_visible(False)\n",
    "\n",
    "    # plt.figure()\n",
    "    #    hm = axs.imshow(gamma_heatmap, cmap = 'viridis')\n",
    "\n",
    "    hm = axs[0].imshow(\n",
    "        gamma_heatmap_ot,\n",
    "        cmap=\"viridis\",\n",
    "        vmin=np.min(gamma_heatmap_ot),\n",
    "        vmax=np.max(gamma_heatmap_ot),\n",
    "    )\n",
    "\n",
    "    axs[0].set_xticks(np.arange(len(relative_levels)), relative_levels)\n",
    "    axs[0].set_yticks(\n",
    "        np.arange(len(groupdata)), [\"<25th %\", \"25-50th %\", \"50-75th %\", \">75th %\"]\n",
    "    )\n",
    "    cb = plt.colorbar(hm, cax=caxs[0]).set_label(\"Mean OT gamma power\", size=10)\n",
    "\n",
    "    axs[0].set_xlabel(\"Relative level [dB]\")\n",
    "    hm2 = axs[1].imshow(\n",
    "        gamma_heatmap_fb,\n",
    "        cmap=\"viridis\",\n",
    "        vmin=np.min(gamma_heatmap_ot),\n",
    "        vmax=np.max(gamma_heatmap_ot),\n",
    "    )\n",
    "\n",
    "    axs[1].set_xticks(np.arange(len(relative_levels)), relative_levels)\n",
    "    axs[1].set_yticks(\n",
    "        np.arange(len(groupdata)), [\"<25th %\", \"25-50th %\", \"50-75th %\", \">75th %\"]\n",
    "    )\n",
    "\n",
    "    cb2 = plt.colorbar(hm2, cax=caxs[1]).set_label(\"Mean nRt gamma power\", size=10)\n",
    "\n",
    "    axs[1].set_xlabel(\"Relative level [dB]\")\n",
    "\n",
    "    hm3 = axs[2].imshow(\n",
    "        gamma_heatmap,\n",
    "        cmap=\"viridis\",\n",
    "        vmin=np.min(gamma_heatmap_ot),\n",
    "        vmax=np.max(gamma_heatmap_ot),\n",
    "    )\n",
    "\n",
    "    axs[2].set_xticks(np.arange(len(relative_levels)), relative_levels)\n",
    "    axs[2].set_yticks(\n",
    "        np.arange(len(groupdata)), [\"<25th %\", \"25-50th %\", \"50-75th %\", \">75th %\"]\n",
    "    )\n",
    "    cb2 = plt.colorbar(hm2, cax=caxs[2]).set_label(\"Mean X-region gamma power\", size=9)\n",
    "    axs[2].set_xlabel(\"Relative level [dB]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "210a60da",
   "metadata": {},
   "source": [
    "## Spike Field Coherence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "dfa3f23f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def within_area_sfc_competition(within_area_sfc_df, axs=None, colors=c_ot):\n",
    "    if axs is None:\n",
    "        fig, axs = plt.subplots(1, 3)\n",
    "\n",
    "    axs.flatten()\n",
    "\n",
    "    lfp_df = within_area_sfc_df.copy()\n",
    "\n",
    "    lfp_df = lfp_df.set_index(\"relative_level\")\n",
    "\n",
    "    frequency = lfp_df.iloc[0][\"frequency\"]\n",
    "\n",
    "    level_plots = [-15, 0, 10]\n",
    "    level_plot_label = [\"-15\", \"0\", \"+10\"]\n",
    "\n",
    "    for kk, level in enumerate(level_plots):\n",
    "\n",
    "        sub = np.mean(np.vstack(lfp_df.loc[level, \"spike_field_coherence\"]), axis=0)\n",
    "        sub_sem = scipy.stats.sem(\n",
    "            np.vstack(lfp_df.loc[level, \"spike_field_coherence\"]), axis=0\n",
    "        )\n",
    "        shuffle = np.mean(np.vstack(lfp_df.loc[level, \"trial_shuffled_sfc\"]), axis=0)\n",
    "        shuffle_sem = scipy.stats.sem(\n",
    "            np.vstack(lfp_df.loc[level, \"trial_shuffled_sfc\"]), axis=0\n",
    "        )\n",
    "\n",
    "        axs[kk].set_title(level_plot_label[kk])\n",
    "        axs[kk].errorbar(frequency, sub, sub_sem, color=colors[0])\n",
    "\n",
    "        shuff = axs[kk].errorbar(\n",
    "            frequency, shuffle, shuffle_sem, color=\"0.4\", label=\"trial shuffled\"\n",
    "        )\n",
    "\n",
    "        axs[kk].set_xlim(0, 100)\n",
    "        if kk == 0:\n",
    "            axs[kk].set_xticks([0, 25, 50, 75], [0, 25, 50, 75])\n",
    "            axs[kk].set_xlabel(\"Frequency (Hz)\")\n",
    "            axs[kk].set_ylabel(\"Spike Field\\nCoherence\")\n",
    "        else:\n",
    "            axs[kk].set_xticks([0, 25, 50, 75], [])\n",
    "            axs[kk].set_yticklabels([])\n",
    "        axs[kk].set_ylim(bottom=0.05, top=0.450)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbb268fc",
   "metadata": {},
   "source": [
    "## Preferred Phase Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "1cae425f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def figure_withinregion_spike_phase(\n",
    "    cross_region_df,\n",
    "    colors=\"b\",\n",
    "    axs=None,\n",
    "):\n",
    "    if axs is None:\n",
    "        fig, axs = plt.subplots(1, 3, sharey=True)\n",
    "\n",
    "    gamma_df = cross_region_df.copy()\n",
    "\n",
    "    print(\"Number of units:\", np.unique(gamma_df.index).size)\n",
    "    print(\n",
    "        \"Number of sessions:\",\n",
    "        len(sorted(set([idx[:2] for idx in gamma_df.index.values]))),\n",
    "    )\n",
    "\n",
    "    gamma_df = gamma_df.set_index(\"relative_level\")\n",
    "    relative_levels = np.unique(gamma_df.index)\n",
    "\n",
    "    # colors = ['slategrey', 'royalblue','blue','navy']\n",
    "\n",
    "    phase_bins = np.linspace(-np.pi, np.pi, 9)\n",
    "\n",
    "    levels_plot = [-15, 0, 10]\n",
    "    for k, relative_level in enumerate(levels_plot):\n",
    "        spike_angles = gamma_df.loc[relative_level, \"gamma_spike_angles\"].values\n",
    "        mean_phases = gamma_df.loc[relative_level, \"gamma_plv_angle\"].values\n",
    "\n",
    "        mean_vs = np.mean(gamma_df.loc[relative_level, \"gamma_plv\"].values)\n",
    "        norm_trace = []\n",
    "        for angles, mean_ang in zip(spike_angles, mean_phases):\n",
    "            transform = angles - mean_ang\n",
    "            transform = wrap_to_pi(transform)\n",
    "            zz = np.histogram(transform, phase_bins)\n",
    "            norm_zz = zz[0] / sum(zz[0])\n",
    "            norm_trace.append(norm_zz)\n",
    "        mean_trace = np.mean(norm_trace, axis=0)\n",
    "        sem_trace = scipy.stats.sem(norm_trace, axis=0)\n",
    "\n",
    "        spike_angles_shuff = gamma_df.loc[\n",
    "            relative_level, \"gamma_spike_angles_shuffled\"\n",
    "        ].values\n",
    "        mean_phases_shuff = gamma_df.loc[\n",
    "            relative_level, \"gamma_plv_angle_shuffled\"\n",
    "        ].values\n",
    "        mean_vs_shuff = np.mean(\n",
    "            gamma_df.loc[relative_level, \"gamma_plv_shuffled\"].values\n",
    "        )\n",
    "        norm_trace_shuff = []\n",
    "        for angles, mean_ang in zip(spike_angles_shuff, mean_phases_shuff):\n",
    "            transform = angles - mean_ang\n",
    "            transform = wrap_to_pi(transform)\n",
    "            zz = np.histogram(transform, phase_bins)\n",
    "            norm_zz = zz[0] / sum(zz[0])\n",
    "            norm_trace_shuff.append(norm_zz)\n",
    "        mean_trace_shuff = np.mean(norm_trace_shuff, axis=0)\n",
    "        sem_trace_shuff = scipy.stats.sem(norm_trace_shuff, axis=0)\n",
    "\n",
    "        # axs[ndata, k].plot(bins[:-1], histogram, color=colors[ndata], ls=\"-\", lw=1)\n",
    "        # print(histogram)\n",
    "        axs[k].errorbar(\n",
    "            phase_bins[:-1],\n",
    "            mean_trace_shuff,\n",
    "            sem_trace_shuff,\n",
    "            color=\"0.4\",\n",
    "            ls=\"-\",\n",
    "            lw=1,\n",
    "        )\n",
    "        axs[k].errorbar(\n",
    "            phase_bins[:-1], mean_trace, sem_trace, color=colors, ls=\"-\", lw=1\n",
    "        )\n",
    "\n",
    "        axs[k].text(\n",
    "            0.5,\n",
    "            0.80,\n",
    "            f\"VS {mean_vs:.3f}\",\n",
    "            transform=axs[k].transAxes\n",
    "            + mpl.transforms.ScaledTranslation(\n",
    "                0 / 72,\n",
    "                +4 / 72,\n",
    "                axs[k].figure.dpi_scale_trans,\n",
    "            ),\n",
    "            ha=\"center\",\n",
    "            fontsize=10,\n",
    "        )\n",
    "        axs[k].set_ylim(bottom=0.0, top=0.50)\n",
    "        axs[k].spines[\"top\"].set_visible(False)\n",
    "        axs[k].spines[\"right\"].set_visible(False)\n",
    "        axs[k].set_xlim(-np.pi, +np.pi)\n",
    "        axs[k].set_xticks(\n",
    "            [-np.pi, 0, +np.pi], labels=[\"-π\", \"0\", \"π\"] if k == 0 else []\n",
    "        )\n",
    "        axs[k].set_xlabel(\"φ - Best φ\")\n",
    "        if k != 0:\n",
    "            axs[k].set_yticks([])\n",
    "            axs[k].set_xlabel(None)\n",
    "\n",
    "    colors = [\"r\", \"orange\", \"k\", \"blue\", \"purple\", \"violet\"]\n",
    "    for k, relative_level in enumerate(relative_levels):\n",
    "        spike_angles = gamma_df.loc[relative_level, \"gamma_spike_angles\"].values\n",
    "        mean_phases = gamma_df.loc[relative_level, \"gamma_plv_angle\"].values\n",
    "        prop = phase_properties(mean_phases)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46e861b7",
   "metadata": {},
   "source": [
    "# Figures "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41eda98e",
   "metadata": {},
   "source": [
    "## Figure 1 Experimental Design"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e44f87de",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(3.30, 2.4 + 2.10 + 2.4 + 1.25))\n",
    "\n",
    "ax1a = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.05,\n",
    "    left=0.8,\n",
    "    width=2.1,\n",
    "    height=2.4,\n",
    "    label=\"A\",\n",
    ")\n",
    "plt.setp(ax1a, frame_on=False, xticks=[], yticks=[], zorder=20)\n",
    "owl_image = matplotlib.image.imread(\n",
    "    r\".\\other_figures\\experimental_design_forebrainmanuscript.png\"\n",
    ")\n",
    "ax1a.imshow(owl_image)\n",
    "\n",
    "subplot_indicator(ax1a, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "ax1b = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3 + 2.4 + 0.15,\n",
    "    left=0.8,\n",
    "    width=1.9,\n",
    "    height=2.0,\n",
    "    label=\"B\",\n",
    ")\n",
    "\n",
    "ax1b_cb = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3 + 2.4 + 0.15,\n",
    "    left=0.8 + 1.9 + 0.05,\n",
    "    width=0.05,\n",
    "    height=2.0,\n",
    ")\n",
    "\n",
    "ax1b.set_title(\"Example OT unit\")\n",
    "\n",
    "figure_spatial_receptive_field(midhighcorr_srf_ot, ax=ax1b, cax=ax1b_cb)\n",
    "subplot_indicator(ax1b, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "ax1c = figure_add_axes_inch(\n",
    "    fig,\n",
    "    bottom=0.5,\n",
    "    left=0.8,\n",
    "    width=1.9,\n",
    "    height=2.0,\n",
    "    label=\"C\",\n",
    ")\n",
    "\n",
    "ax1c_cb = figure_add_axes_inch(\n",
    "    fig,\n",
    "    bottom=0.5,\n",
    "    left=0.8 + 1.9 + 0.05,\n",
    "    width=0.05,\n",
    "    height=2.0,\n",
    ")\n",
    "\n",
    "ax1c.set_title(\"Example nRt unit\")\n",
    "\n",
    "figure_spatial_receptive_field(midhighcorr_srf_fb, ax=ax1c, cax=ax1c_cb)\n",
    "subplot_indicator(ax1c, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "fig.savefig(OUTDIR / \"figure1.pdf\")\n",
    "fig.savefig(OUTDIR / \"figure1.png\")\n",
    "fig.savefig(OUTDIR / \"figure1.eps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f04b754",
   "metadata": {},
   "source": [
    "## Figure 2 Azimuth and Elevation tuning nRt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddd5a4a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "single_srf_df = single_srf[single_srf[\"region\"] == \"Forebrain\"]\n",
    "single_rlf_df = single_rlf[single_rlf[\"region\"] == \"Forebrain\"]\n",
    "\n",
    "fig = plt.figure(figsize=(6.80, 6.40))\n",
    "ax2a = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3,\n",
    "    left=0.6,\n",
    "    width=2.1,\n",
    "    height=3.0,\n",
    "    label=\"A\",\n",
    ")\n",
    "ax2a.set_title(\"nRt Azimuth Tuning\")\n",
    "\n",
    "ax2a_cb = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3,\n",
    "    left=0.6 + 2.1 + 0.05,\n",
    "    width=0.05,\n",
    "    height=3.0,\n",
    ")\n",
    "\n",
    "ax2b = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3 + 3.0 + 0.6,\n",
    "    left=0.7,\n",
    "    width=2.1,\n",
    "    height=2.0,\n",
    "    label=\"B\",\n",
    ")\n",
    "\n",
    "azimuth_heatplot_alt(single_srf_df, axs=[ax2a, ax2b], caxs=ax2a_cb)\n",
    "\n",
    "ax2c = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3,\n",
    "    left=0.5 + 2.4 + 1.2,\n",
    "    width=2.1,\n",
    "    height=3.0,\n",
    "    label=\"C\",\n",
    ")\n",
    "\n",
    "ax2c_cb = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3,\n",
    "    left=0.5 + 2.4 + 2.1 + 1.2 + 0.05,\n",
    "    width=0.05,\n",
    "    height=3.0,\n",
    ")\n",
    "ax2c.set_title(\"nRt Elevation Tuning\")\n",
    "\n",
    "ax2d = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3 + 3.0 + 0.6,\n",
    "    left=0.6 + 2.4 + 1.2,\n",
    "    width=2.1,\n",
    "    height=2.0,\n",
    "    label=\"D\",\n",
    ")\n",
    "\n",
    "elevation_heatplot_alt(single_srf_df, axs=[ax2c, ax2d], caxs=ax2c_cb)\n",
    "\n",
    "subplot_indicator(ax2a, ha=\"left\", va=\"bottom\", pad_inch=0.6)\n",
    "subplot_indicator(ax2b, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(ax2c, ha=\"left\", va=\"bottom\", pad_inch=0.6)\n",
    "subplot_indicator(ax2d, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "# fig.outline()\n",
    "\n",
    "fig.savefig(OUTDIR / \"figure2.pdf\")\n",
    "fig.savefig(OUTDIR / \"figure2.png\")\n",
    "fig.savefig(OUTDIR / \"figure2.eps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33c994c5",
   "metadata": {},
   "source": [
    "## Figure 3 nRt flat noise competition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0f99670",
   "metadata": {},
   "outputs": [],
   "source": [
    "## forebrain specific data\n",
    "forebrain_twostim = twostim_rlf[(twostim_rlf[\"region\"] == \"Forebrain\")]\n",
    "forebrain_singlestim = single_rlf[(single_rlf[\"region\"] == \"Forebrain\")]\n",
    "\n",
    "twostim_ccg_forebrain = twostim_ccg[(twostim_ccg[\"corr_type\"] == \"within_Forebrain\")]\n",
    "single_ccg_forebrain = single_ccg[(single_ccg[\"corr_type\"] == \"within_Forebrain\")]\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(3.30, 5.10 + 2.6 + 0.1))\n",
    "ax3a = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.1,\n",
    "    left=0.8,\n",
    "    width=2.4,\n",
    "    height=1.5,\n",
    "    label=\"A\",\n",
    ")\n",
    "figure_competition_boxplot(\n",
    "    forebrain_twostim,\n",
    "    forebrain_singlestim,\n",
    "    ax=ax3a,\n",
    "    brackets={(0, 5): 0.9, (1, 5): 0.8, (2, 5): 0.7},\n",
    "    colors=c_forebrain,\n",
    ")\n",
    "ax3a.set_ylabel(\"Change in\\nSpike Rate\")\n",
    "ax3a.yaxis.set_label_coords(-0.4 / 2.4, 0.5)\n",
    "subplot_indicator(ax3a, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "axs3b, axg3b = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=2,\n",
    "    group_top=2.4,\n",
    "    group_left=0.8,\n",
    "    individual_width=(2.4 - 0.1) / 2,\n",
    "    individual_height=2.0,\n",
    "    wspace=0.1,\n",
    "    hspace=0.0,\n",
    ")\n",
    "\n",
    "figure_coincident_rasterplot(\n",
    "    spiketrain_data=example_spiketrains_forebrain_flat,\n",
    "    axs=axs3b,\n",
    ")\n",
    "axs3b[0, 0].set_title(\"-15 dB\")\n",
    "axs3b[0, 1].set_title(\"+10 dB\")\n",
    "subplot_indicator(axg3b, \"B\", ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "axs3c, axg3c = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=6,\n",
    "    group_top=2.2 + 2.6 + 0.1,\n",
    "    group_left=0.8,\n",
    "    individual_width=2.4 / 6,\n",
    "    individual_height=0.6,\n",
    "    wspace=0.0,\n",
    "    hspace=0.0,\n",
    ")\n",
    "\n",
    "figure_xcorr_ccg(\n",
    "    twostim_ccg_forebrain,\n",
    "    single_ccg_forebrain,\n",
    "    axs=axs3c.flatten()[:],\n",
    "    colors=c_forebrain,\n",
    ")\n",
    "\n",
    "\n",
    "for ax in axs3c[0, :]:\n",
    "    ax.set_ylim(top=1.0e-5)\n",
    "axs3c[0, 0].set_yticks(np.array([0.5]) * 1e-5, [\"0.5\"], minor=False)\n",
    "axs3c[0, 0].yaxis.set_label_coords(-0.3 / (2.4 / 6), 0.5)\n",
    "\n",
    "subplot_indicator(axg3c, \"C\", ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "ax3d = figure_add_axes_inch(\n",
    "    fig,\n",
    "    left=0.8,\n",
    "    width=2.4,\n",
    "    bottom=0.5,\n",
    "    height=1.5,\n",
    "    label=\"D\",\n",
    ")\n",
    "figure_xcorr_boxplot(\n",
    "    twostim_ccg_forebrain,\n",
    "    single_ccg_forebrain,\n",
    "    ax=ax3d,\n",
    "    colors=c_forebrain,\n",
    "    brackets={},\n",
    "    # brackets={(1,5): 0.9, (2,5):0.8},\n",
    ")\n",
    "ax3d.yaxis.set_label_coords(-0.3 / 2.4, 0.5)\n",
    "subplot_indicator(ax3d, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "ax3d.set_ylim(-0.75, 0.75)\n",
    "\n",
    "# condition_batch(fig, 0.8, 0.3, text=\"flat noise\", fontsize=10, color='b', y_pt=2, x_pt=0, pad_pt=2)\n",
    "condition_batch(\n",
    "    fig,\n",
    "    left=0.8 + 0.2,\n",
    "    top=2.3,\n",
    "    text=\"flat noise\",\n",
    "    color=c_forebrain[0],\n",
    "    fontsize=12,\n",
    "    ha=\"right\",\n",
    ")\n",
    "\n",
    "fig.savefig(OUTDIR / \"figure3.pdf\")\n",
    "fig.savefig(OUTDIR / \"figure3.png\")\n",
    "fig.savefig(OUTDIR / \"figure3.eps\")\n",
    "\n",
    "# fig.outline()\n",
    "fig.get_size_inches()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5073f42",
   "metadata": {},
   "source": [
    "## Figure 4 Tuning Selectivity Index + Competition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da0df2e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "forebrain_twostim = twostim_rlf[twostim_rlf[\"region\"] == \"Forebrain\"]\n",
    "forebrain_singlestim = single_rlf[single_rlf[\"region\"] == \"Forebrain\"]\n",
    "forebrain_singlestim_out = single_rlf_out[single_rlf_out[\"region\"] == \"Forebrain\"]\n",
    "\n",
    "\n",
    "forebrain_single_ccg = single_ccg[single_ccg[\"corr_type\"] == \"within_Forebrain\"]\n",
    "forebrain_ts_ccg = twostim_ccg[twostim_ccg[\"corr_type\"] == \"within_Forebrain\"]\n",
    "\n",
    "fig = plt.figure(figsize=(6.80, 4.5))\n",
    "ax4a = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.1,\n",
    "    left=0.8,\n",
    "    width=5.6,\n",
    "    height=1.0,\n",
    "    label=\"A\",\n",
    ")\n",
    "\n",
    "\n",
    "histogram_selectivity_resubmit(\n",
    "    forebrain_twostim,\n",
    "    forebrain_singlestim,\n",
    "    forebrain_singlestim_out,\n",
    "    colors=c_forebrain,\n",
    "    ax=ax4a,\n",
    ")\n",
    "\n",
    "subplot_indicator(ax4a, ha=\"left\", va=\"top\", pad_inch=0.75)\n",
    "# subplot_indicator(ax4a, ha=\"left\", va=\"top\", pad_inch=0.75)\n",
    "\n",
    "ax4b = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=1.5 + 0.3 + 0.2,\n",
    "    left=0.8,\n",
    "    width=2.4,\n",
    "    height=2.0,\n",
    "    label=\"B\",\n",
    ")\n",
    "\n",
    "ax4c = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=1.5 + 0.3 + 0.2,\n",
    "    left=0.8 + 2.4 + 0.8,\n",
    "    width=2.4,\n",
    "    height=2.0,\n",
    "    label=\"C\",\n",
    ")\n",
    "\n",
    "figure_singlesound_rlf_selectivity_resubmit(\n",
    "    forebrain_singlestim, forebrain_singlestim_out, ax4b\n",
    ")\n",
    "subplot_indicator(ax4b, ha=\"left\", va=\"top\", pad_inch=0.75)\n",
    "\n",
    "# ax2e = figure_add_axes_inch(\n",
    "#     fig,\n",
    "#     bottom = 0.50,\n",
    "#     left=0.8,\n",
    "#     width=5.6,\n",
    "#     height=2.0,\n",
    "#     label=\"D\",\n",
    "# )\n",
    "\n",
    "figure_rlf_selectivity_resubmit(\n",
    "    forebrain_twostim,\n",
    "    forebrain_singlestim,\n",
    "    forebrain_singlestim_out,\n",
    "    colors=c_forebrain,\n",
    "    ax=ax4c,\n",
    ")\n",
    "ax4c.set_ylim(-0.75, 0.75)\n",
    "ax4c.yaxis.set_label_coords(-0.4 / 2.4, 0.5)\n",
    "ax4c.set_title(\"Flat Noise Competition\")\n",
    "subplot_indicator(ax4c, ha=\"left\", va=\"top\", pad_inch=0.75)\n",
    "\n",
    "# fig.outline()\n",
    "fig.savefig(OUTDIR / \"figure4.pdf\")\n",
    "fig.savefig(OUTDIR / \"figure4.png\")\n",
    "fig.savefig(OUTDIR / \"figure4.eps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e75aa70c",
   "metadata": {},
   "source": [
    "## Figure 5 Modeling Figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0558613",
   "metadata": {},
   "outputs": [],
   "source": [
    "forebrain_singlestim = single_rlf[single_rlf[\"region\"] == \"Forebrain\"]\n",
    "forebrain_singlestim_out = single_rlf_out[single_rlf_out[\"region\"] == \"Forebrain\"]\n",
    "\n",
    "sound_level, rate_driver_4b, rate_competitor_4b = rlf_modeling_data(\n",
    "    forebrain_singlestim, forebrain_singlestim_out\n",
    ")\n",
    "\n",
    "forebrain_twostim = twostim_rlf[twostim_rlf[\"region\"] == \"Forebrain\"]\n",
    "forebrain_singlestim = single_rlf[single_rlf[\"region\"] == \"Forebrain\"]\n",
    "forebrain_singlestim_out = single_rlf_out[single_rlf_out[\"region\"] == \"Forebrain\"]\n",
    "\n",
    "\n",
    "relative_level, relative_response_d, relative_response_c = competition_modeling_data(\n",
    "    forebrain_twostim, forebrain_singlestim, forebrain_singlestim_out\n",
    ")\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(6.80, 4.90))\n",
    "\n",
    "ax5a = figure_add_axes_inch(fig, top=0.3, left=0.8, width=2.4, height=1.5, label=\"A\")\n",
    "ax5b = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3,\n",
    "    left=0.8 + 2.4 + 1.1,\n",
    "    width=2.4,\n",
    "    height=1.5,\n",
    "    label=\"B\",\n",
    ")\n",
    "\n",
    "model(relative_response_d, axs=[ax5a, ax5b])\n",
    "ax5b.set_title(\"Invariant nRt Response\")\n",
    "ax5a.set_title(\"OT responses\")\n",
    "\n",
    "ax5c = figure_add_axes_inch(\n",
    "    fig, top=0.3 + 1.5 + 1.0, left=0.8, width=2.4, height=1.5, label=\"C\"\n",
    ")\n",
    "ax5d = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3 + 1.5 + 1.0,\n",
    "    left=0.8 + 2.4 + 1.1,\n",
    "    width=2.4,\n",
    "    height=1.5,\n",
    "    label=\"D\",\n",
    ")\n",
    "model(relative_response_c, axs=[ax5c, ax5d])\n",
    "ax5d.set_title(\"Increasing nRt Response\")\n",
    "ax5c.set_title(\"OT responses\")\n",
    "\n",
    "\n",
    "subplot_indicator(ax5a, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(ax5b, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(ax5c, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(ax5d, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "fig.savefig(OUTDIR / \"figure5.pdf\")\n",
    "fig.savefig(OUTDIR / \"figure5.png\")\n",
    "fig.savefig(OUTDIR / \"figure5.eps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38ee029b",
   "metadata": {},
   "source": [
    "## Figure 6 Extreme Selectivity Resubmission"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ac47e1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "forebrain_twostim = twostim_rlf[twostim_rlf[\"region\"] == \"Forebrain\"]\n",
    "forebrain_singlestim = single_rlf[single_rlf[\"region\"] == \"Forebrain\"]\n",
    "forebrain_singlestim_out = single_rlf_out[single_rlf_out[\"region\"] == \"Forebrain\"]\n",
    "\n",
    "fig = plt.figure(figsize=(6.80, 2.5))\n",
    "ax6a = figure_add_axes_inch(fig, top=0.5, left=0.8, width=2.4, height=1.5, label=\"A\")\n",
    "ax6a.set_title(\"Extreme Frontal\\nDriver Selectivity\")\n",
    "ax6b = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.5,\n",
    "    left=0.8 + 2.4 + 1.1,\n",
    "    width=2.4,\n",
    "    height=1.5,\n",
    "    label=\"B\",\n",
    ")\n",
    "ax6b.set_title(\"Extreme Lateral\\nDriver Selectivity\")\n",
    "\n",
    "extreme_selectivity_competition_resubmit(\n",
    "    forebrain_twostim, forebrain_singlestim, forebrain_singlestim_out, axs=[ax6a, ax6b]\n",
    ")\n",
    "# fig.outline()\n",
    "subplot_indicator(ax6a, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(ax6b, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "\n",
    "fig.savefig(OUTDIR / \"figure6.pdf\")\n",
    "fig.savefig(OUTDIR / \"figure6.png\")\n",
    "fig.savefig(OUTDIR / \"figure6.eps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4a3aff6",
   "metadata": {},
   "source": [
    "## Figure 7 Frontal OT vs. Peripheral OT competition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03e680de",
   "metadata": {},
   "outputs": [],
   "source": [
    "frontal_ot_twostim = twostim_rlf[\n",
    "    (twostim_rlf[\"region\"] == \"OT\") & (abs(twostim_rlf[\"fixedazi\"]) <= 45)\n",
    "]\n",
    "frontal_ot_singlestim = single_rlf[\n",
    "    (single_rlf[\"region\"] == \"OT\") & (abs(single_rlf[\"azimuth\"]) <= 45)\n",
    "]\n",
    "peripheral_ot_twostim = twostim_rlf[\n",
    "    (twostim_rlf[\"region\"] == \"OT\") & (abs(twostim_rlf[\"fixedazi\"]) > 45)\n",
    "]\n",
    "peripheral_ot_singlestim = single_rlf[\n",
    "    (single_rlf[\"region\"] == \"OT\") & (abs(single_rlf[\"azimuth\"]) > 45)\n",
    "]\n",
    "\n",
    "fig = plt.figure(figsize=(6.80, 2.30))\n",
    "ax7a = figure_add_axes_inch(fig, top=0.3, left=0.9, width=2.4, height=1.5, label=\"A\")\n",
    "ax7b = figure_add_axes_inch(\n",
    "    fig, top=0.3, left=0.8 + 2.4 + 1.1, width=2.4, height=1.5, label=\"B\"\n",
    ")\n",
    "\n",
    "figure_competition_boxplot(\n",
    "    frontal_ot_twostim, frontal_ot_singlestim, colors=[\"b\"] * 6, ax=ax7a\n",
    ")\n",
    "ax7a.set_title(\"Frontal OT\")\n",
    "\n",
    "figure_competition_boxplot(\n",
    "    peripheral_ot_twostim, peripheral_ot_singlestim, colors=[\"b\"] * 6, ax=ax7b\n",
    ")\n",
    "ax7b.set_title(\"Lateral OT\")\n",
    "\n",
    "subplot_indicator(ax7a, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(ax7b, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "\n",
    "fig.savefig(OUTDIR / \"figure7.pdf\")\n",
    "fig.savefig(OUTDIR / \"figure7.png\")\n",
    "fig.savefig(OUTDIR / \"figure7.eps\")\n",
    "\n",
    "# fig.outline()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab864822",
   "metadata": {},
   "source": [
    "## Figure 8 AM noise competition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c98be217",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Forebrain Driver 55Hz, Competitor 75Hz\n",
    "am_single_ccg_55 = am_single_ccg[\n",
    "    (am_single_ccg[\"corr_type\"] == \"within_Forebrain\")\n",
    "    & (am_single_ccg[\"modulation_frequency\"] == 55)\n",
    "]\n",
    "am_single_rlf_55 = am_single_rlf[\n",
    "    (am_single_rlf[\"region\"] == \"Forebrain\")\n",
    "    & (am_single_rlf[\"modulation_frequency\"] == 55)\n",
    "]\n",
    "\n",
    "am_twostim_ccg_55 = am_twostim_ccg[\n",
    "    (am_twostim_ccg[\"fixed_modulation_frequency\"] == 55)\n",
    "    & (am_twostim_ccg[\"corr_type\"] == \"within_Forebrain\")\n",
    "]\n",
    "am_twostim_rlf_55 = am_twostim_rlf[\n",
    "    (am_twostim_rlf[\"fixed_modulation_frequency\"] == 55)\n",
    "    & (am_twostim_rlf[\"region\"] == \"Forebrain\")\n",
    "]\n",
    "\n",
    "## Forebrain Driver 75Hz, Competitor 55Hz\n",
    "am_single_ccg_75 = am_single_ccg[\n",
    "    (am_single_ccg[\"corr_type\"] == \"within_Forebrain\")\n",
    "    & (am_single_ccg[\"modulation_frequency\"] == 75)\n",
    "]\n",
    "am_single_rlf_75 = am_single_rlf[\n",
    "    (am_single_rlf[\"region\"] == \"Forebrain\")\n",
    "    & (am_single_rlf[\"modulation_frequency\"] == 75)\n",
    "]\n",
    "\n",
    "am_twostim_ccg_75 = am_twostim_ccg[\n",
    "    (am_twostim_ccg[\"fixed_modulation_frequency\"] == 75)\n",
    "    & (am_twostim_ccg[\"corr_type\"] == \"within_Forebrain\")\n",
    "]\n",
    "am_twostim_rlf_75 = am_twostim_rlf[\n",
    "    (am_twostim_rlf[\"fixed_modulation_frequency\"] == 75)\n",
    "    & (am_twostim_rlf[\"region\"] == \"Forebrain\")\n",
    "]\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(6.80, 5.10 + 2.6 + 0.1))\n",
    "\n",
    "ax8a = figure_add_axes_inch(fig, top=0.1, left=0.8, width=2.4, height=1.5, label=\"A\")\n",
    "figure_competition_boxplot(\n",
    "    am_twostim_rlf_55,\n",
    "    am_single_rlf_55,\n",
    "    ax=ax8a,\n",
    "    colors=c_forebrain_am55,\n",
    "    brackets={(0, 5): 0.9, (1, 5): 0.8, (2, 5): 0.7},\n",
    ")\n",
    "ax8a.set_ylabel(\"Change in\\nSpike Rate\")\n",
    "ax8a.yaxis.set_label_coords(-0.4 / 2.4, 0.5)\n",
    "subplot_indicator(ax8a, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "ax8b = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.1,\n",
    "    left=0.8 + 2.4 + 0.3 + 0.8,\n",
    "    width=2.4,\n",
    "    height=1.5,\n",
    "    label=\"B\",\n",
    ")\n",
    "figure_competition_boxplot(\n",
    "    am_twostim_rlf_75,\n",
    "    am_single_rlf_75,\n",
    "    ax=ax8b,\n",
    "    colors=c_forebrain_am75,\n",
    "    brackets={(0, 5): 0.9, (1, 5): 0.8, (2, 5): 0.7},\n",
    ")\n",
    "\n",
    "ax8b.set_ylabel(\"Change in\\nSpike Rate\")\n",
    "ax8b.yaxis.set_label_coords(-0.4 / 2.4, 0.5)\n",
    "subplot_indicator(ax8b, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "axs8c, axg8c = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=2,\n",
    "    group_top=2.4,\n",
    "    group_left=0.8,\n",
    "    individual_width=(2.4 - 0.1) / 2,\n",
    "    individual_height=2.0,\n",
    "    wspace=0.1,\n",
    "    hspace=0.0,\n",
    ")\n",
    "\n",
    "figure_coincident_rasterplot(\n",
    "    spiketrain_data=example_spiketrains_am55,\n",
    "    axs=axs8c,\n",
    ")\n",
    "axs8c[0, 0].set_title(\"-15 dB\")\n",
    "axs8c[0, 1].set_title(\"+10 dB\")\n",
    "subplot_indicator(axg8c, \"C\", ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "axs8d, axg8d = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=2,\n",
    "    group_top=2.4,\n",
    "    group_left=0.8 + 2.4 + 0.3 + 0.8,\n",
    "    individual_width=(2.4 - 0.1) / 2,\n",
    "    individual_height=2.0,\n",
    "    wspace=0.1,\n",
    "    hspace=0.0,\n",
    ")\n",
    "\n",
    "figure_coincident_rasterplot(\n",
    "    spiketrain_data=example_spiketrains_am75,\n",
    "    axs=axs8d,\n",
    ")\n",
    "axs8d[0, 0].set_title(\"-15 dB\")\n",
    "axs8d[0, 1].set_title(\"+10 dB\")\n",
    "subplot_indicator(axg8d, \"D\", ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "axs8e, axg8e = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=6,\n",
    "    group_top=2.2 + 2.6 + 0.1,\n",
    "    group_left=0.8,\n",
    "    individual_width=2.4 / 6,\n",
    "    individual_height=0.6,\n",
    "    wspace=0.0,\n",
    "    hspace=0.0,\n",
    ")\n",
    "figure_xcorr_ccg(\n",
    "    am_twostim_ccg_55, am_single_ccg_55, axs=axs8e.flatten()[:], colors=c_forebrain_am55\n",
    ")\n",
    "for ax in axs8e[0, :]:\n",
    "    ax.set_ylim(top=1.5e-5)\n",
    "axs8e[0, 0].set_yticks(np.array([0, 1]) * 1e-5, minor=False)\n",
    "axs8e[0, 0].set_yticks(np.array([0.5]) * 1e-5, minor=True)\n",
    "axs8e[0, 0].yaxis.set_label_coords(-0.3 / (2.4 / 6), 0.5)\n",
    "subplot_indicator(axg8e, \"E\", ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "axs8f, axg8f = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=6,\n",
    "    group_top=2.2 + 2.6 + 0.1,\n",
    "    group_left=0.8 + 2.4 + 0.3 + 0.8,\n",
    "    individual_width=2.4 / 6,\n",
    "    individual_height=0.6,\n",
    "    wspace=0.0,\n",
    "    hspace=0.0,\n",
    ")\n",
    "figure_xcorr_ccg(\n",
    "    am_twostim_ccg_75, am_single_ccg_75, axs=axs8f.flatten()[:], colors=c_forebrain_am75\n",
    ")\n",
    "for ax in axs8f[0, :]:\n",
    "    ax.set_ylim(top=1.5e-5)\n",
    "axs8f[0, 0].set_yticks(np.array([0, 1]) * 1e-5, minor=False)\n",
    "axs8f[0, 0].set_yticks(np.array([0.5]) * 1e-5, minor=True)\n",
    "axs8f[0, 0].yaxis.set_label_coords(-0.3 / (2.4 / 6), 0.5)\n",
    "subplot_indicator(axg8f, \"F\", ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "ax8g = figure_add_axes_inch(\n",
    "    fig,\n",
    "    left=0.8,\n",
    "    width=2.4,\n",
    "    bottom=0.5,\n",
    "    height=1.5,\n",
    "    label=\"G\",\n",
    ")\n",
    "figure_xcorr_boxplot(\n",
    "    am_twostim_ccg_55, am_single_ccg_55, ax=ax8g, colors=c_forebrain_am55, brackets={}\n",
    ")\n",
    "\n",
    "ax8g.yaxis.set_label_coords(-0.3 / 2.4, 0.5)\n",
    "subplot_indicator(ax8g, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "ax8h = figure_add_axes_inch(\n",
    "    fig,\n",
    "    left=0.8 + 2.4 + 0.3 + 0.8,\n",
    "    width=2.4,\n",
    "    bottom=0.5,\n",
    "    height=1.5,\n",
    "    label=\"H\",\n",
    ")\n",
    "figure_xcorr_boxplot(\n",
    "    am_twostim_ccg_75, am_single_ccg_75, ax=ax8h, colors=c_forebrain_am75, brackets={}\n",
    ")\n",
    "ax8h.yaxis.set_label_coords(-0.3 / 2.4, 0.5)\n",
    "subplot_indicator(ax8h, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "condition_batch(\n",
    "    fig,\n",
    "    left=0.8,\n",
    "    top=2.3,\n",
    "    text=\"Driver\\n55 Hz\",\n",
    "    color=c_forebrain_am55[0],\n",
    "    fontsize=12,\n",
    "    ha=\"right\",\n",
    ")\n",
    "condition_batch(\n",
    "    fig,\n",
    "    left=0.8 + 2.4 + 0.3 + 0.8,\n",
    "    top=2.3,\n",
    "    text=\"Driver\\n75 Hz\",\n",
    "    color=c_forebrain_am75[0],\n",
    "    fontsize=12,\n",
    "    ha=\"right\",\n",
    ")\n",
    "\n",
    "# fig.savefig(OUTDIR / \"figure3.pdf\")\n",
    "# fig.savefig(OUTDIR / \"forebrain_am_competition.png\")\n",
    "# fig.savefig(OUTDIR / \"figure3.eps\")\n",
    "\n",
    "# fig.outline()\n",
    "print(\n",
    "    f\"Figure size: {' x '.join(fig.get_size_inches().astype('str'))} inches == {' x '.join(f'{l*2.54:.2f}' for l in fig.get_size_inches())} cm\"\n",
    ")\n",
    "\n",
    "fig.savefig(OUTDIR / \"figure8.pdf\")\n",
    "fig.savefig(OUTDIR / \"figure8.png\")\n",
    "fig.savefig(OUTDIR / \"figure8.eps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c3ab54a",
   "metadata": {},
   "source": [
    "## Figure 9 AM noise Switch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79fe4269",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Forebrain Driver 55Hz, Competitor 75Hz\n",
    "am_single_ccg_55 = am_single_ccg[\n",
    "    (am_single_ccg[\"corr_type\"] == \"within_Forebrain\")\n",
    "    & (am_single_ccg[\"modulation_frequency\"] == 55)\n",
    "]\n",
    "am_single_rlf_55 = am_single_rlf[\n",
    "    (am_single_rlf[\"region\"] == \"Forebrain\")\n",
    "    & (am_single_rlf[\"modulation_frequency\"] == 55)\n",
    "]\n",
    "\n",
    "am_twostim_ccg_55 = am_twostim_ccg[\n",
    "    (am_twostim_ccg[\"fixed_modulation_frequency\"] == 55)\n",
    "    & (am_twostim_ccg[\"corr_type\"] == \"within_Forebrain\")\n",
    "]\n",
    "am_twostim_rlf_55 = am_twostim_rlf[\n",
    "    (am_twostim_rlf[\"fixed_modulation_frequency\"] == 55)\n",
    "    & (am_twostim_rlf[\"region\"] == \"Forebrain\")\n",
    "]\n",
    "\n",
    "## Forebrain Driver 75Hz, Competitor 55Hz\n",
    "am_single_ccg_75 = am_single_ccg[\n",
    "    (am_single_ccg[\"corr_type\"] == \"within_Forebrain\")\n",
    "    & (am_single_ccg[\"modulation_frequency\"] == 75)\n",
    "]\n",
    "am_single_rlf_75 = am_single_rlf[\n",
    "    (am_single_rlf[\"region\"] == \"Forebrain\")\n",
    "    & (am_single_rlf[\"modulation_frequency\"] == 75)\n",
    "]\n",
    "\n",
    "am_twostim_ccg_75 = am_twostim_ccg[\n",
    "    (am_twostim_ccg[\"fixed_modulation_frequency\"] == 75)\n",
    "    & (am_twostim_ccg[\"corr_type\"] == \"within_Forebrain\")\n",
    "]\n",
    "am_twostim_rlf_75 = am_twostim_rlf[\n",
    "    (am_twostim_rlf[\"fixed_modulation_frequency\"] == 75)\n",
    "    & (am_twostim_rlf[\"region\"] == \"Forebrain\")\n",
    "]\n",
    "\n",
    "\n",
    "df_twostim_55 = am_twostim_stim_phaselocking[\n",
    "    (am_twostim_stim_phaselocking[\"fixed_modulation_frequency\"] == 55)\n",
    "    & (am_twostim_stim_phaselocking[\"region\"] == \"Forebrain\")\n",
    "]\n",
    "df_twostim_75 = am_twostim_stim_phaselocking[\n",
    "    (am_twostim_stim_phaselocking[\"fixed_modulation_frequency\"] == 75)\n",
    "    & (am_twostim_stim_phaselocking[\"region\"] == \"Forebrain\")\n",
    "]\n",
    "\n",
    "fig = plt.figure(figsize=(6.80, 5.6))\n",
    "\n",
    "axs9a, axg9a = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=6,\n",
    "    group_top=0.3 + 0.15,\n",
    "    group_left=0.8,\n",
    "    individual_width=2.4 / 6,\n",
    "    individual_height=0.6,\n",
    "    wspace=0.030,\n",
    "    hspace=0.0,\n",
    ")\n",
    "figure_xcorr_ccg(\n",
    "    am_twostim_ccg_55, am_single_ccg_55, axs=axs9a, colors=c_forebrain_am55\n",
    ")\n",
    "for ax in axs9a[0, :]:\n",
    "    ax.set_ylim(bottom=-0.10e-5, top=1.0e-5)\n",
    "\n",
    "axs9a[0, 0].set_yticks(np.array([0.5]) * 1e-5, [\".5\"], minor=False)\n",
    "# axs3a[0, 0].set_yticks(np.array([0.5]) * 1e-5, minor=True)\n",
    "axs9a[0, 0].yaxis.set_label_coords(-0.3 / (2.4 / 6), 0.5)\n",
    "subplot_indicator(axg9a, \"A\", ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "axs9b, axg9b = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=6,\n",
    "    group_top=0.3 + 1.0 + 0.15,\n",
    "    group_left=0.8,\n",
    "    individual_width=2.4 / 6,\n",
    "    individual_height=0.6,\n",
    "    wspace=0.030,\n",
    "    hspace=0.0,\n",
    ")\n",
    "figure_xcorr_pspectra(\n",
    "    am_twostim_ccg_55, am_single_ccg_55, axs=axs9b, colors=c_forebrain_am55\n",
    ")\n",
    "subplot_indicator(axg9b, \"B\", ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "axs9c, axg9c = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=6,\n",
    "    group_top=0.4 + 1.0 + 0.15 + 0.6 + 0.4,\n",
    "    group_left=0.8,  #  + 2.4 + 0.3 + 0.8\n",
    "    individual_width=2.4 / 6,\n",
    "    individual_height=0.6,\n",
    "    wspace=0.030,\n",
    "    hspace=0.0,\n",
    ")\n",
    "\n",
    "figure_stimphase_driver_competitor(\n",
    "    df_twostim_55, axs=axs9c[:], colors=[c_forebrain_am55, c_forebrain_am75]\n",
    ")\n",
    "\n",
    "\n",
    "for ax in axs9c[0, :]:\n",
    "    ax.set_title(\"\")\n",
    "axg9c.set_title(\"\")\n",
    "axg9c.set_ylabel(\"VS Competitor\")\n",
    "axg9c.yaxis.set_label_coords(-0.3 / 2.4, 0.5)\n",
    "subplot_indicator(axg9c, \"C\", ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "\n",
    "axs9c[0, 0].set_xlim(0, 0.5)\n",
    "axs9c[0, 0].set_ylim(0, 0.5)\n",
    "\n",
    "ax9d = figure_add_axes_inch(\n",
    "    fig,\n",
    "    bottom=0.5,\n",
    "    left=0.8,\n",
    "    width=2.55,\n",
    "    height=1.5,\n",
    "    label=\"D\",\n",
    ")\n",
    "figure_stimphase_switch(\n",
    "    df_twostim_55, ax=ax9d, colors=[c_forebrain_am55, c_forebrain_am75]\n",
    ")\n",
    "# ax3d.set_ylabel(\"Change in\\nSpike Rate\")\n",
    "# ax3d.yaxis.set_label_coords(-0.4 / 2.4, 0.5)\n",
    "subplot_indicator(ax9d, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "ax9d.set_ylim(0, 0.2)\n",
    "\n",
    "axs9e, axg9e = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=6,\n",
    "    group_top=0.3 + 0.15,\n",
    "    group_left=0.8 + 2.4 + 0.1 + 0.8,\n",
    "    individual_width=2.4 / 6,\n",
    "    individual_height=0.6,\n",
    "    wspace=0.020,\n",
    "    hspace=0.0,\n",
    ")\n",
    "figure_xcorr_ccg(\n",
    "    am_twostim_ccg_75, am_single_ccg_75, axs=axs9e.flatten()[:], colors=c_forebrain_am75\n",
    ")\n",
    "for ax in axs9e[0, :]:\n",
    "    ax.set_ylim(bottom=-0.10e-5, top=1.0e-5)\n",
    "axs9e[0, 0].set_yticks(np.array([0.5]) * 1e-5, [\".5\"], minor=False)\n",
    "# axs3e[0, 0].set_yticks(np.array([0.5]) * 1e-5, minor=True)\n",
    "axs9e[0, 0].yaxis.set_label_coords(-0.3 / (2.4 / 6), 0.5)\n",
    "subplot_indicator(axg9e, \"E\", ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "axs9f, axg9f = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=6,\n",
    "    group_top=0.3 + 1.0 + 0.15,\n",
    "    group_left=0.8 + 2.4 + 0.1 + 0.8,\n",
    "    individual_width=2.4 / 6,\n",
    "    individual_height=0.6,\n",
    "    wspace=0.020,\n",
    "    hspace=0.0,\n",
    ")\n",
    "figure_xcorr_pspectra(\n",
    "    am_twostim_ccg_75, am_single_ccg_75, axs=axs9f, colors=c_forebrain_am75\n",
    ")\n",
    "subplot_indicator(axg9f, \"F\", ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "axs9g, axg9g = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=6,\n",
    "    group_top=0.4 + 1.0 + 0.15 + 0.6 + 0.4,\n",
    "    group_left=0.8 + 2.4 + 0.1 + 0.8,\n",
    "    individual_width=2.4 / 6,\n",
    "    individual_height=0.6,\n",
    "    wspace=0.030,\n",
    "    hspace=0.0,\n",
    ")\n",
    "\n",
    "figure_stimphase_driver_competitor(\n",
    "    df_twostim_75, axs=axs9g[:], colors=[c_forebrain_am75, c_forebrain_am55]\n",
    ")\n",
    "\n",
    "\n",
    "for ax in axs9g[0, :]:\n",
    "    ax.set_title(\"\")\n",
    "axg9g.set_title(\"\")\n",
    "axg9g.set_ylabel(\"VS Competitor\")\n",
    "axg9g.yaxis.set_label_coords(-0.3 / 2.4, 0.5)\n",
    "subplot_indicator(axg9g, \"G\", ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "\n",
    "axs9g[0, 0].set_xlim(0, 0.5)\n",
    "axs9g[0, 0].set_ylim(0, 0.5)\n",
    "\n",
    "\n",
    "ax9h = figure_add_axes_inch(\n",
    "    fig,\n",
    "    bottom=0.5,\n",
    "    left=0.8 + 2.4 + 0.1 + 0.8,\n",
    "    width=2.55,\n",
    "    height=1.5,\n",
    "    label=\"H\",\n",
    ")\n",
    "figure_stimphase_switch(\n",
    "    df_twostim_75, ax=ax9h, colors=[c_forebrain_am75, c_forebrain_am55]\n",
    ")\n",
    "# ax3d.set_ylabel(\"Change in\\nSpike Rate\")\n",
    "# ax3d.yaxis.set_label_coords(-0.4 / 2.4, 0.5)\n",
    "subplot_indicator(ax9h, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "ax9h.set_ylim(0, 0.2)\n",
    "\n",
    "\n",
    "axg9a.set_title(\"Driver 55Hz\", pad=22, fontsize=12)\n",
    "axg9e.set_title(\"Driver 75Hz\", pad=22, fontsize=12)\n",
    "# fig.savefig(OUTDIR / \"forebrain_xcorr_pspectra.png\")\n",
    "# fig.outline()\n",
    "# condition_batch(\n",
    "#     fig,\n",
    "#     left=0.8,\n",
    "#     top=2.3,\n",
    "#     text=\"Driver\\n55 Hz\",\n",
    "#     color=c_forebrain_am55[0],\n",
    "#     fontsize=12,\n",
    "#     ha=\"right\",\n",
    "# )\n",
    "\n",
    "fig.savefig(OUTDIR / \"figure9.pdf\")\n",
    "fig.savefig(OUTDIR / \"figure9.png\")\n",
    "fig.savefig(OUTDIR / \"figure9.eps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2515bf8f",
   "metadata": {},
   "source": [
    "## Figure 10 Spatial tuning Correlation distribution "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "411d325d",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6.80, 10.2))\n",
    "ax10a = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.1,\n",
    "    left=0.8,\n",
    "    width=5.8,\n",
    "    height=1.5,\n",
    "    label=\"A\",\n",
    ")\n",
    "srf_correlation_distribution(single_srf, twostim_ccg, ax10a)\n",
    "subplot_indicator(ax10a, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "ax10b = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3 + 1.5 + 0.7,\n",
    "    left=0.8,\n",
    "    width=2,\n",
    "    height=2,\n",
    "    label=\"B\",\n",
    ")\n",
    "\n",
    "ax10b_cb = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3 + 1.5 + 0.7,\n",
    "    left=0.8 + 2 + 0.05,\n",
    "    width=0.05,\n",
    "    height=2,\n",
    ")\n",
    "\n",
    "figure_spatial_receptive_field(lowcorr_srf_ot, ax=ax10b, cax=ax10b_cb)\n",
    "subplot_indicator(ax10b, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "ax10b.set_title(\"OT\")\n",
    "\n",
    "ax10c = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3 + 1.5 + 0.7,\n",
    "    left=0.8 + 1.4 + 2.0,\n",
    "    width=2,\n",
    "    height=2,\n",
    "    label=\"C\",\n",
    ")\n",
    "\n",
    "ax10c_cb = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3 + 1.5 + 0.7,\n",
    "    left=0.8 + 1.4 + 2.0 + 2 + 0.05,\n",
    "    width=0.05,\n",
    "    height=2,\n",
    ")\n",
    "\n",
    "figure_spatial_receptive_field(lowcorr_srf_fb, ax=ax10c, cax=ax10c_cb)\n",
    "subplot_indicator(ax10c, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "ax10c.set_title(\"nRt\")\n",
    "\n",
    "ax10d = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.03 + 2.0 + 1.0 + 1.5 + 0.7,\n",
    "    left=0.8,\n",
    "    width=2.0,\n",
    "    height=2.0,\n",
    "    label=\"D\",\n",
    ")\n",
    "\n",
    "ax10d_cb = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.03 + 2.0 + 1.0 + 1.5 + 0.7,\n",
    "    left=0.8 + 2.0 + 0.05,\n",
    "    width=0.05,\n",
    "    height=2.0,\n",
    ")\n",
    "\n",
    "figure_spatial_receptive_field(highcorr_srf_ot, ax=ax10d, cax=ax10d_cb)\n",
    "subplot_indicator(ax10d, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "ax10e = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.03 + 2.0 + 1.0 + 1.5 + 0.7,\n",
    "    left=0.8 + 1.4 + 2.0,\n",
    "    width=2.0,\n",
    "    height=2.0,\n",
    "    label=\"E\",\n",
    ")\n",
    "\n",
    "ax10e_cb = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.03 + 2.0 + 1.0 + 1.5 + 0.7,\n",
    "    left=0.8 + 1.4 + 2.0 + 2.0 + 0.05,\n",
    "    width=0.05,\n",
    "    height=2.0,\n",
    ")\n",
    "\n",
    "figure_spatial_receptive_field(highcorr_srf_fb, ax=ax10e, cax=ax10e_cb)\n",
    "subplot_indicator(ax10e, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "ax10f = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.03 + 2.0 + 1.0 + 1.5 + 0.7 + 2.0 + 0.80,\n",
    "    left=0.8,\n",
    "    width=5.8,\n",
    "    height=1.5,\n",
    "    label=\"F\",\n",
    ")\n",
    "ax10f.set_ylim(0, 20)\n",
    "# ax8g = figure_add_axes_inch(\n",
    "#    fig,\n",
    "#     top=0.03 + 2.0 + 1.0 + 1.5+ 0.7 + 2.0 + .80,\n",
    "#     left=0.8 + 1.4 + 2.0,\n",
    "#     width=2.0,\n",
    "#     height=2.0,\n",
    "#     label=\"G\",\n",
    "# )\n",
    "\n",
    "\n",
    "srf_correlation_frontalperiperpheral(single_srf, twostim_ccg, axs=ax10f)\n",
    "subplot_indicator(ax10f, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "# subplot_indicator(ax5b, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "\n",
    "\n",
    "# fig.outline()\n",
    "fig.savefig(OUTDIR / \"figure10.pdf\")\n",
    "fig.savefig(OUTDIR / \"figure10.png\")\n",
    "fig.savefig(OUTDIR / \"figure10.eps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69e36499",
   "metadata": {},
   "source": [
    "## Figure 11 Dual Region Xcorr by SRF Correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "075c7612",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6.80, 7.0))\n",
    "\n",
    "axs11a, axg11a = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=4,\n",
    "    ncols=6,\n",
    "    group_top=0.4,\n",
    "    group_left=0.8,\n",
    "    individual_width=(5.10 / 6) * 0.85,\n",
    "    individual_height=(5.0 / 6) * 0.85,\n",
    "    wspace=0.15,\n",
    "    hspace=0.45,\n",
    ")\n",
    "\n",
    "ax11a_cb = figure_add_axes_inch(fig, top=0.4, left=6.1, width=0.07, height=(5.0 * 0.85))\n",
    "\n",
    "fig_dual_region_xcorr_grid(single_srf, twostim_ccg, axs=axs11a, caxs=ax11a_cb)\n",
    "\n",
    "\n",
    "ax11b = figure_add_axes_inch(\n",
    "    fig, top=0.3 + 3.6 + 0.5 + 0.30 + 0.3, left=0.8, width=2.6, height=1.5, label=\"B\"\n",
    ")\n",
    "csi_boxplot_dualregion(\n",
    "    single_srf,\n",
    "    twostim_ccg,\n",
    "    single_ccg,\n",
    "    ax=ax11b,\n",
    "    anova_align=\"left\",\n",
    "    brackets={(0, 1): 0.62, (1, 2): 0.62, (1, 3): 0.72, (1, 4): 0.82, (1, 5): 0.92},\n",
    ")\n",
    "\n",
    "\n",
    "ax11c = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.3 + 3.6 + 0.5 + 0.3 + 0.3,\n",
    "    left=0.8 + 2.1 + 1.2,\n",
    "    width=2.6,\n",
    "    height=1.5,\n",
    "    label=\"C\",\n",
    ")\n",
    "csi_boxplot_dualregion_srfcorr(\n",
    "    single_srf,\n",
    "    twostim_ccg,\n",
    "    single_ccg,\n",
    "    ax=ax11c,\n",
    "    anova_align=\"left\",\n",
    "    brackets={(0, 3): 0.9, (0, 2): 0.8},\n",
    ")\n",
    "\n",
    "# subplot_indicator(ax3d, ha=\"left\", va=\"top\", pad_inch=0.7)\n",
    "subplot_indicator(axg11a, \"A\", ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(ax11b, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(ax11c, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "\n",
    "# fig.outline()\n",
    "\n",
    "fig.savefig(OUTDIR / \"figure11.pdf\")\n",
    "fig.savefig(OUTDIR / \"figure11.png\")\n",
    "fig.savefig(OUTDIR / \"figure11.eps\")\n",
    "# fig.outline()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91bb1fe3",
   "metadata": {},
   "source": [
    "## Figure 12 Within and Across Region Gamma Power"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fc0f1ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6.80, 7.3))\n",
    "\n",
    "ax12a = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.20,\n",
    "    left=0.8 + 0.3,\n",
    "    width=2.0,\n",
    "    height=2.1,\n",
    "    label=\"A\",\n",
    ")\n",
    "\n",
    "ax12a_cb = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.20 + 0.3,\n",
    "    left=0.8 + 2.0 + 0.05 + 0.3,\n",
    "    width=0.05,\n",
    "    height=1.5,\n",
    ")\n",
    "ax12a.set_title(\"Within OT\")\n",
    "ax12a.set_ylabel(\"Spatial Tuning\\nCorrelation Percentile\", x=-1, fontsize=10)\n",
    "\n",
    "ax12b = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.4 + 0.1 + 2.1,\n",
    "    left=0.8 + 0.3,\n",
    "    width=2.0,\n",
    "    height=2.1,\n",
    "    label=\"B\",\n",
    ")\n",
    "ax12b.set_title(\"Within nRt\")\n",
    "\n",
    "ax12b_cb = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.4 + 0.3 + 0.3 + 1.9,\n",
    "    left=0.8 + 2.0 + 0.05 + 0.3,\n",
    "    width=0.05,\n",
    "    height=1.5,\n",
    ")\n",
    "ax12b.set_ylabel(\"Spatial Tuning\\nCorrelation Percentile\", x=-1, fontsize=10)\n",
    "\n",
    "ax12c = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.60 + 0.1 + 2.1 + 0.1 + 2.1,\n",
    "    left=0.8 + 0.3,\n",
    "    width=2.0,\n",
    "    height=2.1,\n",
    "    label=\"C\",\n",
    ")\n",
    "ax12c.set_title(\"Cross Region\")\n",
    "\n",
    "ax12c_cb = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.60 + 0.3 + 0.3 + 1.9 + 0.3 + 1.9,\n",
    "    left=0.8 + 2.0 + 0.05 + 0.3,\n",
    "    width=0.05,\n",
    "    height=1.5,\n",
    ")\n",
    "ax12c.set_ylabel(\"Spatial Tuning\\nCorrelation Percentile\", x=-1, fontsize=10)\n",
    "\n",
    "\n",
    "ax12d = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.1 + 0.4,\n",
    "    left=0.8 + 2.1 + 0.05 + 1.5,\n",
    "    width=2.0,\n",
    "    height=1.5,\n",
    "    label=\"D\",\n",
    ")\n",
    "\n",
    "ax12e = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.5 + 0.1 + 2.3,\n",
    "    left=0.8 + 2.1 + 0.05 + 1.5,\n",
    "    width=2.0,\n",
    "    height=1.5,\n",
    "    label=\"E\",\n",
    ")\n",
    "\n",
    "ax12f = figure_add_axes_inch(\n",
    "    fig,\n",
    "    top=0.9 + 0.1 + 2.1 + 0.1 + 2.1,\n",
    "    left=0.8 + 2.1 + 0.05 + 1.5,\n",
    "    width=2.0,\n",
    "    height=1.5,\n",
    "    label=\"F\",\n",
    ")\n",
    "\n",
    "lowgamma = twostim_gamma_power[twostim_gamma_power[\"gamma_band\"] == \"lowgamma\"]\n",
    "gamma_heatmap_scatter_resubmit(\n",
    "    single_srf,\n",
    "    twostim_ccg,\n",
    "    lowgamma,\n",
    "    [ax12a, ax12b, ax12c, ax12d, ax12e, ax12f],\n",
    "    [ax12a_cb, ax12b_cb, ax12c_cb],\n",
    ")\n",
    "\n",
    "subplot_indicator(ax12a, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(ax12b, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(ax12c, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(ax12d, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(ax12e, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(ax12f, ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "# fig.outline()\n",
    "\n",
    "fig.savefig(OUTDIR / \"figure12.pdf\")\n",
    "fig.savefig(OUTDIR / \"figure12.png\")\n",
    "fig.savefig(OUTDIR / \"figure12.eps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0f2ea17",
   "metadata": {},
   "source": [
    "## Figure 13 Spike-LFP Coherence "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01fa04bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "## spike field coherence figure\n",
    "fig = plt.figure(figsize=(6.80, 5.7))\n",
    "\n",
    "axs13a, axg13a = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=3,\n",
    "    group_top=0.6,\n",
    "    group_left=0.8,\n",
    "    individual_width=2.4 / 3,\n",
    "    individual_height=1.0,\n",
    "    wspace=0.05,\n",
    "    hspace=0.0,\n",
    ")\n",
    "within_area_sfc_competition(\n",
    "    within_area_sfc[within_area_sfc[\"region\"] == \"OT\"], axs=axs13a[0, :], colors=c_ot\n",
    ")\n",
    "axg13a.set_title(\"OT\", y=1.3)\n",
    "\n",
    "\n",
    "axs13b, axg13b = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=3,\n",
    "    group_top=0.6,\n",
    "    group_left=0.8 + 2.4 + 0.8,\n",
    "    individual_width=2.4 / 3,\n",
    "    individual_height=1.0,\n",
    "    wspace=0.05,\n",
    "    hspace=0.0,\n",
    ")\n",
    "within_area_sfc_competition(\n",
    "    within_area_sfc[within_area_sfc[\"region\"] == \"Forebrain\"],\n",
    "    axs=axs13b[0, :],\n",
    "    colors=c_forebrain,\n",
    ")\n",
    "axs13b[0, 0].set_ylabel(None)\n",
    "axg13b.set_title(\"nRt\", y=1.3)\n",
    "\n",
    "axs13c, axg13c = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=3,\n",
    "    group_top=0.6 + 1.0 + 0.8,\n",
    "    group_left=0.8,\n",
    "    individual_width=2.4 / 3,\n",
    "    individual_height=1.0,\n",
    "    wspace=0.05,\n",
    "    hspace=0.0,\n",
    ")\n",
    "axs13c[0, 0].set_ylabel(\"Normalized\\nSpike Count\")\n",
    "axg13c.set_title(\"Low Gamma 20-50Hz\")\n",
    "gamma_band = twostim_gamma_power[\n",
    "    (twostim_gamma_power[\"gamma_band\"] == \"lowgamma\")\n",
    "    & (twostim_gamma_power[\"region\"] == \"OT\")\n",
    "]\n",
    "figure_withinregion_spike_phase(gamma_band, colors=\"b\", axs=axs13c[0, :])\n",
    "\n",
    "axs13d, axg13d = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=3,\n",
    "    group_top=0.6 + 1.0 + 0.8,\n",
    "    group_left=0.8 + 2.4 + 0.8,\n",
    "    individual_width=2.4 / 3,\n",
    "    individual_height=1.0,\n",
    "    wspace=0.05,\n",
    "    hspace=0.0,\n",
    ")\n",
    "axg13d.set_title(\"Low Gamma 20-50Hz\")\n",
    "gamma_band = twostim_gamma_power[\n",
    "    (twostim_gamma_power[\"gamma_band\"] == \"lowgamma\")\n",
    "    & (twostim_gamma_power[\"region\"] == \"Forebrain\")\n",
    "]\n",
    "figure_withinregion_spike_phase(gamma_band, colors=\"orange\", axs=axs13d[0, :])\n",
    "\n",
    "\n",
    "axs13e, axg13e = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=3,\n",
    "    group_top=0.6 + 1.0 + 0.6 + 1.0 + 1.0,\n",
    "    group_left=0.8,\n",
    "    individual_width=2.4 / 3,\n",
    "    individual_height=1.0,\n",
    "    wspace=0.05,\n",
    "    hspace=0.0,\n",
    ")\n",
    "axs13e[0, 0].set_ylabel(\"Normalized\\nSpike Count\")\n",
    "axg13e.set_title(\"High Gamma 50-75Hz\")\n",
    "gamma_band = twostim_gamma_power[\n",
    "    (twostim_gamma_power[\"gamma_band\"] == \"highgamma\")\n",
    "    & (twostim_gamma_power[\"region\"] == \"OT\")\n",
    "]\n",
    "figure_withinregion_spike_phase(gamma_band, colors=\"b\", axs=axs13e[0, :])\n",
    "\n",
    "axs13f, axg13f = figure_add_axes_group_inch(\n",
    "    fig,\n",
    "    nrows=1,\n",
    "    ncols=3,\n",
    "    group_top=0.6 + 1.0 + 0.6 + 1.0 + 1.0,\n",
    "    group_left=0.8 + 2.4 + 0.8,\n",
    "    individual_width=2.4 / 3,\n",
    "    individual_height=1.0,\n",
    "    wspace=0.05,\n",
    "    hspace=0.0,\n",
    ")\n",
    "gamma_band = twostim_gamma_power[\n",
    "    (twostim_gamma_power[\"gamma_band\"] == \"highgamma\")\n",
    "    & (twostim_gamma_power[\"region\"] == \"Forebrain\")\n",
    "]\n",
    "figure_withinregion_spike_phase(gamma_band, colors=\"orange\", axs=axs13f[0, :])\n",
    "axg13f.set_title(\"High Gamma 50-75Hz\")\n",
    "\n",
    "\n",
    "subplot_indicator(axg13a, \"A\", ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(axg13b, \"B\", ha=\"left\", va=\"bottom\", pad_inch=0.5)\n",
    "subplot_indicator(axg13c, \"C\", ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(axg13d, \"D\", ha=\"left\", va=\"bottom\", pad_inch=0.5)\n",
    "subplot_indicator(axg13e, \"E\", ha=\"left\", va=\"bottom\", pad_inch=0.7)\n",
    "subplot_indicator(axg13f, \"F\", ha=\"left\", va=\"bottom\", pad_inch=0.5)\n",
    "\n",
    "fig.savefig(OUTDIR / \"figure13.pdf\")\n",
    "fig.savefig(OUTDIR / \"figure13.png\")\n",
    "fig.savefig(OUTDIR / \"figure13.eps\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "374.391px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}