{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Process-Oriented Simulation of M/G/1 Queues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The event-driven approach is more than sufficient to handle the single-server queue with the FIFO queuing discipline. However, the situation may get more complicated using the event-driven approach, when we deal with multiple servers with multiple queues and with more complicated queuing disciplines (such as priority-based scheduling with preemption). A strong aspect of simulus that distinguishes simulus from most other simulators is its natural support for process-oriented simulation. Simulus provides high-level modeling abstractions (such as traps, semaphores, resources, stores, and mailboxes) in order for users to create more intuitive and expressive simulation models using the process-oriented approach. We start here with the single-server queue and later we will expand it for more complicated queuing scenarios.\n",
    "\n",
    "Rather than directly scheduling events (one still can, if needs to), one create processes and have the processes to run their own course. Each process is a single thread of control. A process can be suspended, sleeping for time or waiting for resources. A process can resume execution, after sleep or when the wait conditions have been satisfied.\n",
    "\n",
    "Before we start this section, let's first run the following so that we can set up the environment and import the necessary packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for random distributions, random number generators, statistics\n",
    "import random\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "\n",
    "# for simulation\n",
    "import simulus\n",
    "\n",
    "# for data visualization\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# for animation inside the notebook\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact\n",
    "\n",
    "def exp_generator(mean, seed):\n",
    "    rv = stats.expon(scale=mean)\n",
    "    rv.random_state = np.random.RandomState(seed)\n",
    "    while True:\n",
    "        # 100 random numbers as a batch\n",
    "        for x in rv.rvs(100):\n",
    "            yield x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the single-server queue, we use a process to generate the customers. A process starts from a normal function, in this case, `gen_arrivals()`. Within the function, the process enters a forever loop, in which the process sleeps for some time (which is a random value returned from the `inter_arrival_time` generator function), and then creates another process that starts from a function called customer. As its name suggests, the process represents a customer who arrives at the queue. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_arrivals():\n",
    "    while True:\n",
    "        sim.sleep(next(inter_arrival_time))\n",
    "        sim.process(customer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The customer process waits for a server to be available. A server in this case is implemented as a resource. A process trying to acquire a resource may be suspended if the resource is unavailable (i.e., if it has been already acquired by another customer). In that case, the process will be put in a queue (with unlimited space) maintained by the resource. \n",
    "\n",
    "Once the server becomes available, one of the waiting processes will be selected. By default, the resource is FIFO: it selects the waiting process who has entered the queue first. When the selected process acquires the server, it sleeps for some time representing the customer being served. The sleep time is a random value returned from the `service_time` generator function. After the sleep, the process releases the server so that if other customer processes are waiting for the server, the server will select the next customer in queue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def customer():\n",
    "    print('%g: customer arrives (num_in_system=%d->%d)' %\n",
    "          (sim.now, server.num_in_system(), server.num_in_system()+1))\n",
    "    server.acquire()\n",
    "    sim.sleep(next(service_time))\n",
    "    print('%g: customer departs (num_in_system=%d->%d)' %\n",
    "          (sim.now, server.num_in_system(), server.num_in_system()-1))\n",
    "    server.release()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to simulate the queue. Like before, we instantiate a simulator and then create the two generators for the inter-arrival time and for the service time. We create the server as a resource and then starts the gen_arrivals process. Finally we run the simulation for 10 simulated seconds. The results should be the same as those from the event-driven approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.689679: customer arrives (num_in_system=0->1)\n",
      "1.23285: customer departs (num_in_system=1->0)\n",
      "2.23403: customer arrives (num_in_system=0->1)\n",
      "2.56491: customer arrives (num_in_system=1->2)\n",
      "3.00135: customer arrives (num_in_system=2->3)\n",
      "3.79835: customer arrives (num_in_system=3->4)\n",
      "5.02737: customer arrives (num_in_system=4->5)\n",
      "5.40655: customer arrives (num_in_system=5->6)\n",
      "5.50404: customer arrives (num_in_system=6->7)\n",
      "5.56426: customer departs (num_in_system=7->6)\n",
      "6.39297: customer departs (num_in_system=6->5)\n",
      "7.08135: customer departs (num_in_system=5->4)\n",
      "8.00896: customer departs (num_in_system=4->3)\n",
      "8.99307: customer departs (num_in_system=3->2)\n",
      "9.41727: customer arrives (num_in_system=2->3)\n",
      "9.77588: customer arrives (num_in_system=3->4)\n"
     ]
    }
   ],
   "source": [
    "sim = simulus.simulator('ssq')\n",
    "inter_arrival_time = exp_generator(1.2, sim.rng().randrange(2**32))\n",
    "service_time = exp_generator(0.8, sim.rng().randrange(2**32))\n",
    "server = sim.resource()\n",
    "sim.process(gen_arrivals)\n",
    "sim.run(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's change the service time distribution from exponential to a truncated normal distribution, which specifies a particular range between a and b. We create another generator function for this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def truncnorm_generator(a, b, seed):\n",
    "    rv = stats.truncnorm(a, b)\n",
    "    rv.random_state = np.random.RandomState(seed)\n",
    "    while True:\n",
    "        # 100 random numbers as a batch\n",
    "        for x in rv.rvs(100):\n",
    "            yield x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to run the simulation and show the effect of the service time distribution on the customer's wait time *on the fly*. To do this, we disable the print-outs in the original customer() method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_arrivals():\n",
    "    while True:\n",
    "        sim.sleep(next(inter_arrival_time))\n",
    "        sim.process(customer)\n",
    "        \n",
    "def customer():\n",
    "    server.acquire()\n",
    "    sim.sleep(next(service_time))\n",
    "    server.release()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are wrapping around each simulation run using a function called `sim_run()`. In this case, we can simply call this function to have one run of the simulation and return the results. Since in our example both `gen_arrivals()` and `customer()` use module global variables: 'sim', 'inter_arrival_time', 'service_time', and 'server'. We need to make them global from the `sim_run()` function. \n",
    "\n",
    "To gather the statistics, we use DataCollector provided by simulus. When creating the DataCollector instance, we specify the kind of data to be gathered at the resource. In this example, we would like to keep all the data series for `system_times`, which is just another name for the wait time of all the customers using the resource. When we create the resource, we attach the DataCollector instance. After the simulation run, we retrieve the collected data from the DataCollector instance and return them from the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sim_run(a, b):\n",
    "    global sim, inter_arrival_time, service_time, server\n",
    "    sim = simulus.simulator('mg1')\n",
    "    inter_arrival_time = exp_generator(1.2, sim.rng().randrange(2**32))\n",
    "    service_time = truncnorm_generator(a, b, sim.rng().randrange(2**32))\n",
    "    dc = simulus.DataCollector(system_times='dataseries(all)')\n",
    "    server = sim.resource(collect=dc)\n",
    "    sim.process(gen_arrivals)\n",
    "    sim.run(1000)\n",
    "    return dc.system_times.data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are all set for interactive plotting. We fix 'a' and use a slider to choose 'b'. The `plot_hist()` function calls `sim_run()` to get the wait time of all customers in a list and we plot the histogram accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1421af1687b0438dbe5162437799d597",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.0, description='b', max=2.0, min=0.1), Output()), _dom_classes=('wid…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_hist(a, b):\n",
    "    waits = sim_run(a, b)\n",
    "    plt.hist(waits, alpha=0.5, bins='auto', density=True)\n",
    "    plt.xlim(0, 8)\n",
    "    plt.show()\n",
    "interact(plot_hist, a=widgets.fixed(0), b=(0.1, 2.0, 0.1))\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The wait time shown in the histogram is from one run of the simulation. To be statistically meaningful, we need to run many trials of each simulation setting and show the results with confidence intervals.\n",
    "\n",
    "In the following example, we revise the `sim_run()` function to just collect the summary statistics for wait time and return the mean value (which would consume far less time and memory). This is achieved by configuring the DataCollector with the argument 'system_times' set to be 'dataseries' rather than 'dataseries(all)'.\n",
    "\n",
    "We run the simulation for 25 times at each parameter setting. It is important to note that there's a subtle difference from the previous when we create the simulator for each simulation run. Now the simulator is nameless! This is important, because otherwise the simulator will generate the same random sequence regardless of how many times we run the simulation. Without a name, the simulator will be created with a randomly generated name (using the default random number generator in Python's `random` module). Recall that the simulator's attached random number generator is seeded from the namespace (which we do not change since we use the same global random seed) and the simulator's name. In this case, the simulator would create a different random sequence, which is exactly what we'd need. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sim_run(a, b):\n",
    "    global sim, inter_arrival_time, service_time, server\n",
    "    sim = simulus.simulator()\n",
    "    inter_arrival_time = exp_generator(1.2, sim.rng().randrange(2**32))\n",
    "    service_time = truncnorm_generator(a, b, sim.rng().randrange(2**32))\n",
    "    dc = simulus.DataCollector(system_times='dataseries')\n",
    "    server = sim.resource(collect=dc)\n",
    "    sim.process(gen_arrivals)\n",
    "    sim.run(1000)\n",
    "    return dc.system_times.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code will run pretty slow. Be prepared! We are conducting 30*25=750 simulation runs before we can get the plot. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3hb5fXA8e8rWZL3irMTZwJp2BAKlJZRyibsvUcJUPbeJMywR0gYAcLee5RZoBR+zLAhgRCyp+3Ylixb+57fH/caFMd2FFvyPJ/n8WPpzmNJvkf3nUZEUEoppVri6uwAlFJKdW2aKJRSSrVKE4VSSqlWaaJQSinVKk0USimlWqWJQimlVKs0UXRRxpgjjTHvpOE4xxljPk5HTOlijJlsjDk7A8fdxBjzSbqP217GmHJjTNAY4+7sWHoDY8zJxpg70nQsMcaMbmHdeGPMM+k4T1eniaKLEpEnRGTXzo4j3YwxfYFjgPuc515jzPPGmAXOP+WOrezrM8Y8aIxZaIypM8Z8a4zZo3G9iHwP1Bpjxmf671gXIrJIRPJFJNER53Ney390xLm6GmOMF7gcuDlp2WbGmK+MMQ3O783ScS4ReQ3Y0BizSTqO15VpouhkxpisVJb1IMcBb4hIKGnZx8BRwIq17JsFLAZ2AIqwLwjPGmOGJ23zBHBymmJNSXd7v7pLvG2Mc1/gZxFZ6hzDC7wCPA6UAI8ArzjL0+EpYEKajtV1iYj+pPkHuBj4DagDZgH7J607Dvg/4HZgFXBtK8s+dva5B7ilyTleAc5N8Xwft+FvGAO8C1QDvwCHOMu9wLfAGc5ztxP7lc7zScDzwDNOPF8DmyYd933gqBbOuQTYcR3j/B44MOn5YCAE+FLY1ziveQUQAH4ANnLW+YBbgEXASuBeIMdZt6MT60XYye0xYDawd9Kxs4BKYAtgOCBAlrOuFHgIWAbUAC8n7be38/rWAp8Am6zj6/EYYDmvQRC4MOn8Jzp/z/8a/4Ym+y4A/pH0Pj4LPOq8jz8B45K2HQq86PyNq4CpyZ8357WrAeYDeyTtNwh41flczQVOSlrX+Nl53Hk//ukse85ZVue8R+sDlzjv22Jg16RjzAAuT3q+K7AUMEnLFgG7p/h6CnAmMA+owr5TcSWt3w6Y3xnXmY786fQAeuIPcLDzD+ECDgXqgYHOuuOAOHCGczHJaWVZY6LY3vmHMM7zEudCMCjF832cFNv32Beh5n7udrbJc853vBPP5s4/yVhn/UbOReBPwGXAZ4DbWTcJiAEHAR7gfOdi4XHWVwJbtfC6rVOiAPoDYWBMk+UBUrjAArsBXwHF2EnjT0mv2+3YF7RSoAB4DZjsrNvReb9uxE4oOcCVwBNJx94LmO08Hs7qieLf2Im0xHmNdnCWb4598dsaOwEfi33x9jnrX2/lvXs96dwLcC74Tc7/qPPe5pBaoggDezqxTAY+c9a5ge+c1ygPyAb+mvR5iwEnOdudip0QGz+7/wPudvbZzPk8/L3JZ2c/7M9yTlIcu2F/Fh/F/jxd5rx2J5F0oQa+BA5Oen4O8GaTv/N14LwUP2MCfOB8DsqBOcA/k9aXOtsUdvZ1J5M/nR5Ab/jB/oa4r/P4OGBRk/UtLWtMFAb7W9D2zvOTgPfX4XzrdEeBnWw+arLsPmBi0vPzsO80aoD1kpZParygOM9dwHLgb87zGE0u7EnbppwonIvEf4D7mlm3tPG1Wssx/u7842/D6t8SDXayHZW0bNvGCxL2RTYKZCetH439jTfXef4Ef9xlDXcuJlnAQOxv/CXNxHMPcE2TZb/gJJJ1eP8W0HyiGJm0bEfWnij+k7RuLBBKei0qcRJfM5/buUnPc51zD8C+C0kABUnrJwMPJ53zf02ONwl4N+n5eOw7pcYvJgXO8Yud57+SdLcAXAE83eSYTwCTUnwtpcnx/gW81+RzKED5urxH3e1H6ygywBhzjFPRWmuMqcX+Bl6WtMniZnZrbhkAYn8inwYOdxYdgf1hT/V862oYsHXj8ZxjHon9z97oEWe7N0Tk15b+FhGxsBPAIGdRDfY/d5sZY1zYRSxR4PRmNinA/pbdKhF5H5gKTAMqjDHTjTGFQF/sC9xXSX//W87yRpUiEk461lzs4qfxxphcYB/gyWZOOxSoFpGaZtYNA85r8roP5Y/Xrr1a/Iy1ILnOqAHIduoNhgILRSS+tv1EpMF5mI/9d1SLSF3Stguxiwtbi3Fl0uMQUCV/NAxorOvKd343/XwFgcImxyvETuqpSo5pIau/H43nWuvnrTvTRJFmxphhwP3YF7A+IlIM/Ij9LbWRNLNrc8uSPQUc5Bx/a+CFdThfcnw/OU01m/u519lsMfChiBQn/eSLyKlJh7ob+xZ+N2PMX5ucZmjS+VzAEOziB7CLvtZfy9/aImOMAR7ELnY6UERiTdYPxq5H+SWV44nIFBHZEvsb8/rABdjFbCFgw6S/v0hE8pN3beZwT2En832BWU7yaGoxUGqMKW5h3XVNXvdcEXnK+dvebOW9e3MtsTVdXo+dDHGO7Wb1RNiaxUB5Gyqbl2H/7ckX8nLsO8DmYmyLpp+vn4BNnM9No02c5akamvS4nD8+y2AXVy4QkcC6BtqdaKJIvzzsD3slgDHmeOxv+O0iIt9gX8AeAN4WkcZvMOt0PhHZ0LnoN/dzirPZ68D6xpijjTEe52crY8yfnHMcDWyJXcxwJvCIMSb5IrqlMeYA50JyNhDBrscAeAO71dLvnGav2c5TrzEmu8k/drJ7sP85x8vqLaca7YBdLBdxjn2cMWZBcwdy/qatjTEe7AtnGLCcu6D7gduNMf2cbQcbY3ZrIaZGT2NXnp5K83cTiMhy4E3gbmNMifPabu+svh84xYnJGGPyjDF7NV5YRWSPVt67PZJOsxIYuZZY52DfIezl/P2XY9e3pOIL7OLEG5wYs40x261tJxFZjF1BP9nZZxPsCvbHUzxvKpp+vv6LXdx1pvM5a7wDfR9a/3wkucB5r4YCZ2HXLzXaAfv97NE0UaSZiMwCbgU+xf6H3Ri7VVA6PAn8g6SLUCbO5xQN7Aochv3taQVOxa0xphy4AzhGRIIi8iQwE7tis9Er2PUcNcDRwAFJ3/wfBfY0xuQkbf8L9jf4wcDbzuNhAMaYSxu/LTt3TydjV4KuSPo2fWTSsY7EbqHUaGgrr0ch9sW5BrtIYRV/tL+/CLtVzmfGmAB2fcgGLb5o/J4EPgX+wuoXk6aOxq6r+Rm78vpsZ/+Z2PVPU52Y5mIn43U1GbjcKb46v4VY/djl7Q9gf6Ovxy4iXCun2Gc8dr3MIme/Q1OM7XDsOpNlwEvY9V7/SXHfVLwGjDHGDHJijWJXjh+DXTx0ArCfsxxa/3w0egW70cO32A0RHkxadzhOn6CerLElglJpYYyZBIwWkaNa2eZ6oEJE0tJ7Num4m2BXbm+btOwd4CwRmZ3Oc6muyxgzAbuF3lp7/7fn82Hsjp1Hi8ghbQizW9FEodIqlUShlOpetOhJKaVUq/SOQimlVKv0jkIppVSrusXgYMnKyspk+PDhnR2GUkp1K1999VWViKTaV2Y13S5RDB8+nJkzZ3Z2GEop1a0YYxa2dV8telJKKdUqTRRKKaVapYlCKaVUqzKWKIwxQ40xHxhjZjkD0Z3VzDY7GmP8zsin3xpjrsxUPEoppdomk5XZcezJQb52BjX7yhjzrjM2UbKPRGTvDMahlFKqHTJ2RyEiy0Xka+dxHfZY/YNb30sppVRX0yF1FMaY4djTPH7ezOptjTHfOWPtb9jC/hOMMTONMTMrKyszGKlSSqmmMp4onHkKXgDObmZyj6+BYSKyKXAX8HJzxxCR6SIyTkTG9e3bpv4iSiml2iijicKZEOUF7EnnX2y6XkQCIhJ0Hr8BeIwx7ZnCUymluhR/fYRVgfDaN+zCMtnqqXHKytkiclsL2wxonMnMGPNnJ55VmYpJKaU6kmUJ85bV8ctiP8FQbO07dFGZbPW0HfZMXj8YY751ll2KPecsInIvcBBwqjEmjj2r2WGiw9kqpXqIKn+YUCxBni+LnxfVsPHIPvg87s4Oa51lLFGIyMdAS/MeN24zFXvaR6WU6lHiCYtFFXUUZHvwZLmIhSzmLvEzZlgJblerl8YuR3tmK6VUBlT6w8QSgifLvszm53gIhGIsXFFHdys40UShlFJpFos7dxM5ntWWF+d5WV7dwMqaUCdF1jaaKJRSKs1W1jQgFmS5V7/EGmMozvfx2/IA/vpoJ0W37jRRKKVUGkViCRZX1lOY6212vdtlKMj28MviWsLReAdH1zaaKJRSKo2WVdXjdhlcrVRYez1uslwuflnsJ56wOjC6ttFEoZRSaRKKxFle3UB+k7qJ5uRmZxGKxJm3PNDlK7c1USilVJosqarH43bjMqk1fy3O91Hpj7C0qj7DkbWPJgqllEqD+lCMytoQ+Tnr1j2tJM/LwpVBqrvwMB+aKJRSKg0WVQbxedyYFO8mGrlchsJcL3OW+KkPd81hPjRRKKVUOwUaotTURcjLXnvdRHM8WS58Xje/LPITi3e9ym1NFEop1Q4iwoIVQXJ87RsRKcebRUIsvp1bxbzlAWqDkS7TIiqTgwIqpVSPVxuMEgxFKS3IbvexCnK8JCyhOhBhZU0IYwyFOR76FGVTmOsh27vuRVvpoIlCKaXayLKEBSvryPO1rcipOW6X+b15rYgQiVnMXx5AELxZbsqKcijO95Kf4+mwwQU1USilVBtV14UJRRKUFvgycnxjDNleN9lee2jyeMJiZU0Dy6qCGJehb1EOw/rnrzFUSLppolBKqTZIWMKClcE1Bv7LpCy36/ehQSwRqvwhAg1RxgwtbncdSWu0Mlsp1euJCKv8IaoDYYKhGNFYYq29patqQ8Ri1u/DiHc0lzEU5fkQge/nVVNTl7l+GHpHoZTq1eIJi/kr6qioCf0+PpPB7t+Q68siPyeLvGwPPq8bb5Ybn8dFwhIWVQQpyO24u4mW5PqyiLktZi+sZdiAfAb1yUt7hbcmCqVUrxWKxJmzxE8oEqe0wLfaBdYSIRa3qPJHWFGdNH+EAW+Wm4SsOYx4Z/FkuSjO97FwZT314QQjBxakNTZNFEqpXqmmLsycJQE8bvsi25TLGHwe9xpzXIsI8YSQ4+1ac1+7XIbSAh81dWF+jMTZYEhR2uotukY6VEqpDmJZwpLKILMW1pLryyI3e90upsYYPFmuTunPkIqiPB9WQtJab6GJQinVa8TiFr8u9bOoop6SfF+nVURnWm52Frm+LGYvrGVpVbDdw5hr0ZNSqleoD8eYs9hPLGFlrN9DV9JYb7FgRZD6cALsOvo20UShlOrxVvlDzFkaINvrbnGK0p7I5TL0KcymOhDG5cnJbetxNFEopXosyxIWVwRZUlVPUZ63y7RS6mjF+T6My9Xm2ndNFEqpHimesJizxE9tMEJJgS/lWefUmjRRKKV6nIQlzF3qJ1CfnlFde7veeR+mlOqxLEv4bVmAmmC02f4Rat1polBK9Rj2JEIBVvlDlGiSSBtNFEqpHkFEWLQyyPKasN5JpJkmCqVUj7C0qp4lVfWU5Hu7bK/p7koThVKq21te3cDClUFt3ZQhGUsUxpihxpgPjDGzjDE/GWPOamYbY4yZYoyZa4z53hizRabiUUr1TBU1IeYtC1Ccr0kiUzLZPDYOnCciXxtjCoCvjDHvisispG32ANZzfrYG7nF+K6XUWq0KhPl1mZ+iPG+HzR/dG6V0R2GM+asx5njncV9jzIi17SMiy0Xka+dxHTAbGNxks32BR8X2GVBsjBm4Tn+BUqpX8tdH+WWxn6Lc3tvjuqOs9dU1xkwELgIucRZ5gMfX5STGmOHA5sDnTVYNBhYnPV/CmslEKaVWU9cQY9bCGgpyPF06SdTUVHP7TVczb+6czg6lXVJ5hfcH9gHqAURkGVCQ6gmMMfnAC8DZIhJoS5DGmAnGmJnGmJmVlZVtOYRSqoeoD9tJIteX1aWHCQ+HQ5x24uFMn3Y7+++5PfffcwfxeLyzw2qTVF7lqNiDmQuAMSYv1YMbYzzYSeIJEXmxmU2WAkOTng9xlq1GRKaLyDgRGde3b99UT6+U6kFicYuV1Q38OL+m2ZnnuhLLsrjo7FP49usvmXT9bWy/4y7cdsNVHHHAbsz5ZdbaD9DFpJIonjXG3Iddf3AS8B/g/rXtZOyGzA8Cs0XkthY2exU4xmn9tA3gF5HlKcaulOoFGiJxFq6o46s5lcxbUUeuL4vsLjYNaVM3XXsF77z5Khdefi2HHnk8U+57lFunPsiSxQs5aK8duWfKLcRisc4OM2UmlZmPjDG7ALtiT3zxtoi8m8I+fwU+An4ALGfxpUA5gIjc6ySTqcDuQANwvIjMbO2448aNk5kzW91EKdXJ6kMxYgmLbG8WPs+6TxsqIgQaYiytqqc2GMXjNuRle3B1g5ZNj864l8lXXcLRJ5zCpRMnr7auelUV1068kDdfe4k/jd2Ya2+eytiNNumQuDZaf8iv8XBw/bbsm1KiADDGFJLUnFZEqttywvbSRKFU11UfirGkqp5VgTDGGAzgdhkK870U53rJzfaQ43O3WAEdT1isCkRYVlVPKBon22tP6dldvPPma5x96rHsvOte3HHPw7jdzd/5vPvW61x9+XnU1lTzz1PP5tQzzsfrW/uwI6uqKvnmqy/4ZuZnfPfNTGLxOKWlfSjt05fSPn0oLS2jtKzM/t3H+Sktw5edndlEYYw5GbgKCGPfGRhARGRkW07YXpoolOp6fk8Q/jBej5u87Kzf7yIsS4jEEkTjCSyntjMnO4uSPC8FeV5ynGKkitowy6vrEYH8bE+XrqhuzjczP+f4I/ZjzNiNefjpV8jOzml1+9raGiZfdQmvvvgMo9cfw/W3TGPjTf/oc2xZFr/N/YVvZn7ON199wdczP2fRgnkAeLxeNtxoU3Jy86iprqJ61SpqqqtaLM7Kyy8gHLMymih+BbYVkaq2nCDdNFEo1XUEQ3bx0Cp/GJ/XTa4vK6VipljcIhJLEEs4pdICbrchv5sULzW1YP5vHL7/rhQWFfPUi29T2qcs5X3/+97bTLr0XCorVnDMCadSVFzMN199wXdff4nfXwtASWkfthi3NZtvuTVbjNuasRttii979Xk2RIS6QIDq6ipqVlWxalXV74+rq1fxxOOPZjRRvAUcICINbTlBummiUKrzBZ07iJpA4x2Ep7ND6jSrqio54oDdqKsL8NRL7zBs+LoXttQF/Nx03ZU8//SjAIwavQGbj7OTwubjtmbY8JHtHugw00VPmwMPYXeWizQuF5Ez23LC9tJEoVTnCYZiLK4MUlMXwdfLEwRAKNTAcYftw5yfZ/Hw06+y6ebj2nW8RQvnU1hUTHFxSZoi/EN7EkUqtUT3Ae+zeuslpVQvIiIsXBlk2ap6fB63Ti8KJBIJLjjzJH747mum3PdYu5MEQPmwtY6O1ClSSRQeETk345EopbqsFTUhllXVU1Lg65FzPcyfN5dPP/4vg4eUUz58JIOHlOP1elvcXkSYfNUlvPfOG1x+1Y38Y7e9OjDajpdKonjTGDMBeI3Vi546pXmsUqpj1QYjzFtuD+PdE5PEB++9xQVnTqA+WPf7MpfLxaDBQykfPpLyYSMYNnwk5cNHUD5sJEPLh/Pkow/wxCP3c/xJp3PkcRM6MfqOkUqiONz5fUnSMgE6pXmsUqrjNETi/LzYT2FOzxvGW0R44N47uf3Gqxm70aZcf+s06oNBFi2cx6IF81m0YB6LFs7nzddfwl9b8/t+xhhEhN333p/zL72qE/+CjrPWRCEiXbPQTCmVUdFYgl8W1eLLcnW7Pg1rEwmHueLis3jtpWfZY/z+XHfzVHJycgHYfMs/r7F9bW0NixfayWPhgvkAnHjyGbhcPet1aUmLicIY83cRed8Yc0Bz61sY5E8p1QMkLOHXpX7iCYuC3JbL6rujipUrOOPko/n+m5mced5lnHLGeWstUisuLqG4uGS1DnG9SWt3FDtgt3Ya38w6ATRRKNUDiQgLV9QRaIhRkr/2YSW6kx+//4bTTzqKuoCfKfc9xi67793ZIXULLSYKEZnoPLxaROYnr0tlhjulVPe0oibE8uoGSgt6VpJ447UXufS80+jTty9PvPAmY8Zu3NkhdRupFLC90Myy59MdiFKq89XUhZm/LEBJD2rhZFkWd95yHeedfiIbbbIZz77yniaJddRaHcUYYEOgqEk9RSGgvW2U6mHqwzF+XuynINfbLcdbak59fZCLzzmV/7z9OgceejRXXntLq/0jVPNaq6PYANgbKGb1eoo64KRMBqWU6lh2Cyc/2V53t27hFAmHCQT8BPy1rKqq5PqrLubXX2ZzycTJHH38yT3mLqmjtVZH8QrwijFmWxH5tANjUkp1oIQl/LrET0IsCrrwt+0Vy5fyzpuvsWjBPAL+2t8TQsDv//1xJBJebZ+CwkLue/hZ/rrDzp0Udc+QSj8KTRJK9VAiwoIVAQKhrtnCqaqygrffeIU3X3+Jr76wL0WFhUUUFhVTUFhEUVExo9Yb4DwuosBZV+SsHzN2Y/r1H9DJf0X3132mjlJKpd3y6hArqkNdqoVTbW0N7775Gm++/hKff/I/LMti1HpjOPO8y9hj/P4MHzGqs0PsddaaKIwxI5prHtt0mVKqe6lriLFgRV2XaOEUrAvw/rtv8sZrL/J//3ufeDxO+fCRTDjtXPYYvz/rbzC2U+Pr7VK5o3gBaNod8Xlgy/SHo5TqCLG4xS9LasnzZXVaC6dIOMyH77/D6688z4cfvEM0EmHgoMEcc+Kp7Dn+AMZutGmnJ7CuzrKkQ94/bR6rVC8jIsxfUUciIeRnuzv03IlEgi8+/ZjXX36Od956lWBdHWV9+3HI4cey5/gD2HSLrXrN+EntYYngD0ZwuVxYYuF2uTI6jaw2j1Wql6moDVPlD3XY5EMiwqwfv+O1l57jjddepLJiBXn5Beyy+97sve9BbP2X7cnK0urSVMUTFv76KEP75TO4LI9QJE6VP0RFbZh4Qsj2uMnxudN6N6bNY5XqRerDMeYtD1DYAQP9LVwwj9dfeZ5/v/I883/7FY/Hw/Y77cJe+x7ETv/YnezsnIzH0NOEonHC0QRjhhbRp8h+/fJzPOTneBjar4C6higra0JU19lTB+V6s/B523/X2FrR04UichNwhDHm8KbrO2vObKVU28QTFr8uDZDtcZPlzlzxztdffsbN11/Jt19/CcBW22zHcSedxq577JORuaB7i0BDFI/bxSYjS5udq9ztMhTn+yjO9xGNJfDXR1lRbSeN9r7drd3vzXZ+z2zfKZRSXcGilUHCkTjFGeovURfwc9uNV/P04zMYOGgw510yib32OZCBg4Zk5HzdRTxhAbQ5OVsi1AYjlBb4GDWoKKWe816Pm77FOfQtziEUiVMTjCKJeKxNAdB60dNrzu9H2npwpVTXsCoQzuiIsP95+99cc8UFVFWu5JgTT+XM8y4lLy8/I+fqDixLaIjEiSUsslwGjCEej5HldpGXnXpLs1jcItAQpdypj2hLZXWOL4scXxZWPBJa550dqfSj6AtcBIwlqbWTiPy9rSdVSnWccDTO3KV2vUS6m5tWrFzOtVdexLtvvcb6Y8Zy1/TH2GSz3tlyXkQIRRJEYglcbuhblENZUTb5TjFRfTjOqoBd6ZywBG+WixxfFq4W3pNQJE4kbjGmvJg+hZ3b0DSVpgZPAM8AewGnAMcClZkMSimVHpYl/LY0gNtl0jrYn2VZPPfUI9x6wyQikQjnXHgFx084A49nzbLzni4cTRCKxgEoLfAxYmABBbmeNYqaCnI9FOT+Uelc6Q+xyh9BRMj2ZZHt+aOlkr8+itdj10fk+jq/RVgqEfQRkQeNMWeJyIfAh8aYLzMdmFKq/ZZW1RNoiFKSxqaw8+bO4cpLzuarLz5l623/xqTJt3f7YTUaInHC0TjGGJp+vxfAAMYYXMbgctmPo7EEIlCY52FI30KK8nwpJePkSucRAywCTqVzbTCKMZAQ6FvoY+Sgwow2OlgXqSSKxgqQ5caYvYBlQGnmQlJKpYO/PsriiiDFaaqXiEaj3H/PHdw39VZyc3K59qa7OOCQI7t17+lwNEFDJEZRrpfRg0vxZrmwxL4TS1iCZQkiQkKEWFxIJCzilkU8LhT0yaU430d2O5qfZrldlBZmU1qYTSSWoDYYRUToX5LTpV7XVBLFtcaYIuA84C7sntnnZDQqpVS7RGMJfl1SS16Op8Uy8HXx86wfuOCsCcyd8zN7jN+fSyfeQFnffmmItHNEYgnqwzHyc7xsOLyUwlxPp1+YfR43/Uu6Zt+SVBLFf0QkDPiBnTIcj1KqnRqH6BCxLz7t9frLz3HFRWdRUFjE3TOeYqedd09DlJ0jFreoC8XI8bn5U3kJxfnpr+DviVJJFD8aY1YCHzk/H4uIf207GWNmYA8BUiEiGzWzfkfgFaBxFNoXReTqVANXSq3JsoQVNSGq/OF2t5SJxWLcfP2VPDbjXrb887bccffD3fYuIp6wqGuI4vW4WX9IIaUF2T1muteOkMrERaONMeXA37BbPk0zxtSKyGZr2fVhYCrwaCvbfCQie6carFJqTfGERTAUY1UgwqpAmHjCanenuqrKCs49/QS+/Oz/OPqEU7jg0qu7ZYumeMK+g8hyG0YOKqKsKBu3Joh1lko/iiHAdtiJYlPgJ+Djte0nIv8zxgxvZ3xKqWaEowm7iWVtGH9DFAG87vSMIPrdNzM565Rj8dfWcOPt97LPAYemJ+gO1Jg8XS7DsP4F9CvO7jItiLqjVIqeFgFfAteLyClpPv+2xpjvsFtSnS8iPzW3kTFmAjABoLy8PM0hKNX1iQgN4Ti19VEq/WFCEbvdfrbXTXFe+srZn33qEa698kL69RvAky++zZ823Dgtx21NYwujdPTzaLyDcLmgXBNE2qSSKDYH/oo9OODFwK/AhyLyYDvP/TUwTESCxpg9gZeB9ZrbUESmA9MBxo0bJ+08r1LdSjiaYPbCGsKxBC5jyPVlpX1+62gkwjVXXsjzTz/KX3fYmZun3J+xAfwsESJRuwezCLjdhiy3IRiOOn+fZx3znF0AACAASURBVJ2TRmOCcOsdREakUkfxnTHmN+A37OKno4AdgHYlChEJJD1+wxhztzGmTESq2nNcpXoSEWHBigAJS9KeHBotX7aEs089ju+//YqTTzuXM867FLc7fRMaiQiRmEU4GkcQjDGU5PsY3DeP/BwPOV67R3JDJE5tMMrKmgaCdVFcLhe5vqxWk0Zyghg+oIC+RZogMiGVOoqZgA/4BLvV0/YisrC9JzbGDABWiogYY/4MuIBV7T2uUj1JRW2YVYFIxsb6+eLTjznntOOJRCJMue8xdtk9PW1L4gmL+lAcC8EYKMzxMrC0wE4MLUy/muvLIteXxcDSHELRBDV1ESpqQr/XNSQnjeRK6hEDCijTBJFRqRQ97SEi6zy2kzHmKWBHoMwYswSYCHgARORe4CDgVGNMHAgBh4mIFisp5QhF4sxbHqAoL/2TDFmWxYzpd3HHTddQPnwkU6c/zsjR66fl2AlLCNRHKe9fQFGel9zsrHVqaWSc4rVcXxaD+uQ6dxoRKmrCBEMxBLs+Y+SAAvpogugQqRQ9tWkAQBFZY7KjJuunYjefVUo1YVnCvGUBvFmutF8IK1au4OJzT+XTj//Lrnvsw3U330V+QWFajm1Z9twJowcV0r80t93HM8aQl+0hL9vDoD55NETiRGIWhc0Muqcyp/OHJVRKrWFFTQh/QzTt81p/8N5bXHb+6YQaGrjqhjs4+LBj0tZiSkSoCUYY1j8/LUmiqT+SRtoPrdai1ZRsjHEZY/7SUcEopaA+FGPhijqK89JXeR0Jh7lu4kX864TD6T9gIM+//gGHHH5smpNElCFleQwuy0vLMVXX0WqiEBELmNZBsSjV6yUsYe6yANled9qGmJg752cO2XdnHn94OseceCpPv/Quo9bbIC3HblRbH6VfcTbl/fN17KQeKJVCvveMMQcaffeVyrilVfU0ROLkpGGyGhHh6cdncNDeO1FVWcG9Dz3DJVdejy87vWU3tcEIJfleRgws1CTRQ6XyaTwZOBdIGGNC2HN4iIikp/ZLKQVAoCHK4spgWvpL1NRUc8VFZ/Le2/9mu+3/zuRb76Zvv/5piHJ1dQ1RCnI8jB5cpGMo9WCptHoq6IhAlOrN4gmLuUsC5PvaP3/EF59+zIVnn0z1qkouuvxajjnxVFyu9LcQCoZi9misQ4u1BVIPl0qHOwMcCYwQkWuMMUOBgSLyRcajU6qXWLQySDSeWKdRX0WEqsoKFi74jQXzfmPB/LnMmzuH/773NsNGjOKpB95lw403zUi8DZE4LmPYYGhxWufiVl1TKkVPdwMW8HfgGiCIXcG9VQbjUqrXqA1GWFHT0GKRUzgc4tdfZrNg3lwWLpjHgnlzWTD/NxbM/436YN3v23m8XsrLR3DksSdx9oVXkJeXn5F4I9EE8YTFxiNK2zUNqOo+UkkUW4vIFsaYbwBEpMYYk/6uokr1QtFYgl+X+snPbn4qzpUrlnHUQXuyZLE9ao4xhsFDyhk+cjSbbbEVw0eMZvjIUQwfMZqBg4ekdYymluINReNsNKI0LRXuqntI5Z2OGWPcgAAYY/pi32EopdpBRFiwMogIeJuZsrSmehUnHnkANTXV3HzndMaM3Zih5cPT3mopVfGERX04ztjhJeTndL9JjFTbpZIopgAvAf2MMddhj9F0eUajUqoXqA6EqawNUVqwZpFTfbCOk487hMWLF3D/oy/w522264QIbeFogoZonCyXYYOhRRkZe0p1bam0enrCGPMVsDN209j9RGR2xiNTqgerD8WYu6yOwtw1Jx2KhMOc9s8jmfXjd9w1/bFOSRKWJQRDMRJikZ/tZYPBhRTl+7R1Uy+VaiHjSuwhxrOAHGPMFiLydebCUqpnEhEqasPMWxYg2+deo8VQPB7nvDNO5PNPP+KmO+5jp3/s0aHxRaIJGiIxXG7DgNI8yop85GVrMVNvl0rz2GuA47AnLmocBlywW0EppVIUT1jMX1FHRU2IojzvGt/OLcvi8gvP4L133uDyq25k/P6HdEhcliUEwzESCYu8HA/rDSmiWO8eVJJU7igOAUaJSDTTwSjVU9WHYsxZ4icaT1Ba4FujuElEuOHqS3nlhac549xLOfK4CRmPqXHyH2NgQEkOfYtyyNNKatWMVBLFj0AxUJHhWJTqcUSEipoQvy0PkOPLoqiFEWHvnnIzjz10H8ee+C9OPfP8jMZkiVDXEAV0+lCVmlQSxWTgG2PMj0CkcaGI7JOxqJTqAWJxu6ipyh+iKM/X4lhIjz10H1Nvm8z+Bx/BhZdfk9GB9YKhGLG4xaA+uQzsk9tss1ylmkolUTwC3Aj8gPafUColwVCMOYv9xBJWq5MPvfLi01w/6WL+sdveXH3DnRkZkwkgFI0TisTpU5jN0H755GpnObUOUvm0NIjIlIxHolQPICKsqAkxf3kdub4sirJb7nPw/rtvcNn5p7PNdjtwy5T7ycpK/8U7FrfrIfKys9hoRCmFudoHQq27VD6ZHxljJgOvsnrRkzaPVSpJNJZgwYo6Kv1hivNbLmoC+PyTjzjntBPYcOPNmDr9sbT3tk5YQqAhiifLxQZDCikpyE7bREiq90klUWzu/N4maZk2j1XKkbCEqtoQC1YGMQb6FLZ80RcRnnniIW685nLKh43g3oefJS+/9ZH8RYRwLIEIICBOK3WR5G3+WG5ZgjGGYf0L6FesFdWq/VLpmb1TRwSiVHfkr48yb3mAcCROQe6afSOSVVVWcPmFZ/Dh++/w1x12ZvKtd1NSUtrq8RvCccLROMUFPnweF/bgCOBygcFgAONyfhv7x+1yUVLgw6cV1SpNUulwVwRMBLZ3Fn0IXC0i/kwGplRXForEWVQRpMofJj/bQ0krFdYAH7z3FpdfcAbBYB2XTbqRI479Z6sV15FogvpIjOI8L+sP7aOD8KlOlUrR0wzsvhSN3USPBh4CDshUUEp1VfGExbJVDSypDOLNcrdazATQ0FDPTdddwTOPP8SYsRvx8NOvst76f2px+8bK5xyfm7HDSijKW3MsKKU6WiqJYpSIHJj0/CpjzLeZCkiprsiyhCp/mIUr60hYQnG+b61Tlv74/TdccNYEFs7/jRNPOZMzz70Ur6/5DncJy+4E53YbRg8qpKxIK59V15FKoggZY/4qIh8DGGO2A0KZDUupriPQEGX+8joawnHyczxrnfozkUhw/923M+2OGynr24+HnnyFrf/yt2a3tcQepdUSGNovj/4luVr5rLqcVBLFKcCjTl0FQA1wbOZCUqpriCcsFq0MsqKmgRxfFiXNzBvR1JJFC7nonJP5eubn7LnPAVx57a0UFRU3u20wFCMatxhYmsOgsjytfFZdViqJIiAimxpjCgFEJGCMGZHhuJTqVOFonF+XBqgPxSjJX3MQv6YSiQSvvvgM1026GGPg5juns/d+B7e4fU1dhKJ8L2MHFGgvadXlpfIJfQHYQkQCScueB7bMTEhKdS5/fZRfFtfidhmK81u/i6hYuYIXn32c555+lGVLFjNu679ww233MHhIebPbiwg1wQhlRTmMGlTYaqc8pbqKFhOFMWYMsCFQZIxJbuFUCHTOpL1KZZCIsGxVPQtWBCnI8bQ4YJ5lWXz2yf945vGHeP/dN4jH42yz3Q5ccOnV7LL7eNzu5vdrTBL9inMYObBQK6tVt9HaHcUGwN7YQ4yPT1peB5yUyaCU6mjxhMW8ZQGqAhFK8n3NXsSrV1Xx0nNP8uyTD7No4XyKS0o55oRTOPiI4xg+YlSrx7eTRJQBJTkMH6BJQnUvLSYKEXkFeMUYs62IfNqBMSnVoerDzqRCMYvSJhXWIsKXn3/Cs088xDtvvUYsGmXc1n/hjHMvYZfdx6c0RpOIUB2MMrAkh+EDCjRJqG4nlTqK/Y0xP2E3iX0L2AQ4R0Qeb20nY8wM7DuSChHZqJn1BrgT2BNoAI7TgQZVR1vlD/Hr0gA+j5uivNVHVv38k4+45ooL+G3uLxQWFnHYUSdwyOHHMnr9MSkf3xKhpi7C4LI8hvXP185zqltKpcH2rk5F9t7AAmA0cEEK+z0M7N7K+j2A9ZyfCcA9KRxTqbSwLGHRyjp+XuwnP8dDTlLLo0QiwbQ7b+KEI/cjYSWYfOvd/PeLWVw6cXKbksQQTRKqm0vljqJxkJm9gOdExJ/KB15E/meMGd7KJvsCj4qIAJ8ZY4qNMQNFZHkKMSnVZpFYgt+W+fEHY5QUrN7DuqqyggvOmsBn//ch4/c/hInX3UpeXv46n6MxSQztl8/QvnmaJFS3lkqieM0Y8zN20dOpxpi+QDgN5x4MLE56vsRZtkaiMMZMwL7roLy8+WaHSqWipi7M3GUBDGaNDnSff/IR5595EnUBP9fcOIUDDz2qTRd4y7JbN5X3y2eIJgnVA6y16ElELgb+AowTkRhQj3030GFEZLqIjBORcX379u3IU6seIp6wmL88wOxFtWR7slYbjTWRSHD3nTdzwpH7UVBYyLOv/oeDDju6XUliWP98hvbT4ibVM6QyzPgxSY+TVz3aznMvBYYmPR/iLFMqrYKhGL8u8ROJW2v0sq6qrODCs0/m04//y/j9Dmbi9be1qagJ7GTkr48yfEABg8vy0hW+Up0ulaKnrZIeZwM7A1/T/kTxKnC6MeZpYGvAr/UTKp0sy+5At6iinlxfFsVNWjV98enHnHfGP9tV1JSwhIZwnHjCIsttGD24iP4lOen8M5TqdKnMcHdG8nNjTDHw9Nr2M8Y8BewIlBljlmBPfuRxjnkv8AZ209i52M1jj1/H2JVqUSgSZ+6yAMGGKMV5q3egSyQSTJ92O1Nvn8ywEaN48PEXWH/Mhikf27KE+nCcuGWR5TL0Lc6mT2EOedlZ2kdC9UhtGY2sHljroIAicvha1gtwWhvOr1SLRISK2jDzlwfwZrnXmHmuqrKCi845hU8++mCdiposS2iIxInFE7jdLsqKsikryiY/26PJQfV4qdRRvAY0TuPuAsYCz2YyKKXaIhJLMH95gFWBCEV5a85f/cnH/+Wis0+mLhDg6hvuTKnCOhSNE44mcLmgb1EOfQqzyc/x6GB+qldJ5Y7ilqTHcWChiCzJUDxKtUl1wGn2alhjetJYLMZdt13PA/fcychR6/Pg4y+utajJEsFfHyXXl8XYYSWaHFSvlkqiWAQsF5EwgDEmxxgzXEQWZDQypVK0orqB35YFKMz1rjH73NLFizjvjH/y3TdfcvDhx3LJxOvJyclt9XiRaIJgOMbQfvkMLsvTBKF6vVQSxXPY/SgaJZxlWzW/uVIdI3lY8OJ83xoX9Lf+/TJXXnwWIsJt02awx977r/V4/vooniwXG48spTDX2+r2SvUWqSSKLBGJNj4RkagxRv+DVKcSERZXBFlcWb/GsODhcIgbrr6UZ554mE0225Jb73qQIeXDWj1eLG4RaIgysDSX8v75Om+1UklSSRSVxph9RORVAGPMvkBVZsNSqmWWJSxYUcfymtAaYzX9Omc25552AnPn/MyJp5zJWedfjsfjaeVoUNdgfw8aO6x4jVZSSqnUEsUpwBPGmKnO8yXA0ZkLSamWJSzht2UBVvlDlCb1shYRnnvqESZfdSl5+fk88NgLbLf931s9VmNP6r5F2QwfUNDijHZK9XapdLj7DdjGGJPvPA9mPCqlmhFPWMxd6qcmGF3tm3/A72fiJWfz1r9f5i9/24kbbruHvv36t3qsYChGwhLWH1JEWVG2jsmkVCtS7nCnCUJ1pljcYs7iWupCMUry/xj19b133uDqy8+nelUl5148kRNPPhOXq+X6hVA0TkM4TkmBl5EDC8n2tqXPqVK9i/6XqC4vEkvw86JaorEExU6SqKqs4LqJF/HWv19m/TFjmXr/42y86RbN7m+JPR5TLJ6gINfL2GElFOV5tUe1UinSRKG6tHA0zuyFtVgiFOR6ERFeeu5Jbrr2chpCDZx1/mWceMpZzVZYxxMWwVAMC+hXlM3A0lzyclqv2FZKrSmlRGGM+QswPHl7EWnv6LFKtaohEmfWgmpcxkVetofFixZw5cVn89n/fciWW23D1TfcycjR66+xXySaoD4ax+MyDO2XR1lRDj6tqFaqzVIZ6+kxYBTwLXZnO7DHftJEoTKmOhDm16UBvFkuPG54aPpUptx6Pe6sLCZedyuHHHHcanURIvagfZFYgrxsDxsMLqS4IFt7VSuVBqncUYwDxjqjvSqVUfGExaKVdSyvDlGY6+W3OT9xxUVn8eP337DTP3bnymtvYcDAwavt0xCOE44lKCvysX5pMfk5WdqKSak0SiVR/AgMoJm5rJVKp7oGeya6aDxBnkeYdtt1PHjfFIqKS7ht2gx232u/1RJALG5RF4pSlOdlg/Ii8rK1/kGpTEglUZQBs4wxXwCRxoUisk/GolK9SsISlq+qZ2FFkHyfh19/msmVF5/Ngnlz2f/gI7jgsmsoKSn9fXtLhEB9FLfbsMGQIkoLtR+EUpmUSqKYlOkgVO/1x0x0MbISESZPuoRnnniYIUOH8cDjL7Ld33ZabfuGcJxQNMGQslwGleWtMVqsUir9UumZ/WFHBKJ6l6Yz0X396Xtcffn5VFWu5PiTTuf0cy8mNzfv9+0bB+0rzvOywdAibeaqVAdKpdXTNsBdwJ8AL+AG6kWkMMOxqR4qEkuwYEUdqwJhovU1TLrqEt5+4xU2+NOGa3Scs0Soa4jiMoYNhhbRR4uZlOpwqRQ9TQUOw56DYhxwDLBm43WlUlAdCPPbsgCC8MEbz3PTtZcTjoQ5+4LLOeHkM1frONcQsachHdxHi5mU6kwpdbgTkbnGGLeIJICHjDHfAJdkNjTVkwQaoiypDFJTF8VfuZRrLj+Xzz75H1v+eVuuueFORoxa7/dtI7EE9eEYBbleNhiixUxKdbZUEkWDM1HRt8aYm7CbyepXO7VWIkJdKMaSiiC19VGyDLzy1P3cddtksjxZTLr+Ng4+/NjfO87FE3Y9RI43izFDiykp8Gkxk1JdQCqJ4mjsxHA6cA4wFDgwk0Gp7k1ECDTEWFwRxF8fJdvr5rvPP+CuWycze9YP7Lzrnlxxzc30HzAIsJvH1jXYzV1HDbKH/dYe1Up1Ham0elpojMkBBorIVR0Qk+qmGhPEwpVBgqEo2R4333/xAVNvv5FZP37H0PLh3H73Q+y2574YY5yK6hgCDO2XR/+SXJ2CVKkuKJVWT+OBW7BbPI0wxmwGXK0d7lQjEcFfH2VRRT31oRg+r4vvPl89QVx381TG738IHo8HESEYihGLWwwozWFQWZ4O2qdUF5Zqh7s/A/8FEJFvjTEjMhiT6kb89VEWrKijPhwjx5vFt5+/v1qCuP6Waey938G/t2ayWzLFKSvKZkjffHJ9OtK9Ul1dKv+lMRHxN6lU1AECe7lILMHiiiAVtSFyvFlr3EE0TRChSJyGSJziPC/rDe5DQa62ZFKqu0glUfxkjDkCcBtj1gPOBD7JbFiqq7IsodIfZsGKOkD47vP/MvX2G5j143eUDxvR7B1EyEkQIweVUpjr0ZZMSnUzqSSKM4DLsAcEfAp4G7gmk0Gprqk+HGPe8jrqGmIEa1Zw/cQL+fD9d35PEOP3P4SsLPsj1VjEVJznZfTgUgpzvZ0cvVKqrVJp9dSAnSguy3w4qiuKJyyWrWpgSWU9WUZ49ekHmXLr9RhjuOiK6zjquAl/JAhnbojifI8WMSnVQ6TS6mkccClrToW6SebCUl1FbTDC3GUB4nGLpfNnMemSc5j143fsuPNuXHHNzQwaPBSw7zYisQTF+T7WG1KkCUKpHiSVoqcngAuAHwArs+GoriIcTbBoZR2V/jBuK8K9U27isYfupbRPX+64+2F23XMfjDH29KPROMUFPjYYWky+DrehVI+TSqKoFJFX23JwY8zuwJ3YI84+ICI3NFl/HHAzsNRZNFVEHmjLuVR6iAgVNSHmrwzidsEPX37I1Vecz/KlSzj0qOM598KJFBYVOcN+RyjK87L+kD6aIJTqwVJJFBONMQ8A77H6DHcvtraTMcYNTAN2AZYAXxpjXhWRWU02fUZETl+3sFUmhKMJ5i0PUBuMEgmu4qZrLuOtf7/M6PXH8OQLb7H5uK2xRKgNRnC7dNhvpXqLVBLF8cAYwMMfRU8CtJoosDvpzRWReQDGmKeBfYGmiUJ1MhGhyh9m3vI6EIt3X32KW2+YRCQS4azzL+OEk8/E6/Xq7HJK9VKpJIqtRGSDNhx7MLA46fkSYOtmtjvQGLM9MAc4R0QWN93AGDMBmABQXl7ehlBUSyKxBAtX1FFRG+LHrz7mjpuu4pfZP7HNX7Zn4vW3MXzEKGJxi+q6MEW5OrucUr1RKoniE2PM2GaKjNLhNeApEYkYY04GHgH+3nQjEZkOTAcYN26c9gpPk1XOJEK//PQ9d99+LZ/934cMLR/OrVMfZI+990ewWz253YYNhhRRqsVMSvVKqSSKbbDnopiPXUdhAEmheexS7CHJGw3hj0prsA+yKunpA8BNKcSj2ikWt1iwso5vf/iFh++5mTdefZ6S0j5cOukGDj3yeLxer9Pc1WJwWR6Dy3RUV6V6s1QSxe5tPPaXwHrOAIJLsadTPSJ5A2PMQBFZ7jzdB5jdxnOpFPnro3z5wzwevPcOXnjqIdwuNyefdi4nnnIm+QWFNETiBOvCFOV5GVNeTF62FjMp1dulNB9FWw4sInFjzOnYQ364gRki8pMx5mpgptPk9kxjzD5AHKgGjmvLudTaxRMWvy6sZMpdd/HEg3dRXx9k/4OP5IxzL6Zvv4EEwzFqglH6FDod5nJ0TCallM2IdK8i/3HjxsnMmTM7O4xuI56wWFFdz/0PPMS9U26kYsUydvj7rpx78URGjBpDMBzDGBhYmkf/kmyyvTrst1I9kTHmKxEZ15Z99arQQ0VjCRYsrWL6gzN4+pH7WLp4IRtvugU333EfG2+xDaFonGg8wciBhfQp9GkdhFKqRZooephwNMGPcxYybdpUXnzqYQL+GjbdfCsuuOxqtt1+dxKWhc/jYsTAEgpzvbh0bmql1FpooughGiJxPvrsG6bddSdvvfoc8XiMnXfdi+MnnMbosVuQSAh9irIZUJqjFdRKqXWiiaKbCzZEefnf73LPtDv55MN38fmyOfDQozj2n/+i78BywrEEpQU+hvTN0/oHpVSb6JWjmwrUh3ng4Se4/567+Pmn7yjtU8bp517C4UedQE5BMaFInNzsLMYMLdae1EqpdtFE0c1EYwkefOQpJl83icULfmP4yNFcNfl29jngUIzb7ijndrnYcHgpRXk6q5xSqv00UXQTliW89NpbXHnFpcz64VtGrbcBd977KP/YbS8SFtSFouQYYczQYkoKfNoHQimVNpoouoGPPvmCiy66mE8//oABAwdz3c1T2ffAw8C48DdEyXIbRg8qok9RNm5txaSUSjNNFF3YrNm/cNEll/H6Ky9QVFzChZddwxHH/BNcHgINMVxuw7D+BfQrztZ+EEqpjNFE0QUtXbacy66YxOOPziAry8OE087l+Aln4PLk0hCzyMsxrD+kiKJ87SinlMo8TRRdSG2tn2uvv4G7p00hGolw4KFH889/nUd+cRnGbehfnENZcbb2g1BKdShNFF3AgoWLuOW2O3nkoQcI1gXYZY99mXDGRQwuH0FejpfBfXL07kEp1Wk0UXSiL7/8ihtuvoVXXnoesSz+vuveHHH8vxi7yWYMKM7VuwelVJegiaKDiQivvv4GN910M598/CE5ObkcfPgJHHTkiYwcNZKhffMo1rsHpVQXoomig0QiER6c8Sh33HEbv875mb79BvCvcy9n/IFHUj64P4PL8ijM1TkglFJdjyaKDKuqquKOu6Zx7913s6qqgtEbjOXKyXex8+77MrhfkQ7Sp5Tq8jRRZMjs2T9z0y238fSTjxMOh9j2r3/nsGNPZuvtdmBI3zzKinLwedydHaZSSq2VJoo0EhHee+99brrlFt59+y08Hi+7jT+QQ4+ewIYbbcSQsjxKCrT+QSnVvWiiSINIJMKTTz7Frbfdxk8//kBxSR+OP+VcDjzseEaNHMLA0lytf1BKdVuaKNqhqqqKu++5h2nTplGxciXDR23ABRNvYfz+BzNiUCllRToHtVKq+9OrWBvMnj2b226/nccfe4xwOMy4v+zIhRNvZ5ddd2Ww07xVB+dTSvUUmihSlEgkePPNN7njzjt57z//wev18Y+9DuSI407mL1ttRr8Sbb2klOqZNFGsRW1tLTNmzGDq1GnMnz+Psn4DOOG0CzjsyOPZZMxwSgt8eLK0clop1XNpomjBrFmzmDJlCo899hgNDQ1stNlWXHXzRRxwwP6UDyiiIEcrp5VSvYMmiiSJRILXX3+dKVPu4v3338Pr9bHzHvtx2NH/ZOcdtqFPYQ7ZXu37oJTqXTRRYLdeevjhh5k6dRoLFy6gX/9BTDjzEg494hg2HzuCwjyvVk4rpXqtXpsoLMvigw8+YPr0+3n55ZeIRqNsuuU2nHTWZRx84AEM7legldNKKUUvTBTLly/noYce4oEHHmT+/HkUFhazz8FHc+AhR7PT3/5Mn0LtOa2UUsl6RaJIJBK89dZb3H///bz++uskEgk232pbJp1yHgcdeCDDB5WSm52lldNKKdWMHp0oFi5cyIwZM5gxYwZLliyhtE8Zhx93KgceehTbb7MZhbkevXtQSqm16LGJYvr06ZxyyikAbPvXnTjjwqvYb999KR9QrC2XlFJqHfTYRLHDDjtwzvkXs8v4Q9h2i7E6KJ9SSrWREZHMHdyY3YE7ATfwgIjc0GS9D3gU2BJYBRwqIgtaO+a4ceNk5syZmQlYKaV6KGPMVyIyri37ZqyA3hjjBqYBewBjgcONMWObbHYiUCMio4HbgRszFY9SSqm2yWRN7p+BuSIyT0SiwNPAvk222Rd4xHn8PLCz0fIhpZTqUjKZKAYDi5OeRytFCAAABQdJREFUL3GWNbuNiMQBP9Cn6YGMMROMMTONMTMrKyszFK5SSqnmdIu2oSIyXUTGici4vn37dnY4SinVq2QyUSz9//bu7sWKOo7j+PuDrj2gZGCQkLYFIT3Qg4EoRXgTiBd2kYEE1Qbd9EB1GV0Y9Qd0UVESKfQgYVTEJkoICVHQloiWtm3oVYFgD6RGD2B8u5hZWk7nzM45O+fM/KbPCw78zs7vLN8v3935npkz5zfAqjnPr8h/1nWOpMXAJWQfapuZWUMMs1F8CVwj6SpJS4BtwGTHnEnggXy8Ffg4hnkZlpmZ9W1o36OIiPOSHgM+Irs8dldEHJf0HHAoIiaBncCbkk4Av5A1EzMza5ChfuEuIvYB+zp+tn3O+E/gnmHGYGZmCzPUL9wNg6RzwEzdcQzRCuCnuoMYIueXrjbnBu3Pb01ELBvkhSku4TEz6LcLUyDpkPNLV5vza3Nu8P/Ib9DXJnF5rJmZ1ceNwszMCqXYKF6tO4Ahc35pa3N+bc4NnF9PyX2YbWZmo5XiEYWZmY2QG4WZmRVqbKOQtEnSjKQTkp7qsv0CSXvy7VOSxkcf5eBK5Dch6UdJR/LHQ3XEOQhJuySdlnSsx3ZJeiHP/StJa0cd40KUyG+jpDNzare927wmkrRK0kFJ30g6LumJLnOSrV/J/FKu34WSvpB0NM/v2S5z+t93RkTjHmRLfpwErgaWAEeB6zrmPALsyMfbgD11x11xfhPAS3XHOmB+dwBrgWM9tm8G9gMC1gNTdcdccX4bgb11xzlgbiuBtfl4GfBdl7/NZOtXMr+U6ydgaT4eA6aA9R1z+t53NvWIou03PSqTX7Ii4hOytbt6uQt4IzKfA8slrRxNdAtXIr9kRcSpiDicj88B0/z3PjLJ1q9kfsnKa/Jb/nQsf3ResdT3vrOpjaKymx41VJn8AO7OD+3flbSqy/ZUlc0/ZRvyw//9kq6vO5hB5KckbiF7VzpXK+pXkB8kXD9JiyQdAU4DByKiZ/3K7jub2igMPgTGI+JG4AD/vgOw5jsMXBkRNwEvAh/UHE/fJC0F3gOejIizdcdTtXnyS7p+EfF3RNxMdg+gdZJuWOjvbGqjaPtNj+bNLyJ+joi/8qevAbeOKLZRKFPfZEXE2dnD/8hWUB6TtKLmsEqTNEa2E90dEe93mZJ0/ebLL/X6zYqIX4GDwKaOTX3vO5vaKNp+06N58+s457uF7FxqW0wC9+dXz6wHzkTEqbqDqoqky2fP+UpaR/Z/lsSbmDzuncB0RDzfY1qy9SuTX+L1u0zS8nx8EXAn8G3HtL73nY1cPTZaftOjkvk9LmkLcJ4sv4naAu6TpLfJrhxZIekH4BmyD9WIiB1k9yjZDJwAfgcerCfSwZTIbyvwsKTzwB/AtoTexNwG3Ad8nZ/nBngaWA2tqF+Z/FKu30rgdUmLyBrcOxGxd6H7Ti/hYWZmhZp66snMzBrCjcLMzAq5UZiZWSE3CjMzK+RGYWZmhdwozBZI0nivlWTN2sCNwszMCrlRmFVjsaTdkqbzRRwvrjsgs6q4UZhVYw3wckRcC5wlW/PfrBXcKMyq8X1EfJaP3wJurzMYsyq5UZhVo3MtHK+NY63hRmFWjdWSNuTje4FP6wzGrEpuFGbVmAEelTQNXAq8UnM8ZpXx6rFmZlbIRxRmZlbIjcLMzAq5UZiZWSE3CjMzK+RGYWZmhdwozMyskBuFmZkV+gfaMvyVROiksgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0.1, 3.0, 30)\n",
    "y = []\n",
    "e1 = []\n",
    "e2 = []\n",
    "for b in x:\n",
    "    z = []\n",
    "    for _ in range(25):\n",
    "        z.append(sim_run(0, b))\n",
    "    z = np.array(z)\n",
    "    y.append(z.mean())\n",
    "    e1.append(np.percentile(z, 5))\n",
    "    e2.append(np.percentile(z, 95))\n",
    "\n",
    "plt.plot(x, y, 'k-')\n",
    "plt.xlim(0, 3.0)\n",
    "plt.title('arrival=exp(1.2), service=truncnorm(0, b)')\n",
    "plt.ylabel('mean customer wait time')\n",
    "plt.xlabel('b')\n",
    "plt.fill_between(x, e1, e2, color='#386cb0', alpha=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above plot shows that the mean wait time increases as we increase 'b', which increases both the mean and the standard deviation of the service time. The confidence interval also seems to be getting larger as we increase 'b'. \n",
    "\n",
    "As 'b' increases, the service rate decreases and the difference between the arrival rate and service rate decreases as a result. It takes longer for the simulation to reach steady state. We fixed each simulation run to last for 1000 simulated seconds, which may not be sufficient for a larger 'b'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}