assignments/assignment_02/syde556_assignment_02_template.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SYDE 556/750 --- Assignment 2\n",
    "**Student ID: 00000000**\n",
    "\n",
    "*Note:* Please include your numerical student ID only, do *not* include your name.\n",
    "\n",
    "*Note:* Refer to the [PDF](https:/astoeckel/syde556-w20/raw/master/assignments/assignment_02/syde556_assignment_02.pdf) for the full instructions (including some hints), this notebook contains abbreviated instructions only. Cells you need to fill out are marked with a \"writing hand\" symbol. Of course, you can add new cells in between the instructions, but please leave the instructions intact to facilitate marking."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import numpy and matplotlib -- you shouldn't need any other libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Fix the numpy random seed for reproducible results\n",
    "np.random.seed(18945)\n",
    "\n",
    "# Some formating options\n",
    "%config InlineBackend.figure_formats = ['svg']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Generating a random input signal\n",
    "\n",
    "## 1.1 Band-limited white noise\n",
    "\n",
    "**a) Time-domain signals.** Plot $x(t)$ for three randomly generated signals with $\\texttt{limit}$ at $5$, $10$, and $20\\,\\mathrm{Hz}$. For each of these, $\\mathtt{T}=1\\,\\mathrm{s}$, $\\mathtt{dt}=1\\,\\mathrm{ms}$ and $\\mathtt{rms}=0.5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**b) Average power spectrum.** Plot the average $|X(\\omega)|$ (the norm of the Fourier coefficients, or “power spectrum”) over $100$ signals generated with $\\mathtt{T}=1\\,\\mathrm{s}$, $\\mathtt{dt}=1\\,\\mathrm{ms}$, $\\mathtt{rms}=0.5$, and $\\mathtt{limit}=10\\,\\mathrm{Hz}$ (of course, each of these 100 signals should have a different `seed`). The plot should have the $x$-axis labeled “$\\omega$ in radians” and the average $|X|$ value for that $\\omega$ on the $y$-axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 Gaussian power spectrum noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a) Time-domain signals.** Plot $x(t)$ for three randomly generated signals with `bandwidth` at $5$, $10$, and $20\\,\\mathrm{Hz}$. For each of these, $\\mathtt{T}=1\\,\\mathrm{s}$, $\\mathtt{dt}=1\\,\\mathrm{ms}$ and $\\mathtt{rms}=0.5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**b) Average power spectrum.** Plot the average $|X(\\omega)|$ (the norm of the Fourier coefficients, or “power spectrum”) over $100$ signals generated with $\\mathtt{T}=1\\,\\mathrm{s}$, $\\mathtt{dt}=1\\,\\mathrm{ms}$, $\\mathtt{rms}=0.5$, and $\\mathtt{bandwidth}=10$ (of course, each of these 100 signals should have a different `seed`). The plot should have the $x$-axis labeled “$\\omega$ in radians” and the average $|X|$ value for that $\\omega$ on the $y$-axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Simulating a spiking neuron\n",
    "\n",
    "**a) Spike plots for constant inputs.** Plot the spike output for a constant input of $x=0$ over $1$ second. Report the number of spikes. Do the same thing for $x=1$. Use a time step of $\\Delta t = 1\\,\\mathrm{ms}$ for the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**b) Discussion.** Does the observed number of spikes in the previous part match the expected number of spikes for $x=0$ and $x=1$? Why or why not? What aspects of the simulation would affect this accuracy?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "✍ \\<YOUR SOLUTION HERE\\>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**c) Spike plots for white noise inputs.** Plot the spike output for $x(t)$ generated using your function from part 1.1. Use $\\mathtt{T}=1\\,\\mathrm{s}$, $\\mathtt{dt}=1\\,\\mathrm{ms}$, $\\mathtt{rms}=0.5$, and $\\mathtt{limit}=30\\,\\mathrm{Hz}$. Overlay on this plot $x(t)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**d) Voltage over time.** Using the same $x(t)$ signal as in part *c)*, plot the neuron's voltage over time for the first $0.2$ seconds, along with the spikes over the same time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**e) 🌟 Bonus question.** How could you improve this simulation (in terms of how closely the model matches actual equation) without significantly increasing the computation time? $0.5$ marks for having a good idea. Up to $1$ mark for actually implementing it and showing that it works."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "✍ \\<YOUR SOLUTION HERE\\>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Simulating two spiking neurons\n",
    "\n",
    "**a) Spike plots for constant inputs.** Plot $x(t)$ and the spiking output for $x(t)=0$ (both neurons should spike at about $40$ spikes per second), as well as (in a separate plot) $x(t)=1$ (one neuron should spike at $\\approx 150$ spikes per second, and the other should not spike at all)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**b) Spike plots for a sinusodial input.** Plot $x(t)$ and the spiking output for $x(t)=\\frac{1}2 \\sin(10 \\pi t)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**c) Spike plot for a white noise signal.** Plot $x(t)$ and the spiking output for a random signal generated with your function for question 1.1 with $\\mathtt{T}=2\\,\\mathrm{s}$, $\\mathtt{dt}=1\\,\\mathrm{ms}$, $\\mathtt{rms}=0.5$, and $\\mathtt{limit}=5\\,\\mathrm{Hz}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Computing an optimal filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a) Document the code.** Fill in comments where there are `# !`-signs in the Python code. Make sure that your comments (where this makes sense) describe the semantics of the code and do not just repeat what is obvious from the code itself. Run the function with what you wrote for part 3 above, so that it uses the spike signal generated in 3c)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_optimal_filter(\n",
    "        # Signal generated from your white noise generator\n",
    "        x,\n",
    "        # Fourier coefficients from your white noise generator\n",
    "        X,\n",
    "        # Spike train from the previous part\n",
    "        spikes,\n",
    "        # Time step size\n",
    "        dt=1e-3\n",
    "    ):\n",
    "\n",
    "    # x and X should (effectively) be 1D-arrays\n",
    "    assert x.ndim == 1 and X.ndim == 1\n",
    "    assert x.shape[0] == X.shape[0]\n",
    "\n",
    "    # !\n",
    "    Nt = x.size\n",
    "\n",
    "    # Make sure that \"spikes\" is a 2 x Nt array\n",
    "    assert spikes.ndim == 2\n",
    "    assert spikes.shape[0] == 2              \n",
    "    assert spikes.shape[1] == Nt\n",
    "\n",
    "    # !\n",
    "    T = Nt * dt\n",
    "\n",
    "    # !\n",
    "    ts = np.arange(Nt) * dt - T / 2.0\n",
    "\n",
    "    # !\n",
    "    fs = np.arange(Nt) / T - Nt / (2.0 * T)\n",
    "\n",
    "    # !\n",
    "    omega = fs * 2.0 * np.pi\n",
    "\n",
    "    # !\n",
    "    r = spikes[0] - spikes[1]\n",
    "\n",
    "    # !\n",
    "    R = np.fft.fftshift(np.fft.fft(r))\n",
    "\n",
    "    # !\n",
    "    sigma_t = 25e-3\n",
    "\n",
    "    # !\n",
    "    W2 = np.exp(-omega**2*sigma_t**2)\n",
    "\n",
    "    # !\n",
    "    W2 = W2 / sum(W2)\n",
    "\n",
    "    # !\n",
    "    CP = X*R.conjugate()\n",
    "\n",
    "    # !\n",
    "    WCP = np.convolve(CP, W2, 'same')\n",
    "\n",
    "    # !\n",
    "    RP = R*R.conjugate()\n",
    "\n",
    "    # !\n",
    "    WRP = np.convolve(RP, W2, 'same')\n",
    "\n",
    "    # !\n",
    "    XP = X*X.conjugate()\n",
    "\n",
    "    # !\n",
    "    WXP = np.convolve(XP, W2, 'same')\n",
    "\n",
    "    # !\n",
    "    H = WCP / WRP\n",
    "\n",
    "    # !\n",
    "    h = np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(H))).real\n",
    "\n",
    "    # !\n",
    "    XHAT = H*R\n",
    "\n",
    "    # !\n",
    "    xhat = np.fft.ifft(np.fft.ifftshift(XHAT)).real\n",
    "\n",
    "    return ts, fs, R, H, h, XHAT, xhat, XP, WXP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**b) Optimal filter.** Plot the time and frequency plots of the optimal filter for the signal you generated in question 3c). Make sure to use appropriate limits for the $x$-axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**c) Decoded signal.** Plot the $x(t)$ signal, the spikes, and the decoded $\\hat x(t)$ value for the signal from 3c)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**d) Power spectra.** Plot the signal $|X(\\omega)|$, spike response $|R(\\omega)|$, and filtered signal $|\\hat X(\\omega)|$ power spectra for the signal from 3c)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**e) Discussion.** How do these spectra relate to the optimal filter?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "✍ \\<YOUR SOLUTION HERE\\>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**f) Filter for different signal bandwidths.** Plot the optmial filter $h(t)$ in the time domain when filtering spike trains for white noise signals with different `limit` values of $2\\,\\mathrm{Hz}$, $10\\,\\mathrm{Hz}$, and $30\\,\\mathrm{Hz}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**g) Discussion.** Describe the effects on the time plot of the optimal filter as `limit` increases. Why does this happen?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "✍ \\<YOUR SOLUTION HERE\\>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using post-synaptic currents as a filter\n",
    "\n",
    "\n",
    "**a) Plotting the filter for different $n$.** Plot the normalized $h(t)$ for $n=0$, $1$, and $2$, with $\\tau=7\\,\\mathrm{ms}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**b) Discussion.** What two things do you expect increasing $n$ will do to $\\hat{x}(t)$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "✍ \\<YOUR SOLUTION HERE\\>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**c) Plotting the filter for different $\\tau$.** Plot the normalized $h(t)$ for $\\tau=2\\,\\mathrm{ms}$, $\\tau=5\\,\\mathrm{ms}$, $\\tau=10\\,\\mathrm{ms}$, $\\tau=20\\,\\mathrm{ms}$ with $n = 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**d) Discussion.** What two things do you expect increasing $\\tau$ will do to $\\hat{x}(t)$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "✍ \\<YOUR SOLUTION HERE\\>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**e) Decoding a spike-train using the post-synaptic current filter.** Decode $\\hat{x}(t)$ from the spikes generated in question 3c) using an $h(t)$ with $n=0$ and $\\tau=7\\,\\mathrm{ms}$. Do this by generating the spikes, filtering them with $h(t)$, and using that as your activity matrix $A$ to compute your decoders. Plot the time and frequency plots for this $h(t)$. Plot the $x(t)$ signal, the spikes, and the decoded $\\hat{x}(t)$ value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**f) Deocding a spike-train representing a low-frequency signal.** Use the same decoder and $h(t)$ as in part e), but generate a new $x(t)$ with $\\mathtt{limit}=2\\,\\mathrm{Hz}$. Plot the $x(t)$ signal, the spikes, and the decoded $\\hat{x}(t)$ value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ✍ <YOUR SOLUTION HERE>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**g) Discussion.** How do the decodings from e) and f) compare? Explain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "✍ \\<YOUR SOLUTION HERE\\>"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}