tutorial1.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X3ukvoQEumJu"
      },
      "source": [
        "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/jmourabarbosa/lowrank_tutorial/blob/main/tutorial1.ipynb)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5YTGVCJaumJw"
      },
      "source": [
        "### uncomment this if running in google colab"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 82,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "bxW0AuECumJw",
        "outputId": "7253f0b7-03c5-4ee3-fef2-45ecd1d6a7f0"
      },
      "outputs": [],
      "source": [
        "!git clone https://github.com/jmourabarbosa/lowrank_tutorial/\n",
        "import sys\n",
        "sys.path.append('lowrank_tutorial/')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 83,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "64FGBI_yMcUU",
        "outputId": "e70f41d2-f2ca-4b53-f044-56797ab7d585"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<module 'lowrank_helpers' from '/Users/jbarbosa/Dropbox/Neuro/QBio_school/lowrank_tutorial/lowrank_helpers.py'>"
            ]
          },
          "execution_count": 83,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pylab as plt\n",
        "import figures.lowrank_helpers as lrh\n",
        "import imp\n",
        "imp.reload(lrh)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "R4DMYIaMdlKX"
      },
      "source": [
        "\n",
        "# Introduction to (low-rank) recurrent neural networks\n",
        "\n",
        "### Goals of this tutorial:\n",
        "\n",
        "* Introduction to recurrent neural networks (RNN)\n",
        "* basics of low-rank RNN theory (linear)\n",
        "* basics of low-rank RNN theory (non-linear)\n",
        "\n",
        "Note: This tutorial is heavily based on the paper and code provided in [(Mastrogiuseppe & Ostojic, 2019)](https://github.com/fmastrogiuseppe/LowRank). For more advance concepts, please refer to that paper.\n",
        "\n",
        "In a recurrent neural network (RNN), each neuron has a time-dependent membrane voltage $x_i(t)$ and produces a firing rate through the transfer function $\\phi(x_i(t))$.\n",
        "\n",
        "Here, we will use either $\\phi(x) = x$ or $\\phi = tanh(x)$. To declutter the text, we will drop time and use vector notation.\n",
        "\n",
        "As the name suggests, neurons in a RNN are recurrently connected so each pair of neurons ${i, j}$ has a synaptic connexion of strength $J_{ij}$.\n",
        "\n",
        "For the purpose if this tutorial, the recurrent connectivity will be rank 1: $J = \\frac{1}{N} mn^T$.\n",
        "\n",
        "In adition to recurrent inputs, each neuron might be driven by external inputs. For simplicity, we will consider the case of one input $u(t)$ that is fed to each neuron through the weight vector $I$.\n",
        "\n",
        "Finally, each neuron integrates all its inputs into its voltage via the differential equation:\n",
        "\n",
        "$$\\dot{x} = -x + \\frac{1}{N}mn^T\\phi(x) + uI$$\n",
        "\n",
        "Assumming that $m \\perp I$, the activity of this network will be confined to a 2D plane defined by $m$ and $I$:\n",
        "\n",
        "$$x = \\kappa m + vI \\; \\; \\; \\small (1)$$\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# 1. Linear low-rank RNN"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nJAGlTQzoxeG"
      },
      "source": [
        "***Exercise 1.1.***\n",
        "\n",
        "An important insight of the low-rank RNN theory is that the alignment between the vectors definning the RNN (in this case $m,n,I$) fully determine the dynamics performed by the network.\n",
        "\n",
        "\n",
        "<p align=\"center\">\n",
        "<img src=\"https://github.com/jmourabarbosa/lowrank_tutorial/blob/main/figures/mastrogiuseppe.png?raw=1\" alt=\"drawing\" width=\"600\"/>\n",
        "</p>\n",
        "\n",
        "\n",
        "\n",
        "We will gradually develop the intuition (and theory) for why this is the case.\n",
        "\n",
        "* Which scenario are we for the vectors below m, n, IA?\n",
        "\n",
        "* How about IB?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 95,
      "metadata": {
        "id": "Xbrq1-7GumJx"
      },
      "outputs": [],
      "source": [
        "m,n,IA,IB = lrh.m, lrh.n, lrh.IA, lrh.IB\n",
        "\n",
        "N=len(m)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XaRw7ka5umJy"
      },
      "source": [
        "***Exercise 1.2***\n",
        "\n",
        "* Use the following code to run a *linear* rank 1 network with $J = \\frac{1}{N}mn^T$ that recieves inputs from IA\n",
        "\n",
        "* Plot the currents of some neurons.\n",
        "\n",
        "* Plot the mean activity.\n",
        "\n",
        "\n",
        "<p align=\"center\">\n",
        "<img src=\"https://github.com/jmourabarbosa/lowrank_tutorial/blob/main/figures/1.2.png?raw=1\" alt=\"drawing\" width=\"600\"/>\n",
        "</p>\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 96,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 449
        },
        "id": "As3hQBhpZqus",
        "outputId": "72a4703a-2d5b-4975-ac55-3d715fe8022b"
      },
      "outputs": [],
      "source": [
        "def phi(x):\n",
        "  return x\n",
        "\n",
        "dt = 0.1\n",
        "\n",
        "time = np.linspace(0,30,int(30//dt))\n",
        "\n",
        "J = (...)\n",
        "u = np.zeros_like(time)\n",
        "\n",
        "# inputs are on during this period, off otherwise\n",
        "u[50:150] = 1\n",
        "\n",
        "# two trials, one for each input\n",
        "x = np.zeros(...)\n",
        "X = np.zeros(...)\n",
        "\n",
        "for ti, _ in enumerate(time):\n",
        "  x[0] = x[0] + dt*(...)\n",
        "  x[1] = x[1] + dt*(...)\n",
        "  X[:,ti,:] = x\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "plt.plot((...),\"gray\")\n",
        "plt.plot((...),\"r\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "a5GQCwEisaCT"
      },
      "source": [
        "***Exercise 1.3***\n",
        "\n",
        "Instead of looking at each neuron individually, it is more informative to do dimensionality reduction.\n",
        "\n",
        "One way of doing this is using PCA. You can try, but given our theory we know that the actvity of x is going to be in m-I plane (equation 1, above).\n",
        "\n",
        "* Check this fact by projecting the activity on m and I, separatedly, then together.\n",
        "\n",
        "* Can you guess how many PCs would PCA return?\n",
        "\n",
        "\n",
        "<p align=\"center\">\n",
        "<img src=\"https://github.com/jmourabarbosa/lowrank_tutorial/blob/main/figures/1.3a.png?raw=1\" alt=\"drawing\" width=\"600\"/>\n",
        "</p>\n",
        "\n",
        "<p align=\"center\">\n",
        "<img src=\"https://github.com/jmourabarbosa/lowrank_tutorial/blob/main/figures/1.3b.png?raw=1\" alt=\"drawing\" width=\"600\"/>\n",
        "</p>\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 863
        },
        "id": "7ZGr4aJ9XNX9",
        "outputId": "03651e0f-d7e7-4253-a6c6-fadcb4b57601"
      },
      "outputs": [],
      "source": [
        "plt.plot(time,X[0] @ (...) / np.linalg.norm(...)**2)\n",
        "plt.plot(time,X[1] @ (...) / np.linalg.norm(...)**2)\n",
        "plt.xlabel(\"Time\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nHFlTLuDumJy"
      },
      "source": [
        "***Exercise 1.5***\n",
        "\n",
        "* Do the same as above, but make the network receive inputs through IB.\n",
        "* Plot the the activity of some neurons\n",
        "* Plot the mean activity\n",
        "* Plot the m-I plane\n",
        "* What would be the output of PCA, for this case?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# 2. Solving for kappa (Linear case)\n",
        "\n",
        "We have now checked that $x = \\kappa m + v I \\implies \\kappa = \\frac{m}{|m|^2} x $.\n",
        "\n",
        "We can use this fact (and that $m \\perp I$ ) to get an equation for $\\dot \\kappa$ by projecting $\\dot x$ on $\\frac{m}{|m|^2}$. Similarly for $v$:\n",
        "\n",
        "$$ \\dot \\kappa = -k + \\frac{1}{N}n^T\\phi(\\kappa m + v I) \\; \\; \\; \\small (2)$$\n",
        "$$ \\dot v = -v + u $$\n",
        "\n",
        "In the case where $\\phi(x) = x$ and using simple linear algebra we can simply this futher:\n",
        "\n",
        "$$ \\dot \\kappa = -k + \\sigma_{nm}\\kappa + \\sigma_{nI}v \\; \\; \\; \\small (3)$$\n",
        "$$ \\dot v = -v + u$$\n",
        "\n",
        "\n",
        "with $\\sigma_{ab} = a^Tb/N$.\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "***Exercise 2.1.***\n",
        "\n",
        "* Calculate $\\sigma_{mn}$ and $\\sigma_{nI}$ (for both inputs)\n",
        "* Use these overlaps to integrate $\\dot \\kappa$ and $\\dot v$ (tip: reuse the above code to integrate $\\dot x$)\n",
        "* plot on the same plot $\\dot \\kappa$ and $\\dot x$ using $\\phi(x) = x$ and $\\phi(x) = tanh(x)$.\n",
        "* Do they match?\n",
        "\n",
        "<p align=\"center\">\n",
        "<img src=\"https://github.com/jmourabarbosa/lowrank_tutorial/blob/main/figures/2.1.png?raw=1\" alt=\"drawing\" width=\"600\"/>\n",
        "</p>\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# two trials, one for each input\n",
        "kappa = [0,0]\n",
        "v = [0,0]\n",
        "K = np.zeros((2,len(time)))\n",
        "\n",
        "for ti in range(len(time)):\n",
        "\n",
        "    # one input\n",
        "    s_sigma_mn, s_sigma_nI = \n",
        "    kappa_rec = \n",
        "    v[0] = \n",
        "    kappa[0] = kappa[0] + dt*(...)\n",
        "\n",
        "    # second input\n",
        "\n",
        "    K[:,ti] = kappa\n",
        "\n",
        "\n",
        "plt.plot(time,K[0],'k--',label='kappa')\n",
        "plt.plot(...)\n",
        "\n",
        "plt.plot(time,K[1],'r--')\n",
        "plt.plot(...)\n",
        "\n",
        "plt.legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# 3. non-linear low-rank RNN ##"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jIkx2W-rumJy"
      },
      "source": [
        "***Exercise 3***\n",
        "\n",
        "* Use the code provided above to run a *non-linear* network with the same m, n, IA, IB vectors, but with tanh as the transfer function.\n",
        "\n",
        "* Look at the activations and projections on m and I.\n",
        "\n",
        "* What changed?\n",
        "\n",
        "* Generate new set of network vectors in which the overlap between m and n is = 1.2 (the third case in 1.1)\n",
        "\n",
        "* how bad is the linear approximation? Why?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# code generate guassian vectors\n",
        "\n",
        "overlaps = ...\n",
        "m,n,IA,IB  = np.random.multivariate_normal([0,0,0,0], overlaps, size=N).T\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DpoxNv9RumJz"
      },
      "source": [
        "# 4. Solving for kappa (non-linear case)\n",
        "\n",
        "Now that we checked our intutions with the linear case, lets solve $\\dot \\kappa$ also for the non-linear case. \n",
        "\n",
        "There's a way to solve $\\kappa$ just using the overlaps (similar to equation 3; we will try that later), but let reduce the dimensionality of our dynamical system just by integrating equation 2.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8aBCjC27umJz"
      },
      "source": [
        "***Exercise 4.1***\n",
        "* reuse the code you used above for integreating $\\dot x$ and $\\dot \\kappa$ for the linear case to integrate equation 2.\n",
        "* compare what you got with the non-linear RNN.\n",
        "\n",
        "<p align=\"center\">\n",
        "<img src=\"https://github.com/jmourabarbosa/lowrank_tutorial/blob/main/figures/4.png?raw=1\" alt=\"drawing\" width=\"600\"/>\n",
        "</p>\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 88,
      "metadata": {
        "id": "WG6cdVHRumJz"
      },
      "outputs": [],
      "source": [
        "\n",
        "# two trials, one for each input\n",
        "k_ = np.zeros([2])\n",
        "K_ = np.zeros((2,len(time)))\n",
        "v = np.zeros([2])\n",
        "\n",
        "for ti, _ in enumerate(time):\n",
        "  # one input\n",
        "  v[0] = v[0] + (...)\n",
        "  k_[0] = k_[0] + (...)\n",
        "\n",
        "  # second input\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 448
        },
        "id": "KkFigrwXumJz",
        "outputId": "5fd8e11f-8992-4132-85f5-bc534fcb5be4"
      },
      "outputs": [],
      "source": [
        "plt.plot(time,K_[0],'k--',label='kappa')\n",
        "plt.plot(time,K_[1],'r--')\n",
        "\n",
        "plt.plot(...,label='activity proj m')\n",
        "plt.plot(...,'r')\n",
        "\n",
        "plt.legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eV8B2x0humJz"
      },
      "source": [
        "#"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PSNcshOlvlhS"
      },
      "source": [
        "(bonus) ***Exercise 4.2*** As you saw from the previous example, the main difference between linear neurons and tanh neurons is that the latter saturate. Also in this case it is possible to derive $\\dot \\kappa = F(\\Sigma)$, with $\\Sigma$ the overlaps between all network vectors. \n",
        "\n",
        "In the case where $\\phi(x) = tanh(x)$ and assuming $N \\rightarrow \\infty$:\n",
        "\n",
        "$$ \\dot \\kappa = -k + \\tilde{\\sigma}_{nm}\\kappa + \\tilde{\\sigma}_{nI}\\kappa_{I} \\; \\; \\; \\small (3)$$\n",
        "$$ \\dot \\kappa_{I} = -\\kappa_{I} + u$$\n",
        "\n",
        "\n",
        "with:\n",
        "\n",
        "$\\tilde{\\sigma}_{ab} = a^Tb/N <\\phi'(0,\\Delta)>$ \n",
        "\n",
        "$\\Delta  =  \\kappa^2 \\sigma_{m}^2 + \\kappa_{I}^2\\sigma_{I}^2$\n",
        "\n",
        "The derivation of this can be checked in (Mastrogiuseppe & Ostojic, 2019) and (Dubreuil, Valente et al, 2019). \n",
        "\n",
        "Intuitively, the dynamics of $\\kappa$ depend no only on the fixed connectivity (i.e. overlaps), but also on the average gain (i.e. $<\\phi'>$) that will change depending on the network state.\n",
        "\n",
        "* Using phi_prime function below, integrate $\\dot \\kappa$ only using the overlaps.\n",
        "* The match between the network (projected on $m$) and $\\kappa$ might be pretty good but probably not perfect. Think about why that is and how to improve it (hint: think about the two major assumptions we made) \n",
        "\n",
        "\n",
        "<p align=\"center\">\n",
        "<img src=\"https://github.com/jmourabarbosa/lowrank_tutorial/blob/main/figures/4.2.png?raw=1\" alt=\"drawing\" width=\"600\"/>\n",
        "</p>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 90,
      "metadata": {
        "id": "cZRR_fPMZ3X3"
      },
      "outputs": [],
      "source": [
        "gaussian_norm = (1/np.sqrt(np.pi))\n",
        "gauss_points, gauss_weights = np.polynomial.hermite.hermgauss(200)\n",
        "gauss_points = gauss_points*np.sqrt(2)\n",
        "\n",
        "def phi (mu, delta0):\n",
        "    integrand = np.tanh(mu+np.sqrt(delta0)*gauss_points)\n",
        "    return gaussian_norm * np.dot (integrand,gauss_weights)\n",
        "\n",
        "def phi_prime (mu, delta0):\n",
        "    integrand = 1 - (np.tanh(mu+np.sqrt(delta0)*gauss_points))**2\n",
        "    return gaussian_norm * np.dot (integrand,gauss_weights)\n",
        "\n",
        "def calc_effective_conn(kappa, kappa_I,I,n,overlaps):\n",
        "\n",
        "\n",
        "  (...)\n",
        "\n",
        "  return s_sigma_mn, s_sigma_nI\n",
        "\n",
        "# two trials, one for each input\n",
        "kappa = [0,0]\n",
        "kappa_I = [0,0]\n",
        "K = np.zeros((2,len(time)))\n",
        "\n",
        "for ti in range(len(time)):\n",
        "\n",
        "    # one input\n",
        "    s_sigma_mn, s_sigma_nI = calc_effective_conn(kappa[0], kappa_I[0],IA,n,lrh.overlaps)\n",
        "    kappa_rec = (...)\n",
        "\n",
        "    K[:,ti] = kappa\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 448
        },
        "id": "LSX8ImnqaypG",
        "outputId": "2f328c7a-670d-48c4-b906-9d2fc59b98a4"
      },
      "outputs": [],
      "source": [
        "plt.plot(...)"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.8.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}