advanced_python_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Introduction\n",
    "\n",
    "This module will cover intermediate and advanced Python topics -- installing and creating packages, data loading and saving, visualization tools, multiprocessing, functions and classes, and exceptions and debugging.\n",
    "\n",
    "\n",
    "## 1.1 Anaconda & Spyder\n",
    "\n",
    "The great thing about Python is that it has so many packages (also called libraries or modules) written for it. But they're a pain to install yourself. The easiest way to install a \"scientific Python\" installation on your computer is with Anaconda:\n",
    "\n",
    "https://www.anaconda.com/distribution/#download-section\n",
    "\n",
    "Everything you need is preinstalled on the OCNC virtual machines, so if you're using that, you don't need to do anything. However, you are encouraged to set it up on your local computer, since you probably won't keep using the virtual machine after OCNC!\n",
    "\n",
    "Anaconda includes an \"interactive development environment\" (IDE) called Spyder, which provides a Matlab-like interface to Python. For big projects, it's much more flexible than a Jupyter notebook.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Exercise 1.1:</b><br>\n",
    "    Start Spyder.\n",
    "</div>\n",
    "\n",
    "\n",
    "## 1.2 Pylab\n",
    "\n",
    "One useful way of using Python is by including the following at the top of your file:\n",
    "\n",
    "`import pylab as pl`\n",
    "\n",
    "or\n",
    "\n",
    "`from pylab import *`\n",
    "\n",
    "Pylab is a package that contains all the most commonly used scientific functions. The second line is basically just shorthand for\n",
    "\n",
    "```python\n",
    "from numpy import * # numerical functions like arrays, log, exp\n",
    "from scipy import * # scientific functions like interpolation\n",
    "from matplotlib.pyplot import * # plotting functions\n",
    "```\n",
    "\n",
    "Note that `import pylab as pl` is more \"pythonic\" (and is strongly encouraged!), but `from pylab import *` lets you have the most Matlab-like environment. For example, the following is valid Python _or_ Matlab code!\n",
    "\n",
    "```python\n",
    "from pylab import * # comment out in matlab\n",
    "randomnumbers = randn(50,10);\n",
    "cumulativenumbers = cumsum(abs(randomnumbers));\n",
    "plot(cumulativenumbers);\n",
    "```\n",
    "\n",
    "(Caveat: these give slightly different results though!)\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Exercise 1.2:</b><br>\n",
    "    (A) Run the code block in both Matlab and Python. <br>\n",
    "    (B) Read the Python documentation on the \"cumsum\" function (\"help(cumsum)\") and figure out how to make Python match Matlab.\n",
    "</div>\n",
    "\n",
    "\n",
    "## 1.3 Installing packages\n",
    "\n",
    "Numpy, Pylab, and so forth are examples of _packages_. The easiest way to install packages in Python is with `pip` (package installer for python). This searches http://pypi.org for a package with that name, downloads it, and installs it, while you sit back and sip tea. Usually you would just do `pip install the_package_i_want`, but on the OCNC virtual machines you need a slightly different syntax: `python3 -m pip install the_package_i_want`.\n",
    "\n",
    "We will use Sciris a bit later, but as a small preview, Sciris includes a data structure called an `odict`, which is another kind of ordered dictionary that also behaves like a list wherever possible:\n",
    "\n",
    "```python\n",
    "import pylab as pl\n",
    "import sciris as sc\n",
    "\n",
    "citations_dict = {'kenji doya':6363, 'erik de schutter':6366, 'eve marder': 27480}\n",
    "citations_odict = sc.odict({'kenji doya':6363, 'erik de schutter':6366, 'eve marder': 27480})\n",
    "\n",
    "total_citations = pl.array(citations_dict.values()).sum()\n",
    "total_citations = citations_odict[:].sum()\n",
    "\n",
    "alphabetized_dict = {key:citations_dict[key] for key in sorted(list(citations_dict.keys()))}\n",
    "alphabetized_odict = citations_odict.sorted()\n",
    "\n",
    "first_entry = alphabetized_dict[list(alphabetized_dict.keys())[0]]\n",
    "first_entry = alphabetized_odict[0]\n",
    "```\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Exercise 1.3:</b><br>\n",
    "    (A) If you're using the virtual machine: type \"python3 -m pip install sciris\" to install the Sciris package. Otherwise: type \"pip install sciris\".<br>\n",
    "    (B) Run the example code above. Extra credit (warning, difficult!): you can sort an odict by value with \"citations_odict.sorted(sortby='values')\". How would you create a regular dictionary sorted by value?\n",
    "</div>\n",
    "\n",
    "\n",
    "## 1.4 Creating packages\n",
    "\n",
    "If you're working on a big project, you don't want everything in one file. To use functions from one file in another, you can import that file, just as you would `numpy` or `pylab`. For example, if you call your file `myfunction.py`, and inside that file you have defined `def myfunction()`, you would use it as follows:\n",
    "\n",
    "```python\n",
    "import myfunction as myf\n",
    "output = myf.myfunction()\n",
    "```\n",
    "\n",
    "Note: this syntax of `nameofthing.nameofthing()` is frustratingly common in Python, but you will learn to see beauty in it. (Maybe.)\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Exercise 1.4:</b><br>\n",
    "    (A) Convert the block of code from Exercise 1.2 into a function that returns `cumulativenumbers` (but doesn't plot), and save this as a separate file (e.g. \"exercise1.py\".<br>\n",
    "    (B) Create a new file (e.g. \"exercise1plot.py\") that imports this function, calls it, and then plots the results.\n",
    "</div>\n",
    "\n",
    "\n",
    "## 1.5 Classes and methods\n",
    "\n",
    "There are two basic approaches to programming: functional programming and object-oriented programming (OOP). One isn't \"better\" than the other, it depends on what you're trying to do. In general, you have data, and functions that do things on the data. In functional programming, you keep your data together in one place, and you keep your functions together in another place. In OOP, you keep data and functions together in an _object_.\n",
    "\n",
    "Let's make it more concrete. Let's say you have you have a worm, a pigeon, and a cat. If you care about their _shared_ properties and functions that apply to all of them, you want to take a functional programming approach, for example:\n",
    "\n",
    "```python\n",
    "# Functional approach\n",
    "\n",
    "masses = {'worm':0.002, 'pigeon':0.086, 'cat':7.3} # in kg\n",
    "velocities = {'worm':0.01, 'pigeon':9.6, 'cat':4.4} # in m/s\n",
    "\n",
    "def energy(mass, velocity):\n",
    "    return 0.5*mass*velocity**2\n",
    "\n",
    "print('Maximum kinetic energy for:')\n",
    "for animal in masses.keys():\n",
    "    KE = energy(masses[animal], velocities[animal])\n",
    "    print(f'{animal}: {KE}')\n",
    "```\n",
    "\n",
    "However, it might also make more sense to treat each animal as an _object_ with properties/attributes (data) and methods (functions). That would look like this:\n",
    "\n",
    "```python\n",
    "# Object-oriented\n",
    "    \n",
    "class Animal():\n",
    "    \n",
    "    def __init__(self, name, mass, velocity):\n",
    "        self.name = name\n",
    "        self.mass = mass\n",
    "        self.velocity = velocity\n",
    "        \n",
    "    def energy(self):\n",
    "        return 0.5*self.mass*self.velocity**2\n",
    "    \n",
    "    \n",
    "worm   = Animal(name='worm',   mass=0.002, velocity=0.01) # shorthand: Animal('worm', 0.002, 0.01)\n",
    "pigeon = Animal(name='pigeon', mass=0.086, velocity=9.6)\n",
    "cat    = Animal(name='cat',    mass=7.3,   velocity=4.4)\n",
    "\n",
    "print('Maximum kinetic energy for:')\n",
    "for animal in [worm, pigeon, cat]:\n",
    "    print(f'{animal.name}: {animal.energy()}')\n",
    "```\n",
    "\n",
    "Note that these two approaches are **exactly equivalent**, it's just a matter of which one works better for you. OOP is almost always more work to set up initially, but can make life easier later. That said, for most of you most of the time, functional programming probably makes the most sense.\n",
    "\n",
    "\n",
    "## 1.6 Multiprocessing\n",
    "\n",
    "Multiprocessing means running on multiple cores of your computer simultaneously. If you have a big simulation, for example, then if you can run it across 8 cores, it should finish 8 times faster.\n",
    "\n",
    "This code gives an example of a computationally intensive function run in an \"embarrassingly parallel\" way (each run is completely independent). It should take 10-20 seconds in serial (depending on your computer), and be 2-8 times faster in parallel.\n",
    "\n",
    "```python\n",
    "# Imports\n",
    "import pylab as pl\n",
    "import sciris as sc\n",
    "import multiprocessing as mp\n",
    "\n",
    "# Set parameters\n",
    "xmin = 0\n",
    "xmax = 10\n",
    "npts = 200\n",
    "std = 0.1\n",
    "repeats = 100000\n",
    "x = pl.linspace(xmin, xmax, npts)\n",
    "run_series = True\n",
    "run_manual = True\n",
    "run_sciris = True\n",
    "do_plot = True\n",
    "\n",
    "# Define the random wave generator\n",
    "def randgen(std):\n",
    "    pl.seed(int(std*100))\n",
    "    a = 0.2*pl.cos(x)*pl.sqrt(repeats)\n",
    "    b = pl.zeros(npts)\n",
    "    for r in range(repeats):\n",
    "        b += (0.5-pl.rand(npts))*std\n",
    "    return a+b\n",
    "\n",
    "\n",
    "# Set arguments for the parallel run\n",
    "noisevals = pl.linspace(0,1,21)\n",
    "\n",
    "# Run in series\n",
    "if run_series:\n",
    "    sc.tic()\n",
    "    data = []\n",
    "    for noiseval in noisevals:\n",
    "        output = randgen(noiseval)\n",
    "        data.append(output)\n",
    "    sc.toc()\n",
    "\n",
    "# Run in parallel -- manual way\n",
    "if run_manual:\n",
    "    sc.tic() \n",
    "    multipool = mp.Pool(processes=mp.cpu_count())\n",
    "    data = multipool.map(randgen, noisevals)\n",
    "    multipool.close()\n",
    "    multipool.join()\n",
    "    sc.toc()\n",
    "\n",
    "# Run in parallel -- easy way\n",
    "if run_sciris:\n",
    "    sc.tic() \n",
    "    data = sc.parallelize(randgen, noisevals)\n",
    "    sc.toc()\n",
    "\n",
    "# Create 3D plot\n",
    "if do_plot:\n",
    "    sc.surf3d(pl.array(data))\n",
    "```\n",
    "\n",
    "## 1.7 Debugging\n",
    "\n",
    "If you write code, you will also write bugs. Thankfully, debugging is easy in Spyder.\n",
    "\n",
    "Consider this truly awful code:\n",
    "\n",
    "```python\n",
    "def buggy_function1(vector, multiplier):\n",
    "    output = {'max': max(vector*multiplier)}\n",
    "    return output\n",
    "\n",
    "def buggy_function2():\n",
    "    calculation = buggy_function1(vector=[3,4,5], multiplier=3) # should return 15, right?\n",
    "    output = calculation['maximum']\n",
    "    return output\n",
    "\n",
    "output = buggy_function2()\n",
    "print(output)\n",
    "```\n",
    "\n",
    "There are two main kinds of debugging: figuring out why something crashed, and figuring out why something isn't as you expected. For the first, run your code, and after it crashes, type `%debug` in the console; that will put you in your program on the line it crashed. For the second, double-click on the line number you want to pause at; this will put in a _breakpoint_. Then, instead of running your code the normal way (F5), run it in debug mode (Ctrl+F5, or from the menu). It should then stop execution at the breakpoint.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Exercise 1.7:</b><br>\n",
    "    (A) Using the \"%debug\" command, find and fix the bug in buggy_function2().<br>\n",
    "    (B) Using a breakpoint on line 3, find and fix the bug in buggy_function1().\n",
    "</div>\n",
    "\n",
    "\n",
    "## 1.8 Python 2 vs. Python 3\n",
    "\n",
    "Don't use Python 2.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 2. Data saving and loading\n",
    "\n",
    "One of the great things about Python is that it can load/save data in just about any format. Here we will just look at a couple.\n",
    "\n",
    "Except for large files (say >10 MB), it's often easiest just to read/write data files as text. Files in a plain text format can be read by anything and written by anything (more or less).\n",
    "\n",
    "## 2.1 With numpy\n",
    "\n",
    "First, make sure you have the file `data/Okinawa_weather_data.csv`. You would hope you can do something simple, like:\n",
    "\n",
    "```python\n",
    "import pylab as pl\n",
    "filename = 'data/Okinawa_weather_data.csv'\n",
    "data = pl.loadtxt(filename)\n",
    "```\n",
    "\n",
    "But that doesn't work (although it almost does).\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Exercise 2.1A:</b><br>\n",
    "    By reading the documentation for loadtxt(), fix that code so the file will load (hint: you need to add two keyword arguments).\n",
    "</div>\n",
    "\n",
    "Now that we've loaded the data, we want to _do_ something with it. We can plot the 5 variables of interest with:\n",
    "\n",
    "```python\n",
    "pl.figure() # not strictly necessary, python will create a figure automatically\n",
    "pl.plot(data[:,5:10])\n",
    "```\n",
    "\n",
    "We can also for example plot temperature vs. humidity:\n",
    "\n",
    "```python\n",
    "pl.figure() # this is necessary, since otherwise will plot over the other figure\n",
    "pl.scatter(data[:,5],data[:,6])\n",
    "pl.xlabel('temperature (deg C)')\n",
    "pl.ylabel('rel. humidity')\n",
    "```\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Exercise 2.1B:</b><br>\n",
    "    Using the code above, plot the data. For the first figure, modify the plot command so the line width is 3.\n",
    "</div>\n",
    "\n",
    "\n",
    "## 2.2 With Pandas\n",
    "\n",
    "If you are regularly working with tabular files like this, you should have a look at the [`pandas` library](http://pandas.pydata.org/), a powerful library built exactly for this use case:\n",
    "\n",
    "```python\n",
    "import pandas as pd\n",
    "pd_data = pd.read_csv('data/Okinawa_weather_data.csv', sep=';')\n",
    "print(pd_data)\n",
    "```\n",
    "\n",
    "In contrast with `numpy`, `pandas` reads data into a **dataframe**. This is an extremely powerful data structure with a steep learning curve. A full explanation is beyond the scope of this tutorial, but see the [documentation](http://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.html) and [a good tutorial](https://www.datacamp.com/community/tutorials/pandas-tutorial-dataframe-python) for more information.\n",
    "\n",
    "To give a brief example of the possibilities of dataframes, to group the data by days and plot the mean and standard deviation of temperature over a day, run the following:\n",
    "\n",
    "```python\n",
    "grouped_data = pd_data.groupby('Day').aggregate(['mean', 'std'])\n",
    "temp = grouped_data['Temperature  [2 m above gnd]']\n",
    "pl.plot(temp.index, temp['mean'])\n",
    "pl.fill_between(temp.index, temp['mean'] - 2*temp['std'], temp['mean'] + 2*temp['std'], alpha=0.2)\n",
    "pl.xlabel('Day')\n",
    "pl.ylabel('Temperature (C)')\n",
    "```\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Exercise 2.2:</b><br>\n",
    "    Modify the code above for wind speed instead of temperature (hint: you can list the column names in a dataframe with \"dataframe.columns\").\n",
    "</div>\n",
    "\n",
    "\n",
    "\n",
    "## 2.3 Saving/loading objects\n",
    "\n",
    "Sometimes, we want to save more than just text -- for example, the animals we created above. The details are complicated, but we can save arbitrary objects easily using Pylab:\n",
    "\n",
    "```python\n",
    "import sciris as sc\n",
    "import pylab as pl\n",
    "\n",
    "class Animal(sc.prettyobj):\n",
    "\n",
    "    def __init__(self, name, mass, velocity):\n",
    "        self.name = name\n",
    "        self.mass = mass\n",
    "        self.velocity = velocity\n",
    "\n",
    "    def energy(self):\n",
    "        return 0.5*self.mass*self.velocity**2\n",
    "\n",
    "# Create animals\n",
    "worm   = Animal(name='worm',   mass=0.002, velocity=0.01) # shorthand: Animal('worm', 0.002, 0.01)\n",
    "pigeon = Animal(name='pigeon', mass=0.086, velocity=9.6)\n",
    "cat    = Animal(name='cat',    mass=7.3,   velocity=4.4)\n",
    "animals = [worm, pigeon, cat]\n",
    "\n",
    "# Save and load animals\n",
    "pl.save('animals.npy', animals) # file must end .npy (for numpy)\n",
    "saved_animals = pl.load('animals.npy')\n",
    "\n",
    "# Check it worked\n",
    "print(saved_animals[0])\n",
    "```\n",
    "\n",
    "\n",
    "## 2.4 To/from Matlab\n",
    "\n",
    "The `scipy` library comes with support to load matlab files in its `io` package:\n",
    "\n",
    "```python\n",
    "import scipy.io as sio\n",
    "mat_data = sio.loadmat('data/matlab_old.mat')\n",
    "```\n",
    "\n",
    "The data is loaded as a Python dictionary, with the variable names as the keys and the variable values as the, well, values. Note that Matlab does not have the concept of one-dimensional arrays, vectors are either row- or column-vectors and therefore two-dimensional (with one of the dimensions having size one):\n",
    "\n",
    "```python\n",
    "print(mat_data['X'])\n",
    "print(mat_data['X'].shape)\n",
    "```\n",
    "\n",
    "If you want to avoid this, you can automatically \"squeeze\" the dimensions when loading:\n",
    "\n",
    "```python\n",
    "mat_data = sio.loadmat('data/matlab_old.mat', squeeze_me=True)\n",
    "print(mat_data['X'].shape)\n",
    "```\n",
    "\n",
    "You can also access cell arrays or structures:\n",
    "\n",
    "```python\n",
    "print('cell array:', mat_data['c'])\n",
    "print('name elements of structure s:', mat_data['s']['name'])\n",
    "```\n",
    "\n",
    "To save data in the Matlab format, you can use `scipy`'s `savemat` function:\n",
    "\n",
    "```python\n",
    "x = np.arange(10)  # some data\n",
    "y = np.random.rand(10) # some more random data\n",
    "sio.savemat('file.mat', {'x': x}) # variables have to be given as a dictionary of variable names and values\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Advanced plotting\n",
    "\n",
    "There are many ways of plotting in Python. Matplotlib is the de facto standard, and what we'll focus on here, but for some things, other packages might be better.\n",
    "\n",
    "## 3.1 Comparison of plotting libraries\n",
    "\n",
    "The huge variety of plotting libraries is covered by Jake Vanderplas in his 2017 [talk](https://github.com/rougier/python-visualization-landscape). Briefly:\n",
    "\n",
    "**Matplotlib:** By far the most widely used, should be your go-to unless you have specific requirements otherwise.\n",
    "\n",
    "```python\n",
    "import pylab as pl\n",
    "pl.plot([1,4,3,2,5,8])\n",
    "pl.title('Simple example')\n",
    "pl.xlabel('x')\n",
    "pl.ylabel('y')\n",
    "```\n",
    "\n",
    "**Bokeh:** Created by the team behind Anaconda, focuses on visualizing big data interactively. Renders in a browser rather than in Python. A bit more complicated than Matplotlib, but fairly similar.\n",
    "\n",
    "```python\n",
    "import bokeh.plotting as bop\n",
    "p = bop.figure(title=\"Simple example\", x_axis_label='x', y_axis_label='y')\n",
    "p.line(x=range(6), y=[1,4,3,2,5,8])\n",
    "bop.show(p)\n",
    "```\n",
    "\n",
    "**Altair:** Based on the [Grammar of Graphics](https://www.springer.com/gp/book/9780387245447), steep initial learning curve but then allows complicated statistical data to be plotted intuitively. Input more or less has to be a `pandas` dataframe. Only works from a Jupyter notebook. Not very adaptable to custom plot types.\n",
    "\n",
    "```python\n",
    "import altair as alt\n",
    "import pandas as pd\n",
    "data = pd.DataFrame({'x': range(6), 'y': [1,4,3,2,5,8]})\n",
    "alt.Chart(data).mark_line().encode(x='x', y='y').interactive()\n",
    "```\n",
    "\n",
    "**Plotnine:** Also based on the Grammar of Graphics, an attempt at porting `ggplot2` from R. Poor documentation, \"unpythonic\" syntax; input also must be a dataframe.\n",
    "\n",
    "```python\n",
    "import plotnine as pn\n",
    "import pandas as pd\n",
    "data = pd.DataFrame({'x': range(6), 'y': [1,4,3,2,5,8]})\n",
    "pn.ggplot(data) + pn.geoms.geom_line()\n",
    "```\n",
    "\n",
    "**Plotly:** For-profit company, but managed to get a relatively large user base since it was one of the first packages that allowed plotting in the browser. By default, renders plots on their server (which you need an account for, and in some cases a paid subscription).\n",
    "\n",
    "```python\n",
    "import plotly as ply\n",
    "data = ply.graph_objs.Scatter(x=list(range(6)), y=[1,4,3,2,5,8])\n",
    "ply.offline.plot({'data':[data]}, auto_open=True)\n",
    "```\n",
    "\n",
    "**TLDR:** Use `matplotlib`, unless you want to work interactively with large datasets (`bokeh`), or explore complicated statistical relationships within your highly structured data (`altair`).\n",
    "\n",
    "\n",
    "## 3.2 Making publication-quality plots\n",
    "\n",
    "The problem with Matplotlib is that by default the plots it makes are kind of ugly. The simplest way to make them pretty is just to use one of the built-in styles, e.g. adapting our example above:\n",
    "\n",
    "```python\n",
    "import pylab as pl\n",
    "pl.style.use('seaborn')\n",
    "pl.plot([1,4,3,2,5,8])\n",
    "pl.title('Simple example')\n",
    "pl.xlabel('x')\n",
    "pl.ylabel('y')\n",
    "```\n",
    "\n",
    "What makes high-quality plots stand out is typically how much customization is done _for that plot_. Default styles will only get you so far. The following code illustrates some common changes that can make plots look much better: fonts, transparency, spines (whether the plot is in a box or not), legends, and labels. Every aspect of Matplotlib can be customized, so if you can Google it, you can do it!\n",
    "\n",
    "```python\n",
    "\n",
    "import pylab as pl\n",
    "\n",
    "# Simulation parameters\n",
    "ncells = 1000\n",
    "tmax = 500\n",
    "dt = 1\n",
    "period = 30\n",
    "\n",
    "# Calculate spikes\n",
    "pl.seed(1)\n",
    "spikelist = []\n",
    "tvec = pl.arange(0,tmax,dt)\n",
    "firingrate = 5+4*pl.sin(tvec/period)\n",
    "npts = len(tvec)\n",
    "for t in range(npts):\n",
    "    for c in range(ncells):\n",
    "        if pl.rand() < firingrate[t]*dt/1000:\n",
    "            spikelist.append([t,c])\n",
    "spikedata = pl.array(spikelist)\n",
    "\n",
    "# Average firing rate\n",
    "averagerate = pl.zeros(npts)\n",
    "for spike in spikedata:\n",
    "    averagerate[spike[0]] += 1.0/ncells*1000\n",
    "\n",
    "# Trend line\n",
    "trendline = averagerate[:]\n",
    "for repeat in range(5):\n",
    "    trendline = pl.convolve(trendline, [0.25, 0.5, 0.25], 'same')\n",
    "\n",
    "# Set up figure and options\n",
    "pl.rc('font', family='serif', size=14)\n",
    "pl.figure(figsize=(16,12))\n",
    "\n",
    "# Plot the raster\n",
    "pl.axes([0.1,0.35,0.8,0.6])\n",
    "pl.scatter(spikedata[:,0], spikedata[:,1], s=100, alpha=0.3)\n",
    "pl.xlim((0,tmax))\n",
    "pl.gca().spines['top'].set_visible(False)\n",
    "pl.gca().spines['right'].set_visible(False)\n",
    "pl.gca().spines['bottom'].set_visible(False)\n",
    "pl.xticks([])\n",
    "pl.ylabel('Cell ID')\n",
    "\n",
    "# Plot the average firing rate\n",
    "pl.axes([0.1,0.05,0.8,0.25])\n",
    "pl.plot(tvec, averagerate, color=(0.7,0.7,0.7), label='Instantaneous')\n",
    "pl.plot(tvec, trendline, color=(0,0.4,1.0), lw=3, label='Averaged')\n",
    "pl.plot(tvec, firingrate, '--', color=(0,0,0), lw=2, label='Generator')\n",
    "pl.xlim((0,tmax))\n",
    "pl.ylim((0,pl.ylim()[1]))\n",
    "pl.gca().spines['top'].set_visible(False)\n",
    "pl.gca().spines['right'].set_visible(False)\n",
    "pl.xlabel('Time (ms)')\n",
    "pl.ylabel('Average firing rate (Hz)')\n",
    "pl.legend(loc=(0.85,0.93))\n",
    "\n",
    "# Add panel labels\n",
    "pl.axes([0,0,1,1])\n",
    "pl.gca().axis('off')\n",
    "pl.text(x=0.02, y=0.95, s='A', fontsize=35, fontfamily='sans-serif')\n",
    "pl.text(x=0.02, y=0.3, s='B', fontsize=35, fontfamily='sans-serif')\n",
    "\n",
    "# Save the figure\n",
    "pl.savefig('ocnc-example.png', dpi=150)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Exercises\n",
    "\n",
    "\n",
    "## Exercise 4.1: Izhikevich model\n",
    "\n",
    "Convert the following Matlab code into Python. Run both and confirm that they produce the same result.\n",
    "\n",
    "```\n",
    "% Created by Eugene M. Izhikevich, February 25, 2003\n",
    "% Excitatory neurons    Inhibitory neurons\n",
    "\n",
    "% Hint: You will need to import the libraries here (suggestion: import pylab as pl)\n",
    "\n",
    "Ne=800;                 Ni=200;\n",
    "re=rand(Ne,1);          ri=rand(Ni,1);\n",
    "a=[0.02*ones(Ne,1);     0.02+0.08*ri]; % Hint: Use pl.hstack([vector1, vector2])\n",
    "b=[0.2*ones(Ne,1);      0.25-0.05*ri];\n",
    "c=[-65+15*re.^2;        -65*ones(Ni,1)];\n",
    "d=[8-6*re.^2;           2*ones(Ni,1)];\n",
    "S=[0.5*rand(Ne+Ni,Ne),  -rand(Ne+Ni,Ni)];\n",
    "\n",
    "v=-65*ones(Ne+Ni,1);    % Initial values of v\n",
    "u=b.*v;                 % Initial values of u\n",
    "firings=[];             % spike timings % Hint: use pl.empty((2,0))\n",
    "\n",
    "for t=1:1000            % simulation of 1000 ms\n",
    "  I=[5*randn(Ne,1);2*randn(Ni,1)]; % thalamic input\n",
    "  fired=find(v>=30);    % indices of spikes\n",
    "  % Hint: you might need an if statement here\n",
    "  firings=[firings; t+0*fired,fired];\n",
    "  v(fired)=c(fired);\n",
    "  u(fired)=u(fired)+d(fired);\n",
    "  I=I+sum(S(:,fired),2); % Hint: you can specify the axis for sum with .sim(axis=1)\n",
    "  v=v+0.5*(0.04*v.^2+5*v+140-u+I); % step 0.5 ms\n",
    "  v=v+0.5*(0.04*v.^2+5*v+140-u+I); % for numerical\n",
    "  u=u+a.*(b.*v-u);                 % stability\n",
    "end;\n",
    "plot(firings(:,1),firings(:,2),'.'); % Hint: use pl.scatter()\n",
    "% Hint: if the plot doesn't appear, use pl.show() if necessary\n",
    "```\n",
    "\n",
    "## Exercise 4.2: Ant simulation\n",
    "\n",
    "Replace the missing lines with code to run a simulation of ants walking around randomly.\n",
    "\n",
    "```python\n",
    "import pylab as pl\n",
    "\n",
    "class Antsim:\n",
    "    ''' Simulation of ants randomly walking '''    \n",
    "    \n",
    "    def makeants(self, numants=50):\n",
    "        ''' Initialize the ants '''\n",
    "        self.numants = numants # Initialize the number of ants\n",
    "        # Set up the x coordinate: an array of zeros of length numants\n",
    "        # Set up the y coordinate the same way\n",
    "    \n",
    "    def plotants(self, timesteps=150, stepsize=0.03):\n",
    "        ''' Plot the ants '''\n",
    "        pl.figure()\n",
    "        for t in range(timesteps):\n",
    "            pl.clf()\n",
    "            # Increment the x position of each ant by an amount stepsize*randn()\n",
    "            # Increment the y position of each ant by an amount stepsize*randn()\n",
    "            # Scatter x vs y\n",
    "            # Set the x limits to (-1,1)\n",
    "            # Set the y limits to (-1,1)\n",
    "            # Set the plot title to show the current step and the total number of timesteps\n",
    "            pl.pause(1e-3)\n",
    "            \n",
    "\n",
    "# Run the simulation\n",
    "sim = Antsim()\n",
    "sim.makeants(numants=100)\n",
    "sim.plotants()\n",
    "```"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}