{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Sinusoidal gamma generator example\n\nThis script demonstrates the use of the ``sinusoidal_gamma_generator`` and its\ndifferent parameters and modes. The source code of the model can be found in\n``models/sinusoidal_gamma_generator.h``.\n\nThe script is structured into two parts, each of which generates its own\nfigure. In part 1A, two generators are created with different orders of the\nunderlying gamma process and their resulting PST (Peristiumulus time) and ISI\n(Inter-spike interval) histograms are plotted. Part 1B illustrates the effect\nof the ``individual_spike_trains`` switch. In Part 2, the effects of\ndifferent settings for rate, phase and frequency are demonstrated.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, we import all necessary modules to simulate, analyze and\nplot this example.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import nest\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nnest.ResetKernel()   # in case we run the script multiple times from iPython"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We first create a figure for the plot and set the resolution of NEST.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure()\nnest.resolution = 0.01"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then we create two instances of the ``sinusoidal_gamma_generator`` with two\ndifferent orders of the underlying gamma process using ``Create``. Moreover,\nwe create devices to record firing rates (``multimeter``) and spikes\n(``spike_recorder``) and connect them to the generators using ``Connect``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "num_nodes = 2\ng = nest.Create('sinusoidal_gamma_generator', n=num_nodes,\n                params={'rate': 10000.0,\n                        'amplitude': 5000.0,\n                        'frequency': 10.0,\n                        'phase': 0.0,\n                        'order': [2.0, 10.0]})   # note the syntax for different order parameter of the two nodes\n\nm = nest.Create('multimeter', num_nodes, {'interval': 0.1, 'record_from': ['rate']})\ns = nest.Create('spike_recorder', num_nodes)\n\nnest.Connect(m, g, 'one_to_one')\nnest.Connect(g, s, 'one_to_one')\n\nnest.Simulate(200)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "After simulating, the spikes are extracted from the ``spike_recorder`` and\nplots are created with panels for the PST and ISI histograms.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "colors = ['b', 'g']\n\nfor j in range(num_nodes):\n\n    ev = m[j].events\n    t = ev['times']\n    r = ev['rate']\n\n    spike_times = s[j].events['times']\n    plt.subplot(221)\n    h, e = np.histogram(spike_times, bins=np.arange(0., 201., 5.))\n    plt.plot(t, r, color=colors[j])\n    plt.step(e[:-1], h * 1000 / 5., color=colors[j], where='post')\n    plt.title('PST histogram and firing rates')\n    plt.ylabel('Spikes per second')\n\n    plt.subplot(223)\n    plt.hist(np.diff(spike_times), bins=np.arange(0., 0.505, 0.01),\n             histtype='step', color=colors[j])\n    plt.title('ISI histogram')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The kernel is reset and the number of threads set to 4.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.ResetKernel()\nnest.local_num_threads = 4"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, a ``sinusoidal_gamma_generator`` with ``individual_spike_trains`` set to\n`True` is created and connected to 20 parrot neurons whose spikes are\nrecorded by a spike recorder. After simulating, a raster plot of the spikes\nis created.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "g = nest.Create('sinusoidal_gamma_generator',\n                params={'rate': 100.0, 'amplitude': 50.0,\n                        'frequency': 10.0, 'phase': 0.0, 'order': 3.,\n                        'individual_spike_trains': True})\np = nest.Create('parrot_neuron', 20)\ns = nest.Create('spike_recorder')\n\nnest.Connect(g, p)\nnest.Connect(p, s)\n\nnest.Simulate(200)\nev = s.events\nplt.subplot(222)\nplt.plot(ev['times'], ev['senders'] - min(ev['senders']), 'o')\nplt.ylim([-0.5, 19.5])\nplt.yticks([])\nplt.title('Individual spike trains for each target')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The kernel is reset again and the whole procedure is repeated for a\n``sinusoidal_gamma_generator`` with ``individual_spike_trains`` set to `False`.\nThe plot shows that in this case, all neurons receive the same spike train\nfrom the ``sinusoidal_gamma_generator``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.ResetKernel()\nnest.local_num_threads = 4\n\ng = nest.Create('sinusoidal_gamma_generator',\n                params={'rate': 100.0, 'amplitude': 50.0,\n                        'frequency': 10.0, 'phase': 0.0, 'order': 3.,\n                        'individual_spike_trains': False})\np = nest.Create('parrot_neuron', 20)\ns = nest.Create('spike_recorder')\n\nnest.Connect(g, p)\nnest.Connect(p, s)\n\nnest.Simulate(200)\nev = s.events\nplt.subplot(224)\nplt.plot(ev['times'], ev['senders'] - min(ev['senders']), 'o')\nplt.ylim([-0.5, 19.5])\nplt.yticks([])\nplt.title('One spike train for all targets')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In part 2, multiple generators are created with different settings for rate,\nphase and frequency. First, we define an auxiliary function, which simulates\n`n` generators for `t` ms. After `t/2`, the parameter dictionary of the\ngenerators is changed from initial to after.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def step(t, n, initial, after, seed=1, dt=0.05):\n\n    nest.ResetKernel()\n    nest.resolution = dt\n    nest.rng_seed = seed\n\n    g = nest.Create('sinusoidal_gamma_generator', n, params=initial)\n    sr = nest.Create('spike_recorder')\n    nest.Connect(g, sr)\n    nest.Simulate(t / 2)\n    g.set(after)\n    nest.Simulate(t / 2)\n\n    return sr.events"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This function serves to plot a histogram of the emitted spikes.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def plot_hist(spikes):\n    plt.hist(spikes['times'],\n             bins=np.arange(0., max(spikes['times']) + 1.5, 1.),\n             histtype='step')\n\n\nt = 1000\nn = 1000\ndt = 1.0\nsteps = int(t / dt)\noffset = t / 1000. * 2 * np.pi\n\n\n# We create a figure with a 2x3 grid.\n\n\ngrid = (2, 3)\nfig = plt.figure(figsize=(15, 10))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We simulate a ``sinusoidal_gamma_generator`` with default parameter values,\ni.e. ``ac=0`` and the DC value being changed from 20 to 50 after `t/2` and\nplot the number of spikes per second over time.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.subplot(grid[0], grid[1], 1)\nspikes = step(t, n,\n              {'rate': 20.0},\n              {'rate': 50.0, },\n              seed=123, dt=dt)\nplot_hist(spikes)\nexp = np.ones(int(steps))\nexp[:int(steps / 2)] *= 20\nexp[int(steps / 2):] *= 50\nplt.plot(exp, 'r')\nplt.title('DC rate: 20 -> 50')\nplt.ylabel('Spikes per second')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We simulate a ``sinusoidal_gamma_generator`` with the DC value being changed\nfrom 80 to 40 after `t/2` and plot the number of spikes per second over\ntime.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.subplot(grid[0], grid[1], 2)\nspikes = step(t, n,\n              {'order': 6.0, 'rate': 80.0, 'amplitude': 0.,\n               'frequency': 0., 'phase': 0.},\n              {'order': 6.0, 'rate': 40.0, 'amplitude': 0.,\n               'frequency': 0., 'phase': 0.},\n              seed=123, dt=dt)\nplot_hist(spikes)\nexp = np.ones(int(steps))\nexp[:int(steps / 2)] *= 80\nexp[int(steps / 2):] *= 40\nplt.plot(exp, 'r')\nplt.title('DC rate: 80 -> 40')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we simulate a ``sinusoidal_gamma_generator`` with the AC value being\nchanged from 40 to 20 after `t/2` and plot the number of spikes per\nsecond over time.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.subplot(grid[0], grid[1], 3)\nspikes = step(t, n,\n              {'order': 3.0, 'rate': 40.0, 'amplitude': 40.,\n               'frequency': 10., 'phase': 0.},\n              {'order': 3.0, 'rate': 40.0, 'amplitude': 20.,\n               'frequency': 10., 'phase': 0.},\n              seed=123, dt=dt)\nplot_hist(spikes)\nexp = np.zeros(int(steps))\nexp[:int(steps / 2)] = (40. + 40. * np.sin(np.arange(\n    0, t / 1000. * np.pi * 10, t / 1000. * np.pi * 10. / (steps / 2))))\nexp[int(steps / 2):] = (40. + 20. * np.sin(np.arange(\n    0, t / 1000. * np.pi * 10, t / 1000. * np.pi * 10. / (steps / 2)) + offset))\nplt.plot(exp, 'r')\nplt.title('Rate Modulation: 40 -> 20')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we simulate a ``sinusoidal_gamma_generator`` with a non-zero AC value\nand the DC value being changed from 80 to 40 after `t/2` and plot the\nnumber of spikes per second over time.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.subplot(grid[0], grid[1], 4)\nspikes = step(t, n,\n              {'order': 6.0, 'rate': 20.0, 'amplitude': 20.,\n               'frequency': 10., 'phase': 0.},\n              {'order': 6.0, 'rate': 50.0, 'amplitude': 50.,\n               'frequency': 10., 'phase': 0.},\n              seed=123, dt=dt)\nplot_hist(spikes)\nexp = np.zeros(int(steps))\nexp[:int(steps / 2)] = (20. + 20. * np.sin(np.arange(\n    0, t / 1000. * np.pi * 10, t / 1000. * np.pi * 10. / (steps / 2))))\nexp[int(steps / 2):] = (50. + 50. * np.sin(np.arange(\n    0, t / 1000. * np.pi * 10, t / 1000. * np.pi * 10. / (steps / 2)) + offset))\nplt.plot(exp, 'r')\nplt.title('DC Rate and Rate Modulation: 20 -> 50')\nplt.ylabel('Spikes per second')\nplt.xlabel('Time [ms]')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Simulate a ``sinusoidal_gamma_generator`` with the AC value being\nchanged from 0 to 40 after `t/2` and plot the number of spikes per\nsecond over time.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.subplot(grid[0], grid[1], 5)\nspikes = step(t, n,\n              {'rate': 40.0, },\n              {'amplitude': 40.0, 'frequency': 20.},\n              seed=123, dt=1.)\nplot_hist(spikes)\nexp = np.zeros(int(steps))\nexp[:int(steps / 2)] = 40. * np.ones(int(steps / 2))\nexp[int(steps / 2):] = (40. + 40. * np.sin(np.arange(\n    0, t / 1000. * np.pi * 20, t / 1000. * np.pi * 20. / (steps / 2))))\nplt.plot(exp, 'r')\nplt.title('Rate Modulation: 0 -> 40')\nplt.xlabel('Time [ms]')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Simulate a ``sinusoidal_gamma_generator`` with a phase shift at\n`t/2` and plot the number of spikes per second over time.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Phase shift\nplt.subplot(grid[0], grid[1], 6)\nspikes = step(t, n,\n              {'order': 6.0, 'rate': 60.0, 'amplitude': 60.,\n               'frequency': 10., 'phase': 0.},\n              {'order': 6.0, 'rate': 60.0, 'amplitude': 60.,\n               'frequency': 10., 'phase': 180.},\n              seed=123, dt=1.)\nplot_hist(spikes)\nexp = np.zeros(int(steps))\n\nexp[:int(steps / 2)] = (60. + 60. * np.sin(np.arange(\n    0, t / 1000. * np.pi * 10, t / 1000. * np.pi * 10. / (steps / 2))))\nexp[int(steps / 2):] = (60. + 60. * np.sin(np.arange(\n    0, t / 1000. * np.pi * 10, t / 1000. * np.pi * 10. / (steps / 2)) + offset + np.pi))\nplt.plot(exp, 'r')\nplt.title('Modulation Phase: 0 -> Pi')\nplt.xlabel('Time [ms]')\nplt.show()"
      ]
    }
  ],
  "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
}