{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Simulation: Complex social dynamics in a few lines of code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will create a model depicting competition between two cultural traits within a common population. This is a typical cultural dynamics scenario where individuals must adopt one option amongst two or more **mutually exclusive** options (i.e. religion, elections, football teams, ...). In this case we are interested in situations when you have to choose one option (e.g., you cannot practice two religions), but more complex versions with individuals adopting more than one trait can easily be developed.\n",
    "\n",
    "Individuals can change their choice over time. The decision is based on the payoff of each trait. This payoff is a measure of the relative interest of the trait, based on: \n",
    "\n",
    "a) how many people exhibits the trait and \n",
    "\n",
    "b) the attractiveness of the trait.\n",
    "\n",
    "An example of this dynamic could be a competition between two different religions. The number of people practicing a belief makes this belief more appealing. However, some beliefs could be intrinsically more interesting for some individuals so part of the population could adopt them even if they are a minority. Finally, social norms are not static so the attractiveness of specific beliefs can vary over time.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "![general overview of the model](diagram.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time is divided in discrete steps starting at t = 0. \n",
    "\n",
    "At each step t the two populations $A_t$ and $B_t$ are updated as individuals move from A to B and from B to A. Take into accout that the value here can be negative meaning that more people move from B to A.\n",
    "\n",
    "\\begin{equation*}\n",
    "A_{t+1} = A_t + \\Delta AB\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*}\n",
    "B_{t+1} = B_t - \\Delta AB\n",
    "\\end{equation*}   \n",
    "\n",
    "Let's see how we can express it in code. \n",
    "\n",
    "<!--XRC: no needed because they all use python 3.2\n",
    "We start with an auxiliary line of code which prevents truncation when dividing integers. You don't have to worry about it now but just keep in mind to start any code written in Python 2.7 or earlier with this line. If you have installed Python 3.2 or higher you can skip this line.\n",
    "from __future__ import division\n",
    "-->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we need to define the number of individuals in the population. Say we want to start with 100 people. The text after a hash symbol is just a comment so you can skip it for now, but in general it is a good idea to keep documenting the code as you write."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 100   # total population size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we decide on how many believers of each cultural option (religion) we want to start with. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = 65    # initial number of believers A\n",
    "B = N - A # initial number of believers A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we want to update these quantities at every time step depending on some variation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'variation' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-3-3a54da9d0f70>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[0mMAX_TIME\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m100\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;32mwhile\u001b[0m \u001b[0mt\u001b[0m \u001b[1;33m<\u001b[0m \u001b[0mMAX_TIME\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 4\u001b[1;33m     \u001b[0mA\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mA\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mvariation\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      5\u001b[0m     \u001b[0mB\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mB\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mvariation\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      6\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'variation' is not defined"
     ]
    }
   ],
   "source": [
    "t = 0\n",
    "MAX_TIME = 100\n",
    "while t < MAX_TIME:\n",
    "    A = A + variation\n",
    "    B = B - variation\n",
    "    \n",
    "    # advance time to next iteration\n",
    "    t = t + 1    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Type the code in the code tab. Keep in mind that indents are important. \n",
    "\n",
    "Run the code by hitting F5 (or click on a green triangle). What happens? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well, nothing happens or, rather, we get an error. \n",
    "```python\n",
    "NameError: name variation is not defined \n",
    "```\n",
    "Indeed, we have not defined what do we mean by 'variation'. So let's calculate it based on the population switching trait based on a comparison between payoffs. For example if B has higher payoff then A then we should get something like this:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\Delta_{A\\to B} = A · (payoff_{A\\to B} - payoff_{B\\to A})\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the proportion of population A that switches to B is proportional to the difference between payoffs. As we mentioned the payoff of a trait is determined by the population exhibiting the competing trait as well as its intrinsic attractiveness."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To define the payoff we need to implement the following competition equations:\n",
    "\n",
    "\\begin{equation*}\n",
    "Payoff_{B\\to A} = \\frac{A_t}{N}\\ \\frac {T_A} {(T_A + T_B)}\\\n",
    "\\end{equation*} \n",
    "\n",
    "\n",
    "\\begin{equation*}\n",
    "Payoff_{A\\to B} = \\frac{B_t}{N}\\ \\frac {T_B} {(T_A + T_B)}\\\n",
    "\\end{equation*} \n",
    "\n",
    "Let's look at the equations a bit more closely. The first term is the proportion of the entire population N holding a particular cultural trait ($\\frac {A_t}{N}$ for $A$ and $\\frac {B_t}{N}$ for $B$). While the second element of the equations is the balance between the attractiveness of both ideas ($T_A$ and $T_B$) expressed as the attractiveness of the given trait in respect to the total 'available' attractiveness ($T_A + T_B$). \n",
    "\n",
    "You have probably immediately noticed that these two equations are the same in structure and only differ in terms of what is put into them. Therefore, to avoid unnecessary hassle we will create a 'universal' function that can be used for both. Type the code below at the beginning of your script:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def payoff(believers, Tx,Ty):    \n",
    "    proportionBelievers = believers/N\n",
    "    attraction = Tx/(Ty + Tx)\n",
    "    return proportionBelievers * attraction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's break it down a little. First we define the function and give it the input - the number of believers and the two values that define how attractive each cultural option is. \n",
    "```python\n",
    "def payoff(believers, Tx, Ty):\n",
    "```\n",
    "Then we calculate two values:\n",
    "\n",
    "1. percentage of population sharing this cultural option\n",
    "```python\n",
    "    proportionBelievers = believers/N    \n",
    "```    \n",
    "2. how attractive the option is\n",
    "```python\n",
    "    attraction = Tx/(Ty + Tx)\n",
    "```    \n",
    "\n",
    "Look at the equations above, the first element is just the $\\frac {A_t}{N}$ and $\\frac {B_t}{N}$ part and the second is this bit: $\\frac {T_A} {(T_A + T_B)}$. \n",
    "Finally we return the results of the calculations. \n",
    "```python\n",
    "    return proportionBelievers * attraction\n",
    "```\n",
    "\n",
    "Voila! We have implemented the equation into Python code. Now, let's modify the main loop to call the function - we need to do it twice to get the payoff for changing from A to B and from B to A. This is repeated during each iteration of the loop, so each time we can pass different values into it. To get the payoffs for switching from A to B and from B to A we have to add the calls to 'payoff' at the beginning of our loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'Ta' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-5-7c7966597477>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;32mwhile\u001b[0m \u001b[0mt\u001b[0m \u001b[1;33m<\u001b[0m \u001b[0mMAX_TIME\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m     \u001b[0mvariationBA\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpayoff\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mTa\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mTb\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      4\u001b[0m     \u001b[0mvariationAB\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpayoff\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mB\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mTb\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mTa\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'Ta' is not defined"
     ]
    }
   ],
   "source": [
    "while t < MAX_TIME:\n",
    "    \n",
    "    variationBA = payoff(A, Ta, Tb)\n",
    "    variationAB = payoff(B, Tb, Ta)\n",
    "    \n",
    "    A = A + variation\n",
    "    B = B - variation\n",
    "    \n",
    "    # advance time to next iteration\n",
    "    t = t + 1    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is, we pass the number of believers and the attractiveness of the traits. The order in which we pass the input variables into a function is important. In the B to A transmission, B becomes 'believers', while Ta becomes 'Tx' and Tb becomes 'Ty'. In the second line showing the A to B transmission, A becomes '_believers', Tb becomes 'Tx' and Ta becomes 'Ty'."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The obvious problem after this change is that we have not defined the 'attractiveness' of each trait. To do so add their definitions at the beginning of the script around other definitions (N, A, B, MAX_TIME, etc). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Ta = 1.0            # initial attractiveness of option A\n",
    "Tb = 2.0            # initial attractiveness of option B\n",
    "alpha = 0.1         # strength of the transmission process"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now calculate the difference between the perceived payoffs during this time step. To do so, we need to first see which one did better (A or B).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'variationBA' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-7-5ac21a8bcc38>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mdifference\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mvariationBA\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mvariationAB\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'variationBA' is not defined"
     ]
    }
   ],
   "source": [
    "    difference = variationBA - variationAB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, if the difference between the two is negative then we know that during time step B is more interesting than A. On the contrary, if it is positive then A seems better than B. What is left is to see how many people moved based on this difference between payoffs. We can express it in the main while loop, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'difference' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-8-e00b0e1b32ba>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;31m# B -> A\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[1;32mif\u001b[0m \u001b[0mdifference\u001b[0m \u001b[1;33m>\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      3\u001b[0m     \u001b[0mvariation\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdifference\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mB\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;31m# A -> B\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'difference' is not defined"
     ]
    }
   ],
   "source": [
    "    # B -> A\n",
    "    if difference > 0:\n",
    "        variation = difference*B\n",
    "    # A -> B        \n",
    "    else:\n",
    "        variation = difference*A\n",
    "    # update the population    \n",
    "    A = A + variation\n",
    "    B = B - variation    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use an additional term to have a control over how strong the transmission from A to B and back is - we will call it alpha (α). This parameter will multiply change before we modify populations A and B:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'variation' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-9-a661a3be48c1>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mvariation\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0malpha\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mvariation\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'variation' is not defined"
     ]
    }
   ],
   "source": [
    "    variation = alpha*variation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we need to add it to the rest of the parameters of our model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# temporal dimension\n",
    "MAX_TIME = 100\n",
    "t = 0               # initial time\n",
    "\n",
    "# init populations\n",
    "N = 100          # population size\n",
    "A = 65           # initial population of believers A\n",
    "B = N-A             # initial population of believers B\n",
    "\n",
    "# additional params\n",
    "Ta = 1.0            # initial attractiveness of option A\n",
    "Tb = 2.0            # initial attractiveness of option B\n",
    "alpha = 0.1         # strength of the transmission process"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should bear in mind that the main loop code should be located after the definition of the transmission function and the initialization of variables (because they are used here). After all the edits you have done it should look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "while t < MAX_TIME:\n",
    "    # calculate the payoff for change of believers A and B in the current time step       \n",
    "    variationBA = payoff(A, Ta, Tb)      \n",
    "    variationAB = payoff(B, Tb, Ta)     \n",
    "    difference = variationBA - variationAB\n",
    "             \n",
    "    # B -> A   \n",
    "    if difference > 0:\n",
    "        variation = difference*B\n",
    "    # A -> B        \n",
    "    else:\n",
    "        variation = difference*A\n",
    "        \n",
    "    # control the pace of change with alpha\n",
    "    variation = alpha*variation  \n",
    "    \n",
    "    # update the population    \n",
    "    A = A + variation\n",
    "    B = B - variation\n",
    "    \n",
    "    # advance time to next iteration\n",
    "    t = t + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, we have all the elements ready now and if you run the code the computer will churn all the numbers somewhere in the background. However, it produces no output so we have no idea what is actually happening. Let's solve this by visualising the flow of believers from one option to another. \n",
    "First, we will create two empty lists. Second, we will add there the initial populations. Then, at each timestep, we will add the current number of believers to these lists and finally, at the end of the simulation run we will plot them to see how they changed over time. \n",
    "Start with creating two empty lists. Add the following code right after all the variable definitions at the beginning of the code before the while loop: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# initialise the list used for plotting\n",
    "believersA = []\n",
    "believersB = []\n",
    "\n",
    "# add the initial populations\n",
    "believersA.append(A) \n",
    "believersB.append(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The whole initialisation/definition block at the beginning of your code should look like this: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# initialisation \n",
    "MAX_TIME = 100\n",
    "t = 0               # initial time\n",
    "N = 100           # population size\n",
    "A = 65              # initial proportion of believers A\n",
    "B = N-A             # initial proportion of believers B\n",
    "\n",
    "Ta = 1.0            # initial attractiveness of option A\n",
    "Tb = 2.0            # initial attractiveness of option B\n",
    "alpha = 0.1         # strength of the transmission process\n",
    "\n",
    "# initialise the list used for plotting\n",
    "believersA = []    \n",
    "believersB = []\n",
    "\n",
    "# add the initial populations\n",
    "believersA.append(A) \n",
    "believersB.append(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We just added the initial number of believers to their respective lists. However, we also need to do this at the end of each time step. Add the following code at the end of the while loop - remember to align the indents with the previous line! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "    believersA.append(A)\n",
    "    believersB.append(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The whole while-loop block should now look like this: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "while t < MAX_TIME: \n",
    "    # calculate the payoff for change of believers A and B in the current time step       \n",
    "    variationBA = payoff(A, Ta, Tb)      \n",
    "    variationAB = payoff(B, Tb, Ta)     \n",
    "    difference = variationBA - variationAB\n",
    "    \n",
    "        # B -> A   \n",
    "    if difference > 0:\n",
    "        variation = difference*B\n",
    "    # A -> B        \n",
    "    else:\n",
    "        variation = difference*A\n",
    "        \n",
    "    # control the pace of change with alpha\n",
    "    variation = alpha*variation  \n",
    "    \n",
    "    # update the population    \n",
    "    A = A + variation\n",
    "    B = B - variation   \n",
    "    \n",
    "    # save the values to a list for plotting    \n",
    "    believersA.append(A)\n",
    "    believersB.append(B)\n",
    "    \n",
    "    # advance time to next iteration\n",
    "    t = t + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's plot the results. First, we will import Python's plotting library, Matplotlib and use a predifined plotting style. Add these two lines at the beginning of your script:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt      # plotting library\n",
    "plt.style.use('ggplot')              # makes the graphs look pretty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's plot! Plotting in Python is as easy as saying 'please plot this data for me'. Type these two lines at the very end of your code. We only want to plot the results once the simulation has finished so make sure this are not inside the while loop - that is, ensure this block of code is not indented. Run the code!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f6c5c0ec7b8>]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEECAYAAADOJIhPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX5x/HPmewhk5UsZGXfRBFERVABETdcaznWpdXa\nWn9arVtdwFpxqQJarXWpVq21VVuOdnGtu3GtIrgDguwkIQnZ9/We3x8TW7QsIZPJneV5v155mZnM\nZJ7HCfc7595zz1XWWoQQQkQmj9sFCCGEcI+EgBBCRDAJASGEiGASAkIIEcEkBIQQIoJJCAghRASL\n3tMDtNYPA8cDFcaY/XruSwOWAkXAJkAbY+p7fjYfOBfoAi4xxrwcmNKFEEL4qzcjgUeAo7913zXA\nq8aYMcDrwHwArfV4QAPjgGOB+7TWqjeFaK1n9rLmkCT9hTbpL3SFc2/gf397DAFjzDtA7bfuPgl4\ntOf7R4GTe74/EfirMabLGLMJ+Ao4qJe1zOzl40LVTLcLCLCZbhcQYDPdLiDAZrpdQADNdLuAAJvp\nz5P7ekwgyxhTAWCMKQeyeu7PA7bu8LjSnvuEEEIEof46MCxrTwghRAja44HhXajQWmcbYyq01jlA\nZc/9pUDBDo/L77nvf/Tsx5r59W1jzPXA9X2sJ+gZY0D6C1nSX+gK597A15/Wese7io0xxb19fm9D\nQPV8fe0Z4BxgMXA28PQO9z+utb4T326gkcCyXRReDOxY6PVlZWW9LCf0eL1eGhsb3S4jYKS/0BbO\n/YVzbwC5ubkYYxb29fm9mSL6BL5P7Bla6y34EnUR8KTW+lxgM74ZQRhjVmmtDbAK6AQuNMbIriIh\nhAhSKoiWkrYyEghd0l9oC+f+wrk38I0E+Oaemr0iZwwLIUQEkxAQQogI1tfZQUIIEdSstdS2dfNl\nXT2ry2rZUt9OUWocJ4/LcLu0oCIhIIQIedZayps6WV/TxvqaNjbUtLGxth0HGJGRSIE3mnGZiYwd\nnOB2qUFHQkAIEXIa27tZU9XKmqpW1la3sa66lbhoDyPT4xmZHs8JY9MZlhZHekI0ycnJYX1g2F8S\nAkKIoFfR1MHKylZWb29hVWUr1S1djMqIZ8zgBI4bncrojCGkJcjmrC/k/5oQIuhUt3TyWXkLn1U0\n80VFC+3dlglZiYzPSuDYUWkUpcYR5enzrEixAwkBIYTrOrodVla28lFZEx9va6a2tYt9swcxMSeR\nU8dnkJcci1Ky0Q8ECQEhhCtqWrtYXtrEh6VNfF7eQlFqHJNzB/GzqUMYkR4vn/QHiISAEGLAbGvs\n4N0tjXywtZHSxg4mDxnEoYVeLp46hOS4KLfLi0gSAkKIgNrW2MHbmxt4b0sjta1dHFLg5cyJmeyT\nlUhMlHzad5uEgBCi39W1dfH2pgbe3NRAZXMn0wq8nHdANmMzE2Q3T5CREBBC9IvObsuHpY28vqGe\nVZWtHJifxBn7DWZiziDZ8AcxCQEhhF9K6tt5eV0dxRsbKEiJZfaIVK6YnkdCjCxNFgokBIQQe62z\n2/L+1kZeXFdHSX07s4ensPjoIoZ4Y90uTewlCQEhRK/VtXbx0ro6/vVVHXneGI4dncbB+V45wBvC\nJASEEHu0pa6dp7+s4d9bG5le6GXhrHyGpsW7XZboBxICQoidstayqrKVv62qZkNNG8eNTuP+E4aT\nHC+bjXAi76YQ4hustawoa+bJL6qpb+/iO+MzuObwPGKj5EBvOJIQEEIAvo3/stImln5eTZdjmbdP\nBtMKvTK9M8xJCAgR4ay1LNtSz0Pvb6HbWk7bdzAH5yfhkQXbIoKEgBARbFVlC3/6ZDstXZbTJqRz\nSIFXNv4RRkJAiAi0pa6dP32ynU21bZwxMZPj982jpbnJ7bKECyQEhIggda1dPPFZFe9vbeTUfTK4\n6rBcYqM8st8/gkkICBEBOrsdnvmyln+srmHWsGTuPWE4Xlm6WSAhIERYs9ayvLSZhz+qID85jiVH\nFZGbLEs7iP+SEBAiTJU3dvDg8gq2NXXykynZTM5NcrskEYQkBIQIM53dDn9fVcOza2o5eVw61xye\nLmv7iF2SEBAijKysaOG+ZeUM8cZyxzFDyUqKcbskEeQkBIQIA00d3Tz6cSUrSps5b0o2UwuSUDLf\nX/SChIAQIW5ZSSP3L6tgSl4Sdx8/jEGxMutH9J6EgBAhqrG9mweXV7CmqpXLpg9h3+xBbpckQpCE\ngBAhaHlpE/d+UM60Qi93zR1GfLSs8Cn6RkJAiBDS0tnNwysq+ay8hcvl07/oBxICQoSI1dtb+M17\n25iQnchdc4eSGCP7/oX/JASECHLdjmXpF1W89FUdFxyUw9QCr9sliTDiVwhorecDZwHdwOfAD4FB\nwFKgCNgEaGNMvX9lChGZKpo6uOPdbcRHK+48bhjpCfK5TfSvPh9N0loXAecBk4wx++ELlNOBa4BX\njTFjgNeB+f1RqBCR5t3NDVz54mamFiRx/REFEgAiIPz5q2oAOoBBWmsHSABK8W30Z/Q85lGgGF8w\nCCF6oaPb4Q8rKvl4WzPXzcpnVEaC2yWJMNbnkYAxphb4NbAF38a/3hjzKpBtjKnoeUw5kNUfhQoR\nCcoaOrjqpc00tHdzx7FDJQBEwPV5JKC1Hg5chm/ffz3wpNb6TMB+66Hfvv3182cCM7++bYzB6w3f\nA16xsbHSXwgbiP7eXF/Db97ewjkH5nHi+MwBXfYhnN+/cO7ta1rrhTvcLDbGFPf2ucranW6je/Oi\nGphjjDmv5/b3ganAEcBMY0yF1joHeMMYM64Xv9KWlZX1qZZQ4PV6aWxsdLuMgJH++q6z2/Lox5Us\nK23iqkPzGJkRH5DX2Z1wfv/CuTeA3NxcgD5/YvDnmMAa4DqtdTzQDswGPgSagHOAxcDZwNN+vIYQ\nYa26pZPb3iljUIyHO44ZSpJc7UsMMH+OCXwK/AlYAXyKL4l+j2/jP0drvQZfMCzqhzqFCDsrK1q4\n4sXNTB4yiGtn5ksACFf0eXdQAMjuoBAm/fWetZbn19Zivqjm0kOGBMUVv8L5/Qvn3sDd3UFCiL3U\n0e3wu2XlrK9pZ/FRRQzxyvV+hbskBIQYINUtndz6VilZg2JYcnSRrPwpgoKEgBADYE1VK4veKmXu\n6DRO3SddrvolgoaEgBAB9vqGev74USUXTc3hoPzwnq8uQo+EgBAB0u1YHvt0O+9taeTmOYUUpsS5\nXZIQ/0NCQIgAaOns5o53t9Ha2c1tRxeRHC//1ERwkr9MIfrZ9uZObi4uYVRGPFcflkdMlOz/F8FL\nQkCIfrS2qpVb3yrl5HHpnDg2TQ4Ai6AnISBEP3l3SwMPLKuQA8AipEgICOEnay1/X1XD82trWXhE\nAcPTB34BOCH6SkJACD90OZb7l5WzrqaN244uIiMxxu2ShNgrEgJC9FFzRzdL3i4l2qO4dU4RCTFy\nBrAIPRICQvTB9uZObnqjhPFZCZw3JZsojxwAFqFJQkCIvbShpo2b3yzhxLFpnDRWloAQoU1CQIi9\n8FFZE795bxvnH5TN9MJkt8sRwm8SAkL00qvr6/jzJ9uZf3ge47IS3S5HiH4hISDEHlhr+eOHpby8\npppfzSkkP1nWABLhQ0JAiN3ociz3fVBOSWMni48qIjVB/smI8CJ/0ULsQktnN0veLiNKwR0njqWr\nrcXtkoTodzKxWYidqG3t4hevbiFzUDQLZuSTECMXgRfhSUYCQnxLSUM7N75RwuzhKegJGTIFVIQ1\nCQEhdrB6ewuL3irl+/tncuSIVLfLESLgJASE6PH+1kbu/aCcSw8ZwgF5SW6XI8SAkBAQAnhhbS3m\ni2p+OSufURkJbpcjxICREBARzVrLY59W8d6WBhbNKSTHG+t2SUIMKAkBEbE6uy33frCN0oYOFh1V\nRIpcB1hEIPmrFxGppbObxW+VEhPl4eYjC4mLltnSIjJJCIiIU9PaxU1vbGVURgLnHyjLQIvIFlQh\n0H3Lz33fxMSivCmQnALJqZCchkpJg+Q0SEmD5FRUdFCVLkLE1vp2bnxjK3NGpDJPzgEQIrhCwHPa\nj33fdLRjmxqgoc73tXEtTn2t7/v6GmhsgEFJkJoOqRmo9MGQmgHpmaj0TEgfDGmDJSjEN6ysbGHx\n26WcMymLI4anuF2OEEEhqLaSasTY/36/m8dZp9sXCLU1UFeNra2C2mr4YgVOTRXUVPp+npwKGVmo\nwdkwOAcyc1BZQyArF5K88ikwgry7uYEHPqzgsum5TBoyyO1yhAgaQRUCvaU8Ub5P/qkZwKidBobt\n6oK6aqiqwFZVwPYK+Hw5zvZyqCwDC2TnorJzIScPNaQAcvJ990XLxcLDhbWWZ76s5ekva7hhdgHD\n0uLdLkmIoBKSIdAbKjoaBmfD4Oydh0RjA1SWYctLoaIE5/03oXwrVG+HjCzIK0TlFqHyh0L+UN/v\n8cgMklDS7Vj+8FEln5U3s/ioIjIHSbgL8W1hGwJ7orzJ4E3+xi4oANvZCRUl2NItULoZ591XoWQj\ntDRD/jBU4XAoGoEqHAFDClBRsrpkMGrvcvj1u2W0dDrcelQRSbHyPgmxMxEbAruiYmJ8G/v8Yd+4\n3zY1wNaN2C0bYOXHOM8/6TtIXTgcNWw0avgYGD4WlZbhUuXia3VtXdxcXEJecixXHppHTJQc+xFi\nV/wKAa11CvAQMAFwgHOBtcBSoAjYBGhjTL1/ZbpPJSXDuImocRP/c59taYJN67Ab1+K8+xo8dh/E\nxqFGjodR41Gj9vGNFmQ30oDZWt/OTcUlzBiazBn7DZaD/0Lsgb8jgbuAF4wx87TW0cAgYAHwqjFm\nidb6amA+cI2frxOUVGISjN8fNX5/wHcQkopS7LrVsG4Vzsv/hNYWGD2B9v2mYIeP9R2Elg1TQHxW\n3szt75Zx9v6ZzJZloIXolT6HgNY6GTjMGHMOgDGmC6jXWp8EzOh52KNAMWEaAt+mlIKcfFROPhw6\nBwBbsx279gu616/GeeYJcBzUuP1hwmTU+P19Iwzht9c31PPHjyq54tBcJubIFFAhesufkcAwoEpr\n/QgwEVgOXApkG2MqAIwx5VrrLP/LDF0qPRM1dRaJc06kq6EBKrdhV32MXfYW9rH7fKGx3xTUfgdC\nwXAZJewlx1qe+LSKtzY38Ks5hRSkxLldkhAhxZ8QiAYmAz81xizXWt+J7xO//dbjvn07Yiml/ntu\nwqy5vplIX63Efr4c5/7F0N2FmngwatJUGD1BZh7tQXuXw2/f38b25i6WHF1EqqwCKsRe8+dfTQmw\n1RizvOf23/CFQIXWOtsYU6G1zgEqd/ZkrfVMYObXt40xeL1eP8oJbrGxsTvv7+DD4ODDsNbilG6m\n88N36PznYzjby4k+YBoxU2cSPWFy0C+Bscv+AqS6uYNfvrGOIclx3HXyaGIDvAroQPc30MK5v3Du\n7Wta64U73Cw2xhT39rnK2r5/UNdavwmcZ4xZq7W+Hkjs+VGNMWZxz4HhNGNMb44J2LKysj7XEuy8\nXi+NjY29frytrsR+9G/s8negchtq8jTUwYfDyPFBOdtob/vzx4aaNm55s4QjR6Zy2gAtAjeQ/bkh\nnPsL594AcnNzYfcr7eyWvx8vfwY8rrWOATYAPwSiAKO1PhfYDGg/XyMiqYws1JyTYM5J2KoK7Ifv\n4DzxALS2oA6egZp2hO8AdIT599ZG7vugnJ9MyeawoXJQXQh/+TUS6GcyEugFW7IR++83sO8X+xbE\nm34k6sDDUPHuXhc30J+2rLU8ubKaF7+qY/7heQN+HeBw/zQZzv2Fc2/g/khADDCVPww1bxj2lB/A\nyo9w3nkF+9QfUVOmow4/BlU0wu0S+117l8Pd72+jvKmT244uIiNR1gASor9ICIQoFR0NEw8iauJB\n2Lpq7Luv4dx3C6SkoWbN9YVCTOhfNH17cye3vFlCYUocv5LLQArR7yQEwoBKzUDN1dhjT4XPPsR5\n4wXsU4+gZhyDmnksKjnN7RL7ZGVFC7e9U8op4zM4cWyanEMhRABICIQR5YmC/acStf9UbNkW7GvP\n4Vx3IWrSIaijT/FdMyEEWGt5bk0tT66s5rJpchEYIQJJQiBMqdxC1PcvxJ5yFrb4BZzbFsCw0XiO\n/S5q5Di3y9ul9i6H+5aVs7munSVHFZHjDf1dWkIEMwmBMKeSklHHfw971CnY917DeejXMDgbz1wN\nY/cLql0s2xo7WPRWKYWpcSw+qkj2/wsxACQEIoSKjUPNPA576FHYZW/hPHE/JCXjOeF03xLZLofB\nspJG7nm/nNP2Hcxxo1Ndr0eISCEhEGFUdDRq2hHYqTN8J6D95QFISsFzylmo0RMGvJ5ux/LYp9t5\nc1MDC2bkMzbT3fMdhIg0EgIRSnmiUAfPwB54KPaDt3AeuQuyc/Gc8oMBO9eguqWT298pIzbaw53H\nDiVFFoATYsDJv7oIpzxRqENm+cLg7Vdw7r4JNWZf1MlnojJzAva6H5U18dt/b+OY0WnM2yeDKI/s\n/hHCDRICAgAVHYOadRz2kFnYl/+Jc8sVqGlHouZqVGL/TdHsciyPf7qdNzc2cMWhueybLdM/hXCT\nTL8Q36DiE/CceDqehfdAazPOdRfgFL+A7e72+3eXN3aw4JXNbKpt587jhkoACBEEZCQgdkqlpKF+\ncBF260acpQ9h33wRz+k/6fPB4+KN9Ty8opJ5EzI4fkwaHpn9I0RQkBAQu6UKhuG54mbs8ndxHr4D\nNWIcSp+LSs3o1fObOrp54MMK1te0ccMRBQxPjw9wxUKIvSG7g8QeKaXwHHgonhvvg8HZODf8DOe1\nZ/e4i+iz8mYufX4jSbG+2T8SAEIEHxkJiF5TcfGo7/wAe8gsnMfvx773Gp4fXPw/U0rbuxz+/Ml2\n3t3SyEUH53BAXpJLFQsh9kRGAmKvqSEFeK64GTX7BJy7FuI8+QdsWysAa6paufSFTdS3dXPX3GES\nAEIEORkJiD5RSqGmzcbuOwW79CGqrv4/Hp/xf7zVEMdPDsxmeqFc+lGIUCAhIPyivCmsPOEC7nl7\nMyPXrOaulGpSMs90uywhRC/J7iDRZ03t3dz9/jbueK+MC2eN4opzjybZacO54WfYLz9zuzwhRC/I\nSEDsNWstb29u5A8fVTKtIIl7jh9GdnoqjY2NqLMvxn6+AufhO1GTpqJOPQcVF+d2yUKIXZCRgNgr\nJQ3t/PL1rfx9VTXXHJbHTw7MITEm6huPUfsegGfh3dDchHPzpdhNX7lUrRBiT2QkIHqltdPBfFHF\nK+vr0RMymDs6bbeLvqlBSajzrsBZ9hbOb29EHXE86rjv+i6BKYQIGhICYrestby1qYFHP97OvtmJ\n3HXcUDISY3r9fM9Bh2NHjsP5w2+wqz/B86PLUemZAaxYCLE3JATELq2pauXhFZV0OQ5XHprLuKzE\nPv0elZ6J5/IbsS/9A+fmy/GcdQFq8rR+rlYI0RcSAuJ/VDR18NinVXxR0cJZEwcza3iK3wu+KU8U\n6tjvYsfsi/Pg7agvP0PNOxcVIxeSF8JNcmBY/EdDezcPr6jgin9tItcbw30nDGf2iNR+XfFTDR+D\n57o7sQ11OLdeiS0v7bffLYTYexICgpbObv76eRUXPruBjm7L3ccP5/T9MkmICcyfh0pMwnP+1ajD\nj8ZZfDV2+TsBeR0hxJ7J7qAI1tbl8MLaWv65uoZJOYO47egihngHZveMUgo18zjssDE4DyxGrV3p\nW6I6uvcHnYUQ/pORQARq7XT4+6pqzn96Peuq27hpdiGXTc8dsADYkSoagecXd2Brq3CWzMfWVA14\nDUJEMhkJRJCm9m6eW1vLC2tqmZCdyI2zCylKdf9sXpWYhOeC+diX/o5zyxW+aaTjJrpdlhARQUIg\nAmxv7uTZL2t4fUM9B+Z7ueWoQvKT3d/470h5PL7ZQ0NH+a5gNudk1FEno+QylEIElIRAGFtX3caz\nX9awoqyJI4ancOdxw8gcFNz73NW4iXjm347zu1th8zo4+2JUnFyRTIhAkRAIM12O5YOSRp77spbK\n5k7mjknjvAOzSYoNneUaVEYmnqtuxT5+P86iq/BcuACVmeN2WUKEJQmBMFHd0skr6+t5+as6spNi\nmDsmjUMKvLtd3yeYqdg4OOdn8PrzviA47+eosfu5XZYQYUdCIIR1O5aPtzXz8ro6vqhs4dDCZH45\nK5+haeGx+0QphZp9PDa3wHeW8VyNmjVXjhMI0Y/8DgGttQdYDpQYY07UWqcBS4EiYBOgjTH1/r6O\n+K+ShnZeX1/PGxsbyEiM5qiRqVw2LTdgJ3e5TY2biOeaJTj33AylW+D0n6Ci5fOLEP2hP7YalwCr\ndrh9DfCqMWYM8Dowvx9eI+LVtXXx3Joafv7iJq59ZQvdFm44ooDbjxnKUSNTwzYAvqYyc/BcswRb\nV41z5y+xjQ1ulyREWPDr45TWOh84DvgVcHnP3ScBM3q+fxQoxhcMYi81tXfzfkkjb29u5KuqVg7I\nS+KM/QYzMWdQyO7r94dKSMTz0wXYv/8Z59af47n4OtSQArfLEiKk+TumvhO4EkjZ4b5sY0wFgDGm\nXGud5edrRJS6ti6WlTTx3pZG1lS1sl9OInNGpDD/8Dzio8P7035vKE8U6rvn4AwpwLltAZ4fX44a\nP8ntsoQIWX0OAa31XKDCGPOJ1nrmbh5q+/oakcBay9aGDpaXNPFBSRNb69vZf8ggjhyRwtWH5YX9\nbp6+8kyfjc3MxnlgCeqE0/HMPNbtkoQIScravm2jtda3AGcBXUAC4AX+AUwBZhpjKrTWOcAbxphx\nO3n+TGDm17eNMdc3Njb2qZZQEBsbS0dHBwCtnd18UtbIsi31fLClDsfCwYUpTB+axv55XmKjQm/D\nv2N/A6m7vJTmJfOJmTSV+DPPD9jlK93qb6CEc3/h3BuA1+tFa33DDncVG2OKe/v8PofAjrTWM4Ar\nemYHLQGqjTGLtdZXA2nGmN4cE7BlZWV+1xKMuh3LtjYP723YzqflLayrbmNURjyThgzigNxBFKXG\nhfy0R6/Xi1shbpubfGcYJyTi+fEVATnD2M3+BkI49xfOvQHk5uYC9HkDEoh5dosAo7U+F9gM6AC8\nRlDr7Lasr2ljZWULKytbWL29lZzkOCZkxnPKuHT2yUqU3Tz9SA1KwnPpQuyf7/MdJ7j4OlRKmttl\nCRES+mUk0E9CdiRQ39bFl1WtrNneypqqVtbVtJHrjWVcViL7ZiWyT1YCeZlpYf1pJBg+bVlrsc8b\n7Duv+IIgr6jffncw9BdI4dxfOPcGwTkSCGvNHd2sr2ljfU0bX1X7vpo7uhk1OIGxg+M5dZ8MxmYm\nkBgTOmv1hAulFOr403AGZ+P8+he+pSZkSWohdktCYBestVQ2d7K5rp1Nte1sqG1nY20bdW1dDE2N\nZ0RGPFPykjhj4mByvbH9eh1e4R/P1JnYtMG+K5Z994d4ph3hdklCBK2ID4Fux7K9uZOShg5KGtrZ\nWt/BljrffxNiPAxNjaMoNY5phV7O7NngR+KJWqFGjZmA58pbcO66AaemEjX3tJA/+C5EIERECHQ7\nlprWLrY1dlDe1Mm2xg7KGjsoa/Dd9sZFUZAcS15KHCPT4zlieAqFKXF442SXTihTQwrwzL8N5+6b\noHo7nHmBrDkkxLeExb+Izm6H6pYuqlq6qGzupKq5k8odvrY3d5ESF0V2UgxDvLEM8cZwWFEyud5Y\ncpNj5UzcMKZS0vD8/Fc4v78Ne+/NeM6/ChWf6HZZQgSNoA6Bjm6HutZu6tq6qG3rorbV91XT2kVN\nSxfVPf9t7nRIT4hmcGI0mYNiGJwYzciMeKYVeskaFENWUkxInoAl+oeKT8Dz02uxj/+uZwrpL1Gp\n6W6XJURQCKoQ+NnzGwFo73Kob+um03FIiY8mLT6a1Pgo0hKiSUuIZnhaPFPyoslIiCEjMZqU+Cg5\nMCt2S0VFwfd/Cs8b30VqLlmIGpLvdllCuC6oQuDyaUMAiInykBofRWKMRw7miX7znymkaYNxbpuP\n58L5qJHj3S5LCFcFVQiEyxWxRHDzTJ+NTUnDue9WPGddiJp8iNslCeEa2VEuIpKaMBnPJQtx/vIA\nzuvPuV2OEK6REBARSxWNwHPVIuwbz+P87VGs47hdkhADTkJARDSVmYPn6sXYr1ZiH/kNtqvT7ZKE\nGFASAiLiqaRkPJffhG1rxfntjdjWFrdLEmLASAgIAajYODwXXIPKzsVZMh9bV+12SUIMCAkBIXoo\nTxTqjP9DTZmOs+hq7LatbpckRMBJCAixA6UUnrkadeIZOLctwK5d6XZJQgSUhIAQO+GZdgSeH1+O\nc/8i7Ip33S5HiICREBBiF9T4SXguvQHnrw/R9sJTbpcjREBICAixG6pwOJ5rFtPx2nM4Sx+ScwlE\n2JEQEGIPVEYWSTfcjd2yHuf3S7CdHW6XJES/kRAQohc8SV48l96A8kTh3HEdtqnB7ZKE6BcSAkL0\nkoqJRf34CtTI8Ti3XoWt3OZ2SUL4TUJAiL2gPB48p56NmnMSzpJrsOu/dLskIfwiISBEH3hmHovn\n7Itx7rkZu+I9t8sRos8kBIToI7XvFDyX3YDz1wdxXvo71lq3SxJir0kICOEHVTgCz/wl2PeLsY/d\nh+3qcrskIfaKhIAQflLpmXiuXoStrca5+0ZsS7PbJQnRaxICQvQDFZ+I56fXorLzcBZdhd1e7nZJ\nQvSKhIAQ/URFReE543zUrONwFl+N/WqV2yUJsUcSAkL0M8+suXjOuQTnd7fivPea2+UIsVsSAkIE\ngJowGc+Vt2CfW4rz1CNYp9vtkoTYKQkBIQJEDSnAs+B27KZ1OPfeIpetFEFJQkCIAFJJyb41h9Iy\ncG69EltR5nZJQnyDhIAQAaaio/GcdSHqiON9B4y/+MjtkoT4j2i3CxAiUnhmHovNLcT5/W2oI09A\nHf0dlFJulyUinIwEhBhAavQ+eBbchl3xHvb3t2Hb29wuSUS4Po8EtNb5wJ+AbMABHjTG/FZrnQYs\nBYqATYA2xtT3Q61ChAWVnonnqluxj/8O55af47lgPionz+2yRITyZyTQBVxujNkHOAT4qdZ6LHAN\n8KoxZgz46fEuAAALpklEQVTwOjDf/zKFCC8qJhZ19s98xwmWXIP96N9ulyQiVJ9HAsaYcqC85/sm\nrfVqIB84CZjR87BHgWJ8wSCE2IFSCjXjGGzhcJz7F6M2rEGd8n1UVJTbpYkI0i/HBLTWQ4H9gfeB\nbGNMBfwnKLL64zWECFdq2Gg8v7gTW7IR59fXYuuq3S5JRBC/ZwdprZOAp4BLekYE315UfaeLrGut\nZwIzv75tjMHr9fpbTtCKjY2V/kJYwPvzerHX3k77Px+n/VdXkHDhfGL2mxK41/uWcH7/wrm3r2mt\nF+5ws9gYU9zb5yp/LoShtY4GngP+ZYy5q+e+1cBMY0yF1joHeMMYM64Xv86WlYXviTRer5fGxka3\nywgY6a//2NWf4vzhTtS02agTzxiQ3UPh/P6Fc28Aubm5AH2ea+zv7qA/AKu+DoAezwDn9Hx/NvC0\nn68hRERR4ybiue43vuUmbl+Ara50uyQRxvo8EtBaTwfeAj7Ht8vHAguAZYABCoDN+KaI1vXiV8pI\nIIRJf/3POg72lX9iX/oH6vTz8Rx4aMBeK5zfv3DuDfwfCfi1O6ifSQiEMOkvcOymr3AevB01ajzq\ne+eh4hP7/TXC+f0L597A/d1BQogAU0NH4bnuTlAenBsvxa6Ti9WI/iMhIEQIUPGJeM6+GM+8c3F+\ntwjn73/Cdna6XZYIAxICQoQQNWkqnuvvwm7bivOry7FbNrhdkghxEgJChBiVnIbnwgWoo7+D85vr\ncZ55AtslowLRNxICQoQgpRSeQ2b5ppJu2YBz8+XYjV+5XZYIQRICQoQwlZaB56fXoo79Ls49N+Es\nfRjb1up2WSKESAgIEeKUUngOnoFn4T3Q3ICz8GLsZx+6XZYIEXJlMSHChPImo869DLvqE5zHfwfv\nvILntPNQGZlulyaCmIwEhAgzavz+eBbejSoYjnPzpTj/ekqmk4pdkhAQIgypmFg8J3wPz/zbseu/\nxFl4EfbTZQTRCgEiSMjuICHCmMoaQtRFv8B+sQJn6UPw+vN49LmovCK3SxNBQkJAiAigJhyAZ+x+\n2DdfxPn1L1CTpvqWqU5Jc7s04TLZHSREhFDRMXhmn4DnpvsgNh7n+otwnn4c29ridmnCRRICQkQY\nNciL57Qf+Ralq67EufZ82p4z2I52t0sTLpDdQUJEKJWR5ZtSWrqZ7hcMznNLUcd+F3XYUaiYWLfL\nEwNERgJCRDiVV8Sgy2/Ec9EvsF98hLPgfJzXnpWRQYSQkYAQAgBVNJKon/0Su3kdzrN/xf7rKdTs\nE1AzjkUlDnK7PBEgMhIQQnyDKhpJ1EW/wHPpDVC2BWfBT3CefARbs93t0kQAyEhACLFTKn8o6keX\nY6sqsK89i3PDJah9D0AdeSJq6Ci3yxP9REJACLFbanA26rQfY0/4Hvbtl3HuXwwpaagjjkdNnoaK\niXG7ROEHCQEhRK+oxCTU0d/BHnkSfLoM543nsUsfQk0/0jejKGuI2yWKPpAQEELsFRUVBZMPIWry\nIdjyUuzbL+HceiXkFfnCYNJUVGyc22WKXpIQEEL0mcrJQ807F3vy9+HTD3DeeQX7xAOoA6ahps6C\nUeNRSrldptgNCQEhhN9UTAxMOZSoKYdia6qwy97Eeew+6GhHHXgY6qDDIX+oBEIQkhAQQvQrlT4Y\ndcyp2KO/AyWbsMvewrn3VxAd4xshHDAdCoZJIAQJCQEhREAopXwb+4Jh2O/8ADatw654F+f+ReA4\nvmMH+x8MI8ahomVT5Bb5Py+ECDilFAwbhRo2Cnvq2VC6GfvJ+zhPPgLbt6HGT4J9D0CNn4RKTXe7\n3IgiISCEGFBKKd/xgfyhcPz3sHU12C9WYD/7ELv0IUjPQo3fHzVuou/Acly82yWHNQkBIYSrVGo6\n6tA5cOgcbHc3bFyLXf0pzgsGtmyAguGoMRNQoyfA8DGo+AS3Sw4rEgJCiKChoqJg5DjUyHFwwvew\n7W2wfjV2zUqcZ/8CWzdCTr7v5yPGooaPhfTBcpDZDxICQoigpeLiYfwk3zEDwHZ2+A4wr1uNXfY2\n9q8PglIwdBRq6EhU0UgoGoFKlstm9paEgBAiZKiYWN9xglHjAbDWQk0VbFqL3fgVzitPw+b1EBsL\n+b6ZSR2jxmHTsyE7V2Yh7YT8HxFChCylFGRkQkam7/wDeoKhqsJ3jsLWjXT+uxhn83qorULNPgHP\nqWe7XHVwkRAQQoQVpRRk5kBmDmrSVAZ5vTQ2NmLb26G91e3ygo6EgBAiIqi4OIiThe2+LWAhoLU+\nBvgNvquXPWyMWRyo1xJCCNE3Abm8pNbaA9wDHA3sA5yutR4biNcSQgjRd4G6xvBBwFfGmM3GmE7g\nr8BJAXotIYQQfRSoEMgDtu5wu6TnPiGEEEEkUCEghBAiBATqwHApULjD7fye+/5Daz0TmPn1bWMM\nubm5ASonOHi9XrdLCCjpL7SFc3/h3BuA1nrhDjeLjTHFvX6ytbbfv+bNmxc1b968dfPmzSuaN29e\n7Lx58z6ZN2/euD08Z2EgagmWL+kvtL+kv9D9Cufe+qO/gOwOMsZ0AxcBLwMrgb8aY1YH4rWEEEL0\nXcDOEzDGvAiMCdTvF0II4b9gOjBc7HYBAVbsdgEBVux2AQFW7HYBAVbsdgEBVOx2AQFW7M+TlbW2\nn+oQQggRaoJpJCCEEGKASQgIIUQEC4pVRMNpsTmtdT7wJyAbcIAHjTG/1VqnAUuBImAToI0x9a4V\n6qee9aGWAyXGmBPDqT+tdQrwEDAB33t4LrCW8OlvPnAW0A18DvwQGESI9qe1fhg4HqgwxuzXc98u\n/x57+j8X6AIuMca87EbdvbWL/pYAJwDtwHrgh8aYhp6f7VV/ro8EwnCxuS7gcmPMPsAhwE97+rkG\neNUYMwZ4HZjvYo394RJg1Q63w6m/u4AXjDHjgInAl4RJf1rrIuA8YFLPBiUaOJ3Q7u8RfNuPHe20\nH631eEAD44Bjgfu01sF+geKd9fcysI8xZn/gK/zoz/UQIMwWmzPGlBtjPun5vglYje+M6ZOAR3se\n9ihwsjsV+q9ntHMcvk/LXwuL/rTWycBhxphHAIwxXT2fIMOiP6AB6AAGaa2jgQR8Z/OHbH/GmHeA\n2m/dvat+TsR33lKXMWYTvg3oQQNRZ1/trD9jzKvGGKfn5vv4tjHQh/6CIQTCdrE5rfVQYH98b1K2\nMaYCfEEBZLlYmr/uBK4EdpxaFi79DQOqtNaPaK0/0lr/XmudSJj0Z4ypBX4NbMG38a83xrxKmPS3\ng6xd9PPt7U0pob+9ORd4oef7ve4vGEIgLGmtk4Cn8O2Ta+KbG0x2cjskaK3n4ts3+Qmwu2FmSPaH\nb/fIZOBeY8xkoBnfroVwef+GA5fh21eei29EcCZh0t9uhFs/AGitrwU6jTF/6evvCIYQ2ONic6Gm\nZ5j9FPBnY8zTPXdXaK2ze36eA1S6VZ+fpgMnaq03AH8BjtBa/xkoD5P+SoCtxpjlPbf/hi8UwuX9\nmwK8a4yp6Vne5R/ANMKnv6/tqp9SoGCHx4Xs9kZrfQ6+3bJn7HD3XvcXDCHwITBSa12ktY4Fvgc8\n43JN/voDsMoYc9cO9z0DnNPz/dnA099+UigwxiwwxhQaY4bje69eN8Z8H3iW8OivAtiqtR7dc9ds\nfOtfhcX7B6wBpmqt43sOGM7Gd4A/1PtTfHNkuqt+ngG+p7WO1VoPA0YCywaqSD98o7+eGZVXAica\nY9p3eNxe9xcUZwz3NHQX/50iusjlkvpMaz0deAvf1Dvb87UA3xth8KX0ZnxT1urcqrM/aK1nAFf0\nTBFNJ0z601pPxHfQOwbYgG8KZRTh09+V+DaQ3cDHwI8BLyHan9b6CXzL0mcAFcD1wD+BJ9lJPz1T\nKH8EdBIaU0R31t8CIBao7nnY+8aYC3sev1f9BUUICCGEcEcw7A4SQgjhEgkBIYSIYBICQggRwSQE\nhBAigkkICCFEBJMQEEKICCYhIIQQEUxCQAghItj/A0jLaCFGNaliAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6c5c0ec6d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "# plot the results    \n",
    "plt.plot(believersA)\n",
    "plt.plot(believersB)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see how over time one set of believers increases while the other decreases. This is not particularly surprising as the attractiveness Tb is higher than Ta. However, can you imagine a configuration where this is not enough to sway the population? Have a go at setting the variables to different initial values to see what combination can counteract this pattern. \n",
    "1. set the initial value of A to 5, 10, 25, 50, 75\n",
    "2. set the MAX_TIME to 1000,\n",
    "3. set the Ta and Tb to 1.0, 10.0, 0.1, 0.01, etc.,\n",
    "4. set alpha to 0.01, and 1.0\n",
    "\n",
    "You can try all sorts of configurations to see how quickly the population shifts from one option to another or what are the minimum values of each variable that prevent it. \n",
    "\n",
    "However, we can make the model more interesting if we allow the attractiveness of each option to change through time. To do so let's define a new function. Add the following line at the beginning of the while loop (remember indentation!). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'attractiveness' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-18-c7b0a9890c8d>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mTa\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mTb\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mattractiveness\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mTa\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mTb\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'attractiveness' is not defined"
     ]
    }
   ],
   "source": [
    "    Ta, Tb = attractiveness(Ta, Tb) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ta and Tb will then be modified based on some dynamics we want to model. Let's define the 'attractiveness' function. We have already done it once for the 'payoff' function so it should be a piece of cake. At each time step we will slightly modify the attractiveness of each trait using a kernel K that we will define. This can be expressed as: \n",
    "\n",
    "\\begin{equation*}\n",
    "T_{A, t+1} = {T_A} + {K_a}\n",
    "\\end{equation*} \n",
    "\n",
    "\\begin{equation*}\n",
    "T_{B, t+1} = {T_B} + {K_b}\n",
    "\\end{equation*} \n",
    "\n",
    "K can have several shapes such as:\n",
    "* Fixed traits with $K_a = K_b = 0$\n",
    "* A gaussian stochastic process such as $K = N (0, 1)$\n",
    "* A combination (e.g., $K_a = N (0, 1)$ and $K_b = 1/3)$\n",
    "\n",
    "Let's start with a simple case scenario, such as a) $T_a$ will increase each step by a fixed $K_a$ and b) $K_b$ is equal to zero (so $T_b$ will be fixed over the whole simulation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def attractiveness(Ta, Tb):\n",
    "\n",
    "    Ka = 0.1  \n",
    "    Kb = 0\n",
    "    \n",
    "    Ta = Ta + Ka\n",
    "    Tb = Tb + Kb\n",
    "    return Ta, Tb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we define the function and give it the input values. \n",
    "```python\n",
    "def attractiveness(Ta, Tb):\n",
    "```\n",
    "\n",
    "Then, we establish how much the attractiveness of each trait changes (i.e., define $K_a$ and $K_b$)\n",
    "```python \n",
    "Ka = 0.1         \n",
    "Kb = 0\n",
    "```\n",
    "\n",
    "And plug them into the equations:\n",
    "\n",
    "```python \n",
    "Ta = Ta + Ka\n",
    "Tb = Tb + Kb\n",
    "```\n",
    "\n",
    "Finally, we return the new values:\n",
    "```python \n",
    "return Ta, Tb\n",
    "```\n",
    "\n",
    "This is how the function is defined. The main loop will now look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "while t < MAX_TIME: \n",
    "    # update attractiveness\n",
    "    Ta, Tb = attractiveness(Ta, Tb)\n",
    "    # calculate the payoff for change of believers A and B in the current time step           \n",
    "    variationBA = payoff(A, Ta, Tb)      \n",
    "    variationAB = payoff(B, Tb, Ta)     \n",
    "    difference = variationBA - variationAB\n",
    "    \n",
    "        # B -> A   \n",
    "    if difference > 0:\n",
    "        variation = difference*B\n",
    "    # A -> B        \n",
    "    else:\n",
    "        variation = difference*A\n",
    "        \n",
    "    # control the pace of change with alpha\n",
    "    variation = alpha*variation  \n",
    "    \n",
    "    # update the population    \n",
    "    A = A + variation\n",
    "    B = B - variation   \n",
    "    \n",
    "    # save the values to a list for plotting    \n",
    "    believersA.append(A)\n",
    "    believersB.append(B)\n",
    "    \n",
    "    # advance time to next iteration\n",
    "    t = t + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you plot this you will get a different result than previously:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f6c52920b38>]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEECAYAAADOJIhPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcHHWd//FX1XTP3XNkzkyOScKRxFzcp5AsKJccuku+\ngsIqqD91FfH4uYKP3R947O8hP9dVVt11XVmEXQ++osh9xTByQ4BwJSEhhJxDZnLNfXd/f39UJz3E\nkEymZ6Z7pt7Px6Mf3VVd3f2Zbzr17vpW1bc85xwiIhJOfqYLEBGRzFEIiIiEmEJARCTEFAIiIiGm\nEBARCTGFgIhIiEUOtYAx5hbgQqDJWrswOa8cuAOoBzYCxlrbmnzueuBqYAC41lr7yOiULiIi6RrK\nlsCtwLn7zbsOWGatnQ0sB64HMMa8DzDAXOB84N+MMd5QCjHGLBlizROe2iJFbZGitkhRW6Sk2xaH\nDAFr7ZPAnv1mXwLclnx8G/Dh5OOLgd9aawestRuBN4GThljLkiEuFwZLMl1AFlmS6QKyyJJMF5BF\nlmS6gCyyJJ0XD3efQLW1tgnAWrsdqE7OnwJsGbTctuQ8ERHJQiO1Y1hjT4iIjEOH3DH8HpqMMTXW\n2iZjTC3QnJy/DZg2aLmpyXl/IdmPtWTvtLX2BuCGYdYzoVhrQW0BqC0GU1ukqC1SrLUYYwbParDW\nNgz19UMNAS952+se4JPATcAngLsHzf+VMeaHBN1ARwLPv0fhDcDgQm9obGwcYjkTWywWo729PdNl\nZAW1RcpItYWLx6GrAzrboasTujpxXR3Q3RVMd3dCTxd0d+G6u6CnO3iupxt6u6GnB/p6IBqF3HzI\nG3TLzUvecvH2Pc6DaDCPaPKWmwuRKF40GjwXiQTvFxl02zcdgZxI8DgnB8/z9L0YpK6uDmvtjcN9\nvXeoUUSNMb8m+MVeATQRpO8fgd8R/OrfRHCIaEty+euBTwH9HN4hok4hENAXPEVtkfJebeF6uqGt\nBdpboa0F194C7W3BdEcbrqMNOtqhoy1Y+ff2QGERFBYnb0V4BUXJeUVQUAT5hVBQiFdQCPkFwXR+\nQeqWl4fn52SgFQL6XqTU1dXBu3+kH5ZDhsAYUggk6QueEua2cP390LIL9uzE7dlFXlc7vU3vQOse\nXMtuaN0drPxdAkrKIVYKJWV4sVKIlUBxKcRK8YpLoDgGRbHgPr8Qzx/f54mG+Xuxv3RDYLj7BEQk\nTa6vF3Y1w47tuF3NsLMZt6sJdu+E3TuC7prSSVBeiVdegaueDBXVcMQc/NLy4LnSMsgrwPOGvQ6Q\nkFMIiIwiN9APO7bDO1txzY3Q1BjcN28PumcmVUFlDV5lTXBffwReRTVUVAW/6gd1uRTEYgzo16+M\nMIWAyAhw/f2wfStu2yZ4ZwuucTO8szX4pT+pCmqn4NXUQf2R+CedCdWTobwio/3qIqAQEDlsrr0V\ntmzAbd4AW97Gbd0Y/NqvrMGbUg+Tp+GfsgRqp0H15OAIGJEspRAQOQjX1Qkb1+E2rsdtfBM2rYfu\nbpg2E2/6LJh3LP65fw2Tp+JFczNdrshhUwiIJDnngi6d9WvgrTW4DeuCHbTTZ+HNPBrvxDPwLr0K\nqmq1I1YmDIWAhJZLxIPunLWv49a9DuvXBMfGHzk3OALnrItgSj1ejvrtZeJSCEhoOOegcTPujVdx\na16BN1dB6SS82fPxTl6Md8Xn8coqMl2myJhSCMiE5trbcKtXwqqVuNUvQzSKN3cR3kln4v3tF/BK\nyjNdokhGKQRkQnHOwdaNuFdX4F57AbZtgtkL8OYdi3+hwauuy3SJIllFISDjnovHYd3ruFeex738\nHHge3qKT8C+6HI6er0M0RQ5CISDjkhvohzdexb34dLDir6jGO+Zk/Gv+Eeqm6+gdkSFSCMi44RJx\nWPs6bsUTuJXPQHUd3vGn43/IBMMuiMhhUwhIVnPO4Ta/hXumAbfi8eBonpPOwP+HHwZj7IhIWhQC\nkpVc6x7csw20P9dAoqsT7+Ql+F/7J7zJUzNdmsiEohCQrOHicXhtBYknl8Gbq/COPYXCq6+lq27G\nuB//XiRbKQQk49yuZtyTj+KefDTYwXvGOXif/hpefgGRWAxPwyeLjBqFgGSESyRg9UoSDQ/C+jV4\nJy/G//K3glE4RWTMKARkTLnuLtxTj+IeewDy8vH+6kN4n/nfeHn5mS5NJJQUAjImXPM7uOX34Z55\nLDh796ovwxFzdDy/SIYpBGRUubfeIPHIXbDudbz3n4N/w814k6oyXZaIJCkEZMS5RAJee4HEQ7+H\nPbvwPvhhvKu+jJdfkOnSRGQ/CgEZMS4eD87mfej34Pt451+Kd9xpGo9fJIspBCRtbqAf9/TyYOVf\nXol/6VUw71j194uMAwoBGTbX3x8c3//QnVA7Df+qL+Md9b5MlyUih0EhIIfNDfTjnlyGe/B3MGUG\n/me/gTdrdqbLEpFhUAjIkLl4HPfMctx9d8DkqVr5i0wACgE5JJdIBDt87/kNlFfgf/prwcXYRWTc\nUwjIe3LOwesvkfjD7RCN4l/xeby5izJdloiMIIWAHJDbsJbE72+Dthb8j1wJx56io31EJiCFgLyL\na2okcdftsGEd3kWX4Z12to7zF5nAFAICgGtvw933W9zzf06e4fsVvLy8TJclIqNMIRByrr8ft/xe\n3EN/wDvxDPxv/xterDTTZYnIGFEIhJRzDl56msSdv4Qp9fjf+B5erS7dKBI2CoEQcpvWk7jjF9Dd\njf+Ja/DmLMx0SSKSIQqBEHGte3B/uB236iW8Sz6Od/rZeL52+oqEWVohYIy5HrgCiAOvAVcBRcAd\nQD2wETDW2tb0ypR0uP5+3J/uwT38B7zTP4D/nX/HKyjMdFkikgX84b7QGFMPfAY41lq7kCBQLgeu\nA5ZZa2cDy4HrR6JQOXzOOdwrK0jc+EXc+jX4130f/9KrFAAisk86WwJtQB9QZIxJAAXANoKV/uLk\nMrcBDQTBIGPIbd8W9Pvv3I5/+Wfx5h+X6ZJEJAsNe0vAWrsH+AGwmWDl32qtXQbUWGubkstsB6pH\nolAZGtfTReLOX5K46e/x5i7Ev+FfFQAi8p6GvSVgjJkFfIWg778V+J0x5uOA22/R/af3vn4JsGTv\ntLWWWCw23HImlNzc3MNuC+cc/U/9ie5f/weR+cdR8M+/xC+bNEoVjp3htMVEpbZIUVu8mzHmxkGT\nDdbahqG+1nPugOvooXyoAT5orf1McvpK4BTgLGCJtbbJGFMLPGatHcqQk66xsXFYtUw0sViM9vb2\nIS/vtr5N4tf/Ab09QdfPBBrh83DbYiJTW6SoLVLq6uoAhj2wVzr7BNYC/2iMyQd6gbOBFUAH8Eng\nJuATwN1pfIYchOvqwN39a9yKJ/Au/hjemefokE8ROSzp7BN4BbgdeBF4hSCJfk6w8v+gMWYtQTB8\nbwTqlEFcIkHiqWUk/s8XYGAA/9s/xV9yvgJARA7bsLuDRoG6g5IOtqnrNq0Pun6cw//YZ/FmHDW2\nxY0xbfanqC1S1BYpmewOkjHkOtpwf/wf3Mpn8T5yZTDEsz/sDTkREUAhkPVcIo574lHc3b/CO+H9\nwSifRcWZLktEJgiFQBZzb70RdP3k5eF/5dt402ZmuiQRmWAUAlko0bKbxG0/wa15Fe9vPoF38mJd\n2lFERoVCIIu4gX7cn+6j/eHfw2kfwP/OT/HyNc6PiIwehUCWcK+9GIz1U1VL7Fs/pitWnumSRCQE\nFAIZ5rZvI2FvgaZt+B/9NN7CE8mJxUCHv4nIGFAIZIjr6sDddwfumeV45/0N3uevx4tGM12WiISM\nQmCMuXgc98TDuHt/i7foJPxv/QSvRF0/IpIZCoEx5F57kcTv/gtKyvCvvRFv+qxMlyQiIacQGANu\ny9sk7rwVdjbjL70KFp2kQz5FJCsoBEaR270jGOph1Uq8Cz+Kd8a5eBE1uYhkD62RRoHraMM98Dvc\n08vxFp+P/92f6bq+IpKVFAIjyPV04Zbdi/vTPcE4Pzf+GG8CXN1LRCYuhcAIcL29uIb7cQ/fhTf3\nGPzrv49XXZfpskREDkkhkAbX24t74iHcQ3fBEXPwv/ZPeFOmZ7osEZEhUwgMg+vpxv35Qdwjf4Qj\n5+J/6R/xph+R6bJERA6bQuAwuPY23PL7cH9+EO/o+cHwzlNnZLosEZFhUwgMgWt+B/ene3HPNuAd\nfxr+338Pr3ZKpssSEUmbQuA9OOfgzVUkHr0H1q/Ce/85+N/6MV5ZRaZLExEZMQqB/bieLtyzf8b9\n+UHo78f7wEV4n/4qXl5+pksTERlxCgGSv/o3vol7ahluxZMwZwH+0qthzkJdzF1EJrRQh4Dbswv3\n/OO4p5bBQD/eaWcHJ3iVq8tHRMIhdCHg2lpwLz2DW/E4bN2Ed+zJ+Fd8Ho6ap0HdRCR0JnwIOOeg\naRvuledxLz8H2zbjzT8O/4OXwLzjdSEXEQm1CRkCrnUPbt0qWPMybtVKSCTwFp6If4EJ+vm14hcR\nASZACLjeXtj6Nm7LBtiwDrd+NXS2wxFz8d53TPCLv3aqunpERA4g60PAOQddHbBnF7Tsxu3eAc2N\nuKZG2L4NdjcHK/lps2DWbPxzPwKTp+moHhGRIciqEIh/96uQiEMiAT3d0N0FPV2QVwBlk6BsEl55\nJdTU4Z/6V1BdB7VT8CLq3hERGY6sCgH/458H3w9ueflQUAQFhboal4jIKMmqtas386hMlyAiEirq\nOBcRCTGFgIhIiCkERERCTCEgIhJiae0YNsaUAr8A5gMJ4GpgHXAHUA9sBIy1tjW9MkVEZDSkuyVw\nM/CAtXYusAh4A7gOWGatnQ0sB65P8zNERGSUDDsEjDElwBnW2lsBrLUDyV/8lwC3JRe7Dfhw2lWK\niMioSKc7aCaw0xhzK8FWwAvAl4Eaa20TgLV2uzGmOv0yRURkNKTTHRQBjgN+aq09Dugk6Apy+y23\n/7SIiGSJdLYEtgJbrLUvJKd/TxACTcaYGmttkzGmFmg+0IuNMUuAJXunrbXEYrE0ypk4cnNz1RZJ\naosUtUWK2uLdjDE3DppssNY2DPW1nnPD/6FujPkz8Blr7TpjzA1AYfKp3dbam4wx3wDKrbXXDeHt\nXGNj47BrmUhisRjt7e2ZLiMrqC1S1BYpaouUuro6gGGPlZ/u2EFfAn5ljIkCG4CrgBzAGmOuBjYB\nJs3PEBGRUZJWCFhrXwFOPMBTH0jnfUVEZGzojGERkRBTCIiIhJhCQEQkxBQCIiIhphAQEQkxhYCI\nSIgpBEREQkwhICISYgoBEZEQUwiIiISYQkBEJMQUAiIiIaYQEBEJMYWAiEiIKQREREJMISAiEmIK\nARGREFMIiIiEmEJARCTEFAIiIiGmEBARCTGFgIhIiCkERERCTCEgIhJiCgERkRBTCIiIhJhCQEQk\nxBQCIiIhphAQEQkxhYCISIgpBEREQkwhICISYgoBEZEQUwiIiISYQkBEJMQUAiIiIRZJ9w2MMT7w\nArDVWnuxMaYcuAOoBzYCxlrbmu7niIjIyBuJLYFrgdWDpq8DlllrZwPLgetH4DNERGQUpBUCxpip\nwAXALwbNvgS4Lfn4NuDD6XyGiIiMnnS3BH4IfB1wg+bVWGubAKy124HqND9DRERGybBDwBjzIaDJ\nWvsy4B1kUXeQ50REJIM854a3jjbG/F/gCmAAKABiwF3ACcASa22TMaYWeMxaO/cAr18CLNk7ba29\nob29fVi1TDS5ubn09fVluoysoLZIUVukqC1SYrEYxphvDZrVYK1tGOrrhx0CgxljFgNfSx4d9P+A\nXdbam4wx3wDKrbXXDeFtXGNjY9q1TASxWAwFYkBtkaK2SFFbpNTV1cHBe2MOajTOE/ge8EFjzFrg\n7OS0iIhkoRHZEhgh2hJI0q+cFLVFitoiRW2Rko1bAiIiMk4oBEREQkwhICISYgoBEZEQUwiIiISY\nQkBEJMTSHkp6JH35gbfxPfA9j4KIT2GuT2E0h5K8HCoKI1QURKgsijI5lktJXk6myxURGfeyKgSu\nOWUyCeeIJ6BnIEFnf5yuvgStvXG2t/exurmLHZ0DNLb3keNBXUkuM8rymVmex6xJwX1ujjZuRESG\nKqtC4IhJ+UNazjlHa2+cba19bGzpZf3uHh5Z38K2tj5mlOfzvqoC5lYXsKCmkMKothhERN5LVoXA\nUHmeR1l+hLL8CPNqCvfN7xlIsG5nN2/s6Ob+tXv4l6feYVZ5HosmF3HilGJmlefhecM+sU5EZMIZ\nlyHwXvIjPgtri1hYW4QBegcSrN7RzcrGDv75yW30xR0nTS3m1Gkx5lUXkuMrEEQk3CZUCOwvL+Jz\n7OQijp1cxFXHVbO1rY/ntnTwy5XN7O4a4PT6Es6cUcLRFfnaQhCRUJrQITCY53lMK81jWmkel86v\nYGtbL09uaudHT7+D58HZs0pZMrOEisJopksVERkzoQmB/U0tyeOyBXl8dH4Fb+zs5k9vtXLN/W8z\nr7qQ848q45jJRfjaOhCRCU5DSQ/SM5Dg8Y1tPLhuD139Cc47qowPHllGce7YHmGkYXJT1BYpaosU\ntUVKukNJh3ZL4EDyIz7nHFnGB48oZd2uHu5fu4fP3v0Wi2eUcOHsSdSV5Ga6RBGREaUQOADP85hd\nWcDsygJ2dfXzwLoWvvHIJuZVF/DX76vg6MqCTJcoIjIiFAKHUFEY5cpjqrh0XgXL3mrhpie2MTmW\ny6XzKlhUW6ijikRkXFMIDFFB1OeiOZM4/+hyHt/Yxs9faKI4N4fLFlRw7OQihYGIjEsKgcMU8T3O\nmlXK4hklPLW5nVtfaubXEZ+PLaxUGIjIuKMQGKYc3+PMGSW8vz7G05vbueXFZkrycvjYokoW1BRl\nujwRkSFRCKTJ9zzeX1/CqdNiPL6xjZ88u53aWC5/e0zVkAfEExHJFIXACMnxPf5qVinvry/h0bda\n+M5jW5hXU8gVi6qYHNOhpSKSnTT4/giL5nhccHQ5P7vkCOrL8vj6w5v4+QtNtPUMZLo0EZG/oBAY\nJfkRHzO/kp9cOBPnHH9339vcuWoXffFEpksTEdlHITDKyvIjfPbEWm46p543d3XzhXs38PjGNrJo\nuA4RCTGFwBiZUpLL9WdO5UunTuau1bv4xiObWLuzO9NliUjIKQTG2IKaIn5w/gzOPbKM7z2+jR8+\n3ciurv5MlyUiIaUQyADf8zj7iDJ+etFMKgujXPvARuzrO7W/QETGnEIggwqjOVx5TBU/OK+eDbt7\nuOa+t3l2S7v2F4jImNF5AlmgpjiX686cysvvdPKLF5t4ZEM7Vx1TwbTSvEyXJiITnLYEssgxk4v4\n0QUzOWV6Kd98dDO3vtRMV38802WJyASmEMgyEd/jbxbW8uMLZ9LRF+fv7n2bxza0qotIREaFuoOy\nVFl+hGtOmczand38fEUTD69v4X+dUMMsjUckIiNIWwJZbnZlAd8/r56zZpVy42Nb+Nnz22nvVReR\niIwMhcA44Hse5xxZxk8vnAXAF+/bwCPrW0ioi0hE0jTs7iBjzFTgdqAGSAD/aa39V2NMOXAHUA9s\nBIy1tnUEag29WF4OnzuplnOOLONnK5p4ZH0Lnz2xhqMqdM1jERmedLYEBoCvWmvnAacCXzDGzAGu\nA5ZZa2cDy4Hr0y9TBps1KZ/vnTOdC44u558atvKTZ9/RKKUiMizDDgFr7XZr7cvJxx3AGmAqcAlw\nW3Kx24APp1uk/CXfCy5z+dOLZpEf8fnifW9z/9o9xBPqIhKRoRuRfQLGmBnAMcCzQI21tgmCoACq\nR+Iz5MCKcnP49Ak1fOcD03l6cxtfe2gjq5u7Ml2WiIwTaR8iaowpBu4ErrXWdhhj9v8pesCfpsaY\nJcCSvdPWWmKxWLrlTAi5ubmH3RbzYzFunlpBw1u7+Zent7CoLsZnTp5GVfH4vqrZcNpiolJbpKgt\n3s0Yc+OgyQZrbcNQX+ulcxKSMSYC3Ac8aK29OTlvDbDEWttkjKkFHrPWzh3C27nGxsZh1zKRxGIx\n2tvbh/367v4Ed67axcPrW/jI3ElcPKecaM74PBAs3baYSNQWKWqLlLq6OgBvuK9Pd83wX8DqvQGQ\ndA/wyeTjTwB3p/kZcpgKoj5XHlPF98+tZ82Obq65/22e36qB6UTkLw17S8AYczrwOPAaQZePA74J\nPA9YYBqwieAQ0ZYhvKW2BJJG+lfOS40d3PJiM1VFUT51fPW4GphOv/hS1BYpaouUdLcE0uoOGmEK\ngaTR+IIPJBwPrtuDfX0XZ84o4bIFlcTyckb0M0aD/rOnqC1S1BYpme4OknEi4ntcNGcSP7lwJvGE\n4wv3buC+tbsZ0CGlIqGmEAiZ0vwInzuplu98YDrPb+3g2vvf5oVtHdpfIBJSGkU0pOrL8vjWWdN4\nsbGTW19q5u43Ilx9XDUzyzVKqUiYKARCzPM8TphSzDGTi3h0fQs3Lt/CcXVFfGxhFVVF0UyXJyJj\nQN1BQsT3OP/ocv794llUFET5ygNvc/vKZjr6NGS1yESnEJB9CqM5XHFMFT/60Exae+P83T0buGv1\nLnoHEpkuTURGibqD5C9UFka55pTJbG7t5X9e3sG9a/dw+YJKzppVSo4/7CPRRCQLKQTkPU0vzeOb\ni6fyxo5u/vuVHfxh9S4uX1jF++tj+J7CQGQiUHeQHNKcqgK+e/Y0PndSLfe+sZsvP7CRZza368pm\nIhOAtgRkSDzPY1FtEQvPLeSFbZ385rUd3PH6Ti5bUMnJU4vxtGUgMi4pBOSweJ7HiVOLOWFKEc9v\n7eA3r+3kt6/tZOn8Ck6dpm4ikfFGISDD4nkeJ0+LcdLUYlZs68C+vovfvLqTS+dVcEZ9iXYgi4wT\nCgFJi+d5nDQ1xolTiln5Tid3rtrFr1/dyUfmTuLsI0rJHafXMRAJC4WAjAjP8ziurpjj6opZs6OL\n36/azR2v7eSC2eWcd1Q5JeNgxFKRMFIIyIibW1XIPywpZHNLL39cs5vP3/MWi2eUcNGcSUyOje/L\nXYpMNAoBGTXTy/L40qmT2dVVyf1r9/D3D29idmUBF88pZ0FNoY4oEskCCgEZdRWFUf722Go+uqCS\nx95u5T9WNOF7cMHR5SyeWUJhVF1FIpmiEJAxkxfxOe+ocs49soxXm7p4YN0e/ueVHSyeUcK5R5VT\nXzZ+LnspMlEoBGTM7T3xbFFtETs6+3kkOYx1VVGU844q4/TpMfIiOqpIZCwoBCSjqoqifHxRFZct\nqOSFbR08vL6FW15s4rTpMT5wRBnHFxdnukSRCU0hIFkhxw9OPjt5WoxdXf08tqGNHz3dSOS5Js6s\nL2bxjFKqi3WhG5GR5mXRtWVdY2NjpmvICrFYjPb29kyXkXHOOTZ3eTy4ajtPbm5nemkuZ84o4bRp\nMUryw/f7Rd+LFLVFSl1dHcCwD7UL3/8kGTc8z2N+bYz6IvjU8TW81NjBE5vauH3lDmZXFnB6fYyT\np8aI6UQ0kWFTCMi4EM1JdRf1DCR4fmsHT29u45YXmzm6soDTpsU4cWoxkwr0lRY5HPofI+NOfsTn\nzBklnDmjhJ6BBC82dvDM5nZue7mZqSW5+8Yyml6aqxPSRA5BISDjWn7E5/TpJZw+vYT+uGNVcxfP\nbW3nuw1bADi+rpjj64qZX1NIQVSHnYrsTzuGs5B2eqUMty2cc2xp7eOFbR289E4nb+7q4ciKfI6d\nXMSi2kJmleePu+Gu9b1IUVukaMewyAF4nsf0sjyml+Xx1/Mq6O5P8HpTFyu3d3LzM++wu3uA+dWF\nLKgpZF51IfVleeMuFERGgkJAQqEg6nPi1GJOnBqcfLane4BXt3eyqrmbB99soaVngLmVBcytKmRu\nVQFHVuTrrGUJBYWAhFJ5QYTFM0tZPLMUCEJh9Y4u1uzo5taVzWxu6WVaaR5HV+YzuzIIhbpYri6f\nKROOQkCEIBT27mAG6B1IsGF3D2t3dfPc1g5+9cpOOvriHDEpnyMm5TOzPI9Zk/KZEstVN5KMawoB\nkQPIi/jMrS5kbnXhvnltPQOs393DW7t7eH5rB799bSe7uwaYVppHfVkeM8rzmF4a7Icoz8/R4aky\nLigERIaoJD+y7xKae3X1x9nc0sfGlh427unluS3tbG7tI+EcU0vymFqay5SS5C2WS01xLtEchYNk\nD4WASBoKoznMqSpgTlXBu+a39AywpbWXbW19bGvr4/WmLhrb+9jZOUBFYYTaWC61xdHkLZfq4ig1\nRVGKcn1tQciYUgiIjIKy/Ahl+REW1BS9a/5AwtHU0c/29j62d/SzvaOPNTu6ae7sp6mjH4CqwiiV\nRRGqiqJUFkaoKAzup8ej5CUSOulNRtSohYAx5jzgR4AP3GKtvWm0PktkvIj43r7uof055+joS7Cz\nq58dnf3s6BxgZ1c/W97pZFdXP7t7mtnZ2UeO51FeEGFSQQ5lBRHKCyLJ0MmhLD9CSV4Opfk5lORF\nyI942rKQgxqVM4aNMT6wDjgbaARWAJdZa984yMt0xnCSzoZMUVukxGIx2tra6OxPsLt7gJbuAfZ0\nD7CnZ4DWnjgtPXFaewZo6YnT3hvMc0AsN4dYXg4lecF9cXK6KNenODe4L4qm7gtzfYqiPtGc7N3i\n0PciJVvPGD4JeNNauwnAGPNb4BLgYCEgIofgeR7FucGKfHrpoa/J3DuQoK03TntvnLbeOB19weOO\nvjitPXG2tfXR2Renoy9BV3+czr4Enf0JuvrieJ5HYdSnMOpTEPUpiAT3+ZHUdF7EJz/ikb/vsU9e\njkfeoPvciEdejk80xyM3J3isw2qzx2iFwBRgy6DprQTBICJjKC/iUxXxqSo6vKuyOefoTzi6+hN0\n9SXoHkjQ3Z+8DSToHUjQM2jenp44Pf0J+uIJeuOOnoEEfQOO3niCvrijdyC470tOA+TmeERzfHJ9\nj2iOR8QPQiKSnI4Omr93Xo4XPC7M30NioH/fczl+cHW6iJ9aJseHHM8jx/fI8cBP3u+d5ycf+z74\n3qDp5LKDpz0GLQfBPM/D84LH/uDlPPC84Kf5eOiK045hEfkLnheskHNzfMryR/794wm3LxT6E47+\neHDrizuxXAUdAAAEsklEQVT6EwkGkvMGEm6/x8HO9ZxoLp3dCeKJIKz6+h1xF7zvQMIRd8Gy8YQj\n4VLPJZwjnoAB50g4SCSC5xLOJW97l2Pf9LsfO5yDBIMeD5rvSL0G9gYB+wIiFQ6DpgcFRuq5QdPJ\nN/II3of93ueez9Wl9W8xWiGwDZg+aHpqct4+xpglwJK909bavX1bQtDnKQG1RYraQg7EGHPjoMkG\na23DkF/snBvx29KlS3OWLl26funSpfVLly7NXbp06ctLly6de4jX3DgatYzHm9pCbaG2UFuMVVuM\nyu5/a20c+CLwCLAK+K21ds1ofJaIiAzfqO0TsNY+BMwerfcXEZH0ZdOBwA2ZLiCLNGS6gCzSkOkC\nskhDpgvIIg2ZLiCLNKTz4my6vKSIiIyxbNoSEBGRMaYQEBEJsaw4WSzMg80ZY6YCtwM1BOeg/Ke1\n9l+NMeXAHUA9sBEw1trWjBU6RpLjTr0AbLXWXhzWdgAwxpQCvwDmE3w3riYYkytU7WGMuR64AogD\nrwFXAUWEpB2MMbcAFwJN1tqFyXnv+f8i2V5XAwPAtdbaRw72/hnfEkj+p/8JcC4wD7jcGDMns1WN\nqQHgq9baecCpwBeSf/91wDJr7WxgOXB9BmscS9cCqwdNh7UdAG4GHrDWzgUWEYy9Far2MMbUA58B\njk2uACPA5YSrHW4lWD8OdsC/3xjzPsAAc4HzgX8zxhx07IqMhwCDBpuz1vYDewebCwVr7XZr7cvJ\nxx3AGoIzrC8Bbksudhvw4cxUOHaSW0UXEPz63St07QBgjCkBzrDW3gpgrR1I/tILW3u0AX1AkTEm\nAhQQjD4Qmnaw1j4J7Nlv9nv9/RcTnJc1YK3dCLzJIcZty4YQONBgc1MyVEtGGWNmAMcAzwI11tom\nCIICqM5gaWPlh8DXgcGHrIWxHQBmAjuNMbcaY14yxvzcGFNIyNrDWrsH+AGwmWDl32qtXUbI2uEA\nqt/j799/fbqNQ6xPsyEEBDDGFAN3EvThdfDuFSEHmJ5QjDEfIujzfJmDj40+odthkAhwHPBTa+1x\nQCdBF0DYvhezgK8Q9H3XEWwRfJyQtcMQDPvvz4YQOORgcxNdcjP3TuC/rbV3J2c3GWNqks/XAs2Z\nqm+MnA5cbIzZAPwGOMsY89/A9pC1w15bgS3W2heS078nCIWwfS9OAJ6y1u5ODkdzF3Aa4WuH/b3X\n378NmDZouUOuT7MhBFYARxpj6o0xucBlwD0Zrmms/Rew2lp786B59wCfTD7+BHD3/i+aSKy137TW\nTrfWziL4Diy31l4J3EuI2mGv5Kb+FmPM0clZZxOMwxWq7wWwFjjFGJOf3MF5NsGBA2Frh32jSie9\n199/D3CZMSbXGDMTOBJ4/qBvnA1nDCcPEb2Z1CGi38twSWPGGHM68DjBoW8uefsmwT+cJUj1TQSH\ngLVkqs6xZIxZDHwteYjoJMLbDosIdpJHgQ0Eh0bmELL2MMZ8nWCFFwdWAp8GYoSkHYwxvyYYdr8C\naAJuAP4I/I4D/P3JQ0Q/BfQzhENEsyIEREQkM7KhO0hERDJEISAiEmIKARGREFMIiIiEmEJARCTE\nFAIiIiGmEBARCTGFgIhIiP1/8kgiLUt+9b8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6c52920b00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(believersA)\n",
    "plt.plot(believersB) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Have a go at changing the values of $K_a$ and $K_b$ and see what happens. Can you see any equilibrium where both traits coexist?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a number of functions we can use to dynamically change the 'attractiveness' of each trait. Try the following ones:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np # stick this line at the beginning of the script alongside other 'imports'\n",
    "\n",
    "\n",
    "def attractiveness2(Ta, Tb):\n",
    "    # temporal autocorrelation with stochasticity (normal distribution)\n",
    "    # we get 2 samples from a normal distribution N(0,1)\n",
    "    Ka, Kb = np.random.normal(0, 1, 2)\n",
    "    # compute the difference between Ks\n",
    "    diff = Ka-Kb\n",
    "    # apply difference of Ks to attractiveness\n",
    "    Ta += diff\n",
    "    Tb -= diff\n",
    "    return Ta, Tb\n",
    "\n",
    "def attractiveness3(Ta, Tb):\n",
    "    #  anti-conformism dynamics (more population means less attractiveness)\n",
    "    \n",
    "    # both values initialized at 0\n",
    "    Ka = 0\n",
    "    Kb = 0\n",
    "    \n",
    "    # first we sample gamma with mean=last popSize of A times relevance    \n",
    "    diffPop = np.random.gamma(believersA[t])\n",
    "    # we sustract from this value the same computation for population B\n",
    "    diffPop = diffPop - np.random.gamma(believersB[t])\n",
    "    \n",
    "    # if B is larger then we need to increase the attractiveness of A\n",
    "    if diffPop < 0:\n",
    "        Ka = -diffPop\n",
    "    # else A is larger and we need to increase the attractiveness of B\n",
    "    else:\n",
    "        Kb = diffPop\n",
    "    \n",
    "    # change current values\n",
    "    Ta = Ta + Ka\n",
    "    Tb = Tb + Kb\n",
    "    \n",
    "    return Ta, Tb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the functions above (just change 'attractiveness' to, e.g., 'attractiveness2' in the main code) to see what happens when we add small dynamic variations to the process. What happens when the initial conditions are changed? What if we look at a much longer time scale? "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}