{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from collections import namedtuple, deque\n",
    "\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Route:\n",
    "Numbers from Google maps at 5pm on a Wednesday\n",
    "\n",
    "* Bus station\n",
    "    * 5 minutes / 0.9 miles / 4 junctions\n",
    "* Oldfield Park\n",
    "    * 2 minutes / 0.4 miles / 0 junctions\n",
    "* Broughton Hayes\n",
    "    * 3 minutes / 0.6 miles / 2 junctions\n",
    "* Corn Street\n",
    "    * 5 minutes / 0.7 miles / 3 junctions\n",
    "* Pulteney Court\n",
    "    * 3 minutes / 0.5 miles / 1 junction\n",
    "* St Marys Church\n",
    "    * 5 minutes / 1.7 miles / 6 junctions\n",
    "* University of Bath\n",
    "\n",
    "\n",
    "\n",
    "## Timetable\n",
    "Times acording to timetable\n",
    "* Oldfield Park -> Corn Street\n",
    "    * 6, 9, 11 minutes at deferent times\n",
    "* Corn street -> Uni\n",
    "    * 12, 16, 11 mintues at different times\n",
    "* Every 6 minutes during the day\n",
    "* U2 is every 15 to factor in percentage it carrys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.25       0.54545455]\n",
      "[0.48780488 0.51282051 0.47244094 0.51282051 0.51282051 0.44444444]\n"
     ]
    }
   ],
   "source": [
    "to_minutes = numpy.array([5, 2, 3, 5, 3, 5], dtype=float)\n",
    "to_miles = numpy.array([0.9, 0.4, 0.6, 0.7, 0.5, 1.7])\n",
    "to_junctions = numpy.array([4, 0, 2, 3, 1, 6], dtype=float)\n",
    "\n",
    "base = 1.0\n",
    "a = numpy.concatenate((to_miles[:,None], to_junctions[:,None]), axis=1)\n",
    "answer = numpy.linalg.lstsq(a, to_minutes - base, rcond=None)\n",
    "print(answer[0])\n",
    "\n",
    "to_correction = to_minutes / (base + 1.45 * to_minutes + 0.5 * to_junctions)\n",
    "print(to_correction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.08490807 1.6641328  3.08908929 5.75467603 4.10598185 4.81128028]\n",
      "22.510068320518467\n"
     ]
    }
   ],
   "source": [
    "def draw_travel_time():\n",
    "    \"\"\"Returns a draw of the travel time, in (float) minutes to get from the previous\n",
    "    stop to the given stop, as an array indexed by each stop - effectively drive time\n",
    "    for a bus along the route. For stop 0 that is time to drive there from the bus station;\n",
    "    note that this is calculated as for rush hour!\"\"\"\n",
    "    base = 1.0\n",
    "    drive = (1.45 + 0.2 * numpy.random.standard_normal(to_minutes.shape[0])) * to_minutes\n",
    "    wait = numpy.random.exponential(0.5 * to_junctions)\n",
    "    return  (base + drive + wait) * to_correction\n",
    "\n",
    "\n",
    "\n",
    "# Quick test of the above...\n",
    "tt = draw_travel_time()\n",
    "print(tt)\n",
    "print(tt.sum())\n",
    "\n",
    "\n",
    "\n",
    "# Histogram...\n",
    "drive_time = numpy.empty(1024)\n",
    "for i in range(drive_time.shape[0]):\n",
    "    drive_time[i] = draw_travel_time().sum()\n",
    "\n",
    "#plt.figure(figsize=(12,12))\n",
    "#plt.hist(drive_time)\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Student populations\n",
    "Using 2011 census, as updated by ONS to 2017 based on some unknown model, giving population by ward for 18-22:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bus stop usage weightings:\n",
      "  Oldfield Park:\t 1516\n",
      "  Broughton Hayes:\t 2693\n",
      "  Corn Street:\t 2163\n",
      "  Pulteney Court:\t 732.0\n",
      "  St Marys Church:\t 2120.5\n",
      "  University of Bath:\t 0.0\n",
      "Total = 9224.5\n",
      "\n",
      "Expected staff percentage each day: 0.928\n",
      "Expected student percentage each day: 0.9535353535353536\n",
      "Staff bus usage per day = 332.192448\n",
      "Bus usage per day = 9789.845248000001\n"
     ]
    }
   ],
   "source": [
    "ward_pop = {'Abbey' : 787,\n",
    "            'Bathavon North' : 274,\n",
    "            'Bathavon South' : 215,\n",
    "            'Bathavon West' : 801,\n",
    "            'Bathwick' : 3530,\n",
    "            'Chew Valley North' : 105,\n",
    "            'Chew Valley South' : 134,\n",
    "            'Clutton' : 117,\n",
    "            'Combe Down' : 403,\n",
    "            'Farmborough' : 107,\n",
    "            'High Littleton' : 117,\n",
    "            'Keynsham East' : 306,\n",
    "            'Keynsham North' : 284,\n",
    "            'Keynsham South' : 266,\n",
    "            'Kingsmead' : 1376,\n",
    "            'Lambridge' : 221,\n",
    "            'Lansdown' : 490,\n",
    "            'Lyncombe' : 373,\n",
    "            'Mendip' : 153,\n",
    "            'Midsomer Norton North': 299,\n",
    "            'Midsomer Norton Redfield': 280,\n",
    "            'Newbridge': 495,\n",
    "            'Odd Down': 471,\n",
    "            'Oldfield': 1516,\n",
    "            'Paulton': 268,\n",
    "            'Peasedown': 407,\n",
    "            'Publow and Whitchurch': 147,\n",
    "            'Radstock': 315,\n",
    "            'Saltford': 207,\n",
    "            'Southdown': 477,\n",
    "            'Timsbury': 92,\n",
    "            'Twerton': 620,\n",
    "            'Walcot': 506,\n",
    "            'Westfield': 394,\n",
    "            'Westmoreland': 2693,\n",
    "            'Weston': 330,\n",
    "            'Widcombe': 1464}\n",
    "\n",
    "\n",
    "\n",
    "# Bus stop usage estiamted by assigning each central ward to a stop...\n",
    "bs_usage = {}\n",
    "bs_usage['Oldfield Park'] = ward_pop['Oldfield']\n",
    "bs_usage['Broughton Hayes'] = ward_pop['Westmoreland']\n",
    "bs_usage['Corn Street'] = ward_pop['Abbey'] + ward_pop['Kingsmead']\n",
    "bs_usage['Pulteney Court'] = 0.5 * ward_pop['Widcombe']\n",
    "bs_usage['St Marys Church'] = 0.5 * ward_pop['Widcombe'] + 0.25 * ward_pop['Bathwick'] + ward_pop['Walcot']\n",
    "bs_usage['University of Bath'] = 0.0 # 0.75 * ward_pop['Bathwick']\n",
    "\n",
    "print('Bus stop usage weightings:')\n",
    "bs_usage_total = 0\n",
    "for key, value in bs_usage.items():\n",
    "    print('  {}:\\t {}'.format(key, value))\n",
    "    bs_usage_total += value\n",
    "print('Total = {}'.format(bs_usage_total))\n",
    "print()\n",
    "\n",
    "\n",
    "\n",
    "# Total bus usage...\n",
    "staff_per_day = 0.4 * 0.03 +  0.6 * 0.09 + 0.8 * 0.09 + 1.0 * (0.72 + 0.04 + 0.03)\n",
    "staff_renorm = 0.03 + 0.09 + 0.09 + 0.72 + 0.04 + 0.03\n",
    "staff_per_day /= staff_renorm\n",
    "print('Expected staff percentage each day: {}'.format(staff_per_day))\n",
    "\n",
    "students_per_day = 0.6 * 0.05 + 0.8 * 0.13 + 1.0 * (0.55 + 0.19 + 0.07)\n",
    "students_renorm = 0.05 + 0.13 + 0.55 + 0.19 + 0.07\n",
    "students_per_day /= students_renorm\n",
    "print('Expected student percentage each day: {}'.format(students_per_day))\n",
    "\n",
    "\n",
    "\n",
    "# Total bus usage...\n",
    "staff_bus_usage = 2841 * 0.126 * staff_per_day\n",
    "student_bus_usage = 15819 * 0.627 * students_per_day\n",
    "bus_usage = staff_bus_usage + student_bus_usage\n",
    "print('Staff bus usage per day = {}'.format(staff_bus_usage))\n",
    "print('Bus usage per day = {}'.format(bus_usage))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Usage\n",
    "Useful URLs:\n",
    "* http://www.bath.ac.uk/estates/docs/University_2014_15_Travel_Survey_Report.pdf\n",
    "* http://www.bath.ac.uk/student-records/statistics/stats-20141201-all.pdf\n",
    "* http://www.bath.ac.uk/student-records/statistics/stats-20161201-all.pdf\n",
    "* http://www.bath.ac.uk/student-records/statistics/stats-20171201-all.pdf\n",
    "\n",
    "* Full time students:\n",
    "    * 2014: 13136\n",
    "    * 2017: 15819\n",
    "* Full time staff equivalent:\n",
    "    * 2017: 2841\n",
    "\n",
    "## 2014 Travel Survey\n",
    "* 12.6% staff travel by bus, based on survey\n",
    "* 62.7% students travel by bus, based on survey\n",
    "* 55.3% of 18,695 arrived by bus, counting as people arrive\n",
    "\n",
    "* Days of week travelled on (reading from graph, so numbers not accurate):\n",
    "    * Staff: 2: 3%, 3: 9%, 4: 9%, 5: 72%, 6: 4%, 7: 3%\n",
    "    * Student: 3: 5%, 4: 13%, 5: 55%, 6: 19%, 7: 7%\n",
    "\n",
    "# Bus sizes\n",
    "* 85 (double decker), 70 (single decker) and 120 (bendy)\n",
    "* Numbers from TFL model, i.e. not correct numbers but what their data suggests the true carrying capacity is.\n",
    "\n",
    "* At new route intro: U1 runs on a core fleet of eight Euro VI buses and up to 25 U1 buses may be in operation at peak times. But the additional 17 are not of the environmentally-friendly Euro VI spec, which runs on a diesel engine but emits ultra low levels of harmful emissions.\n",
    "* After chaos: 13 extra buses = 1000 extra seats per hour\n",
    "\n",
    "* Euro 6/VI is the double decker\n",
    "\n",
    "# How bad it was\n",
    "All first day of term:\n",
    "* “This morning I set off at 8.45 for my 10.15 lecture. Of the 5 U2 buses that were supposed to arrive in the time that I stood at the bus stop (8.53, 9.09, 9.28, 9.48, 10.08), only two turned up, and both drove by without stopping. I ended up missing the first lecture of the year.” – Alice\n",
    "* Corn street, 9:24: 90 waiting\n",
    "* Waited at the Fox Hill U2 bus stop this morning from 8:30 until 10:10, one bus was a no show, the other four drove straight by. Two of the “full” buses did seem packed, so I don’t mind that they drove by, but the other 2 were definitely not full, and I feel that we could have been let on without any problem.” – Oliver\n",
    "* Waited 1hr 45mins to get on a bus in Southdown. Turned up at 8.15am and didn’t get on a bus until 10am, four U2 services in a row were full and one didn’t even show up” – Sarah"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notes\n",
    "* Mention waiting time paradox\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Journeys generated = 9784\n"
     ]
    }
   ],
   "source": [
    "bus_stops = ['Oldfield Park',\n",
    "             'Broughton Hayes',\n",
    "             'Corn Street',\n",
    "             'Pulteney Court',\n",
    "             'St Marys Church',\n",
    "             'University of Bath']\n",
    "bus_stop_ratios = numpy.array([bs_usage[stop] / bs_usage_total for stop in bus_stops])\n",
    "\n",
    "staff_travel_chance = numpy.array([0.4, 0.6, 0.8, 1.0])\n",
    "staff_travel_dist = numpy.array([0.03, 0.09, 0.09, 0.72 + 0.04 + 0.03])\n",
    "staff_travel_dist /= staff_travel_dist.sum()\n",
    "\n",
    "students_travel_chance = numpy.array([0.6, 0.8, 1.0])\n",
    "students_travel_dist = numpy.array([0.05, 0.13, 0.55 + 0.19 + 0.07])\n",
    "students_travel_dist /= students_travel_dist.sum()\n",
    "\n",
    "\n",
    "\n",
    "def draw_journeys(scale = 1.0):\n",
    "    \"\"\"Generates a random set of journeys for an entire day at the university. Loops the entire\n",
    "    body of staff and students who take buses, decides if each one will be travling that day. Return\n",
    "    is a numpy array, indexed [journey, 0] for when they arrive at the bus stop (minutes since midnight),\n",
    "    [journey, 1] for the time they would be late (minutes since midnight), [journey, 2] for which\n",
    "    bus stop. Array is all integers. scale is the size of the university as a function of current population.\"\"\"\n",
    "    # Staff...\n",
    "    staff = numpy.zeros((int(scale * 2841 * 0.126), 3), dtype=int)\n",
    "    staff[:,1] = 15 * numpy.random.randint(33, 44, size=staff.shape[0]) # 8:15 to 11:15 in 15 minute incriments\n",
    "    staff[:,0] = staff[:,1] - numpy.random.randint(30, 60, size=staff.shape[0]) #  Leave 30-60 minutes before\n",
    "    staff[:,2] = numpy.random.choice(len(bus_stops), size=staff.shape[0], p=bus_stop_ratios)\n",
    "    \n",
    "    staff_travel = numpy.random.choice(staff_travel_chance, size=staff.shape[0], p=staff_travel_dist)\n",
    "    staff_keep = staff_travel >= numpy.random.rand(staff.shape[0])\n",
    "    staff = staff[staff_keep,:]\n",
    "    \n",
    "    # Students...\n",
    "    students = numpy.zeros((int(scale * 15819 * 0.627), 3), dtype=int)\n",
    "    students[:,1] = (15 + 60 * numpy.random.randint(8, 18, size=(students.shape[0],4))).min(axis=1)\n",
    "    # Above = Standard lecture times, uniform distribution, 4 per day\n",
    "    students[:,0] = students[:,1] - numpy.random.randint(15, 45, size=students.shape[0])\n",
    "    students[:,2] = numpy.random.choice(len(bus_stops), size=students.shape[0], p=bus_stop_ratios)\n",
    "    \n",
    "    students_travel = numpy.random.choice(students_travel_chance, size=students.shape[0], p=students_travel_dist)\n",
    "    students_keep = students_travel >= numpy.random.rand(students.shape[0])\n",
    "    students = students[students_keep,:]\n",
    "    \n",
    "    # Return...\n",
    "    return numpy.concatenate((staff, students), axis=0)\n",
    "\n",
    "\n",
    "\n",
    "# Quick test of journeys...\n",
    "journeys = draw_journeys()\n",
    "\n",
    "print('Journeys generated = {}'.format(journeys.shape[0]))\n",
    "\n",
    "#plt.figure(figsize=(12,12))\n",
    "#plt.hist(journeys[:,0].astype(float) / 60)\n",
    "#plt.hist(journeys[:,1].astype(float) / 60, alpha=0.5)\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Students at each stop:\n",
      "Oldfield Park:\t 1611\n",
      "Broughton Hayes:\t 2873\n",
      "Corn Street:\t 2319\n",
      "Pulteney Court:\t 771\n",
      "St Marys Church:\t 2228\n",
      "University of Bath:\t 0\n"
     ]
    }
   ],
   "source": [
    "def draw_queues(scale = 1.0):\n",
    "    \"\"\"Returns a list of queues (deque objects), where they are aligned with the bus stops.\n",
    "    Each student is represented as a tuple (arrived at bus stop, late time)\"\"\"\n",
    "    journeys = draw_journeys(scale)\n",
    "    \n",
    "    order = numpy.argsort(journeys[:,0])\n",
    "    journeys = journeys[order,:]\n",
    "\n",
    "    ret = [deque() for _ in bus_stops]\n",
    "    \n",
    "    for row in journeys:\n",
    "        ret[row[2]].append((row[0], row[1]))\n",
    "    \n",
    "    return ret\n",
    "\n",
    "\n",
    "\n",
    "# Quick test...\n",
    "print('Students at each stop:')\n",
    "for i, q in enumerate(draw_queues()):\n",
    "    print('{}:\\t {}'.format(bus_stops[i], len(q)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_bus_size():\n",
    "    return 85 if numpy.random.randint(25)<8 else 120\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "Journey = namedtuple('Journey', ['interval', # Timetabled gap between buses in minutes\n",
    "                                 'buses', # Number of buses covering route\n",
    "                                 'scale', # Scale factor for population of Bath uni\n",
    "                                 'start', # Start time\n",
    "                                 'end', # End time\n",
    "                                 'passengers', # List of passengers who boarded at each stop\n",
    "                                 'arrive', # Arrival time at each stop\n",
    "                                 'depart', # Departure time at each stop\n",
    "                                 'free_space', # Free space\n",
    "                                 'left_behind', # Passengers left at bus stop\n",
    "                                 'period' # How much time the bus 'covered' for each stop\n",
    "                                ])\n",
    "\n",
    "\n",
    "\n",
    "def sim(interval = 6, buses = 25, scale = 1.0):\n",
    "    \"\"\"Runs a simulation, recording loads of stuff. Parameters:\n",
    "    interval - interval, in minutes, between buses starting the journey\n",
    "    buses - number of buses the company has\n",
    "    scale - multiplier of university population\"\"\"\n",
    "    # Queue of people for each stop...\n",
    "    queue = draw_queues()\n",
    "    \n",
    "    # For recording when a stop is blocked - tuples of (start, end)...\n",
    "    # (kept for final stats)\n",
    "    blocked = [list() for _ in bus_stops]\n",
    "    \n",
    "    # For detecting when they have run out of buses...\n",
    "    busy = [] # List of times when buses will return\n",
    "    \n",
    "    # For recording the statistics...\n",
    "    stats = []\n",
    "    \n",
    "    # Loop and run each bus through the system...\n",
    "    for start_time in range(7 * 60, 22 * 60, interval):\n",
    "        # Return all buses and move to next release if they are all busy...\n",
    "        busy = [return_time for return_time in busy if return_time > start_time]\n",
    "        if len(busy)>=buses:\n",
    "            continue\n",
    "        \n",
    "        # Get buses size and drive times between stops...\n",
    "        size = draw_bus_size()\n",
    "        travel_time = draw_travel_time()\n",
    "        \n",
    "        # Stats...\n",
    "        passengers = [deque() for _ in bus_stops]\n",
    "        arrive = [None] * len(bus_stops)\n",
    "        depart = [None] * len(bus_stops)\n",
    "        free_seats = [None] * len(bus_stops)\n",
    "        left_behind = [None] * len(bus_stops)\n",
    "        block_index = [None] * len(bus_stops)\n",
    "        free_space = size\n",
    "        \n",
    "        # Run the route...\n",
    "        now = start_time\n",
    "        for stop in range(len(bus_stops)):\n",
    "            now += travel_time[stop]\n",
    "            arrive[stop] = now\n",
    "            \n",
    "            pickup = True\n",
    "            for low, high in blocked[stop]:\n",
    "                if low < now and now < high:\n",
    "                    pickup = False\n",
    "                    break\n",
    "\n",
    "            if pickup:\n",
    "                while free_space > 0 and len(queue[stop])>0 and queue[stop][-1][0]<=now:\n",
    "                    passengers[stop].append(queue[stop].popleft())\n",
    "                    free_space -= 1\n",
    "                    now += max(numpy.random.exponential(5 / 60), 1.0)\n",
    "                    \n",
    "            free_seats[stop] = free_space\n",
    "            depart[stop] = now\n",
    "            left_behind[stop] = len(queue[stop])\n",
    "            if stop+1<len(bus_stops) and pickup:\n",
    "                block_index[stop] = len(blocked[stop])\n",
    "                blocked[stop].append((arrive[stop], depart[stop]))\n",
    "\n",
    "        busy.append(depart[-1] + 15) # 15 minutes to get back to depo - not realistic but meh\n",
    "        \n",
    "        period = [None] * len(bus_stops)\n",
    "        for stop, i in enumerate(block_index):\n",
    "            if i is not None:\n",
    "                period[stop] = blocked[stop][i][1] - blocked[stop][i][0]\n",
    "                if i>0:\n",
    "                    period[stop] += 0.5 * (blocked[stop][i][0] - blocked[stop][i-1][1])\n",
    "                if i+1<len(blocked[stop]):\n",
    "                    period[stop] += 0.5 * (blocked[stop][i+1][0] - blocked[stop][i][1])\n",
    "                    \n",
    "            else:\n",
    "                period[stop] = 0.0\n",
    "        \n",
    "        stats.append(Journey(interval, buses, scale, start_time, now,\n",
    "                             passengers, arrive, depart, free_seats, left_behind, period))\n",
    "    \n",
    "    # Return, after adding the time range covered by each bus...\n",
    "    return stats\n",
    "\n",
    "\n",
    "\n",
    "def sub_sim(start = 8, end = 10, interval = 6, buses = 25, scale = 1.0):\n",
    "    \"\"\"Filters the results to only include buses that perform a pickup\n",
    "    in the given time range (hours)\"\"\"\n",
    "    stats = sim(interval, buses, scale)\n",
    "    \n",
    "    start *= 60\n",
    "    end *= 60\n",
    "    \n",
    "    def keep(journey):\n",
    "        if journey.end < start:\n",
    "            return False\n",
    "        if journey.start > end:\n",
    "            return False\n",
    "        return True\n",
    "    \n",
    "    return [journey for journey in stats if keep(journey)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "29969\n"
     ]
    }
   ],
   "source": [
    "def to_output(journey):\n",
    "    \"\"\"Converts a journey to the output format:\n",
    "    interval between buses\n",
    "    for each stop (5x)\n",
    "      excess capacity per minute when it departs (goes negative for people left behind)\n",
    "      average waiting time\n",
    "      passengers on time\n",
    "      passengers late\n",
    "    \"\"\"\n",
    "    ret = [journey.interval]\n",
    "    \n",
    "    for stop in range(5):\n",
    "        ret.append((journey.free_space[stop] - journey.left_behind[stop]) / max(journey.period[stop], 1.0))\n",
    "        \n",
    "        waiting = [max(journey.arrive[stop] - p[0], 0) for p in journey.passengers[stop]]\n",
    "        if len(waiting)>0:\n",
    "            ret.append(numpy.mean(waiting))\n",
    "            \n",
    "        else:\n",
    "            ret.append(0.0)\n",
    "        \n",
    "        before = 0\n",
    "        after = 0\n",
    "        for passenger in journey.passengers[stop]:\n",
    "            if passenger[1]<journey.end:\n",
    "                before += 1\n",
    "            else:\n",
    "                after += 1\n",
    "        ret.append(before)\n",
    "        ret.append(after)\n",
    "\n",
    "    return ret\n",
    "\n",
    "\n",
    "\n",
    "data = []\n",
    "for interval in [1,2,3,4,5,6,7,8,9]:\n",
    "    for _ in range(100):\n",
    "        stats = sub_sim(interval = interval, buses=60)\n",
    "        data += [to_output(journey) for journey in stats]\n",
    "\n",
    "print(len(data))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stats:\n",
      "  Interval: mean=4.52, sd=2.24\n",
      "  Oldfield Park:\n",
      "    Capacity per minute: mean=-1248.71, sd=343.92\n",
      "                E[wait]: mean=0.00, sd=0.00\n",
      "                on time: mean=0.61, sd=5.67\n",
      "                   late: mean=5.38, sd=22.97\n",
      "  Broughton Hayes:\n",
      "    Capacity per minute: mean=-2312.81, sd=777.23\n",
      "                E[wait]: mean=0.00, sd=0.00\n",
      "                on time: mean=0.01, sd=0.85\n",
      "                   late: mean=5.91, sd=24.94\n",
      "  Corn Street:\n",
      "    Capacity per minute: mean=-1758.08, sd=666.63\n",
      "                E[wait]: mean=0.00, sd=0.00\n",
      "                on time: mean=0.09, sd=2.84\n",
      "                   late: mean=5.95, sd=25.04\n",
      "  Pulteney Court:\n",
      "    Capacity per minute: mean=-435.35, sd=191.81\n",
      "                E[wait]: mean=0.05, sd=1.07\n",
      "                on time: mean=2.40, sd=15.98\n",
      "                   late: mean=3.72, sd=19.42\n",
      "  St Marys Church:\n",
      "    Capacity per minute: mean=-1617.65, sd=717.53\n",
      "                E[wait]: mean=0.00, sd=0.00\n",
      "                on time: mean=0.23, sd=4.81\n",
      "                   late: mean=5.99, sd=25.00\n"
     ]
    }
   ],
   "source": [
    "data_mean = numpy.mean(data, axis=0)\n",
    "data_std = numpy.std(data, axis=0)\n",
    "\n",
    "print('Stats:')\n",
    "print('  Interval: mean={:.2f}, sd={:.2f}'.format(data_mean[0], data_std[0]))\n",
    "for stop in range(5):\n",
    "    print('  {}:'.format(bus_stops[stop]))\n",
    "    print('    Capacity per minute: mean={:.2f}, sd={:.2f}'.format(data_mean[1+4*stop], data_std[1+4*stop]))\n",
    "    print('                E[wait]: mean={:.2f}, sd={:.2f}'.format(data_mean[2+4*stop], data_std[2+4*stop]))\n",
    "    print('                on time: mean={:.2f}, sd={:.2f}'.format(data_mean[3+4*stop], data_std[3+4*stop]))\n",
    "    print('                   late: mean={:.2f}, sd={:.2f}'.format(data_mean[4+4*stop], data_std[4+4*stop]))\n"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
