{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Auto- and crosscorrelation functions for spike trains\n\nA time bin of size `tbin` is centered around the time difference it\nrepresents. If the correlation function is calculated for `tau` in\n`[-tau_max, tau_max]`, the pair events contributing to the left-most\nbin are those for which `tau` in `[-tau_max-tbin/2, tau_max+tbin/2)` and\nso on.\n\nCorrelate two spike trains with each other assumes spike times to be ordered in\ntime. `tau > 0` means spike2 is later than spike1\n\n* tau_max: maximum time lag in ms correlation function\n* tbin:    bin size\n* spike1:  first spike train [tspike...]\n* spike2:  second spike train [tspike...]\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import nest\nimport numpy as np\n\n\ndef corr_spikes_sorted(spike1, spike2, tbin, tau_max, resolution):\n    tau_max_i = int(tau_max / resolution)\n    tbin_i = int(tbin / resolution)\n\n    cross = np.zeros(int(2 * tau_max_i / tbin_i + 1), 'd')\n\n    j0 = 0\n\n    for spki in spike1:\n        j = j0\n        while j < len(spike2) and spike2[j] - spki < -tau_max_i - tbin_i / 2.0:\n            j += 1\n        j0 = j\n\n        while j < len(spike2) and spike2[j] - spki < tau_max_i + tbin_i / 2.0:\n            cross[int(\n                (spike2[j] - spki + tau_max_i + 0.5 * tbin_i) / tbin_i)] += 1.0\n            j += 1\n\n    return cross\n\n\nnest.ResetKernel()\n\nresolution = 0.1    # Computation step size in ms\nT = 100000.0        # Total duration\ndelta_tau = 10.0\ntau_max = 100.0     # ms correlation window\nt_bin = 10.0        # ms bin size\npc = 0.5\nnu = 100.0\n\nnest.local_num_threads = 1\nnest.resolution = resolution\nnest.overwrite_files = True\nnest.rng_seed = 12345\n\n# Set up network, connect and simulate\nmg = nest.Create('mip_generator')\nmg.set(rate=nu, p_copy=pc)\n\ncd = nest.Create('correlation_detector')\ncd.set(tau_max=tau_max, delta_tau=delta_tau)\n\nsr = nest.Create('spike_recorder', params={'time_in_steps': True})\n\npn1 = nest.Create('parrot_neuron')\npn2 = nest.Create('parrot_neuron')\n\nnest.Connect(mg, pn1)\nnest.Connect(mg, pn2)\nnest.Connect(pn1, sr)\nnest.Connect(pn2, sr)\n\nnest.Connect(pn1, cd, syn_spec={'weight': 1.0, 'receptor_type': 0})\nnest.Connect(pn2, cd, syn_spec={'weight': 1.0, 'receptor_type': 1})\n\nnest.Simulate(T)\n\nn_events_1, n_events_2 = cd.n_events\n\nlmbd1 = (n_events_1 / (T - tau_max)) * 1000.0\nlmbd2 = (n_events_2 / (T - tau_max)) * 1000.0\n\nspikes = sr.get('events', 'senders')\n\nsp1 = spikes[spikes == 4]\nsp2 = spikes[spikes == 5]\n\n# Find crosscorrelation\ncross = corr_spikes_sorted(sp1, sp2, t_bin, tau_max, resolution)\n\nprint(\"Crosscorrelation:\")\nprint(cross)\nprint(\"Sum of crosscorrelation:\")\nprint(sum(cross))"
      ]
    }
  ],
  "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.8.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}