modelling_spatial_representations.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# **Modelling Spatial Representations** [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/TomGeorge1234/ModellingSpatialRepresentations/blob/master/modelling_spatial_representations.ipynb)\n",
    "### **TReND CaMinA, Computational Neuroscience and Machine Learning Summer School, 2024**\n",
    "#### made by: **Tom George, UCL**\n",
    "\n",
    "---\n",
    "\n",
    "Ιn this tutorial we write code to model spatial behaviour and spatial representations in a variety of environments. To do so we'll use the `RatInABox` simulation toolkit (full documentation can be found [here](https://github.com/TomGeorge1234/RatInABox)).\n",
    "\n",
    "In summary, the goal of this package is to make it easy to generate spatial behaviour and spatial representation neuronal data in complex worlds. \n",
    "\n",
    "Researchers use it to: \n",
    "* Generate data to test theories.\n",
    "* Build models of behaviour and Hippocampal networks.\n",
    "* Compare or augment neural data with theoretical models. \n",
    "\n",
    "<center><img src=\"./figures/display_figures/ratinabox.gif\" width=1000></center>\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## **Learning Objectives**\n",
    "* Install and Import the `RatInABox` package for modelling spatial behaviour and representations.\n",
    "* Learn its basic API: be able to build  2D Environments and generate locomotion trajectories for one or multiple Agents _within_ these Environments.\n",
    "* Learn how to generate single cell representations: place cells, grid cells, head direction cells and boundary cells\n",
    "* Built a simple feed-forward threshold model for how grid cells could drive place cells.   "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## **Contents** \n",
    "0. [Import dependencies](#dependencies)\n",
    "1. [Generating trajectories in 2D Environments](#trajectories)\n",
    "    1. [Square environment and random trajectories](#square)\n",
    "    2. [Building complex environments](#complex-environments)\n",
    "    3. [Varying the trajectory parameters](#changing-trajectories)\n",
    "    4. [Controlling the Agent for non-random trajectories](#non-random)\n",
    "    5. [Animating trajectories](#animating)\n",
    "    6. [Multiple Agents](#multiagent)\n",
    "    7. [Simulating the behaviour of TReND students](#trend)\n",
    "2. [Spatial representations](#spatial-reps)\n",
    "    1. [Place cells](#placecells)\n",
    "    3. [Grid cells](#gridcells)\n",
    "    4. [Head direction cells](#headdirectioncells)\n",
    "    5. [Boundary vector cells](#boundaryvectorcells)\n",
    "3. [Grid cells drive place cells: A model](#gc2pc)\n",
    "4. [Advanced exercises](#advanced)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## **0 Install and import `RatInABox`** <a name=\"dependencies\"></a>\n",
    "\n",
    "Below we install and import `RatInABox`. \n",
    "\n",
    "Note the line where we set the `ratinabox.figure_directory = './figures/'`. This is where any figures we make will be automatically saved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "#@title Run this cell to install ratinabox and import the necessary libraries {display-mode: \"form\" }\n",
    "!pip install ratinabox>=1.12.4 pandas\n",
    "\n",
    "import ratinabox\n",
    "from ratinabox.Environment import Environment\n",
    "from ratinabox.Agent import Agent \n",
    "from ratinabox.Neurons import *\n",
    "ratinabox.autosave_plots = True; ratinabox.stylize_plots(); ratinabox.figure_directory = './figures/'\n",
    "import pandas as pd "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = 0.1 # throughout this tutorial we will use a time step of 0.1 seconds"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## **1.0 Generating random trajectories in 2D Environments** <a name=\"trajectories\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "RatInABox consists of three core classes: \n",
    "\n",
    "* `Environment()`: a 1D or 2D arena where the actions take place. \n",
    "* `Agent()`: These live inside the `Environment` and move around.\n",
    "* `Neurons()`: Belong to an `Agent` and report its \"state\" (e.g. place cells or grid cells)\n",
    "\n",
    "<center><img src=\"./figures/display_figures/ratinabox_structure.png\" width=400></center>\n",
    "\n",
    "For now we'll focus on the first two: "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **1.1 Basic square environment with random trajectory** <a name=\"square\"></a>\n",
    "\n",
    "If initialised with no parameters then, by default\n",
    "* The `Environment` will be square\n",
    "* The `Agent` will move under a smooth random motion policy "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Env = Environment(params={})\n",
    "Ag = Agent(Env, params={})"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the environment to see its shape (it's a square!) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Env.plot_environment()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may notice this figure was automatically saved so you can find it later if you want. \n",
    "\n",
    "`Agent.update()` progresses the simulation by a small amount of time `Agent.dt`. So for 60 seconds of simulation..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(int(60/dt)):\n",
    "    Ag.update(dt=dt)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data is saved into a dictionary `Agent.history`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "agent_data = pd.DataFrame(Ag.history)\n",
    "agent_data"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_📝TASK_** Starting from this `pandas` dataframe calculate the average speed of the `Agent`.\n",
    "\n",
    "$$speed = \\left< \\| velocity\\| \\right>$$\n",
    "\n",
    "Hint: you might like to loop over each element in the dataframe: `for vel in dataframe['vel']: ...` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "#@title Double click for the answer {display-mode: \"form\" }\n",
    "speeds = []\n",
    "for vel in agent_data['vel']:\n",
    "    speed = np.linalg.norm(vel)\n",
    "    speeds.append(speed)\n",
    "average_speed = np.mean(speeds)\n",
    "print(f\"The average speed is: {average_speed : .3f} m/s\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1.1 Use `Agent.plot_trajectory()` to visualise the path the Agent took \n",
    "Thankfully `RatInABox` has in-built plotting functions which means we can quickly visualise this data with relatively little code. \n",
    "\n",
    "The most useful one for the `Agent` class is `Agent.plot_trajectory()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ag.plot_trajectory() # plot the trajectory of the agent"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **1.2 Building more _complex_ 2D environments** <a name=\"complex-environments\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_📝TASK_**: Use `Environment.add_wall([[x0,y0],[x1,y1])` to add a single wall to the `Environment`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code goes here\n",
    "Env = Environment()\n",
    "Env.add_wall([[0.3,0.6],[0.3,0]]) #<-- add you wall in here\n",
    "Env.plot_environment()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice we left the `params` dictionaries empty giving us a default square `Environment`. Here are the default values: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Env.get_all_default_params()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_📝TASK_**: By passing the `boundary` param into a new `Environment` make a non-square `Environment`\n",
    "\n",
    "The `boundary` parameter should be a list of corners definin the non-square `Environment`\n",
    "Ideas: \n",
    "* Make a T-Maze\n",
    "* Octagon `Environment`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Env = Environment(\n",
    "    params={\n",
    "        'boundary': [[0.25,0],[0,0.25],[0,0.75],[0.25,1],[0.75,1],[1,0.75],[1,0.25],[0.75,0]]#You list goes here\n",
    "        }\n",
    ")\n",
    "Env.plot_environment()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_📝TASK_**: Pass the `holes` parameter (this is a list of lists, each inner list element giving the corners of a \"hole\") to make an `Environment` with a hole"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code goes here \n",
    "Env = Environment(params={'holes':[]}) #<-- add one or more holes into this list \n",
    "Env.plot_environment()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "#@title Double click for the answer {display-mode: \"form\" }\n",
    "# your code goes here \n",
    "Env = Environment(params={'holes':[[[0.2,0.2],[0.4,0.2],[0.4,0.4],[0.2,0.4]]]}) #<-- add one or more holes into this list \n",
    "Env.plot_environment()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **1.3 Varying the trajectory parameters** <a name=\"changing-trajectories\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice we left the `params` dictionaries empty giving us a default `Agent`. Here are the default values: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Agent.get_all_default_params()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, the `Agent` parameter `rotational_velocity_std` controls how curvy the trajectory is. \n",
    "\n",
    "**_📝TASK_**: By passing the `rotational_velocity_std` (or others) param into a new `Agent` vary the trajectory statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Env = Environment()\n",
    "Ag = Agent(Env,params={}) #<--- add parameters here\n",
    "for i in range(int(60/dt)):\n",
    "    Ag.update(dt=dt)\n",
    "Ag.plot_trajectory()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **1.4 Non-random trajectories** <a name=\"non-random\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, it's possible to \"guide\" the `Agent` trajectory, rather than just having a random motion.\n",
    "\n",
    "To do this we set pass a `drift_velocity` to the Agent on each update set: `Agent.update(drift_velocity=[vx,vy])`\n",
    "\n",
    "In this example we'll set a drift velocity so that an `Agent` in a T-Maze moves up the central arm then either left or right. Upon reaching the end of the arm it teleports back to the bottom of the central arm. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Env = Environment(\n",
    "    params={\n",
    "        'boundary': \n",
    "        [[-0.1,0],\n",
    "         [0.1,0],\n",
    "         [0.1,0.8],\n",
    "         [0.7,0.8],\n",
    "         [0.7,1.0],\n",
    "         [-0.7,1.0],\n",
    "         [-0.7,0.8],\n",
    "         [-0.1,0.8]]})\n",
    "Ag = Agent(Env)\n",
    "Env.plot_environment()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def drift_velocity_from_position(pos : np.ndarray):\n",
    "    \"\"\"Returns a drift velocity that will move the agent up the central arm and then along the T-arms\"\"\"\n",
    "    if pos[1] < 0.8: return 0.2*np.array([0,1])\n",
    "    else: \n",
    "        if pos[0] < 0: return 0.2*np.array([-1,0])\n",
    "        elif pos[0] >= 0: return 0.2*np.array([1,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(int(100/dt)):\n",
    "    Ag.update(dt=dt,drift_velocity=drift_velocity_from_position(Ag.pos)) # drift along the arms\n",
    "    if np.abs(Ag.pos[0]) > 0.65: # teleport back and reset the velocity to zero\n",
    "        Ag.pos = np.array([0,0.05])\n",
    "        Ag.velocity = np.array([0,0])\n",
    "\n",
    "Ag.plot_trajectory(color='changing')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_📝TASK_**: Make a default square environment and set a drift velocity in such a way that the Agent moves around in circles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code goes here "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "#@title Double click for the answer {display-mode: \"form\" }\n",
    "Env = Environment()\n",
    "Ag = Agent(Env)\n",
    "for i in range(int(100/dt)):\n",
    "    [x,y] = Ag.pos - np.array([0.5,0.5]) #distance from the centre of the Environment\n",
    "    [vx,vy] = [-0.3*y,0.3*x]\n",
    "    Ag.update(dt=dt, drift_velocity=np.array([vx,vy]))\n",
    "Ag.plot_trajectory()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **1.5 Animating trajectories** <a name=\"animating\"></a>\n",
    "Finally, we can animate the trajectory to see the `Agent` moving around the `Environment` in real time. \n",
    "Note: this works best if running locally and can be a little slow on Google colab. Feel free to skip this! \n",
    "\n",
    "To do this use the `Agent.animate_trajectory()` function. I recommend passing some parameters: \n",
    "* `speed_up`: How many times real speed the animation will play\n",
    "* `t_end`: when, in Agent time coordinates, to plot until\n",
    "* `fps`: frame per second of resulting animation (set lower for faster run times)\n",
    "* `color`: The color of  the trajectory (any `matplotlib` color). If its passed as `\"changing\"` it will shade the trajectory.\n",
    "* `autosave`: Whether the animation is automatically saved to the figure directory (recommend `False` to save time) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ag.animate_trajectory(speed_up=10,t_end=60,fps=15,color=\"changing\",autosave=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_📝TASK_**: Go back up to the T-Maze example we built and animate the trajectory of the `Agent`. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **1.6 Multiple agents** <a name=\"multiagent\"></a>\n",
    "Finally we can add multiple Agents to the environment, no problem. \n",
    "When ploting, specify `plot_all_agents=True` to see them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Env = Environment()\n",
    "Ag1 = Agent(Env)\n",
    "Ag2 = Agent(Env)\n",
    "Ag3 = Agent(Env)\n",
    "for i in range(int(60/dt)):\n",
    "    Ag1.update(dt=dt)\n",
    "    Ag2.update(dt=dt)\n",
    "    Ag3.update(dt=dt)\n",
    "\n",
    "Ag1.plot_trajectory(plot_all_agents=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_📝OPTIONAL TASK_**: Make an Environment with multiple `Agents` who \"interact\".\n",
    " \n",
    "Consider how we could set a `drift_velocity` for each Agent which might make them approach or run away from each another. If you feel comfortable, give it a go! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "#@title Double click to see my solution to this problem (there might be many ways of doing it) \n",
    "\n",
    "# here the algorithm I use is as follows: \n",
    "# for each agent, loop over all the _other_ agents and set a drift velocity which points towards them \n",
    "# you drift more strongly the closer you are to them\n",
    "\n",
    "Env = Environment()\n",
    "Ag1 = Agent(Env)\n",
    "Ag2 = Agent(Env)\n",
    "Ag3 = Agent(Env)\n",
    "Ag4 = Agent(Env)\n",
    "\n",
    "All_Agents = [Ag1,Ag2,Ag3,Ag4]\n",
    "for _ in range(int(120/dt)):\n",
    "    for i in range(len(All_Agents)):\n",
    "        drift_vel = np.array([0.0,0.0])\n",
    "        Me = All_Agents[i]\n",
    "        other_Ags = All_Agents.copy(); other_Ags.pop(i) #remove yourself from the list of other agents\n",
    "        for other_Ag in other_Ags:\n",
    "            vec_to_other_agent = np.array(other_Ag.pos) - np.array(Me.pos)\n",
    "            dist_to_other_agent = np.linalg.norm(vec_to_other_agent)\n",
    "            norm_vec = vec_to_other_agent / dist_to_other_agent\n",
    "            attract_speed = 0.005*np.exp(-(dist_to_other_agent/0.2)**2/2)\n",
    "            drift_vel += attract_speed*norm_vec\n",
    "        Me.update(dt=dt,drift_velocity=drift_vel)\n",
    "\n",
    "Ag1.animate_trajectory(speed_up=10,fps=10,plot_all_agents=True,autosave=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **1.7 Simulating the behaviour of TReND students...** <a name=\"trend\"></a>\n",
    "### **...when they hear about an awesome new summer school happening in Accra**\n",
    "\n",
    "Let's piece this all together and make one simulation showing the great migration of students from all over Africa to Kigali for the TReND summer school \n",
    "\n",
    "We'll go through it step-by-step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1 get the coordinates of the boundary of Africa (\n",
    "# These are stored as a massive list in the very last cell of this notebook\n",
    "# Go and run that cell and then come back here\n",
    "\n",
    "kigali = np.array([0.4,-0.3]) # coordinates of Kigali\n",
    "\n",
    "# Step 2: create an Environment object with the boundary of Africa as the boundary\n",
    "Africa = Environment(params={\n",
    "    'dx': 0.05,\n",
    "    'boundary': africa_boundary\n",
    "})\n",
    "\n",
    "\n",
    "# Step 3: create many Agents (or students) within the Environment you just created\n",
    "# put them all in a list called Students\n",
    "Ag0 = Agent(Africa)\n",
    "Ag1 = Agent(Africa)\n",
    "Ag2 = Agent(Africa)\n",
    "Ag3 = Agent(Africa)\n",
    "Ag4 = Agent(Africa)\n",
    "Ag5 = Agent(Africa)\n",
    "Ag6 = Agent(Africa)\n",
    "Ag7 = Agent(Africa)\n",
    "Ag8 = Agent(Africa)\n",
    "Students = [Ag0,Ag1,Ag2,Ag3,Ag4,Ag5,Ag6,Ag7,Ag8]\n",
    "\n",
    "\n",
    "# Step 4: Simulate!\n",
    "for i in range(int(60/dt)):\n",
    "    for Student in Students:\n",
    "        if Student.t < 30: # initially students just randomly drift around Africa\n",
    "            Student.update(dt=dt)\n",
    "        else: # after 30 seconds students drift towards kigali\n",
    "            vector_to_kigali = kigali - Student.pos\n",
    "            vector_to_kigali = vector_to_kigali / np.linalg.norm(vector_to_kigali)\n",
    "            drift_velocity = 0.1 * vector_to_kigali\n",
    "            Student.update(dt=dt,drift_velocity=drift_velocity)\n",
    "\n",
    "# Step 5: Plot the trajectories of all the students\n",
    "Students[0].plot_trajectory(plot_all_agents=True)\n",
    "\n",
    "\n",
    "# or animate but...this may be VERY slow to run on Google Colab.\n",
    "# Students[0].animate_trajectory(speed_up=5,fps=10,plot_all_agents=True,autosave=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## **2.0 Spatial Representations** <a name=\"spatial-reps\"></a>\n",
    "\n",
    "Great! So far we've\n",
    "* Successfully been able to make `Environment`s\n",
    "* We then added `Agent`s to these environments\n",
    "* We simulate these `Agent`s moving around under complex, realistic motion policies.\n",
    "\n",
    "But what about what goes on inside the brain? We're neuroscientists after all!\n",
    "\n",
    "Now we'll add `Neurons` to `Agents`. Whereas `Agents` have continuously updating _positions_ and _velocities_ in an Environment, `Neurons` have continuously updating _firing rates_. \n",
    "\n",
    "As we learnt during the lecture, certain types of neurons e.g. place cells or grid cells, have localised receptive field (they only fire in specific regions of the environment), others e.g. head direction cells, are specific to aspects of how the agent is moving. \n",
    "\n",
    "Here are some example `Neurons` subclasses we can add to `Agents`\n",
    "\n",
    "* `PlaceCells`\n",
    "* `GridCells`\n",
    "* `BoundaryVectorCells`\n",
    "* `HeadDirectionCells`\n",
    "\n",
    "Just like `Agent`s, each is initialised with a dictionary of `params` and can be updated by calling `Neuron.update()`. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **2.1 Place cells** <a name=\"placecells\"></a>\n",
    "\n",
    "Here we'll make a `PlaceCell` which has a simple `gaussian` receptive field: \n",
    "\n",
    "$$ f_i(\\vec{x}) = \\textrm{fr}_{\\textrm{max}}\\exp \\big( -\\frac{(\\vec{x}(t) - \\vec{x}_{i})^2}{2 \\sigma^2} \\big)$$ \n",
    "\n",
    "$\\vec{x}_i$ is the centre of the Gauassian and $\\sigma$ is the width. \n",
    "\n",
    "Let's see how we'd go about making this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Env = Environment()\n",
    "Ag = Agent(Env)\n",
    "PCs = PlaceCells(Ag, params={'n':5, #<--- 5 place cells\n",
    "                             'max_fr':5, #<--- max firing rate = 5Hz\n",
    "                             'widths':0.1, #<--- sigma = 0.1 m\n",
    "                                     }) "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note we didn't specify the centre of the place fields, $x_i$, so they are chosen randomly. \n",
    "\n",
    "Now lets simulate 10 mins seconds of exploration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(int(10*60/dt)):\n",
    "    Ag.update(dt=dt)\n",
    "    PCs.update() # this line updates the firing rate of the place cells"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the data is saved into a dictionary, `Neurons.history`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "neuron_data = pd.DataFrame(PCs.history)\n",
    "neuron_data"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_📝TASK_**: Calculate the average firing rate of one of the neurons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "#@title Double click to see solution {display-mode: \"form\" }\n",
    "\n",
    "# To calculate the firing rate we can just see how many spikes we observed in the 10 minutes and divide by 10 minutes\n",
    "# we'll do this for the first neuron only. \n",
    "spikes_bool = [spike_list[0] for spike_list in neuron_data['spikes']]\n",
    "num_spikes = np.sum(spikes_bool)\n",
    "av_firing_rate = num_spikes / (10*60)\n",
    "\n",
    "print(f\"Total number of spikes: {num_spikes:g}\\nAverage firing rate (spikes / second): {av_firing_rate:.3f}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`RatInABox` has in-built plotting functions so we can explore this data easily and without much code. These include: \n",
    "\n",
    "| Function | Description | \n",
    "| --- | --- |\n",
    "| `Neurons.plot_rate_timeseries()` | show the firing rate of the neurons over time |\n",
    "| `Neurons.plot_rate_map()` | shows the ground truth firing rate of the neurons as a function of the position in the Environment |\n",
    "| `Neurons.plot_rate_map(method='neither',spikes=True)` | shows where the neurons spiked, from the data already collected |\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = Ag.plot_trajectory()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PCs.plot_rate_timeseries(t_end=60)\n",
    "fig.suptitle(\"The observed spiking rate of the neuron over time \")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PCs.plot_rate_map()\n",
    "fig.suptitle(\"This is the ground-truth receptive field of the place cell (gaussian)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = PCs.plot_rate_map(method='neither',spikes=True)\n",
    "fig.suptitle(\"The is the data (spikes) after 10 minutes of simulation\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To animate this we can call `Agent.animate_trajectory` but passing an additional function which displays the location of the spikes as they appear. \n",
    "\n",
    "I've written this function below. Expect this to take about a minute to run, feel free to skip."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#you can ignore this function for now\n",
    "def plot_spikes_on_animation(fig, ax, t, **kwargs):\n",
    "    kwargs['cells'].plot_rate_map(fig=fig,ax=ax,method='neither',spikes=True,t_end=t,autosave=False,chosen_neurons=kwargs['chosen_neurons'])\n",
    "    return fig, ax\n",
    "\n",
    "Ag.animate_trajectory(speed_up=30,\n",
    "                      fps=10,\n",
    "                      additional_plot_func=plot_spikes_on_animation,\n",
    "                      cells=PCs,\n",
    "                      chosen_neurons='1',\n",
    "                      autosave=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_📝TASK_**: Start playing around with these cells, what do you see?\n",
    "\n",
    "Good ideas include\n",
    "\n",
    "* Number: Try simulating more than place cells, e.g., 10, using `params['n']=10`\n",
    "* Widths: Vary the widths of the place cells using `params['widths']`\n",
    "* Walls: Add a single wall to the Environment using `Env.add_wall([[x0,y0],[x1,y1]])`, how do the place cells respond to the wall? \n",
    "* Models: The parameter `params['description']` gets the specific place cell model being used. Try varying this and see what you observe, options include: \n",
    "    * `\"gaussian\"` (default)\n",
    "    * `\"one_hot\"`\n",
    "    * `\"gaussian_threshold\"`\n",
    "* Firing rate: Vary the firing rate of the place cells `params[\"max_fr\"]`"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's some basic code to get you started"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Env = Environment()\n",
    "Ag = Agent(Env)\n",
    "PCs = PlaceCells(Ag,params={\n",
    "\n",
    "    # insert your params here\n",
    "})\n",
    "\n",
    "# 10 min simulation \n",
    "for i in range(int(10*60/dt)):\n",
    "    Ag.update(dt=dt)\n",
    "    PCs.update()\n",
    "\n",
    "#display the data\n",
    "PCs.plot_rate_timeseries()\n",
    "PCs.plot_rate_map(method='neither',spikes=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **2.2 Grid cells** <a name=\"gridcells\"></a>\n",
    "\n",
    "Other Spatial Representations exist, for example grid cells. The class `GridCells` handles this. \n",
    "\n",
    "The default model for grid cells is to sum three cosine wave together \n",
    "\n",
    "$$ f_{i}(\\vec{x}) = \\frac{1}{3} \\textrm{fr}_{\\textrm{max}} \\bigg[ \\textrm{cos}\\bigg( 2\\pi\\frac{\\vec{x}(t)\\cdot \\vec{e}_{\\theta_{i}}}{\\lambda_i} + \\phi_{i}\\bigg) + \\textrm{cos}\\bigg( 2\\pi\\frac{\\vec{x}(t)\\cdot \\vec{e}_{\\theta_{i} + \\pi/3}}{\\lambda_i} + \\phi_{i}\\bigg) + \\textrm{cos}\\bigg( 2\\pi\\frac{\\vec{x}(t)\\cdot \\vec{e}_{\\theta_{i}+2\\pi/3}}{\\lambda_i} + \\phi_{i}\\bigg)  \\bigg]_{+} $$\n",
    "\n",
    "<center><img src=\"./figures/display_figures/gridcell_oim.png\" width=800></center>\n",
    "\n",
    "In this model, each grid cell, $i$, has an \n",
    "* Orientation, $\\theta_i$\n",
    "* Phase offset, $\\phi_i$\n",
    "* Grid scale, $\\lambda_i$\n",
    "\n",
    "**_📝TASK_**: Initialise, simulate and then visualise some grid cells. To see the available parameters do `GridCells.get_all_default_params()`\n",
    "\n",
    "Extensions once you've done this: \n",
    "* Howdo grid cells respond to adding a wall to the Environment?\n",
    "* By setting `param['gridscale_distribution']` to `'delta'` and `param['gridscale']` to whatever-value-you-like make a set of grid cells which all have the same gridscale.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code goes here "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "#@title Double click to see solution {display-mode: \"form\" }\n",
    "Env = Environment()\n",
    "Ag = Agent(Env)\n",
    "GCs = GridCells(Ag,params={'n':5,'max_fr':5})\n",
    "\n",
    "for i in range(int(10*60/dt)):\n",
    "    Ag.update(dt=dt)\n",
    "    GCs.update()\n",
    "\n",
    "GCs.plot_rate_timeseries()\n",
    "GCs.plot_rate_map()\n",
    "GCs.plot_rate_map(method='neither',spikes=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **2.3 Head Direction Cells** <a name=\"headdirectioncells\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Head direction cells are the brains internal compass. They allow us to know where we're facing and are therefore crucial to the brains ability to path integrate. \n",
    "\n",
    "**_📝TASKS_** \n",
    "1. Make some `HeadDirectionCells` and explore a 2D environment for 10 minutes\n",
    "2. Visualise their receptive fields in polar coordinates using `HeadDirectionCells.plot_HDC_receptive_field()`\n",
    "3. Plot their rate maps using `HeadDirectionCells.plot_rate_map(method='history')`...what do you see\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "#@title Double click to see solution {display-mode: \"form\" }\n",
    "Env = Environment()\n",
    "Ag = Agent(Env)\n",
    "HDCs = HeadDirectionCells(Ag,params={'n':10,'max_fr':5})\n",
    "\n",
    "for i in range(int(10*60/dt)):\n",
    "    Ag.update(dt=dt)\n",
    "    HDCs.update()\n",
    "\n",
    "\n",
    "HDCs.plot_HDC_receptive_field()\n",
    "HDCs.plot_rate_map(method=\"history\")\n",
    "HDCs.plot_rate_timeseries(t_end=60)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **2.4 Boundary Vector Cells** <a name=\"boundaryvectorcells\"></a>\n",
    "\n",
    "The last cell type we discussed was boundary vector cells. These fire at a fixed angle and distance from the walls.\n",
    "\n",
    "`RatInABox` contains these as `BoundaryVectorCells()`\n",
    "\n",
    "**_📝TASK_** \n",
    "Remember the Africa `Environment` we built. Let's add some BVCs to that environment and see how they look. Let's start with 5 cells "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Africa = Environment(params={'boundary':africa_boundary,\n",
    "                             'dx':0.05 #<--- 5cm resolution (ignore this, it just speeds up the code)\n",
    "                             })\n",
    "Ag = Agent(Africa)\n",
    "BVCs = BoundaryVectorCells(Ag,params={'n':5})\n",
    "\n",
    "BVCs.plot_rate_map()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## **3 Grid cells to place cells** <a name=\"gc2pc\"></a>\n",
    "### **A model of the place cell formation**\n",
    "\n",
    "As we discussed in the lecture, it is thought that grid cells provide a primary source of input to place cells. \n",
    "\n",
    "Instead of using the `RatInABox` `PlaceCells` class, can we instead sum together `GridCells` in a way which makes place cells??\n",
    "\n",
    "<center><img src=\"./figures/display_figures/ec_to_hpc.png\" width=400></center>\n",
    "\n",
    "To do this we'll using the `RatInABox` class called `FeedForwardLayer`. \n",
    "\n",
    "<center><img src=\"./figures/display_figures/feedforwardlayer.png\" width=800></center>\n",
    "\n",
    "`FeedForwardLayer` is quite unique. Unlike `PlaceCells`, `GridCells` etc. who's firing rates are directly determined by ther `Agents` state, this layer takes as input one (or many) other layer, $g_j(t)$ and calculates its firing rate $p_i(t)$$ by summing and \"activating\" these inputs (like a layer in a deep neural network if you're familiar with those) according to: \n",
    "\n",
    "$$\n",
    "p_{i}(t) = \\sigma \\big(\\sum_j w_{ij} g_{j}(t) \\big)\n",
    "$$\n",
    "\n",
    "The input layer can be _any other `RatInABox` `Neurons` class_, for example `GridCells`. In which case the new layer calculates an activated linear sum of GridCell inputs. Perhaps this is a good model for hippocampus? Lets see..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Env = Environment()\n",
    "Ag = Agent(Env)\n",
    "GCs = GridCells(Ag,params={'gridscale':np.logspace(np.log10(0.05),np.log10(1),5),\n",
    "                           'description':'three_shifted_cosines',\n",
    "                           'orientation_distribution':'delta',\n",
    "                           'phase_offset_distribution':'delta'})\n",
    "PCs = FeedForwardLayer(Ag, params={\n",
    "    'n':1,\n",
    "    'input_layers':[GCs],\n",
    "    'activation_params':{'activation':'relu'},\n",
    "})\n",
    "\n",
    "fig, ax = GCs.plot_rate_map(chosen_neurons='10')\n",
    "fig.suptitle(\"Grid cells: the inputs to our 'place cells'\")\n",
    "\n",
    "fig, ax = PCs.plot_rate_map(chosen_neurons='5')\n",
    "fig.suptitle(\"Place cells as a random linear sum of the grid cells. Do these look like place cells? (They shouldn't!)\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe we can do better by thresholding the placec ells near their maximum values .This could be done in the brain "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratemaps = PCs.get_state(evaluate_at='all')\n",
    "max_of_ratemaps = np.max(ratemaps,axis=1)\n",
    "PCs.biases = -0.5*max_of_ratemaps\n",
    "\n",
    "PCs.plot_rate_map()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok...not great!! Some of them look \"place cell-y\" but most don't. Perhaps this is because we just chose random weights. It's reasonable to assume the brain wouldn't actually set weights randomly. \n",
    "\n",
    "In fact, if you read Solstad et al. (2006) they derived a model for how strong a grid cell should connect to a place cell. It turns out to only depend on the \"grid scale\" $\\lambda$ of the grid cell and the desired size $\\sigma_{pc}$ of the place cells. Roughly, their model suggest weights should be set proportional to: \n",
    "\n",
    "$$w_{ij} \\propto \\frac{\\exp\\bigg(-\\frac{4\\pi^{2}\\sigma_{pc,i}^{2}}{3\\lambda_j^{2}}\\bigg)}{\\lambda_{j}^2}$$\n",
    "\n",
    "where $\\lambda_{j}$ is the gridscale and $\\sigma_{pc,i}$ is the size of the place cell. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Env = Environment(params={'boundary':[[-0.5,-0.5],[0.5,-0.5],[0.5,0.5],[-0.5,0.5]]})\n",
    "Ag = Agent(Env,params = {'dt':0.1})\n",
    "GCs = GridCells(Ag,params={'gridscale':np.logspace(np.log10(0.05),np.log10(1),50),\n",
    "                           'description':'three_shifted_cosines',\n",
    "                           'orientation_distribution':'delta',\n",
    "                           'phase_offset_distribution':'delta'})\n",
    "PCs = FeedForwardLayer(Ag, params={\n",
    "    'n':1,\n",
    "    'input_layers':[GCs],\n",
    "    'activation_params':{'activation':'relu'},\n",
    "})"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_📝TASK_** Complete the following code to implement to Solstad model. Use the equation above to help (note there's only 1 place cell) so you can ignore the index i. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gridscales = GCs.gridscales\n",
    "sigma = 0.2\n",
    "weights = #your code goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "#@title Double click to see solution {display-mode: \"form\" }\n",
    "gridscales = GCs.gridscales\n",
    "sigma = 0.2\n",
    "weights = np.exp(-(4*np.pi**2*sigma**2)/(3*gridscales**2))/gridscales**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Here we set the weights by accessing the weight matrix. \n",
    "PCs.inputs['GridCells']['w'] = weights\n",
    "\n",
    "#Threshold at 50% the max \n",
    "ratemaps = PCs.get_state(evaluate_at='all')\n",
    "max_of_ratemaps = np.max(ratemaps,axis=1)\n",
    "PCs.biases = -0.5*max_of_ratemaps\n",
    "\n",
    "fig, ax = PCs.plot_rate_map()\n",
    "fig.suptitle(\"A place cell, calculated as a weighted linear sum of the grid cells.\\nDoes it look like a place cell? (It should!)\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## **4.0 Advanced exercises** <a name=\"advanced\"></a>\n",
    "\n",
    "In any remaining time here's a list of task you might like to try. Pick one you think looks fun and give it a go (the content above should be enough to get you there):\n",
    "\n",
    "**_📝TASK A_** Make a rectangular environment with two square holes so that it forms a figure-of-eight. Then, add an agent to the environment and by designing an appropriate `drift_velocity` make it move around the arms of the figure-of-eight.\n",
    "\n",
    "**_📝TASK B_** Initialise a population of `GridCells` and a population of `HeadDirectionCells`. Using a `FeedForwardLayer` to create neurons which \"mix\" these two cell types. How might this be related to the role of the entorhinal cortex in path integration (i.e. dead reckoning). \n",
    "\n",
    "**_📝TASK C_** Add one or many objects to the `Environment` using `Env.add_object(position=[x,y])` (I didn't tell you about this functionality). Then initialise a population of `ObjectVectorCells` and plot their receptive field. What do you observe. \n",
    "\n",
    "**_📝TASK D_** (Very hard) Generate data from a population of 10 place cells as an `Agent` randomly explores a 2D Environment. Extract the firing rate data from the `PlaceCells.history` data frame. Extract the position data from the `Agent.history` data frame. Trai na regression decoder to _predict_ the position from the firing rate. You may like to check out [scikit-learns linear regression model](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html) to do this. We will do something like this,and more, next Thursday. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## **Congrats!** \n",
    "Well done. You made it to the end of the tutorial! \n",
    "\n",
    "### **Re-use**\n",
    "Feel free you to adapt and use this tutorial for your own teaching needs! \n",
    "\n",
    "### **About the author: Tom George**\n",
    "* Feel free to get in touch at tom.george.20@ucl.ac.uk \n",
    "* Links: [Twitter](https://twitter.com/TomNotGeorge), [Github](https://github.com/TomGeorge1234), [Google Scholar](https://scholar.google.com/citations?user=AG49j3MAAAAJ&hl=en)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "africa_boundary = np.array([\n",
    "[-0.83935318,  0.30607098],\n",
    "[-0.83532515,  0.18512068],\n",
    "[-0.84741901,  0.06820623],\n",
    "[-0.75063883, -0.06080795],\n",
    "[-0.61756729, -0.16160638],\n",
    "[-0.48852379, -0.1414467 ],\n",
    "[-0.43207275, -0.13338087],\n",
    "[-0.29093539, -0.10919315],\n",
    "[-0.27076592, -0.16160638],\n",
    "[-0.17398574, -0.16160638],\n",
    "[-0.16591991, -0.25836702],\n",
    "[-0.20221126, -0.31480827],\n",
    "[-0.09736524, -0.41559693],\n",
    "[-0.0651117 , -0.58493048],\n",
    "[-0.12560053, -0.70991663],\n",
    "[-0.07720555, -0.82682718],\n",
    "[-0.01671671, -1.02034844],\n",
    "[ 0.05183795, -1.1251749 ],\n",
    "[ 0.03536411, -1.19568492],\n",
    "[ 0.08340712, -1.21970643],\n",
    "[ 0.18979787, -1.19225328],\n",
    "[ 0.26873546, -1.18195835],\n",
    "[ 0.40944265, -1.0344172 ],\n",
    "[ 0.40944265, -0.95893081],\n",
    "[ 0.47464389, -0.91776086],\n",
    "[ 0.47464389, -0.7907998 ],\n",
    "[ 0.62222414, -0.67071181],\n",
    "[ 0.57760299, -0.40650454],\n",
    "[ 0.69428866, -0.25210011],\n",
    "[ 0.76293132, -0.21779345],\n",
    "[ 0.88991193, -0.00505106],\n",
    "[ 0.87618536,  0.03955053],\n",
    "[ 0.71144688, -0.00848271],\n",
    "[ 0.67026715,  0.0292556 ],\n",
    "[ 0.59132956,  0.13905648],\n",
    "[ 0.56387641,  0.20768056],\n",
    "[ 0.52955997,  0.22483683],\n",
    "[ 0.53642326,  0.30375487],\n",
    "[ 0.40944265,  0.51992107],\n",
    "[ 0.46434895,  0.47531556],\n",
    "[ 0.47464389,  0.52678436],\n",
    "[ 0.44719073,  0.57138987],\n",
    "[ 0.34423163,  0.57482151],\n",
    "[ 0.33050506,  0.5542336 ],\n",
    "[ 0.13488179,  0.61599635],\n",
    "[ 0.08683877,  0.59884009],\n",
    "[ 0.09027041,  0.54737129],\n",
    "[ 0.05938562,  0.53021503],\n",
    "[-0.02641526,  0.58511449],\n",
    "[-0.15682752,  0.65030791],\n",
    "[-0.12937436,  0.68118879],\n",
    "[-0.14996423,  0.75324452],\n",
    "[-0.21517524,  0.74295057],\n",
    "[-0.42108367,  0.72579528],\n",
    "[-0.46570482,  0.68118879],\n",
    "[-0.54120099,  0.68462044],\n",
    "[-0.56866392,  0.70863901],\n",
    "[-0.59611707,  0.64687724],\n",
    "[-0.66474996,  0.61599635],\n",
    "[-0.67162302,  0.54394063],\n",
    "[-0.7402559 ,  0.46502259],\n",
    "[-0.79174034,  0.40669148],\n",
    "[-0.82948843,  0.32777344]])\n",
    "\n",
    "accra = np.array([-0.4,-0.1])\n",
    "kigali = np.array([0.2,0.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}