From 35c1cee4b726f1d9ec06cf5642610b365a1e5b27 Mon Sep 17 00:00:00 2001 From: TomMonks Date: Mon, 27 Apr 2026 16:18:26 +0100 Subject: [PATCH 1/3] feat(nspp): nb 16 data for NSPP --- content/data/nspp_example1.csv | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 content/data/nspp_example1.csv diff --git a/content/data/nspp_example1.csv b/content/data/nspp_example1.csv new file mode 100644 index 0000000..69a19be --- /dev/null +++ b/content/data/nspp_example1.csv @@ -0,0 +1,10 @@ +ļ»æt,mean_iat +0,15 +60,12 +120,7 +180,5 +240,8 +300,10 +360,15 +420,20 +480,20 \ No newline at end of file From dc4c0b22dee91e8a04080a5155b3023b0d9ef54b Mon Sep 17 00:00:00 2001 From: TomMonks Date: Mon, 27 Apr 2026 16:18:45 +0100 Subject: [PATCH 2/3] feat(nspp): nb 16 code for NSPP --- content/16_time_dependent_arrivals.ipynb | 900 +++++++++++++++++++++++ 1 file changed, 900 insertions(+) create mode 100644 content/16_time_dependent_arrivals.ipynb diff --git a/content/16_time_dependent_arrivals.ipynb b/content/16_time_dependent_arrivals.ipynb new file mode 100644 index 0000000..951bcba --- /dev/null +++ b/content/16_time_dependent_arrivals.ipynb @@ -0,0 +1,900 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Modelling a non-stationary poisson process\n", + "\n", + "> A non-stationary Poisson process (NSPP) is an arrival process where inter-arrival times are Exponentially distributed with a mean that varies by time.\n", + "\n", + "One of the limitations of queuing theory is the difficulty of modelling time-dependent arrivals. Computer simulation offers a number of ways of modelling non-stationary arrivals. \n", + "\n", + "šŸŽ“ **In this notebook you will learn:**\n", + "\n", + "* āœ… How to implement the thinning algorithm to model a non-stationary poisson process (NSPP)\n", + "* āœ… How to check that that thinning is working correctly.\n", + "* šŸŽ **BONUS**: How to improving the efficiency of your thinning algorithm code in Python. \n", + "\n", + "---" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Imports" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'4.1.1'" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.ticker as ticker\n", + "import itertools\n", + "import simpy\n", + "\n", + "# please use simpy version 4\n", + "simpy.__version__" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "----\n", + "## An example NSPP\n", + "\n", + "The table below breaks an arrival process down into 60 minutes intervals.\n", + "\n", + "\n", + "| t(min) | Mean time between arrivals (min) | Arrival Rate $\\lambda(t)$ (arrivals/min) |\n", + "|:------:|:--------------------------------:|:--------------------------------------:|\n", + "| 0 | 15 | 1/15 |\n", + "| 60 | 12 | 1/12 |\n", + "| 120 | 7 | 1/7 |\n", + "| 180 | 5 | 1/5 |\n", + "| 240 | 8 | 1/8 |\n", + "| 300 | 10 | 1/10 |\n", + "| 360 | 15 | 1/15 |\n", + "| 420 | 20 | 1/20 |\n", + "| 480 | 20 | 1/20 |\n", + "\n", + "> **Interpretation**: In the table above the fastest arrival rate is 1/5 customers per minute or 5 minutes between customer arrivals." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Thinning\n", + "\n", + "Thinning is a acceptance-rejection sampling method and is used to generate inter-arrival times from a NSPP. \n", + "\n", + "> A NSPP has arrival rate $\\lambda(t)$ where $0 \\leq t \\leq T$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**The thinning algorithm**\n", + "\n", + "A NSPP has arrival rate $\\lambda(t)$ where $0 \\leq t \\leq T$\n", + "\n", + "Here $i$ is the arrival number and $\\mathcal{T_i}$ is its arrival time.\n", + "\n", + "1. Let $\\lambda^* = \\max_{0 \\leq t \\leq T}\\lambda(t)$ be the maximum of the arrival rate function and set $t = 0$ and $i=1$\n", + "\n", + "2. Generate $e$ from the exponential distribution with rate $\\lambda^*$ and let $t = t + e$ (this is the time of the next entity will arrive)\n", + "\n", + "3. Generate $u$ from the $U(0,1)$ distribution. If $u \\leq \\dfrac{\\lambda(t)}{\\lambda^*}$ then $\\mathcal{T_i} =t$ and $i = i + 1$\n", + "\n", + "4. Go to Step 2." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exercise 1: simulation **without thinning**\n", + "\n", + "**Task:**\n", + "* Build a simple `simpy` model that simulates time-dependent arrivals\n", + "* For this exercise please **IGNORE** the need for a thinning process.\n", + "\n", + "**Optional task:**\n", + "* It is useful to set the sampling of arrivals using a random seed. This will allow you to compare the number of arrivals before and after adding thinning. **Remember that an issue with DES without thinning occurs when moving from a period $t$ with a low arrival rate to $t+1$ that has a high one.**\n", + "\n", + "**Hints:**\n", + "* Build your model up gradually. \n", + "* Start by building a model that simulates exponential arrivals using a single mean inter-arrival time then add in logic to change which mean you use depending on the simulation time.\n", + "* The logic to decide the time period is equivalent to asking yourself \"given `env.now()` and that arrival rates are split into 60 minute chunks which row of my dataframe should I select\".\n", + "* To simplify the task you set the run length of the simulation to no more than 540 minutes. For an extra challenge think about how you would run the model for longer than 480 minutes and loop back round to the first period (the code to do this is surprising simple).\n", + "\n", + "The data are stored in a file `data/nspp_example1.csv`. " + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# your code here ..." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
tmean_iatarrival_rate
00150.066667
160120.083333
212070.142857
318050.200000
424080.125000
5300100.100000
6360150.066667
7420200.050000
8480200.050000
\n", + "
" + ], + "text/plain": [ + " t mean_iat arrival_rate\n", + "0 0 15 0.066667\n", + "1 60 12 0.083333\n", + "2 120 7 0.142857\n", + "3 180 5 0.200000\n", + "4 240 8 0.125000\n", + "5 300 10 0.100000\n", + "6 360 15 0.066667\n", + "7 420 20 0.050000\n", + "8 480 20 0.050000" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# example answer\n", + "\n", + "# read in the data and calculate lambda\n", + "arrivals = pd.read_csv('data/nspp_example1.csv')\n", + "arrivals['arrival_rate'] = 1 / arrivals['mean_iat']\n", + "arrivals" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/IAAAEmCAYAAADBUqsoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABZ70lEQVR4nO3dB3gVVfrH8V96CISQQkgCCTWUQEKRIkVElKKAIsW+Ytld/atrWwWxggUU67Lquq4uthVUqoqKsChFkCIdBIEACRAINYGE9Pyfc2KyCQIGSHJLvp/nGZKZWzJ3Dvfeeec95z0ehYWFhQIAAAAAAC7B09E7AAAAAAAAyo9AHgAAAAAAF0IgDwAAAACACyGQBwAAAADAhRDIAwAAAADgQgjkAQAAAABwIQTyAAAAAAC4EAJ5AAAAAABciLejd8AZFRQUaO/evQoMDJSHh4ejdwcAAAAA4OYKCwt17NgxRUVFydPzzDl3AvlTMEF8dHR0ZbUPAAAAAACnlJycrAYNGuhMCORPwWTiiw9g7dq1z3gAAQAAAAA4X+np6TahXByPngmB/CkUd6c3QTyBPAAAAACgqpRneDfF7gAAAAAAcCEE8gAAAAAAuBACeQAAAAAAXIhDA/nx48erU6dOdjB/eHi4Bg8erC1btvymBP+YMWNsCf4aNWqoV69e2rhx4+8+97Rp0xQXFyc/Pz/7c8aMGZX4SgAAAAAAqAaB/IIFC3T33Xfrxx9/1Ny5c5WXl6e+ffsqIyOj5D4TJkzQK6+8otdff10rVqxQRESE+vTpY+fXO52lS5fq2muv1R/+8AetXbvW/rzmmmu0bNmyKnplAAAAAABUDo9Ck/J2EgcOHLCZeRPg9+zZ02bjTSb+/vvv16hRo+x9srOzVa9ePb3wwgu64447Tvk8Jog3pfu//vrrkm39+/dXcHCwJk+e/Lv7YR4bFBSktLQ0qtYDAAAAACrd2cShTjX9nNlhIyQkxP7csWOH9u3bZ7P0xUxX+YsvvlhLliw5bSBvMvIPPPBAmW39+vXTa6+9dsr7m4sDZil9AAEAAADgdHLzC/T2wkR9tT5F+QVOkxvFGcy4q7tq+HrJHThNIG+y7w8++KB69OihNm3a2G0miDdMBr40s75r167TPpd53KkeU/x8pxqrP3bs2Ap4FQAAAADc3aa96Xp46lpt3EsC0JUUOE9ndPcJ5O+55x6tW7dOixcv/s1tHh4evwn6T952Po8ZPXq0vYhQOiMfHR19lq8AAAAAgDvLySvQG99ts0teQaHqBPjo4X4t1DCkpqN3DeXg5+0+k7Y5RSD/l7/8RZ9//rkWLlyoBg0alGw3he0Mk0mPjIws2Z6amvqbjHtp5nEnZ9/P9BjTXd8sAAAAAHAqG/ak6aHP1mrzvqKi2/1a19Mzg9soPNCfA4Yq59BLEiZLbjLx06dP1/z589W4ceMyt5t1E5SbivbFcnJybDG8bt26nfZ5u3btWuYxxrfffnvGxwAAAADAybLz8vXyt1t01Rs/2CA+OMBHf7++vd666QKCeFTPjLyZeu7jjz/WrFmz7FzyxVl0U6nPzBlvusKbivXjxo1TbGysXczvAQEBuuGGG0qe5+abb1b9+vXtWHfjvvvus1XvTWX7q666yj7/vHnzTtltHwAAAABOZd3uo3r4s3Xasr8oCz8gPlJjr2qtsFr05kU1DuT/8Y9/2J+9evUqs33SpEm65ZZb7O8jR47UiRMndNddd+nIkSPq0qWLza6bwL9YUlKSPD3/17nAZN6nTJmixx9/XE888YSaNm2qTz75xD4WAAAAAM4kKzdfE/+7Vf9cmGgr0ofW9NXTV7XRgIT/DfcFHMmp5pF3FswjDwAAAFRPq5OO6OGp67Qt9bhdH9Q2SmMGxSmULDwqmcvOIw8AAAAAjsrCvzrvF/1rYaLMtPCm+/yzg9uof5uiAtyAMyGQBwAAAFCt/bTLZOHXKvFAhl0f3C5KTw1qreCavo7eNeCUCOQBAAAAVEsncooq0r/7ww6ZAcd1A/007up49Yk7/VTXgDMgkAcAAABQ7azYeVgjp67TjoNFWfghHerryYFxqhNAFh7Oj0AeAAAAQLWRmZOnF+ds0XtLdtosfL3afho/JF69W5KFh+sgkAcAAABQLfyYeEijpq3TrkOZdv2ajg302IA4BdXwcfSuAWeFQB4AAACAW8vIztOEbzbr/aW77HpkkL/NwvdqEe7oXQPOCYE8AAAAALe1ZNtBjZq+TsmHT9j16zvH6NErWirQnyw8XBeBPAAAAAC3czw7T89//bM++jHJrtevU0PPD43XRbF1Hb1rwHkjkAcAAADgVhZvPWjHwu85WpSFv+nCGD1yeSvV8iP8gXvgfzIAAAAAt3AsK1fjvvpZk5cn2/XokBp6YUiCujULc/SuARWKQB4AAACAy1vwywGNnrZOe9Oy7PqIrg01sn9L1SQLDzdEIA8AAADAZaWdyNVzszfp05W77XrD0AC9MDRBFzYJdfSuAZWGQB4AAACAS5q/eb8enb5B+9Kz5OEh3dKtkR7u10IBvoQ5cG/8DwcAAADgUtIyc/X0l5s0bVVRFr5xWE1NGJagTo1CHL1rQJUgkAcAAADgMuZu2q/HZqxX6rFsm4X/Y4/GerBPC9Xw9XL0rgFVhkAeAAAAgNM7kpGjsV9s1Mw1e+16k7o19eKwtrqgYbCjdw2ocgTyAAAAAJzaNxv26fGZG3TweLY8PaQ/9WyiBy5rLn8fsvCongjkAQAAADilwxk5eurzjfpibVEWPja8lh0L3z6GLDyqNwJ5AAAAAE7nq/UpemLmBh3KyJGXp4fu6NlE914aSxYeIJAHAAAA4ExM9/mnZm3U7PUpdr1FvUC9ODxBCQ3qOHrXAKdBRh4AAACAwxUWFurLdSm2K/3hX7Pwd/dqqrt7N5OfN2PhgdII5AEAAAA4VOqxLNuNfs7G/Xa9VWRtvTgsQW3qB9EywCkQyAMAAABwWBZ+1pq9GvPFRh3NzJW3p4fu6d1Md/VqJl9vT1oFOA2HvjsWLlyoQYMGKSoqSh4eHpo5c2aZ2822Uy0vvvjiaZ/zvffeO+VjsrKyquAVAQAAACiP1PQs/emDn3T/J2tsEN86qrY+v6eH7r+sOUE84MwZ+YyMDLVt21a33nqrhg4d+pvbU1KKClwU+/rrr3X77bef8r6l1a5dW1u2bCmzzd/fv4L2GgAAAMD5ZOGnr9qjsV9sVHpWnny8PHRv71jd2aupfLzIwgNOH8hffvnldjmdiIiIMuuzZs3SJZdcoiZNmpzxeU0G/uTHAgAAAHCsfWlZenTGes3fnGrX4+sH2Yr0LSNq0zSAO46R379/v2bPnq3333//d+97/PhxNWzYUPn5+WrXrp2eeeYZtW/f/rT3z87Otkux9PT0CttvAAAAoLozWfjPftqtZ77cpGNZefL18tT9fWL154uayJssPOC+gbwJ4AMDAzVkyJAz3q9ly5Z2nHx8fLwNyP/2t7+pe/fuWrt2rWJjY0/5mPHjx2vs2LGVtOcAAABA9bX36Ak9Mn29Fv5ywK63ja6jl4YlKLZeoKN3DXBZHoXm8pgTMN3hZ8yYocGDB582QO/Tp4/+/ve/n9XzFhQUqEOHDurZs6cmTpxY7ox8dHS00tLS7Hh7AAAAAGfHhBlTViTrudk/63h2ni1g99c+zXV7j8Zk4YFTMHFoUFBQueJQl8jIL1q0yBav++STT876sZ6enurUqZO2bt162vv4+fnZBQAAAMD5230kU6Onr9eirQfteoeYOpowrK2ahdfi8AIVwCUC+XfffVcXXHCBrXB/LlcC16xZY7vaAwAAAKg8BQWF+nh5ksZ/9bMycvLl5+2ph/u10K3dG8vL04NDD7hDIG+K0m3btq1kfceOHTboDgkJUUxMTEn3gs8++0wvv/zyKZ/j5ptvVv369e04d8OMdb/wwgvteHjzWNOd3jznG2+8UUWvCgAAAKh+kg9natS0dVqy/ZBd79QoWC8MTVCTumThAbcK5FeuXGmnkyv24IMP2p8jRoywBeuMKVOm2Kz69ddff8rnSEpKst3nix09elR//vOftW/fPju+wFSrX7hwoTp37lzprwcAAACojln4j5bt0vNfb1ZmTr78fTw1qn9LjejaSJ5k4QH3LnbnqkUGAAAAgOpq16EMjZy6Tst2HLbrnRuHaMLQBDUKq+noXQNcjtsVuwMAAADgXFn495fu1IRvtuhEbr4CfL30yOUtdVOXhmThgSpAIA8AAACg3HYcNFn4tVqx84hd79okVBOGJSg6JICjCFQRAnkAAAAAvyu/oFCTftihF+dsUXZegWr6emn0Fa10Q+cYsvBAFSOQBwAAAHBG2w8c18OfrdWqpKN2vUezMD0/NF4NgsnCA45AIA8AAADgtFn4dxYl6uW5vygnr0C1/Lz1+IBWurZTtDw8mBcecBQCeQAAAAC/sXX/MT08dZ3WJBdl4Xs2r6vxQ+JVv04NjhbgYATyAAAAAErk5Rfo7UWJem3uVuXkFyjQ31tPDIzT8AsakIUHnASBPAAAAABryz6ThV+rdbvT7HrvluEad3W8IoL8OUKAEyGQBwAAAKq53PwC/XPBdv3tv1uVm1+o2v7eempQaw3pUJ8sPOAugfzOnTu1aNEi+zMzM1N169ZV+/bt1bVrV/n7c7UOAAAAcBU/p6Troc/WauPedLt+Wat6Gnd1G4XX5rwecItA/uOPP9bEiRO1fPlyhYeHq379+qpRo4YOHz6s7du32yD+xhtv1KhRo9SwYcPK22sAAAAA58VUoX/z+216ff425RUUqk6Aj8Ze2VpXto0iCw+4SyDfoUMHeXp66pZbbtGnn36qmJiYMrdnZ2dr6dKlmjJlijp27Kg333xTw4cPr4x9BgAAAHAeNuxJsxXpTTbe6Ne6np4Z3EbhgWThAVfgUVhYWFieO86ePVsDBgwo15MePHhQO3bsUKdOneSK0tPTFRQUpLS0NNWuXdvRuwMAAABUWBb+9flb9eb3220WPqSmr83CD0yIJAsPuFAcWu6MfHmDeCMsLMwuAAAAAJzD+t0mC79Wm/cds+sD4iM19qrWCqvl5+hdA1AVxe5WrVolHx8fxcfH2/VZs2Zp0qRJiouL05gxY+Tr63suTwsAAACggmXn5Wvif7fqrQWJyi8oVGhNX9uN/or4SI414KI8z+VBd9xxh3755Rf7e2Jioq677joFBATos88+08iRIyt6HwEAAACcgzXJRzVw4mK98d12G8QPahuluQ9eTBAPVMeMvAni27VrZ383wXvPnj1tRfsffvjBBvWvvfZaRe8nAAAAgHLKys3Xq/N+0b8WJqqgULb7/LOD26h/mwiOIVBdA3lTH6+goMD+Pm/ePA0cOND+Hh0dbQvdAQAAAHCMn3Yd0cipa7X9QIZdv7p9fT05ME7BNRn+ClTrQN5ML/fss8/qsssu04IFC/SPf/zDbjeV6uvVq1fR+wgAAACgHFn4l7/doncW75CZlyo80E/PXR2vPnGcnwPu5pwCedN1/sYbb9TMmTP12GOPqVmzZnb71KlT1a1bt4reRwAAAABnsHLnYY2cuk6JB4uy8EM7NLBZ+KAAH44bUJ3nkS+PrKwseXl52Yr2rox55AEAAOAKTuTk68U5WzRpSVEWPqK2v8YNaaPeLcnCA66mUuaRP5WcnBylpqaWjJcvFhMTcz5PCwAAAOB3LEs8pJHT1mnXoUy7fk3HBnpsQJyCarh2Ug1AJVatv/3227VkyZIy201y38PDQ/n5+efytAAAAAB+R0Z2niZ8s1nvL91l1yOD/DV+SLx6tQjn2AHVxDkF8rfeequ8vb315ZdfKjIy0gbvAAAAACrXku0HNWraOiUfPmHXr+8co0evaKlAf7LwQHXieS4PWrNmjf75z3/q8ssvt/PJt23btsxSXgsXLtSgQYMUFRVlLwaY4nml3XLLLXZ76eXCCy/83eedNm2a4uLi5OfnZ3/OmDHjXF4mAAAA4BSOZ+fp8ZnrdcO/ltkgvn6dGvrw9s42E08QD1Q/5xTIm+C4IuaLz8jIsIH/66+/ftr79O/fXykpKSXLV199dcbnXLp0qa699lr94Q9/0Nq1a+3Pa665RsuWLTvv/QUAAACq2uKtB9Xv1YX66Mcku37ThTGa80BPXRRbl8YAqqlzqlo/f/58Pf744xo3bpzi4+N/U6X+9yrsnXJHPDxs5nzw4MFlMvJHjx79Tab+TEwQb6r9ff3112UuBgQHB2vy5MluVbW+oKBQ+YWF8vE6p+sxAAAAcGLHsnI17qvNmry8KICPDqmhF4YkqFuzMEfvGgBXrFp/2WWX2Z+XXnpppRe7+/777xUeHq46dero4osv1nPPPWfXz5SRf+CBB8ps69evn1577bXTPiY7O9supQ+gK1iw9YCdL3T4BQ10XacYxYQGOHqXAAAAUAEW/HJAo6et0960LLs+omtDjezfUjX9zmvSKQBu4pw+Cb777jtVBTMGf/jw4WrYsKF27NihJ554Qr1799ZPP/1kx7+fyr59+1SvXtl5M8262X4648eP19ixY+Vqvli7VweOZevN77fb5aLYMFvw5LJW9eTrTZYeAADA1aRn5eq5L3/WJyuT7XrD0AC9MDRBFzYJdfSuAXD1QN5kxquC6SZfrE2bNurYsaMN6mfPnq0hQ4ac9nEnV9Ev7ilwOqNHj9aDDz5YJiMfHR0tZ2c+1PvG1dPHy5O1aOsBLdp60C6hNX01rGNRlr5xWE1H7yYAAADK4bvNqRo9fb32pWfJnLre0q2RHu7XQgG+ZOEBlFXuT4V169bZYNrT09P+fiYJCQmqDGaqOxPIb9269bT3iYiI+E32PTU19TdZ+tJMdv90GX5nZsbG928TaZfkw5n6dGWyPlmRrNRj2frngkS7dGsaqus6x6hf63ry8/Zy9C4DAADgJGmZuXr6y02atmq3XTeJmAnDEtSpUQjHCsD5BfJmmjkTIJvx6eZ3k+E+VZ28ih4jX9qhQ4eUnJxsA/rT6dq1q+bOnVtmnPy3336rbt26yZ1FhwTor31b6L5LYzV/c6qmrEjWd1tStWT7IbsEB/homBlL3zlGTevWcvTuAgAAQNK8Tfv16Iz1NhFjsvB/7NFYD/ZpoRq+JGAAVEAgb8ao161bt+T3inD8+HFt27atzN8wc9SHhITYZcyYMRo6dKgN3Hfu3KlHH31UYWFhuvrqq0sec/PNN6t+/fp2nLtx3333qWfPnnrhhRd01VVXadasWZo3b54WL16s6sDby1N9W0fYZc/RE/p0RbLN1KekZelfi3bYpXPjEN3QOUb920TI34cvCQAAgKp2NDNHY7/YpBmr99j1JnVr6sVhbXVBw2AaA0DlTD9XkRXpL7nkkt9sHzFihP7xj3/YqehWr15tp6Azwby57zPPPFNm/HqvXr3UqFEjvffeeyXbpk6daqfHS0xMVNOmTW2l+zONqXfV6efKKy+/wFY+nbw8WfM371fBry0eVMNHQzrUtwXymtcLdPRuAgAAVAtzNu7TYzM26ODxbHl6SH/q2UQPXNacBAtQzaWfRRx6zoH8nj179MMPP9jx5wUFBWVuu/fee+XK3C2QLy0l7YQ+W7nbjqU3GftiHRsG24D+ivhIunIBAABUgsMZORrz+UZ9vnavXY8Nr2XHwrePIQsPQJUfyE+aNEl33nmnfH19FRoaWqYivPndZMJdmTsH8sXyCwptpfvJy5M07+dUu24E+ntrSPv6ur5LjFpGuOdrBwAAqGpfrU/Rk7NMFj5HXp4euqNnE917aSxZeABVF8ibru0mkDfTtpkq9u6mOgTypaWmZ+mzn3ZryookJR/+X5a+fUwdXd8pRgPbRjLtCQAAwDkw3eefmrVRs9en2PUW9QL14vAEJTSow/EEULWBvMnCL1++3I4/d0fVLZAvVlBQqB+2H7RZ+m837ldecZbez1tXtY+yXe9bRwU5ejcBAACcnjnF/nJdip76fKPtUm+y8Hf3aqq7ezdjSmAAjgnkR44caavKP/LII3JH1TWQL+3AsWw7l+mU5UnaeSizZHtCgyAb0A9qG6VafuWe9AAAAKBanUc9MXODvtm4z663iqytF4clqE19EiIAHBjIm3niBw4cqBMnTig+Pl4+Pj5lbn/llVfkygjky2bpf0w8pMkrkvXNhhTl5hf9d6np66Ur2xVl6ePrB5WpkwAAAFAdmdNqU8jOZOGPZubK29ND9/Ruprt6NZOvt/sNRwXguDj0nFKq48aN05w5c9SiRQu7fnKxO7gPT08PdWsWZpdDx+M0fdUe2/U+8WCGnc7OLK2jatuA/qp2UQr0L3tRBwAAoDowNYcem7lBczftt+vm/MjMCx8XVT17dwKoXOeUkQ8ODtarr76qW265Re6IjPyZmf8yy3cctgH9Vxv2KSevaPrBGj5eGtQ20gb17aLrcFEHAABUi/OiGav3aOwXm5R2Ilc+Xh66t3es7uzVVD5eZOEBOFHX+oiICC1atEixsbFyRwTy5XckI0fTV++xY+m3ph4v2d4yItAG9IPb11dQDbL0AADA/exLy9KjM9Zr/uZUu26GG5qK9EzhC8ApA/nx48crJSVFEydOlDsikD975r/RT7uO6OPlSZq9LkXZv2bp/X08NSDejKWP1gUNg8nSAwAAtzjvMVP3PvPlJh3LypOvl6fu7xOrP1/URN5k4QE4ayB/9dVXa/78+XYautatW/+m2N306dPlygjkz09aZq5mrikaS79537GS7bHhtWyWfkiH+qoT4Hve7QQAAFDV9h49odHT12vBLwfsetvoOnppWIJi6wXSGACcO5C/9dZbz3j7pEmT5MoI5CuG+a+1Ovmo7Xb/xdoUncjNt9tN1dYr2kTYoL5z4xCy9AAAwCXOaz5ZkaznZv+sY9l59nzmr32a6/YejcnCA3CNQN7dEchXwjHNytWsNXs1eVmSNqWkl2xvUrembrBZ+gYKqUmWHgAAOJ89R0/okWnrtGjrQbveIaaOJgxrq2bhtRy9awDcCIF8FR5AnB1z3Wj9njTb7d4E9pk5v2bpvTzVz2bpo9W1SShZegAA4BTnLab+z7jZPysjJ19+3p56uF8L3dq9sbw8mXIZgAsE8v3799eTTz6pbt26nfF+x44d05tvvqlatWrp7rvvlisikK8ax7Pz9LnJ0i9PssF9scZhNXVdp2gNvaCBwmr5VdHeAAAA/E/y4UyNmrZOS7YfsusdGwZrwrAENalLFh6ACwXy7777rp566ikFBgbqyiuvVMeOHRUVFSV/f38dOXJEmzZt0uLFi/XVV19p4MCBevHFFxUdHS1XRCBf9TaUytKbAN8w87D2jYvQdZ2j1b1pmDy58g0AACpZQUGhPlq2S89/vdn2HDQz8Izs11IjujUiCw/ANbvW5+TkaOrUqfrkk0/sPPJHjx4tehIPD8XFxalfv37605/+pBYtWsiVEcg7TkZ2np2+znRjW5Nc9P/LiA6poes6xWh4xwYKD/R34B4CAAB3tetQhkZOXadlOw7bdVOUd8LQBDUKq+noXQNQDaRXVbE78wdOnDhhp6E7eQo6V0Yg7xx+Tkm3Fe+nr95j52g1vD09dGmrcFvx/qLYulwZBwAAFZKFf3/pTk34ZoudZSfA10uj+rfUHy5sSI9AAFWGYndVeABR+U7k5Gv2+hTb9f6nXUdKttevY7L00RreMVoRQWTpAQDA2dtxMEOjpq7T8p1FWXhTdPeFoQmKCQ3gcAKoUgTyVXgAUbV+2X/MBvTTV+1R2olcu81Ujb2kRbhu6BKti5uHk6UHAAC/K7+gUJN+2KGXvt2irNwC1fT10ugrWtlpcanLA8ARCOSr8ADCMbJy8/X1hhRNXpZccgXdiAry1zWdonVNx2hF1alB8wAAgN/YfuC4Hv5srVYlFdXj6dEsTM8PjVeDYLLwAByHQL4KDyAcb1vqMU1Znqxpq3brSGZRlt4UuO/Vomgs/SUt6srby9PRuwkAAJwgC//OokS9MvcXZecVqJaftx4b0MoO1TPFmwHALQP5hQsX2nnkvb295c4I5F1Tdl6+5mzcr8nLkrQ0sWjOV6NebT+boTdLdAhX2gEAqK4X/h/6bF3JrDg9m9fV+CHxtuYOALh1IO/l5aWUlBSFh4fLnRHIu77EA8f1yYpkTf1ptw5l5Nht5kJ7z9i6ur5ztC5tVU8+ZOkBAHB7efkFentRol6bt1U5eQUK9PfWEwPi7JS2ZOEBVItA3tPTU/v27SOQh8swX9hzN+23BfIWbztYsj2slsnSN7Bz01OVFgAA97Rl3zE9PHWt1u1Os+tmuN24IfGKDCILD8C1A/mzHjhckVcuTVf9QYMGKSoqyj7vzJkzS27Lzc3VqFGjFB8fr5o1a9r73Hzzzdq7d+8Zn/O9996zz3XykpWVVWH7Ddfh6+2pAQmR+uiPXbTg4V66q1dTG8QfPJ6tN7/frp4vfqeb3lmm2etSbNAPAABcX25+gV6fv1UD/77IBvG1/b318vC2+vctnQjiAbiFsx7s/sQTTygg4MzjjF955ZVyPVdGRobatm2rW2+9VUOHDi1zW2ZmplatWmX/nrnPkSNHdP/99+vKK6/UypUrz/i85urFli1bymzz92ee8equYWhNjezfUg/0aa7//rxfHy9P1qKtB2ym3iyhNX017IIGuq5zjBqH1XT07gIAgHPwc0q6HvpsrTbuTbfrl7UK13NXx6tebc4FAVTjQH79+vXy9fWtkIz95ZdfbpdTMV0K5s6dW2bb3//+d3Xu3FlJSUmKiYk54z5ERESUez9QvZix8f3bRNol+XCmPl2ZbMfTpx7L1j8XJtqla5NQXd8lRv1a15Oft5ejdxkAAPwO07Puze+36Y3vtik3v1B1Anw0ZlBrXdWuqOcnAFTrQH7GjBkOGyNvxgqYD+I6deqc8X7Hjx9Xw4YNlZ+fr3bt2umZZ55R+/btT3v/7Oxsu5Qem4DqwVSx/2vfFrrv0ljN35yqKSuS9f2WVFv13izBAT4a2qEoS98svJajdxcAAJzCxr1ptiK9ycYbfePq6dmr2yg8kCw8APd0VoG8I69mmjHujzzyiG644YYzDvxv2bKlHSdvxtabgPxvf/ubunfvrrVr1yo2NvaUjxk/frzGjh1biXsPZ2fmme/bOsIue46e0Kcrkm2mPiUtS+8s3mGXzo1DbMX7y9tEyt+HLD0AAM6QhTdj4U3dm7yCQnsBfuxVbTQoIZIsPAC3VqFV6w8dOqQPP/zQjmU/6x3x8LDZ/sGDB//mNlP4bvjw4bZL/ffff/+7FfxKKygoUIcOHdSzZ09NnDix3Bn56OjoclULhPvKLyjUgl9S9fGyZM3fvF8Fv75Tgmr4aEiH+rq+c4ya1wt09G4CAFAtrd+dZivSb953zK5fER+hp69qY4vaAoC7V60/q4z8pEmT7BOXZq4DfPvtt3r33Xc1a9Ys+wfPJZA/HRPEX3PNNdqxY4fmz59/1oG1ufjQqVMnbd269bT38fPzswtQmpenh3q3rGeXfWlZ+mxlsu16bzL2k37YaZcLGgbbgH5AfKRq+JKlBwCgsmXn5Wvif7fqrQWJ9qK7KVZrAngzSw0AVBdnNf3ciBEjSgLenTt36sknn7Rj0a+44gpbFX727Nk2Y1/RQbwJwufNm6fQ0NCzfg5zoWHNmjWKjOTDHecuIshff7k0VgtHXqL3bu1ki+CZQP+nXUdsZdzO4+bpyVkbSsbmAQCAircm+agGTlysN77bboP4QW2j9O0DPQniAVQ7Z5WRN93Pp0+frnfeeUdLliyxFefNVHPXX3+9Hb8eFxd3Vn/cFKXbtm1bybrJupugOyQkxM4bP2zYMDsF3ZdffmkL1xVfJDC3F1fON3PL169f345zN8xY9wsvvNCOhzddE0x3evOcb7zxxlntG3AqJnjv1SLcLqnpWfrsp92asiJJyYdP6IOlu+zSLrqObugco4FtIxXge9b1JAEAwEmycvP16rxf9K+FiXaoW1gtXz07uI2dgQYAqqOzGiMfFhZmg/WbbrrJjlkPDg622318fGwxubMN5M1490suueSUmf8xY8aocePGp3zcd999p169etnfzc9GjRrZAnfGAw88YC82mKDfDAMw1erNc3Xt2rVSxiYABQWF+mH7QU1Znqw5G/fZYjtGLT9vO+WN6Xrfpn7ZISkAAKB8TO+3kVPXavuBDLs+uF2UnhrUWsE1Tz8dMgC4orOJQ88qkDeBe0JCgg3kr7322pInP9dA3lkRyONcHTiWrWmrdmvK8iTtPJRZst1k6Z++qrUSGpx56kQAAPC/LPzL326xM8eYs9W6gX4ad3W8+sTV4xABcEtnE4ee1Rj5lJQU/fnPf9bkyZMVERGhoUOH2krzjpyWDnAm5iTjzoubav5fe+njP3WxY/d8vTztmL6r31yiF77ZbE9MAADA6a3ceVhX/G2R/rWoKIg3s8XMfaAnQTwAnEtGvrTt27fbKvbvv/++9uzZY8fJ33LLLerdu7e8vFy7ejcZeVR0lv6ZLzfp87V77Xqz8Fp6cViC2scUDU0BAABFTuTk68U5WzRpSVEAX6+2n8YPibczyACAu0uvrK71p5unfc6cOXb6uS+++EKBgYE6ePCgXBmBPCqDGT//2IwNOng8W54e0p8uaqIH+jSXv49rX/gCAKAiLEs8pJHT1mnXr0PTrunYQI8NiFNQDR8OMIBqIb0qA/nSDhw4oA8//FAPPvigXBmBPCrL0cwcjf1ik2as3mPXm9StabPzFzQM4aADAKqljOw8Tfhms95fusuuRwb52yy8mSEGAKqTdEcF8u6CQB6Vbd6m/Xp0xnqlHsuWKTFxe/fG+mvfFqrhS3YeAFB9LNl+UKOmrbPTuBrXd47W6CtaqbY/WXgA1U96ZVatL09hu8OHD8uVEcijKqRl5uqZ2Zs09afddr1RaIAmDGurzo3JzgMA3Nvx7Dw9//XP+ujHJLtev04NPT80XhfF1nX0rgGAS8Sh3mfzxK+99tr57huAXwUF+Oil4W01ICFSo6ett9PVXfv2Uo3o2kgj+7dQgO9ZvT0BAHAJi7cWZeH3HC3Kwt/YJcZm4Wv58b0HAOVF1/pTICOPqpaelavnvvxZn6xMtusxISY7n6ALm4TSGAAAt3AsK1fjvtqsycuLsvANgmtowtAEdWsW5uhdAwCnwBj5KjyAQEVa8MsBjZ62TnvTsuz6zV0balT/lqpJlgIA4ML4fgOA30cgf54I5OFIZCwAAO7c4+yFoQnq2pQeZwBwMgL580QgD2fAGEIAgCv7bnOqRk9fr33pWXaGllu6NdLD/agBAwCnQyB/ngjk4UxVfV/4erM+/LFobl2q+gIAXGFWlqe/3KRpq4pmZWkcVtPWfenUiFlZAOBMCOTPE4E8nA3z7AIAXMG8Tfv16Iz1Sj2WbbPwt3dvrL/2baEavl6O3jUAcHqVHsjn5+frvffe03//+1+lpqaqoKCgzO3z58+XKyOQhzPKyM7Ti3O26L0lO+16ZJC/xg+JV68W4Y7eNQBANXc0M0djv9ikGav32PUmYTX14vAEXdCQLDwAOHwe+WL33XefDeQHDBigNm3ayMNccgVQqUzl+jFXttblbSI0cto67TqUqVsmrdDwCxro8YFxCqrhQwsAAKrcnI379NiMDTp4PFueHtKfLmqiB/o0l78PWXgAqCznlJEPCwvTBx98oCuuuELuiIw8nN2JnHybnZ+0ZIfMO7hebT+bne/dsp6jdw0AUE0czsjRmM836vO1e+16s/BaenFYgtrHBDt61wDA7eNQz3P5A76+vmrWrNm57h+A82TGGj45KE6f3dHVFhHan56t295bqQc/XWOLDAEAUJm+Wp+ivq8usEG8ycL/X6+m+vIvPQjiAcCZM/Ivv/yyEhMT9frrr7tlt3oy8nAlWbn5emXuL3pnUaIKCqW6gX4ad3W8+sSRnQcAVCzTff6pWRs1e32KXW9ez2Th26ptdB0ONQA4e7G7q6++Wt99951CQkLUunVr+fiUHZs7ffp0uTICebiiVUlH9PBna7X9QIZdH9wuSk8Naq3gmr6O3jUAgIszp4tfrkvRU59vtF3qvTw9dFevprqndzP5eTMWHgBcothdnTp1bDAPwHl0iAnW7Hsv0mvzturthds1c81eLd52SM8Obq3+bSIdvXsAABd14Fi2npi5Qd9s3GfXW0YE6qXhbdWmfpCjdw0Aqq1zysi7OzLycHVrko/a7PzW1ON2fWBCpMZe2VqhtfwcvWsAABdhThHNGHiThT+amStvTw/dfUkzu/h6n1OZJQCAI7vWuzsCebiD7Lx8TfzvVr21IFH5BYUKremrp69qowEJZOcBAGeWmp6lx2Zu0NxN++16XGRtOy986yiy8ADg0oH81KlT9emnnyopKUk5OTllblu1apVcGYE83Mn63Wl66LO12rL/mF2/Ij7CBvRhZOcBACcxp4UzVu/R2C82Ke1Erny8PPSX3rG2Kr2PF1l4AHDp6ecmTpyoW2+9VeHh4Vq9erU6d+6s0NBQW8n+8ssvP9f9BlAJ4hsE6Yu/9NC9l8babpFfrd+nPq8UTRlEhxwAQLF9aVm6/X0zlelaG8S3qV+75PuDIB4AnMs5BfJvvvmm3n77bTv9nJlTfuTIkZo7d67uvfdee/WgvBYuXKhBgwYpKirKTmM3c+bMMrebIGPMmDH29ho1aqhXr17auHHj7z7vtGnTFBcXJz8/P/tzxowZ5/IyAbdhxjI+2Ke5Zt7dXa0ia+tIZq7unbxad370k1KPZTl69wAADmTOtz5dmaw+ry7Q/M2p8vXy1MP9WmjGXd3VMuLMGSEAgAsF8qY7fbdu3ezvJsA+dqyoy+4f/vAHTZ48udzPk5GRobZt29oLAqcyYcIEvfLKK/b2FStWKCIiQn369Cn5e6eydOlSXXvttXZf1q5da39ec801WrZs2Vm/TsDdmArDs+7urgcua26z83M27lffVxdq5uo9ZOcBoBrae/SEbpm0QiOnrtOxrDy1bRCkL+/tYQvakYUHAOd1TmPkmzRpYsfId+jQQZ06ddIf//hH3XHHHfr222913XXX6fDhw2e/Ix4eNnM+ePBgu252y2Ti77//fo0aNcpuy87OVr169fTCCy/Yv3cqJog3Ywu+/vrrkm39+/dXcHBwuS8yMEYe1cHPKel6eOpabdiTbtcvaxWu566OV73a/o7eNQBAJTPnWZ+sSNZzs3/Wsey8kp5bf+zRWN6MhQcA9xwj37t3b33xxRf299tvv10PPPCAzZSbILqi5pffsWOH9u3bp759+5ZsM13lL774Yi1ZsuSMGfnSjzH69et3xseYCwTmoJVeAHdnutibbpMP9W1uixnN+znVjp2f+tNusvMA4IbMDCZrk4/qje+26eo3l+iR6ettEN8+po6+uvci3XlxU4J4AHAR3ufyIDM+vqCgwP5+5513KiQkRIsXL7bj3c16RTBBvGEy8KWZ9V27dp3xcad6TPHzncr48eM1duzY895nwNWYbpP39I5Vn7gIm51f92uF+9nr9mr8kARFBJGdBwBXzrrvOJihH7Yd1OJtB7V0+yGlZ+WV3O7n7amH+rbQbT0ay8vTw6H7CgCogkDe09PTLsXMGHSzVAbT5f7kL6WTt53vY0aPHq0HH3ywZN1k5KOjo895nwFX0yIiUNP/r5veXpSo1+Zu1XdbDtiiR08MiNPwjg1+9z0HAHAOpoCpCdgXbz1oA/i9aWULmgb6eevCpqHq0SxMl8XVU/06NRy2rwCAKg7kjUWLFumf//yntm/fbsfL169fXx9++KEaN26sHj166HyZwnaGyaRHRkaWbE9NTf1Nxv3kx52cff+9x5gu+2YBqjMzJvKuXs3Up1U9PTx1ndYkH9XIaev05foUPT8kXlGc7AGA0zmenaflO0zgfsgG7lv2ly0IbCrQd2hYxwbu3ZuFKb5+EN3nAaC6BvJmejdTDf7GG2+088ibMeaGqSY/btw4ffXVV+e9Y+aCgAnKzbR27du3t9tycnK0YMECW+zudLp27WofY8btFzNF+Iqr7AM4s9h6gZr2f9307uJEvfTtL1r4ywFb2f6xAa10XadosvMA4EC5+QX2QqvJuC/ZflCrk44qr6Bs3eLWUbVt4N6tWZg6NQpWgO85520AAE7qnD7Zn332Wb311lu6+eabNWXKlJLtJlh++umny/08x48f17Zt28oUuFuzZo0dcx8TE2Mr1psLA7GxsXYxvwcEBOiGG24oeYzZB9MbwIxzN+677z717NnTBvtXXXWVZs2apXnz5tkx/ADKx4yV/HPPpurdsp5GTl2rVUlHNXr6es1el6Lnh8arQXAAhxIAqoAZHmiy7D9sK8q4L0s8pIyc/DL3iQkJUPdmoTbj3rVJqEJr0csQANzdOQXyW7ZsscHyyUyJ/KNHj5b7eVauXKlLLrmkZL14nPqIESP03nvvaeTIkTpx4oTuuusuHTlyRF26dLHZ9cDAwDJz2pcer28uJpiLC48//rieeOIJNW3aVJ988ol9LICz0yy8lj67s5sm/bBDL87ZYosl9Xt1oR65opVu7BwjT4ojAUCF23P0hA3ai5ZDOni8qOdjseAAH5ttt93lm4YpJpSLqwBQ3ZzTPPImODbj4y+77DIbVK9du9bOLf/BBx/o+eef16ZNm+TKmEce+C1T+dhk51fsPGLXTdZnwrAERYdwAgkA5yMtM1dLE4sqyy/ZdkiJBzPK3O7v46nOjU2BulB1axqmuMjaXEgFADd0NnHoOWXk77jjDtuF/d///rcdL7t37147f/tDDz2kJ5988lz3G4ATaxxWU5/8uaveX7pTE77ZoqWJh9TvtYUa1b+l/nBhQ04qAaCcsnLztWrXERu4m6z7+j1pKj3M3XR2SmjwvwJ1plidn7cXxxcAcH4ZeeOxxx7Tq6++qqysomlNTNV3E8g/88wzcnVk5IEz23XIZOfXadmOw3a9c+MQTRiaoEZhNTl0AHCS/IJCbdqbXhK4r9h5WNl5BWXu07RuzZLAvUuTUAXV8OE4AkA1k34WGflzDuSNzMxM242+oKBAcXFxqlWrltwBgTzw+woKCvWfZbs0/uvNyszJt10/R/ZrqVu6NSI7D6BaM6dWuw5lFnWV326WQzqamVvmPuGBfiWBu1kigvwdtr8AgGoWyLsrAnmg/JIPZ2rUtHX2RNXo2DDYjp1vUtc9LuwBQHmYgnTmc/CHrUVj3U3ButJq+XnrwiYhNmg3AbwpJmqGJwIAUOmB/G233Vau+5mx866MQB44O+Zj5OPlSRo3+2c7LZKft6ce6ttCt/VobKeyAwB3k5Gdp+U7D5cE7pv3HStzu4+Xh9rHBJdk3RMaBMnH63+z7AAAUGWBvJnmrWHDhmrfvr09cT+dGTNmyJURyAPnZveRTDvf/KKtB+16+5g6enFYW5t5AgBXlptfoHW7j9rp4EzgvjrpiHLzy54LtYqsre5NQ9U9NkydG4Wopt851RQGAFRT6ZUVyJv53M0c7TExMTY7f9NNNykkJETuhkAeOHfmI+WTFcl6dvbPOp6dJ19vTz3Yp7n+2KOxvMlGAXChz7JtqcdLCtT9mHjYfqaVVr9OjaKMe2yYujUNVVgtP4ftLwDA9VXqGPns7GxNnz7ddp9fsmSJBgwYoNtvv119+/Z1m7FeBPLA+dt79ITNzi/45YBdb9sgSC8Ob6vm9QI5vACcUkraCZtxN4G7WVKPZZe5vU6Ajw3Yi8e5x4QEuM25DwCgGhW727Vrl9577z198MEHys3NtRXs3aFyPYE8UDHMx8vUn3br6S836VhWnny9PHXfZbG6o2cTsvMAHC7tRK5+TDykJduKxrlvP5BR5nZT78NMr9mtaVHgHhdVm7ofAACniEPPa/CWuQptFnOybqagA4CTPyOGd4zWRbF19eiM9Zq/OVUvztmibzbs04vDE9Qy4swfUABQkbLz8rVqlxnnXhS4mzHvBaXSGSa5nlA/qCTj3qFhsPx9vGgEAIDTOa+u9YsXL9bAgQN16623qn///rYYnjsgIw9UPPNRM2P1Ho35fKPSs/JsRee/9I7V//VqSiVnAJWioKBQm1LSi7rKbz+k5TsOKSu3bOKhSVjNkrncuzYJVVCAD60BAHDfYncmeDfF7kJDQ+VuCOSBypOanqVHZ2zQvJ/32/XWUbVtZXvTZRUAzlfy4UybbTeL6TJ/JDO3zO2mIF2PZqHq9mvwbgrWAQDg9tPPmSDeTD93puIuJmPvygjkgcplPnY+X7tXT32+UUczc+Xt6aG7L2lmF1PlHgDK63BGjpZsLypOZ4L35MMnytwe4OulC5v8r0Bd83q1KFAHAKheY+RvvvlmvvwAnDdzIfCqdvXVtWmonpy5Ud9s3Ke//Xer5mzcp5eGt1Wb+kEcZQCndCInX8t3Hi6pLL9xb3qZ282FwfYxdUoC97bRdRi+AwBwO+dVtd5dkZEHqo75CJq9PkVPztpoM2tenh66q1dT3dO7mfy8KTIFVHd5+QVavyetJONuitXl5Jcd596iXmBR4B4bqs6NQ1XL77xq+QIA4N7Tz7krAnmg6h08nq2nZm20Qb1hur+a7HxCgzo0B1CNmNMSMw1cceD+4/ZDOpadV+Y+UUH+vwbuYbZnT3igv8P2FwCAikIgX4UHEEDF+mp9ip6YuUGHfs3O/7lnE913aSxTQAFubH961q9d5Q/Zn/vSs8rcXtvf287l3j22qLt8o9AAhvoBANwOgXwVHkAAFc90sTeF8L5Yu9euNwuvpReHJah9TDCHG3ADx7JytSzxsM24m8B9a+rxMrebopcdGwaXjHM3dTPMhT0AANxZOl3rq+4AAqg832zYp8dnbrDd7s05/J8uaqIH+jQnOw+4mJy8Aq1OOlLSXX7t7jTlF/xvZJ+ZCKdNVFBJ4N6xUTDvcwBAtZNOIF91BxBA5TqSkaOnv9ykGav32PUmdWva7PwFDUM49ICTKigo1OZ9x4q6y28/aLPvJ3Lzy9zHdI83gbtZujYJVXBNX4ftLwAAzoBAvgoPIICqMW/Tfj06Y71Sj2Xb7N3t3Rvrr31bqIYvle0BZ7FhT5o+Xp6kORv22ToXpYXW9FU3m3EPtePdo0MCHLafAAA4IwL5KjyAAKpOWmauzc5PW7W7JKM3YVhbdW5Mdh5wlOPZefp8zV5NXp5kp4krVsPHS12ahNiu8ibrbqaI82ScOwAAp0Ugf54I5AHn9t3mVI2evt5WtjbZ+RFdG2lk/xYK8GXuaKCqpohbtzvNBu+fr92rzJyibvO+Xp7q1yZC13RsoC6NQ23ROgAAUD4E8ueJQB5wfmkncvXc7E36dGVRdj4mxGTnE3Rhk1BH7xrgttKzcjVr9R5NXp6sTSnpJdtN7YrrO8VoSIf6Cq3l59B9BADAVblVIN+oUSPt2rXrN9vvuusuvfHGG7/Z/v333+uSSy75zfaff/5ZLVu2LNffJJAHXMeCXw7okWnrlJJWNO/0zV0balT/lqrpR3YeqAjmNGFV0lFNWZ6kL9ellBStM9n2K9pE6PrOMXZ4i4fpHgMAAM7Z2cShTn+mu2LFCuXn/6/S7YYNG9SnTx8NHz78jI/bsmVLmRdft27dSt1PAI5xcfO6+vaBnhr31WbbzfeDpbs0f3OqJgxNsIW1AJx7TYoZq3fb7PuW/cdKtseG17LBu8m+1wmg0jwAAI7g9IH8yQH4888/r6ZNm+riiy8+4+PCw8NVp06dSt47AM4g0N9H44fEa0B8pEZNW6fdR07ohneW6cYuMRp9RSvVIjsPlDv7vnLXEU1elqTZ61OUnVdgt/t5e2pgQpRu6BKtDjHBZN8BAHAwpw/kS8vJydFHH32kBx988HdPItq3b6+srCzFxcXp8ccfP2V3+2LZ2dl2Kd2lAYDr6REbpjkP9NTzX/+sj35M0n+WJen7LQf0/NB4XRRLrxzgdI5k5Gi6HfuepG2px0u2t4wItNn3we3qKyjAhwMIAICTcPox8qV9+umnuuGGG5SUlKSoqKjTdqlfuHChLrjgAhucf/jhh3rrrbfs2PmePXue8jFjxozR2LFjf7Od6ecA17Vk+0GbnU8+fMKuX9852mbna/sTjACG+fpftuOwDd6/3rBPOb9m3820cYPaRtoAvl10HbLvAABUEbcqdldav3795Ovrqy+++OKsHjdo0CB7IvL555+XOyMfHR1NIA+4uIzsPE34ZrPeX1pUMDMyyN92we/VItzRuwY4zKHj2Zq+qij7nngwo2R7XGRtXd8lRle1i+KCFwAADuBWxe6Kmcr18+bN0/Tp08/6sRdeeKHtkn86fn5+dgHgXkzl+rFXtdHlv46d33UoU7dMWqHhFzTQ4wPjFFSD7Dyqh4KCQv2YeEgfL0/SnI37lJtfdA2/pq+XrmwXZbPv8fWDyL4DAOAiXCaQnzRpki1gN2DAgLN+7OrVqxUZGVkp+wXA+Zm55b++7yK9NOcXTVqyQ5/9tFsLtx6w2fneLes5eveASnPgWLam/rRbU1Yk2QtZxRIaBNngfVDbKIpBAgDgglwikC8oKLCB/IgRI+TtXXaXR48erT179uiDDz6w66+99pqde75169YlxfGmTZtmFwDVV4Cvt54cFKcr4iP08NR12nEwQ7e9t9JOofXUwNYU8oJbZd8Xbztou87P3bRfeQVF2Xcze8Pg9lG6rlOM2tQPcvRuAgAAdw/kTZd6U+Dutttu+81tKSkp9rZiJnh/6KGHbHBfo0YNG9DPnj1bV1xxRRXvNQBn1LFRiL669yK9MneL3lm8w44VXrT1oMZdHa8+cWTn4bpS07NsbxOTfS8u8miYgnU3dI7RwLaR9oIWAABwfS5V7M4ZiwwAcF0/7Tqih6euVeKBooJfF8WG2UJ43ZuFqkW9QMYLw+nlFxTaYSJm3vf/bk6160agv7eGtK+v6zrHqFUk32MAALgCt61aX1UI5IHqIys3X6/O+0X/WpioX2MgK6yWnw3ouzcLs0v9OjUcuZtAGfvSsvTpymR9siJZe47+L/vesWGwDd4HxEeqhq8XRw0AABdCIF+FBxCAe9iWelzzN+/X4m2HtHzHIWXlFs2pXaxJWE11axaqHs3C1LVJGGPqUeXy8gu04JcDduz7/M2pJReezOwLptaDKV7XvF4gLQMAgIsikK/CAwjA/WTn5WvVrqP6YdtBWzRs3e6jZbL1nh6yU3WZTL0J7Ds0DJa/D9lPVA6TcTeZ989WJislLatke+fGIXbse/82Efz/AwDADRDIV+EBBOD+0k7kalnioZLAfvuvY+qL+Xl7qlOjkJLAPi6qtrxMtA+cR/bdZN1N9v37Xw6oeBBccICPhl3QQNd2ilGz8FocXwAA3AiBfBUeQADVT0raCf2w7ZCW/BrYpx7LLnN7nQAfdWv66/j6pmFqGBpA4TyUS/LhTJt9N+PfS/+/6tokVNd3iVG/1vXk503vDwAA3BGBfBUeQADVm6kXasbXm4DeZOx/TDys49l5Ze5jCuWZTH332DAb4JtCekCx3PwCzdu0Xx8vT7L/j4qz76E1fTWsYwM773vjsJocMAAA3Fw6Veur7gACwMlB2brdaSXd8FcnHVFuftnJQcx0YN1Nxj42TF0ahzC3dzW161CGptix77t18Pj/su9mGkQTvPeJqydfb0+H7iMAAKg6BPJVeAAB4EwysvO0fOdh/bD1oH7Yfkg/p6SXud3Hy0PtY4KLMvbNQpXQoI58vAje3FVOXoG+3bTPjn03wzOK1Q3003A79j1aDUPJvgMAUB2lk5GvugMIAGfDZF6XbD9kA3uTsS89B7hRy89bFzb5X+E8U9DMw4PCea4u8cBxm32f9tNuHcrIsdtMs/aMrWunjbu0VTgXcAAAqObSCeSr7gACwPmMr086nFkyvt4E+Eczc8vcJzzQzwb03X7N2EcG1eCAu4is3HzN2ViUfTe1E4rVq+2naztGa3jHaEWHBDh0HwEAgPMgkK/CAwgAFSW/oFCb9qbbwH7J9oNavuOwsvMKytynad2av3bDD9OFTUNV29+HBnAy21KPafLyZE1ftVtHfr0wY2Yj7NUi3GbfL2lRV94MnwAAACchkD9PBPIAnCWju2rXkZKM/bo9aSUVzYuDQzOmvjiw79CwDlOTObCtvt6QosnLkm1NhGKRQf523Ps1HaMVVYfeFAAA4PQI5M8TgTwAZ5SWmauliYdsUG+WxIMZZW739/FU58ah6tEsVN2ahikusrY8TbSPSrNln8m+J9nse3pW0bSDXp4euqRFuG7oEq2Lm4fbdQAAgN9DIH+eCOQBuAJTKK84qDcV0EtPYWaE1PRV16YmsC8qnMd47IpxIidfX67bawP4VUlHS7bXr1ND13UqGvseEeRfQX8NAABUF+kUu6u6AwgAzlI475f9x0u64f+YeEiZOfll7hMTEmC74Hf/NWNvAn2Un6lfMGVFkmas3qNjv2bfvT09dFmrerq+S4y9WEL2HQAAnCsC+fNEIA/A1eXmF2ht8tGSwH510lHlFZQaYC+pdVTtkor4nRuFqIavl8P211llZOfZ7PvHy5Pt8Sx9UeS6ztEadkEDhQeSfQcAAOePQL4KDyAAuILj2XlavuOQFm8tGmO/Zf+xMrf7ennaYnnFhfPi6wdV68rqG/ak6ePlSfp8zV577AwfLw/1jYuwlee7NQ2l/gAAAKhQBPJVeAABwBWlHsvS0u0msC/K2O9Nyypze6C/t7o2Cf21K36YnfbOw8O9i7aZgN0E7mbs+/o9aSXbG4UG2OB96AUNFFbLz6H7CAAA3BeBfBUeQABwh/H1Ow5mlBTNM3PYF1dgLxZR298G9D1iQ9W9aZjCa/u7zWtftzvNBu+fr91bUlfA9FDo3ybCdp83FzTc/SIGAABwPAL5KjyAAOBu8gsKbdfy4vH1K3cdUU5eQZn7xIbXKgrsm4WpS5MQBfr7yJWkZ+Vq1uo9mrw8WZtS0ku2m54HJvs+pEMDigECAIAqRSBfhQcQANxdVm6+Vu48UhLYb9ibpsJSdfNMpfZ20XWKuuE3DVX7mGD5ens6ZfbdTBc3ZXmSvli3V1m5RRcnzL4OiI+0AXynRsFk3wEAgEMQyFfhAQSA6uZIRo6WJhYVzTPLzkOZZW4P8PVS58YhJYXzWtQLdGhhuLTMXM1Yvdtm30sX+TO9Coqy7/VVJ4Cp+AAAgGMRyFfhAQSA6i75cKYdV7/YjK/fdlCHMnLK3B5Wy9fOW2/mrzeBfYPggCrJvpshAZOXJWn2+hRl/zo0wM/bUwMTonRDl2h1iCH7DgAAnAeBfBUeQADA/xQUFNqst8nUm674yxIP60RuUQG50lXgi8fXd20aWqHZcNNbYLod+56kbanHS7a3jAjUDV1idFW7+gqq4Vrj+QEAQPWQfhZxqEehSVs4qTFjxmjs2LFlttWrV0/79u077WMWLFigBx98UBs3blRUVJRGjhypO++886z+LoE8AFQMUyRvddIR/bC9qCv+muSjtpheMVMMvk1UUElg37FRsPx9vM7qb5ivsWU7Dtvg/esN+0oK89Xw8dKVbaNs5Xkzhp/K8wAAwJmdTRzqLSfXunVrzZs3r2Tdy+v0J3g7duzQFVdcoT/96U/66KOP9MMPP+iuu+5S3bp1NXTo0CraYwBAMVNIrkuTULs82Ke5jmXl2ix9ceG8ranH7ZztZnlrwXZ7f1NwznTFN4F9m/pBtpjeqRw6nq3pq4qy74kHM0q2t46qbce+X9UuyuWq6QMAAJSH0wfy3t7eioiIKNd933rrLcXExOi1116z661atdLKlSv10ksvEcgDgBMwgfVlcfXsYuxPzyoaX7+1KGO/Lz3LzmVvlhfnbFFtf++i8fWxRYF9w5AA/Zh4SB8vT9KcjfuUm1+U3a/p66Ur29XXDZ1jFN8gyMGvEgAAoJoH8lu3brVd5P38/NSlSxeNGzdOTZo0OeV9ly5dqr59+5bZ1q9fP7377rvKzc2Vj8+pMzPZ2dl2Kd2lAQBQ+erV9tfV7RvYxXSR334go6QavqmMn56Vp2827rNLccCekfO/MfdtGwTpus4xGtQ2SrX8nP4rDQAAoEI49VmPCdw/+OADNW/eXPv379ezzz6rbt262fHvoaGhv7m/GTtvxtCXZtbz8vJ08OBBRUZGnvLvjB8//jdj8QEAVcuMYW8WXssuI7o1Ul5+ge1yX1w4b9WuozaINwH74PZRuq5TjO16DwAAUN04dSB/+eWXl/weHx+vrl27qmnTpnr//fdtQbtTObmYUXEtvzMVORo9enSZ5zMZ+ejo6Ap4BQCAc+Xt5an2McF2uad3rE7k5Gtr6jEb6Af4OvXXFwAAQKVyqTOhmjVr2oDedLc/FTOW/uSK9qmpqXac/aky+MVMt32zAACcVw1fLyU0qOPo3QAAAHA4T7kQM479559/Pm0XeZOxnzt3bplt3377rTp27Hja8fEAAAAAALgSpw7kH3roITsvvJlWbtmyZRo2bJjt9j5ixIiSLvE333xzyf3NfPG7du2y3eRNwP/vf//bFrozzwMAAAAAgDtw6q71u3fv1vXXX28L1Zm54C+88EL9+OOPatiwob09JSVFSUlJJfdv3LixvvrqKz3wwAN64403bLX7iRMnMvUcAAAAAMBteBQWV4NDCZP1DwoKUlpammrXrs2RAQAAAAA4TRzq1F3rAQAAAABAWQTyAAAAAAC4EAJ5AAAAAABciFMXu3OU4rIBZowCAAAAAACVrTj+LE8ZOwL5Uzh27Jj9GR0dXdFtAwAAAADAGeNRU/TuTKhafwoFBQXau3evAgMD5eHhIWe/amMuOCQnJ1Nh30nRRq6BdnINtJPzo41cA+3kGmgn50cbuYZ0F4mZTCbeBPFmGnVPzzOPgicjfwrmoDVo0ECuxPyHdOb/lKCNXAXvJddAOzk/2sg10E6ugXZyfrSRa6jtAjHT72Xii1HsDgAAAAAAF0IgDwAAAACACyGQd3F+fn566qmn7E84J9rINdBOroF2cn60kWugnVwD7eT8aCPX4OeGMRPF7gAAAAAAcCFk5AEAAAAAcCEE8gAAAAAAuBACeQAAAAAAXAiBPAAAAAAALoRA3oW9+eabaty4sfz9/XXBBRdo0aJFjt6lam3Pnj266aabFBoaqoCAALVr104//fRTye2FhYUaM2aMoqKiVKNGDfXq1UsbN2506D67u4ULF2rQoEH2mHt4eGjmzJklt+Xm5mrUqFGKj49XzZo17X1uvvlm7d27t8xzZGdn6y9/+YvCwsLs/a688krt3r3bAa+meraTcfz4cd1zzz1q0KCBfe+0atVK//jHP8rch3aqXOPHj1enTp0UGBio8PBwDR48WFu2bDnt/e+44w7blq+99hrtVIXM+yIhIUG1a9e2S9euXfX111+f1fcQ7yXHtpHx888/2++aoKAg+5678MILlZSURBs5+DPQfKbdf//9dp1zCOdvo+pw/kAg76I++eQT+x/1scce0+rVq3XRRRfp8ssvL/NBj6pz5MgRde/eXT4+PvYLedOmTXr55ZdVp06dkvtMmDBBr7zyil5//XWtWLFCERER6tOnj44dO0ZTVZKMjAy1bdvWHvOTZWZmatWqVXriiSfsz+nTp+uXX36xH+ClmffZjBkzNGXKFC1evNh+KQwcOFD5+fm0WxW0k/HAAw/om2++0UcffWRPcM26+dKdNWsW7VRFFixYoLvvvls//vij5s6dq7y8PPXt29e23cnMhZhly5bZYPFkvJ8qlzlZff7557Vy5Uq79O7dW1dddVVJsF6e7yHayLFttH37dvXo0UMtW7bU999/r7Vr19rvKZO0oY0cw7xX3n77bXsBphjnEM7fRtXi/KEQLqlz586Fd955Z5ltLVu2LHzkkUcctk/V2ahRowp79Ohx2tsLCgoKIyIiCp9//vmSbVlZWYVBQUGFb731VhXtZfVmPu5mzJhxxvssX77c3m/Xrl12/ejRo4U+Pj6FU6ZMKbnPnj17Cj09PQu/+eabSt/n6uhU7dS6devCp59+usy2Dh06FD7++OP2d9qp6qWmptq2WrBgQZntu3fvLqxfv37hhg0bChs2bFj46quvltxGOzlGcHBw4TvvvFOu7yHayLFtZFx77bWFN91002nvSxtVrWPHjhXGxsYWzp07t/Diiy8uvO+++057X84hnK+NWrv5+QMZeReUk5Nju2ybbEhpZn3JkiUO26/q7PPPP1fHjh01fPhw2+20ffv2+te//lVy+44dO7Rv374ybebn56eLL76YNnMiaWlptltWcU8K8z4z3edKt5vJMrZp04Z2q0ImO2XeY2b4ion1v/vuO9t7ol+/frSTA98rRkhISMm2goIC/eEPf9DDDz+s1q1b/+YxvJ+qlskmmQyT6TVhum+X53uINnJsG5n30OzZs9W8eXP7+WbOJ7p06VJmuBFtVLVMT6QBAwbosssu+937cg7hfG3Uw83PHwjkXdDBgwfth3+9evXKbDfr5ksaVS8xMdGOuYmNjdWcOXN055136t5779UHH3xgby9uF9rMeWVlZemRRx7RDTfcYMctFrebr6+vgoODy9yX91rVmjhxouLi4myXVNMe/fv3tzVCzBc07VT1zMnQgw8+aI+/Odkp9sILL8jb29t+9p0K76eqsX79etWqVcsG6ea7yHQZNe+f8nwP0UaObaPU1FTbrdd0vTefc99++62uvvpqDRkyxA5voY2qlrnIYobembHXv4dzCOdso4lufv7g7egdwLkzmcOTT65O3oaqYa6im4z8uHHj7LrJyJvxbia4NwXUitFmzslcjb3uuutsO5oP+N/De61qmS9iMzbbXFVv2LChLY531113KTIy8oxZEtqpcpjCQevWrbNjCYuZrMbf/vY3e0J1tt9DtFPFatGihdasWaOjR49q2rRpGjFiREkQeK7fQ7RR1bRRcW8wM2bejOU1TOFckxl86623bO8J2qhqJCcn67777rMXU0rXJzgVziGct40muvn5Axl5F2SqKnp5ef3mSpG5knvylXZUDfOBYK74lWYqYxYXHzQFhQzazPmYL+BrrrnGdjs1RbyKs/HF7WaGsphihqXxXqs6J06c0KOPPmoLdJnK9qaQjQkkr732Wr300ku0UxUzRYLMCZHpnmgyHMXMrCnmfRETE2Oz8mbZtWuX/vrXv6pRo0a0UxUy2aVmzZrZi8smS2UKSZqLLOX5HuIzz7FtZM7vzHvn984n+F6qfObipHlvmFmhij/TzMUWExia34sLoXEO4bxtlJGR4fbnDwTyLvoFYP7TmqCjNLPerVs3h+1XdWYq1p88FZMZg2Ou/hlmmkDzYVG6zcwHh/nAoc0cp/gLeOvWrZo3b56dOrA08z4zMxGUbreUlBRt2LCBdqvCNjKLp2fZrytzMdP0oKCdqobJTpgTIDO7w/z58+1nWmlmbLzJ0pssY/Fixhma8fJmuBHt5Ni2M9Mrled7iM88x7aROb8z0zye6XyCNqoal156qR0CUfozzVx4ufHGG+3v5juIcwjnbqP8/Hz3P39wdLU9nBtTXdFUWXz33XcLN23aVHj//fcX1qxZs3Dnzp0cUgcwlUq9vb0Ln3vuucKtW7cW/uc//ykMCAgo/Oijj0ruYyoFm+rA06dPL1y/fn3h9ddfXxgZGVmYnp5Om1ViJdPVq1fbxXzcvfLKK/Z3U5U+Nze38Morryxs0KBB4Zo1awpTUlJKluzs7JLnMLNDmPvMmzevcNWqVYW9e/cubNu2bWFeXh7tVgXtZJgqtKby7HfffVeYmJhYOGnSpEJ/f//CN998k3aqIv/3f/9nP7++//77Mu+VzMzM0z7m5Kr1Bu+nyjV69OjChQsXFu7YsaNw3bp1hY8++qitvvztt9+W+3uINnJsG5m2Med3b7/9tj2f+Pvf/17o5eVVuGjRItrIwUpXROccwjldfFLVenc/fyCQd2FvvPGGPVHy9fW1UymcPA0QqtYXX3xR2KZNm0I/Pz87FaD5Ei7NTP3z1FNP2el/zH169uxpT6RQecwHtwkMT15GjBhhT6JOdZtZzOOKnThxovCee+4pDAkJKaxRo0bhwIEDC5OSkmi2KmonwwSMt9xyS2FUVJT9Am7RokXhyy+/bN9TtFPVON17xZwUnU0gz/upct12220l5wV169YtvPTSS0sCxPJ+D9FGjm0jwyRpmjVrZj/vTEAxc+ZM2sjJgkTOIVwjkE9x8/MHD/OPo3sFAAAAAACA8mGMPAAAAAAALoRAHgAAAAAAF0IgDwAAAACACyGQBwAAAADAhRDIAwAAAADgQgjkAQAAAABwIQTyAAAAAAC4EAJ5AAAAAABcCIE8AACVaMyYMWrXrh3H+Bw1atRIr7322hnvk5OTo2bNmumHH344r+N8yy23aPDgwaoonTp10vTp0yvs+QAAKEYgDwDAOfLw8DjjYgLDhx56SP/973+r/Bjv3LnT7sOaNWvc/qLD22+/rYYNG6p79+7n9Tx/+9vf9N5771XYfj3xxBN65JFHVFBQUGHPCQCAQSAPAMA5SklJKVlM1rh27dpltpnAsFatWgoNDa1Wxzg3N7dK/97f//53/fGPfzzv5wkKClKdOnVUUQYMGKC0tDTNmTOnwp4TAACDQB4AgHMUERFRspgg0GTAT952cpa7uPv2uHHjVK9ePRs4jh07Vnl5eXr44YcVEhKiBg0a6N///neZv7Vnzx5de+21Cg4OthcGrrrqKpt1L6/vv//e7p/pHdCxY0cFBASoW7du2rJli73dZKLNfqxdu7akR0FxdtoEo3/+858VHh5uL1b07t3b3q9Y8Ws0+9ykSRP5+fnpn//8p+rXr/+bbPSVV16pESNG2N+3b99uX4c5DuaCh+mKPm/evLNqg1WrVmnbtm02aD65N8Knn36qiy66SDVq1LDP/csvv2jFihX29Zu/179/fx04cOA3bVOsV69euvfeezVy5EjbLqZNzWstzazHxMTY1xwVFWXvX8zLy0tXXHGFJk+efFavCQCA30MgDwBAFZs/f7727t2rhQsX6pVXXrHB4MCBA22QvmzZMt155512SU5OtvfPzMzUJZdcYoNP85jFixeXBKJmfPjZeOyxx/Tyyy9r5cqV8vb21m233Wa3m4sEf/3rX9W6deuSHgVmW2FhoQ2S9+3bp6+++ko//fSTOnTooEsvvVSHDx8ueV4TTJvAedq0abY7/7Bhw3Tw4EF99913Jfc5cuSIzU7feOONdv348eM20DXB++rVq9WvXz8NGjRISUlJ5X495ng0b97cXmA42VNPPaXHH3/cBvvmtV5//fU2KDc9JRYtWmQvJDz55JNnfP73339fNWvWtO0yYcIEPf3005o7d669berUqXr11VftRYutW7dq5syZio+PL/P4zp07278FAEBF8q7QZwMAAL/LZHcnTpwoT09PtWjRwgaIJlh/9NFH7e2jR4/W888/b4u3XXfddZoyZYq97zvvvGMzzcakSZNsNt9k2vv27Vvuo/7cc8/p4osvtr+b8dsmSM/KyrJZa3NxwAS8JvNc+qLD+vXrlZqaarPOxksvvWSDVhPImky9YS4ofPjhh6pbt27JY82Fho8//tgG/cZnn31mX3vxetu2be1S7Nlnn9WMGTP0+eef65577inX6zHZd5MJPxVTn8BcHDDuu+8+G8ibHgnFY+lvv/323x0Tn5CQYC8IGLGxsXr99dftc/Tp08decDDH6rLLLpOPj4/NzJvAvTTTK8Hcz/RMMG0IAEBF4BsFAIAqZrLepYM607W8dCbXdMk23edN8GyYLLjJeAcGBtpg2ywmIDYBuMkqnw0TmBaLjIy0P4v/zqmYv20y52Z/iv+2WXbs2FHmb5tic6WDeMNk3k2GPjs7267/5z//sRcmzOszMjIybIY8Li7OXpQwz7t58+azysifOHFC/v7+v/tazTE2Sh9ns+1Mr/3k5yg+ZsWPGT58uP37ZjjBn/70J3sRwgyRKM1cIDFBfPExAACgIpCRBwCgipnsbWkmy36qbcXjy83PCy64wAbCJzs5eD6bv12c3T9TVXVzmwleTeb/ZKULw5nu5ycz3eTN42fPnm3HqJsu5mYoQTFTE8B0tTcZfjN9nAl6TZf8sxkuEBYWZnsMlPe1nrzt9yrKn6ldoqOjbY0B09XeDA+466679OKLL2rBggUljzPDD0w9AvPaAACoKATyAAA4OTMm/ZNPPikpNldZfH19lZ+f/5u/bcbHmy73Zk73s2GC1yFDhtgLEKZHgRnLbi5IFDOBvSkwd/XVV9t1k/k/mwJ+Rvv27fWPf/zDjuUvDtarknmNpoCfWe6++261bNnSXlgwx83YsGFDye8AAFQUutYDAODkTBd1k3k2Fd5N8Gu6tZusrxn3vXv37gr7OyZQN89titWZQnWmO7gZ/921a1dbzd1kz02gvWTJEltEzhTMK8++m4y8qWh/0003lbnNZOGnT59u/56pgn/DDTec9Zzrpgig6aK/ceNGVTUzvv7dd9+1wXpiYqKtEWACezPMoJhpr7OpYQAAQHkQyAMA4ORM12xTnd0UUzMZ7latWtlq82Z8dkVm6IcOHWoL1Jng2HTZN9OmmSy3qVbfs2dP+zdNVt2MczcBffG48zMxU9WZ8fymC7oJ1EszFd9NpX4zDZ7phm8K051t9tqM3S/O+lc1M7TgX//6ly2eZ8bSmyJ4X3zxhd2n4ikDzUWPW2+9tcr3DQDg3jwKTV80AAAAF2W6spueA8UFAZ2FqQGQlpamt99+29G7AgBwM2TkAQCASzOV6M0Ufmc7vr6ymZoGzzzzjKN3AwDghsjIAwAAAADgQsjIAwAAAADgQgjkAQAAAABwIQTyAAAAAAC4EAJ5AAAAAABcCIE8AAAAAAAuhEAeAAAAAAAXQiAPAAAAAIALIZAHAAAAAMCFEMgDAAAAACDX8f+QtkgS7g88cAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# plot the data\n", + "fig = plt.figure(figsize=(12,3))\n", + "ax = fig.add_subplot()\n", + "_ = ax.plot(arrivals[\"t\"], arrivals[\"mean_iat\"])\n", + "_ = ax.xaxis.set_major_locator(ticker.MultipleLocator(60))\n", + "_ = ax.set_xlabel(\"Time Interval (mins)\")\n", + "_ = ax.set_ylabel(\"Mean IAT (mins)\")" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# This example answer shows how to build up your functionality gradually.\n", + "\n", + "# iteration 1: ignoring time dependence. \n", + "# I just use the first mean in the dataframe and build up the basic structure of \n", + "# the code. A few lines of code and have an incorrect(!) but running model.\n", + "\n", + "def arrivals_generator(env, means, random_seed=None):\n", + " rng = np.random.default_rng(random_seed)\n", + " \n", + " for n in itertools.count():\n", + " interarrival_time = rng.exponential(means['mean_iat'].iloc[0])\n", + " yield env.timeout(interarrival_time)\n", + " print(f'arrival {n} at {env.now}')" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "arrival 0 at 36.06312905948992\n", + "arrival 1 at 71.10597389685671\n", + "arrival 2 at 106.87738889497054\n", + "arrival 3 at 111.07430324274428\n", + "arrival 4 at 112.37086423821994\n", + "arrival 5 at 134.1607719738122\n", + "arrival 6 at 155.31018238767578\n", + "arrival 7 at 202.17462173591824\n", + "arrival 8 at 203.36403469367087\n", + "arrival 9 at 219.06244739181835\n", + "arrival 10 at 220.11899199378547\n", + "arrival 11 at 236.45434641414488\n", + "arrival 12 at 262.424256707402\n", + "arrival 13 at 268.22767907625064\n", + "arrival 14 at 286.7014640161742\n", + "arrival 15 at 289.0080628235676\n", + "arrival 16 at 290.3817218344205\n", + "arrival 17 at 295.1094098664064\n", + "arrival 18 at 308.6273989969588\n", + "arrival 19 at 314.82217894303614\n", + "arrival 20 at 333.5329681902405\n", + "arrival 21 at 336.8866137734253\n", + "arrival 22 at 364.45616773559226\n", + "arrival 23 at 382.86245568380616\n", + "arrival 24 at 392.6993979932512\n", + "arrival 25 at 398.95569858792845\n", + "arrival 26 at 405.7554939274024\n", + "arrival 27 at 406.9116511332372\n", + "arrival 28 at 409.6061349488555\n", + "arrival 29 at 419.88594172291926\n", + "arrival 30 at 425.7161452653166\n", + "arrival 31 at 444.6792481287932\n", + "arrival 32 at 455.30660754003503\n", + "arrival 33 at 458.87555685450525\n", + "arrival 34 at 465.7918818983673\n", + "arrival 35 at 475.415232956225\n", + "arrival 36 at 480.5727333777919\n", + "arrival 37 at 485.4014016913866\n", + "arrival 38 at 498.58512192427474\n", + "arrival 39 at 503.03933005820676\n", + "arrival 40 at 523.044866433685\n" + ] + } + ], + "source": [ + "RUN_LENGTH = 540\n", + "env = simpy.Environment()\n", + "env.process(arrivals_generator(env, arrivals, random_seed=42))\n", + "env.run(RUN_LENGTH)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "# iteration 2. I've now added an index t used to select a mean IAT.\n", + "\n", + "def arrivals_generator(env, means, random_seed=None):\n", + " rng = np.random.default_rng(random_seed)\n", + " \n", + " for n in itertools.count():\n", + " \n", + " # this give us the index of dataframe to use\n", + " # I've used mod 9 so that run_length can be > 540\n", + " t = int(env.now // 60) % 9\n", + " interarrival_time = rng.exponential(means['mean_iat'].iloc[t])\n", + " yield env.timeout(interarrival_time)\n", + " print(f'arrival {n} at {env.now}')" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "arrival 0 at 36.06312905948992\n", + "arrival 1 at 71.10597389685671\n", + "arrival 2 at 99.72310589534777\n", + "arrival 3 at 103.08063737356677\n", + "arrival 4 at 104.1178861699473\n", + "arrival 5 at 121.54981235842111\n", + "arrival 6 at 131.41953721822412\n", + "arrival 7 at 153.2896089140706\n", + "arrival 8 at 153.8446682943552\n", + "arrival 9 at 161.17059422015734\n", + "arrival 10 at 161.66364836774198\n", + "arrival 11 at 169.2868137639097\n", + "arrival 12 at 181.40610523409634\n", + "arrival 13 at 183.3405793570459\n", + "arrival 14 at 189.49850767035377\n", + "arrival 15 at 190.2673739394849\n", + "arrival 16 at 190.72526027643588\n", + "arrival 17 at 192.30115628709785\n", + "arrival 18 at 196.80715266394864\n", + "arrival 19 at 198.8720793126411\n", + "arrival 20 at 205.1090090617092\n", + "arrival 21 at 206.2268909227708\n", + "arrival 22 at 215.4167422434931\n", + "arrival 23 at 221.55217155956439\n", + "arrival 24 at 224.8311523293794\n", + "arrival 25 at 226.91658586093848\n", + "arrival 26 at 229.1831843074298\n", + "arrival 27 at 229.56857004270807\n", + "arrival 28 at 230.46673131458084\n", + "arrival 29 at 233.8933335726021\n", + "arrival 30 at 235.8367347534012\n", + "arrival 31 at 242.15776904122674\n", + "arrival 32 at 247.82569406055575\n", + "arrival 33 at 249.7291336949399\n", + "arrival 34 at 253.41784038499966\n", + "arrival 35 at 258.5502942825238\n", + "arrival 36 at 261.30096117402616\n", + "arrival 37 at 263.8762509412766\n", + "arrival 38 at 270.907568398817\n", + "arrival 39 at 273.28314607024737\n", + "arrival 40 at 283.9527654705024\n", + "arrival 41 at 295.07967736046965\n", + "arrival 42 at 303.7523407841912\n", + "arrival 43 at 304.48262426399555\n", + "arrival 44 at 315.8232992812231\n", + "arrival 45 at 329.3666791598748\n", + "arrival 46 at 340.5870756744491\n", + "arrival 47 at 343.38756687681337\n", + "arrival 48 at 346.5969287871737\n", + "arrival 49 at 348.25428786018966\n", + "arrival 50 at 351.7979229871043\n", + "arrival 51 at 352.00932055120217\n", + "arrival 52 at 353.6828845123593\n", + "arrival 53 at 366.83345295783084\n", + "arrival 54 at 426.76979702560345\n", + "arrival 55 at 438.04055887337614\n", + "arrival 56 at 443.7213529053927\n", + "arrival 57 at 451.6764779835374\n", + "arrival 58 at 458.90332244223566\n", + "arrival 59 at 461.4959945067077\n", + "arrival 60 at 483.16249720413987\n", + "arrival 61 at 510.6794742744887\n" + ] + } + ], + "source": [ + "RUN_LENGTH = 540\n", + "env = simpy.Environment()\n", + "env.process(arrivals_generator(env, arrivals, random_seed=42))\n", + "env.run(RUN_LENGTH)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exercise 2: Thinning the arrivals\n", + "\n", + "**Task:**\n", + "* Update your exercise 1 code to include an implementation of thinning\n", + "* What do you notice about the total number of arrivals compared to the previous example? Why has the changed occurred?\n", + " * If you are not controlling your sampling with random seeds you will need to run each implementation a few times.\n", + "\n", + "**Hints:**\n", + "* You will need a second distribution - Uniform(0, 1) to do the thinning. If you are controlling random sampling through seeds that means you will need a second seed.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "# your code here ..." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "# Example answer ...\n", + "\n", + "# I've added an extra bit of code here to report the number of rejections for \n", + "# each arrival. This can be useful for debugging.\n", + "\n", + "def arrivals_generator_with_thinning(\n", + " env: simpy.Environment, \n", + " means: pd.DataFrame, \n", + " audit: list | None = None, \n", + " seed1: int | np.random.SeedSequence | None = None, \n", + " seed2: int | np.random.SeedSequence | None = None\n", + "):\n", + " \"\"\"\n", + " Simulate interā€‘arrival times from a NSPP using thinning.\n", + " (acceptance/rejection)\n", + "\n", + " Parameters\n", + " ----------\n", + " env : simpy.Environment\n", + " Used to access the simulation clock with env.now\n", + " means : pandas.DataFrame\n", + " Must have an 'arrival_rate' column (or be the rate series)\n", + " audit : list, optional\n", + " For logging details of accepted arrivals\n", + " seed1, seed2 : int | SeedSequence, optional (default = None)\n", + " Seeds for the exponential RNG (interā€‘arrival) and uniform RNG (thinning).\n", + " \"\"\"\n", + "\n", + " # exponential sampling for next IAT\n", + " exp_rng = np.random.default_rng(seed1)\n", + "\n", + " # uniform sampling - for accept/reject\n", + " unif_rng = np.random.default_rng(seed2)\n", + " \n", + " # maximum arrival rate (smallest time between arrivals)\n", + " lambda_max = means['arrival_rate'].max()\n", + " \n", + " # Main thinning loop\n", + " for n_arrivals in itertools.count():\n", + "\n", + " interarrival_time = 0.0 \n", + " rejects = 0\n", + " \n", + " while True:\n", + " \n", + " # 1. Sample exponential interā€‘arrival time using. Ī»_max\n", + " # note we use += here as we may reject an arrival \n", + " interarrival_time += exp_rng.exponential(1 / lambda_max)\n", + "\n", + " # 2. Advance simulation clock to the candidate arrival time\n", + " cand_time = env.now + interarrival_time\n", + "\n", + " # 3. Look up the rate at that clock time of candidate arrival\n", + " # note that the 60 (interval width) and 9 (no. intervals) are specific to this problem.\n", + " t_idx = int(cand_time // 60) % 9\n", + " lambda_t = means['arrival_rate'].iloc[t_idx]\n", + "\n", + " # 4. Accept with probability Ī»(t) / Ī»_max\n", + " u = unif_rng.uniform(0, 1)\n", + " if u <= lambda_t / lambda_max:\n", + " break # accepted ā†’ exit inner while\n", + " # else: reject and draw another candidate\n", + " rejects += 1\n", + "\n", + " # yield the accepted interā€‘arrival time\n", + " yield env.timeout(interarrival_time)\n", + "\n", + " # debugging\n", + " if audit is not None:\n", + " audit[-1] += 1\n", + " else:\n", + " # Show how many candidates were rejected before this acceptance\n", + " print(\n", + " f'arrival {n_arrivals} at {env.now:.2f}. '\n", + " f'Rejected samples = {rejects}'\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "arrival 0 at 37.46. Rejected samples = 4\n", + "arrival 1 at 67.39. Rejected samples = 2\n", + "arrival 2 at 73.02. Rejected samples = 1\n", + "arrival 3 at 89.41. Rejected samples = 3\n", + "arrival 4 at 96.79. Rejected samples = 2\n", + "arrival 5 at 98.37. Rejected samples = 0\n", + "arrival 6 at 104.94. Rejected samples = 1\n", + "arrival 7 at 112.30. Rejected samples = 1\n", + "arrival 8 at 121.49. Rejected samples = 0\n", + "arrival 9 at 127.62. Rejected samples = 0\n", + "arrival 10 at 135.25. Rejected samples = 2\n", + "arrival 11 at 135.64. Rejected samples = 0\n", + "arrival 12 at 141.91. Rejected samples = 2\n", + "arrival 13 at 148.23. Rejected samples = 0\n", + "arrival 14 at 155.26. Rejected samples = 2\n", + "arrival 15 at 158.47. Rejected samples = 0\n", + "arrival 16 at 160.19. Rejected samples = 0\n", + "arrival 17 at 167.68. Rejected samples = 2\n", + "arrival 18 at 174.35. Rejected samples = 0\n", + "arrival 19 at 181.30. Rejected samples = 0\n", + "arrival 20 at 186.72. Rejected samples = 0\n", + "arrival 21 at 187.09. Rejected samples = 0\n", + "arrival 22 at 192.76. Rejected samples = 0\n", + "arrival 23 at 199.53. Rejected samples = 0\n", + "arrival 24 at 205.14. Rejected samples = 0\n", + "arrival 25 at 206.54. Rejected samples = 0\n", + "arrival 26 at 208.15. Rejected samples = 0\n", + "arrival 27 at 208.97. Rejected samples = 0\n", + "arrival 28 at 210.75. Rejected samples = 0\n", + "arrival 29 at 210.85. Rejected samples = 0\n", + "arrival 30 at 211.69. Rejected samples = 0\n", + "arrival 31 at 218.26. Rejected samples = 0\n", + "arrival 32 at 238.24. Rejected samples = 0\n", + "arrival 33 at 244.47. Rejected samples = 2\n", + "arrival 34 at 246.92. Rejected samples = 1\n", + "arrival 35 at 259.22. Rejected samples = 1\n", + "arrival 36 at 266.76. Rejected samples = 0\n", + "arrival 37 at 293.31. Rejected samples = 0\n", + "arrival 38 at 320.47. Rejected samples = 4\n", + "arrival 39 at 328.36. Rejected samples = 1\n", + "arrival 40 at 346.04. Rejected samples = 1\n", + "arrival 41 at 347.95. Rejected samples = 1\n", + "arrival 42 at 380.10. Rejected samples = 3\n", + "arrival 43 at 393.26. Rejected samples = 4\n", + "arrival 44 at 402.71. Rejected samples = 2\n", + "arrival 45 at 406.88. Rejected samples = 1\n", + "arrival 46 at 428.54. Rejected samples = 3\n", + "arrival 47 at 435.37. Rejected samples = 3\n", + "arrival 48 at 435.75. Rejected samples = 0\n", + "arrival 49 at 485.80. Rejected samples = 5\n", + "arrival 50 at 505.93. Rejected samples = 3\n", + "arrival 51 at 517.75. Rejected samples = 1\n" + ] + } + ], + "source": [ + "RUN_LENGTH = 540\n", + "\n", + "env = simpy.Environment()\n", + "env.process(arrivals_generator_with_thinning(env, arrivals, \n", + " seed1=42, seed2=101))\n", + "env.run(RUN_LENGTH)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exercise 3: Validate the total number of arrivals in 540 minutes.\n", + "\n", + "Here we will repeat the same 10,000 times and then explore the distribution of the number of arrivals. If all has gone to plan this should be a Poisson distribution with mean ~53." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "# %%timeit # <- uncomment this line to run the code multiple times and evaluate efficiency\n", + "\n", + "RUN_LENGTH = 540\n", + "REPLICATIONS = 10_000\n", + "audit = []\n", + "\n", + "for i in range(REPLICATIONS):\n", + " # set up audit for replication.\n", + " audit.append(0)\n", + " env = simpy.Environment()\n", + "\n", + " seed_sequence = np.random.SeedSequence(i)\n", + " seeds = seed_sequence.spawn(2)\n", + " \n", + " env.process(arrivals_generator_with_thinning(\n", + " env, arrivals, audit, seeds[0], seeds[1]))\n", + " env.run(RUN_LENGTH)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjEAAAGdCAYAAADjWSL8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAjf0lEQVR4nO3dC7BV5Xk/4Bc5gIBAuQhIRSGJEhS0GUy5xBEUBKlIjU4h0lAciZpqIFSogulMaMeI2omYDlNLrCMRsTidiLGFGLAqlkEu0qECQasjtFC5GMM9BBD3f77V7P3nACIH0cN3eJ6ZJXut9Z591vnYh/3zu6xdr1QqlQIAIDNn1PYFAACcCCEGAMiSEAMAZEmIAQCyJMQAAFkSYgCALAkxAECWhBgAIEtVUUd99NFH8d5770WzZs2iXr16tX05AMBxSPfg3bVrV3To0CHOOOOM0zPEpADTsWPH2r4MAOAEbNiwIc4999zTM8SkHphyIzRv3ry2LwcAOA47d+4sOiHK7+OnZYgpDyGlACPEAEBejmcqiIm9AECWhBgAIEtCDACQJSEGAMiSEAMAZEmIAQCyJMQAAFkSYgCALAkxAECWhBgAIEtCDACQJSEGAMiSEAMAZEmIAQCyVFXbFwDUvk4T50Zu1j9wbW1fAlDL9MQAAFkSYgCALAkxAECWhBgAIEtCDACQJSEGAMiSEAMAZEmIAQCyJMQAAFkSYgCALAkxAECWhBgAIEtCDACQJSEGAMiSEAMAZEmIAQCyJMQAAHU/xDz66KNxySWXRPPmzYutd+/e8fOf/7xyvlQqxeTJk6NDhw7RuHHj6NevX6xZs6bac+zbty/GjBkTbdq0iaZNm8bQoUNj48aN1Wq2bdsWI0eOjBYtWhRberx9+/ZP+7MCAKdriDn33HPjgQceiNdff73YrrrqqvjjP/7jSlB56KGH4uGHH45p06bF8uXLo3379nH11VfHrl27Ks8xbty4mDNnTsyePTsWLVoUu3fvjiFDhsTBgwcrNSNGjIiVK1fGCy+8UGzpcQoyAABl9Uqp++RTaNWqVfzt3/5t3HLLLUUPTAop99xzT6XXpV27dvHggw/G7bffHjt27Iizzz47Zs6cGcOHDy9q3nvvvejYsWPMmzcvBg0aFGvXro2LLroolixZEj179ixq0uPU6/Pmm29Gly5djuu6du7cWfTipO+Zeo2Aj9dp4tzsmmf9A9fW9iUAn4GavH+f8JyY1HOSelP27NlTBIx169bF5s2bY+DAgZWaRo0aRd++fWPx4sXF/ooVK+LAgQPValLw6datW6XmtddeKy6+HGCSXr16FcfKNUeTAlP6wQ/dAIC6q8YhZtWqVXHWWWcVAeXb3/52MTSUek5SgElSz8uh0n75XPqzYcOG0bJly2PWtG3b9ojvm46Va45mypQplTk0aUu9OwBA3VXjEJOGc9IclTTE8+d//ucxatSo+OUvf1k5X69evWr1abTq8GOHO7zmaPWf9DyTJk0qup7K24YNG2r4kwEAdTrEpJ6UL33pS3HZZZcVvR+XXnpp/OhHPyom8SaH95Zs3bq10juTavbv31+sPjpWzZYtW474vu+///4RvTyHSj1D5VVT5Q0AqLs+9X1iUg9Jmo/SuXPnIoAsWLCgci4FloULF0afPn2K/R49ekSDBg2q1WzatClWr15dqUnza1JPyrJlyyo1S5cuLY6VawAAqmrSBPfee28MHjy4mG+Slk2nib2vvPJKsQw6DfWklUn3339/XHDBBcWWHjdp0qRYMp2kuSqjR4+O8ePHR+vWrYuVTRMmTIju3bvHgAEDipquXbvGNddcE7feemtMnz69OHbbbbcVy7CPd2USAFD31SjEpGGedL+W1HuSAkm68V0KMOleMMndd98de/fujTvuuKMYMkorjObPnx/NmjWrPMfUqVOjqqoqhg0bVtT2798/ZsyYEfXr16/UzJo1K8aOHVtZxZRuiJfuPQMAcNLuE3Oqcp8YOH7uEwOcVveJAQCoTUIMAJAlIQYAyJIQAwBkSYgBALIkxAAAWRJiAIAsCTEAQJaEGAAgS0IMAJAlIQYAyJIQAwBkSYgBALIkxAAAWRJiAIAsCTEAQJaEGAAgS0IMAJAlIQYAyJIQAwBkSYgBALIkxAAAWRJiAIAsCTEAQJaEGAAgS0IMAJAlIQYAyJIQAwBkSYgBALIkxAAAWRJiAIAsCTEAQJaEGAAgS0IMAJAlIQYAyJIQAwBkSYgBALIkxAAAWRJiAIAsCTEAQJaEGAAgS0IMAJAlIQYAyJIQAwDU/RAzZcqU+OpXvxrNmjWLtm3bxvXXXx9vvfVWtZqbb7456tWrV23r1atXtZp9+/bFmDFjok2bNtG0adMYOnRobNy4sVrNtm3bYuTIkdGiRYtiS4+3b9/+aX5WAOB0DTELFy6MO++8M5YsWRILFiyIDz/8MAYOHBh79uypVnfNNdfEpk2bKtu8efOqnR83blzMmTMnZs+eHYsWLYrdu3fHkCFD4uDBg5WaESNGxMqVK+OFF14otvQ4BRkAgKSqJs2QwsShnnjiiaJHZsWKFXHFFVdUjjdq1Cjat29/1OfYsWNHPP744zFz5swYMGBAceypp56Kjh07xosvvhiDBg2KtWvXFt8rhaWePXsWNY899lj07t276Pnp0qWLvz0AOM19qjkxKZAkrVq1qnb8lVdeKcLNhRdeGLfeemts3bq1ci4FngMHDhQ9OGUdOnSIbt26xeLFi4v91157rRhCKgeYJA1JpWPlmsOlIaqdO3dW2wCAuuuEQ0ypVIq77rorLr/88iKAlA0ePDhmzZoVL730Uvzwhz+M5cuXx1VXXVWEjGTz5s3RsGHDaNmyZbXna9euXXGuXJNC0OHSsXLN0ebrlOfPpC317AAAdVeNhpMO9Z3vfCfeeOONYk7LoYYPH155nMLNZZddFueff37MnTs3brjhhmOGojQJuOzQxx9Xc6hJkyYVoaos9cQIMgBQd51QT0xaWfT888/Hyy+/HOeee+4xa88555wixLz99tvFfpors3///mL10aHSkFPqjSnXbNmy5Yjnev/99ys1h0vzcJo3b15tAwDqrhqFmNQTknpgnn322WK4qHPnzp/4NR988EFs2LChCDNJjx49okGDBsXqprK0gmn16tXRp0+fYj9N4E3zbZYtW1apWbp0aXGsXAMAnN5qNJyUllc//fTT8bOf/ay4V0x5fkqag9K4ceNiqfTkyZPjxhtvLELL+vXr49577y3uB/P1r3+9Ujt69OgYP358tG7dupgUPGHChOjevXtltVLXrl2LZdppUvD06dOLY7fddluxDNvKJACgxiHm0UcfLf7s16/fEUut003u6tevH6tWrYonn3yyuDFdCjJXXnllPPPMM0XoKZs6dWpUVVXFsGHDYu/evdG/f/+YMWNG8fVlaXLw2LFjK6uY0g3xpk2b5m8NACjUK6UxojooTexNvT5pCMr8GDi2ThPnZtdE6x+4trYvAajl92+fnQQAZEmIAQCyJMQAAFkSYgCALAkxAECWhBgAIEtCDACQJSEGAMiSEAMAZEmIAQCyJMQAAFkSYgCALAkxAECWhBgAIEtCDACQJSEGAMiSEAMAZEmIAQCyJMQAAFkSYgCALAkxAECWhBgAIEtVtX0BUNd0mji3ti8B4LSgJwYAyJIQAwBkSYgBALIkxAAAWRJiAIAsCTEAQJaEGAAgS0IMAJAlIQYAyJIQAwBkSYgBALIkxAAAWRJiAIAsCTEAQJaEGAAgS0IMAJAlIQYAyJIQAwBkSYgBALIkxAAAWRJiAIC6H2KmTJkSX/3qV6NZs2bRtm3buP766+Ott96qVlMqlWLy5MnRoUOHaNy4cfTr1y/WrFlTrWbfvn0xZsyYaNOmTTRt2jSGDh0aGzdurFazbdu2GDlyZLRo0aLY0uPt27d/mp8VADhdQ8zChQvjzjvvjCVLlsSCBQviww8/jIEDB8aePXsqNQ899FA8/PDDMW3atFi+fHm0b98+rr766ti1a1elZty4cTFnzpyYPXt2LFq0KHbv3h1DhgyJgwcPVmpGjBgRK1eujBdeeKHY0uMUZAAAknql1HVygt5///2iRyaFmyuuuKLohUk9MCmk3HPPPZVel3bt2sWDDz4Yt99+e+zYsSPOPvvsmDlzZgwfPryoee+996Jjx44xb968GDRoUKxduzYuuuiiIiz17NmzqEmPe/fuHW+++WZ06dLlE69t586dRQ9O+n7Nmzf3t83nptPEuVr7c7D+gWu1M9RBNXn//lRzYtI3SFq1alX8uW7duti8eXPRO1PWqFGj6Nu3byxevLjYX7FiRRw4cKBaTQo+3bp1q9S89tprxQ9QDjBJr169imPlmsOlsJR+8EM3AKDuOuEQk3pd7rrrrrj88suLAJKkAJOknpdDpf3yufRnw4YNo2XLlsesST08h0vHyjVHm69Tnj+TttSzAwDUXSccYr7zne/EG2+8Ef/0T/90xLl69eodEXgOP3a4w2uOVn+s55k0aVLRM1TeNmzYUIOfBgA4LUJMWln0/PPPx8svvxznnntu5XiaxJsc3luydevWSu9Mqtm/f3+x+uhYNVu2bDnqHJzDe3kOHbZKY2eHbgBA3VWjEJN6QlIPzLPPPhsvvfRSdO7cudr5tJ8CSFq5VJYCS5r426dPn2K/R48e0aBBg2o1mzZtitWrV1dq0gTe1JuybNmySs3SpUuLY+UaAOD0VlWT4rS8+umnn46f/exnxb1iyj0uaQ5KuidMGupJK5Puv//+uOCCC4otPW7SpEmxZLpcO3r06Bg/fny0bt26mBQ8YcKE6N69ewwYMKCo6dq1a1xzzTVx6623xvTp04tjt912W7EM+3hWJgEAdV+NQsyjjz5a/JluYHeoJ554Im6++ebi8d133x179+6NO+64oxgySiuM5s+fX4SesqlTp0ZVVVUMGzasqO3fv3/MmDEj6tevX6mZNWtWjB07trKKKd0QL917BgDgU98n5lTmPjHUFveJ+Xy4TwzUTZ/bfWIAAGqLEAMAZEmIAQCyJMQAAHV/dRLAqSLHCdQmI8PJpScGAMiSEAMAZEmIAQCyJMQAAFkSYgCALAkxAECWhBgAIEtCDACQJSEGAMiSEAMAZEmIAQCyJMQAAFkSYgCALAkxAECWhBgAIEtCDACQJSEGAMiSEAMAZEmIAQCyJMQAAFkSYgCALAkxAECWhBgAIEtCDACQJSEGAMiSEAMAZEmIAQCyJMQAAFkSYgCALAkxAECWhBgAIEtCDACQJSEGAMiSEAMAZEmIAQCyJMQAAFkSYgCALAkxAECWhBgA4PQIMa+++mpcd9110aFDh6hXr14899xz1c7ffPPNxfFDt169elWr2bdvX4wZMybatGkTTZs2jaFDh8bGjRur1Wzbti1GjhwZLVq0KLb0ePv27Sf6cwIAp3uI2bNnT1x66aUxbdq0j6255pprYtOmTZVt3rx51c6PGzcu5syZE7Nnz45FixbF7t27Y8iQIXHw4MFKzYgRI2LlypXxwgsvFFt6nIIMAEBSVdNmGDx4cLEdS6NGjaJ9+/ZHPbdjx454/PHHY+bMmTFgwIDi2FNPPRUdO3aMF198MQYNGhRr164tgsuSJUuiZ8+eRc1jjz0WvXv3jrfeeiu6dOnibw8ATnOfyZyYV155Jdq2bRsXXnhh3HrrrbF169bKuRUrVsSBAwdi4MCBlWNpaKpbt26xePHiYv+1114rhpDKASZJQ1LpWLnmcGmIaufOndU2AKDuOukhJvXSzJo1K1566aX44Q9/GMuXL4+rrrqqCBnJ5s2bo2HDhtGyZctqX9euXbviXLkmhaDDpWPlmsNNmTKlMn8mbalnBwCou2o8nPRJhg8fXnmcelcuu+yyOP/882Pu3Llxww03fOzXlUqlYhJw2aGPP67mUJMmTYq77rqrsp96YgQZAKi7PvMl1uecc04RYt5+++1iP82V2b9/f7H66FBpyCn1xpRrtmzZcsRzvf/++5Wao83Dad68ebUNAKi7PvMQ88EHH8SGDRuKMJP06NEjGjRoEAsWLKjUpBVMq1evjj59+hT7aQJvmgC8bNmySs3SpUuLY+UaAOD0VuPhpLQc+p133qnsr1u3rlj+3KpVq2KbPHly3HjjjUVoWb9+fdx7773F/WC+/vWvF/Vpvsro0aNj/Pjx0bp16+JrJkyYEN27d6+sVuratWuxTDtNCp4+fXpx7LbbbiuWYVuZBACcUIh5/fXX48orr6zsl+ehjBo1Kh599NFYtWpVPPnkk8WN6VKQSbXPPPNMNGvWrPI1U6dOjaqqqhg2bFjs3bs3+vfvHzNmzIj69etXatLk4LFjx1ZWMaUb4h3r3jQAwOmlXinNlq2D0sTe1OuThqDMj+Hz1GniXA3OUa1/4FotAyfx/dtnJwEAWRJiAIAsCTEAQJaEGAAgS0IMAJAlIQYAyJIQAwBkSYgBALIkxAAAWRJiAIAsCTEAQJaEGAAgS0IMAJAlIQYAyJIQAwBkSYgBALIkxAAAWRJiAIAsCTEAQJaEGAAgS0IMAJAlIQYAyJIQAwBkSYgBALIkxAAAWRJiAIAsCTEAQJaEGAAgS0IMAJAlIQYAyFJVbV8AHEuniXM1EABHpScGAMiSEAMAZEmIAQCyJMQAAFkSYgCALAkxAECWhBgAIEtCDACQJSEGAMiSEAMAZEmIAQCyJMQAAFkSYgCA0yPEvPrqq3HddddFhw4dol69evHcc89VO18qlWLy5MnF+caNG0e/fv1izZo11Wr27dsXY8aMiTZt2kTTpk1j6NChsXHjxmo127Zti5EjR0aLFi2KLT3evn37if6cAMDpHmL27NkTl156aUybNu2o5x966KF4+OGHi/PLly+P9u3bx9VXXx27du2q1IwbNy7mzJkTs2fPjkWLFsXu3btjyJAhcfDgwUrNiBEjYuXKlfHCCy8UW3qcggwAQFKvlLpOTlDqiUlh5Prrry/201OlHpgUUu65555Kr0u7du3iwQcfjNtvvz127NgRZ599dsycOTOGDx9e1Lz33nvRsWPHmDdvXgwaNCjWrl0bF110USxZsiR69uxZ1KTHvXv3jjfffDO6dOnyide2c+fOogcnfb/mzZv7285Up4lza/sS4KRZ/8C1WhNO4vv3SZ0Ts27duti8eXMMHDiwcqxRo0bRt2/fWLx4cbG/YsWKOHDgQLWaFHy6detWqXnttdeKH6AcYJJevXoVx8o1h0thKf3gh24AQN11UkNMCjBJ6nk5VNovn0t/NmzYMFq2bHnMmrZt2x7x/OlYueZwU6ZMqcyfSVvq2QEA6q7PZHVSGmY6VBpmOvzY4Q6vOVr9sZ5n0qRJRddTeduwYcMJXz8AcJqFmDSJNzm8t2Tr1q2V3plUs3///mL10bFqtmzZcsTzv//++0f08hw6bJXGzg7dAIC6q+pkPlnnzp2LALJgwYL4yle+UhxLgWXhwoXFxN6kR48e0aBBg6Jm2LBhxbFNmzbF6tWri5VNSZrAm3pTli1bFn/4h39YHFu6dGlxrE+fPifzkgE+NzlOVDcZmToVYtJy6HfeeafaZN60/LlVq1Zx3nnnFSuT7r///rjggguKLT1u0qRJsWQ6SfNVRo8eHePHj4/WrVsXXzdhwoTo3r17DBgwoKjp2rVrXHPNNXHrrbfG9OnTi2O33XZbsQz7eFYmAQB1X41DzOuvvx5XXnllZf+uu+4q/hw1alTMmDEj7r777ti7d2/ccccdxZBRWmE0f/78aNasWeVrpk6dGlVVVUVPTKrt379/8bX169ev1MyaNSvGjh1bWcWUboj3cfemAQBOP5/qPjGnMveJqRty7H6HusRwEqfNfWIAAD4vQgwAkCUhBgDIkhADAGRJiAEAsiTEAABZEmIAgCwJMQBAloQYACBLQgwAkCUhBgDIkhADAGRJiAEAsiTEAABZEmIAgCwJMQBAloQYACBLQgwAkCUhBgDIkhADAGRJiAEAsiTEAABZEmIAgCwJMQBAloQYACBLQgwAkCUhBgDIkhADAGRJiAEAsiTEAABZEmIAgCwJMQBAloQYACBLQgwAkCUhBgDIkhADAGRJiAEAsiTEAABZEmIAgCwJMQBAloQYACBLQgwAkCUhBgDIkhADAGTppIeYyZMnR7169apt7du3r5wvlUpFTYcOHaJx48bRr1+/WLNmTbXn2LdvX4wZMybatGkTTZs2jaFDh8bGjRtP9qUCABn7THpiLr744ti0aVNlW7VqVeXcQw89FA8//HBMmzYtli9fXgScq6++Onbt2lWpGTduXMyZMydmz54dixYtit27d8eQIUPi4MGDn8XlAgAZqvpMnrSqqlrvy6G9MI888kh873vfixtuuKE49pOf/CTatWsXTz/9dNx+++2xY8eOePzxx2PmzJkxYMCAouapp56Kjh07xosvvhiDBg36LC4ZAMjMZ9IT8/bbbxfDRZ07d45vfOMb8e677xbH161bF5s3b46BAwdWahs1ahR9+/aNxYsXF/srVqyIAwcOVKtJz9WtW7dKzdGkIaidO3dW2wCAuuukh5iePXvGk08+Gb/4xS/iscceK0JLnz594oMPPigeJ6nn5VBpv3wu/dmwYcNo2bLlx9YczZQpU6JFixaVLfXcAAB110kPMYMHD44bb7wxunfvXgwHzZ07tzJsVJYm+x4+zHT4scN9Us2kSZOKoajytmHDhk/9swAAp/ES67S6KAWaNMRUnidzeI/K1q1bK70zqWb//v2xbdu2j605mjQs1bx582obAFB3feYhJs1VWbt2bZxzzjnFHJkUUhYsWFA5nwLLwoULiyGnpEePHtGgQYNqNWmF0+rVqys1AAAnfXXShAkT4rrrrovzzjuv6D257777ikm2o0aNKoaD0vLp+++/Py644IJiS4+bNGkSI0aMKL4+zWcZPXp0jB8/Plq3bh2tWrUqnrM8PAUA8JmEmHRTuptuuil+9atfxdlnnx29evWKJUuWxPnnn1+cv/vuu2Pv3r1xxx13FENGaSLw/Pnzo1mzZpXnmDp1arFMe9iwYUVt//79Y8aMGVG/fn1/awBAoV4pzZitg1LvT+rVSZN8zY/JV6eJ/zcxHKgd6x+4VtNzyr5/++wkACBLQgwAkCUhBgDIkhADAGRJiAEAsiTEAABZEmIAgCwJMQBAloQYACBLQgwAkCUhBgDIkhADAGRJiAEAsiTEAABZEmIAgCwJMQBAloQYACBLQgwAkCUhBgDIkhADAGRJiAEAsiTEAABZEmIAgCxV1fYFAHDq6jRxbuRm/QPX1vYl8DkRYk4jOf5jBAAfx3ASAJAlIQYAyJIQAwBkSYgBALIkxAAAWRJiAIAsCTEAQJaEGAAgS0IMAJAlIQYAyJIQAwBkSYgBALIkxAAAWRJiAIAsCTEAQJaEGAAgS0IMAJAlIQYAyFJVbV8AAJxMnSbOzbJB1z9wbW1fQnZO+Z6Yv//7v4/OnTvHmWeeGT169Ih///d/r+1LAgBOAad0T8wzzzwT48aNK4LM1772tZg+fXoMHjw4fvnLX8Z5551Xq9eWa9IHgLrilO6Jefjhh2P06NHxrW99K7p27RqPPPJIdOzYMR599NHavjQAoJadsj0x+/fvjxUrVsTEiROrHR84cGAsXrz4iPp9+/YVW9mOHTuKP3fu3PmZXN9H+37zmTwvAKenz+r9Ktd2KJVK+YaYX/3qV3Hw4MFo165dteNpf/PmzUfUT5kyJf76r//6iOOp5wYATnUtHqntKzi17Nq1K1q0aJFniCmrV69etf2UzA4/lkyaNCnuuuuuyv5HH30Uv/71r6N169ZHrf+sU2QKTxs2bIjmzZt/rt+7LtGO2vFU4vWoHU8ldfn1WCqVigDToUOHT6w9ZUNMmzZton79+kf0umzduvWI3pmkUaNGxXao3/u934valF5Yde3FVRu0o3Y8lXg9asdTSfM6+j7zST0wp/zE3oYNGxZLqhcsWFDteNrv06dPrV0XAHBqOGV7YpI0PDRy5Mi47LLLonfv3vHjH/84/ud//ie+/e1v1/alAQC17JQOMcOHD48PPvgg/uZv/iY2bdoU3bp1i3nz5sX5558fp7I0rPX973//iOEttGNt8HrUjqcSr0fteDLVKx3PGiYAgFPMKTsnBgDgWIQYACBLQgwAkCUhBgDIkhBzgtKHUF5yySWVGw2lJeA///nPK+fTfOnJkycXdxxs3Lhx9OvXL9asWXOy/t7qrPTxEekOy+nTy8u05SdLr7XUbodu7du314Yn4H//93/jm9/8ZnG37yZNmsQf/MEfFJ/j5vV4/Dp16nTE6zFtd955p9/pGvjwww/jr/7qr6Jz587F+8gXvvCFYrVuuiO91+PvpNVJ1Nzzzz9fmjt3bumtt94qtnvvvbfUoEGD0urVq4vzDzzwQKlZs2aln/70p6VVq1aVhg8fXjrnnHNKO3fu1NwfY9myZaVOnTqVLrnkktJ3v/vdynFt+cm+//3vly6++OLSpk2bKtvWrVu1YQ39+te/Lp1//vmlm2++ubR06dLSunXrSi+++GLpnXfe0ZY1kF57h74WFyxYkFbBll5++WW/0zVw3333lVq3bl3613/91+K1+M///M+ls846q/TII494Pf6OEHMStWzZsvSP//iPpY8++qjUvn374s237Le//W2pRYsWpX/4h384md+yzti1a1fpggsuKP6x69u3byXEaMvjDzGXXnrpUc9pw+N3zz33lC6//PKPPa8tT0z6ff7iF79YtJ82PH7XXntt6ZZbbql27IYbbih985vf9Hr8HcNJJ0H6tO3Zs2fHnj17imGldevWFZ/5NHDgwGo3eOrbt28sXrz4ZHzLOid1M1977bUxYMCAase15fF7++23i+HL1PX8jW98I959911tWEPPP/98cYfwP/mTP4m2bdvGV77ylXjssce8Hj+F/fv3x1NPPRW33HJLMaTkd/r4XX755fFv//Zv8V//9V/F/n/+53/GokWL4o/+6I+K/XXea07tO/ae6latWlWElt/+9rdx1llnxZw5c+Kiiy6qBJXDP6gy7f/3f/93LV3tqSsFwP/4j/+I5cuXH3Gu/AGg2vLYevbsGU8++WRceOGFsWXLlrjvvvuKzxhL87C04fFLwS/Nd0sfeXLvvffGsmXLYuzYscX/hPzZn/2ZtjwBzz33XGzfvj1uvvlmv9M1dM8998SOHTviy1/+cvGByOl/mH/wgx/ETTfdpC1/R4j5FLp06RIrV64sfkF/+tOfxqhRo2LhwoWV8+n/Og6Vhu8OP3a6Sx8j/93vfjfmz58fZ5555sfWactjGzx4cOVx9+7di3D9xS9+MX7yk59Er169tOFxShMmU0/M/fffX+ynnpgUBFOwSSHG67HmHn/88eL1mXoJD+V3+pM988wzRS/W008/HRdffHHxfpMWPaS2TO83ZadzWxpO+pSftP2lL32p+Ecvraq59NJL40c/+lFlVUj5/4DLtm7dekSPwukurfpI7ZI+sbyqqqrYUhD8u7/7u+Jxub20Zc00bdq0CDNpiMnr8fidc845RW/qobp27Vp88GyiLWsm9Ty/+OKL8a1vfatyTBsev7/8y7+MiRMnFsPD6fc5fSDyX/zFXxTvN9ry/wgxJ1FKv/v27SvmJKRf1AULFlQbF05vzqmLn/+vf//+xbBc+j+M8pZC4Z/+6Z8Wj9OSQm1Zc+l1uHbt2uJN2evx+H3ta1+Lt956q9qxNB+h/KGz2rJmnnjiiWJuUZrvVqYNj99vfvObOOOM6m/TaVipvMS6s/caS6xP1KRJk0qvvvpqseztjTfeKJZYn3HGGaX58+cX59PKpLQa6dlnny2WWN90002WWB+nQ1cnacvjM378+NIrr7xSevfdd0tLliwpDRkypFjiv379em1Yw2X+VVVVpR/84Aelt99+uzRr1qxSkyZNSk899ZTXYw0dPHiwdN555xUrvg7n38fjM2rUqNLv//7vV5ZYp/eTNm3alO6++25t+TuWWJ+gtOwt3U+iYcOGpbPPPrvUv3//SoBJ0jLCtOw1LbVu1KhR6YorrijCDDUPMdryk5XvQ5TuVdShQ4diGeaaNWu04Qn4l3/5l1K3bt2K39svf/nLpR//+MfVzns9Hp9f/OIXxb1h0n20DqcNj0+6r1j6tzCFwTPPPLP0hS98ofS9732vtG/fPm35O/XSf2rYQwgAUOvMiQEAsiTEAABZEmIAgCwJMQBAloQYACBLQgwAkCUhBgDIkhADAGRJiAEAsiTEAABZEmIAgCwJMQBA5Oj/AU+tNR3tPCIlAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# distribution\n", + "plt.hist(audit);" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "53.15\n" + ] + } + ], + "source": [ + "# mean\n", + "print(np.array(audit).mean().round(2))" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "53.07" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# expected arrivals from data.\n", + "round(sum(arrivals['arrival_rate'] * 60), 2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Optional Exercise 4: Pre-calculate the acceptance probabilities\n", + "\n", + "At the moment the code calculates the probability of accepting an arrival on each iteration. We can make the code more efficient by pre-calculating the acceptance probabilities before we run the function." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "def arrivals_generator_with_thinning2(\n", + " env: simpy.Environment, \n", + " means: pd.DataFrame, \n", + " audit: list | None = None, \n", + " seed1: int | np.random.SeedSequence | None = None, \n", + " seed2: int | np.random.SeedSequence | None = None\n", + "):\n", + " \"\"\"\n", + " Simulate interā€‘arrival times from a NHPP using thinning.\n", + " This version pre-calculates the probabilities.\n", + "\n", + " Parameters\n", + " ----------\n", + " env : simpy.Environment\n", + " Used to access the simulation clock with env.now\n", + " means : pd.DataFrame\n", + " Must have an 'arrival_rate' column (or be the rate series)\n", + " audit : list, optional (default = None)\n", + " For logging details of accepted arrivals\n", + " seed1, seed2 : int | SeedSequence, optional (default = None)\n", + " Seeds for the exponential RNG (interā€‘arrival) and uniform RNG (thinning).\n", + " \"\"\"\n", + "\n", + " # exponential sampling for next IAT\n", + " exp_rng = np.random.default_rng(seed1)\n", + "\n", + " # uniform sampling - for accept/reject\n", + " unif_rng = np.random.default_rng(seed2)\n", + " \n", + " # maximum arrival rate (smallest time between arrivals)\n", + " lambda_max = means['arrival_rate'].max()\n", + "\n", + " # -- MODIFICATION : pre-calculate the acceptance probabilities --\n", + " accept_probs = (means['arrival_rate'] / lambda_max).to_numpy()\n", + " \n", + " # Main thinning loop\n", + " for n_arrivals in itertools.count():\n", + "\n", + " interarrival_time = 0.0 \n", + " rejects = 0\n", + " \n", + " while True:\n", + " \n", + " # 1. Sample exponential interā€‘arrival time using. Ī»_max\n", + " interarrival_time += exp_rng.exponential(1 / lambda_max)\n", + "\n", + " # 2. Advance simulation clock to the candidate instant\n", + " cand_time = env.now + interarrival_time\n", + "\n", + " # 3. Look up clock time of candidate arrival\n", + " # note that the 60 and 9 are specific to this problem.\n", + " t_idx = int(cand_time // 60) % 9\n", + " \n", + " # 4. Accept with probability Ī»(t) / Ī»_max\n", + " u = unif_rng.uniform(0, 1)\n", + "\n", + " # -- MODIFICATION - use pre-calc accept prob in interval t\n", + " if u <= accept_probs[t_idx]:\n", + " break # accepted ā†’ exit inner while\n", + " # else: reject and draw another candidate\n", + " rejects += 1\n", + "\n", + " # yield the accepted interā€‘arrival time\n", + " yield env.timeout(interarrival_time)\n", + "\n", + " # debugging\n", + " if audit is not None:\n", + " audit[-1] += 1\n", + " else:\n", + " # Show how many candidates were rejected before this acceptance\n", + " print(\n", + " f'arrival {n_arrivals} at {env.now:.2f}. '\n", + " f'Rejected samples = {rejects}'\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "# %%timeit # <- uncomment this line to run the code multiple times and evaluate efficiency\n", + "\n", + "RUN_LENGTH = 540\n", + "REPLICATIONS = 10_000\n", + "audit = []\n", + "\n", + "for i in range(REPLICATIONS):\n", + " # set up audit for replication.\n", + " audit.append(0)\n", + " env = simpy.Environment()\n", + "\n", + " seed_sequence = np.random.SeedSequence(i)\n", + " seeds = seed_sequence.spawn(2)\n", + " \n", + " env.process(arrivals_generator_with_thinning2(\n", + " env, arrivals, audit, seeds[0], seeds[1]))\n", + " env.run(RUN_LENGTH)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "53.15\n" + ] + } + ], + "source": [ + "# mean is the same, but we cut down on computations!\n", + "print(np.array(audit).mean().round(2))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.15" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From af0cfa3b1639b252782f1d13038079791dfe04bf Mon Sep 17 00:00:00 2001 From: TomMonks Date: Mon, 27 Apr 2026 16:21:03 +0100 Subject: [PATCH 3/3] docs(changes): added nb 16 nspp notes --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index fc7fa8c..0d755ba 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,6 +19,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * `14_initial_conditions.ipynb`: contains explanation of manually setting up processes before a model is run * `15_resource_store.ipynb`: introduction to `Store` and `FilterStore` for advanced resource modelling * `sim_utility.py`: added `trace`, `set_trace`, and `spawn_seeds` functions to use across notebooks. +* `16_time_dependent_arrivals.ipynb` to introduce the thinning algorithm and using it with `simpy`. ## [v0.2.0 - 11/02/2024](https://github.com/pythonhealthdatascience/intro-open-sim/releases/tag/v0.2.0) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.14849934.svg)](https://doi.org/10.5281/zenodo.14849934)