x - single_notebook_version_2020.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Creative approaches to problem solving in neuroscience using Python SOLUTIONS",
      "provenance": [],
      "collapsed_sections": [
        "4nDVu2jSRzvT",
        "O3rgE1q5qPhu",
        "9yDCHHQbVxmc",
        "qng-l0-A_FI6",
        "KXVGVwY3Y19V",
        "XHswTKdwRu8Q",
        "X1oiD8xuR3r-",
        "yk_PX4DORoAD",
        "b8BKoP-ZhGPB",
        "tnbG0uPNm8-J",
        "VKG9gt_fnMrz",
        "AEMalT1aIRAZ",
        "a1J02DGdYJjJ",
        "RYSJnOLiYhX9",
        "EBLPGMNhgLqr",
        "i_haO-qkjJEa",
        "yssyMI-lGQ4w",
        "8SMrEq7Fi8nw",
        "9FxuTzrHm8gc",
        "Yjq2P8HAnIk-",
        "Gvaswh34tGRg",
        "I7XSnGBGqTuC",
        "SQV6nUl4To0P",
        "h9lYLV7FtaZp",
        "gKHVmt32vqqp",
        "rE41i9p1wyNh",
        "k0mbotKJvpQ3",
        "cZ5OMIIsJwrG",
        "z9IZefnbyMsS",
        "ahdTrDriMYyM",
        "QVgyo7a4wbsQ",
        "za14qhk3HyvE",
        "gCZK0xtAYwd2",
        "OPalZrijY0gz",
        "9_US2e-VEQFF",
        "sncniINvE-sq",
        "Z80HwVN9PcY_",
        "Km1ccE7fKahR",
        "FEoRvdsUtiK6",
        "AAb-KyPxZQ34",
        "uLXiFlW7ZAeh",
        "FW13rGMfr0PQ",
        "WhwrRDnw-PBy",
        "TMTZyPVwKhhS",
        "i5gzV5i0RdKr",
        "JicSGsIJRjyy",
        "EIx1XJlv5OId",
        "y827oGSuRgCw",
        "iDZ2rKjaoknX",
        "l1FfFW7Czr8v",
        "Dmz4UP-7jw81",
        "JDFoHI6Jj2_N",
        "RX-vtmE3Lmv0",
        "BscX0nrLDjX2",
        "8yR9-qLBIPCW",
        "_KPPX3i8K0KZ",
        "CjhhUvzoK3R7",
        "eVC5T5GrlEbm",
        "twbhAnd3toQ3",
        "ukHE8H6mYhtZ",
        "4ypvA18bZqxp",
        "cNsMUWK_z6Yc",
        "DshsdNisGOTv",
        "as77YuNwHO0e",
        "5N-t6_pEHNrQ",
        "DGU7fv0eblZz",
        "kRiiQcWJSjMF",
        "xZH2g33NpnxT",
        "eFr7kXmQ9U2i"
      ],
      "toc_visible": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rUoEJ0qJ8DTf",
        "colab_type": "text"
      },
      "source": [
        "# Creative approaches to problem solving[$^1$](https://qz.com/192071/how-one-college-went-from-10-female-computer-science-majors-to-40/) in neuroscience using Python\n",
        "## Introduction to (Scientific) Progamming\n",
        "### featuring Neurons and Maths\n",
        "\n",
        "by [Christopher Brian Currin](https://chriscurrin.github.io)\n",
        "\n",
        "A *functional approach* to learning code. *Why-based* in contrast to *fact-based*.\n",
        "\n",
        "![Cape Town](https://images.unsplash.com/photo-1563656157432-67560011e209?ixlib=rb-1.2.1&ixid=eyJhcHBfaWQiOjEyMDd9&auto=format&fit=crop&w=900&h=300&q=80)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "22iZQjozpExG",
        "colab_type": "text"
      },
      "source": [
        "**Imbizo 2020** - version 1.0"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8FbCIjIsBNa3",
        "colab_type": "text"
      },
      "source": [
        "> **what's in a name**\n",
        "\\\n",
        "Python is named after the creator's favourite TV show at the time...\n",
        "\\\n",
        "![what's in a name - monty python](https://external-content.duckduckgo.com/iu/?u=https%3A%2F%2Fm.media-amazon.com%2Fimages%2FM%2FMV5BOTFmYTc3ZWEtNTYxNi00OTA4LTk2NjEtNTI2MTJlNzkyMDdlXkEyXkFqcGdeQWpybA%40%40._V1_UX477_CR0%2C0%2C477%2C268_AL_.jpg&f=1&nofb=1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4nDVu2jSRzvT",
        "colab_type": "text"
      },
      "source": [
        "# The Notebook\n",
        "We'll go over the basics of using Python within the Google Colaboratory (\"Colab\") notebook. [Overview of features](https://colab.research.google.com/notebooks/basic_features_overview.ipynb).\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "O3rgE1q5qPhu",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "## Cells \n",
        "A notebook is a list of cells. Cells contain either explanatory text or executable code and its output. Click a cell to select it. \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9yDCHHQbVxmc",
        "colab_type": "text"
      },
      "source": [
        "## Text cells\n",
        "This is a **text cell**. You can **double-click** to edit this cell. Text cells\n",
        "use markdown syntax. To learn more, see our [markdown\n",
        "guide](/notebooks/markdown_guide.ipynb).\n",
        "\n",
        "You can also add math to text cells using [LaTeX](http://www.latex-project.org/)\n",
        "to be rendered by [MathJax](https://www.mathjax.org). Just place the statement\n",
        "within a pair of **\\$** signs. For example `$\\sqrt{3x-1}+(1+x)^2$` becomes\n",
        "$\\sqrt{3x-1}+(1+x)^2.$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qng-l0-A_FI6",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "## Code cells\n",
        "Below is a **code cell**. Once the toolbar button indicates CONNECTED, click in the cell to select it and execute the contents in the following ways:\n",
        "\n",
        "* Click the **Play icon** in the left gutter of the cell;\n",
        "* Type **Cmd/Ctrl+Enter** to run the cell in place;\n",
        "* Type **Shift+Enter** to run the cell and move focus to the next cell (adding one if none exists); or\n",
        "* Type **Alt+Enter** to run the cell and insert a new code cell immediately below it.\n",
        "\n",
        "There are additional options for running some or all cells in the **Runtime** menu.\n",
        "\n",
        "> **NB Lesson**\n",
        "> \n",
        "> ---\n",
        "> learn keyboard shortcuts for things you do often! \n",
        "\n",
        "Use the `Ctrl/Cmd + ]` shortcut to collapse sections if you're feeling the notebook get a bit slow.\n",
        "\n",
        "*it may take a while the first time you run any cell in a notebook while it CONNECTS to a Python runtime*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JPH6-yM6V8X3",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "to_output = \"this is a code cell with output below\"\n",
        "to_output"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KXVGVwY3Y19V",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "### Tasks\n",
        "- Try adding some text to this cell.\n",
        "- \\[**advanced**] Add a *formatted* equation to this cell.\n",
        "- Run the code cell\n",
        "\n",
        "Note: If you want to clear the notebook's workspace (all the currently stored variables) use the menu option Runtime --> Restart Runtime...\n",
        "\n",
        "Okay, let's get going!\n",
        "\n",
        "*remember that along the way we will be encountering some weird and wonderfully new terms + concepts that may not complete sense at first...* **Be patient with a) yourself and a) others** "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XHswTKdwRu8Q",
        "colab_type": "text"
      },
      "source": [
        "# Why are we coding?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "urVNdk-tT-e6"
      },
      "source": [
        "\n",
        "(A) To make lots of money?\n",
        "\n",
        "(B) To boast to our friends?\n",
        "\n",
        "(C) To solve a problem?\n",
        "\n",
        "(D) To run statistics on our results?\n",
        "\n",
        "(E) To visualise our results?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pXPe_yJIR8ep",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title Answer { display-mode: \"form\" }\n",
        "print(\"Coding is merely a tool within our large repertoire of possible approaches to solving a problem.\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X1oiD8xuR3r-",
        "colab_type": "text"
      },
      "source": [
        "#### But before you even start coding\n",
        "\n",
        "1. identify what you *do* know\n",
        "1. identify what you *do* ***not*** know\n",
        "1. form an hypothesis about what you *expect* to find\n",
        "1. have a plan\n",
        "\n",
        "This is called software **engineering**. \n",
        "\n",
        "A **programmer** knows how to code and may have the technical skills needed to solve problems. *This is hard.*\n",
        "\n",
        "A **software engineer** follows a process of understanding requirements, working with others and developing a solution that solves a problem (or answers a question). *This is harder.*\n",
        "\n",
        "But we are **scientists**, so we want to answer the *how* and *why* something occurs, perhaps with no clear end goal. We use a systematic approach to problem solving that often involves creative and innovate solutions. *This is hardest.*\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yk_PX4DORoAD",
        "colab_type": "text"
      },
      "source": [
        "## Coding is *formal* Problem Solving"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "s0p1jzZGRkT9",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "**Step 1: Start with a problem**\n",
        "\n",
        "*Remember that a problem is informed by our current knowledge and limitations thereof*\n",
        "\n",
        "E.g. how do I best visualise my results?\n",
        "\n",
        "![example data](https://miro.medium.com/max/458/1*Y0-OMtzCuC4H5ldl_rC5Yw.png)\n",
        "\n",
        "**Step 2: Formulate an hypothesis** <-- this is often skipped but very N.B.\n",
        "\n",
        "E.g. a pie chart is the best\n",
        "\n",
        "**Step 3: Identify and evaluate solutions**\n",
        "\n",
        "E.g. compare a pie chart to a text table\n",
        "\n",
        "E.g. try a bunch of charts\n",
        "\n",
        "E.g. make a donut pie chart with an inner and outer portion\n",
        "\n",
        "**Step 4: Implement solution**\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YV6t2oT_r2_v",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title **[Run me]** Donut pie chart { display-mode: \"form\" }\n",
        "import matplotlib.pyplot as plt\n",
        " \n",
        "# Data to plot\n",
        "labels = ['Python', 'C++', 'Ruby', 'Java']\n",
        "sizes = [504, 337, 415, 280]\n",
        "labels_gender = ['Man','Woman']*len(labels)\n",
        "sizes_gender = [315,189,125,212,270,145,190,90]\n",
        "colors = ['#ff6666', '#ffcc99', '#99ff99', '#66b3ff']\n",
        "colors_gender = ['#c2c2f0','#ffb3e6']*len(labels)\n",
        " \n",
        "# Plot\n",
        "fig, ax = plt.subplots()\n",
        "ax.pie(sizes, labels=labels, colors=colors, startangle=90, frame=True)\n",
        "ax.pie(sizes_gender, colors=colors_gender, radius=0.75, startangle=90)\n",
        "\n",
        "# create donut hole\n",
        "centre_circle = plt.Circle((0,0),0.5,color='black', fc='white',linewidth=0)\n",
        "ax.add_artist(centre_circle)\n",
        "\n",
        "ax.axis('equal')\n",
        "fig.tight_layout()\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uaDrvQwGr1pz",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "\n",
        "**Step 5: Evaluate results**\n",
        "\n",
        "[the issue with pie chart](https://www.data-to-viz.com/caveat/pie.html)\n",
        "\n",
        "**Step 6: You probably encountered a problem (i.e. step 6 is actually step 1)**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "T0LP_FTQxtQ0",
        "colab_type": "text"
      },
      "source": [
        "# The Neuron\n",
        "![neuron](https://images.unsplash.com/photo-1564325724739-bae0bd08762c?ixlib=rb-1.2.1&ixid=eyJhcHBfaWQiOjEyMDd9&auto=format&fit=crop&w=1366&&h=500&q=80)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "I13G6P4oVzG_",
        "colab_type": "text"
      },
      "source": [
        "**Step 1: Start with a problem**\n",
        "\n",
        "How does inhibition influence a neuron's transfer function (input-output [I-O], I-F)?\n",
        "\n",
        "*what do we know already?*\n",
        "\n",
        "**Step 2: Formulate an hypothesis**\n",
        "\n",
        "Inhibition causes a divisive modulation of the transfer function (decreases the maximum firing rate and changes the I-O slope).\n",
        "\n",
        "Alternative hypotheses?\n",
        "\n",
        "**Step 3: Identify and evaluate solutions**\n",
        "\n",
        "some options\n",
        "1. record a neuron while it is receiving inhibition.\n",
        "  \n",
        "    where we know:\n",
        "    1. the amount of *inhibition* being received at each receptor\n",
        "    1. the amount of *excitation* being received at each receptor \n",
        "    1. the density of *passive* channels in the dendrites\n",
        "    1. the density of *active* channels in the dendrites\n",
        "    1. the concentration of ions that change the strength of channels\n",
        "    1. ... more\n",
        "1. simulate a neuron where we *define* all$^1$ of these\n",
        "\n",
        "**Step 4: Implement solution**\n",
        "\n",
        "We will be doing (part of) this now\n",
        "\n",
        "**Step 5: Evaluate results**\n",
        "\n",
        "- Plot voltage trace over time.\n",
        "- Plot firing rate as a function of (excitatory) input frequency\n",
        "- Plot firing rate as a function of (excitatory + inhibitory) input frequency\n",
        "\n",
        "**Step 6: You probably encountered a problem (i.e. step 6 is actually step 1)**\n",
        "\n",
        "What about stats?\n",
        "\n",
        "---\n",
        "\n",
        "Remember that **a problem is informed by our current knowledge and limitations thereof**\n",
        "\n",
        "\n",
        "$1$: most$^2$\n",
        "\n",
        "$2$: some\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "K7r7fyJxuK1m",
        "colab_type": "text"
      },
      "source": [
        "**Let's start with what we know**\n",
        "\n",
        "![alt text](https://images.unsplash.com/photo-1483706571191-85c0c76b1947?ixlib=rb-1.2.1&ixid=eyJhcHBfaWQiOjEyMDd9&auto=format&fit=crop&w=1366&h=500&q=80)\n",
        "\n",
        "> note: we will be implementing this from first principles instead of using a library like [NEURON](https://neuron.yale.edu/neuron/), [Brian2](https://brian2.readthedocs.io/en/stable/), or the many other available.\n",
        "> This will help teach us some of the **fundamentals of programming**. Plus we get to appeciate some of the math (*yes I know, scary right!*🙀) and its implementation.\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "b8BKoP-ZhGPB",
        "colab_type": "text"
      },
      "source": [
        "## Neurons are electrical ⚡ \n",
        "![alt text](https://images.unsplash.com/photo-1507413245164-6160d8298b31?ixlib=rb-1.2.1&ixid=eyJhcHBfaWQiOjEyMDd9&auto=format&fit=crop&w=1366&h=500&q=80)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fblsVmY0vnV7",
        "colab_type": "text"
      },
      "source": [
        "**Membranes are capacitors.**\n",
        "\n",
        "> Q: What is the equation for a capacitor?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yURmHTrivlIf",
        "colab_type": "text"
      },
      "source": [
        "![capacitance](https://upload.wikimedia.org/wikipedia/commons/0/0c/Cable_theory_Neuron_RC_circuit_v3.svg)\n",
        "\n",
        "The membrane potential ($V$ or $V_m$ or even sometimes $u$) of a neuron **changes** in time according to all the electrical input it receives (in the form of currents, $I$) and its *membrane capacitance* $C_m$.\n",
        "\n",
        "In the absence of any input, the potential is at its resting value $V_{rest} = -65 mV$. What is this set by?\n",
        "\n",
        "If the neuron receives synaptic input from other neurons, the *time-dependent*, $t$, potential $V(t)$ will be deflected from its resting value.\n",
        "\n",
        "$C_m \\frac{dV}{dt} = \\sum_{k}{I_k(t)} \\tag{1}$\n",
        "\n",
        "where $k$ are the sources of current, such as synaptic input.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tnbG0uPNm8-J",
        "colab_type": "text"
      },
      "source": [
        "### Neurons are lazy\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4Gzvf0nm-ejX",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "However, the neuron will want to return to it's resting membrane potential. Because everyone likes to rest.\n",
        "\n",
        "$C_m \\frac{dV}{dt} = -\\frac{(V(t) - V_{rest})}{R_m} + \\sum_{k}{I_k(t)} \\tag{2}$\n",
        "\n",
        "where $R_m$ is the *membrane resistance*. This new term $-\\frac{(V(t) - V_{rest})}{R_m}$ represents the **leak current** and is often formulated as \n",
        "$I_{leak}(t) = \\frac{(V(t) - V_{rest})}{R_m} = \\bar g_{leak}(V(t) - \\bar E_{leak}) \\tag{3}$\n",
        "\n",
        "Where $g$ is the **conductance** and $E$ is the **reversal potential**. The bar symbol $\\bar g$ mean it is a constant. \n",
        "\n",
        "> Much of the complexity and richness of neuronal dynamics arises because membrane conductances change over time. However, some of the factors that contribute to the total membrane current can be treated as relatively constant (i.e. $\\bar g$ instead of $g(t)$), and these are typically grouped together into a single term called the **leakage current**. The currents carried by ion pumps that maintain the concentration gradients that make equilibrium potentials nonzero typically fall into this category.\n",
        "\n",
        "Together, $C_m$ and $R_m$ therefore set the time it takes for input current to shift the voltage away from rest (aka the rise time as the capacitor charges) and the time it takes to return to rest (aka the decay time as the capacitance dessipates through the leak). This is known as the **membrane time constant** $\\tau_m = C_m \\cdot R_m$. \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VKG9gt_fnMrz",
        "colab_type": "text"
      },
      "source": [
        "### Experimentalists suck"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "U2L4kok0-fs3",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "Of course, experimentalists love poking neurons with electrodes, so we also need to consider externally injected current $I_e$. \n",
        "\n",
        "$C_m \\frac{dV}{dt} = -I_{leak}(t) - \\sum_{k}{I_k(t)} + \\frac{I_e(t)}{A}\n",
        "\\tag{4}$\n",
        "\n",
        "where $A$ is the surface area of the injected current (sometimes $A$ is ignored when dealing with neurons with no size - $\\rm{point}$ $\\rm{neurons}$)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FsXuTvOsS0kc",
        "colab_type": "text"
      },
      "source": [
        "**How do we even go about putting this into a computer!?**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4wNwmtEXS5L5",
        "colab_type": "text"
      },
      "source": [
        "What do we need to know in order to figure out what we don't?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "mJC8NnOVUHmp",
        "colab": {}
      },
      "source": [
        "# Parameters for neuronal voltage\n",
        "V_rest = -65  # resting membrane potential (mV)\n",
        "C_m = 1       # membrane capacitance (nF)\n",
        "R_m = 100     # membrane resistance (MOhm)\n",
        "I_e = 0       # external current (nA)\n",
        "A = 10        # surface area of electrode (um^2)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8mgYr9c5Vvfo",
        "colab_type": "text"
      },
      "source": [
        "Our first block of code!\n",
        "\n",
        "Did you run the cell?\n",
        "\n",
        "what does it do?\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AEMalT1aIRAZ",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "## Assignment"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kIDvh6aGYUtO",
        "colab_type": "text"
      },
      "source": [
        "In general terms, these are examples of **variable assignment** - assigning a value to a name (because we can change the value later)\n",
        "\n",
        "As in many programming languages, a Python assignment statement associates a symbolic *name* on the left-hand side with a *value *on the right-hand side. In Python, we say that names refer to values, or a name is a reference to a value. This is done using the equals '`=`' operator. (In R this can be '`<-`' and in Go this can be '`:=`').\n",
        "\n",
        "> *if you're very confuzzled by variables and assigning values to names, check out this [post](https://nedbatchelder.com/text/names.html)*"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "a1J02DGdYJjJ",
        "colab_type": "text"
      },
      "source": [
        "## Comments\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-RG4j4UbIHAy",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "**Comments** in Python are preceded by a '`#`' and indicate extra information *that is magically ignored when you run the cell*. Comments *in other languages* can be '`// single line comment`' or '\n",
        "```\n",
        "/* multi\n",
        "line\n",
        "comment */```\n",
        "\n",
        "We use comments because **most** of the time we write code in lots of files that end in '`.py`' (aka a python file) instead of in a notebook (which allows nicely formatted text cells).\n",
        "\n",
        "Comments are necessary for writing good code. \n",
        "They are there for a few reasons.\n",
        "1. To describe what the code below does (best to describe blocks of code at a time).\n",
        "1. To provide documentation that can be quickly previewed when developing.\n",
        "1. To provide context for why values or algorithms were chosen, especially for\n",
        "obscure cases where someone may try an 'obvious' optimisation that you've tried and doesn't work. Save human time. When variable names are short symbols (c) because that is what your equations use, a comment next to them are vital. Also, if you're not using a library with units, include the units in comments.\n",
        "\n",
        "    Example:\n",
        "    ```\n",
        "    # Speed of light (in a vacuum)\n",
        "    # This value was found on https://en.wikipedia.org/wiki/Speed_of_light\n",
        "    c = 299792458 # m/s. \n",
        "    mass = 100000 # g\n",
        "    E = mass*c**2 # python uses ** instead of ^ for powers\n",
        "    ```\n",
        "\n",
        "**Who** are they for?\n",
        "1. Future you. \n",
        "    \n",
        "    That person in 6+ months that needs to run or edit the code because\n",
        "    someone asked you to add *just one* more feature.\n",
        "1. Collaborators\n",
        "    \n",
        "    Those who work on the same code as you, or guide what you're coding and \n",
        "    need to understand what your code does and why. Plus, it builds confidence\n",
        "    in your ability to write reproducible code.\n",
        "1. Scientists building on your work\n",
        "    \n",
        "    Open science is important for generating knowledge. Sharing code that only \n",
        "    a handful of people know how it works is counter to the ethos of open \n",
        "    science."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RYSJnOLiYhX9",
        "colab_type": "text"
      },
      "source": [
        "## Showing the result"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "u3JjcAnUYoxe",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# we 'print out' the value of a variable to the *console*\n",
        "# Python has the 'print' *method* that does this for us. There is more than one \n",
        "# way to *call* the method, which are mostly equivalent.\n",
        "print('V_rest =', V_rest) # 'default'\n",
        "print('C_m = %d' % C_m) # very old school\n",
        "print('R_m = {}'.format(R_m)) # old\n",
        "print(f'I_e = {I_e}') # newest and best (note the f before the ')\n",
        "print(f\"The surface area of the neuron is {A} mm^2\")\n",
        "\n",
        "display(\"display is notebook specific\")\n",
        "\"the last line in a code cell is displayed\""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gyrMuYgNZhKs",
        "colab_type": "text"
      },
      "source": [
        "Notice the **comments** are ignored in the output.\n",
        "\n",
        "The 'stuff' in between the apostrophes `'` or `\"` are **strings** of characters.\n",
        "\n",
        "Strings are just another value we can assign to a name.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cOcoep2Kau_g",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "nrn_name = 'my first neuron'\n",
        "print(nrn_name)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2euL4PxvE1KS",
        "colab_type": "text"
      },
      "source": [
        "a ***formatted* string** starts with `f\"` (or `f'`) and ends with `\"` (or `'`). \n",
        "\n",
        "Formatted strings easily allow us to create a line of text along with values. Values can be added to the string by wrapping them in curly brackets `{` `}`. \n",
        "\n",
        "*use these by default* but be aware you may encounter the other forms shown above."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZjSa8MWvF4Nf",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "n = 10\n",
        "nth_neuron_lit = \"my {n}th neuron\"  # note the value doesn't start with an 'f'\n",
        "print(nth_neuron_lit)\n",
        "nth_neuron = f\"my {n}th neuron\"     # we assign the formatted string to a name\n",
        "print(nth_neuron)                   # print out the variable\n",
        "ten_more = f\"my {n+10}th neuron\"    # we can do some basic operations too\n",
        "print(ten_more)                   "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LKCJGa7qbIB9",
        "colab_type": "text"
      },
      "source": [
        "**Task**: How do you display\n",
        "\n",
        "`Chris says, \"I have just create my 100th neuron!\"`?\n",
        "\n",
        "*challenge*:\n",
        "\n",
        "`Chris says, \"I've just create my 100th neuron!\"`? (hint: characters can be *escaped* by preceding it with a `\\`)\n",
        "\n",
        "`Chris says, \"I've just create my {100}th neuron!\"`? (hint: use triple braces `{{{nrn_num}}}` instead of escaping with `\\`)\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cxtXN6c7Hv_k",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title TODO: display text above { run: \"auto\" }\n",
        "nrn_num =  99#@param {type:\"integer\"}\n",
        "\n",
        "# implement your solution here\n",
        "# ...\n",
        "print('Chris says, \"I\\'ve just create my 100th neuron!\"')\n",
        "print(\"Chris says, \\\"I've just create my 100th neuron!\\\"\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "U9wVlGfabUbX",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title solution { display-mode: \"form\" }\n",
        "print(f\"Chris says, \\\"I've just create my {nrn_num}th neuron!\\\"\")\n",
        "print(f'Chris says, \"I have just create my {nrn_num}th neuron!\"')\n",
        "print(f'Chris says, \"I have just create my {{{nrn_num}}}th neuron!\"')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "bRflH_QlgG63"
      },
      "source": [
        "if you're thinking, \"wait, this is _way_ too easy\"\n",
        "\n",
        "1. display `γ=1.20 t/m³` without unicode characters (how many ways can solve this?)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "MZucuNDBgG6x"
      },
      "source": [
        "still too easy?\n",
        "\n",
        "2. write a preprocessor to remove all comments including a multiline comment that starts with `#/` and ends with `/#`\n",
        "\n",
        "   so this\n",
        "   ```Python\n",
        "   a = 3 # single line comment\n",
        "   # comment for line below\n",
        "   # b = 2 # this line is commented out from the beginning\n",
        "   ## double hashes are valid comments\n",
        "   #/ start of multiline comment\n",
        "   print(\"this should not display\")\n",
        "   /#\n",
        "   c = \"👋\"      #/ start of multiline comment\n",
        "                   print(\"this should not display\")\n",
        "   d,e = 'd',\"e\" /#\n",
        "   ```\n",
        "   becomes\n",
        "   ```Python clean\n",
        "   a = 3\n",
        "   c = \"👋\"\n",
        "   ```\n",
        "   (note: you must raise an error if there is anything after `/#` on the same line) "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EBLPGMNhgLqr",
        "colab_type": "text"
      },
      "source": [
        "## From text to code\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "paY6po8E-qX1",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "Reminder\n",
        "\n",
        "$C_m \\frac{dV}{dt} = -I_{leak}(t) - \\sum_{k}{I_k(t)} + \\frac{I_e(t)}{A}\n",
        "\\tag{4}$\n",
        "\n",
        "$I_{leak}(t) = \\frac{(V(t) - V_{rest})}{R_m} \\tag{3}$\n",
        "\n",
        "\n",
        "We can't write equation 4 like this in python\n",
        "```\n",
        "C_m * dV/dt = -I_leak + sum_I_k + I_e/A\n",
        "```\n",
        "because '=' is an **assignment operator** so the left-hand side needs to be a **name**.\n",
        "\n",
        "Note that we **can** write this to *check if assignment was done correctly*\n",
        "```\n",
        "C_m * dV/dt == -I_leak + sum_I_k + I_e/A\n",
        "```\n",
        "\n",
        "To fix this, we need to \n",
        "- re-arrange the equation to have a single term on the left-hand side\n",
        "- define all symbols as variables (names with values)\n",
        "\n",
        "let's ignore $I_k$ and $I_e$ for now\n",
        "\n",
        "___\n",
        "*Ultimatey, we want to implement how voltage changes by updating and keeping track of the voltage*\n",
        "\n",
        "\n",
        "$$\n",
        "\\begin{aligned}\n",
        "  V(t+1) &= V(t) + dV_1 \\\\\n",
        " V(t+2) &= V(t+1) + dV_2 \\\\\n",
        " V(t+3) &= V(t+2) + dV_3 \\\\\n",
        " \\cdots\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "We therefore need to calculate $dV$ at each time step and update $V$ accordingly.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i_haO-qkjJEa",
        "colab_type": "text"
      },
      "source": [
        "#### 1. re-arrange equation 4\n",
        "From (4), \n",
        "\n",
        "$\\begin{align}\n",
        "C_m\\frac{dV}{dt} &= -I_{leak}(t) - \\sum_{k}{I_k(t)} + \\frac{I_e(t)}{A} \\\\\n",
        "\\frac{dV}{dt} &= (-I_{leak}(t) - \\sum_{k}{I_k(t)} + \\frac{I_e(t)}{A})\\cdot \\frac{1}{C_m} \\\\\n",
        "dV &= (-I_{leak}(t) - \\sum_{k}{I_k(t)} + \\frac{I_e(t)}{A}) \\cdot \\frac{dt}{C_m}\n",
        "\\end{align}$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yssyMI-lGQ4w",
        "colab_type": "text"
      },
      "source": [
        "#### 2. Define all symbols as variables\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vk053GZALVb9",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "For equations above:\n",
        "\n",
        "Math symbol | Variable name | Defined | Comment\n",
        "---    | ---        | --- | ---\n",
        "$V(t)$  | `V_t`     | N | Membrane potential (mV)\n",
        "$dV$  | `dV`      | N | Change in membrane potential (mV)\n",
        "$I_{leak}$ | `I_leak` | N | Leakage current (nA)\n",
        "$I_e$  | `I_e` | N/A | External current (nA)\n",
        "$A$    | `A`   | Y | Surface area of electrode (mm$^2$)\n",
        "$I_k$  | `I_syn` | N/A | Synaptic current (nA)\n",
        "$dt$  | `dt`      | N | Time step (ms)\n",
        "$C_m$  | `C_m`      | Y | Membrane capacitance (nF)\n",
        "\n",
        "For equation 3: \n",
        "$I_{leak}(t) = \\frac{(V(t) - V_{rest})}{R_m} = \\bar g_{leak}(V(t) - \\bar E_{leak}) \\tag{3}$\n",
        "\n",
        "Math symbol | Variable name   | Defined | Comment\n",
        "---    | ---         | --- | ---\n",
        "$I_{leak}$ | `I_leak` | N | Leakage current (nA)\n",
        "$V(t)$  | `V_t`      | N | Membrane potential (mV)\n",
        "$R_m$  | `R_m`        | Y | Membrane resistance (MΩ) \n",
        "$V_{rest}$  | `V_rest`| Y | Resting membrane potential (mV)\n",
        "$\\bar g_{leak}$  | `g_leak` | N | Leakage conductance (μS)\n",
        "$\\bar E_{leak}$  | `E_leak` | N | Leakage reversal potential (mV)\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BatV0q8Xe7wq",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "V_t = V_rest + 2 # this is the initial value of V(t)\n",
        "I_leak = (V_t - V_rest)/R_m # according to equation 3 (we divide using /)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HT10pHdecJmG",
        "colab_type": "text"
      },
      "source": [
        "We could implement `I_leak` like this too"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vWOQ_9ZLcXhN",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "E_leak = V_rest # these are the same thing and now have the same value! (this only works if V_rest exists)\n",
        "g_leak = 1/R_m # we divide using /\n",
        "I_leak_alt = g_leak * (V_t - E_leak) # according to equation 3 (we multiply using *)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ABwcGNs7d8cR",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# we can check if things are equal to each other using ==\n",
        "print(f\"are V_rest and E_leak the same? Answer: {V_rest == E_leak}\")\n",
        "# remember this last line is outputted to the console automatically\n",
        "g_leak*R_m == 1"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xqelV7k4dJnp",
        "colab_type": "text"
      },
      "source": [
        "What if you only have a single '`=`'?\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1k4Hk4-qR59U",
        "colab_type": "text"
      },
      "source": [
        "Equation 3 is now complete\n",
        "\n",
        "Math symbol | Variable name   | Defined\n",
        "---    | ---         | ---\n",
        "$I_{leak}$ | `I_leak` | Y\n",
        "$V(t)$  | `V_t`      | Y\n",
        "$R_m$  | `R_m`        | Y\n",
        "$V_{rest}$  | `V_rest`| Y\n",
        "$\\bar g_{leak}$  | `g_leak` | Y\n",
        "$\\bar E_{leak}$  | `E_leak` | Y\n",
        "\n",
        "Now for equation 4:\n",
        "\n",
        "Math symbol | Variable name   | Defined\n",
        "---    | ---        | ---\n",
        "$V(t)$  | `V_t`      | Y\n",
        "$dV$  | `dV`      | N\n",
        "$I_{leak}$ | `I_leak` | Y\n",
        "$I_e$  | `I_e` | N/A\n",
        "$A$    | `A`   | Y\n",
        "$I_k$  | `I_k` | N/A\n",
        "$dt$  | `dt`      | N\n",
        "$C_m$  | `C_m`      | Y\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8SMrEq7Fi8nw",
        "colab_type": "text"
      },
      "source": [
        "#### 3. Discretisation\n",
        "i.e. making some discrete (in time) instead of continous\n",
        "\n",
        "*because life is full of holes*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "M7iTco2bjGI8",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# we want to perform a simulation *over time*\n",
        "dt = 0.01 # ms --> note that this has decimals\n",
        "# we could also write it this way where after the 'e' is the exponent part\n",
        "0.01 == 1e-2 == 10**-2 == pow(10,-2) # note that python uses ** for 'to the power of'\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-rlTUZ8CjITF",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "dV = (-I_leak + I_e) * dt/C_m\n",
        "print(f\"dV = {dV} mV\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LZdUMQr0ljbL",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title What about different dts? { run: \"auto\", vertical-output: true, display-mode: \"form\" }\n",
        "dt = 0.445 #@param {type:\"slider\", min:0, max:1, step:0.005}\n",
        "dV = (-I_leak + I_e) * dt/C_m\n",
        "print(f\"dV = {dV:.4f} mV\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CmOb3YfpmbeP",
        "colab_type": "text"
      },
      "source": [
        "## Ok, we have some problems"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9FxuTzrHm8gc",
        "colab_type": "text"
      },
      "source": [
        "#### 1. We have not updated `V_t`. Solution : **re-assignment**\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "J6XVOgf2nBgY",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "print(f\"V_t before is {V_t:.4f}\")\n",
        "V_t = V_t + dV # overwrite the old value with the new value\n",
        "print(f\"V_t after is {V_t:.4f}\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Yjq2P8HAnIk-",
        "colab_type": "text"
      },
      "source": [
        "#### 2. we only do this once. Solution: **loops**\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3Z76aHWvngMR",
        "colab_type": "text"
      },
      "source": [
        "If you want to do the same thing over and over again we use **loops**.\n",
        "\n",
        "The most common loop is the **for loop** because it starts with the special name **`for`** (you cannot assign anything to `for`. Go ahead, try it).\n",
        "\n",
        "The general syntax (way to write valid code) is...\n",
        "> note that the `<stuff between>` the less than `<` and greater than `>` symbols are placeholders\n",
        "\n",
        "```\n",
        "for <temp var name> in range(<number of iterations>):\n",
        "    <statement on new indented line using tabs or spaces>\n",
        "```\n",
        "\n",
        "This is rather long so we use shorthand var names such as `i`, `j`, `k` etc. \n",
        "\n",
        "```\n",
        "N = 10\n",
        "for i in range(N):\n",
        "    print(i)\n",
        "```\n",
        "\n",
        "Note: if you have previously defined `i` (you have assigned a value to it), it will be overwritten and lost forever (or until you run the program again)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vUpp_BSbBuTg",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "i = 100\n",
        "print(f\"value of i BEFORE for loop: i={i}\")\n",
        "N = 10\n",
        "for i in range(N):\n",
        "    print(i)\n",
        "print(f\"value of i AFTER for loop: i={i}\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7k0ZUgRenMO1",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "print(f\"V_t at the start is {V_t}\")\n",
        "for i in range(10):\n",
        "  #print(f\"iteration number: {i}\") # you may want to comment this out if you're doing a bunch of interations\n",
        "  V_t = V_t + dV\n",
        "print(f\"V_t at the end is {V_t}\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jb0M8tyjpJNC",
        "colab_type": "text"
      },
      "source": [
        "Notice that `i` starts at 0!\n",
        "\n",
        "This is because loops are often used for retrieving an index of a list.\n",
        "\n",
        "99% of languages start indexing from 0, which language(s) don't?\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oFkg0Ir-syfO",
        "colab_type": "text"
      },
      "source": [
        "**Q1: How many times do you have to run the cell until you reach V_rest?** \n",
        "\n",
        "**Q2: What happens once you reach V_rest?**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "83hM3nh0MXRk",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "`range` actually has a number of ways to use it.\n",
        "\n",
        "1. `range(stop)`\n",
        "1. `range(start, stop)`\n",
        "1. `range(start, stop, step)` \n",
        "\n",
        "where start is *inclusive* and stop is *exclusive*"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Gvaswh34tGRg",
        "colab_type": "text"
      },
      "source": [
        "##### How did I know this?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pZ5wrGtCtLK3",
        "colab_type": "text"
      },
      "source": [
        "Someone told me!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "I7XSnGBGqTuC",
        "colab_type": "text"
      },
      "source": [
        "###### What if no-one is around?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dyO8R8Ugqp74",
        "colab_type": "text"
      },
      "source": [
        "1. Read the [documentation](https://docs.python.org/3/)\n",
        "1. search [Google](https://www.google.com), (or [Bing](https://www.bing.com) or [DuckDuckgo](https://duckduckgo.com/) if that's your thing) and [StackOverflow](https://stackoverflow.com/)\n",
        "1. use an IDE (integrated development environment, like this one or [PyCharms](https://www.jetbrains.com/pycharm/)) that performs *code inspection*.\n",
        "1. play around with it"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SQV6nUl4To0P",
        "colab_type": "text"
      },
      "source": [
        "###### **Tasks**\n",
        "- try different values in `range`.\n",
        "- what if `start` is bigger than `stop`?\n",
        "- what if values are negative?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BPn6rnv7NBDa",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# play around with for loops and range here\n",
        "for i in range(1, 10):\n",
        "  print(i)\n",
        "print() # prints a blank line\n",
        "for i in range(1, 10, 2):\n",
        "  print(i)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "h9lYLV7FtaZp",
        "colab_type": "text"
      },
      "source": [
        "#### 3. `V_t` depends on `I_leak`, which depends on `V_t`. Solution: **blocks of code** for *dependencies*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "feu-mB51tnxe",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title I_leak is time-dependent { run: \"auto\", vertical-output: true }\n",
        "num_iter = 4060 #@param {type:\"slider\", min:10, max:10000, step:10}\n",
        "dt = 0.5 #@param {type:\"slider\", min:0, max:1, step:0.005}\n",
        "\n",
        "# because we only define I_leak once, but it depends on V_t, \n",
        "# we need to reassign it in the loop too\n",
        "\n",
        "# TODO: re-define V_t\n",
        "V_t = V_rest + 2\n",
        "print(f\"V_t at the start is {V_t:.4f} mV\")\n",
        "for i in range(num_iter):\n",
        "  # TODO: equation 3\n",
        "  # I_leak = ...\n",
        "  # TODO: calculate dV\n",
        "  # dV = ...\n",
        "  V_t = V_t + dV\n",
        "print(f\"V_t at the end is {V_t:.4f} mV\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gKHVmt32vqqp",
        "colab_type": "text"
      },
      "source": [
        "##### Q: How many iterations for a dt of 0.5 gets us to $V_{rest}$ (within an error margin of 0.0001)?\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HiduBm_jwfiK",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title Solution { run: \"auto\", vertical-output: true, display-mode: \"form\" }\n",
        "decimal_places = 8 #@param {type:\"slider\", min:0, max:8, step:1}\n",
        "num_iter = 10000 # choose max\n",
        "dt = 0.5\n",
        "\n",
        "# because we only define I_leak once, but it depends on V_t, \n",
        "# we need to reassign it in the loop too\n",
        "\n",
        "V_t = V_rest + 2 # let's start V_t again\n",
        "print(f\"V_t at the start is {V_t:.4f} mV\")\n",
        "for i in range(num_iter):\n",
        "  I_leak = (V_t - V_rest)/R_m\n",
        "  dV = (-I_leak + I_e) * dt/C_m\n",
        "  V_t = V_t + dV\n",
        "  if round(V_t, decimal_places) <= V_rest:\n",
        "    # stop the loop early\n",
        "    break\n",
        "print(f\"V_t at the end is {V_t:.4f} mV\")\n",
        "print(f\"the numer of iterations it took is {i+1}\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pMnn9-rrQT3z",
        "colab_type": "text"
      },
      "source": [
        "**Challenge** implement the above as a `while` loop"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PKN2G9s4QS_P",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# while ...:\n",
        "  # do stuff"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EDomwRF9ocNM",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "\n",
        "Good work, we now have a program that simulates neuronal voltage over time\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MjaCvqMWoFb6",
        "colab_type": "text"
      },
      "source": [
        "I've plotted an example below, but we will go through how to plot in Python a bit later"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5DH17807l3Xs",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title Voltage over time with injected current{ run: \"auto\", vertical-output: true, display-mode: \"form\" }\n",
        "\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "num_iter = 5000\n",
        "dt = 0.5 \n",
        "I_e = 0.05 #@param {type:\"slider\", min:0, max:0.1, step:0.01}\n",
        "start = 500 #@param {type:\"slider\", min:0, max:2500, step:10}\n",
        "end = 1000 #@param {type:\"slider\", min:0, max:2500, step:10}\n",
        "I_i1 = start/dt\n",
        "I_i2 = end/dt\n",
        "V_t = V_rest\n",
        "_v = []\n",
        "for i in range(num_iter):\n",
        "  _v.append(V_t)\n",
        "  I_leak = (V_t - V_rest)/R_m\n",
        "  _i_e = I_e  if I_i1 < i < I_i2 else 0\n",
        "  dV = (-I_leak + _i_e) * dt/C_m\n",
        "  V_t = V_t + dV\n",
        "\n",
        "plt.plot(np.arange(0,num_iter*dt,dt),_v, c='k')\n",
        "plt.xlabel(\"Time (ms)\")\n",
        "plt.ylabel(\"Voltage (ms)\")\n",
        "y = plt.ylim()[1]\n",
        "\n",
        "plt.annotate(\"$I_e$\",xy=(end,y),xytext=(start,y), \n",
        "             ha='right', va='center', fontsize=14, c='b',\n",
        "             arrowprops=dict(arrowstyle='-', ec='b', linewidth=5))\n",
        "plt.annotate(\"$V_m$\",xy=(num_iter*dt,_v[-1]), \n",
        "             ha='left', va='center', fontsize=14)\n",
        "ax = plt.gca()\n",
        "for spine in [ax.spines['top'],ax.spines['right']]:\n",
        "  spine.set_visible(False)\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rE41i9p1wyNh",
        "colab_type": "text"
      },
      "source": [
        "## Connecting the pieces: Synapses\n",
        "![lights](https://images.unsplash.com/photo-1556691421-cf15fe27a0b6?ixlib=rb-1.2.1&ixid=eyJhcHBfaWQiOjEyMDd9&auto=format&fit=crop&w=1366&&h=500&q=80)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pP2v_syUupPs",
        "colab_type": "text"
      },
      "source": [
        "Let's provide some *synaptic input* to our neuron.\n",
        "\n",
        "i.e. Let's implement $\\sum_{k}{I_k(t)}$\n",
        "\n",
        "$I_{syn}(t) = g_{syn}(t) (V(t)−E_{syn}) \\tag{5}$\n",
        "\n",
        "Note that $g_{syn}$ has 'lost' it's flat cap $\\bar{g}$ and gained a dependence on time $(t)$\n",
        "\n",
        "*not a problem*\n",
        "\n",
        "A simple choice for the time course of the synaptic conductance in is an exponential decay\n",
        "\n",
        "$g_{syn}(t) = \\sum_{s} \\bar g_{max} e^{−(t−t_{s})/τ} \\cdot Θ(t−t_{s}) \\tag{6}$\n",
        "\n",
        "which can be viewed as conductance given a time point and all the incoming spike times. We can split this up into \n",
        "1. conductance given a time point and a presynaptic spike time\n",
        "1. The sum of 1. for each spike time, which gives us the conductance given a time point\n",
        "\n",
        "$\\begin{aligned}\n",
        "&g_{syn}(t, t_{s}) = \\bar g_{max} e^{−(t−t_{s})/τ} \\cdot Θ(t−t_{s}) \\\\\n",
        "g_{syn}(t) = \\sum_{s}^{S} &g_{syn}(t, t_{s})\n",
        "\\end{aligned}$\n",
        "\n",
        "**don't freak out if this doesn't make sense**\n",
        "\n",
        "- $\\tau$ is the **decay [time constant](https://en.wikipedia.org/wiki/Time_constant)** of the synapse (e.g. 5 ms).\n",
        "  > after 1 $\\tau$, a value will have decreased to $1/e \\approx 36.8 %$ of its initial value.\n",
        "\n",
        "  ![time constant](https://upload.wikimedia.org/wikipedia/commons/thumb/3/39/Series_RC_resistor_voltage.svg/230px-Series_RC_resistor_voltage.svg.png)\n",
        "- $t_{s}$ denotes the arrival time of a presynaptic action potential - $s$ is each presynaptic *spike*.\n",
        "- Θ(x) or $H(x)$ is the [Heaviside step function](https://en.wikipedia.org/wiki/Heaviside_step_function).\n",
        "\n",
        "  ![Heaviside step function](http://mathworld.wolfram.com/images/eps-gif/HeavisideStepFunction_900.gif)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZCXRyW380ted",
        "colab_type": "text"
      },
      "source": [
        "Let's implement the Heaviside step function and introduce Python abilities, *viz.* **functions** (also known as **methods** but there is a small difference) and **conditionals**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "x7hGNZGQuEje",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def heaviside(x):\n",
        "  \"\"\" a python method starts with 'def' (definition) followed by the method name \n",
        "  'heaviside', the *parameters* (here just 'x') for the method are in '(' ')'\n",
        "   followed by a ':' then an **indented** newline where the *method body* starts.\n",
        "  \n",
        "  note this comment just below the method definition is in quotations and can be\n",
        "   multi-line\n",
        "  \"\"\"\n",
        "  if x < 0:\n",
        "    print(0)\n",
        "  elif x > 0: # we could also write this as 0 < x because we are doing a *comparison*\n",
        "    print(1)\n",
        "  else: # only thing left is x being 0. could also write \"elif x==0\":\n",
        "    print(1/2)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gyXOcldd0PxS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# test our method\n",
        "heaviside(-1000)  # should display 0\n",
        "heaviside(0)      # should display 0.5\n",
        "heaviside(0.1)    # should display 1"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "k0mbotKJvpQ3",
        "colab_type": "text"
      },
      "source": [
        "### Functions\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "P7aVd4RNuaF-",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "If there is a computation during your program's flow that you will be doing a bunch, functions help *abstract* the logic into more manageable bits. This *abstraction* reduces complexity and increases efficiency. \n",
        "\n",
        "\n",
        "So we want to change our logic we implemented from this\n",
        "\n",
        "![Half-max convention](https://wikimedia.org/api/rest_v1/media/math/render/svg/34f164d5bf42583f4f09a2871a3f589ff0a89d43)\n",
        "\n",
        "to this\n",
        "\n",
        "![regular](https://wikimedia.org/api/rest_v1/media/math/render/svg/f1783c84465f7a602fae566c34efa63f48c84212)\n",
        "\n",
        "\n",
        "Imagine we used the logic of the Heaviside function in a bunch of places. We'd have to find and replace each of these. But because we used a **function**, we can change our definition in one place.\n",
        "\n",
        "While we are changing things, it's not a good idea to have `print` statements inside methods because they don't perform any useful logic. Instead, we want the **method** to provide us with a **value** back, which we can then display *if we want*. \n",
        "E.g. \n",
        "\n",
        "```\n",
        "print(heaviside(1e9))\n",
        "# or\n",
        "h_offset = heaviside(1e-2) + 2 # offset of 2\n",
        "print(h_offset)\n",
        "```\n",
        "\n",
        "Try this **before** running the next code block **and after**.\n",
        "\n",
        "---\n",
        "\n",
        "If you want to use any names or values from inside the function, you need to **`return`** these names or a value directly (*a method does **not** return anything, functions do*)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6Na015JMyLYr",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def heaviside(x):\n",
        "  \"\"\"Heaviside step function which returns 1 when x is 0\"\"\"\n",
        "  if x < 0:\n",
        "    return 0\n",
        "  elif x >= 0: # or could use else:\n",
        "    return 1\n",
        "\n",
        "# test the method\n",
        "x = 0.1\n",
        "y=2\n",
        "answer = heaviside(y)\n",
        "print(f\"H({y}) = {answer}\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PWm7z851j7ai",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def heaviside(x):\n",
        "  \"\"\"Heaviside step function which returns True when x is 0 or greater, and \n",
        "  False otherwise\"\"\"\n",
        "  return x >= 0     # does this make sense to you? True == 1 and False == 0\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cZ5OMIIsJwrG",
        "colab_type": "text"
      },
      "source": [
        "#### *parameters* and *arguments*\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zVXU7Lvhuc8H",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "When you define a method, you can have a number of variables in its definition (aka *declaration*) which can be used in the **method body**. These variables in the definition are known as **parameters** (`x` was a parameter in our `heaviside` declaration). \n",
        "\n",
        "These **parameters** only exist within the method. Outside the method, no-one cares about them. The concept of variables being local to where they are defined (e.g. in a method) is known as **scope**.\n",
        "\n",
        "You may be asking the question then, why do we have `x` outside the function/method then?\n",
        "\n",
        "Well, this `x` is an **argument** we *passed* to the function so it can use its value. The name of the argument could be anything (or you can pass a value directly).\n",
        "\n",
        "For our use case they had the same name which helps with *readability*, but this may not always be the case.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0Oy4p6x6jr6L",
        "colab_type": "text"
      },
      "source": [
        "**meta tips**\n",
        "\n",
        "hovering over a method reveals its signature - it's parameters and description (the text in `\"\"\" \"\"\"` at the top of a method.)\n",
        "\n",
        "`dir(<object>)` returns a list of methods which the object can call\n",
        "\n",
        "`help(<object>)` similarly works and combines the above 2"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xPFik7xaCqQS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "help(heaviside)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "z9IZefnbyMsS",
        "colab_type": "text"
      },
      "source": [
        "### Conditional expressions\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "M7kLiXblubcO",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "If you want to branches in your logic, you can use Python's `if`, `elif`, `else` keywords with **comparison expressions** (`<expr>`).\n",
        "\n",
        "The syntax is as follows:\n",
        "```\n",
        "if <expr1>:\n",
        "  <statement if expr1 is True>\n",
        "  <other statements>\n",
        "  <as many as you like>\n",
        "elif <expr2>:\n",
        "  <statement if expr1 is False but expr2 is True>\n",
        "elif <expr3>:\n",
        "  <statement if expr1 is False AND expr2 is False but expr3 is True>\n",
        "...\n",
        "else:\n",
        "  <statement if nothing in the logic branch is True> \n",
        "```\n",
        "\n",
        "Things to keep in mind:\n",
        "1. You must **start** with `if`. A new `if` is a new logic tree (not a branch).\n",
        "1. You can have **many** `elif`s. `elif`is short for `else if` and means that the previous expressions have been been satisfied.\n",
        "1. You can end in `else`, but there must be only one if you do.\n",
        "1. `else if` won't work. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ahdTrDriMYyM",
        "colab_type": "text"
      },
      "source": [
        "### Implementing $g_{syn}$\n",
        "\n",
        "We are going to combine a bit of what we learnt to write the $g_{syn}$ equation (6).\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3Eo9hIOCNhAe",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "e = 2.71828  # Euler's constant\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Dwm8iEdQDpmx",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "\n",
        "def g_syn(t, spk_t, g_max=50, tau=5):\n",
        "  \"\"\"Single exponential conductance decay using equation 6 \n",
        "  but limited to a single incoming signal (spk_t).\n",
        "  \"\"\"\n",
        "  # note that g_max is in nS\n",
        "  return g_max * e**(-(t-spk_t)/tau) * heaviside(t-spk_t)\n",
        "\n",
        "# test our function\n",
        "t = 15\n",
        "print(g_syn(t, 14, 100, 6))   # call with all arguments passed\n",
        "print(g_syn(t, 14, 100))      # argument 'tau' omitted, the default value is used\n",
        "print(g_syn(t, 14, tau=6))    # argument 'g_max' omitted and 'tau' explicitly provided\n",
        "print(g_syn(t, 14))           # both g_max and tau defaults used\n",
        "print(g_syn(t, spk_t=14, g_max=100, tau=6))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tqDzgXZ-STPe",
        "colab_type": "text"
      },
      "source": [
        "The main issue with this definition is that it only works for a single incoming signal (`t_f`).\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QVgyo7a4wbsQ",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "### Types\n",
        "\n",
        "Up to now we have been dealing with\n",
        "- **`int`egers** (numbers with no decimal point) like `-1`, `0`, `2`, or `1000`\n",
        "- **`float`**ing points numbers (aka real numbers with decimal point) like `e = 2.71828`\n",
        "- **`str`**ing of characters like `'this'` or `\"also this\"` or ```\"\"\"this which can span multiple lines\"\"\"```\n",
        "\n",
        "- [**`bool`ean**](https://techterms.com/definition/boolean) values which are either `True` or `False`. These arise from **comparisons** such as in **conditionals**. \n",
        "  \n",
        "  E.g. \n",
        "  \n",
        "  ```\n",
        "  a=1\n",
        "  b=2\n",
        "\n",
        "  if a>b:\n",
        "    print(\"a is bigger than b\")\n",
        "  else:\n",
        "    print(\"a is not bigger than b)\n",
        "\n",
        "  # the below is the same\n",
        "  c = a>b   # assign our comparison to a variable\n",
        "  print(c)  # displays 'False' in the console\n",
        "  if c:\n",
        "    print(\"a is bigger than b\")\n",
        "  else:\n",
        "    print(\"a is not bigger than b)\n",
        "  ```"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "za14qhk3HyvE",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "#### Aside: Boolean Logic\n",
        "Aka multiple comparisons.\n",
        "\n",
        "So far, we have been doing single comparisons at a time.\n",
        "\n",
        "If we want to perform multiple comparisons we use **`and`** or **`or`**. If we want to **negate** some logic we use **not**. Logic can be grouped with parentheses `(` `)`.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3g2RtYmqHn_m",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title Try display each of the 8 possibilities { run: \"auto\", vertical-output: true, display-mode: \"both\" }\n",
        "is_horse_like = True #@param [\"True\", \"False\"] {type:\"raw\"}\n",
        "has_horn = True #@param [\"True\", \"False\"] {type:\"raw\"}\n",
        "has_stripes = True #@param [\"True\", \"False\"] {type:\"raw\"}\n",
        "\n",
        "# with 3 booleans there are 2^3 = 8 possible combinations\n",
        "#   2 is the base because it is either True or False\n",
        "\n",
        "if is_horse_like:\n",
        "  print(\"is_horse_like, so either 🐴 or 🦓 or 🦄 or striped unicorn. It's a...\")\n",
        "  # 'nest' some logic\n",
        "  if has_horn and not has_stripes:\n",
        "    print(\"🦄\")\n",
        "  if not (has_horn or has_stripes): # could also be `elif not has_stripes` --> why?\n",
        "    print(\"🐴 1\")\n",
        "  elif has_stripes and not has_horn:\n",
        "    print(\"🦓\")\n",
        "  else:\n",
        "    print(\"striped unicorn!\")\n",
        "elif has_horn:\n",
        "  print(\"not horse_like, has a horn, but could have stripes...\")\n",
        "  if not has_stripes:\n",
        "    print(\"🦏\")\n",
        "  if has_stripes: # better to just use 'else'\n",
        "    print(\"a tatooed narwhal\")\n",
        "elif has_stripes:\n",
        "  print(\"neither horselike, nor has_horn, but has stripes\")\n",
        "  print(\"🐯\")\n",
        "else:\n",
        "  print(\"could be anything else, even a 🦑\")\n",
        "\n",
        "if is_horse_like and not has_horn and not has_stripes:\n",
        "  print(\"🐴 1\") \n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "g-M0HDrfHIxH",
        "colab_type": "text"
      },
      "source": [
        "it is useful to read the Python docs on [Types and Comparisons](https://docs.python.org/3.8/library/stdtypes.html)\n",
        "\n",
        "but these are the 2 important bits (haha) right now\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gCZK0xtAYwd2",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "##### Boolean Operations — `and`, `or`, `not`\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "L6WsnWDGYt4m",
        "colab_type": "text"
      },
      "source": [
        "These are the Boolean operations, ordered by ascending priority:\n",
        "\n",
        "Operation | Result | Notes\n",
        "---    | --- | ---\n",
        "`x or y` | if x is false, then y, else x | (1)\n",
        "`x and y`| if x is false, then x, else y | (2)\n",
        "`not x` | if x is false, then True, else False | (3)\n",
        "\n",
        "Notes:\n",
        "1. This is a short-circuit operator, so it only evaluates the second argument if the first one is false.\n",
        "1. This is a short-circuit operator, so it only evaluates the second argument if the first one is true.\n",
        "1. `not` has a lower priority than non-Boolean operators, so `not a == b` is interpreted as `not (a == b)`, and `a == not b` is a syntax error.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OPalZrijY0gz",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "##### Comparisons\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ue-Wf7RgY2a8",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "There are eight comparison operations in Python. They all have the same priority (which is higher than that of the Boolean operations). Comparisons can be chained arbitrarily; for example, `x < y <= z` is equivalent to `x < y and y <= z`, except that `y` is evaluated only once (but in both cases `z` is not evaluated at all when `x < y` is found to be false).\n",
        "\n",
        "This table summarizes the comparison operations:\n",
        "\n",
        "Operation | Meaning\n",
        "---|---\n",
        "`<`|strictly less than\n",
        "`<=`|less than or equal\n",
        "`>` | strictly greater than\n",
        "`>=` | greater than or equal\n",
        "`==` | equal\n",
        "`!=` | not equal\n",
        "`is` | object identity\n",
        "`is not` | negated object identity\n",
        "\n",
        "Objects of different types, except different numeric types, never compare equal. The `==` operator is always defined but for some object types (for example, class objects) is equivalent to `is`. The `<`, `<=`, `>` and `>=` operators are only defined where they make sense; for example, they raise a `TypeError` exception when one of the arguments is a complex number."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cjrtrzsWFBB5",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# play around with Booleans and conditionals here\n",
        "n = 20\n",
        "n_less = n - 10\n",
        "n_more = n + 10\n",
        "\n",
        "n_less < n < n_more"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Vu_RmvFdFxKj",
        "colab_type": "text"
      },
      "source": [
        "for collections such as `str` strings we also have access to **`in`**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sCBvpfoZFyin",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "s = \"string of characters\"\n",
        "small = \"small\"\n",
        "big = \"big\"\n",
        "\n",
        "print(small < big)\n",
        "print(small is not big)\n",
        "print(small != big)\n",
        "print('c' in 'character')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9_US2e-VEQFF",
        "colab_type": "text"
      },
      "source": [
        "##### [*Advanced*] Bitwise operators"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mSFGFKzxEe6B",
        "colab_type": "text"
      },
      "source": [
        "many languages use `&&` for `and` and `||` for `or`.\n",
        "\n",
        "Python is different$^1$. However, it *also* has `&` and `|` which does **bitwise** comparison. \n",
        "\n",
        "Bits are 1s and 0s and can represented using `0b` at the start or by converting a number using `bin(<number>)`"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6gD5hKxREkW_",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "s = 0b1001\n",
        "c = 67\n",
        "print(f\"\\t...8421\")\n",
        "print(f\"-\"*20)\n",
        "print(f\"\\t{s:0>7b} ({s:3g})\\n\\t{c:0>7b} ({c:3g})\")\n",
        "print(f\"s|c\\t{s|c:0>7b} ({s|c:3g}) bitwise or\")\n",
        "print(f\"s&c\\t{s&c:0>7b} ({s&c:3g}) bitwise and\")\n",
        "print(f\"s^c\\t{s^c:0>7b} ({s^c:3g}) bitwise xor (exclusive or)\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6xkLfwP3sPBv",
        "colab_type": "text"
      },
      "source": [
        "These operators are useful when dowing element-wise comparisons in arrays (see below)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6YoKo2H2s-3o",
        "colab_type": "text"
      },
      "source": [
        "> Q: Implement exclusive or using logical (instead of bitwise) operator(s)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qEfWWEiIIpvu",
        "colab_type": "text"
      },
      "source": [
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sncniINvE-sq",
        "colab_type": "text"
      },
      "source": [
        "#### Lists\n",
        "\n",
        "A **List** is an ordered$^1$ collection of values (or names which point to values).\n",
        "\n",
        "They are created with hard brackets `[ ]` or `list( )`.\n",
        "  \n",
        "  E.g. `[1,2,3]` or `[3, 2, 'hi', 4.1]` or  even\n",
        "  \n",
        "  `[\"a\", \"list\", [\"inside\", \"a\"], \"list\"]`\n",
        "\n",
        "Because lists are ordered, we can retrieve a *list element* by **indexing** with integers. \n",
        "\n",
        "$^1$ lists are ordered by *insertion order* - the order in which elements are put in the list\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XSM1touAJm7J",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "spks = [1.9, 22, 40] # arriving action potentials [spike times] (ms)\n",
        "print(f\"spks looks like this: {spks}\")\n",
        "print(f\"the first element (index 0) is {spks[0]}\")\n",
        "print(f\"the second element (index 1) is {spks[1]}\")\n",
        "print(f\"the third element (index 2) is {spks[2]}\")\n",
        "print(f\"-1 is a shortcut for the last element!: spks[-1] == {spks[-1]}\")\n",
        "two_spks = [10, 20]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Z80HwVN9PcY_",
        "colab_type": "text"
      },
      "source": [
        "##### Arrays\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SMV0zl-oR8vI",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "A similar concept to lists is **arrays**. \n",
        "\n",
        "*Lists* are useful for keeping a collection of values that stay the same.\n",
        "\n",
        "*Arrays* are useful for performing computations on a collection of values, e.g. matrices.\n",
        "\n",
        "*Arrays* must have the same data type (all numbers or all strings). *Lists* can combine types.\n",
        "\n",
        "We create arrays in Python using a `library` known as *numerical python* (`numpy`).\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3Ef1rC7rL7vP",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "arr = np.array([1.9, 22, 40])\n",
        "print(f\"spks looks like this: {arr}\")\n",
        "print(f\"the first element (index 0) is {arr[0]}\")\n",
        "print(f\"the second element (index 1) is {arr[1]}\")\n",
        "print(f\"the third element (index 2) is {arr[2]}\")\n",
        "print(f\"-1 is a shortcut for the last element!: arr[-1] == {arr[-1]}\")\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DNL5QY_GJ4nR",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "help(np.array)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Km1ccE7fKahR",
        "colab_type": "text"
      },
      "source": [
        "For a 500 ms simulation (spike times can be between 0 ms and 500 ms)\n",
        "\n",
        "#### Tasks\n",
        "1. Create a list of 100 spike times - `one_hundy_spks`\n",
        "1. Create a list of 100 spike times from a random process - `rand_spks`\n",
        "1. Sort `rand_spks` by time"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9p3ljCsOmuTk",
        "colab_type": "text"
      },
      "source": [
        "**Task 1**:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nAPInwP5myKF",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import random\n",
        "one_hundy_spks = [] # start with an empty list\n",
        "rand_spks = []\n",
        "for i in range(100):\n",
        "  rand_spks.append(random.random()*400 + 100) # between 100 and 500 ms\n",
        "sorted(rand_spks)\n",
        "\n",
        "rand_spks = np.random.uniform(low=0, high=500, size=100)\n",
        "# implement your solution here\n",
        "# hint: we can add a value to the end of a list using one_hundy_spks.append(<value>)\n",
        "# TODO: put a number in the list\n",
        "# TODO: 100 times..."
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "USALH5lsndfs",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title run this to test your code\n",
        "import numpy as np\n",
        "try:\n",
        "  assert len(one_hundy_spks) == 100, \"you don't have 100 spike times\"\n",
        "  assert np.all(np.array(one_hundy_spks) <= 500), \"there's a spike time more than 500\"\n",
        "  assert np.all(np.array(one_hundy_spks) >= 0), \"there's a spike time less than zero\"\n",
        "except AssertionError as ae:\n",
        "  print(f\"failed because {ae}\")\n",
        "else:\n",
        "  print(\"passed\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZbBoqBD4nLn3",
        "colab_type": "text"
      },
      "source": [
        "Solution:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sU7Y3lpynK7v",
        "colab_type": "code",
        "cellView": "both",
        "colab": {}
      },
      "source": [
        "#@title Solution\n",
        "one_hundy_spks = []\n",
        "for i in range(100):\n",
        "  one_hundy_spks.append(i*5)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ue1Vy3Roowmp",
        "colab_type": "text"
      },
      "source": [
        "**Task 2**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IL0lNFruKehs",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# option 1: using the random library\n",
        "import random # <-- we have now gained access to extra functionality!"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nozwuhIWxRCQ",
        "colab_type": "text"
      },
      "source": [
        "Imagine if someone (or a bunch of people) had implemented useful functions already and made them available to us...\n",
        "\n",
        "Enter **libraries** or *packaged bits of functionality*.\n",
        "\n",
        "We have just **imported** the `random` library. \n",
        "\n",
        "Read the [documentation](https://docs.python.org/3/library/random.html) to find out what methods it has!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fn8GTdqPxi6j",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "rand_spks = [] # start with an empty list\n",
        "# implement your solution here\n",
        "# hint: we can use the randint function in the random library\n",
        "# TODO: get a random number\n",
        "# TODO: put the number in the list\n",
        "# TODO: 100 times..."
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZpZI78IwpQM_",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title run this to test your code\n",
        "import numpy as np\n",
        "try:\n",
        "  assert len(rand_spks) == 100, \"you don't have 100 spike times\"\n",
        "  _test_arr = np.array(rand_spks)\n",
        "  assert np.all(_test_arr <= 500), \"there's a spike time more than 500\"\n",
        "  assert np.all(_test_arr >= 0), \"there's a spike time less than zero\"\n",
        "  assert sorted(rand_spks) != rand_spks, \"you didn't use a random process (sneaky sneaky)\"\n",
        "except AssertionError as ae:\n",
        "  print(f\"failed because {ae}\")\n",
        "else:\n",
        "  print(\"passed\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "neEvtpZ4o4SQ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title Solution (integers) { display-mode: \"form\" }\n",
        "rand_spks=[]\n",
        "for i in range(100):\n",
        "  r = random.randint(0, 500) # the randint function in the random library\n",
        "  rand_spks.append(r)\n",
        "  \n",
        "print(f\"number of items in `rand_spks` is {len(rand_spks)}\")\n",
        "print(f\"here are the first ten: {rand_spks[0:10]}\")\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AwueEuykrAn2",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title Solution (floating point numbers) { display-mode: \"form\" }\n",
        "rand_spks=[]\n",
        "for i in range(100):\n",
        "  # the random function in the random library returns a value between 0 and 1\n",
        "  r = random.random()\n",
        "  r_scale = r * 500 + 0 # upper bound is 500 and lower bound is 0\n",
        "  rounded = round(r_scale, 2)\n",
        "  rand_spks.append(rounded)\n",
        "  \n",
        "print(f\"number of items in `rand_spks` is {len(rand_spks)}\")\n",
        "print(f\"here are the first ten: {rand_spks[0:10]}\")\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "unkvcWOVsQyT",
        "colab_type": "text"
      },
      "source": [
        "**Task 3**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "f3iqPwsLsS3b",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# implement your solution here\n",
        "# TODO: write some code ..."
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dQw3-GKEroWG",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title Solution (using *list comprehension*) { display-mode: \"form\" }\n",
        "rand_spks = sorted([round(random.random()*500+0, 2) for i in range(100)])\n",
        "print(f\"here are the first ten: {rand_spks[0:10]}\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Im4maxi5reNi",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title Solution (using arrays) { display-mode: \"form\" }\n",
        "import numpy as np\n",
        "rand_spks = np.sort(np.random.rand(100) * 500)\n",
        "rand_spks"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FEoRvdsUtiK6",
        "colab_type": "text"
      },
      "source": [
        "### Tying it all together\n",
        "Creating `g_syn` that takes a list of spikes"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6l_qfz7Etlat",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from math import e # instead of defining e manually, we use the math library which has a more accurate value\n",
        "\n",
        "def g_syn(t, spks, g_max=50, tau=5):\n",
        "  \"\"\"Single exponential conductance decay using equation 6.\n",
        "  see https://colab.research.google.com/drive/1xFXWT3ScBHAJgaszYKpUzhcFh6jBwxoi#scrollTo=pP2v_syUupPs\n",
        "\n",
        "  :param t: time point at which to evaluate (ms)\n",
        "  :param spks: list of time points (ms) of incoming signals (presynaptic action potentials)\n",
        "  :param g_max: maximum synaptic conductance (nS) [default: 50]\n",
        "  :param tau: synapse time constant (ms) [default: 5]\n",
        "  :return: value of synaptic conductance (nS) at time t\n",
        "\n",
        "  Tests\n",
        "  ------\n",
        "  >>> g_syn(1,[1],50,5) # command starting with '>>>' and expected output below\n",
        "  50.0\n",
        "  >>> g_syn(1,1,50,5) # a number instead of list is passed\n",
        "  50.0\n",
        "  >>> g_syn(1,2,50,5) # a number instead of list is passed\n",
        "  0.0\n",
        "  >>> g_syn(1,np.array([1]),50,5) # a numpy array\n",
        "  50.0\n",
        "  >>> g_syn(1,[10000],50,5) # a spk is waaaay ahead of t\n",
        "  0.0\n",
        "  >>> g_syn(1,[1,10000],50,5) # a spk is waaaay ahead of t\n",
        "  50.0\n",
        "  >>> g_syn(10000,[1],50,5) # t is waaaay after spk\n",
        "  0.0\n",
        "  \"\"\"\n",
        "  total = 0.0  # we iterate from 0 until the length of the spks list (e.g. 100)\n",
        "  for i in range(len(spks)):\n",
        "    total += g_max * e**(-(t-spks[i])/tau) * heaviside(t-spks[i])\n",
        "  return total\n",
        "\n",
        "# test our function\n",
        "t = 40                # time point we want to investigate\n",
        "spks = [1.9, 22, 40]  # a list of spike times (ms)\n",
        "print(f\"the first spike comes at t = {spks[0]} ms\") # lists start indexing at 0!\n",
        "print(f\"the last spike comes at t = {spks[-1]} ms\") # -1 is shortcut for the last element"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XLC4Wcj1DgLD",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "t = 40\n",
        "spks = np.array([1.9, 22, 40])\n",
        "spks+10"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AAb-KyPxZQ34",
        "colab_type": "text"
      },
      "source": [
        "#### **Testing**\n",
        "run basic tests in docstrings (method comments)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sW2BZVboztQ9",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import doctest\n",
        "# test the entire module (in this case every cell that has been run)\n",
        "#doctext.testmod() # this becomes burdensome if we just want to test g_syn\n",
        "doctest.run_docstring_examples(g_syn, globals(), name=g_syn.__name__, verbose=False)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uLXiFlW7ZAeh",
        "colab_type": "text"
      },
      "source": [
        "#### **Tasks**\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WuQ3nOrgupiO",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "You will have (hopefully) noticed some doctsring tests failed. Let's fix that 🛠\n",
        "\n",
        "of course we can change the tests themselves, but *shortcuts **now** lead to broken results **later**.*\n",
        "\n",
        "**Task 1**\n",
        "\n",
        "Fix the errors by examining what error was \"raised\" or \"thrown\" (terminology for when an error occurs) and changing the `g_syn` method appropriately.\n",
        "\n",
        "Hints:\n",
        "1. you can check a variable's type using `type(<var name>) is int` or `isinstance(<var name>, int)`. `int` is an example of a type you can check against. The others include `float`, `str`, `list`, `dict`, `bool`, and `np.array`. numpy also has the `np.iterable` method\n",
        "1. you can **cast** to a type (convert from one to another) by using the type name and parenthesis. E.g. `int(3.2) == 3`\n",
        "1. `list` and `[]` are **not** synonymous. Try this:\n",
        "```\n",
        "d = {'a':1,'b':2} \n",
        "list(d) == [d] # is this True or False?\n",
        "```\n",
        "  \n",
        "**Task 2**\n",
        "\n",
        "Make the method more efficient by using the numpy library. \n",
        "\n",
        "Hints:\n",
        "1. \"Vectorise\" the `for` loop (i.e. use a `np.array` which enables linear algebra operations)\n",
        "1. `np.exp(...)` is shorthand for `e**(...)` or `power(e,...)`(and more efficient too!)\n",
        "1. you can create **boolean masks** to only select certain elements of an array\n",
        "```\n",
        "a = np.array([1,2,3,4,5,6,7,8,9])\n",
        "mask = a>5\n",
        "print(a[mask]) # outputs [6 7 8 9]\n",
        "```\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "z13jtxbLfyLX",
        "colab": {}
      },
      "source": [
        "#@title Solution { display-mode: \"form\" }\n",
        "def g_syn_solution(t, spks, g_max=50, tau=5):\n",
        "  if not np.iterable(spks) or type(spks) is list:\n",
        "      spks = np.array(spks)\n",
        "  return np.sum(g_max * np.exp(-(t-spks[t>=spks])/tau))\n",
        "g_syn_solution.__doc__ = g_syn.__doc__.replace(\"g_syn\",\"g_syn_solution\")\n",
        "doctest.run_docstring_examples(g_syn_solution, globals(), \n",
        "                               name=g_syn_solution.__name__, verbose=True)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4RypqroqS0Fb",
        "colab_type": "text"
      },
      "source": [
        "## Neuron with synapses"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OPEngNMkvmgk",
        "colab_type": "text"
      },
      "source": [
        "<img src=\"https://4pyffg.am.files.1drv.com/y4m3icQdXi7J2szFDJRQvhu7bgff9sXokUqmKwgHxqs0G88zqGf82V_CTHD5pAff3KgVbVp7GI_thNPL5fqkM_AQ8cXFH9pbNu-fWZue5QnJyuTuimYhUDEoOtN-vUO7NIRwVx12VEQeCLsvBZEMUUXcWKf-bUiORFfkduEI2UvPswTEgpbcrUZSV0AmOgKuNQnOadJ7lYyA8KI3R0-7NY61w?width=1494&height=1079&cropmode=none\" width=\"700\" />"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "F2mcIJDSPDj2",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Excitatory synapse\n",
        "def I_AMPA(t, V_t, spks, g_max=50, tau=4, E=0):\n",
        "  # note we convert from nS and mV to mA \n",
        "  # --> motivate why 1e-3 instead of having the reader guess\n",
        "  return 1e-3 * g_syn(t, spks, g_max, tau) * (V_t-E)\n",
        "\n",
        "# Inhibitory synapse\n",
        "def I_GABA(t, V_t, spks, g_max=50, tau=8, E=-70):\n",
        "  return 1e-3 * g_syn(t, spks, g_max, tau) * (V_t-E)\n",
        "\n",
        "T = 100       # duration\n",
        "t = 10        # time point\n",
        "V_t = -64.2   # V_t at t\n",
        "# spike times\n",
        "E_spks = [1.9, 22, 40]                    # defined\n",
        "I_spks = np.random.randint(0, t, size=2)  # random\n",
        "print(f\"I_AMPA{t, V_t, E_spks} = {I_AMPA(t, V_t, E_spks):>8.4f} mA\")\n",
        "print(f\"I_GABA{t, V_t, I_spks} = {I_GABA(t, V_t, I_spks):>8.4f} mA\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "33fT0mJ9SMHS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "V_t = V_rest + 2 # let's start V_t again\n",
        "print(f\"V_t at the start is {V_t:.2f} mV\")\n",
        "spks = [1.9] # a list of spike times (ms)\n",
        "num_iter = 1000\n",
        "dt = 0.1\n",
        "for i in range(num_iter):\n",
        "  t = i*dt\n",
        "  I_leak = (V_t - V_rest)/R_m\n",
        "  I_syn = I_GABA(t, V_t, I_spks) + I_AMPA(t, V_t, E_spks)\n",
        "  dV = (-I_leak - I_syn + I_e) * dt/C_m\n",
        "  V_t = V_t + dV\n",
        "print(f\"V_t after {dt*num_iter} ms is {V_t:.2f} mV\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Zz17JYMCTRdk",
        "colab_type": "text"
      },
      "source": [
        "___"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4qlKKIdjfjPU",
        "colab_type": "text"
      },
      "source": [
        "## Simulation"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Wjt1ql5YVgts",
        "colab_type": "text"
      },
      "source": [
        "Like how libraries provide packaged functionality, we can create a **Class** that holds *attributes*  (or *properties*) and *methods* (and *functions*). \n",
        "\n",
        "This helps **abstract** a lot of detail by **ENCAPSULATING** *implementation* details and allowing us to re-use blocks of related code easily (i.e. collated pieces of logic forming a concept / object). \n",
        "\n",
        "You should aim to abstract away as much detail as necessary, but not more. \n",
        "\n",
        "> The process of re-visiting old code and changing it (*hopefully an improvement!*) is known as **refactorisation** (noun) or **refactoring** (verb).\n",
        "\n",
        "You can define an object by creating a **CLASS** definition using the syntax\n",
        "\n",
        "```\n",
        "class <MyClassName>(object):\n",
        "  \"\"\"some comment about the class\"\"\"\n",
        "\n",
        "  def __init__(self):\n",
        "    <things to do when object is created>\n",
        "\n",
        "  def <other class method>(self, a):\n",
        "    \"\"\"a short comment about the method\"\"\"\n",
        "    self.a  = a # this class now has property a!\n",
        "```\n",
        "\n",
        "> Note how the first parameter for a class method must be `self`. This provides the method with access to *other* class methods and properties. It's also a way to store arguments provided to methods as object properties. \n",
        "\n",
        "> `__init__` is a special method name (hence it starts *and* ends with `__`). It is called when a class is *initialised* or created (also known as **instantiation**) using `obj = MyClassName()`). `obj` is the **instantiated OBJECT** of the **class**. Hence, \"Object Orientated Programming\", or OOP.\n",
        "\n",
        "> Classes define the attributes and methods, Objects are what you play around with when you run your program.\n",
        "\n",
        "> There are a number of other special methods that classes have that can be overwritten!\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UwheocVbF5UP",
        "colab_type": "text"
      },
      "source": [
        "#### Encapsulate the neuron logic so far in a single entity"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tHZ-tL6Gfh35",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import math\n",
        "\n",
        "# create constants with uppercase variable names\n",
        "DEFAULT_V_REST = -65  # resting membrane potential (mV)\n",
        "DEFAULT_C_M = 1       # membrane capacitance (nF)\n",
        "DEFAULT_R_M = 100     # membrane resistance (MOhm)\n",
        "\n",
        "# define a 'point' neuron's core properties\n",
        "class MyFirstNeuron(object):\n",
        "  \"\"\" A point neuron that can receive spike times for GABAA and AMPA synapses, \n",
        "  as well as constant external current. Note this only models subthreshold \n",
        "  membrane potentials (no action potentials are generated).\n",
        "\n",
        "  Note: although a point neuron, the length and radius of the neuron must be\n",
        "  provided (in um).\n",
        "  \"\"\"\n",
        "  def __init__(self, length, radius, name=\"my first neuron\", \n",
        "               R_m=DEFAULT_R_M, C_m=DEFAULT_C_M, V_0=DEFAULT_V_REST):\n",
        "    \"\"\"\n",
        "    \"\"\"\n",
        "    print(\"creating a neuron\")\n",
        "    # initialise variables\n",
        "    self.name = name          # name now belongs, and is a property of \"MyFirstNeuron\"\n",
        "    self.length = length*1e-5 # um to dm\n",
        "    self.r = radius*1e-5      # um to dm\n",
        "    self.R_m = R_m\n",
        "    self.C_m = C_m\n",
        "    self.V_rest = V_0\n",
        "    self.I_e = 0\n",
        "    \n",
        "    # surface area of electrode placement (cm2)\n",
        "    self.A = (np.pi * self.length * self.r)*100              \n",
        "    # volume of sphere (dm3)     \n",
        "    self.volume = 4/3 * np.pi * (self.length/2) * self.r**2      \n",
        "    # surface area of open cylinder (cm2)\n",
        "    self.SA = (self.length*np.pi*self.r + 2*np.pi*self.r**2)*100  \n",
        "\n",
        "    # create lists, even if they start empty\n",
        "    self.V_t = []               # we are going to record the voltage\n",
        "    self.I_spks = []            # list for spikes that arrive at GABA synapses\n",
        "    self.E_spks = []            # list for spikes that arrive at AMPA synapses\n",
        "\n",
        "    # initialise variables for a simulation\n",
        "    self.reinit()\n",
        "\n",
        "  def add_inh_input(self, I_spks):\n",
        "    self.I_spks = I_spks\n",
        "\n",
        "  def add_exc_input(self, E_spks):\n",
        "    self.E_spks = E_spks\n",
        "\n",
        "  def add_external_input(self, I_e):\n",
        "    # current is in μA\n",
        "    self.I_e = I_e\n",
        "\n",
        "  def _I_syn(self, t, V_t):\n",
        "    \"\"\" Returns the sum of all synaptic currents.\n",
        "    \"\"\"\n",
        "    return I_GABA(t, V_t, self.I_spks) + I_AMPA(t, V_t, self.E_spks)\n",
        "\n",
        "  def step(self, t, dt):\n",
        "    V_t = self.V_t[-1] # get last value\n",
        "    I_leak = (V_t - self.V_rest)/self.R_m\n",
        "    dV = (-I_leak - self._I_syn(t, V_t) + self.I_e*1e-5/self.A) * dt/self.C_m\n",
        "    self.V_t.append(V_t + dV)\n",
        "\n",
        "  def reinit(self, V_0=None):\n",
        "    \"\"\"Re-initialise the membrane potential list, starting at V_0 (or V_rest \n",
        "    if no argument is provided).\n",
        "    \"\"\"\n",
        "    if V_0 is None:   # None is a special value that provides a good placeholder\n",
        "      V_0 = self.V_rest # use original value\n",
        "    self.V_t = [V_0]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rnwdjN-GEVKu",
        "colab_type": "text"
      },
      "source": [
        "> Note that methods that start with a single underscore *discourage* use    outside the class. E.g. `_I_syn`\n",
        "\n",
        "> Double underscore makes the method private (difficult to 'see' outside the  class and child classes)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CFL3IYqYF-ff",
        "colab_type": "text"
      },
      "source": [
        "#### Add some logic that allows us to 'step through' time.\n",
        "\n",
        "Objects in the simulator move through time through their `step` method."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oh5cDnIr_s6t",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class Simulator(object):\n",
        "  \"\"\" A Simulator can be used to step through time by calling each of its\n",
        "  encapsulated neuron's ``MyFirstNeuron.step`` method.\n",
        "  The duration and the time step can be chosen when call ``run`` or \n",
        "  ``continuefor`` (which does does initialise objects)\n",
        "  The time is recorded in the ``t_list`` variable for retrieval later.\n",
        "  \n",
        "  Note that only one Simulator ever exists.\n",
        "  We call this a Singleton class because efforts to construct a new object \n",
        "  instead retrieve the existing object and so there only exists a single object.\n",
        "  Here, this is implemented here by \n",
        "  - overriding the __new__ method (takes in `cls` - shorthand for 'class' - \n",
        "    instead of `self`)\n",
        "  - having a _sim class variable (it is at the top of the class definition)\n",
        "  other approaches can be found\n",
        "  https://python-3-patterns-idioms-test.readthedocs.io/en/latest/Singleton.html\n",
        "  \"\"\"\n",
        "  _sim = None\n",
        "\n",
        "  def __new__(cls):\n",
        "    if not cls._sim:\n",
        "      print(\"new Simulator\")\n",
        "      cls._sim = super(Simulator, cls).__new__(cls)\n",
        "      cls._sim.neurons = []\n",
        "    print(\"existing Simulator\")\n",
        "    return cls._sim\n",
        "  \n",
        "  def __init__(self):\n",
        "    self.init()\n",
        "\n",
        "  def add_neuron(self, neuron):\n",
        "    # ensure only neurons with unique names are in the list\n",
        "    for nrn in self.neurons:\n",
        "      if neuron.name == nrn.name:  # found a previous version\n",
        "        sim.neurons.remove(nrn)    # remove from simulator\n",
        "        del nrn                    # delete object\n",
        "    self.neurons.append(neuron)\n",
        "\n",
        "  def _run(self, duration, time_step):\n",
        "    \"\"\"private method to run the simulation from t_list[-1] until duration.\n",
        "    \"\"\"\n",
        "    print()\n",
        "    start = self.t_list[-1]\n",
        "    stop = self.t_list[-1]+duration+time_step\n",
        "    for t in np.arange(start, stop, time_step):\n",
        "      self.t_list.append(t)\n",
        "      for neuron in self.neurons:\n",
        "        neuron.step(t, time_step)\n",
        "      # if i%log_every==0:\n",
        "      #   print(\"#\"*int(100*i/num_iter) + f\"> {t} ms\", end='\\r')\n",
        "    print()\n",
        "\n",
        "  def run(self, duration=100, time_step=0.1):\n",
        "    \"\"\"Run a simulation over time, calling the step method of each neuron in the\n",
        "    neurons list. \n",
        "    \n",
        "    Values at t_0 are defined in all contained objects' reinit() methods\n",
        "    The first time step is therefore t_1 = dt.\n",
        "    \"\"\"\n",
        "    num_iter = int(duration/time_step)\n",
        "    print(f\"running a simulation for {duration} ms \"\n",
        "          f\"using a dt of {time_step} ms.\"\n",
        "          f\"\\nThere will be {num_iter} iterations for each of \"\n",
        "          f\"{len(self.neurons)} neurons\")\n",
        "    self.init()\n",
        "    self._run(duration, time_step)\n",
        "    print(\"done!\")\n",
        "  \n",
        "  def continuefor(self, duration=100, time_step=0.1):\n",
        "    num_iter = int(duration/time_step)\n",
        "    print(f\"continuing a simulation from {self.t_list[-1]} ms \"\n",
        "          f\"for {duration} ms \" \n",
        "          f\"using a dt of {time_step} ms.\"\n",
        "          f\"\\nThere will be {num_iter} iterations for each of \"\n",
        "          f\"{len(self.neurons)} neurons\")\n",
        "    self._run(duration, time_step)\n",
        "    print(\"done!\")\n",
        "\n",
        "  def init(self):\n",
        "    self.t_list = [0]\n",
        "    # reinitialise all neurons (for repeated simulations)\n",
        "    for neuron in self.neurons:\n",
        "      neuron.reinit()\n",
        "\n",
        "  def clean(self):\n",
        "    while len(self.neurons)>0:\n",
        "      neuron = self.neurons.pop()\n",
        "      del neuron\n",
        "    self.t_list = []    "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "s3f96XNg_3PF",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "T = 1000 # duration (ms)\n",
        "dt = 0.1 # time step (ms)\n",
        "\n",
        "# create a simulator\n",
        "sim = Simulator()\n",
        "\n",
        "# create a neuron with input\n",
        "nrn = MyFirstNeuron(23, 11, name='my 1st neuron!')\n",
        "nrn.add_inh_input(np.random.randint(0, T, size=10))\n",
        "nrn.add_exc_input(np.random.randint(0, T, size=10))\n",
        "\n",
        "# TODO: create ANOTHER neuron with external input\n",
        "# ...\n",
        "\n",
        "# adding all neurons we want to update\n",
        "sim.add_neuron(nrn)\n",
        "\n",
        "# run a simulation of a neuron\n",
        "sim.run(T, dt)\n",
        "\n",
        "# visualise\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# create a figure with an axis for plotting\n",
        "fig, ax = plt.subplots()\n",
        "\n",
        "# plot voltage over iterations\n",
        "ax.plot(sim.t_list, nrn.V_t)\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FW13rGMfr0PQ",
        "colab_type": "text"
      },
      "source": [
        "#### **Tasks**\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UHawHpxYxJKL",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "**Task 1**\n",
        "\n",
        "Run the cell below with different time steps. \n",
        "\n",
        "How does this affect the voltage trace?\n",
        "\n",
        "**Task 2**\n",
        "\n",
        "Create another neuron with external input and compare the traces\n",
        "\n",
        "we will use this neuron in future traces"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pIqdvVssCnIG",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title Task 2 Solution { display-mode: \"form\" }\n",
        "T = 1000 # duration (ms)\n",
        "dt = 0.1 # time step (ms)\n",
        "\n",
        "# create a simulator\n",
        "sim = Simulator()\n",
        "\n",
        "# create a neuron with input\n",
        "nrn = MyFirstNeuron(23, 11, 'my 1st neuron!')\n",
        "E_input = np.random.randint(0, T, size=10)\n",
        "I_input = np.random.randint(0, T, size=10)\n",
        "\n",
        "nrn.add_inh_input(I_input)\n",
        "nrn.add_exc_input(E_input)\n",
        "\n",
        "# TODO: create ANOTHER neuron with external input\n",
        "nrnE = MyFirstNeuron(23, 11, 'my 1st neuron with external input')\n",
        "nrnE.add_inh_input(I_input)\n",
        "nrnE.add_exc_input(E_input)\n",
        "nrnE.add_external_input(0.02)\n",
        "\n",
        "# adding all neurons we want to update\n",
        "sim.add_neuron(nrn)\n",
        "sim.add_neuron(nrnE)\n",
        "\n",
        "# run a simulation of a neuron\n",
        "sim.run(T, dt)\n",
        "\n",
        "# visualise\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# create a figure with an axis for plotting\n",
        "fig, ax = plt.subplots()\n",
        "\n",
        "# plot voltage over iterations\n",
        "names = []\n",
        "for n in sim.neurons:\n",
        "  ax.plot(sim.t_list, n.V_t)\n",
        "  names.append(n.name)\n",
        "\n",
        "names = [n.name for n in sim.neurons]\n",
        "ax.legend(names, frameon=False)\n",
        "\n",
        "# add some axis labels\n",
        "ax.set_ylabel(\"Voltage [$V_m$] (mV)\")\n",
        "ax.set_xlabel(\"Time (ms)\")\n",
        "\n",
        "# don't forget a title\n",
        "ax.set_title(\"How does external input affect voltage?\")\n",
        "\n",
        "# make the borders prettier\n",
        "ax.spines['right'].set_visible(False)\n",
        "ax.spines['top'].set_visible(False)\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8zUXuf54lDWY",
        "colab_type": "text"
      },
      "source": [
        "**Task 3**\n",
        "\n",
        "Instead of the **analytical** solution (can compute any time `t`) for `g_syn`, implement the **numerical** solution (need to know `t-1` to calculate `t`).\n",
        "\n",
        "$\\begin{aligned}\n",
        "\\frac{dg_{syn}}{dt} &= -\\frac{g_{syn}}{\\tau_{syn}} \\\\\n",
        "g_{syn} &= g_{syn} + g_{syn_{\\rm{max}}}   (\\textrm{when } t = t_{s})\n",
        "\\end{aligned}$\n",
        "\n",
        "Further questions\n",
        "1. What are the advantages + disadvantages between the 2 solutions?\n",
        "1. Derive the numerical solution from the analytical solution (ignore the $\\sum$)\n",
        "1. Integrate the numerical solution to get the analytical solution"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vJvdzD4DlCPu",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title Task 3 Solution { display-mode: \"form\" }\n",
        "class SynapseG(object):\n",
        "  \"\"\"Implemented synaptic conductance using Object Orientated Programming and\n",
        "  numerical solution to g_syn\n",
        "  \"\"\"\n",
        "  def __init__(self, spks, g_max=50, tau=5):\n",
        "    self.g_syn = 0\n",
        "    self.spks = sorted(spks) # sort the spikes (makes it a list)\n",
        "    self.g_max = g_max\n",
        "    self.tau = tau\n",
        "\n",
        "  def step(self, t, dt):\n",
        "    # decay any conductance\n",
        "    self.g_syn -= self.g_syn/self.tau * dt\n",
        "    # add conductance for incoming spikes\n",
        "    while len(self.spks) and t>=self.spks[0]:\n",
        "      self.spks.pop(0) # remove this spike from the list\n",
        "      self.presyn_spike()\n",
        "\n",
        "  def presyn_spike(self):\n",
        "    self.g_syn += self.g_max\n",
        "\n",
        "class Synapse(object):\n",
        "  def __init__(self, *args, E=0, **kwargs):\n",
        "    self.syn_g = SynapseG(*args,**kwargs)\n",
        "    self.E = E\n",
        "\n",
        "  def step(t, dt, V_t):\n",
        "    self.syn_g.step(t, dt)\n",
        "    self.g = self.syn_g.g_syn\n",
        "    self.I_syn = self.g*(V_t-self.E)\n",
        "\n",
        "\n",
        "class NeuronWithNumericalSynapse(object):\n",
        "\n",
        "  def __init__(self, length, radius, name=\"numerical synapses\", \n",
        "               R_m=DEFAULT_R_M, C_m=DEFAULT_C_M, V_0=DEFAULT_V_REST):\n",
        "    print(\"creating a neuron\")\n",
        "    # initialise variables\n",
        "    self.name = name\n",
        "    self.length = length\n",
        "    self.r = radius        \n",
        "    self.A = math.pi * length * radius\n",
        "    self.R_m = R_m\n",
        "    self.C_m = C_m\n",
        "    self.V_rest = V_0\n",
        "    self.I_e = 0\n",
        "\n",
        "    self.V_t = []               \n",
        "    self.reinit()\n",
        "\n",
        "  def add_inh_input(self, I_spks):\n",
        "    self.I_synapse = Synapse(I_spks, g_max=50, tau=8, E=-70)  # <-- changed\n",
        "\n",
        "  def add_exc_input(self, E_spks):\n",
        "    self.E_synapse = Synapse(I_spks, g_max=50, tau=4, E=0)  # <-- changed\n",
        "\n",
        "  def add_external_input(self, I_e):\n",
        "    self.I_e = I_e\n",
        "\n",
        "  def _I_syn(self, t, dt, V_t): # <-- method changed (call in `step` changed)\n",
        "    \"\"\" Returns the sum of all synaptic currents.\n",
        "    \"\"\"\n",
        "    self.I_synapse.step(t, dt, V_t)\n",
        "    self.E_synapse.step(t, dt, V_t)\n",
        "    return self.I_synapse.I_syn + self.E_synapse.I_syn\n",
        "\n",
        "  def step(self, t, dt):\n",
        "    V_t = self.V_t[-1] # get last value\n",
        "    I_leak = (V_t - self.V_rest)/self.R_m\n",
        "    dV = (-I_leak - self._I_syn(t, dt, V_t) + self.I_e/self.A) * dt/self.C_m \n",
        "    self.V_t.append(V_t + dV)\n",
        "\n",
        "  def reinit(self, V_0=None):\n",
        "    \"\"\"Re-initialise the membrane potential list, starting at V_0 (or V_rest \n",
        "    if no argument is provided).\n",
        "    \"\"\"\n",
        "    if V_0 is None:\n",
        "      V_0 = self.V_rest\n",
        "    self.V_t = [V_0]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WhwrRDnw-PBy",
        "colab_type": "text"
      },
      "source": [
        "## The Leaky Integrate-and-Fire Neuron"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IXpTUVkT-Ors",
        "colab_type": "text"
      },
      "source": [
        "Thus far, we have implemented the change in voltage given inputs, but there's no concept of an action potential going on!\n",
        "\n",
        "The simplest approach is to define a *threshold* $V_{thresh}$ for the neuron that once reached at time $t_i$\n",
        "\n",
        "1. $V_m$ goes to `0` (or `20`) at $t_{i+1}$\n",
        "1. After hyperpolarises to $V_{reset}$  at $t_{i+2}$\n",
        "\n",
        "Because we want to add extra functionality, we'll use this opportunity to introduce another 2 (of the 3) tenets of programming. **Inheritance** and **Polymorphism** (the other tenet was **encapsulation**)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TMTZyPVwKhhS",
        "colab_type": "text"
      },
      "source": [
        "### Inheritance\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Lch7V1p0MdCE",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "We already have a bunch of logic in `MyFirstNeuron` that has the same applicability to a leaky integrate-and-fire (LIF) neuron (setting initial values in `__init__`, adding inputs, summing synaptic inputs, and re-initialising the neuron's voltage). \n",
        "\n",
        "We can *inherit* this logic using \n",
        "\n",
        "```\n",
        "class <NewClass>(<ClassToInheritFrom>):\n",
        "```\n",
        "\n",
        "> note that previously we *inherited from object*. \n",
        "\n",
        "All methods that are in `<ClassToInheritFrom>` (e.g. `MyFirstNeuron`) are accessible in this new class!\n",
        "\n",
        "```\n",
        "class LIFNeuron(MyFirstNeuron):\n",
        "```\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RaniDxe-Qnif",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "But what if we want to *change* a method.\n",
        "\n",
        "For example, the `step` method in `MyFirstNeuron` doesn't have the logic to handle our (artificial) action potentials  -the \"fire\" 🔥 in \"integrate-and-fire\". It only has the sub-threshold \"integrate\" part. \n",
        "\n",
        "> Q: Where does the \"leaky\" come from?\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i5gzV5i0RdKr",
        "colab_type": "text"
      },
      "source": [
        "### Polymorphism"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QvDPEhL3QsII",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "To change this, we create a `step` method in this new class with the new logic we want. Using the same name to implement new functionality is known as **polymorphism**! Python determines which `step` method you want by the object that calls it (i.e. is it a `MyFirstNeuron` or `LIFNeuron`).\n",
        "\n",
        "In fact, the way we have implemented things, we do need to explicitly *repeat* some logic. We could go back to `MyFirstNeuron` and **refactor** `step` to make things smoother by re-using the same logic. But we're only repeating logic once so we'll accept this for now. If this happens again, go back and optimise to reduce sources of error and the number of places you would need to alter something.\n",
        "\n",
        "```\n",
        "def step(self, t, dt):\n",
        "  V_t = self.V_t[-1] # get last value\n",
        "  if V_t == 0:\n",
        "    # neuron has spiked at previous time step\n",
        "    V_t = self.V_rest\n",
        "  elif V_t >= V_thresh: # what if you change the order of the if statements?\n",
        "    # voltage has reached threshold, spike!\n",
        "    V_t = 0\n",
        "  else:\n",
        "    # this is repeated from MyFirstNeuron\n",
        "    I_leak = (V_t - self.V_rest)/self.R_m\n",
        "    dV = (-I_leak - self.__I_syn(t, V_t) + self.I_e/self.A) * dt/self.C_m\n",
        "    V_t = V_t + dV\n",
        "  self.V_t.append(V_t)\n",
        "```\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JicSGsIJRjyy",
        "colab_type": "text"
      },
      "source": [
        "### Polymorphism with Inheritance"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LFWDxqUNRVcT",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "We also have new parameters, `V_thresh` and `V_reset`, we want to introduce that only needs to be defined once when we create (**instantiate**) our neuron. We overwrite `__init__` by  creating a method with the same name, but then we **call** the **parent** (`MyFirstNeuron`) using the super special `super()` function, along with the method we want: `super().__init__()`. Because our parent method has parameters it accepts, we need to pass these along the hierarchy. We take in arbitrary *arguments* (e.g. `radius`) and *keyword arguments* (e.g. `R_m=1`) that other methods want with `*args` and `**kwargs`, represpectively. \n",
        "\n",
        "> We can explicitly accept arguments and pass them on too, e.g. maybe we want a different default value!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EIx1XJlv5OId",
        "colab_type": "text"
      },
      "source": [
        "### Another thing: [Logging](https://docs.python.org/3/library/logging.html)\n",
        "We have been using `print` to display text, but a better (more informative and robust) way is using the `logging` module."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XLLDq2xh5ydZ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title example\n",
        "import logging\n",
        "logging.basicConfig(format='%(asctime)s %(name)-15s %(levelname)-8s %(message)s',\n",
        "                    datefmt='%H:%M:%S',\n",
        "                    level=logging.INFO)\n",
        "\n",
        "logger = logging.getLogger('MyLogger')\n",
        "\n",
        "logger.setLevel(logging.DEBUG)\n",
        " \n",
        "logger.debug(\"Only show when debugging\")\n",
        "logger.warning(\"this displays a warning\")\n",
        "logger.error(\"there was an error\")\n",
        "logger.info(\"replace print statements with this instead\")\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y827oGSuRgCw",
        "colab_type": "text"
      },
      "source": [
        "### Implementation"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lDh9lxxn98ph",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class LIFNeuron(MyFirstNeuron):\n",
        "  \"\"\"\n",
        "  Leaky integrate-and-fire neuron that inherits from MyFirstNeuron.\n",
        "  \"\"\"\n",
        "  def __init__(self, *args, name=\"LIF\", \n",
        "               V_thresh=-55, V_reset=-68, **kwargs):\n",
        "    logging.debug(f\"other arguments passed = {args}\")\n",
        "    logging.debug(f\"other keyword arguments passed = {kwargs}\")\n",
        "    super().__init__(*args, name=name, **kwargs)\n",
        "    self.V_thresh = V_thresh\n",
        "    self.V_reset = V_reset\n",
        "\n",
        "  def step(self, t, dt):\n",
        "    V_t = self.V_t[-1] # get last value\n",
        "    if V_t == 0:\n",
        "      # neuron has spiked at previous time step\n",
        "      V_t = self.V_reset\n",
        "    elif V_t >= self.V_thresh: # what if you change the order of the if statements?\n",
        "      # voltage has reached threshold, spike!\n",
        "      V_t = 0\n",
        "    else:\n",
        "      # this is repeated from MyFirstNeuron\n",
        "      I_leak = (V_t - self.V_rest)/self.R_m\n",
        "      dV = (-I_leak - self._I_syn(t, V_t) + self.I_e*1e-5/self.A) * dt/self.C_m\n",
        "      V_t = V_t + dV\n",
        "    self.V_t.append(V_t)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Fg10qLgVFVLY",
        "colab_type": "text"
      },
      "source": [
        "Create a shortcut for visualisations"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WlSbtab1FXdU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# visualise\n",
        "%matplotlib inline\n",
        "import matplotlib.pyplot as plt\n",
        "from matplotlib.cbook import flatten\n",
        "\n",
        "def plot_v(sim, attrs=None, title=\"\", **kwargs):\n",
        "  # note that list parameters should have defaults as None in the method and \n",
        "  #   defined in the body, even if it's empty\n",
        "  if attrs is None:\n",
        "    attrs = []\n",
        "  \n",
        "  # create a figure with an axis for plotting\n",
        "  fig, axs = plt.subplots(nrows=1+len(attrs), squeeze=False, sharex=True,\n",
        "                          gridspec_kw={'height_ratios':[5]+[1]*len(attrs)})\n",
        "  ax = axs[0, 0]\n",
        "  # plot voltage over iterations\n",
        "  for neuron in sim.neurons:\n",
        "    ax.plot(sim.t_list, neuron.V_t, **kwargs)\n",
        "    for i, attr in enumerate(attrs):\n",
        "      try:\n",
        "        axs[i+1,0].plot(sim.t_list, getattr(neuron, attr))\n",
        "      except AttributeError:\n",
        "        logger.debug(f\"attribute {attr} did not exist in {neuron.name}. Continuing...\")\n",
        "        logger.debug(f\"\"\"do this next time `val = getattr(neuron, attr, None)`\n",
        "                          [see `help(getattr)`]\"\"\")\n",
        "\n",
        "  names = [n.name for n in sim.neurons]\n",
        "  ax.legend(names, loc='upper left', bbox_to_anchor=(1,1), frameon=False)\n",
        "\n",
        "  # add some axis labels\n",
        "  ax.set_ylabel(\"$V_m$ (mV)\") # <-- math symbols in-between $ $\n",
        "  axs[-1, -1].set_xlabel(\"Time (ms)\")\n",
        "\n",
        "  # don't forget a title\n",
        "  ax.set_title(title)\n",
        "\n",
        "  # make the borders prettier\n",
        "  for ax in flatten(axs):\n",
        "    ax.spines['right'].set_visible(False)\n",
        "    ax.spines['top'].set_visible(False)\n",
        "\n",
        "  return fig, axs"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BGuB-c_JJIoK",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "T = 1000 # duration (ms)\n",
        "dt = 0.1 # time step (ms)\n",
        "\n",
        "# create a neuron with input\n",
        "lif_nrn = LIFNeuron(23, 11)\n",
        "lif_nrn.add_inh_input(np.random.randint(0, T, size=10))\n",
        "lif_nrn.add_exc_input(np.random.randint(0, T, size=10))\n",
        "\n",
        "sim.add_neuron(lif_nrn)\n",
        "\n",
        "# run a simulation of a neuron\n",
        "sim.run(T, dt)\n",
        "\n",
        "plot_v(sim, title=\"Let there be spikes!\")\n",
        "print(f\"there were {np.sum(np.array(lif_nrn.V_t)>=0)} spikes\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9eGzo8DP5VOz",
        "colab_type": "text"
      },
      "source": [
        "### How does *external input current* affect the *firing rate* of a neuron?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iwIaHGNa7LNp",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "sim = Simulator()\n",
        "sim.clean()\n",
        "\n",
        "# create neurons each with their own external current value\n",
        "for current in np.arange(0,0.5,0.05):\n",
        "  nrn = LIFNeuron(23, 11, name=f\"{current:.2f}\")\n",
        "  nrn.add_external_input(current)\n",
        "  sim.add_neuron(nrn)\n",
        "sim.run(T, dt)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YcvzeRzcEsOU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "fig, ax = plot_v(sim, title=\"Let there be spikes!\", alpha=0.7, lw=1)\n",
        "ax[0,0].get_legend().set_title(\"External\\nCurrent (μA)\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "F_W775u4CXQp",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# convert Voltage for each neuron to an array\n",
        "v = np.empty(shape=(len(sim.t_list),len(sim.neurons)))\n",
        "for i,nrn in enumerate(sim.neurons):\n",
        "  v[:,i] = nrn.V_t\n",
        "# calculate the instantaneous firing rate\n",
        "ifr = np.sum(v>=0,axis=0)/(sim.t_list[-1]/1000)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Smtlt98vHl9S",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "fig, ax = plt.subplots()\n",
        "ax.plot([float(nrn.name) for nrn in sim.neurons],ifr)\n",
        "ax.set_xlabel(\"Current (μA)\")\n",
        "ax.set_ylabel(\"Firing rate (Hz)\")\n",
        "ax.set_title(\"I-F 'curve'\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "06LikNcZLeNj",
        "colab_type": "text"
      },
      "source": [
        "### **Homework 1**: How does _**excitatory** and **inhibitory** input_ affect the *firing rate* of a neuron?\n",
        "\n",
        "Similar to above for external current, determine\n",
        "- how does excitatory drive affect firing rate? \n",
        "\n",
        "  you should produce an input-output (I-O) curve where excitation is on the x-axis and firing rate is on the y-axis\n",
        "\n",
        "- how does inhibition affect this I-O relationship?\n",
        "\n",
        "- how does timing of the excitatory + inhibitory inputs affect the relationship? \n",
        "  - what if they both have the same input frequency (\"Balanced input\")\n",
        "  - what if they have the same frequency, but are offset?\n",
        "  - what does noise add or take away from this I-O function?\n",
        "\n",
        "- create a circuit with 2 excitatory neurons and 1 inhibitory neuron. Connect them to each other. \\\n",
        "  *hint when one fires, add a spike time to its connected pair*\n",
        "  - Increase the number of neurons (try keep a 80:20 percentage of excitatory: inhibitory neurons).\n",
        "  - Implement Hebbian learning"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "btB5clIQ9fM7",
        "colab_type": "text"
      },
      "source": [
        "## Ions 🏋️\n",
        "![ions](https://images.unsplash.com/photo-1531748049999-7217dc4cfbd1?ixlib=rb-1.2.1&ixid=eyJhcHBfaWQiOjEyMDd9&auto=format&fit=crop&w=1275&h=500&q=80)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iDZ2rKjaoknX",
        "colab_type": "text"
      },
      "source": [
        "### Equilibrium potentials\n",
        "\n",
        "\n",
        "In general, we create models that ignore (\"abstract away\") unncessary complexity but capture what we're investigating [1]. \n",
        "\n",
        "1. Herz AVM, Gollisch T, Machens CK, Jaeger D. Modeling single-neuron dynamics and computations: A balance of detail and abstraction. Science (80- ). 2006;314(5796):80–5. doi: [10.1126/science.1127240](http://doi.org/10.1126/science.1127240) [[DOWNLOAD PDF](https://drive.google.com/open?id=1WcUPYYdBcDpXbFcR5JSLd5sp1FvmsWo0)]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "l1FfFW7Czr8v",
        "colab_type": "text"
      },
      "source": [
        "#### Where does the resting membrane potential come from?\n",
        "\n",
        "The resting membrane potential is determined by the uneven distribution of ions (charged particles) between the inside and the outside of the cell, and by the different permeability of the membrane to different types of ions.\n",
        "\n",
        "<div>\n",
        "<figure style=\"float:left\"><img src=\"https://www.news-medical.net/image.axd?picture=2018%2f10%2fshutterstock_480412786.jpg&ts=20181025104435&ri=673\" width=\"600px\"></figure>\n",
        "<img style=\"float:left\" src=\"https://upload.wikimedia.org/wikipedia/commons/f/fb/Basis_of_Membrane_Potential2.png\" width=\"300px\">\n",
        "<img style=\"float:right\" src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/7/73/Action_potential_ion_sizes.svg/1920px-Action_potential_ion_sizes.svg.png\" width=\"300px\">\n",
        "</div>\n",
        "\n",
        "- [Wikipedia](https://en.wikipedia.org/wiki/Membrane_potential)\n",
        "- [Khan Academy](https://www.khanacademy.org/science/biology/human-biology/neuron-nervous-system/a/the-membrane-potential)\n",
        "\n",
        "Despite the small differences in their radii, ions rarely go through the \"wrong\" channel."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ONtYUXG__TBz",
        "colab_type": "text"
      },
      "source": [
        "#### Implementation\n",
        "The way we are going to implement ions is going to be simple and doesn't require the complexity of class and objects (which we could certainly do). We are going to using a **Dictionary**."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bBxk9ccNJiNr",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "A **Dictionary** is an unordered collection of `key:value` pairs. This means the index is the **key** instead of an integer as in a `list`.\n",
        "They are created with curly brackets `{ }` or the `dict` keyword.\n",
        "\n",
        "```\n",
        "# create\n",
        "my_dict = {'key': 'value',\n",
        "           'key2': 'another value, \n",
        "           'b': 2, \n",
        "           10: 'j',\n",
        "           'list_val': my_list\n",
        "           }\n",
        "# add another key: value pair\n",
        "my_dict['new key'] = {'a': 'dict',\n",
        "                      'within': True}\n",
        "\n",
        "# overwrite a key (reassign)\n",
        "my_dict['key'] = dict(another_dict_key=\"another dict's value\")\n",
        "\n",
        "# retrieve a key\n",
        "print(my_dict['key2'])\n",
        "\n",
        "# iterate\n",
        "for k, v in my_dict.items():\n",
        "  print(f\"my_dict[{k}] = {v}\")\n",
        "```\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HxrEf0bjBGQ-",
        "colab_type": "text"
      },
      "source": [
        "> the `dict` method of creation looks *suspiciously* like keyword arguments in methods 😉"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uFbHPcoVa0k9",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#-----------------------------\n",
        "# Variables\n",
        "#-----------------------------\n",
        "C = concentrations = {'K':    {'i':140,'o':  5},\n",
        "                      'Na':   {'i': 15,'o':125},\n",
        "                      'Cl':   {'i':  5,'o':135},\n",
        "                      'HCO3': {'i': 10,'o': 25},\n",
        "                      }\n",
        "\n",
        "currents = list(concentrations.keys()) # get the keys and convert to a list\n",
        "currents  # the exact order of keys is *not* guaranteed!\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZIuZX-fO3vMN",
        "colab_type": "text"
      },
      "source": [
        "#### Reversal Potential"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2nlGi_mpRn1F",
        "colab_type": "text"
      },
      "source": [
        "> In a biological membrane, the reversal potential (also known as the Nernst potential) of an ion is the membrane potential at which there is no net (overall) flow of that particular ion from one side of the membrane to the other. In the case of post-synaptic neurons, the reversal potential is the membrane potential at which a given neurotransmitter causes no net current flow of ions through that neurotransmitter receptor's ion channel.\n",
        "  https://en.wikipedia.org/wiki/Reversal_potential\n",
        "\n",
        "For a single ion species, the Nernst potential is used,\n",
        "\n",
        "$ E = \\frac{R \\cdot T}{z \\cdot F} \\cdot \\ln(\\frac{[ion]_{out}}{[ion]_{in}})$\n",
        "\n",
        "Where \n",
        "- $[ion]_{out}$ is the extracellular concentration of that ion (in moles per cubic meter, to match the other SI units, though the units strictly don't matter, as the ion concentration terms become a dimensionless ratio),\n",
        "- $[ion]_{in}$ is the intracellular concentration of that ion (in moles per cubic meter),\n",
        "- R is the ideal gas constant (joules per kelvin per mole),\n",
        "- T is the temperature in kelvins,\n",
        "- F is Faraday's constant (coulombs per mole).\n",
        "- z is the charge of the ion species (e.g. +1 for $Na^+$ and -1 for $Cl^-$)\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4SRoeGGGTyk0",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#-----------------------------\n",
        "# Constants\n",
        "#-----------------------------\n",
        "R = 8.314   # Real gas constant (J K^-1 mol^-1)\n",
        "F = 96485   # Faraday constant (C mol^-1)\n",
        "T = 310.25  # Temperature (K)\n",
        "RTF = R*T/F # All treated as constant for the simulations\n",
        "\n",
        "z = valence = {'K': 1, 'Na': 1, 'Cl': -1,'HCO3': -1}"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MtHGBMdYQrDM",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def nernst(ion, conc=None):\n",
        "  \"\"\"Calculate the Nernst potential (in mV) for a given ion.\n",
        "     The ion concentrations are retrieved from the conc dict.\n",
        "     The ion valences are retrieved from the valence dict.\n",
        "     >>> round(nernst('Cl', conc=dict(Cl=dict(i=5, o=135))),2)\n",
        "     -88.11\n",
        "  \"\"\"\n",
        "  if conc is None:\n",
        "    conc=concentrations\n",
        "  return RTF/valence[ion] * np.log(conc[ion]['o']/conc[ion]['i'])*1000 # note that log is the natural log (base e not base 10)\n",
        "\n",
        "def nernst_ind(C_out, C_in, sign):\n",
        "  \"\"\"Calculate the Nernst potential (in mV) for given C out, C in, z\n",
        "  >>> round(nernst_ind(135, 5, -1),2)\n",
        "  -88.11\n",
        "  \"\"\"\n",
        "  return RTF/sign * np.log(C_out/C_in)*1000"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8ZNoD4k04Y1V",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# some simple tests\n",
        "for ion in concentrations:\n",
        "  _e = nernst(ion)\n",
        "  print(f\"E{ion:<4} = {_e:>8.2f} mV\")\n",
        "  assert _e == nernst_ind(C[ion]['o'], C[ion]['i'], valence[ion]), \"Nernst error\""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FKzRR6ag4Ng6",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# remember to test more explicitly\n",
        "import doctest\n",
        "doctest.run_docstring_examples(nernst, globals(), \n",
        "                               name=nernst.__name__, verbose=True)\n",
        "doctest.run_docstring_examples(nernst_ind, globals(), \n",
        "                               name=nernst_ind.__name__, verbose=True)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ERvH9pfOT8z8",
        "colab_type": "text"
      },
      "source": [
        "#### Resting potential\n",
        "\n",
        "To calculate the equilibrium potential for section of membrane that is permeable to *multiple* ion species, the Goldman-Hodgkin-Katz (GHK) equation is used.\n",
        "\n",
        "This is true for some receptors like the $GABA_A$ receptor which is permeable to both Chloride and Bicarbonate ions. \n",
        "\n",
        "Because the neuronal membrane has passive channels and active ion pumps that make it permeable to ions even when not \"doing anything\" (spiking), the GHK equation is used to calculate the resting membrane potential.\n",
        "\n",
        "![](https://wikimedia.org/api/rest_v1/media/math/render/svg/e290b78cf2f968293d8dde2a4732c6d22c6d1226)\n",
        "\n",
        "Where \n",
        "- $P$ is the relative permeability of that ion (or selectivity in meters per second)\n",
        "- $M$ are monovalent cations (*positive* ions with a single charge)\n",
        "- $A$ are monovalent anions (*negative* ions with a single charge)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kbzOwA_3T6yX",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#-----------------------------\n",
        "# Parameters\n",
        "#-----------------------------\n",
        "# GABAA receptor\n",
        "pcl = 0.8   # Cl permeability\n",
        "phco3 = 0.2 # HCO3 permeability\n",
        "\n",
        "# membrane permeability to ions (multiple values have been found in the lit)\n",
        "pK = 1\n",
        "pNa = 0.05\n",
        "pCl = 0.45\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tNiZtMN44Kq-",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def ghk(C_outs, C_ins, ps, zs):\n",
        "  \"\"\"Calculate the potential of given ions using the Goldman–Hodgkin–Katz (GHK)\n",
        "  equation.\n",
        "  >>> round(ghk([4,125,110], [150,15,10], [pK, pNa, pCl], [1,1,-1]),2)\n",
        "  -69.73\n",
        "  \"\"\"\n",
        "  dividend = 0\n",
        "  divisor = 0\n",
        "  for cin, cout, p, z  in zip(C_ins, C_outs, ps, zs):\n",
        "    assert abs(z) == 1, \"only monovalent ions supported\"\n",
        "    if z>0:\n",
        "      dividend += p*cout\n",
        "      divisor += p*cin\n",
        "    else:\n",
        "      dividend += p*cin\n",
        "      divisor += p*cout\n",
        "  return RTF*np.log(dividend/divisor)*1000\n",
        "\n",
        "# some simple tests\n",
        "for ion in concentrations:\n",
        "  _e = nernst(ion)\n",
        "  print(f\"E{ion:<4} = {_e:>8.2f} mV\")\n",
        "  _e_ghk = ghk([C[ion]['o']], \n",
        "             [C[ion]['i']], \n",
        "             [1], [valence[ion]])\n",
        "  assert _e == nernst_ind(C[ion]['o'], C[ion]['i'], valence[ion]), \"Nernst error\"\n",
        "  assert round(_e,8) == round(_e_ghk,8), f\"Nernst ({_e}) != GHK ({_e_ghk})\""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cl5wmEYLwQoI",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "doctest.run_docstring_examples(ghk, globals(), \n",
        "                               name=ghk.__name__, verbose=True)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8-AkipZCwfFQ",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "> ### **Homework 2** Investigate the reversal potential for $GABA_A$ synapses ($E_{GABA}$) as function of its permeability to its ions: $Cl^-$ and $HCO_3^-$. \n",
        "  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "g0xrjInYtkiQ",
        "colab_type": "text"
      },
      "source": [
        "![GABAA synapse](https://n6adig.dm.files.1drv.com/y4m1rMjfg2Tg0NwOpdM_nzOAULV_h6UxibvnoWXYliG8Sg-A8--d_smNQeA4UiprJbjkuLMSZZG2Ir6X84gUo_8ofgDQA4HXf3b8AaR3YP-AQfGyWaiwbkEsl7SJCFCYaADmezTUA3ZlvCNWbrN7U1ZyQz6DGMyfFRW1o8vWEXFo1UB4cflv2ATi-MTF60PRZWQ9enjYHUcYTl2lzCNqOiy6A?width=720&height=1040&cropmode=none) - _C.B. Currin, 2020_"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Dmz4UP-7jw81",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "#### **2.1**\n",
        "  Experiments have shown that $GABA_A$ can be excitatory ($V_m - E_{GABA} < 0$ ). \n",
        "  \n",
        "  Given \n",
        "  $\\begin{align}\n",
        "  [Cl^-]_o &= 135 mM\\\\\n",
        "  [HCO_3^-]_o &= 25 mM \\\\\n",
        "  [HCO_3^-]_i &= 10 mM \\\\\n",
        "  V_m &= -65 mV\n",
        "  \\end{align}$, what value for $[Cl^-]_i$ elicits a negative *driving force* $(V_m - E_{GABA})$? \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JDFoHI6Jj2_N",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "#### **2.2**\n",
        "The change in concentration for an ion based on its current, can be formulated as $\\frac{d[X]}{dt} = I_{[X]} \\cdot \\frac{1}{F \\cdot \\rm{volume}}$. \n",
        "  \n",
        "  Given \n",
        "  $\\begin{align}\n",
        "  \\rm{radius} &= 5 \\mu m \\\\\n",
        "  \\rm{length} &= 12 \\mu m \\\\\n",
        "   g_{GABA_{\\rm{max}}} &= 500 nS\n",
        "  \\end{align}$\n",
        "\n",
        " what frequency of spikes from an interneuron (i.e. activation of $GABA_A$ synapses) causes $E_{GABA} > V_m$ within $1000 ms$.\n",
        "\n",
        " *hint: the GHK equation won't work (why not?)*\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E0nU8r_ej6Ik",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "#### **2.3**\n",
        "Implement chloride extrusion\n",
        "\n",
        " $\\nu = \\frac{[Cl^-]_i - [Cl^-]_i^\\infty}{\\tau_{Cl^-}}$; where $\\nu$ is in mM/s. Do **1.2** again but with $\\tau_{Cl^-} = 3 s$ or $\\tau_{Cl^-} = 30 s$\n",
        "\n",
        "What $\\tau_{Cl^-}$ keeps $E_{GABA} < V_m$ for a 5 Hz inhibitory input?\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NvnlAVdgscXf",
        "colab_type": "text"
      },
      "source": [
        "![Example solution](https://3uazqw.dm.files.1drv.com/y4mQWEZlEfWeIsKIFokTXSY-J9mt1cxv7UhzM1hTGl4h5jDWII1ggBLAUlGDOqRQ2yXxKDaNWwuRPaxD4d9k1gVL0UrOCtaVnmZYY13rkijBmafqu1U-RX6g1sOoSJkwzc1gHkNPft8cnuymDkd6-zCIbNHG5GkQOJ0nuz2CdSkOGFHKD9gzYWBDebyD5omSzCNqbgTsC5WUuFFrmdSv3UC9A?width=798&height=360&cropmode=none) - _C.B. Currin, 2020_"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B1UMC3vNn_EH",
        "colab_type": "text"
      },
      "source": [
        "# More\n",
        "\n",
        "1. Take a breath. Take another. And many more.\n",
        "\n",
        "1. High-five your neighbour: you've made it past the first hurdle and it's ok to not understand everything. \n",
        "> *pssst, programmers are great at using Google!*\n",
        "\n",
        "1. coding is just problem solving using a formal *language* (in this case python) to convert instructions we give into some form of processing. \n",
        "1. **most** of programming is figuring out which instructions are possible/allowed (aka the **syntax**) and which are appropriate for your problem (aka the **semantics**) - just because something compiles (no syntax errors) does not mean it is correct in solving the problem (or the best way to solve a problem). \n",
        "\n",
        "  \"A compiler or interpreter could complain about syntax errors. Your co-workers will complain about semantics.\"\n",
        "\n",
        "1. Test often. Test well.\n",
        "1. [read and contribute to **The role of scientific code:**](https://docs.google.com/document/d/11NcnnCllrHD19EpHc2wSjdD_5mimqClDkdmweZkOPrk/edit)\n",
        "How can we write better code, i.e. use our laboratory more effectively?\n",
        "1. [The Zen of Python](https://www.python.org/dev/peps/pep-0020/#id3)\n",
        "1. Use a good file structure (https://drivendata.github.io/cookiecutter-data-science)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2aHDCIdl3h1i",
        "colab_type": "text"
      },
      "source": [
        "```\n",
        "├── LICENSE\n",
        "├── Makefile           <- Makefile with commands like `make data` or `make train`\n",
        "├── README.md          <- The top-level README for developers using this project.\n",
        "├── data\n",
        "│   ├── external       <- Data from third party sources.\n",
        "│   ├── interim        <- Intermediate data that has been transformed.\n",
        "│   ├── processed      <- The final, canonical data sets for modeling.\n",
        "│   └── raw            <- The original, immutable data dump.\n",
        "│\n",
        "├── docs               <- A default Sphinx project; see sphinx-doc.org for details\n",
        "│\n",
        "├── models             <- Trained and serialized models, model predictions, or model summaries\n",
        "│\n",
        "├── notebooks          <- Jupyter notebooks. Naming convention is a number (for ordering),\n",
        "│                         the creator's initials, and a short `-` delimited description, e.g.\n",
        "│                         `1.0-jqp-initial-data-exploration`.\n",
        "│\n",
        "├── references         <- Data dictionaries, manuals, and all other explanatory materials.\n",
        "│\n",
        "├── reports            <- Generated analysis as HTML, PDF, LaTeX, etc.\n",
        "│   └── figures        <- Generated graphics and figures to be used in reporting\n",
        "│\n",
        "├── requirements.txt   <- The requirements file for reproducing the analysis environment, e.g.\n",
        "│                         generated with `pip freeze > requirements.txt`\n",
        "│\n",
        "├── setup.py           <- Make this project pip installable with `pip install -e`\n",
        "├── src                <- Source code for use in this project.\n",
        "│   ├── __init__.py    <- Makes src a Python module\n",
        "│   │\n",
        "│   ├── data           <- Scripts to download or generate data\n",
        "│   │   └── make_dataset.py\n",
        "│   │\n",
        "│   ├── features       <- Scripts to turn raw data into features for modeling\n",
        "│   │   └── build_features.py\n",
        "│   │\n",
        "│   ├── models         <- Scripts to train models and then use trained models to make\n",
        "│   │   │                 predictions\n",
        "│   │   ├── predict_model.py\n",
        "│   │   └── train_model.py\n",
        "│   │\n",
        "│   └── visualization  <- Scripts to create exploratory and results oriented visualizations\n",
        "│       └── visualize.py\n",
        "│\n",
        "└── tox.ini            <- tox file with settings for running tox; see tox.testrun.org\n",
        "```"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TUmsKMHPv1bX",
        "colab_type": "text"
      },
      "source": [
        "## The beauty of the brain's complexity\n",
        "\n",
        "Whatever we've learnt so far, it's **always** _slighty_ more complicated than that.\n",
        "\n",
        "You get so many types of neurons of all shapes and sizes. \n",
        "\n",
        "![Neurons](http://www.scholarpedia.org/w/images/thumb/e/e4/Neuron_Llinas_photographs.jpeg/356px-Neuron_Llinas_photographs.jpeg)\n",
        "\n",
        "Thus far we haven't considered a single \"type\" *per se*, but generally when neuroscientists model neurons they model pyramidal cells (**PC**) which are commonly found in the cortex (layer 2/3 & layer 5) and hippocampus (CA1, CA3). \n",
        "\n",
        "Inhibition is normally modelled as Parvalbumin+ (PV; the neurons stain positive for the protein parvalbumin when exposed to the right antibodies) interneurons (**IN**).\n",
        "[Many others exist, even just in the cortex](https://knowingneurons.com/2014/11/05/inhibitory-neurons-keeping-the-brains-traffic-in-check/): ![](https://knowingneurons.files.wordpress.com/2014/11/interneuron_650.jpg?resize=610%2C809)\n",
        "\n",
        "\n",
        "Many parameters are chosen based on experimental recordings from non-human mammals. Keep in mind there are known differences between the human cortex and a rat cortex, for example.\n",
        "\n",
        "Just remember:\n",
        "\n",
        "> All models are wrong, but some are useful. - George E. Box\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BABSiFs91hLY",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "# Analysis"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RX-vtmE3Lmv0",
        "colab_type": "text"
      },
      "source": [
        "## Towards better data management and plotting: Pandas and Seaborn\n",
        "\n",
        "Pandas allows the creation of `DataFrames` that neatly capture multi-dimensional arrays that can be more easily understood and visualised. \n",
        "\n",
        "Seaborn neatly plots tidy (“long-form”) dataframes where each column is a variable and each row is an observation. [Tidy Data in Python](https://www.jeannicholashould.com/tidy-data-in-python.html)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B_UhkcuQsO5K",
        "colab_type": "text"
      },
      "source": [
        "Pandas has some awesome ways to load data.\n",
        "\n",
        "E.g. read from an excel spreadsheet (can even reference an online document)\n",
        "\n",
        "[`df = pd.read_excel(...)`](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.read_excel.html)\n",
        "\n",
        "E.g. save to a \"comma-separated value\" file\n",
        "\n",
        "[`df.to_csv(\"file name.csv\")`](https://dev.pandas.io/docs/user_guide/io.html#io-store-in-csv)\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kzsOHGff1fuf",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import seaborn as sns\n",
        "import pandas as pd\n",
        "\n",
        "# seaborn has some built-in datasets\n",
        "dots = sns.load_dataset(\"dots\")\n",
        "# display a pandas dataframe, loaded from data\n",
        "print(type(dots))\n",
        "dots"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BscX0nrLDjX2",
        "colab_type": "text"
      },
      "source": [
        "### Converting between long and wide data\n",
        "[Reshaping](https://pandas.pydata.org/pandas-docs/stable/user_guide/reshaping.html)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YpvZT63xJp5J",
        "colab_type": "text"
      },
      "source": [
        "- to wide:\n",
        "  - [`pivot`](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.pivot.html#pandas.DataFrame.pivot) or [`pivot_table`](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.pivot_table.html#pandas.DataFrame.pivot_table)\n",
        "\n",
        "  <img src=\"https://pandas.pydata.org/pandas-docs/stable/_images/reshaping_pivot.png\" width=\"500px\">\n",
        "\n",
        "  - [`stack`](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.unstack.html?highlight=unstack#pandas.DataFrame.stack)\n",
        "\n",
        "  <img src=\"https://pandas.pydata.org/pandas-docs/stable/_images/reshaping_stack.png\" width=\"500px\">\n",
        "\n",
        "- to long\n",
        "  - [`melt`](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.melt.html#pandas.DataFrame.melt)\n",
        "\n",
        "  <img src=\"https://pandas.pydata.org/pandas-docs/stable/_images/reshaping_melt.png\" width=\"500px\">\n",
        "\n",
        "  - [`unstack`](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.unstack.html?highlight=unstack#pandas.DataFrame.unstack)\n",
        "\n",
        "  <img src=\"https://pandas.pydata.org/pandas-docs/stable/_images/reshaping_unstack.png\" width=\"500px\">\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8yR9-qLBIPCW",
        "colab_type": "text"
      },
      "source": [
        "#### Example"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AtV5lSYAaRye",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "wide_df = dots.pivot_table(index='time', columns=['align','coherence','choice'], values=['firing_rate'])\n",
        "wide_df"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mcHRnetLFOi-",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# first make a copy so we can manipulate it before converting without affecting the original\n",
        "time_as_col = wide_df.copy()\n",
        "# make the index an explicit column (melt tends to ignore index)\n",
        "time_as_col['time'] = time_as_col.index\n",
        "# reshape\n",
        "long_df = time_as_col.melt(id_vars=['time'], value_vars=['firing_rate'], value_name='firing_rate')\n",
        "# show a sample\n",
        "print(long_df.head())\n",
        "print(f\"(rows, columns) = {long_df.shape}\")\n",
        "# clean up by removing NaNs and removing the NaN (None) column\n",
        "long_df = long_df.dropna().drop(columns=[None])\n",
        "long_df"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_KPPX3i8K0KZ",
        "colab_type": "text"
      },
      "source": [
        "### Plotting long and wide data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CjhhUvzoK3R7",
        "colab_type": "text"
      },
      "source": [
        "##### Plot wide data using Pandas and Matplotlib"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "STDU-_rwLFqk",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#simplest version (with 2 columns)\n",
        "fig, ax = plt.subplots(1, 2, sharey=True)\n",
        "for _ax, align in zip(ax,['dots','sacc']):\n",
        "  dropped_df = wide_df['firing_rate',align].dropna(how='all')\n",
        "  style=['-',':']*int(dropped_df.columns.get_level_values('choice').size/2)\n",
        "  dropped_df.plot(style=style, ax=_ax, cmap='viridis', legend=align=='sacc')\n",
        "  _ax\n",
        "\n",
        "# annotate\n",
        "ax[0].set_ylabel('firing_rate')\n",
        "ax[1].legend(title=ax[1].get_legend().get_title().get_text(), frameon=False,\n",
        "             loc='upper left', bbox_to_anchor=(1,1))\n",
        "# clean up\n",
        "for _ax in ax:\n",
        "  for s,spine in _ax.spines.items():\n",
        "    if s=='top' or s=='right':\n",
        "      spine.set_visible(False)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JTH0et0uNwWd",
        "colab_type": "text"
      },
      "source": [
        "more control over the plotting\n",
        "\n",
        "what we want is \n",
        "- each `align` to be on different axes\n",
        "- each ` choice` to be a different line width\n",
        "- each `coherence` to be a different colour\n",
        "\n",
        "we iterate over slices of the `pd.DataFrame` (where each slice is a `pd.Series` object with a flattened index (`name` in the for loop) \n",
        "\n",
        "`name=(<align>,<coherence>,<choice>)` and plot these slices.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "50wwpG-2goQE",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import matplotlib\n",
        "from matplotlib import cm\n",
        "\n",
        "# Change ax to plot on by 'align'\n",
        "# Get all alignments (set is used to only have unique values)\n",
        "align = set(wide_df.columns.get_level_values('align'))\n",
        "ncols = len(align)\n",
        "fig, ax = plt.subplots(1,ncols, sharey=True)\n",
        "# zip combines 2 iterables (like list or set) such that \n",
        "#   zip([a,b],[1,2]) == [(a,1),(b,2)]\n",
        "#   (zip returns anobject generator, but wrap list around it for above)\n",
        "# when we pass zip to dict, it converts the list of tuples () to key:value pairs\n",
        "align_to_ax = dict(zip(sorted(align),ax))\n",
        "\n",
        "# Change linewidth by 'choice' - `[::-1]` reverses a list\n",
        "choices = sorted(set(wide_df.columns.get_level_values('choice')))[::-1]\n",
        "lw_min = 1\n",
        "lw_scale = 5\n",
        "choice_to_lw = dict(zip(choices,\n",
        "                        np.arange(lw_min,len(choices)*lw_scale+lw_min,lw_scale)))\n",
        "\n",
        "# Change color by 'coherence'\n",
        "coherence = sorted(set(wide_df.columns.get_level_values('coherence')))\n",
        "norm = matplotlib.colors.Normalize(vmin=min(coherence), vmax=max(coherence))\n",
        "# like list comphrehension, we an do dict comprehension\n",
        "coh_to_color = {coh:cm.magma_r(norm(coh)) for coh in coherence}\n",
        "\n",
        "# iterate over slices of the dataframe\n",
        "for name, series in wide_df['firing_rate'].iteritems():\n",
        "  # print(name)\n",
        "  series.dropna(how='all').plot(ax=align_to_ax[name[0]], \n",
        "              lw=float(choice_to_lw[name[2]]), \n",
        "              c=coh_to_color[name[1]],legend=False)\n",
        "  \n",
        "# annotate\n",
        "for align,_ax in align_to_ax.items():\n",
        "  _ax.set_title(f'align = {align}')\n",
        "ax[0].set_ylabel('firing_rate')\n",
        "\n",
        "from matplotlib.lines import Line2D\n",
        "# create legend for choice - list of tuples (Line, label)\n",
        "leg_choice = [(Line2D([],[],c='k', lw=choice_to_lw[choice]), choice) for choice in choices][::-1]\n",
        "\n",
        "# *zip(*leg_choice) works by converting [(line1,label1), (line2, label2)] to \n",
        "#   [[line1, line2], [label1, label2]] - this is zip(*leg_choice)\n",
        "# legend expects the first argument to be lines and the 2nd labels, so we spread\n",
        "#   the list of 2 lists [[],[]] using the * operator to only have 2 lists [],[].\n",
        "leg = ax[1].legend(*zip(*leg_choice), title='choice',\n",
        "                   loc='lower left', bbox_to_anchor=(1,-0.1), ncol=2, frameon=False)\n",
        "# use the existing dict, which can be problematic if ordering matters\n",
        "#   (see comment below)\n",
        "leg_coh = [(Line2D([],[],c=c), coh) for coh,c in coh_to_color.items()]\n",
        "# note that dicts are unsorted\n",
        "# we can use something like `from collections import OrderedDict` and some \n",
        "#   settings to maintaing an ordering.\n",
        "# here we use an anonymous function on the leg_coh list\n",
        "leg_coh = sorted(leg_coh, key=lambda linelabel: (linelabel[1],linelabel[0]))\n",
        "# we lose the previous legend when we call legend() again\n",
        "ax[1].legend(*zip(*leg_coh), title='coherence',\n",
        "             loc='upper left', bbox_to_anchor=(1,1.1), frameon=False)\n",
        "# but we can have both by explicitly adding the object back\n",
        "ax[1].add_artist(leg)\n",
        "\n",
        "# clean up axes\n",
        "for _ax in ax:\n",
        "  for s,spine in _ax.spines.items():\n",
        "    if s=='top' or s=='right':\n",
        "      spine.set_visible(False)\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lfDI66y6Y_3N",
        "colab_type": "text"
      },
      "source": [
        "we learnt some **new** Python above\n",
        "> `tuple` - an immutable iterable akin to a list\n",
        "  \\\n",
        "  *immutable* means you cannot change the data after you have created the object. You can reassign the name, though. \n",
        "  \\\n",
        "  this will fail\n",
        "  ```\n",
        "  t = (1,5,10)\n",
        "  t[0] = 2\n",
        "  ```\n",
        "  hint: you can use tuples as default arguments in methods directly (instead of using `None` like with `list`s and `dict`s\n",
        "\n",
        "---\n",
        "> `set` - an iterable that will only have unique values (no duplicates)\n",
        "   \\\n",
        "   sets are created with `set()` or `{}` and are **unsorted** but very fast - O(1) lookup\n",
        "   \\\n",
        "   instead of `.append` you use `.add`\n",
        "   \\\n",
        "   ```\n",
        "   s = {1,10,10,100,100,100,1000,1000,1000,1000}\n",
        "   print(s) # {1000, 1, 10, 100}\n",
        "   ```\n",
        "\n",
        "---\n",
        "> `zip` - combines multiple iterables (e.g. `list`, `array`, `set`, `tuple`, `str`)\n",
        "  \\\n",
        "  ```\n",
        "  chars = 'abcde'\n",
        "  nums = range(1,10)\n",
        "  zipped = zip(s, range(10))\n",
        "  print(zipped)       # <zip object at 0x7ff505aabc88>\n",
        "  print(list(zipped)) # [('a', 0), ('b', 1), ('c', 2), ('d', 3), ('e', 4)]\n",
        "  print(*zipped)      # ('a', 0), ('b', 1), ('c', 2), ('d', 3), ('e', 4)\n",
        "  ```\n",
        "  iterables are truncated to the shortest one\n",
        "\n",
        "---\n",
        "> `*` - spread operator\n",
        "  \\\n",
        "  we have used `*` in methods as a **rest operator** to 'collect' other arguments (and `**` to collect other *keyword* arguments). \n",
        "  \\\n",
        "  here, the operator 'flattens' an iterable. This can be useful for passing to methods, for copying data, and can be combined with `zip` to adjust the format of 2 related lists.\n",
        "  \\\n",
        "  ```\n",
        "  unzipped = zip(*zipped) # unzip and get original 2 lists\n",
        "  chars, nums, *rest = unzipped # note *rest is recent in Python\n",
        "  print(f\"chars={chars}, nums={nums}, rest={rest}\")\n",
        "  ```\n",
        "\n",
        "---\n",
        "> `lambda` - anonymous functions\n",
        "  \\\n",
        "  when we're feeling lazy$^1$ and don't want to name a function explicitly to implement a \"simple\" expression. Actually it can be named, go figure 🤷‍♂️.\n",
        "  \\\n",
        "  ```\n",
        "  anon = lambda x: x*2\n",
        "  for i in range(10):\n",
        "    print(f\"{i} - {anon(2):>2g\"})\n",
        "  ```\n",
        "  we passed an annoymous function to the `sorted` function (as part of the `key` param) in order to sort by `tuple`'s 2nd value (index `1`). The `key` param expects a single arg so we access the 2nd value of the passed `tuple`. We (silently) return `(x[1],x[0])` but it would sufficent in this case to just return `x[1]` because we don't care about additional sorting.\n",
        "  \\\n",
        "  [read more on sorting by value](https://stackoverflow.com/questions/613183/how-do-i-sort-a-dictionary-by-value)\n",
        "\n",
        "\n",
        "1. efficient\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eVC5T5GrlEbm",
        "colab_type": "text"
      },
      "source": [
        "#### Plot long data using Seaborn"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TRRpvxnor2El",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# try changing style and context to see how it affects the plot\n",
        "sns.set(style=\"ticks\", context=\"talk\")\n",
        "\n",
        "# Define a palette to ensure that colors will be\n",
        "# shared across the facets\n",
        "palette = dict(zip(dots.coherence.unique(),\n",
        "                   sns.color_palette(\"rocket_r\", 6)))\n",
        "\n",
        "# Plot the lines on two facets\n",
        "sns.relplot(x=\"time\", y=\"firing_rate\",\n",
        "            hue=\"coherence\", size=\"choice\", col=\"align\",\n",
        "            size_order=[\"T1\", \"T2\"], palette=palette,\n",
        "            height=5, aspect=.75, facet_kws=dict(sharex=False),\n",
        "            kind=\"line\", legend=\"full\", data=dots)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "twbhAnd3toQ3",
        "colab_type": "text"
      },
      "source": [
        "## Analysing someone else's code\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "shnGsTGCX-oB",
        "colab_type": "text"
      },
      "source": [
        "- can you make sense of the code below?\n",
        "- what happens when you run it?\n",
        "- are there inconsistencies in the variable names, function calls, or comments?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IRRfpYiR4E0R",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import cmath\n",
        "import time\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "from numpy.random import uniform, gamma, seed\n",
        "from scipy.sparse import csr_matrix\n",
        "\n",
        "seed(86664)\n",
        "\n",
        "## Parameter values\n",
        "#   Time\n",
        "dt = 0.5  # ms\n",
        "timeEnd = 5000  # ms\n",
        "T = int(np.ceil(timeEnd / dt))  # Number of time_array points\n",
        "\n",
        "#   Neurons\n",
        "n = 1000\n",
        "inh_ratio = 0.2\n",
        "n_inh = int(n * inh_ratio)\n",
        "n_exc = int(n * (1 - inh_ratio))\n",
        "k_inh = np.append(np.ones(n_inh, dtype='bool'), np.zeros(n_exc, dtype='bool'))\n",
        "# uniform(size=n) < 0.2  # 20 % inhibitory\n",
        "k_exc = np.logical_not(k_inh)  # opposite -> could also use (k_inh == False)\n",
        "\n",
        "a = k_inh.choose(0.02, 0.1)  # inhibitory cells a=0.1, excitatory a=0.02\n",
        "b = 0.2\n",
        "c = -65\n",
        "d = k_inh.choose(8, 2)\n",
        "tau_s = 10  # decay of synapses (ms)\n",
        "tau_d = 500  # synaptic depression (ms)\n",
        "stp_u = 0.5  # STP parameter\n",
        "ErevExc = 0.0\n",
        "ErevInh = -75.0\n",
        "\n",
        "# Input\n",
        "Itime = np.array([100, timeEnd-500]) / dt\n",
        "n_in = 100  # number of synaptic inputs\n",
        "w_in = 0.2\n",
        "pConnection = 0.01  # connection probability 'spontaneous activity'\n",
        "seed(int(pConnection * 100))\n",
        "C = uniform(size=(n, n_in)) < pConnection\n",
        "W_in = C.choose(0, w_in)\n",
        "f_rate = 0.010  # ms^-1)\n",
        "prate = f_rate * dt\n",
        "\n",
        "# Recurrent Input\n",
        "w_avg = 0.002  # average recurrent weight\n",
        "pConnectionRecurrent = 0.2  # connection probability\n",
        "W = np.zeros((n, n))\n",
        "C = uniform(size=(n, n))\n",
        "idx = np.nonzero(C < pConnectionRecurrent)\n",
        "g_sc = 0.002\n",
        "g_sh = w_avg / g_sc\n",
        "W[idx] = gamma(g_sh, g_sc, size=idx[0].size)  # gamma distributed array (shape,scale,size)\n",
        "scaleEI = 2  # weights from I->E different\n",
        "W[np.ix_(k_exc, k_inh)] *= scaleEI\n",
        "W = csr_matrix(W)  # make row sparse\n",
        "\n",
        "## Memory\n",
        "time_array = np.arange(1, timeEnd + 1, dt)\n",
        "v = np.zeros((T, n))\n",
        "u = np.zeros((T, n))\n",
        "v[0] = -70  # resting potential\n",
        "# Leak channels\n",
        "EL = -70\n",
        "gL = 0.0001\n",
        "u[0] = -14  # steady state\n",
        "\n",
        "S_in = np.zeros(n_in)  # synaptic input from external source\n",
        "s = np.zeros(n)  # recurrent synapses\n",
        "E_in = np.zeros(n_in)\n",
        "H_in = np.ones(n_in)  # synaptic depression term\n",
        "H = np.ones(n)  # synaptic depression term\n",
        "prevSpike = -np.infty * np.ones(n_in)\n",
        "prevSpikeRec = -np.infty * np.ones(n)\n",
        "E = k_inh.choose(ErevExc, ErevInh)\n",
        "\n",
        "r = uniform(size=n_in)\n",
        "# v = (r<p).choose(b,a)\n",
        "\n",
        "for t in np.arange(T-1):\n",
        "    if t * dt % 1000 == 0:\n",
        "        print(f'{t * dt / 1000}  s simulated.')\n",
        "    if Itime[0] < t < Itime[1]:\n",
        "        I = 0                           # background current\n",
        "        # Get input spikes (each neuron receives ~10% of input spikes)\n",
        "        P = uniform(size=n_in) < prate  # background input\n",
        "        tmp = np.exp(dt * (prevSpike[P] - t) / tau_d)\n",
        "        H_in[P] = 1 - (1 + (stp_u - 1) * H_in[P]) * tmp\n",
        "        prevSpike[P] = t\n",
        "    else:\n",
        "        I = 2\n",
        "        P = 0\n",
        "\n",
        "    # Calculate input current (all excitatory input)\n",
        "    S_in = (1 - dt / tau_s) * S_in + P * H_in\n",
        "    I += W_in.dot(S_in * E_in) - (W_in.dot(S_in) * v[t])\n",
        "\n",
        "    # handle all neurons\n",
        "    fired = v[t] >= 35\n",
        "    tmp = np.exp(dt * (prevSpikeRec[fired] - t) / tau_d)\n",
        "    H[fired] = 1 - (1 + (stp_u - 1) * H[fired]) * tmp\n",
        "    prevSpikeRec[fired] = t\n",
        "\n",
        "    # recurrent input\n",
        "    s = (1 - dt / tau_s) * s + fired\n",
        "    Irec = W * s * (E - v[t])\n",
        "    I += Irec\n",
        "\n",
        "    # Update ODE\n",
        "    dv = (0.04 * v[t] + 5) * v[t] + 140 - u[t]\n",
        "    v[t + 1] = v[t] + (dv + I) * dt\n",
        "    du = a * (b * v[t] - u[t])\n",
        "    u[t + 1] = u[t] + du * dt\n",
        "    # Spike!\n",
        "    v[t][fired] = 35\n",
        "    v[t + 1][fired] = c\n",
        "    u[t + 1][fired] = u[t][fired] + d[fired]\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "M_a4l83iYdZq",
        "colab_type": "text"
      },
      "source": [
        "#### Tasks"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ukHE8H6mYhtZ",
        "colab_type": "text"
      },
      "source": [
        "##### **Task 1**\n",
        "Write the system of equations that govern the neurons' behaviour\n",
        "\n",
        "Extra questions:\n",
        "\n",
        "- what does `np.ix_` do?\n",
        "- what does `csr_matrix` do? Is the comment useful?\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "b5DeUKIUY7Qr",
        "colab_type": "text"
      },
      "source": [
        "**Space for solution**\n",
        "\n",
        "[EDIT ME]\n",
        "\n",
        "eq.1: $\\frac{dv}{dt} = ... \\tag{neuronal voltage}$\n",
        "\n",
        "..."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4ypvA18bZqxp",
        "colab_type": "text"
      },
      "source": [
        "##### **Task 2**: plotting\n",
        "1. plot the **membrane potential** over time for the excitatory \"pyramidal cell\" neurons (**PC**) in red and the inhibitory \"interneuron\" neurons (**IN**) in blue. \n",
        "1. use a **raster plot** to show the spikes times (x-axis) for each neuron (y-axis). \n",
        "\n",
        "  ![example](https://upload.wikimedia.org/wikipedia/en/5/58/Sample_raster_plot_from_Brian_neural_network_simulator.jpg)\n",
        "\n",
        "1. plot a **histogram** of firing rates for each population\n",
        "\n",
        "Extra questions:\n",
        "- are PCs always excitatory? If not, under what condition(s)?\n",
        "- are INs always inhibitory? If not, under what condition(s)?\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0Y_bbzlE5Mvk",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title **Task 2** Solution { run: \"auto\", vertical-output: true, display-mode: \"form\" }\n",
        "time_width = 530 #@param {type:\"slider\", min:100, max:1000, step:50}\n",
        "# Plot\n",
        "fig, ax = plt.subplots(3, 1, sharex=True, figsize=(4,8),\n",
        "                       gridspec_kw={'height_ratios':[5,10,7]})\n",
        "\n",
        "# voltage vs time\n",
        "ax[0].plot(time_array, v[:, k_exc], 'r', alpha=0.1, lw=0.1)\n",
        "ax[0].plot(time_array, v[:, k_inh], 'b', alpha=0.1, lw=0.1)\n",
        "ax[0].plot(time_array, np.mean(v[:, k_exc],axis=1), 'r',  label='PC')\n",
        "ax[0].plot(time_array, np.mean(v[:, k_inh],axis=1), 'b',  label='IN')\n",
        "\n",
        "\n",
        "# get time (rows) and neuron index (columns) of spikes from v\n",
        "tspk, nspk = np.nonzero(v >= 35)\n",
        "idx_i = np.in1d(nspk, np.nonzero(k_inh)[0])  # inh spikes\n",
        "idx_e = np.logical_not(idx_i)                # exc spikes\n",
        "\n",
        "# nrn index vs time\n",
        "ax[1].plot(tspk[idx_e] * dt, nspk[idx_e], 'r.', label='Exc.', markersize=1)\n",
        "ax[1].plot(tspk[idx_i] * dt, nspk[idx_i], 'b.', label='Inh.', markersize=1)\n",
        "\n",
        "# plot histogram of firing rate\n",
        "bins = int(timeEnd / time_width)\n",
        "# plot histogram\n",
        "for i in range(0, bins - 1):\n",
        "    spikes_in_bin_indices = (tspk > time_width * i) & (tspk < time_width * (i + 1))\n",
        "    if np.any(spikes_in_bin_indices):\n",
        "        spikes_in_bin = tspk[spikes_in_bin_indices]\n",
        "        inh_spikes_in_bin = tspk[(idx_i) & (spikes_in_bin_indices)]\n",
        "        exc_spikes_in_bin = tspk[(idx_e) & (spikes_in_bin_indices)]\n",
        "        exc = ax[2].bar((time_width * i), len(exc_spikes_in_bin) / n_exc, \n",
        "                width=0.5*time_width,align='edge',color='r')\n",
        "        inh = ax[2].bar((time_width * i)+0.5*time_width, len(inh_spikes_in_bin) / n_inh, \n",
        "                width=0.5*time_width,align='edge',color='b')\n",
        "        avg = ax[2].bar((time_width * i), 0.0, bottom=len(spikes_in_bin) / n,\n",
        "                width=time_width, align='edge',color='k', alpha=1,edgecolor='k')\n",
        "    else:\n",
        "        ax[2].bar(time_width * i, 0, width=time_width)\n",
        "\n",
        "# annotate\n",
        "ax[0].set_ylabel('Membrane voltage \\n(mV)')\n",
        "ax[1].set_ylabel('Neuron index \\n(#)')\n",
        "ax[2].set_ylabel('Firing rate \\n(Hz)')\n",
        "\n",
        "ax[-1].set_xlim(0, timeEnd)\n",
        "ax[-1].set_xlabel('Time (ms)')\n",
        "ax[0].legend(loc='upper left', bbox_to_anchor=(1,1), frameon=False)\n",
        "ax[2].legend([avg], ['$\\mu (PC + IN)$'], loc='best', frameon=False)\n",
        "\n",
        "# clean\n",
        "sns.despine() # remove top + right borders"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wuV9TFmQg0kQ",
        "colab_type": "text"
      },
      "source": [
        "##### **Task 3**: plotting with Seaborn\n",
        "\n",
        "1. create a *long form* `DataFrame` with \"Neuron Index\" \"Time\" and \"Type\" columns (\"Time\" is the spike time) and each row is an observation.\n",
        "\n",
        "1. plot a raster plot\n",
        "\n",
        "1. plot a raster plot with *kernel density estimates* (run, but don't look at, the solution to see what I mean)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_DosMAikhMEs",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import seaborn as sns"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9x2GthFbf0k1",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title **Task 3.1** Solution { run: \"auto\", vertical-output: true, display-mode: \"form\" }\n",
        "import pandas as pd\n",
        "spikes = pd.DataFrame({\"Neuron Index\":nspk, \"Time\":tspk*dt})\n",
        "spikes[\"Type\"] = spikes[\"Time\"]\n",
        "spikes.loc[idx_i,\"Type\"] = \"IN\"\n",
        "spikes.loc[idx_e,\"Type\"] = \"PC\"\n",
        "spikes"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "eSHR-zGgaTSd",
        "colab": {}
      },
      "source": [
        "#@title **Task 3.2** Solution { run: \"auto\", vertical-output: true, display-mode: \"form\" }\n",
        "sns.scatterplot(x=\"Time\", y=\"Neuron Index\",\n",
        "             hue=\"Type\",  palette=\"Dark2\", s=2, ec='None',\n",
        "             data=spikes)\n",
        "plt.xlim(0, timeEnd)\n",
        "plt.ylim(0, n)\n",
        "plt.title(\"Raster plot\")\n",
        "plt.legend(loc=(1,0.5), frameon=False)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IJILvQnri0EK",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title **Task 3.3** Solution { run: \"auto\", vertical-output: true, display-mode: \"form\" }\n",
        "g = sns.jointplot(x=\"Time\", y=\"Neuron Index\", s=0.2, data=spikes)\n",
        "g.ax_joint.axhline(y=200, lw=2, c='k')\n",
        "g.ax_joint.set_xlim(0, timeEnd)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jhnvIARJ9XgD",
        "colab_type": "text"
      },
      "source": [
        "# Homework/Projects\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "364OraHiz1om",
        "colab_type": "text"
      },
      "source": [
        "> ## **Homework 3** Implement a balanced spiking neural network with 1000 neurons, 20% interneurons, 2% general connectivity (including recurrent). Implement $AMPA$ and $GABA_A$ synapses, and a *low* background current. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cNsMUWK_z6Yc",
        "colab_type": "text"
      },
      "source": [
        "Build from scratch or use [Brian2](https://brian2.readthedocs.io/en/stable/) (or another simulator!)\n",
        "\n",
        "1. How do you scale the synaptics conductances with the number of neurons? \n",
        "1. Change the $AMPA$ and $GABA_A$ time constants, is the scaling rule still the same?\n",
        "1. Implement $NMDA$ and/or $GABA_B$ synapses. How has the scaling rule changed now?\n",
        "1. Restrict interneuron to only connect with neurons within 10 mm. Restrict pyramidal neurons to only connect with neurons at least 50 mm away. *hint: you may need to set up a virtual \"grid\" (2D/3D) on which to place neurons*.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Yzzz4H86gIyc",
        "colab_type": "text"
      },
      "source": [
        "> ## **Morphology**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "23YlP_qzgrIh",
        "colab_type": "text"
      },
      "source": [
        "A hippocampal pyramidal neuron can be considered a two-layer neural network. \n",
        "\\\n",
        "Build a simple dendritic structure and see if you can repeat the results found in [Pyramidal neuron as two-layer neural network. Poirazi P, Brannon T, Mel BW. Neuron. 2003 Mar 27;37(6):989-99.](https://doi.org/10.1016/S0896-6273(03)00149-1).\n",
        "\\\n",
        "  *Use [NEURON](https://neuron.yale.edu/neuron/) or [Brian2](https://brian2.readthedocs.io/en/stable/) to remain sane.*\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "L0whNOumgfMI",
        "colab_type": "text"
      },
      "source": [
        "Implement a spiking neural network in PyTorch that repeats this experiment: [Towards deep learning with segregated dendrites](https://elifesciences.org/articles/22901)\n",
        "![Towards deep learning with segregated dendrites](https://iiif.elifesciences.org/lax:22901%2Felife-22901-fig3-v1.tif/full/1500,/0/default.jpg)\n",
        "\\\n",
        "*see [BindsNET](https://groups.cs.umass.edu/binds/bindsnet-spiking-neural-networks-in-pytorch/)*"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lWfg0ElYpXUj",
        "colab_type": "text"
      },
      "source": [
        "___\n",
        "# Appendix\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C8eXA81-AvOd",
        "colab_type": "text"
      },
      "source": [
        "## Eigen$[x]$\n",
        "where $ [x]  \\in  \\{\\rm{vector, value, faces,...}\\}$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DshsdNisGOTv",
        "colab_type": "text"
      },
      "source": [
        "### **What** are eigenvectors and eigenvalues?\n",
        "> One of the most widely used kinds of matrix decomposition is called eigendecomposition, in which we decompose a matrix into a set of eigenvectors and eigenvalues. *Page 42, Deep Learning, 2016.*\n",
        "\n",
        "Eigendecomposition of a matrix is a type of decomposition that involves decomposing a **square** ($n \\times n$) matrix into a **set** of *eigenvectors* and *eigenvalues*.\n",
        "\n",
        "A vector $\\rm{x}$ is an eigenvector of a matrix $A$ if it satisfies the following equation.\n",
        "\n",
        "$A\\rm{x} = \\lambda \\rm{x} \\tag{E.1}$\n",
        "\n",
        "This is called the eigenvalue equation, where $A$ is the parent square matrix that we are decomposing, $x$ is the eigenvector of the matrix, and $\\lambda$ is the lowercase Greek letter lambda and represents the eigenvalue scalar.\n",
        "\n",
        "This can also be rearranged into\n",
        "\n",
        "$(A-\\lambda I)\\rm{x} = 0 \\tag{E.2}$\n",
        "\n",
        "Where $I$ is the square identity matrix ($1$s in its main diagonal and $0$s in every other entry) with the same dimensions as $A$. \n",
        "\n",
        "The corresponding linear map is the identity map. For any $n \\times n$ matrix $A$, we have $AI = IA = A$. The inverse $A^{-1}$ of $A$ is the unique matrix such that $AA^{-1} = A^{-1}A = I$.\n",
        "\n",
        "$$\\mathbf{I} = \\left.\\left( \n",
        "                  \\vphantom{\\begin{array}{c}1\\\\1\\\\1\\\\1\\\\1\\end{array}}\n",
        "                  \\smash{\\underbrace{\n",
        "                      \\begin{array}{ccccc}\n",
        "                             1&0&0&\\cdots &0\\\\\n",
        "                             0&1&0&\\cdots &0\\\\\n",
        "                             0&0&1&\\cdots &0\\\\\n",
        "                             \\vdots&&&\\ddots&\\\\\n",
        "                             0&0&0&\\cdots &1\n",
        "                      \\end{array}\n",
        "                      }_{n\\text{ columns}}}\n",
        "              \\right)\\right\\}\n",
        "              \\,n\\text{ rows}\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "as77YuNwHO0e",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "### But **why** should I care?\n",
        "> Almost all vectors change direction, when they are multiplied by A. Certain exceptional vectors x are in the same direction as Ax. Those are the “eigenvectors”. Multiply an eigenvector by A, and the vector Ax is the number lambda times the original x. […] The eigenvalue lambda (λ) tells whether the special vector x is stretched or shrunk or reversed or left unchanged – when it is multiplied by A.\n",
        "*— Page 289, Introduction to Linear Algebra, Fifth Edition, 2016.*\n",
        "\n",
        "![](https://upload.wikimedia.org/wikipedia/commons/thumb/5/58/Eigenvalue_equation.svg/2560px-Eigenvalue_equation.svg.png)\n",
        "> Matrix A acts by stretching the vector x, not changing its direction, so x is an eigenvector of A.\n",
        "\n",
        "A matrix could have one eigenvector and eigenvalue for each dimension of the parent matrix. Not all square matrices can be decomposed into eigenvectors and eigenvalues, and some can only be decomposed in a way that requires complex numbers. \n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5N-t6_pEHNrQ",
        "colab_type": "text"
      },
      "source": [
        "### Okay, show me **how**\n",
        "In order to determine the eigenvectors of a matrix, you must first determine the eigenvalues. \n",
        "\n",
        "Non-trivial solutions exist only if the matrix ($A-\\lambda I$) is singular which means $\\det(A-\\lambda I) = 0$. Where $\\det(M)$ or $|M|$ is the [determinant](https://en.wikipedia.org/wiki/Determinant) of the matrix.\n",
        "\n",
        "Therefore eigenvalues of $A$ are roots of the [characteristic polynomial](https://en.wikipedia.org/wiki/Characteristic_polynomial)\n",
        "\n",
        "$p(\\lambda) = \\det(A-\\lambda I) \\tag{E.3}$\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DGU7fv0eblZz",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "#### To find $\\lambda$ such that $\\det(A-\\lambda I) = 0$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xsUwUROHbm9W",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "   if $A = \\begin{bmatrix}\n",
        "                             a_{11}&a_{12}\\\\\n",
        "                             a_{21}&a_{22}\n",
        "                      \\end{bmatrix}$\n",
        "   \n",
        "   i.e. $A$ is a $2 \\times 2$ matrix, \n",
        "   \n",
        "   then\n",
        "\n",
        "   $\\lambda I = \\begin{bmatrix}\n",
        "                             \\lambda&0\\\\\n",
        "                             0&\\lambda\n",
        "                      \\end{bmatrix}$ \n",
        "\n",
        "   therefore \n",
        "\n",
        "   $\\begin{align} A &- \\lambda I & &= \\\\\n",
        "   \\begin{bmatrix}\n",
        "                             a_{11}&a_{12}\\\\\n",
        "                             a_{21}&a_{22}\n",
        "                      \\end{bmatrix} &- \n",
        "                      \\begin{bmatrix}\n",
        "                             \\lambda&0\\\\\n",
        "                             0&\\lambda\n",
        "                      \\end{bmatrix} & &= \n",
        "                      \\begin{bmatrix}\n",
        "                             a_{11}-\\lambda&a_{12}\\\\\n",
        "                             a_{21}&a_{22}-\\lambda\n",
        "                      \\end{bmatrix}\n",
        "                      \\end{align}$\n",
        "\n",
        "   the determinant of which is \n",
        "\n",
        "   $\\begin{align}\\begin{vmatrix}\n",
        "                             a_{11}-\\lambda&a_{12}\\\\\n",
        "                             a_{21}&a_{22}-\\lambda\n",
        "                      \\end{vmatrix} &= (a_{11}-\\lambda) \\cdot  (a_{22}-\\lambda)- a_{12} \\cdot a_{21}\\\\\n",
        "   &= \\lambda^2 - \\lambda (a_{11} + a_{22}) + a_{11} \\cdot a_{22} - a_{12} \\cdot a_{21}\n",
        "   \\end{align}$\n",
        "\n",
        "   finding the roots of the determinant (characteristic polynomial) gives us the eigenvalues. In this case for a $2 \\times 2 $ matrix, there are 2 roots $\\lambda_1$ and $\\lambda_2$. \n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CA6aiZzSKe0V",
        "colab_type": "text"
      },
      "source": [
        "Once the eigenvalues of a matrix (A) have been found, we can find the eigenvectors by Gaussian Elimination.\n",
        "\n",
        "- STEP 1: For each eigenvalue λ, we have\n",
        "$(A − λI)\\rm{x} = 0$,\n",
        "where $\\rm{x}$ is the eigenvector associated with eigenvalue λ.\n",
        "- STEP 2: Find $\\rm{x}$ by [Gaussian elimination](https://en.wikipedia.org/wiki/Gaussian_elimination). That is, convert the augmented matrix\n",
        "$(A − λI: 0)$\n",
        "to [row echelon form](https://en.wikipedia.org/wiki/Row_echelon_form), and solve the resulting linear system by back substitution.\n",
        "\n",
        "see [this](https://www.scss.tcd.ie/Rozenn.Dahyot/CS1BA1/SolutionEigen.pdf) resource for more information.\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FSfMaCsnkvWD",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "Q:\n",
        "Determine the **eigenvalues** and **eigenvectors** for \\begin{bmatrix}\n",
        "                             2&-3&0\\\\\n",
        "                             2&-5&0\\\\\n",
        "                             0&0&3\n",
        "                      \\end{bmatrix}\n",
        "\n",
        "*hint*\n",
        "\n",
        "1. There are 3 eigenvalues.\n",
        "1. The determinant for a $3 \\times 3$ matrix is\n",
        "\n",
        "   ![](https://wikimedia.org/api/rest_v1/media/math/render/svg/14f2f2a449d6d152ee71261e47551aa0a31c801e)\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ghfaE05VkybJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title Solution { display-mode: \"form\" }\n",
        "print(f\"Eigenvalues (λ) are [3,1,-4]\")\n",
        "print(\"\"\"Eigenvectors are \n",
        "λ= 3, x=[ 0  0 1].T\n",
        "λ= 1, x=[ 3  1 2].T\n",
        "λ=-4, x=[0.5 1 0].T\n",
        "\"\"\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kRiiQcWJSjMF",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "##### *Determinant*\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PYXJfU2Ofgdp",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "The determinant encodes certain properties of the linear transformation described by the matrix. \n",
        "\n",
        "Geometrically, it can be viewed as the area/volume scaling factor of the linear transformation described by the matrix\n",
        "\n",
        "![](https://upload.wikimedia.org/wikipedia/commons/thumb/a/ad/Area_parallellogram_as_determinant.svg/440px-Area_parallellogram_as_determinant.svg.png)\n",
        "\n",
        "**$2 \\times 2$ matrix**\n",
        "\n",
        "![](https://wikimedia.org/api/rest_v1/media/math/render/svg/84a06831f59ecb179385f947303f3ef3b4a6a509)\n",
        "\n",
        "\n",
        "**$3 \\times 3$ matrix**\n",
        "\n",
        "![](https://wikimedia.org/api/rest_v1/media/math/render/svg/1e2b66aee07eda072fc4996c89915b0554220226)\n",
        "\n",
        "**$n \\times n$ matrix**\n",
        "\n",
        "The determinant of a matrix of arbitrary size can be defined by the [Leibniz formula](https://en.wikipedia.org/wiki/Leibniz_formula_for_determinants) or the [Laplace formula](https://en.wikipedia.org/wiki/Laplace_expansion).\n",
        "\n",
        "![](https://wikimedia.org/api/rest_v1/media/math/render/svg/6c70973175f00a8a18b80784dbeff28808093898)\n",
        "\n",
        "[*read more*](https://en.wikipedia.org/wiki/Determinant#n_%C3%97_n_matrices)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RyC2thOWO0P5",
        "colab_type": "text"
      },
      "source": [
        "### Python to the rescue\n",
        "The function [`numpy.linalg.eig`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html) computes eigenvalues and eigenvectors of a square matrix. `scipy.linalg.eig` **also** works.\n",
        "\n",
        "note: the returned vectors are *normalised* (<=1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uCbf9eOzPV6-",
        "colab_type": "text"
      },
      "source": [
        "Let's consider a simple example with a diagonal matrix:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MK2Rf3ZfPrBr",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "from numpy.linalg import eig\n",
        "import scipy.linalg as la\n",
        "# define matrix\n",
        "A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])\n",
        "print(A)\n",
        "# calculate eigendecomposition\n",
        "values, vectors = eig(A)\n",
        "print(values)\n",
        "print(vectors)\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pji7mgCaQC-y",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "A = np.array([[1, -3, 3],\n",
        "              [3, -5, 3],\n",
        "              [6, -6, 4]])\n",
        "values, vectors = eig(A)\n",
        "print(values)\n",
        "print(vectors)\n",
        "print(np.allclose(A@vectors[:,0], values[0]*vectors[:,0]))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xZH2g33NpnxT",
        "colab_type": "text"
      },
      "source": [
        "## Style"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eXCkHD1D-Tus",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "\"\"\"\n",
        "Files should start with a multi-string comment to describe what it does and any \n",
        "appropriate licensing/authorship information. Try to limit lines in your code to a fixed width (e.g. 80) otherwise it starts becoming awkward for others to read the comment or code.\n",
        "\n",
        "This file is a Notebook that makes programming in Python (and other supported \n",
        "languages) more interactive and sometimes easier to develop.\n",
        "\n",
        "Authors: \n",
        "Christopher Brian Currin (https://www.github.com/chriscurrin)\n",
        "\"\"\""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zH5X3bov9kUO",
        "colab_type": "text"
      },
      "source": [
        "### Object-Orientated Programming (OOP)\n",
        "\n",
        "Some consider OOP to be an oops because it is easy to be abused (everything becomes a class and is inherited, poluting a bunch of empty parent abstract classed) and can take longer to understand. However, once the code is understood, it *can* be easier to extend the programme. \n",
        "The OOP style is generally preferred for developing *libraries* as they encapsulate a lot of complexity from the user."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OfwLH3oXjGfq",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "\n",
        "class ObjectOrientatedObject(object):\n",
        "  \"\"\"This parent object should be inherited in objects that wish to have a \n",
        "  'foo'method.\n",
        "  \"\"\"\n",
        "  def __init__(self, passed_variable, *args, **kwargs):\n",
        "    \"\"\"To create an object (this line can be tested by the `doctest` module):\n",
        "    >>> ooo = ObjectOrientatedObject(\"ahh\")\n",
        "    ooo.object_variable == \"ahh\"\n",
        "    \n",
        "    Including *args and **kwargs are useful for class methods so that arguments\n",
        "    can be passed to parent and child methods more robustly.\n",
        "    \"\"\"\n",
        "    self.object_variable = passed_variable\n",
        "    self.__init()\n",
        "    \n",
        "  def __init(self):\n",
        "    \"\"\"\n",
        "    Method names starting with '__' are not autocompleted by created (also known\n",
        "    as 'instantiated') objects, and thus are useful names for private methods\n",
        "    that are only to be used within the class by other class methods.\n",
        "    This method initialises some objects to remove some logic from the main\n",
        "    __init__ method, making it clearer what relies on passed variables. \n",
        "    \"\"\"\n",
        "    self.class_list = []\n",
        "      \n",
        "  def foo(self):\n",
        "    # generally we don't want to print in an object.\n",
        "    print(\"foo has been called\")\n",
        "    \n",
        "  def bar(self):\n",
        "    \"\"\"Provide a stub for methods that *must* be implemented by child classes.\n",
        "    \"\"\"\n",
        "    raise NotImplementedError(\"implementation for future classes or versions\")\n",
        "    \n",
        "  def zoo(self):\n",
        "    \"\"\"Provide a stub for methods that *should* be implemented by child classes.\n",
        "    \"\"\"\n",
        "    pass\n",
        "    \n",
        "  def manipulate_data(self, more_str: str):\n",
        "    \"\"\"Object orientated programming often hides data manipulation in an object's\n",
        "    method such that it's not always clear what occurs when running the object.\n",
        "    \"\"\"\n",
        "    self.object_variable = more_str*2 + self.object_variable\n",
        "    self.class_list.append(more_str)\n",
        "    \n",
        "  def init_list(self, other_list: list) -> list:\n",
        "    \"\"\" The `other_list` is assigned to this classes `class_list` property.\n",
        "    However, the objects become linked such that changes to one affects the \n",
        "    other. To prevent this, we can use the `copy` module.\n",
        "    \"\"\"\n",
        "    self.class_list = other_list\n",
        "  \n",
        "# in notebooks, it can be useful to have some small test to check the class \n",
        "# works as expected.\n",
        "ooo = ObjectOrientatedObject(\"AHH!\")\n",
        "# is this clear?\n",
        "ooo.foo()\n",
        "\n",
        "print(ooo.object_variable)\n",
        "ooo.manipulate_data(\"OOO\")\n",
        "print(ooo.object_variable)\n",
        "\n",
        "try:\n",
        "  ooo.bar()\n",
        "except NotImplementedError as err:\n",
        "  print(err)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wsEn7DjvJ8Lv",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class BetterOOP(ObjectOrientatedObject):\n",
        "\n",
        "  def __init__(self, *args, **kwargs):\n",
        "    super().__init__(*args,**kwargs)\n",
        "    \n",
        "  def foo(self) -> str:\n",
        "    return \"better foo has been called\"\n",
        "  \n",
        "  def bar(self):\n",
        "    return \"implemented function!\"\n",
        "  \n",
        "  def manipulate_data(self, *args) -> bool:\n",
        "    old_object = self.object_variable\n",
        "    super().manipulate_data(*args)\n",
        "    new_object = self.object_variable\n",
        "    return old_object\n",
        "  \n",
        "boop = BetterOOP(\"required by parent\")\n",
        "print(boop.foo())"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eFr7kXmQ9U2i",
        "colab_type": "text"
      },
      "source": [
        "## Version Control\n",
        "![](https://madusudanan.com/images/xkcd-vc.jpg)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "go4ZwjkrZXa-",
        "colab_type": "text"
      },
      "source": [
        "1. Go to https://github.com/ or https://bitbucket.org/ or https://gitlab.com/ and sign up for an account.\n",
        "\n",
        "1. create a new repository (\"repo\")\n",
        "\n",
        "1. Go to https://git-scm.com/downloads and follow the instructions to download git\n",
        "\n",
        "1. Open a Terminal (Windows 10: download [Windows Terminal Preview](https://www.microsoft.com/en-us/p/windows-terminal-preview/9n0dx20hk701))\n",
        "\n",
        "1. \"clone\" (download and create a local version) your repo from (2) using `git clone <repo url>`\n",
        "  > what is the difference between \n",
        "  \\\n",
        "  `git clone https://github.com/ChrisCurrin/comp_neuro_count_countries.git` \n",
        "  \\\n",
        "  and \n",
        "  \\\n",
        "  `git clone git@github.com:ChrisCurrin/comp_neuro_count_countries.git`"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GKPQ6Lizb0-5",
        "colab_type": "text"
      },
      "source": [
        "Here are some git commands for creating a repo and commit a file.\n",
        "\n",
        "> Note: the commands start with a bang (!) because we are in a notebook and git is a command prompt tool and not Python. Lines start with ! are passed to the underlying terminal/console for execution. In your own terminal, ignore the '!'s."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1W9Zxlz5dF4m",
        "colab_type": "text"
      },
      "source": [
        "If you've down `git clone ...` you don't need to `git init`"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wwLxyWWb9XF9",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!git init"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "y5vshgEL9XUB",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!git status"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "c87sFpwrXwrp",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!cat > Readme.md"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_FLHUMt8aZnE",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!echo \"Readme\" >> Readme.md\n",
        "!git add Readme.md\n",
        "!git status"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "z4--r-D3X9Pu",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!cat Readme.md"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GVtsuo81ZOVG",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!git commit -m \"commiting files which are staged\""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KUsmtHeRZXhB",
        "colab_type": "text"
      },
      "source": [
        "Pushing and pulling won't work until you've set up an \"upstream\" repository (repo) on a website like https://github.com\n",
        "\n",
        "- `git push` = upload local files and changes to another location (e.g. a cloud service)\n",
        "- `git pull` = download (or `fetch`) changes from the cloud to your local machine and `merge` the changes (if there are conflicts where the same line has been changed, you will need to fix these before a merge is completed)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jItPSHnRaTwA",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!echo \"this is a version controlled file \" >> Readme.md"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QBnjRQfoakjc",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!git status"
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}