{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Random Number Generation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Random number generation is essential to simulation. Before we discuss how to simulate different queuing models, we need to first describe how to generate random numbers in simulation, particularly in simulus.\n",
"\n",
"In this tutorial, we will use the `scipy.stats` module to generate random numbers to simulate the queuing models. Python's native `random` module can also be used for this purpose. (In fact, `simulus` uses the `random` module as its core random number generator.) But the `scipy.stats` module contains a huge number of probability distributions (more than 120 of them!), as well as many useful statistical functions, which are extremely useful for simulation tasks, and of course it would be more than sufficient for our tutorial.\n",
"\n",
"The `scipy.stats` module uses the pseudo-random number generator provided by `numpy`. Both `numpy` and Python's `random` module use the Mersenne Twister as the core generator. This generator produces 53-bit precision floats and has a period of $2^{19937}-1$. It has been proven that this generator has very good random properties.\n",
"\n",
"Before we start this section of the tutorial, run the following cell so that we set up the environment by importing 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 and for qmodels\n",
"import simulus, qmodels\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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Random Variables in scipy.stats"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `scipy.stats` module contains a large number of probability distributions. Let's first examine a few. One of the most common distribution used in simulation is the exponential distribution. The function `expon(scale)` creates an exponential random variable and the parameter 'scale' is the mean of the distribution. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"x = stats.expon(scale=2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can easily find the mean, the median, and the standard deviation of the random variable, using the following functions:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(2.0, 1.3862943611198906, 2.0)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x.mean(), x.median(), x.std()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also use the `stats(moments)` method to find the mean(‘m’), variance (‘v’), skew (‘s’), and kurtosis(‘k’), where the argument 'moments' specifies which we would like to have. In the following, we print them all:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"mean=2, var=4, skew=2, kurtosis=6\n"
]
}
],
"source": [
"m, v, s, k = x.stats(moments='mvsk')\n",
"print('mean=%g, var=%g, skew=%g, kurtosis=%g' % (m, v, s, k))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can plot the probability density function (pdf) and the cumulative density function (cdf) of the random variable. In the following, we first define a couple of functions to get ready for plotting:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def plot_rv(rv, title, xmin=None, xmax=None):\n",
" # find the range (ppf() is the inverse cdf)\n",
" if xmin is None: xmin = rv.ppf(0.01)\n",
" if xmax is None: xmax = rv.ppf(0.99)\n",
" xs = np.linspace(xmin, xmax, 100)\n",
"\n",
" # get the data points for pdf and cdf\n",
" ys = rv.pdf(xs)\n",
" zs = rv.cdf(xs)\n",
" \n",
" plt.fill_between(xs, ys, color='#7fc97f', alpha=0.7)\n",
" plt.plot(xs, ys, color='#7fc97f', lw=3, alpha=0.9, label='pdf')\n",
"\n",
" plt.fill_between(xs, zs, color='#beaed4', alpha=0.7)\n",
" plt.plot(xs, zs, color='#beaed4', lw=3, alpha=0.9, label='cdf')\n",
" \n",
" plt.title(title)\n",
" plt.xlim(xmin, xmax)\n",
" plt.ylim(0)\n",
" plt.legend()\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The above function takes a random variable, a string title of the plot, and the optional min and max data range for plotting. If min and max are ignored, we choose the 1% and 99% of the distribution as the range for plotting. The following function creates an exponentially distributed random variable with the given scale (mean) and uses the above function to plot it. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def plot_expon(mean):\n",
" rv = stats.expon(scale=mean)\n",
" plot_rv(rv, \"exponential (mean=%g)\" % mean)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we are ready to plot the exponential distribution:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZRcd33n/ff33tp7U6slWbIkW8LIBsMQjBUSxknwIc6MgRw8OZk8gMcsIc8wJxOykQWY5DAchhmGJ37IQOLJ4ECSx2HxYfAQBBhMzGbwgi0sb9psbW51t1q9d1fXfu/9Pn/c29XVrW6pW6pWLf19+ZTqblX167b60z/97m8RVcUYY0zrcxpdAGOMMfVhgW6MMW3CAt0YY9qEBboxxrQJC3RjjGkTFujGGNMmLNDNuiAiV4nIrIi4K7j2ZhEZuMA1XxKRf1O/EjYPEfldEflEo8thVs8C3bQlETklIrfM7atqv6p2qqpfh/d+FfAzwNcu9b0aQUSuFZGvicioiEyIyAMicl3NJX8L/DsR2dKoMpqLY4FuzOr9B+AL2rqj8jYA+4DrgCuAx6n55aSqReBbwDsbUjpz0SzQDQAicqWI3BfV2k6KyO/VnLtfRP7fmv17ReTvou13i8jDIvLXIjItIkdE5JcXve++qCZ4TET+fc25j4jIl0XkHhHJishBEdm7wjIt+1oR+UfgKuDrUTPLn4rILhFREYlF1/ymiByOXntCRP7DKr5dbwR+WFOWue/BX4rIVPR+/zI6flpERkTkXTXXJ0XkThHpF5GzIvK/RCQdnesVkW9EX/NktL2j5rU/EJH/En1eVkS+IyKbVlF2VPVxVf2cqk6oagX4S+A6EemruewHwJtX876mCaiqPdb5g/AX+0+BDwMJ4CXACeBfR+e3AiPAG4B/F53ris69G/CAPwTiwFuBaWBjdP4h4H8CKeDVwCjwhujcR4Ai8CbABT4OPLbCMi372uj8KeCWmv1dgAKxaP/NwDWAAK8H8sBronM3AwPLfK86ovfZXHNs7nvwm1FZPgb0A3cBSeBfAVmgM7r+LwlryBuBLuDrwMejc33ArwOZ6Nz/Bv6p5rN+ABwHrgXS0f5/rzk/dZ7HB5f5mv4NcGbRsdcAE43+u2mPVf4sN7oA9mj8A/g5oH/RsQ8Bf1+z/+vAaWAM+IWa4+8GhgCpOfY48A5gJ+AThX907uPAP0TbHwEerDl3PVBYSZnO99po/7yBvsT34J+A34+2zxfo26P3SS36HrxQs/8vomuuqDk2TvgLTYAccE3NudcBJ5f5vFcDkzX7PwD+vGb/PwLfvoT/9zuAQeDti47vAfxG/920x+oeMYyBq4ErRWSq5pgL/Khm/+vAXwFHVfXHi14/qFEKRF4EroweE6qaXXRub83+cM12HkhFzSIrKdOSr1VVb5mvs0pE3gj8Z8KarkNYI372Qq8jrOlCWHsu1hw/W7NdAFDVxcc6gc3RZ/1URKrFIfzaEJEMYQ3+VqB37rNExNX5G7qLv+7OFZT7HCKyGfgO8D9V9UuLTncR/kvLtBBrQzcQ1rxPquqGmkeXqr6p5pr/ChwGtonI2xe9frvUpBNh+/VQ9NgoIl2Lzg3WqUzns+wNSxFJAvcBdxLWojcA9xMG6/nfVDXHfJPHxRgjDPdX1HxdPao6F8p/RHiz8udUtRv4pblir+TNo3sGyz3+U811vYRhvk9V/+sSb/Vy4OmL/BpNg1igGwibSLIi8gERSYuIKyKvFJGfBRCRXyJsH34n8C7gr0Rke83rtwC/JyJxEfkNwjC4X1VPA48AHxeRVNTd77eAz19qmVbgLGG7+1IShG3bo4AX1db/1QrfF8Lwf/0qrq9S1YCwW+BfznULFJHtIvKvo0u6CAN/SkQ2Ev4rYjXv33mex3+LPq8beAB4WFU/uMxbvZ6wp4tpIRbohuif8r9K2F57krAW+VmgJ/rhvwd4n6oOquqPgM8Bf19TK/8JYZvrGGFN/t+q6nh07u2E7ddDwFeB/6yqD15KmVb4ZX0c+POo18kfL3rvLPB7wJeBSeB2wpuUK3U3YT/tFdWal/AB4BjwmIjMAA8S1soB/gfhzc4x4DHg2xf5Gefza8DPAr+5qAZ/FYCIpAhvNv9/a/DZZg3JwqZPY1ZHRN4N/N+q+guNLsvlJCJfBL6sqv/U6LLUm4j8LrBTVf+00WUxq2M3RY25CKp6e6PLsFZU9a8aXQZzcazJxRhj2oQ1uRhjTJuwGroxxrSJhrWhb9q0SXft2tWojzfGmJb005/+dExVNy91rmGBvmvXLvbv39+ojzfGmJYkIi8ud86aXIwxpk1YoBtjTJuwQDfGmDZxwTZ0CRcy+FVgRFVfucR5AT5FOFQ4D7xbVZ+8mMJUKhUGBgYoFosXvrjJpVIpduzYQTweb3RRjDHrxEpuiv4D8NeE83ks5Y2E83jsIZzD+m+i51UbGBigq6uLXbt2cfHTZDSeqjI+Ps7AwAC7d+9udHGMMevEBZtcVPUhYOI8l9wG3KOhx4ANIrLtYgpTLBbp6+tr6TAHEBH6+vra4l8axpjWUY829O2Ec1fPGYiOXZRWD/M57fJ1GGNax2Xthy4i7wXeC3DVVVddzo82xpjLQlUJAiXwlcAPwucgwJ/bD2qP1xybe00Q7mt1Pzym0TXnU49AHyRcO3LO3BqFS32hdxPOJc3evXtbdhKZH/zgB9x555184xvfoFQq8eY3v5mxsTE+9KEP8da3vrXRxTPGrMBc8PpeED78oLod+AG+p/PH/IBg7tnX6DkMad+LAtkL8IPw/KUVbO5JzzkGkIgnEsu9tB6Bvg94n4jcS3gzdFpVz9ThfVvCgQMHAHjqqacaXBJj1pcgCPAqAV7Fx68EeJ4f7Qf4no/nBdXjc0HtefOh7XsBq5qcUMM/dMH+glMLDi54az33dQs+eRXl2LJp287lzq2k2+KXCFdB3yQiA4RLYsXDMuj/IlyO602EK7DkCZcqa1mnTp3i1ltv5cYbb+TJJ5/kFa94Bffccw8PPfQQf/AHf0Amk+EXfiFcy2FkZIQ77riD0dFRXv3qV3PfffdxzTXXNPgrMKa1qCpexadS9vHKPpWKj1cOqFQ8vHIY2HOPSiXAK/t4nr+ymrDWRHDtU034VrN00bXzAXxxNW4RAYkWgxXBEcJ9keoxccJtif4QWfi6xe+hQcD5Vsq6YKCr6uIFgRefV+B3VvxVrtCdB+6s91tW/fENf3ze80ePHuVzn/scN910E+95z3v45Cc/yWc+8xm+973v8dKXvrTarLJlyxY++9nPVptfjDEh3w+olDzKJZ9KyaNS9ilHz17Zp1IOz80F+HnNBe0SgTyXtdWa9txlqwjhamiK4IiEgSuEYSuCOFEILw7cpfahur0WAv/872srFi1h586d3HTTTQDccccdfPrTn2b37t3s2bOneuzuu+9uZBGNaQjfDygXvfBRWuK55FEp+fh+sPQbLA7npYJZz61FL2cuTJ1q+IYh6zhSE8jzwbsgqNcweBvFAn0Ji/8nT09PN6gkxlw+QaCUix6lQoVSsUKp4FEuVigVPUpRiHvL1aa1NpRhLqQXBvR5wlkWBbEjODWPBcHshDXptawJt6qmDfQLNYuspf7+fh599FFe97rX8cUvfpFbbrmFz3zmMxw/fpxrrrmGL33pSw0rmzEXSzUM7GK+QrFQoTT3HD3KRe/c+nBNjXoulFca1EsFs+MKjuPMn6upMZtL17SB3kjXXXcdd911F+95z3u4/vrr+fSnP82NN97Im9/8ZjKZDL/4i79INpttdDGNOYeqUip4FPJlirkyxXyFQvRcKlTO7cdcG9grDOtqSLvOfEjPbTsLmzjM5WWBvoRYLMbnP//5BcduvfVWjhw5cs61N998MzfffPNlKpkxIc/zKcyWKeTK889RcC8O7WpA14Z1tYlk4bULwtoV3CXD2oK6WVmgG9PEfC8gP1siP1sOn7NlCrMlSkVv4YWLati1wV1LRKKgdqNnBzc2X9O2sG5tFuiL7Nq1i+eee67RxTDrjKpSKnrkZkrkZorksyVy2RLFfOWc6xaEd3BuTXs+tB3cuINbrW2HbdemfVmgG3OZqSrFfIXZ6SK5mSKz02GIe15NV7+amvZ8cC+sbTuu4MYdYrGwlu3GnGrziFmfLNCNWWPlokd2qkB2uhiG+PTC8K6tdQdzNe6a8HZcIRZ3cONuNbytecQsxQLdmDoKAiU3UyQ7WSA7VSQ7VVjQ3j0X3kFU49YAqs0lEt6EXBzeFtxmpSzQjbkEnueTnSwyM1kgO5knO1Wc72VS02yyuOYtIrgxh1jKIZZwqzcmLbzNpbBAvwidnZ3Mzs4C8Cd/8ifcf//9vOlNb+Iv/uIvGlwys9Z8L2BmMs/0eIHpiTy5mWLN5E5KMNfevaDNW3BjQizuEks4xOKuhbdZExbol+juu+9mYmIC13UbXRSzBoJAyU4VmB7PMz2WJztdrAb1XJPJXA28Wvt2hFjCJR6FtzWbmMvFAn0J99xzD3feeSciwqte9So++tGPcvvttzM7O8ttt91Wve4tb3kLs7Oz3Hjjjba4RRsp5stMjuaYGsszPZ6fn2iqpgZeG+CO65BIucTjLvGEi+PWY2VHY1avaQP94W8dXbP3vumN1y177uDBg3zsYx/jkUceYdOmTUxMTPDud7+b3/7t3+ad73wnd911V/Xaffv20dnZaYtbtLggUGYm8kyO5pgczVHIlcMTtT1PappQqgGecIknYtZN0DSNpg30Rvne977Hb/zGb7Bp0yYANm7cyMMPP8x9990HwDve8Q4+8IEPNLKIpg68is/kSI6JkVkmx3L4UTfCsBllPsQhbEKJJ6MAT1qAm+Zlgb5C1gba+spFj/GzWSbOzjI9UTin/3e1GUWEWMwh0RHVwK3Pt2kRTRvo52sWWUtveMMb+LVf+zXe//7309fXx8TEBDfddBP33nsvd9xxB1/4whcaUi5zccpFj7HhLOPDWWYmC+FBhUAVjVZYh7AWnki6JFIx4gnXAty0pKYN9EZ5xStewZ/92Z/x+te/Htd1ueGGG/jUpz7F7bffzic+8YkFN0VNc6qUPcaHZxk7k2V6Ih8ejAbzBEFQbUpxXIdUR4Jkymrhpj1YoC/hXe96F+9617sWHHv00Uer2x/72Meq23P90U1j+X7A5EiO0cFpJsfyC0Zk1oa4G3NIZuIkUjHrjWLajgW6aVmqyuxUkbOD04ydyYY3NufaxGuaU9yYQ7IjTiIVtxuapq1ZoJuWUy55jA7OcHZgutrFcK5nSuAroDgxh0wmYTVxs640XaCralu0ZS6e6tRcGlVlejzP8OlpJs7ORk0q8yGuqogjpDIxkuk4bsxC3Kw/TRXoqVSK8fFx+vr6WjrUVZXx8XFSqVSji9LyKmWfkcFphvunqos9aKD4ftQuLkI84ZLKxIjFrXeKWd+aKtB37NjBwMAAo6OjjS7KJUulUuzYsaPRxWhZuWyJM6cmGR2aifqHQxAE+H5YM3diDunOsIeKhbgxoaYK9Hg8zu7duxtdDNMgqsrkaI6hU5NMj4fdDefbxgNEwr7i6Y6ENakYs4SmCnSzPgV+wOjQDIMnJ8ObnHPdDf0AVcVxHTJdSZJpq40bcz4W6KZhvIrP8Okphk5OUin75zSrxBIu6Y64tY0bs0IW6Oayq5Q9hk5NMfziJJ4XzPcb9wMQIZmKkeqI41p3Q2NWxQLdXDaVssfgyUmGX5wKe6nUBLk4QqojQSpjg3+MuVgW6GbNVco+QycnOFMb5F44klMcsfZxY+rEAt2sGd8LGDo1ydDJifmmlSjIHdehoztJwrodGlM3Fuim7oJAOXt6itPHxqmU/QVNK47r0NGTJJG0IDem3lZ010lEbhWRoyJyTEQ+uMT5q0Tk+yJyQESeEZE31b+optmpKmNnshz40UlOHBqhUvLwvQCv7IMqHd1JevrSJFNxC3Nj1sAFa+gi4gJ3Ab8CDABPiMg+VT1Uc9mfA19W1b8RkeuB+4Fda1Be06SyUwVOHh4hO1UEBT8ICDxrIzfmclpJk8trgWOqegJARO4FbgNqA12B7mi7BxiqZyFN8yoVK7x4dIzRoRkgHCTk+4oA6Y44qY6EBbkxl8lKAn07cLpmfwD4uUXXfAT4joj8LtAB3LLUG4nIe4H3Alx11VWrLatpIkGgnDk1yelj49WJsvzoxmcyHSfdmbDuh8ZcZvUaufF24B9UdQfwJuAfReSc91bVu1V1r6ru3bx5c50+2lxuU2M5nvrxKU4dHcXzfLxK+HBjDt19aTq6kxbmxjTASmrog8DOmv0d0bFavwXcCqCqj4pICtgEjNSjkKY5lEseJw+PMHYmOz9M34t6rmwIe64YYxpnJT+BTwB7RGQ3YZC/Dbh90TX9wC8D/yAiLwdSQOvPgWuAsPfKyMA0p46Mhv3JA8XzwuXd0p3h6E5rJzem8S4Y6Krqicj7gAcAF/g7VT0oIh8F9qvqPuCPgL8VkT8kvEH6brUle9pCIVfm2LPDzEwWwt4rXkAQBMQTbti0YvOtGNM0VvRvZFW9n7ArYu2xD9dsHwJuqm/RTCOpKkOnJul/fqy6zJvvB4hA54aUNa8Y04Tsp9Kco5Ar88IzZ8hOFVGNeq8EYe+VTJd1QzSmWVmgmypVZbh/ilNHRqsLTMzd9OzcmCIedxtdRGPMeVigGyAcIHTsmWGmxvMLa+WZOJlOq5Ub0wos0A3jZ7Mce/YsXsUP28q9AMcVujamiFmt3JiWYYG+jvl+wKnDIwyfnl7Qg8Xayo1pTRbo61Q+W+LoU0PkZ8vVfuUi0LUhRdx6sBjTkuwndx0aGZzm+MGz1TnKfS8glnDp7EnZkH1jWpgF+jri+wEnD41wdiBsYvE8Hw0g3ZUkZdPbGtPyLNDXiWKhwpEnB8nNlKpNLI5gNz6NaSMW6OvA1FiOo0+diXqxhE0s8aiJRayJxZi2YYHexuaG7794ZBSt6cViE2oZ054s0NtUEAQcf+4sI4Mz4UChSjg7YldvinjC/rcb047sJ7sNlUseR54cIjtVqLaXu67Q2ZPCjdnsiMa0Kwv0NpPLlji8f4BS0ZtvL09G7eXWxGJMW7NAbyNTYzmOHBjC98IgD/zAFmo2Zh2xQG8TZwemOf7c2aiJxUcVOntSJFL2v9iY9cJ+2lucqnL62Dinj42DKl5087O71/qXG7PeWKC3MFXlxKERhvun0CjMHUfo2pDEjVmYG7PeWKC3qCAIeOHpYcaGszU9WRy6e22wkDHrlQV6C/K9gMNPDjI9ng/DvOKHIz83WE8WY9YzC/QW41V8Du0fJDtViBaj8EmkYnR0Jy3MjVnnLNBbSKXscfCJAXIzpWof81QmTtqWiDPGYIHeMsolj4OPnyY/W66GebojQboz0eiiGWOahAV6CygXPZ57/DSFXBnfDwi8gExXklQm3uiiGWOaiAV6k1sQ5tHoz0x3klTawtwYs5DN1NTELMyNMathgd6kyiULc2PM6ligNyGv4nPoiQELc2PMqligNxnPC8M8lw27JlqYG2NWygK9ifh+wOH9g2Sni9WuiZkuC3NjzMpYoDeJIFCePzDEzOTcCND5tT+NMWYlLNCbgKpy/LlhJkZzBEE4nD/VESfdYYOGjDErZ4HeBF58fixczDlQ/IpPMh0n05lsdLGMMS1mRYEuIreKyFEROSYiH1zmmv9LRA6JyEER+WJ9i9m+hk5OMHhiojqfeTzpkumymrkxZvUuOFJURFzgLuBXgAHgCRHZp6qHaq7ZA3wIuElVJ0Vky1oVuJ2MD2c5eWS0GuaxuGOLORtjLtpKauivBY6p6glVLQP3ArctuubfA3ep6iSAqo7Ut5jtJztZ4Pmnz4CCXwlwHaGr18LcGHPxVhLo24HTNfsD0bFa1wLXisjDIvKYiNy61BuJyHtFZL+I7B8dHb24EreBQq7M4ScHCfxwQWcEujamLcyNMZekXjdFY8Ae4Gbg7cDfisiGxRep6t2quldV927evLlOH91avIrP4f2DVMo+vh+AQldvCseWjTPGXKKVBPogsLNmf0d0rNYAsE9VK6p6EnieMOBNjSBQjh4YopCfH9LfuSFJzBZ0NsbUwUoC/Qlgj4jsFpEE8DZg36Jr/omwdo6IbCJsgjlRx3K2hVOHR5gazxP4Wh3SH0/YDMbGmPq4YKCrqge8D3gAOAx8WVUPishHReQt0WUPAOMicgj4PvAnqjq+VoVuRcP9U5zpnwr7mns+yUzchvQbY+pqRdVDVb0fuH/RsQ/XbCvw/uhhFpkez3Pi0Nmwe6IXEE+4ZGzpOGNMndlI0TVWKlQ4+tQQGoTdEx1H6Nxg3RONMfVngb6GAj/gyIGh+R4tQNeGpIW5MWZNWKCvEVXl+MERZqOpcAM/oLMniWs9Wowxa8QCfY2cPT3NyOB0dBM0IN2RIJ60Hi3GmLVjgb4GZqeLnDg0Mn8TNOmS6rAeLcaYtWWBXmdexefogaFqzdxxxCbcMsZcFhbodaSqvPDMMMVCBd8PUIVOuwlqjLlMLNDr6MypSSZGZgmCcCRoR1fChvUbYy4bC/Q6yU4VOHV0DFXFrwQkUjGSNhLUGHMZWaDXgVfxef6pM/Pt5q7Q0W1LyBljLi8L9Do4cWjE2s2NMQ1ngX6JRgdnGB0KF3gO/ICMtZsbYxrEAv0SFPNljh88C8p8f3NrNzfGNIgF+kVSVZ5/+gy+H+B5PiLQ2Z1qdLGMMetYwwI9W8kSzrrbmgZPTJCdCudp0UDp7EkitoycMaaBGhbo+UqegxMHG/XxlyQ3U6T/hfGwi6IXkMrEbeUhY0zDNSzQFeXBgQeZKE40qggXJfADnn96uBrmjuuQtsUqjDFNoHGBrkqukuObp76JH/iNKsaq9b8wRn62FDa1KGFTi3VRNMY0gYbeFPUDn8HcID8+8+NGFmPFspMFBk9O1kyJGycWty6Kxpjm0LBAdx0XEaEclHn87OO8mH2xUUVZkcAPeOHZYVDwvQA35pDKWBdFY0zzaFigOzikY2kASn6Jb576JrOV2UYV54L6j41TyJXxgwDFmlqMMc2ncU0uAptTm4lJjICAqdIU95+6n0CDhhVpObPTRYZOTqCqBFFTiy0lZ4xpNg1tQ3cdly2ZLQiCpx4nZ07y6PCjjSzSOYIgnONcA/Ar1tRijGleDR8pmoll6E32IgiloMQjZx7h1MypRherauD4eNirJWpq6bCmFmNMk2p4oAP0JnsXtqe/+E2y5WyDSwWF2TIDxyeqfc7TmbhNvGWMaVpNEegiwpb0FuISJyBgujTN109+HS/wGlYmVeX4wbMLBhDZQs/GmGbWFIEOEHNibMlswcHBU4/+2X6+P/j9hpVnZHCG6Yk8ga/zc7VYU4sxpok1TaADpGNpNqY2IoT90w+MHuDZ8WcvezkqJY9TR0bDPud+QDJtA4iMMc2vqQIdoCfRQ2e8M7xJ6pf459P/zJncmctahlNHR/EqPr4fIAKZLpurxRjT/Jou0EWEzenNJNwEilLwCnzt5Ncu26Cj6Yk8I4PzKxB1dFtTizGmNTRdoAM44rA1vRVXXHz1mShO8LUTX6MSVNb0c4NAOXFopDq8P5ZwSSRtWlxjTGtoykAHiLtxrshcUb1Jenr2NN/p/86aLoox3D9FPjvf57yzO7lmn2WMMfXWtIEO4aCjvlRf9SbpwYmDPH728TX5rHLJo/+FsZpFK2I4blN/e4wxZoEVJZaI3CoiR0XkmIh88DzX/bqIqIjsrVcBuxPddCe6qzdJHxp6iKOTR+v19lWnjoziewGBpziOkO6wG6HGmNZywUAXERe4C3gjcD3wdhG5fonruoDfB35SzwKKCJtSm0jH0ihaHUk6ODtYt8+YmSwwOhTdCA0CMnYj1BjTglZSQ38tcExVT6hqGbgXuG2J6/4L8AmgWMfyAWGoX5G+gqSTJCCg4BX46omvMlmavOT3VlVOHrYbocaY1reSQN8OnK7ZH4iOVYnIa4CdqvrN872RiLxXRPaLyP7ZidV1Q3Qdl60dW4lJDF99Zsoz3HfsPvJeflXvs9jo0Ayz08XwRqiq9Tk3xrSsS77rJyIO8Engjy50rarerap7VXVv58bOVX9W3ImzNbO12vNltDjKV49/lbJfvoiShzXyF4+ORSNClaRNvmWMaWErCfRBYGfN/o7o2Jwu4JXAD0TkFPDzwL563hitlYqlqnOoV4IKp2dP8/VTX7+ohaYHTkxQLnnzI0LtRqgxpoWtJNCfAPaIyG4RSQBvA/bNnVTVaVXdpKq7VHUX8BjwFlXdvyYlBjrjnWxKbap2Zzw2fYwH+h9Y1WpHxUJlfhUiPyDdkUAcuxFqjGldFwx0VfWA9wEPAIeBL6vqQRH5qIi8Za0LuJyeZM/8whh+iecmnuOHgz9c8cCjF4+OEgTzU+Mm03Yj1BjT2laUYqp6P3D/omMfXubamy+9WCvTm+yt3iAt+SX2j+wn4Sa4adtN533d7HSRsTNZNIimxu21borGmNbX0kMh5/qod8bDG6xFv8gjZx4572hSVQ2nxiXqphh3iSesdm6MaX0tHegwv9pRJpZBUYp+kR8O/pADoweWvH5yNDe/cIV1UzTGtJGWD3SIBh5lrlgQ6g+efpBnxp5ZcJ2q8uLR+YUrEqmYLVxhjGkbbRHoEE25m9lK2k1XQ/2B/gcWhPrI4Az52TJBEPaGyXTZbIrGmPbRNoEOUah3bCXlpqrzvjzQ/wDPjj+L7wf0Px8NIvKUZDqGY90UjTFtpO3uBrrisq1jG2dyZyj6xWqoF4agUkrYICJjTNtqqxr6nLlQT7kpAgIqZY+xF3NU/DKBH5DqiNsgImNM22nLQIco1DNhqG/JXY3jxyj5ZcpuiVQm3ujiGWNM3bVtoEM4Q+OVqR1cmXtJ9diprud43nkWZe2WsjPGmEZouzb0xbontxAnnEc9784ymhxkXM5QpsQr9UakvX+nGWPWkbZOM8eLkZ7cCAqxIM5UehhECdSnX47zpDyCz+pnaTTGmGbU1oGemehD1EECB3UCNvtXsrGyFYBAfYZlgCecH1Lh4uZTN8aYZtK2ge54MdJTvaAgKniJIo44XF28ji3lcHr3QH3GOcujzgLI3nEAABNWSURBVHcpkGtwiY0x5tK0baAvrp0Hbti0IgjbSy/hylJ4ozTQgBmmeMR5kGkufY1SY4xplLYM9HNr5yWo6XYuCFvLV7Gr8HIEB9WAguZ4zPkuIww1ruDGGHMJ2jLQz62de0tet9G7gpfm/wWuxlCUipbZ7/yIk3LUujUaY1pO2wX6hWrni3X5vVybv4FEEM7/4qvPITnAs/IEAStf0s4YYxqt7QJ9pbXzWumgg+vyr6HD7wbCbo2n5QQ/cb5PieLaF9oYY+qgrQJdPHdh7Tx+/tp5rbgm2JN/9YJujROM8GPnAaYYX8NSG2NMfbRVoGemNoa1c3VQRwliF66d13IIuzWGPWCEQAMKmudR57uclhNrU2hjjKmTtgl08R0y0ahQCQQvXlxx7XzB+0Q9YK6p3iwN8NTjGXmcZ+RxfFb3S8IYYy6Xtgn09HQvErhh7VxWXztfrMffyMtyN5IOOqltV3/EeZAc2foU2hhj6qg9Aj0QMhN987XzC/RsWamkprk2dwO9lSvCj1GfaSZ52PkOQ/Rf+gcYY0wdtUWgp2c24PixqHYOQaxSt/d2cdlVfBk7i9dWByGVtcRTziM8I4/jWROMMaZJtH6ga9hVEcLaub+Kni0rJQibK1dyXe4GkkE66q8ecFpO8LDzHWZsygBjTBNo+UBPZrtxKwkkEFTAj6/dzImZoIvrcjfSW9nCXLt6Vqd52Plnjsth1AYiGWMaqLUDXaFjYhMAog5+rFz32vliMWLsKr6cq4svw1G32gvmiDzNY873yTO7tgUwxphltHSgxwsZYqUUEoQp7icuz7zmgtBX2crLcjeSqRldOsEIP3K+Tb8cs7lgjDGXXUsHerXtXB18twJyeUM0pRmuzb+abaVdSDQQqaJlnpX9PC4/sNq6MeayatlAd0tJkrkuUAEFP1FqSDkcHLaVd3Ft7gZSQQaNautjMsyPnG9zSp63tnVjzGXRsoGemZzr2eIQuD7qNLaJoyPo5mW5G7mifBXU1NYPypM86nyXGaYaWj5jTPtryUB3vBipmR5QQZRwEq4m4OCyvfQSrsvdQCroqNbWJxnjYec7HJGnrd+6MWbNtGSgpyd7EZWwq6IToDG/0UVaYK62HratOwQa9oQ5Lod5yPkWwwzYTVNjTN2tKNBF5FYROSoix0Tkg0ucf7+IHBKRZ0TkuyJydf2LGgmE9PSiKXKb0Fzb+stze+n0NzDXEyavszzp/Jj98pDNCWOMqasLBrqIuMBdwBuB64G3i8j1iy47AOxV1VcBXwH+n3oXdE4q21MzzP/SJ+Faa6kgw578z3B18WXENIES4KvPiAzxkPMtDssBKlye7pbGmPa2khr6a4FjqnpCVcvAvcBttReo6vdVNR/tPgbsqG8x5z6IcIpc5ob5r/1AonqY67d+/ezPsrm8nbmbpp56nJCj/ND5Ji/KMVvyzhhzSVYS6NuB0zX7A9Gx5fwW8K2lTojIe0Vkv4jsn51YfR/tBQOJ1niY/1qIEWdnaQ8vy72GTr+HuWaYohZ4TvbzI+fbnGXQ2teNMRelrjdFReQOYC/wF0udV9W7VXWvqu7t3Ni56vevdlWsDiS6lNI2TiboYk/+1ewuXF9dnDpQnyzT/NT5EY8532OSsUYX0xjTYmIruGYQ2FmzvyM6toCI3AL8GfB6Va37nUqnHCc5WzuQqLVq54sJQq+3hR6vj5HEIMOJfgI8fGBCRnjUeZAteiXX6qvoZkOji2uMaQErCfQngD0ispswyN8G3F57gYjcAHwGuFVVR+peSsL1QiFsOw8HErVHe7ODy9byVfRVtjKceJGxxBkCDQDhrAwyIkNs06t4qb6CLnoaXVxjTBO7YKCrqici7wMeAFzg71T1oIh8FNivqvsIm1g6gf8tIgD9qvqWupUyENLTG5q+q+KliGuCnaU9bCnv4EzyFBPxkWqwD0k/Z+Q023QnL9Xr6bIauzFmCSupoaOq9wP3Lzr24ZrtW+pcrgVS2Z6F64W6zTWQqJ6SmmZX8eVcUb6KoeRJpmNjBOozH+z9XKE7uEZfzgb6Gl1cY0wTWVGgN9Siror1Wi+02aWDDq4pvJK8k+VM8sUFwT4sA5yVATbpVl6iL6ePLch6+KYYY86r6QM9Vky3dFfFS5UJupYN9lE5w5gM000vu4OXsY2dOK05m4Mxpg6aPtAzU71A1FXR8dZF7Xwpc8FecHIMJ/qZrGljn5YJnnYe5ShPc7XuYae+hATJRhfZGHOZNXWgi+eSzPbMd1Vsw5uhq5UOOthdfDlXlnYxkhhgPD5MQHhPIU+OI/I0L8hzbNddXK17rMujMetIUwd6emZDOKuizs2q2B5dFeshqWl2lvawrbyL0fgQo4lBPCmjCkpAvxzntBynl01cHexhKztwcBtdbGPMGmreQFdIT0WzKgZCJWm186XENM628tVcUd7JZHyEkfgABXcWjdrZJ2SUSWeMBEm262526kvopLvRxTbGrIGmDfREvgO3kkA0vBkaxCqNLlJTc3Doq2xlY+UKZt1pxuJDTMVHo3Z2KFLkhBzhpByhl83s1N1s1Z3EiDe45MaYemnaQE/X3gxt4XlbLjdB6PI30OVvoFIqMRYfZiw+RMUpoRpeMSEjTMooB+VJtulOrtRd1vXRmDbQlIHueLFFC0Cvr66K9RLXJNvKV7O1fBUz7gTjiTNMxcartfaAgNNykgE5SYoMV+rVXKlX00WPhbsxLagpAz01Ha0XOrfEXJvM29IogtDj99FT6KMiJSbiI4zHhyk6uWpbe54cx+UwJ+QwnXSzTa/mSr2KDroaXXxjzAo1X6ArC5aYqySKjS5RW4lrkivKO9lS3kHeyTIRP8tkfARPKtUmmRmZJivP8oI8Sxcb2Ko72aY76LTJwYxpak0X6PFzboY29xJzrUoQOoJuOkrdbC9dQzY2yUTsLNOxcQL8aIkNYUYmmZEpXpBn6aCbK3Q7W3UHPWy0ZhljmkzTBXp6OhwIYzdDLx8Hhx6vjx6vDx+fmdg4k/FRpt3xmmXxhKxMMysznJDDJEmzRbexRbeziStwm++vkjHrTlP9FIrnkpzttpuhDeTi0uttodfbgo/HTGyCyfgoM+7Egpp7gTz9coLTcgIHl426hS1sY7Nus3Z3YxqkqQI9PdOzcGSo3QxtKJdYNdwDfGZik0zHxpiKjeNX29zD3jKjMsQYZ0AgQyebdCubdCt9bCFOorFfiDHrRPMEukJqJmpusZGhTcfBZYO3iQ3eJq5CybkzTMfGmY6NR71l5ha2FmYlS05m6ZdjgNDDRvp0C5v0CnrZZM0zxqyRpvnJipVSC6bJtZGhzUsQOv0eOv0etpdeQkmKzMQmmI6NM+tO1TTNhNdOyRhTMs4JOYzgsIE+NupmNuoWeumz0arG1EnTBPqCm6HreJrcVpTUFJsrV7K5ciUBATl3hmxskhl3grw7i+r8jVUhYEJGmJBRRA6FI1vpZaNuplf76GUTKTIN/XqMaVXNEeiBkJrpASWaJtduhrYqB6c69cCV7MajwmxsiqwbPoruwuYZAaZlnGmZ4FT0SzxFB73axwb62KB9dNOLazNFGnNBTRHoyVzXgjVDtY3XDF1vYsTZ4G1mg7cZgIqUmXWnmXWnmI1NU3ByNTX4MOLzMktBcpyhHwQEhy566NGN9LCRHu2liw22OpMxizRFoM83twherGzNLW0srgl6vc30epuhBB4eOXeanDvDbGyGvDOzoA0+jHhlWiaYlkmE4wtCvlt76aaXbt1AFxuIW3u8WccaHuiOFyOR66hZlciaW9aTGLFwnhm/D8rh4hwFJ0fOnSHnZsm5M5ScfLWLZGhxyFOtBKTppEt76GYDXdpDFz1k6LLavFkXGh7oqekeYH4iLhy94GtM+xIcMkEXmaCLzVFHJw+Pgpsl52bJu1nyzixlp7BEyENesuRllhEGqyEvOHTSRaf20Ek3ndpNB9100GVt86atNDbQ5/qe20Rc5jxixOjye+nye6vHPCoU3Fny7iwFJ0feyVJy89WpgeeFtfkZmWJGphfU5kHI0EGHdtFB9NAuMnSSJoNYrd60mIYGeqyUJlZO2kRcZtVixM8J+YCAopOn4MxSdPMUnBwFZ7ZmcY9agiDkokFQcGZB2AsOGTrp0E7SdNJBJ2ntIEMnGTpscJRpSg39W5mesb7npn4cHDJBJ5mgE2rqBh4eJSdP0c1RcPIUo0fZKdZ0oZwzN4ekMhtNRhYeZcHfzwQp0nSQ0Q7SZEjRQVozpOkgRZo4CZuN0lx2jQt0FdJZ63tu1l6MGLGgm45g4eLYAQElJ0/JKVB0ChSj7ZJTwJPyMrX68LkoBUoUmJLx6AgLAt8lRooMac2QIk2K8Dmp6fCZNElSFvqmrhoW6E7g4Pgx63tuGsbBIR10kg46zznn40XhXqTkFCg7BUpSpOQUKTtFUOXc2/dS/TOQCh6LavgLNsJfDwlSJEmT0lQ15JOkSGhqfpsUMWIW/uaCGhjo4Udb33PTjFxi1d42iykBFSlTdoqUpEjZKVF2ipSj7YqUFvWlryUL/ixKniJ5pqX2OOf8PDi4JEiGD01Vt5MkiZMkoUkSJMJtksRJWFfNdaiBge5Wm1uCuE3EZVqH4ISh6qc4t24PiuLjVcO94pQoS4myUw73o2O+eEs068x/yvyfEEj4ngVyILLoKpasEMWIEycRPjQxv119xInp/HacRPSaOI5152xJDWxDD2+GqqM277lpK4IQI04siAOdsExrYhDV9D0pU3FK0XaFslMKj0klei6jBMvU+Oc/df7PcKssZSpE96aW+iVwzs48Bzf8GqJHXOe2Y8SIR2fjuHPHNFbddqOrwutcHFxrLrpMGtrLRQLBS9i852Z9cnBIaoqkpuACdRofHy8KeE8qeM7ctkdFKvhzx8XDkwq+eMu08y+2VNALgXjhfQQK0SFZ5lXLHVh4ci7oXdzoOYar89sxXJzqebfmWhdHa7ajR+3+/HFn3f/iaHhnWt/mPTfmglxcXHXD8Idla/1zwmYfH18q1YD3xMOPtmv3PfzqsbnHyn4Z1FrmXwAIoJSljFBeeHxR9i4bxavI6Llgd2pC3q0ec3DUrbnGqbn23G1ZdJXoEseqz3LO/tLnhLkxEGthRYEuIrcCnwJc4LOq+t8XnU8C9wA3AuPAW1X11HnfEyGwof7GrAmJasUxjZHU9KpeqyhKEIW7j0/4HFT3F2774hHUHA8IouvDhxKAsspfEAu/muX3Fh4JxDv/NefJ0QtGbB0zWBYEvSy7P/9ftC/nL8QFA11EXOAu4FeAAeAJEdmnqodqLvstYFJVXyoibwM+Abz1Qu9tfc+NaT5hgLg46hKvQ30rmPtPfALmgj6Y366eC2qOL/UcoDXXafX6AI32535rrE01cfkwPX/MLvVPEX8Fr1vi9QKILNt9aSU19NcCx1T1BICI3AvcBtQG+m3AR6LtrwB/LSKi5w7Dq9Lorj3++m7zMqbdSbXFe22nNtbqf0H1l8b8frDgWUVrfhHUXqPVY+G/LnTB9aA176vhGJoF77PoeLRdW7raXzer/sWjEEu4y/6TayWBvh04XbM/APzccteoqici00AfMFZ7kYi8F3gvQDyR8P/j79xho4mWkc3OOF1d3db95zzse3Rh9j26sFb7Ho2MDC9bC76sN0VV9W7g7sv5ma1KRPaPj43ubXQ5mpl9jy7MvkcX1k7fo5UMJRsEdtbs74iOLXmNiMSAHsKbo8YYYy6TlQT6E8AeEdktIgngbcC+RdfsA94Vbf9b4Hvnaz83xhhTfxdsconaxN8HPEDYbfHvVPWgiHwU2K+q+4DPAf8oIseACcLQN5fGmqYuzL5HF2bfowtrm++RWEXaGGPag03HZowxbcIC3Rhj2oQFepMRkZ0i8n0ROSQiB0Xk9xtdpmYkIq6IHBCRbzS6LM1IRDaIyFdE5IiIHBaR1zW6TM1GRP4w+hl7TkS+JCKpRpfpUlmgNx8P+CNVvR74eeB3ROT6BpepGf0+cLjRhWhinwK+raovA34G+14tICLbgd8D9qrqKwk7fLR8Zw4L9CajqmdU9cloO0v4g7i9saVqLiKyA3gz8NlGl6UZiUgP8EuEvc9Q1bKqTjW2VE0pBqSjsTMZYKjB5blkFuhNTER2ATcAP2lsSZrO/wD+lAvOIr5u7QZGgb+PmqU+KyIdjS5UM1HVQeBOoB84A0yr6ncaW6pLZ4HepESkE7gP+ANVnWl0eZqFiPwqMKKqP210WZpYDHgN8DeqegOQAz7Y2CI1FxHpJZxUcDdwJdAhInc0tlSXzgK9CYlInDDMv6Cq/6fR5WkyNwFvEZFTwL3AG0Tk840tUtMZAAZUde5fdl8hDHgz7xbgpKqOqmoF+D/Av2xwmS6ZBXqTkXAG+88Bh1X1k40uT7NR1Q+p6g5V3UV4E+t7qtryNat6UtVh4LSIXBcd+mUWTndtwqaWnxeRTPQz98u0wY3jhi9BZ85xE/AO4FkReSo69p9U9f4Glsm0nt8FvhDNv3QC+M0Gl6epqOpPROQrwJOEPcsO0AZTANjQf2OMaRPW5GKMMW3CAt0YY9qEBboxxrQJC3RjjGkTFujGGNMmLNCNMaZNWKAbY0yb+P8BBGXKHE/O+6MAAAAASUVORK5CYII=\n",
"text/plain": [
"