{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "503cb687",
   "metadata": {},
   "source": [
    "# Hands-on session 3: Elephant\n",
    "\n",
    "These exercises build on concepts introduced in Tutorial 3\n",
    "\n",
    "In this exercise, we will use two methods available in the Elephant library, namely the Unitary Event analysis to identify significant spike synchrony in a pair of spiketrains across all correct trials and SPADE method to find recurring patterns in the spike data of trial 1.\n",
    "\n",
    "\n",
    "## Imports and Preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "253c2101-0e2e-4e29-8019-d4b4eb050620",
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip install git+https://github.com/INM-6/viziphant\n",
    "#!pip install nixio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fc1d993",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import quantities as pq\n",
    "import neo.utils\n",
    "import elephant\n",
    "import viziphant"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d1cd9a4",
   "metadata": {},
   "source": [
    "## Data loading\n",
    "\n",
    "\n",
    "1. Load the data `reach_to_grasp_material/i140703-001.nix` file using the `neo.io.NixIO` into a neo Block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2e41fc5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "d86077ac",
   "metadata": {},
   "source": [
    "## Data preparation for Unitary Event (UE) analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a15b721",
   "metadata": {},
   "source": [
    "2. At this step, extract the start events of all trial types that are correct. We will cut data from the trial start `TS-ON` to `2105*pq.ms`. Then, in the call to `neo.utils.add_epoch()`, supply these as correct parameters (Hint: call `neo.utils.add_epoch?` to view function documentation). The result should be an `Epoch` from `TS-ON` to 2105 ms *after* `TS-ON`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fff561d5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "047e6aee",
   "metadata": {},
   "source": [
    "3. Now you will need to find the best spiketrains to use for the UE analysis. These should: \n",
    "- be `'sua'`, \n",
    "- not be classified as `'noise'`,\n",
    "- have a high signal-to-noise ratio (e.g. above 10)\n",
    "- have a high spike count (e.g. above 10000)\n",
    "- not be captured by electrodes rejected due to low, intermediate or high frequency components.\n",
    "\n",
    "Create a new segement that contains the list of selected spiketrains that match these criteria. \n",
    "\n",
    "Hint: You may want to look into `segment.spiketrain[0].annotations` and list comprehensions to take the best spiketrains from all the spiketrains present in the segment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a96e8797",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "235933b7",
   "metadata": {},
   "source": [
    "4. Print the number of spiketrains you have found. You will need a pair of spiketrains for the UE analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dffbd02d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3dfaefbc",
   "metadata": {},
   "source": [
    "5. Now that we have isolated the good spiketrains, cut the segment according to the defined epoch using `neo.utils.cut_segment_by_epoch()`. From the result, prepare the input for the unitary event analysis, which is a list, where the `n`-th element is a list of spiketrains for the `n`-th trial, or, in other words, where the `n`-th element is the `spiketrains` attribute of the `n`-th segment.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ee5a825",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "298d5185",
   "metadata": {},
   "source": [
    "## UE analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05b4a995",
   "metadata": {},
   "source": [
    "The Unitary Event (UE) analysis method is very useful in detecting spiking coincidences between a two or several spike trains recorded simulaneously and determining their significance (Grün et al 2002a, Grün et al 2002b, Grün et al 2010). We are going to apply this method to look at a pair of spike trains from the Reach2Grasp dataset you have been exploring previously. The aim is to look for correlated spiking activity in the clean spike trains during all correct trials of all trial types.\n",
    "\n",
    "The dataset should be a list of trials, containing a list of spike trains, which, in this case, is going to contain only two spike train objects.\n",
    "\n",
    "You may want to consider the following parameters for the Unitary Event analysis:\n",
    "* `bin_size`: The size of bins for discretizing spike trains.\n",
    "* `win_size`: The size of the window of analysis\n",
    "* `win_step`: The size of the window step.\n",
    "* `pattern`: The binary pattern, indicating which neurons should participate in the pattern. E.g. [1, 1] for two out of two neurons, [1, 1, 0] for the first two out of three neurons.\n",
    "* `pattern_hash`: The list of interesting patterns in hash values. Generated from `pattern` using `hash_from_pattern` function.\n",
    "* `significance_level`: The level of significance for p-value evaluation. \n",
    "* `method` Determines the approach for unitary event computation.\n",
    "    * `analytic_TrialByTrial`: calculates the expectancy analytically on each trial, then sums over all of them\n",
    "    * `analytic_TrialAverage`: calculates the expectancy analytically by averaging over trials\n",
    "    * `surrogate_TrialByTrial`: calculates the distribution of expected coincidences by spike time randomization in each different trial, then sums over all of them\n",
    "    \n",
    "6. Run UE analysis using `elephant.unitary_event_analysis.jointJ_window_analysis()` and pass the result together with additional parameters to the viziphant function provided below. You may want to use a `plot_params` dictionary to add unit names and `TS-ON` trigger marker to the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59d18c64",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "44bf3ef6",
   "metadata": {},
   "source": [
    "Optional: Explore the parameters of the functions. Can you add the missing triggers between `TS-ON` and 2105 ms after `TS-ON` to the generated plot?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e7c6c1e",
   "metadata": {},
   "source": [
    "## Data preparation for SPADE"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a42cbb4e",
   "metadata": {},
   "source": [
    "7. In close analogy to the beginning of tutorial 3 and what we have done above, prepare a Neo `Block`, containing one `Segment` of SGHF data of the first correct trial. Name this `Segement` by the variable `trial`, as in tutorial 3.  In contrast to the lecture and the UE analysis above, we will cut data from the trial start `TS-ON` to reward administration indicated by event `RW-ON`. To this end, first find also the `RW-ON` events similar to how we found `TS-ON` in tutorial 3. Then, in the call to `neo.utils.add_epoch()`, supply these as a second event `event2=` *instead* of giving `t_pre=0` and `t_post=2*pq.s`. This will cut from event 1 to event 2, instead of a fixed amount of 2 s around event 1 (as in tutorial 3 or the UE exercise above)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "943d26af",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7564a093",
   "metadata": {},
   "source": [
    "## SPADE analysis\n",
    "\n",
    "Note that patterns in this data are not easily to spot by eye in the rasterplot we developed in tutorial 3.\n",
    "We use SPADE (Torre et al., 2013; Stella et al. 2019) as a technique that does that for us by finding all patterns and checking for each pattern if it occurs more often than expected given its complexity (number of neurons participating) and frequency. Before going directly to the analysis we briefly explain SPADE's most important parameters:\n",
    "\n",
    "- `binsize`: temporal precision of the method. The smaller the binsize is, the more precisely we expect each single pattern to repeat. This raises an important question: which is the temporal precision that are you interested in? It depends on the scientific question! We often use 5ms, based on a number of studies on the minimal neuronal temporal precision.\n",
    "- `winlen`: window length, or maximal length allowed to each pattern, expressed in bin units. SPADE will detect patterns with a temporal duration up to the window length. If winlen=1, then only synchronous patterns are detected. Are you interested in synchronous higher-order correlations? Are you interested in patterns with delays? Note: the higher the winlen parameter is, the more expensive (memory and time) the pattern search is!\n",
    "- `min_spikes` and `min_neu`: minimum number of spikes and minimum number of neurons allowed in a pattern. These parameters are important, too. Do you want to look for pattern with a particular size? Are you interested in patterns with multiple spikes coming from the same neuron?\n",
    "- `n_surr`: number of surrogates used for the statistical testing. If the number of surrogates is set to zero, then all patterns are retrieved by SPADE without further testing. If, instead, the number of surrogates is different from zero (typically we suggest to use a large number - the more the better!) then only the significant patterns are retrieved from the analysis. For simplicity, we will set this parameter to `n_surr=0` in the following.\n",
    "\n",
    "Next steps:\n",
    "8. As in tutorial 3, select only good quality neurons using the annotations (i.e., SUA units, more than 10000 spikes,...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d7cad57",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "10d32abe",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "9. Run the SPADE analysis by executing the code given below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23263b16",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Fix parameters for SPADE\n",
    "binsize = 5 * pq.ms\n",
    "winlen = 1\n",
    "min_spikes = 2\n",
    "n_surr = 0\n",
    "min_neu = 2\n",
    "min_occ = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23f7d6b4",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "result = elephant.spade.spade(\n",
    "        spiketrains=spiketrains, binsize=binsize, winlen=winlen,\n",
    "        min_spikes=min_spikes, n_surr=n_surr, min_neu=min_neu, min_occ=min_occ)\n",
    "patterns = result['patterns']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c3d05a6",
   "metadata": {},
   "source": [
    "10. Let's look at the output in `patterns`. How many patterns were found? How does one pattern look as an output?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "994e3df5",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Print the length of the list of patterns found, and the first pattern as an example\n",
    "len(patterns), patterns[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11523576",
   "metadata": {},
   "source": [
    "11. Display five detected patterns in the data and their statistics by executing the Viziphant function calls given below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5dfda9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "viziphant.patterns.plot_patterns(spiketrains, patterns[:5], circle_sizes=(1, 50, 70))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2a28cb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "viziphant.patterns.plot_patterns_statistics_all(patterns)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "079c0771",
   "metadata": {},
   "source": [
    "# References\n",
    "1. Grün et al (2002a) DOI: 10.1162/089976602753284455\n",
    "2. Grün et al (2002b) DOI: 10.1162/089976602753284464\n",
    "3. Grün et al (2010) DOI: 10.1007/978-1-4419-5675-0_10\n",
    "4. Torre et al (2013) DOI: 10.3389/fncom.2013.00132\n",
    "5. Stella et al (2019) DOI: 10.1016/j.biosystems.2019.104022"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "458adf0e-6b24-4c7a-95a5-2896a7247d0d",
   "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}