{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Abi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "abis = np.asarray([trial['abi'] for trial in abifile.trials])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spiketrains = np.empty((len(abifile.events), ), dtype='object')\n",
    "for kt, tevents in enumerate(abifile.events):\n",
    "    spiketimes = np.asarray([t for t, unit in tevents])\n",
    "    spiketrains[kt] = spiketimes / 10000\n",
    "order = np.argsort(abis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "abispikecounts = np.zeros((abifile.params['abi']['range'].size, abifile.params['reps']))\n",
    "for kabi, abi in enumerate(abifile.params['abi']['range']):\n",
    "    for krep, spiketimes in enumerate(spiketrains[abis == abi]):\n",
    "        abispikecounts[kabi, krep] = within(spiketimes, bffile.params['delay'], bffile.params['delay'] + \\\n",
    "                                        bffile.params['dur']*2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ILD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ilds = np.asarray([trial['iid'] for trial in ildfile.trials])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spiketrains = np.empty((len(ildfile.events), ), dtype='object')\n",
    "for kt, tevents in enumerate(ildfile.events):\n",
    "    spiketimes = np.asarray([t for t, unit in tevents])\n",
    "    spiketrains[kt] = spiketimes / 10000\n",
    "order = np.argsort(ilds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ildspikecounts = np.zeros((ildfile.params['iid']['range'].size, ildfile.params['reps']))\n",
    "for kabi, iid in enumerate(ildfile.params['iid']['range']):\n",
    "    for krep, spiketimes in enumerate(spiketrains[ilds == iid]):\n",
    "        ildspikecounts[kabi, krep] = within(spiketimes, bffile.params['delay'], bffile.params['delay'] + \\\n",
    "                                        bffile.params['dur']*2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ildx = ildfile.params['iid']['range']\n",
    "ildx = np.round(ildx, decimals=2)\n",
    "ildy = ildspikecounts.mean(axis=1)\n",
    "ildy = np.round(ildy, decimals=2)\n",
    "ildz = sps.stats.sem(ildspikecounts, axis=1)\n",
    "ildz[ildz==0] = 0.01\n",
    "so = [np.max(ildy), 1, 0, 1, np.min(ildy)]\n",
    "go = [np.min(ildy), np.max(ildy), ildx[np.argmax(ildy)], 10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    gopt, gcov = spo.curve_fit(Gaussian, ildx, ildy, go, sigma = ildz, absolute_sigma = True)\n",
    "    gcorr = np.corrcoef(ildy, Gaussian(ildx, *gopt))\n",
    "    Rsqg = gcorr[0,1]**2\n",
    "except RuntimeError:\n",
    "    print('No good Gaussian fit')\n",
    "    Rsqg = np.nan\n",
    "try:\n",
    "    sopt, scov = spo.curve_fit(Sigmoidal, ildx, ildy, so, sigma = ildz, absolute_sigma = True)\n",
    "    scorr = np.corrcoef(ildy, Sigmoidal(ildx, *sopt))\n",
    "    Rsqs = scorr[0,1]**2\n",
    "except RuntimeError:\n",
    "    print('no good Sigmoidal fit')\n",
    "    Rsqs = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "peak = np.where(ildy == ildy.max())\n",
    "ildhalfheight = np.where(ildy == find_nearest(ildy,(ildy.max() - ildy.min())/2 ) )\n",
    "if ildhalfheight[0].size == 0:\n",
    "    ildhalfwidth = np.empty([0])\n",
    "    ildhalfwidth = np.nan    \n",
    "elif ildhalfheight[0].size == 1:\n",
    "    ildhalfwidth = np.abs(ildx[peak] - ildx[ildhalfheight])\n",
    "elif ildhalfheight[0].size == 2: \n",
    "    ildwidth = ildhalfheight[0][1] - ildhalfheight[0][0]\n",
    "best_ild = ildx[peak]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### BF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bfs = np.asarray([trial['stim'] for trial in bffile.trials])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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": {
    "scrolled": true
   },
   "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']*2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bfx = bffile.params['bf']['range']\n",
    "bfy = bfspikecounts.mean(axis=1)\n",
    "bfse = sps.stats.sem(bfspikecounts, axis=1)"
   ]
  },
  {
   "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",
    "\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[0]\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])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**ITD**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "itds = np.asarray([trial['itd'] for trial in itdfile.trials])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "itdspiketrains = np.empty((len(itdfile.events), ), dtype='object')\n",
    "for kt, tevents in enumerate(itdfile.events):\n",
    "    spiketimes = np.asarray([t for t, unit in tevents])\n",
    "    itdspiketrains[kt] = spiketimes / 10000\n",
    "itdorder = np.argsort(itds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "itdspikecounts = np.zeros((itdfile.params['itd']['range'].size, itdfile.params['reps']))\n",
    "for kabi, stim in enumerate(itdfile.params['itd']['range']):\n",
    "    for krep, spiketimes in enumerate(itdspiketrains[itds == stim]):\n",
    "        itdspikecounts[kabi, krep] = within(spiketimes, itdfile.params['delay'], itdfile.params['delay'] + \\\n",
    "                                        itdfile.params['dur']*2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "itdx = itdfile.params['itd']['range']\n",
    "itdy = itdspikecounts.mean(axis=1)\n",
    "itdz = sps.stats.sem(itdspikecounts, axis=1)\n",
    "itdz[itdz==0] = 0.01\n",
    "H = lambda n: (bfx[n+1]-bfx[n]) * ((halfmax_y-bfy[n]) / (bfy[n+1]-bfy[n])) + bfx[n]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "peak_idx = np.argmax(itdy)\n",
    "halfmax_y = (np.max(itdy) - np.min(itdy))/2 + np.min(itdy)\n",
    "#halfmax_y = np.mean(itdy) + np.std(itdy)\n",
    "if np.where(itdy[:peak_idx] >= halfmax_y)[0].size > 0:\n",
    "    l = np.where(itdy[:peak_idx] >= halfmax_y)[0][0]\n",
    "else: #if no values besides peak is greater than halfmax_y\n",
    "    l = peak_idx - 1          \n",
    "if np.where(itdy[peak_idx+1:] >= halfmax_y)[0].size > 0: #peak__idx + 1 due to indexing issues\n",
    "    r = np.where(itdy[peak_idx+1:] >= halfmax_y)[0][-1] + peak_idx +1\n",
    "    if r == itdx.argmax():\n",
    "        r = r-1\n",
    "else:\n",
    "    r = peak_idx + 1 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if np.size(spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 150/10, prominence = 1)[0]) == 3:\n",
    "    [peak, locs] = spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 150/10, prominence = 1)\n",
    "    mainpk = np.arange(locs['left_bases'][1],locs['right_bases'][1])\n",
    "    po = np.array([np.min(itdy), locs['peak_heights'][1]- np.min(itdy[mainpk]), itdx[peak[1]], 50],dtype = object)\n",
    "elif np.size(spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 150/10, prominence = 1.4)[0]) == 1: \n",
    "    [peak, locs] = spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 150/10, prominence = 1.4)\n",
    "    mainpk = np.arange(locs['left_bases'][0],locs['right_bases'][0])\n",
    "    po = np.array([np.min(itdy), locs['peak_heights'][0]- np.min(itdy[mainpk]), itdx[peak[0]], 50],dtype = object)\n",
    "elif np.size(spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 150/10, prominence = 1.4)[0]) > 0:\n",
    "    [peak, locs] = spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 150/10, prominence = 1.4)\n",
    "    ppk = locs['peak_heights'].argmax()\n",
    "    mainpk = np.arange(locs['left_bases'][ppk],locs['right_bases'][ppk])\n",
    "    po = np.array([np.min(itdy), locs['peak_heights'][ppk]- np.min(itdy[mainpk]), itdx[peak[ppk]], 50],dtype = object)\n",
    "if np.size(spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 150/10, prominence = 1.4)[0]) > 0:\n",
    "    try:\n",
    "        bnds = [np.min(itdy) - itdz[np.argmin(itdy)], np.max(itdy[mainpk]) - np.min(itdy[mainpk]) - itdz[np.argmin(itdy[mainpk])], -550, 1], \\\n",
    "               [np.min(itdy) + itdz[np.argmin(itdy)], np.max(itdy[mainpk]) - np.min(itdy[mainpk]) + itdz[np.argmax(itdy[mainpk])], 550, 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",
    "else: \n",
    "    ITDRsq = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "H = lambda n: (itdx[n+1]-itdx[n]) * ((halfmax_y-itdy[n]) / (itdy[n+1]-itdy[n])) + itdx[n]\n",
    "if ITDRsq >.5:\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)])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if ~np.isnan(ITDRsq):\n",
    "    print('best ITD:')\n",
    "    display([best_itd.round(), ITD_lo.round(), ITD_hi.round(), ITDRsq.round(2)])\n",
    "    print('best freq:')\n",
    "    display([best_freq.round(),freq_lo.round(),freq_hi.round()])\n",
    "    print('best ILD:')\n",
    "    display(best_ild)\n",
    "    print('Gaussian or Sigmoidal')\n",
    "    display(np.round([Rsqg, Rsqs], decimals = 2))\n",
    "else:\n",
    "    print('best ITD:')\n",
    "    display([best_itd.round(), ITD_lo.round(), ITD_hi.round(), ITDRsq])\n",
    "    print('best freq:')\n",
    "    display([best_freq.round(),freq_lo.round(),freq_hi.round()])\n",
    "    print('best ILD:')\n",
    "    display(best_ild)\n",
    "    print('Gaussian or Sigmoidal')\n",
    "    display(np.round([Rsqg, Rsqs], decimals = 2))"
   ]
  }
 ],
 "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
}