{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "Created on Mon Aug  6 13:16:08 2018\n",
    "\n",
    "@author: kperks\n",
    "\"\"\"\n",
    "# INITIALIZATION\n",
    "from neo import Spike2IO\n",
    "import os\n",
    "import numpy as np\n",
    "import scipy.io as sio\n",
    "import scipy \n",
    "import h5py\n",
    "from scipy import stats\n",
    "from scipy import signal\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.axes as ax\n",
    "from pathlib import Path\n",
    "import random\n",
    "import sys\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from scipy.stats import spearmanr, ks_2samp "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_range_bounds(win_times,trigger):\n",
    "    start = np.min(np.where(trigger>win_times[0])[0])\n",
    "    stop= np.max(np.where(trigger<win_times[1])[0])\n",
    "    trials = list(range(start,stop))\n",
    "    return trials\n",
    "    \n",
    "def TrialEventTimes(trigger,analevent,trial,win):\n",
    "    t1 = trigger[trial]\n",
    "    t2 = t1+win\n",
    "    EventTimes = analevent[np.where((analevent>t1)&(analevent<t2))]-t1\n",
    "    return(EventTimes)\n",
    "    \n",
    "def get_response(trigger_times,data_func,trial_duration,tau):\n",
    "    response = [data_func(np.linspace(t,t+trial_duration,int(trial_duration/dt))) for t in trigger_times]\n",
    "    return np.asarray(response)\n",
    "    \"\"\"\n",
    "    Creates a spike density function from a spike train to do continuous analysis with (such as trial correlations)\n",
    "    \"\"\" \n",
    "    \n",
    "def filtered_response(spk_times, tau):\n",
    "    \"\"\"\n",
    "    Creates a function with a gaussian centered at every spike time (the mean) with standard deviation tau\n",
    "    \"\"\"\n",
    "    spk_times = spk_times.reshape((-1, 1))\n",
    "    norm_factor = tau * np.sqrt(2. * np.pi)\n",
    "    return lambda t: np.sum(np.exp(-(spk_times - t.reshape((1, -1))) ** 2 / (2 * tau * tau)), 0) / norm_factor\n",
    "    \n",
    "def raster(trigger, data, **kwargs):\n",
    "    ax = plt.gca();\n",
    "    \n",
    "    for ith, t in enumerate(trigger):\n",
    "        inds = np.where((data>t)&(data<t+md.trial_duration[0]))[0]\n",
    "        trialdata = data[inds]-t\n",
    "        plt.scatter(trialdata, np.repeat(ith,len(trialdata)), **kwargs);\n",
    "        \n",
    "    plt.ylim(.5, len(trigger) + .5);\n",
    "\n",
    "###########\n",
    "#functions to estimate optimal bandwidth for spike density:\n",
    "def CostFunction(y_hist, N, w, dt):\n",
    "    \"\"\"\n",
    "    References\n",
    "    ----------\n",
    "    .. [1] H. Shimazaki and S. Shinomoto, \"Kernel Bandwidth Optimization in \n",
    "           Spike Rate Estimation,\" in Journal of Computational Neuroscience \n",
    "           29(1-2): 171–182, 2010 http://dx.doi.org/10.1007/s10827-009-0180-4\n",
    "    \"\"\"\n",
    "    # build normal smoothing kernel\n",
    "    yh = fftkernel(y_hist, w / dt)\n",
    "\n",
    "    # formula for density\n",
    "    C = np.sum(yh**2) * dt - 2 * np.sum(yh * y_hist) * dt + 2 \\\n",
    "        / (2 * np.pi)**0.5 / w / N\n",
    "    C = C * N**2\n",
    "\n",
    "    return C, yh\n",
    "\n",
    "\n",
    "def fftkernel(x, w):\n",
    "    \"\"\"\n",
    "    References\n",
    "    ----------\n",
    "    .. [1] H. Shimazaki and S. Shinomoto, \"Kernel Bandwidth Optimization in \n",
    "           Spike Rate Estimation,\" in Journal of Computational Neuroscience \n",
    "           29(1-2): 171–182, 2010 http://dx.doi.org/10.1007/s10827-009-0180-4\n",
    "    \"\"\"\n",
    "    L = x.size\n",
    "    Lmax = L + 3 * w\n",
    "    n = 2 ** np.ceil(np.log2(Lmax))\n",
    "\n",
    "    X = np.fft.fft(x, n.astype(np.int))\n",
    "\n",
    "    f = np.linspace(0, n-1, n) / n\n",
    "    f = np.concatenate((-f[0: np.int(n / 2 + 1)],\n",
    "                        f[1: np.int(n / 2 - 1 + 1)][::-1]))\n",
    "\n",
    "    K = np.exp(-0.5 * (w * 2 * np.pi * f) ** 2)\n",
    "\n",
    "    y = np.real(np.fft.ifft(X * K, n))\n",
    "\n",
    "    y = y[0:L]\n",
    "\n",
    "    return y\n",
    "\n",
    "def sskernel(x, tin, W):\n",
    "    \"\"\"\n",
    "    Generates a kernel density estimate with globally-optimized bandwidth.\n",
    "    The optimal bandwidth is obtained as a minimizer of the formula, sum_{i,j}\n",
    "    \\int k(x - x_i) k(x - x_j) dx  -  2 sum_{i~=j} k(x_i - x_j), where k(x) is\n",
    "    the kernel function.\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : array_like\n",
    "        The one-dimensional samples drawn from the underlying density\n",
    "    tin : array_like\n",
    "        The values where the density estimate is to be evaluated in generating\n",
    "        the output 'y'.\n",
    "    W : array_like\n",
    "        The kernel bandwidths to use in optimization. Should not be chosen\n",
    "        smaller than the sampling resolution of 'x'.\n",
    "    nbs : int, optional\n",
    "        The number of bootstrap samples to use in estimating the [0.05, 0.95]\n",
    "        confidence interval of the output 'y'\n",
    "    Returns\n",
    "    -------\n",
    "    y : array_like\n",
    "        The estimated density, evaluated at points t / tin.\n",
    "    t : array_like\n",
    "        The points where the density estimate 'y' is evaluated.\n",
    "    optw : double\n",
    "        The optimal global kernel bandwidth.\n",
    "    C : array_like\n",
    "        The cost functions associated with the bandwidths 'W'.\n",
    "\n",
    "    References\n",
    "    ----------\n",
    "    .. [1] H. Shimazaki and S. Shinomoto, \"Kernel Bandwidth Optimization in \n",
    "           Spike Rate Estimation,\" in Journal of Computational Neuroscience \n",
    "           29(1-2): 171–182, 2010 http://dx.doi.org/10.1007/s10827-009-0180-4\n",
    "    \"\"\"\n",
    "\n",
    "    # set argument 't' if not provided\n",
    "    #always provide xtime input so that output is aligned to trial duration\n",
    "    T = np.max(x) - np.min(x)\n",
    "    x_ab = x[(x >= np.min(tin)) & (x <= np.max(tin))]\n",
    "    dx = np.sort(np.diff(np.sort(x)))\n",
    "    dt_samp = dx[np.nonzero(dx)][0]\n",
    "    if dt_samp > np.min(np.diff(tin)):\n",
    "        t = np.linspace(np.min(tin), np.max(tin), np.min([int(np.ceil(T / dt_samp)), int(1e3)]))\n",
    "    else:\n",
    "        t = tin\n",
    "\n",
    "    # calculate delta t\n",
    "    dt = min(np.diff(t))\n",
    "\n",
    "    # create the finest histogram\n",
    "    thist = np.concatenate((t, (t[-1]+dt)[np.newaxis]))\n",
    "    y_hist = np.histogram(x_ab, thist-dt/2)[0]\n",
    "    N = sum(y_hist).astype(np.float)\n",
    "    y_hist = y_hist / N / dt\n",
    "\n",
    "    # always provide W\n",
    "    C = np.zeros((1, len(W)))[0]\n",
    "    C_min = np.Inf\n",
    "    ymat=[]\n",
    "    for k, w in enumerate(W):\n",
    "        C[k], yh = CostFunction(y_hist, N, w, dt)\n",
    "        ymat.append(yh)\n",
    "        if((C[k] < C_min).any()):\n",
    "            C_min = C[k]\n",
    "            optw = w\n",
    "            y = yh\n",
    "    ymat = np.asarray(ymat)       \n",
    "    \n",
    "    return y, t, optw, C, C_min\n",
    "\n",
    "def get_spikes_kde(trigger_times,data_times,trigger,trial_duration):\n",
    "    aligned_times = []\n",
    "    for t in trigger_times:\n",
    "        thesetimes = data_times[\n",
    "          np.where((data_times>(trigger[t]))&(data_times<(trigger[t]+trial_duration)))[0]\n",
    "          ]-(trigger[t])\n",
    "        aligned_times = np.concatenate((aligned_times,thesetimes)) \n",
    "\n",
    "    return np.asarray(aligned_times)\n",
    "\n",
    "def get_optw(these_trials,spikes,xtime,trigger,trial_duration):\n",
    "    W=np.arange(0.006,0.2,0.002)\n",
    "    aligned_times = get_spikes_kde(these_trials,spikes,trigger,trial_duration)\n",
    "    y, t, optw, C, C_min = sskernel(aligned_times, xtime, W)\n",
    "    \n",
    "    return y, t, optw, C, C_min"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "meta = pd.read_csv('MetaData.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are the summary statistics for the trial duration in paired experiments? (for Methods section)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "count    33.000000\n",
       "mean      2.757576\n",
       "std       1.565603\n",
       "min       1.000000\n",
       "25%       1.500000\n",
       "50%       2.000000\n",
       "75%       4.000000\n",
       "max       6.000000\n",
       "Name: trial_duration, dtype: float64"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "meta[((meta['trigger_type']=='paired')\n",
    "      &(meta['phase']=='phase1')\n",
    "      &((meta['condition']=='aen_vent') \n",
    "        | (meta['condition']=='aen_proprio') \n",
    "       |(meta['condition']=='aff_proprio')))].trial_duration.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 20050727_3224_SPK.mat\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.6/site-packages/neo-0.8.0.dev0-py3.6.egg/neo/rawio/spike2rawio.py:609: RuntimeWarning: overflow encountered in short_scalars\n",
      "  info['time_per_adc']) * 1e-6\n",
      "/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:61: DeprecationWarning: object of type <class 'numpy.float64'> cannot be safely interpreted as an integer.\n",
      "/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:96: DeprecationWarning: object of type <class 'numpy.float64'> cannot be safely interpreted as an integer.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 20050728_3404_latency1_SPK.mat\n",
      "2 20050728_3404_latency2_SPK.mat\n",
      "3 20050808_3629_latency1_SPK.mat\n",
      "4 20050808_3629_latency2_SPK.mat\n",
      "5 20050811_3221_SPK.mat\n",
      "6 20050818_2874_SPK.mat\n",
      "7 20050819_2601_SPK.mat\n",
      "8 20050820_2960_trial1_SPK.mat\n",
      "9 20050824_2627_control_SPK.mat\n",
      "10 20050824_2627_trial1_SPK.mat\n",
      "11 20060626_2815_SPK.mat\n",
      "12 20060627_2673_SPK.mat\n",
      "13 20060708_3004_SPK.mat\n",
      "14 20060710_2767_SPK.mat\n",
      "15 20060710_3196_SPK.mat\n",
      "16 20060804_2107_SPK.mat\n",
      "17 20060812_2921_SPK.mat\n",
      "18 20060812_3066_SPK.mat\n",
      "19 20060812_3205_SPK.mat\n",
      "20 20060816_3893_SPK.mat\n",
      "21 20070718_598_coupled_SPK.mat\n",
      "22 20070719_91_coupled_SPK.mat\n",
      "23 20070720_3352_SPK.mat\n",
      "24 20070723_287_SPK.mat\n",
      "25 20070723_311_SPK.mat\n",
      "26 20070723_339_SPK.mat\n",
      "27 20070723_415_SPK.mat\n",
      "28 20070727_3798_SPK.mat\n",
      "29 20070727_3798_control_SPK.mat\n",
      "30 20070816_3375_SPK.mat\n",
      "31 20070816_3509_SPK.mat\n",
      "32 20071025_3335_SPK.mat\n",
      "33 20071025_3335_control_SPK.mat\n",
      "using pre pairing expt to get baseline (recovery end)\n",
      "34 20071026_3108_SPK.mat\n",
      "35 20071026_3108_control_SPK.mat\n",
      "using pre pairing expt to get baseline\n",
      "36 20170810_expt5_vent_SPK.mat\n",
      "37 20180124_expt1_swim_SPK.mat\n",
      "38 20170815_expt8_vent_SPK.mat\n",
      "39 20170815_expt9_vent_SPK.mat\n",
      "40 20180124_expt2_swim_SPK.mat\n",
      "41 20180407_freerun2_SPK.mat\n",
      "42 20180407_freerun3_SPK.mat\n",
      "43 20180407_freerun4_SPK.mat\n",
      "44 20180407_freerun5_SPK.mat\n",
      "45 20170815_expt10_swim_SPK.mat\n",
      "46 20170810_expt4_vent_SPK.mat\n",
      "47 20050808_3629_latency2b_SPK.mat\n",
      "48 20050808_3629_latency1b_SPK.mat\n",
      "49 20070801_3838_SPK.mat\n",
      "50 20070801_3838_control_SPK.mat\n",
      "51 20070808_3727_SPK.mat\n"
     ]
    }
   ],
   "source": [
    "######MAIN##########\n",
    "\n",
    "meta = pd.read_csv('MetaData.csv')\n",
    "\n",
    "for iexpt,thisexpt_name in enumerate(meta.exptname):\n",
    "    dt = 1/1000 \n",
    "    ntrials_per_period = 75 #results did not depend on choice of 50, 75, or 100\n",
    "\n",
    "    print(iexpt,thisexpt_name)\n",
    "\n",
    "    thisexpt = meta.loc[meta['exptname'] == thisexpt_name]\n",
    "\n",
    "    trial_duration = thisexpt.trial_duration.values[0]\n",
    "\n",
    "    basename = thisexpt_name[0:-8]\n",
    "    data_folder = Path.cwd() / basename\n",
    "    file_to_open = data_folder / Path(basename + '.smr')\n",
    "    bl = Spike2IO(file_to_open,try_signal_grouping=False).read_block()\n",
    "\n",
    "    trigger = []\n",
    "    for sublist in np.asarray([seg.events[[s.name for \n",
    "       s in seg.events].index(thisexpt.chan_trigger.values[0])].magnitude \n",
    "       for seg in bl.segments]):\n",
    "        for item in sublist:\n",
    "            trigger.append(item)\n",
    "    trigger = np.asarray(trigger)\n",
    "\n",
    "    if thisexpt.spike_times_from.values[0] == 'Spyking-Circus':\n",
    "        spikes = np.load(data_folder / 'spikes.npy')\n",
    "\n",
    "    if thisexpt.spike_times_from.values[0] == 'Spike2':\n",
    "        spikes = []\n",
    "        for sublist in np.asarray([seg.events[[s.name for \n",
    "           s in seg.events].index(thisexpt.chan_spikes.values[0])].magnitude \n",
    "           for seg in bl.segments]):\n",
    "            for item in sublist:\n",
    "                spikes.append(item)\n",
    "        spikes = np.asarray(spikes)\n",
    "\n",
    "    postpairing_range = get_range_bounds([thisexpt.post_start.values[0],\n",
    "                                          thisexpt.post_stop.values[0]],trigger)\n",
    "    stim_range = get_range_bounds([thisexpt.stim_start.values[0],\n",
    "                                   thisexpt.stim_stop.values[0]],trigger)\n",
    "    base_range = get_range_bounds([thisexpt.base_start.values[0],\n",
    "                                   thisexpt.base_stop.values[0]],trigger)\n",
    "\n",
    "    if thisexpt_name == '20060812_3066_SPK.mat':\n",
    "        #removing a few trials in which spike detection was clearly mis-triggered due to movement artifacts\n",
    "        removeinds = np.arange(360,np.max(postpairing_range)-1,1)\n",
    "        [postpairing_range.remove(x) for x in removeinds]\n",
    "        removeinds = [263, 264, 265, 326, 327] \n",
    "        [postpairing_range.remove(x) for x in removeinds]\n",
    "\n",
    "    xtime = np.linspace(0,trial_duration,trial_duration/dt)\n",
    "\n",
    "    iei = [(TrialEventTimes(trigger,trigger,trial,10))[0] if len(TrialEventTimes(trigger,trigger,trial,10))>0 else 0 for trial in list(range(0,len(trigger)-1))]\n",
    "\n",
    "    these_trials = base_range\n",
    "    y, t, optw, C, C_min= get_optw(these_trials,spikes,xtime,trigger,trial_duration)\n",
    "\n",
    "    data_func = filtered_response(spikes, optw)\n",
    "\n",
    "    baseline_r = get_response(trigger[these_trials],data_func,trial_duration,optw)\n",
    "    baseline_mean = np.mean(baseline_r,0)\n",
    "\n",
    "    n_trials = len(base_range)\n",
    "    base_spkrt = [np.shape(np.where((spikes>(trigger[t]))&(spikes<(trigger[t]+trial_duration)))[0])[0]\n",
    "                  /trial_duration for t in base_range]\n",
    "\n",
    "    spikemat_base = [spikes[np.where((spikes>trigger[t])&\n",
    "                    (spikes<trigger[t]+trial_duration))[0]\n",
    "                    ]-(trigger[t]) for t in base_range]\n",
    "\n",
    "    if thisexpt_name == '20071025_3335_control_SPK.mat': \n",
    "    #not enough baseline recorded so get estimate of baseline from test cell clock trigger baseline\n",
    "        print('using pre pairing expt to get baseline (recovery end)')\n",
    "        suppname = '20071025_3335_SPK.mat'\n",
    "        suppexpt = meta.loc[meta['exptname'] == suppname]\n",
    "        basename = suppname[0:-8]\n",
    "        data_folder = Path.cwd() / basename\n",
    "        supptrigger = np.arange(suppexpt.post_stop.values[0]-250,suppexpt.post_stop.values[0],trial_duration)\n",
    "        suppbase_range = np.arange(0,len(supptrigger),1)\n",
    "        base_range = suppbase_range\n",
    "\n",
    "        suppspikes = np.load(data_folder / 'spikes.npy')\n",
    "\n",
    "        these_trials = suppbase_range\n",
    "        y, t, optw, C, C_min= get_optw(these_trials,suppspikes,xtime,supptrigger,trial_duration)\n",
    "\n",
    "        data_func = filtered_response(suppspikes, optw)\n",
    "\n",
    "        baseline_mean = np.mean(get_response(\n",
    "            supptrigger[these_trials],data_func,trial_duration,optw),0)\n",
    "\n",
    "        base_spkrt = [np.shape(np.where((suppspikes>(supptrigger[t]))&\n",
    "                                        (suppspikes<(supptrigger[t]+trial_duration)))[0])[0]\n",
    "                      /trial_duration for t in suppbase_range]\n",
    "        spikemat_base = [suppspikes[np.where((suppspikes>supptrigger[t])&\n",
    "                (suppspikes<supptrigger[t]+trial_duration))[0]\n",
    "                ]-(supptrigger[t]) for t in suppbase_range]\n",
    "\n",
    "\n",
    "    if thisexpt_name == '20071026_3108_control_SPK.mat': \n",
    "    #not enough baseline recorded so get estimate of baseline from test cell clock trigger baseline\n",
    "        print('using pre pairing expt to get baseline')\n",
    "        suppname = '20071026_3108_SPK.mat'\n",
    "        suppexpt = meta.loc[meta['exptname'] == suppname]\n",
    "        basename = suppname[0:-8]\n",
    "        data_folder = Path.cwd() / basename\n",
    "        file_to_open = data_folder / Path(basename + '.smr')\n",
    "        bl = Spike2IO(file_to_open,try_signal_grouping=False).read_block()\n",
    "        supptrigger = np.arange(0,suppexpt.base_stop.values[0],trial_duration)\n",
    "        suppbase_range = np.arange(0,len(supptrigger),1)\n",
    "        base_range = suppbase_range\n",
    "\n",
    "        suppspikes = []\n",
    "        for sublist in np.asarray([seg.events[[s.name for \n",
    "           s in seg.events].index(suppexpt.chan_spikes.values[0])].magnitude \n",
    "           for seg in bl.segments]):\n",
    "            for item in sublist:\n",
    "                suppspikes.append(item)\n",
    "        suppspikes = np.asarray(suppspikes)\n",
    "\n",
    "        these_trials = suppbase_range\n",
    "        y, t, optw, C, C_min= get_optw(these_trials,suppspikes,xtime,supptrigger,trial_duration)\n",
    "\n",
    "        data_func = filtered_response(suppspikes, optw)\n",
    "\n",
    "        baseline_mean = np.mean(get_response(\n",
    "            supptrigger[these_trials],data_func,trial_duration,optw),0)\n",
    "\n",
    "        base_spkrt = [np.shape(np.where((suppspikes>(supptrigger[t]))&\n",
    "                                        (suppspikes<(supptrigger[t]+trial_duration)))[0])[0]\n",
    "                      /trial_duration for t in suppbase_range]\n",
    "        spikemat_base = [suppspikes[np.where((suppspikes>supptrigger[t])&\n",
    "                (suppspikes<supptrigger[t]+trial_duration))[0]\n",
    "                ]-(supptrigger[t]) for t in suppbase_range]\n",
    "\n",
    "    data_func = filtered_response(spikes, optw)\n",
    "\n",
    "    npair = len(stim_range)\n",
    "    ntrials = np.min([npair,ntrials_per_period*2])\n",
    "    stiminitial_trials = stim_range[0:int(ntrials/2)]\n",
    "    these_trials = stiminitial_trials\n",
    "    stimi_spkrt = [np.shape(np.where((spikes>(trigger[t]))&(spikes<(trigger[t]+trial_duration)))[0])[0]\n",
    "                   /trial_duration for t in these_trials]\n",
    "    y, t, optw_stim, C, C_min= get_optw(these_trials,spikes,xtime,trigger,trial_duration)\n",
    "    stiminitial_r = get_response(trigger[these_trials],data_func,trial_duration,optw)\n",
    "    stiminitial_response = np.mean(stiminitial_r,0)\n",
    "\n",
    "    stimfinal_trials = stim_range[-int(ntrials/2):]\n",
    "    these_trials = stimfinal_trials\n",
    "    stimf_spkrt = [np.shape(np.where((spikes>(trigger[t]))&(spikes<(trigger[t]+trial_duration)))[0])[0]\n",
    "                   /trial_duration for t in these_trials]\n",
    "    stimfinal_r = get_response(trigger[these_trials],data_func,trial_duration,optw)\n",
    "    stimfinal_response = np.mean(stimfinal_r,0)\n",
    "\n",
    "    npost = len(postpairing_range)\n",
    "    ntrials = np.min([npost,ntrials_per_period])\n",
    "    postpairing_trials = postpairing_range[0:ntrials]\n",
    "    these_trials = postpairing_trials\n",
    "    post_spkrt = [np.shape(np.where((spikes>(trigger[t]))&(spikes<(trigger[t]+trial_duration)))[0])[0]\n",
    "                  /trial_duration for t in these_trials]\n",
    "    postpairing_r = get_response(trigger[these_trials],data_func,trial_duration,optw)\n",
    "    postpairing_response = np.mean(postpairing_r,0)\n",
    "\n",
    "    recovery_spkrt = np.nan\n",
    "    recovery_response = np.full_like(np.empty(len(baseline_mean)),np.nan)\n",
    "    recovery_trials = [np.nan]\n",
    "    if npost>=(300+ntrials_per_period):\n",
    "        recovery_trials = postpairing_range[300:300+ntrials_per_period]\n",
    "        recovery_spkrt = [np.shape(np.where((spikes>(trigger[t]))&(spikes<(trigger[t]+trial_duration)))[0])[0]\n",
    "                          /trial_duration for t in recovery_trials]\n",
    "        these_trials = recovery_trials\n",
    "        recovery_r = get_response(trigger[these_trials],data_func,trial_duration,optw)\n",
    "        recovery_response = np.mean(recovery_r,0)\n",
    "\n",
    "\n",
    "    spikemat_stim = [spikes[np.where((spikes>trigger[t])&\n",
    "                    (spikes<trigger[t]+trial_duration))[0]\n",
    "                    ]-(trigger[t]) for t in stim_range]\n",
    "    spikemat_post = [spikes[np.where((spikes>trigger[t])&\n",
    "                    (spikes<trigger[t]+trial_duration))[0]\n",
    "                    ]-(trigger[t]) for t in postpairing_range]\n",
    "\n",
    "    ccpost = scipy.signal.correlate(postpairing_response,stiminitial_response,mode='same')/(np.var(postpairing_r)*np.var(stiminitial_r))/len(stiminitial_response)\n",
    "    lags = np.linspace(-(0.5*len(postpairing_response)*dt),(0.5*len(postpairing_response)*dt),len(ccpost))\n",
    "\n",
    "    datamat = pd.DataFrame({'exptname' : thisexpt_name,\n",
    "                  'ntrials_base' : [len(base_range)],\n",
    "                  'ntrials_per_period' : [ntrials_per_period],\n",
    "                  'ntrials_post' : [len(postpairing_trials)],\n",
    "                  'ntrials_post_total' : [len(postpairing_range)],\n",
    "                  'ntrials_recovery' : [len(recovery_trials)],\n",
    "                  'ntrials_pairing' : [len(stim_range)],\n",
    "                  'iei_median' : [np.median(iei)],\n",
    "                  'spkrt_base' : [np.mean(base_spkrt)],\n",
    "                  'spkrt_stimi' : [np.mean(stimi_spkrt)],\n",
    "                  'spkrt_stimf' : [np.mean(stimf_spkrt)],\n",
    "                  'spkrt_post' : [np.mean(post_spkrt)],\n",
    "                  'spkrt_recovery' : [np.mean(recovery_spkrt)],\n",
    "                  'latency_minpost' : [xtime[np.where(postpairing_response==np.min(postpairing_response))][0]],\n",
    "                  'latency_maxstim' : [xtime[np.where(stiminitial_response==np.max(stiminitial_response))][0]],\n",
    "                  'sp_poststim_mean' : [stats.spearmanr(postpairing_response,stiminitial_response,nan_policy='propagate')[0]],\n",
    "                  'sp_recoverstim_mean' : [stats.spearmanr(recovery_response,stiminitial_response,nan_policy='propagate')[0]],\n",
    "                  'sp_basestim_mean' : [stats.spearmanr(baseline_mean,stiminitial_response,nan_policy='propagate')[0]],\n",
    "                  'ccneg' : [np.min(ccpost)],\n",
    "                  'ccneg_lag' : [(dt*(np.argmin(ccpost)-len(postpairing_response)))],\n",
    "                  'ccpos' : [np.max(ccpost)],\n",
    "                  'ccpos_lag' : [(dt*(np.argmax(ccpost)-len(postpairing_response)))],\n",
    "                  'cell' : thisexpt.cell.values[0],\n",
    "                  'condition' : thisexpt.condition.values[0],\n",
    "                  'trigger_type' : thisexpt.trigger_type.values[0],\n",
    "                  'phase' : thisexpt.phase.values[0],\n",
    "                  'optw' : optw,\n",
    "                  'optw_stim' : optw_stim,\n",
    "                  'control' : thisexpt.control.values[0]\n",
    "                    })\n",
    "\n",
    "    path_to_file = Path.cwd()\n",
    "    file_to_open = path_to_file / 'DataMat.csv'\n",
    "    f = open(file_to_open)\n",
    "    dfimport = pd.read_csv(f.name)\n",
    "    dfimport = dfimport.loc[:, ~dfimport.columns.str.contains('^Unnamed')]\n",
    "\n",
    "    dfimport = dfimport.append(datamat,sort=True)\n",
    "    dfimport = dfimport.drop_duplicates('exptname','last')\n",
    "    dfimport.to_csv(f.name)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}