ball-and-stick-4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the final part of the series where we build a multicompartment cell and evolve it into a network of cells running on a parallel machine (which is basically all computers made within the last decade). On this page, we translate the classes we have previously constructed so that they operate in either a parallel or serial mode."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Note:</b> If you followed the <a href=\"https://neuron.yale.edu/ftp/neuron/2019umn/neuron-quickstart.pdf\">Quick Start</a> instructions on the NEURON website, you should have everything you need to run parallel simulations. If not, if you do not already have an MPI installation, go to that link and follow the instructions in \"Step 4: Install MPI to enable parallel simulations.\" If you compiled NEURON yourself instead of using an installer (this is almost-never necessary), this part of the tutorial requires you to have used the `--with-paranrn` flag at configure time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Parallel communication in NEURON</h2>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parallel communication takes place via logical events in network connection objects known as <a href=\"https://neuron.yale.edu/neuron/static/py_doc/modelspec/programmatic/network/netcon.html#NetCon\">NetCon</a>. NetCon sources are threshold detectors. They monitor some variable, say the membrane potential of a cell, and when the variable reaches some threshold, it triggers an event sent to the targets. Targets are frequently synapses on other cells. When they receive the event, they activate.\n",
    "\n",
    "In a parallel context across several machines, communication between hosts can be computationally inefficient when the frequency of events is high and when the message being sent is large. NEURON uses an efficient priority queueing mechanism to deliver events to targets after the delay specified by the NetCon. The message passed is succinct. It is an integer, the unique global identifier (gid) of the source. The following two figures illustrate these ideas and come from <a href=\"http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2430920/\">Hines M.L. and Carnevale N.T, Translating network models to parallel hardware in NEURON, Journal of Neuroscience Methods 169 (2008) 425–455</a>. Users should also consult the <a href=\"https://neuron.yale.edu/neuron/static/py_doc/modelspec/programmatic/network/parcon.html#ParallelContext\">ParallelContext</a> reference."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><img src=\"https://neuron.yale.edu/neuron/static/docs/neuronpython/_images/parfig2.1.png\"/></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><img src=\"https://neuron.yale.edu/neuron/static/docs/neuronpython/_images/parfig2.2.png\"/></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main step involved in making a parallel implementation is to assign the global identifiers across the various hosts. Care should also be taken to assign cells to the various hosts such that the system is load balanced. For example, in a network with computationally complex and simple cells, several simple cells may be assigned to a host while few complex cells may be assigned to another host."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Test MPI and Parallel NEURON</h2>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again, parallel NEURON requires MPI support. If this is your first time using it in a while, you should test your computer setup."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a file called `testmpi.py`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron import h\n",
    "h.nrnmpi_init()       # initialize MPI\n",
    "pc = h.ParallelContext()\n",
    "print('I am {} of {}'.format(pc.id(), pc.nhost()))\n",
    "h.quit()              # necessary to avoid a warning message on parallel exit on some systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the script in parallel via e.g. `mpiexec -n 4 python testmpi.py` from the command line in a terminal. You should see output resembling:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<pre>\n",
    "numprocs=4\n",
    "I am 1 of 4\n",
    "I am 2 of 4\n",
    "I am 3 of 4\n",
    "I am 0 of 4\n",
    "</pre>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These could appear in any order since in theory they are running simultaneously but must print out in some order. If instead you see four processes claiming to be 0 of 1, then your copy of NEURON was not compiled with support for parallel simulation. Reconfigure with the `--with-paranrn` flag, recompile, and try again. If you get an error saying that `mpiexec` is an unknown command, then MPI is either not installed or not on your PATH; correct your MPI setup and try again."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Back to the example</h2>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Begin by downloading ballandstick.py into your working directory. This is equivalent to the classes we created in the previous part of the tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will construct a `ring.py` based on the previous `Ring` class. Changes are indicated with `###`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron import h\n",
    "from ballandstick import BallAndStick\n",
    "\n",
    "### MPI must be initialized before we create a ParallelContext object\n",
    "h.nrnmpi_init()\n",
    "pc = h.ParallelContext()\n",
    "\n",
    "class Ring:\n",
    "    \"\"\"A network of *N* ball-and-stick cells where cell n makes an\n",
    "    excitatory synapse onto cell n + 1 and the last, Nth cell in the\n",
    "    network projects to the first cell.\n",
    "    \"\"\"\n",
    "    def __init__(self, N=5, stim_w=0.04, stim_t=9, stim_delay=1, syn_w=0.01, syn_delay=5, r=50):\n",
    "        \"\"\"\n",
    "        :param N: Number of cells.\n",
    "        :param stim_w: Weight of the stimulus\n",
    "        :param stim_t: time of the stimulus (in ms)\n",
    "        :param stim_delay: delay of the stimulus (in ms)\n",
    "        :param syn_w: Synaptic weight\n",
    "        :param syn_delay: Delay of the synapse\n",
    "        :param r: radius of the network\n",
    "        \"\"\" \n",
    "        self._N = N\n",
    "        self.set_gids()                   ### assign gids to processors\n",
    "        self._syn_w = syn_w\n",
    "        self._syn_delay = syn_delay\n",
    "        self._create_cells(r)             ### changed to use self._N instead of passing in N\n",
    "        self._connect_cells()\n",
    "        ### the 0th cell only exists on one process... that's the only one that gets a netstim\n",
    "        if pc.gid_exists(0):\n",
    "            self._netstim = h.NetStim()\n",
    "            self._netstim.number = 1\n",
    "            self._netstim.start = stim_t\n",
    "            self._nc = h.NetCon(self._netstim, pc.gid2cell(0).syn)   ### grab cell with gid==0 wherever it exists\n",
    "            self._nc.delay = stim_delay\n",
    "            self._nc.weight[0] = stim_w\n",
    "    \n",
    "    def set_gids(self):\n",
    "        \"\"\"Set the gidlist on this host.\"\"\"\n",
    "        #### Round-robin counting.\n",
    "        #### Each host has an id from 0 to pc.nhost() - 1.\n",
    "        self.gidlist = list(range(pc.id(), self._N, pc.nhost()))\n",
    "        for gid in self.gidlist:\n",
    "            pc.set_gid2node(gid, pc.id())\n",
    "    \n",
    "    def _create_cells(self, r):\n",
    "        self.cells = []\n",
    "        for i in self.gidlist:    ### only create the cells that exist on this host\n",
    "            theta = i * 2 * h.PI / self._N\n",
    "            self.cells.append(BallAndStick(i, h.cos(theta) * r, h.sin(theta) * r, 0, theta))\n",
    "        ### associate the cell with this host and gid\n",
    "        for cell in self.cells:\n",
    "            pc.cell(cell._gid, cell._spike_detector)\n",
    "\n",
    "    def _connect_cells(self):\n",
    "        ### this method is different because we now must use ids instead of objects\n",
    "        for target in self.cells:\n",
    "            source_gid = (target._gid - 1 + self._N) % self._N\n",
    "            nc = pc.gid_connect(source_gid, target.syn)\n",
    "            nc.weight[0] = self._syn_w\n",
    "            nc.delay = self._syn_delay\n",
    "            target._ncs.append(nc)        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The call to `h.nrnmpi_init()` must happen before any use of the <a href=\"https://www.neuron.yale.edu/neuron/static/py_doc/modelspec/programmatic/network/parcon.html\">ParallelContext</a> class -- which forms a key part of any NEURON parallel simulation.\n",
    "\n",
    "The only conceptually new method here is the `set_gids` method where each process specifies which cells it will simulate. Here we use what is known as a round-robin approach, where the `pc.id()`th process starts at `pc.id()` and skips by however many processes are running (`pc.nhost`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In `_create_cells`, we now associate the cell `gid` with the NetCon `_spike_detector`. This allows the `_connect_cells` to make connections based on gids instead of objects, using `pc.gid_connect`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Test the model</h2>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a basic test, `test_ring1.py` that loads the `Ring` class and plots cell 0's membrane potential timeseries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron import h\n",
    "from neuron.units import ms, mV\n",
    "import matplotlib.pyplot as plt\n",
    "from ring import Ring\n",
    "\n",
    "cell_to_plot = 0\n",
    "\n",
    "ring = Ring()\n",
    "\n",
    "pc = h.ParallelContext()\n",
    "pc.set_maxstep(10 * ms)\n",
    "\n",
    "t = h.Vector().record(h._ref_t)\n",
    "h.finitialize(-65 * mV)\n",
    "pc.psolve(100 * ms)\n",
    "\n",
    "if pc.gid_exists(cell_to_plot):\n",
    "    plt.figure()\n",
    "    plt.plot(t, pc.gid2cell(cell_to_plot).soma_v)\n",
    "    plt.show()\n",
    "\n",
    "pc.barrier()\n",
    "pc.done()\n",
    "h.quit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code above should look very familiar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The conceptually new pieces are:\n",
    "<ul>\n",
    "    <li><tt>pc.set_maxstep(10 * ms)</tt> -- sets an upper bound on how far MPI can simulate without communicating, here a simulated 10 ms. This <i>must</i> be called before attempting a parallel simulation.</li>\n",
    "    <li><tt>pc.psolve(100 * ms)</tt> -- a parallel version of <tt>h.continuerun</tt>, but does not support updating NEURON graphics during the simulation.</li>\n",
    "    <li><tt>pc.gid_exists</tt>  -- only the process that owns the specified cell should make the plot.</li>\n",
    "    <li><tt>pc.gid2cell</tt> -- lookup a cell by gid</li>\n",
    "    <li><tt>pc.barrier()</tt> -- wait until all processes reach this point; used to make sure processes don't shut down before the graph is closed</li>\n",
    "</ul>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start by testing this without using MPI:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<tt>python test_ring1.py</tt>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And once you are satisfied that works, you can try MPI, e.g."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<tt>mpiexec -n 2 python test_ring1.py</tt>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Gathering spikes</h2>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our above test runs the simulation successfully, but in the end, no single process knows when all the spikes occurred. There are a number of ways to deal with this: one solution is to have each process write its data to a file. Instead, we will use <a href=\"https://www.neuron.yale.edu/neuron/static/py_doc/modelspec/programmatic/network/parcon.html#ParallelContext.py_alltoall\">pc.py_alltoall</a> to send all the data to node 0, at which point node 0 can plot the raster, save data, or whatever."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Store this in `test_ring2.py`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron import h\n",
    "from neuron.units import ms, mV\n",
    "import matplotlib.pyplot as plt\n",
    "from ring import Ring\n",
    "\n",
    "ring = Ring()\n",
    "\n",
    "pc = h.ParallelContext()\n",
    "pc.set_maxstep(10 * ms)\n",
    "\n",
    "t = h.Vector().record(h._ref_t)\n",
    "h.finitialize(-65 * mV)\n",
    "pc.psolve(100 * ms)\n",
    "\n",
    "# send all spike time data to node 0\n",
    "local_data = {cell._gid: list(cell.spike_times) for cell in ring.cells}\n",
    "all_data = pc.py_alltoall([local_data] + [None] * (pc.nhost() - 1))\n",
    "\n",
    "if pc.id() == 0:\n",
    "    # combine the data from the various processes\n",
    "    data = {}\n",
    "    for process_data in all_data:\n",
    "        data.update(process_data)\n",
    "    # plot it\n",
    "    plt.figure()\n",
    "    for i, spike_times in data.items():\n",
    "        plt.vlines(spike_times, i + 0.5, i + 1.5)\n",
    "    plt.show()\n",
    "\n",
    "pc.barrier()\n",
    "pc.done()\n",
    "h.quit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running this via e.g. `mpiexec -n 2 python test_ring2.py` displays the familiar raster plot. If you are wondering why node 0 was the one chosen to make the plot, it is because that is the only node that is guaranteed to exist (nothing else exists if there is only one process being used for the simulation)."
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}