{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyxdphys\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as sps\n",
    "import scipy.optimize as spo\n",
    "import scipy.signal as spg\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# matplotlib settings\n",
    "mpl.rcParams['mathtext.default'] = 'regular'\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['axes.labelsize'] =  18 # default: 'medium' == 10\n",
    "mpl.rcParams['xtick.labelsize'] = 18 # default: 'medium' == 10\n",
    "mpl.rcParams['ytick.labelsize'] = 18 # default: 'medium' == 10\n",
    "mpl.rcParams['legend.fontsize'] = 18\n",
    "\n",
    "# Other texts:\n",
    "in_panel_fontsize = 18\n",
    "bc_indicator_size = in_panel_fontsize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_nearest(array, value):\n",
    "    array = np.asarray(array)\n",
    "    idx = (np.abs(array - value)).argmin()\n",
    "    return [idx]\n",
    "\n",
    "def Gaussian(t, ga, gb, gc, gd):\n",
    "    return ga + np.abs(gb) * np.exp( - ( (t-gc)**2 / gd**2) )\n",
    "\n",
    "def Sigmoidal(t, sa, sb, sc, sd, se):\n",
    "    return sa / (1 + np.exp(-sb*(t - sc)/sd)) + se\n",
    "\n",
    "def within(tval, tmin, tmax):\n",
    "    return np.sum((tval > tmin) & (tval < tmax))\n",
    "\n",
    "def between(val, array):\n",
    "    return array[0] <= val <= array[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Select neuron. d = date, m = owl.neuron\n",
    "d = '2021-0929'\n",
    "m =  '44.02'\n",
    "n = '0' + m\n",
    "\n",
    "#alt, in cases where abi was recorded using using monaural and binuaral (L,R,B):\n",
    "#abifile = pyxdphys.XDdata(fr'F:\\xdphys\\{d}\\{m}\\{n}.0.1.abi')\n",
    "abifile  = pyxdphys.XDdata(fr'F:\\xdphys\\{d}\\{m}\\{n}.0.abi')\n",
    "ildfile  = pyxdphys.XDdata(fr'F:\\xdphys\\{d}\\{m}\\{n}.1.iid')\n",
    "bcfile   = pyxdphys.XDdata(fr'F:\\xdphys\\{d}\\{m}\\{n}.2.gen')\n",
    "itdfile  = pyxdphys.XDdata(fr'F:\\xdphys\\{d}\\{m}\\{n}.3.itd')\n",
    "bffile   = pyxdphys.XDdata(fr'F:\\xdphys\\{d}\\{m}\\{n}.4.bf')\n",
    "#sptfile  = pyxdphys.XDdata(r'D:\\xdphys\\2021-0811\\046.00\\046.00.5.iid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best ITD:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[20.0, -21.0, 61.0, 0.99]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best freq:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[4900.0, 4000.0, 5800.0]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best ILD:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([4.])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gaussian or Sigmoidal\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([0.99, 0.08])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "##Running this will do all analysis automatically. The rest of this script is in case there are issues,\n",
    "## or to see rasterplots/plots\n",
    "%run ./Get_Spike_Counts.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#run next cell if best ITD isn't found. This usually fixes the issue. xlsx notes when this is used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "siteITD = 0\n",
    "[peak, locs] = spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 100/10, prominence = 1.)\n",
    "\n",
    "tt = find_nearest(itdx[peak], siteITD)\n",
    "\n",
    "mainpk = np.arange(locs['left_bases'][tt],locs['right_bases'][tt])\n",
    "po = np.array([np.min(itdy), locs['peak_heights'][tt]- np.min(itdy[mainpk]), itdx[peak[tt]], 50],dtype=object)\n",
    "\n",
    "try:\n",
    "    bnds = [[np.min(itdy) - itdz[np.argmin(itdy)], np.max(itdy[mainpk]) - np.min(itdy[mainpk]) - itdz[np.argmin(itdy[mainpk])], itdx[mainpk].min() - 100, 1], \\\n",
    "           [np.min(itdy) + itdz[np.argmin(itdy)], np.max(itdy[mainpk]) - np.min(itdy[mainpk]) + itdz[np.argmax(itdy[mainpk])], itdx[mainpk].max() + 100, 300]]\n",
    "    ITDopt, ITDcov = spo.curve_fit(Gaussian, itdx[mainpk], itdy[mainpk], po, sigma = itdz[mainpk],\n",
    "                                   bounds = bnds, absolute_sigma = True)\n",
    "    ITDcorr = np.corrcoef(itdy[mainpk], Gaussian(itdx[mainpk], *ITDopt))\n",
    "    ITDRsq = ITDcorr[0,1]**2\n",
    "except RuntimeError:\n",
    "    print('no good Gaussian fit')\n",
    "    ITDRsq = np.nan\n",
    "    \n",
    "H = lambda n: (itdx[n+1]-itdx[n]) * ((halfmax_y-itdy[n]) / (itdy[n+1]-itdy[n])) + itdx[n]\n",
    "if ITDRsq >.25:\n",
    "    gITD = Gaussian(itdx, *ITDopt)\n",
    "    xint = np.arange(itdx[0],itdx[-1])\n",
    "    gint = np.interp(np.arange(itdx[0],itdx[-1]), itdx, gITD)\n",
    "    gpeak_idx = np.argmax(gint)\n",
    "    ghalfpeak = (np.max(gint) - np.min(gint))/2 + np.min(gint)\n",
    "    best_itd = xint[gpeak_idx]\n",
    "    ghalfwidth = np.abs(xint[np.abs(gint-ghalfpeak).argmin()] - xint[gpeak_idx])\n",
    "    ITD_lo = xint[gpeak_idx] - ghalfwidth\n",
    "    ITD_hi = xint[gpeak_idx] + ghalfwidth\n",
    "    display([np.mean((H(l), H(r))).round(2),xint[gpeak_idx]], \\\n",
    "            [H(l).round(2), xint[gpeak_idx] - ghalfwidth], [H(r).round(2),xint[gpeak_idx] + ghalfwidth])\n",
    "else:\n",
    "    best_itd = np.mean((H(l), H(r)))\n",
    "    ITD_lo = H(l)\n",
    "    ITD_hi = H(r)\n",
    "    display([np.mean((H(l), H(r)))], [H(l)], [H(r)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Best ITD Graph\n",
    "fig, ax = plt.subplots(figsize = [5,4])\n",
    "ax.errorbar(itdfile.params['itd']['range'], itdspikecounts.mean(axis=1), sps.stats.sem(itdspikecounts, axis=1))\n",
    "ax.plot(itdx, Gaussian(itdx, *ITDopt))\n",
    "plt.ylabel('Spike Count')\n",
    "plt.xlabel('ITD (us)')\n",
    "\n",
    "#fig.savefig('ExampleITD.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ITD raster plot\n",
    "fig, ax = plt.subplots()\n",
    "ax.eventplot(itdspiketrains[itdorder], color='k')\n",
    "ax.set_yticks(np.arange(itdfile.params['itd']['range'].size)[::4] \\\n",
    "              * itdfile.params['reps'] + itdfile.params['reps']/2)\n",
    "ax.set_yticklabels(itdfile.params['itd']['range'][::4])\n",
    "\n",
    "ax.axvspan(itdfile.params['delay'], itdfile.params['delay']+itdfile.params['dur']+10, color='yellow')\n",
    "ax.axvspan(itdfile.params['delay'], itdfile.params['delay']+ 15, color='pink')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#run the next four cells to change bffile analysis window to just the sound duration\n",
    "bfs = np.asarray([trial['stim'] for trial in bffile.trials])\n",
    "bfspiketrains = np.empty((len(bffile.events), ), dtype='object')\n",
    "for kt, tevents in enumerate(bffile.events):\n",
    "    spiketimes = np.asarray([t for t, unit in tevents])\n",
    "    bfspiketrains[kt] = spiketimes / 10000\n",
    "bforder = np.argsort(bfs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bfspikecounts = np.zeros((bffile.params['bf']['range'].size, bffile.params['reps']))\n",
    "\n",
    "for kabi, stim in enumerate(bffile.params['bf']['range']):\n",
    "    for krep, spiketimes in enumerate(bfspiketrains[bfs == stim]):\n",
    "        bfspikecounts[kabi, krep] = within(spiketimes, bffile.params['delay'], bffile.params['delay'] + \\\n",
    "                                        bffile.params['dur']*1)\n",
    "#print(bfspikecounts.shape)\n",
    "#fig, ax = plt.subplots(figsize = [5,4])\n",
    "#\n",
    "#ax.errorbar(bffile.params['bf']['range'], bfspikecounts.mean(axis=1), sps.stats.sem(bfspikecounts, axis=1))\n",
    "#plt.xlabel('Frequency (Hz)')\n",
    "#plt.ylabel('Spike Count')\n",
    "##fig.savefig('ExampleFreq.png')\n",
    "bfy = bfspikecounts.mean(axis=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "peak_idx = np.argmax(bfy)\n",
    "halfmax_y = (np.max(bfy) - np.min(bfy))/2 + np.min(bfy)\n",
    "lup = np.where(bfy[:peak_idx] > halfmax_y)[0]\n",
    "ldown = np.where(bfy[:peak_idx] < halfmax_y)[0]\n",
    "rup = np.where(bfy[peak_idx+1:] > halfmax_y)[0]\n",
    "rdown = np.where(bfy[peak_idx+1:] < halfmax_y)[0]\n",
    "pm = np.array([-1,1])\n",
    "\n",
    "\n",
    "if lup.size > 0 and ldown.size > 0:\n",
    "    if between(halfmax_y, bfy[lup[0]] +  pm*bfse[lup[0]]):\n",
    "        l = lup[1]\n",
    "    elif between(halfmax_y, bfy[ldown[0]] +  pm*bfse[ldown[0]]):\n",
    "        l = ldown[-1] \n",
    "    elif lup[0] == ldown[-1] + 1:\n",
    "        l = lup[0] \n",
    "    else:\n",
    "        l = lup[0]\n",
    "elif lup.size > 0:\n",
    "    l = lup[0]\n",
    "elif ldown.size > 0: \n",
    "    l = ldown[-1] \n",
    "else:\n",
    "    l = peak_idx\n",
    "\n",
    "if rup.size > 0 and rdown.size > 0:\n",
    "    if between(halfmax_y, bfy[rup[-1] + peak_idx + 1] +  pm*bfse[rup[-1] + peak_idx + 1]):\n",
    "        r = rup[-1] + peak_idx + 1\n",
    "    elif between(halfmax_y, bfy[rdown[0] + peak_idx + 1] +  pm*bfse[rdown[0] + peak_idx + 1]):\n",
    "        r = rdown[0] + peak_idx +1\n",
    "    elif rup[-1] == rdown[0] - 1:\n",
    "        r = rup[-1] + peak_idx + 1\n",
    "    else:\n",
    "        r = rup[-1] + peak_idx + 1\n",
    "elif rup.size > 0:\n",
    "    r = rup[-1] + peak_idx + 1\n",
    "elif rdown.size > 0: \n",
    "    r = rdown[0] + peak_idx +1\n",
    "else:\n",
    "    r = peak_idx   \n",
    "\n",
    "    \n",
    "freq_lo, lofr = bfx[l], bfy[l]\n",
    "freq_hi, hifr = bfx[r], bfy[r]\n",
    "best_freq = np.mean([freq_lo, freq_hi])\n",
    "display([best_freq, freq_lo, freq_hi])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Best Freq\n",
    "\n",
    "fig, ax = plt.subplots(figsize = [5,4])\n",
    "\n",
    "ax.errorbar(bffile.params['bf']['range'], bfspikecounts.mean(axis=1), sps.stats.sem(bfspikecounts, axis=1))\n",
    "plt.xlabel('Frequency (Hz)')\n",
    "plt.ylabel('Spike Count')\n",
    "\n",
    "ax.scatter([freq_lo, freq_hi, best_freq], [lofr, hifr, np.max(bfy)], color = 'black')\n",
    "plt.hlines((np.max(bfy) - np.min(bfy))/2 + np.min(bfy), 400, 10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Frequency raster plot\n",
    "fig, ax = plt.subplots()\n",
    "ax.eventplot(bfspiketrains[bforder], color='k')\n",
    "ax.set_yticks(np.arange(bffile.params['bf']['range'].size)[::4] \\\n",
    "              * bffile.params['reps'] + bffile.params['reps']/2)\n",
    "ax.set_yticklabels(bffile.params['bf']['range'][::4])\n",
    "\n",
    "ax.axvspan(bffile.params['delay'], bffile.params['delay']+bffile.params['dur']+10, color='yellow')\n",
    "ax.axvspan(bffile.params['delay'], bffile.params['delay']+ 10, color='pink')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ILD graph\n",
    "fig, ax = plt.subplots(figsize= [5,4])\n",
    "ax.errorbar(ildx, ildy, sps.stats.sem(ildspikecounts, axis=1), linewidth=2, color='black',marker = 'o')\n",
    "ax.plot(ildx, Gaussian(ildx, *gopt), 'green')\n",
    "ax.plot(ildx, Sigmoidal(ildx, *sopt), 'brown')\n",
    "\n",
    "plt.xlabel('ILD (dB)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.12"
  },
  "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": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}