{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Rate neuron decision making\n\nA binary decision is implemented in the form of two rate neurons\nengaging in mutual inhibition.\n\nEvidence for each decision is reflected by the mean input\nexperienced by the respective neuron.\nThe activity of each neuron is recorded using multimeter devices.\n\nIt can be observed how noise as well as the difference in evidence\naffects which neuron exhibits larger activity and hence which\ndecision will be made.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import nest\nimport matplotlib.pyplot as plt\nimport numpy"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First,the function ``build_network`` is defined to build the network\nand return the handles of two decision units and the ``multimeter``\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def build_network(sigma, dt):\n    nest.ResetKernel()\n    nest.SetKernelStatus({'resolution': dt, 'use_wfr': False})\n    Params = {'lambda': 0.1, 'sigma': sigma, 'tau': 1., 'rectify_output': True}\n    D1 = nest.Create('lin_rate_ipn', params=Params)\n    D2 = nest.Create('lin_rate_ipn', params=Params)\n\n    nest.Connect(D1, D2, 'all_to_all', {\n        'synapse_model': 'rate_connection_instantaneous', 'weight': -0.2})\n    nest.Connect(D2, D1, 'all_to_all', {\n        'synapse_model': 'rate_connection_instantaneous', 'weight': -0.2})\n\n    mm = nest.Create('multimeter')\n    mm.set(interval=dt, record_from=['rate'])\n    nest.Connect(mm, D1, syn_spec={'delay': dt})\n    nest.Connect(mm, D2, syn_spec={'delay': dt})\n\n    return D1, D2, mm"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The function ``build_network`` takes the noise parameter sigma\nand the time resolution as arguments.\nFirst, the Kernel is reset and the ``use_wfr`` (waveform-relaxation)\nis set to false while the resolution is set to the specified value\n`dt`.  Two rate neurons with linear activation functions are created\nand the handle is stored in the variables `D1` and `D2`. The output\nof both decision units is rectified at zero.  The two decisions\nunits are coupled via mutual inhibition.  Next the multimeter is\ncreated and the handle stored in mm and the option ``record_from``\nis set. The multimeter is then connected to the two units in order\nto 'observe' them.  The ``Connect`` function takes the handles as\ninput.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The decision making process is simulated for three different levels\nof noise and three differences in evidence for a given decision. The\nactivity of both decision units is plotted for each scenario.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig_size = [14, 8]\nfig_rows = 3\nfig_cols = 3\nfig_plots = fig_rows * fig_cols\nface = 'white'\nedge = 'white'\n\nax = [None] * fig_plots\nfig = plt.figure(facecolor=face, edgecolor=edge, figsize=fig_size)\n\ndt = 1e-3\nsigma = [0.0, 0.1, 0.2]\ndE = [0.0, 0.004, 0.008]\nT = numpy.linspace(0, 200, 200 / dt - 1)\nfor i in range(9):\n    c = i % 3\n    r = int(i / 3)\n    D1, D2, mm = build_network(sigma[r], dt)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First using build_network the network is build and the handles of\nthe decision units and the multimeter are stored in `D1`, `D2` and `mm`\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.Simulate(100.0)\n    D1.mu = 1. + dE[c]\n    D2.mu = 1. - dE[c]\n    nest.Simulate(100.0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The network is simulated using ``Simulate``, which takes the desired\nsimulation time in milliseconds and advances the network state by\nthis amount of time. After an initial period in the absence of evidence\nfor either decision, evidence is given by changing the state of each\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "senders = mm.get('events', 'senders')\n    voltages = mm.get('events', 'rate')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The activity values ('voltages') are read out by the multimeter\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax[i] = fig.add_subplot(fig_rows, fig_cols, i + 1)\n    ax[i].plot(T, voltages[numpy.where(senders == D1.global_id)],\n               'b', linewidth=2, label=\"D1\")\n    ax[i].plot(T, voltages[numpy.where(senders == D2.global_id)],\n               'r', linewidth=2, label=\"D2\")\n    ax[i].set_ylim([-.5, 12.])\n    ax[i].get_xaxis().set_ticks([])\n    ax[i].get_yaxis().set_ticks([])\n    if c == 0:\n        ax[i].set_ylabel(r\"activity ($\\sigma=%.1f$) \" % (sigma[r]))\n        ax[i].get_yaxis().set_ticks([0, 3, 6, 9, 12])\n\n    if r == 0:\n        ax[i].set_title(r\"$\\Delta E=%.3f$ \" % (dE[c]))\n        if c == 2:\n            plt.legend(loc=0)\n    if r == 2:\n        ax[i].get_xaxis().set_ticks([0, 50, 100, 150, 200])\n        ax[i].set_xlabel('time (ms)')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The activity of the two units is plotted in each scenario.\n\nIn the absence of noise, the network will not make a decision if evidence\nfor both choices is equal. With noise, this symmetry can be broken and a\ndecision wil be taken despite identical evidence.\n\nAs evidence for `D1` relative to `D2` increases, it becomes more likely that\nthe corresponding decision will be taken. For small differences in the\nevidence for the two decisions, noise can lead to the 'wrong' decision.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.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.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}