{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Class 20: Monte Carlo simulations III\n",
"\n",
"---\n",
"\n",
"\n",
"\n",
"This notebook is licensed under a [Creative Commons Attribution-ShareAlike 4.0 International License](http://creativecommons.org/licenses/by-sa/4.0/)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load packages"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"from pathlib import Path\n",
"\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib import animation, rc\n",
"import numpy as np\n",
"import pandas as pd\n",
"import seaborn as sns\n",
"import scipy\n",
"from IPython.display import HTML\n",
"\n",
"rc(\"animation\", html=\"html5\")\n",
"\n",
"np.random.seed(686500035)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Random numbers from various distributions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Definitions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Distribution of numbers\n",
"\n",
"A **distribution** of numbers is a description of portion of times each possible outcome or each possible range of outcomes occurs on the average."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Discrete distribution\n",
"\n",
"A **discrete distribution** is a distribution with discrete values."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Continuous distribution\n",
"\n",
"A **continuous distribution** is a distribution with continuous values."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Probability function\n",
"\n",
"* For a discrete distribution, a **probability function** (or **density function** or **probability density function**) returns the probability of occurrence of a particular argument.\n",
"\n",
"* For a continuous distribution, a **probability function** (or **density function** or **probability density function**) indicates the probability that a given outcome falls inside a specific range of values."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### To Generate Random Numbers in Discrete Distribution with Equal Probabilities for Each of *n* Events\n",
"\n",
"Generate a uniform random integer from a sequence of *n* integers, where each integer corresponds to an event."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Uniform distribution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Just as likely to return value in any interval\n",
"\n",
"* In list of many random numbers, on the average each interval contains same number of generated values"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Using `numpy`\n",
"\n",
"There are two methods available in `numpy` for sampling integers from a uniform distribution.\n",
"There is `np.random.randint`, which let's you specify the integer range directly, and works like so:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([4, 3, 1, 1, 2, 4, 1, 1, 3, 1])"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.random.randint(low=1, high=5, size=10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And then there's `np.random.choice`, which is more general and allows you to randomly sample from an array-like list.\n",
"To replicate the above behavior, you would use the following inputs like so:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([4, 2, 1, 3, 4, 4, 2, 2, 2, 2])"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.random.choice(a=[1, 2, 3, 4], size=10, replace=True, p=None)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Non-uniform discrete distributions "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Manual method\n",
" \n",
"* **To generate random numbers in discrete distribution with probabilities** $p_{1}, p_{2}, \\ldots, p_{n}$ **for events** $e_{1}, e_{2}, \\ldots, e_{n}$, **respectively, where** $p_{1} + p_{2} + \\ldots + p_{n} = 1$\n",
"\n",
" Generate `rand`, uniform random floating point number in `[0, 1)` \n",
" If `rand` < $p_{1}$, then return $e_{1}$ \n",
" else if `rand` < $p_{1} + p_{2}$, then return $e_{2}$ \n",
" ... \n",
" else if `rand` < $p_{1} + p_{2} + \\ldots + p_{n} - 1$, then return $e_{n} - 1$ \n",
" else return $e_{n}$\n",
" \n",
"As an example, consider the sample distribution `[1, 2, 3, 4]` where value `1` has weight `0.10`, value `2` has weight `0.30`, value `3` has weight `0.40`, and value `4` has weight `0.20`.\n",
"Following the approach recipe, this results in the following Python code:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3\n"
]
}
],
"source": [
"x = None\n",
"rand = np.random.uniform()\n",
"p = {\n",
" 1: 0.10,\n",
" 2: 0.30,\n",
" 3: 0.40,\n",
" 4: 0.20\n",
"}\n",
"\n",
"if rand < p[1]:\n",
" x = 1\n",
"elif rand < p[1] + p[2]:\n",
" x = 2\n",
"elif rand < p[1] + p[2] + p[3]:\n",
" x = 3\n",
"else:\n",
" x = 4\n",
" \n",
"print(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Using `numpy`\n",
"\n",
"A general approach to the discrete distribution again uses `np.random.choice`.\n",
"The different weights for different values are specified using the input parameter `p`.\n",
"To replicate the prior example, you would use the following inputs:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([3, 3, 2, 4, 4, 3, 3, 3, 1, 4, 2, 3, 2, 2, 3, 3, 3, 3, 2, 3, 3, 3,\n",
" 3, 3, 3, 3, 4, 2, 2, 2, 2, 4, 4, 2, 3, 4, 3, 2, 3, 4, 3, 3, 3, 1,\n",
" 3, 3, 4, 3, 3, 3, 4, 3, 3, 3, 2, 1, 2, 2, 2, 3, 4, 2, 3, 2, 2, 3,\n",
" 4, 4, 1, 4, 3, 2, 3, 2, 3, 3, 3, 4, 2, 4, 3, 3, 2, 1, 1, 4, 3, 4,\n",
" 2, 1, 2, 3, 4, 3, 2, 3, 2, 4, 4, 2])"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.random.choice([1, 2, 3, 4], size=100, p=[0.10, 0.30, 0.40, 0.20])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Continuous distributions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Normal or Gaussian distribution\n",
"\n",
"Probability density function, where \\\\(\\mu\\\\) is mean and \\\\(\\sigma\\\\) is standard deviation:\n",
"\n",
"\\begin{equation}\n",
"\\dfrac{1}{\\sqrt{2\\pi{}\\sigma{}}}e^{-(x-\\mu{})^{2}/(2\\sigma{})}\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"mean_norm = 0\n",
"sd_norm = 1\n",
"x_norm = np.linspace(-3, 3, 200)\n",
"normal_df = pd.DataFrame({\n",
" \"x\": x_norm,\n",
" \"y\": (\n",
" (1 / np.sqrt(2 * np.pi * sd_norm)) *\n",
" np.exp(-(x_norm - mean_norm)**2 /(2 * sd_norm))\n",
" )\n",
"})"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAwoAAAHMCAYAAAB1Fl+hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAASdAAAEnQB3mYfeAAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3Xl8FdX9//HXJyFhCTsSZN9FATdUZJG64NYq1rXaVgWttta1m6222tpWW9tqaxe1iwpu/brvP21dWhWBIEpRQEXZd5B9DQnJ+f1x5obJcO/NTUgyucn7+XjcB2Tu3Lln5s72mXPO55hzDhERERERkbCcuAsgIiIiIiINjwIFERERERHZiwIFERERERHZiwIFERERERHZiwIFERERERHZiwIFERERERHZiwIFERERERHZiwIFERERERHZiwIFERERERHZiwIFERERERHZiwIFERERERHZS7O4CyCyL8ysO3B18Oc059wLcZanITOz9sB3gebAx865B2MukoiIiDRgNa5RMLNpZnZH6PU7M/uDmd1pZr1rs5Ch77y1DpZ5nJmdmOG8I8zsb6G/+5jZZaner29mdpmZ9Ynr+8PMbLqZtazj78gDfgn8yjl3Y0MPEszseTPrG2MRfgPc45y7QUHCvjOzCWY2YB8+f2vk73rbP6LHZ12cW6v4/r+Y2Zj6/M7akOl5zcxuMTM9iEvCzLoE9wu/NbPbzezPZnZE3OVKpiFdU2uLmTUzs1vq+TtrfH4xswFmNqEWixNdfq3eq5hZnpm9V1vL2xdmdoyZHb6vy9mXE9kbzrmbohPNrDVwp5nd7pxbtA/Lb3Ccc0Vm9r/qvB8ckLc653bXdfkamC8453Yl/ghOtic65+6rxe8YDTzhnNtai8vcZ8FJ7R3n3PzIW18Jb5P6FGz//znn1sTx/ZKRGu0fNTy2vlAf+6KZHQc0c869Hnnr+3EdC/uoPs5rjZaZFQC34X//zcG0HOA2Myt2zs2NtYB1KM2xIPGq8bnQzG6N3gc750rNbHTtFK1GZeoHXAB0BOYDu4GU962ZqPU+Cs65bcB3gCtre9kNQVU7VJZe/GpdPW2HnkDWBKMx7xtZta2aovrcP+I+T8X9/TWVreVuQK7A1wBvTkxwzpUDNwOXpfyUSB2pi2M6xgeCvYAzgfuccz8A/lUby62TqlHn3E4za2pP0KX+5QJlcRciS2hbiUjc9nfOLYxOdM7tNrM1ZtbBObcxjoKJZDvn3FLg97W93LpsQ1kRKCSaYgAHA8cDu4CHnHOzg/eHAt8IpjcDCoCXnHP/L7pQM2sF3AC0Db4jH3jTOfdMZL4cYDwwBHD4GyUH3OGcW5VkuWOBLwfLdPibqtudcxtC8/QhTTVz+H0zOxYYB4wC2piZA54FpgF/cM5dl2IZHYFrnXO3JHs/mKc3vlPqbvbUCk1KM+9lwXytgte9zrlZoXludc7dZGanAqcAJcE2WA/83jlXFpq3BfB9oEMwTw5QjN+uG5MssxXwC/zv1dvMDgSWOef+aGY3Af9I1RzGzP4IfDd44hSe3gNfazUEOMbMNgE7nXM3h787xTIr3kv8XsBTwfLaBuveEvh7tBo8tO6dAAvmmwf80TlXbmYHs2efS5TrTefcS6nKlem+H9qeXweOAEqDbT8vXbOH0Pbv6f+0U4O3bnLOFQfzjAa+AuwE8oL1esQ5NzVaBvxJ6If4426Fc+7OFN97C3BrUNbz8ftqebCetzvndkbm74XvlJ44VlsC051zDyUpw0/xv0O3YJk/cc4VB+/9ErgWaB+sSwvgZefcv8ysOXAV/vfLCd7/MMl3tAGuA9oEZckDVuGP25o0DToEuDxY99xgW/w5yXyV9o+qjrV0x1bw+Vvwv8EE/Lm3FH88r0yxL3YGrsf/truDbfeEc+7tdOVM9l7oGA3vdxXbOlmTzOrsA5mer5KU72/At5OcU34JzIj2cwq26WmJ/TyT81rk80eTwf6fbDvin7JfEWzDHPzv8Z5z7hEzM/w+1Rt/LsrF76N/SrJuxwCnBX8WBOW4zTn3eWiefHwrgG7B8nLw54PfhZ/+h67lzfDXlcS10oL1SnmTb75P2fY0q/0+cAjwVpplZHRs1uD8k/E1NUmZOgLfwx83FpRpDXBX8NA07bEQLONM4AvB9yfW68/Ouc9C8/SheterXPz5o2uw7gbMBfbqn2a+b9S3gz8T9wpvOecej8x3K3673gC0A7ZHzlnfBAazZzsuByodE6kEx9SP8L9t4nz3NjA7ybwt8eeKjvjt3gZ41Tn3VPD+ZUCRc25Oiu+6k9B1I7IOVR4Lwe91DDA6OE8D3O+c+zixnZKcX08HTmLPdbY5/l4s+rslrmPfAbqw5zd+0jn3TtqNWFecczV64U/y6d7/bej/E/AH0gVJ5hsB/A7ID00z/A/17ci8v8bvdF0i078BnBGZ9htgRGRaF+B5oEVo2nHAHcBlkXk7Af8AckPT+oTnq+rvYNot+HaJ4Wk3AN1SbLergEFptms74B6gZWhaDv7m7S6gT2h6X+DG8Pfjd9A/AgeHf0v8yXR85LuOBq6KTPtV+DuCab2BJ4CcVPtHim3TC99WNdl6dgRuqWIfmwAMqM6+GX4vKNM1wN34J12J6c2CbRzeTwz4U5J1PxZflV7tclVz378VfwN8amT6udF9P8V6H4cPYqPTz8SfnMO/XS7wM+CsyLy/DvaxThl83y3AWPyNZ3jZ/aK/K9Af+CvQNjL9bOAXSbbDj4GhybZvsD17JCnLkKDs0fcuBL4Y+rs5MJHI8QkMx1+Mq/ydI/P0BX5L5fNIPr6d9v1V7B+ZHmt9iBxbofW+Ejgug33xrmD/jv4GNwFH1uT4qmK/u4XK56Xq7gMZna+SfO/lwFFJpt9G5DgOpl8b3teSrF+6bZ/R/p9qO+IDhcFJyn9CsG9E3xsDfCsy7TTgnMi0QvzNYkFo2t+IXHeC3+SpJPv8uOD7w+etzsCdVazT/sDVad4/GN9PJ9X71Tk2M97+VOOamqJcfwY6RqYdhm/+kcmx8D3gy5FpbfABQWFkX8voehXaBodGph0F/CC8DYBB+HuClpF5rwG+kWS/vB3olWQ9royuX7AP/YIq7heDef8YXS7wJfy1aEJoWl6w/0W3+aXAJcH/26b6Tnyw/LvwOkXez+hYSPbZNMu8IlG20LQWwb4TvU+9LdjvBoamGf460iXZ96XZpn1Icn6q7qtOxlEws8JgxcJ6O+cei8xn+KcSP3TOlSSmO+8eoIdVzgJyBP7pbaUn0M65+/EXlMRyewLrnHNFkfnWAA/hT7RhHVzkqaxzbj3wOPDVqta3Bh7C36Akc6Bzbl6az16Jv3BWPBFxzpU7534LDI3M+3V8wLY7NG8p/kRxeWi+VkA/F8mE45ybDvSILLPMObc4Mt8S4GXg0DTl3ovz1WS9Urz9VeD/qrO8GhqPv0FYHSrXbvyF9NTQfGcDTydZ97eA0uBpbMZqsO/3BIqdc5XaHDr/BGVkdb47VIY2wAnOud+40FNI51yZc+7nwKnBPAmn4p9Yrs/wK05yzv0usuyF+BNk2PfxNUdbwhOdryVca2ajQpM7Aktc8idFZ+JrKpdHpt+J35f+En3POfcI/gKecB7+Kc/KyHzvAjsi2yMTVwA/c6Gn3MHv/VN8Z/x0auNYa+6cezOD+b6Ar2XaEpn+a/Y8aaxL1dkHqnO+ivo3cHJ4QtD5bxnQOjguw4ak2Ncyken+n8yR+NqDjyLT7wd+Dvwr+p5zbjL+hgaoSCwyyDn3dGS+tfibkcuD+Y7BJyeZF5lvATDNzAZHynA+cHPkvPU5sDp4up5Kc3yNWCo78U9OU6nusZnp9q/ONbUSM9sPmO9CLQ+Cz88CllWxPQhqObs6556PfH4r/tg7N/KRjK5XwT69yTn3QWS5M4DVVPZN4HoXqWVxzv0Z3yIi7GDgleDaHV6P1kB3F+moHexDr+NrNVIKWmC8FV2uc+5l9m4yOwFf2xLd5g8Ah5tZXnAeaRPUYkV9BXg6yfSaHAtVCq7jXZxzEyPLLMbXjl0ZOe+cA0x0odok5+/676Bu7kerVOuBQnBCvxsf/YS9lGT20cD/CzZCMn+h8oaZ6ZK0bwy8FfoRt+NvxpOZw9477ZPJZgx2+mrd/GYiONF1jV6UzGwISarZIjqETxIRFethfnyBJS5JNXwQLGwLTWoHPJximaWRv/NTzPcw8GGK99J5KzhJRA2qImCqLTOdcyuSTP8U/zQ44aggKEjmFfwTteqo7r7fnsx/o0x9Gf+ELpX7g3kSPkpz/CWTKl3tGjPrABXNyOZFL1Ihfyf0EAB/zno+xbzLXNCcMSy46P7P7Z2FKiHcn+pdfBOIZD7GP42tDku2bsFx+WwVn62NYy3ZeTeZV5MECYlyLq1BgJSxGuwD1TlfVZLi4cQ44EX8735UqFz5pL+xrUqV+38azjn3SpKJ5cBCF2kOFhJe/zNIcUPknPuUPfvyYvw5LJlk18uXk11X8BlWeqdYDuz98LC681T32Mx0+2d0TU1hU5LvTbgNqKq/heFrA5JJtu0zvV6dRorfnr3X6f5w0BcRfSi0NcV18HhSdJwN9tW1KZafcBKpz1XR8u7nkjQfD7yND2bAP1AZl2SeUS7yEDlkMdU7FjJxAf6eeC/B8fwSlR8aLY0GeMG8a/EPyurdvvRRON7Mbg/9nYuv7pmLryaKtkVclGQZh5HmQHTOrQoi1YR0HaQX4KvQPopGmhG72PtklG65ddWP4w18zcYboWnn459+ppOurOGT9+HAsUHwkUwvM8sJdtRVSZ7CJkSDyZfM7C78hXpm4kY3xYUjEy/gn5xUnHwyDJhqy7sppm/F35AkpFy/NCeddKq7789OcyNV04B/IPBomjK8a2ZnhCYlO4bTqWrbbsRvh5TbzzlXYmbhi9g65zOrJTMzTVmWpXmvIlgLbqBSSXbuSCl4mpguHW1Vx0xtHGuZ/mbpziuLgAHsY4q9NKq7D1TnfJXMZjNrGwqMujnnVpjZi/jmFon99hh8e/yaymT/T2Wf92XgQOCQvStJKpQDpNmWsKfvVFim58x0Zav2PDU4NjPd/pleU5OVabeZvW9mvwIedkEb9cR76T4bzFNMimPUObcrSQ1Xptu+q3Mu1X5SaZ2S1FqFRQOIxSnmG0iSvg+pvjOJnDTBSriPZAH+3jPVftaOPdvzNXwQVtF/1cwG4u8Vk6rBsZCJVs65dWnefwf/ICRxrkn1G0MdPNzPxL7cBP/XpejQVg0dqTriLq/i/YS1RJ5mmNlI4HT8iaB5MLkdqZ9KJLOjGvNWxyv4fhRvQEXHowIX6jiWQiYnW/Db9g7n3Cc1L2KSL3duspnNwLcdPD2o2vuvc+6NKj6aanm7zWyrmbV3zm0KJmcSMNW3TLd7pmpz36+p3DQ1GvWlI/5pWDqZbodaWRfzA2V9Gd/8Yze+TSz4J1XVuXEsIH3nzbRq+1jbB2uB1lXOVXO1uQ9k4j/4hzTPmR+tfAOAc25D8HfCWHx1fxxqY19u5py7IZMZgxvSk/BPhkvYc73sBTxQC2UBXzuTrtlVS3zzo5Rq8dgM26dt7Zx7xsxeB8aZ2fn4m7kXnHMZD7plPqnFeUFZmrMn6Klp9shqrVNw430hvh9JDv7BL2TeomKfznVkXt5OwDPOub9WuUCfYGSdmXUN1UCMx/eFSqkOjoWqHi6twyesaLDiHjlyI/4ikaraDzK/QLTBR9WJm+578W3jfh6OVEOZAzKV8RPE6nDOlZnZRjPr5Hyb71Px7WerkukIguvxvfZrNVCAiqcgzwQvzOw8M7vD+by9NfEovj/F3dUImNIWcR8+W1/LrM19v04EJ8y6Tqma2A7p1Nt2MLNu+A7Rdzvnboy8N6Gai9tG9ZukVVIHx1pNtKHyflpXx0I6tbkPvIVvFvIcPggLZxhbaGa9giZK7V12p+rMMbP8NE9qgYr25ffiz8M/caE2/RYMElZL5dkA7Jfm/W7A56nerOVjM2yfR+UNaqceDcqSC3zTzI53zv2uqs+a2W/wT7l/45zbEXmvpiMaNw+1GKjq+4/HP5z7jYsMkluN70+c69LVdqVT3fuaTD0KXAT8NvhduqWrNaijY6Gqe8gOBA8rGqpYqjFCPiBNR8yg6j7tE4aQwfhUleB3+qecc09UdZLMQF1uo0fwN8jgn15lMmJjpu3RZwXLrHPOuSeBFWY2qIafX8CeTninUkuDhKTQvOpZkkq5H5hZLzOrbtvF2tz3a+pTSz+8+yHAZ2nerw0f4rM/JRU8QUzWIa2u/AC4xkVSw9ZEcJNZq0/i9/VYq6FB+PbnmajJ8VWv+0DQhC9xYxLtrPw8vvamCz7tZjb7AP9UtCrXAL90zv0rkxvLmgr6xhWkmeUI0ve9qbVjM6KmfbyScj4ZxL1AYaT56F6CPp0LnHN/jwYJ+2gBPgFGJi7CZ9mrbtPSsE/wGaVqamvQsTutoEl7un4w0fk/Zk+H9FOo+h6rLo6FHZY+2ckI/LHaYMUdKLyDr65LFXF9G595KOHANG3TDgm1D+xC6jaeyWoTkl6kzGw4+94udxcpLp7OZzDpFdwUbsxwx9xglbPhhB0RWvYKfOacpDuo+XEmqsXMTjWzVBf0FaR/KphyOwSmm9lR7N1voyaSPp0Ibopr0hkJYKaZfSHFe2dTed2qWleo/r5fF54HLknz/gRSdwisFcExcJD5vNjJXETlJ751rThNH6fq1EQm5KTpCHxEiunVPdYy2d+qckRwQx4thwHtIv1Cmic7f5jPdndkZHKVZYtpH/jMzPoT6awc1CT0wJ+HXs1gObWx7evKC8A3Up1jQr9ha1I/EKjJPp/OiqCdeLQsiae96TKq1faxmZDRNTUZM7vcfOKQZD7H18YlJNtXUt6rmNm+rNNLpM6sGF2n1cmaoAbHY1WZ2RLexGed24v58VHS1SSBb01xfor3ouX9MNW2SXFfU2R+rKBzqDqBRF0cC4+RInNccGyeBtR28FurYg0Ughvj+4Dfmc8wUcH8wFIbXOUsK6uAW6IXXjO7kMqR4qsk+WHM7Kv46rFolD/IzMZF5m2PP9CqynpQlfmkb+c3FZ+//LE084TdA3zH9gzyAYCZXYxvvxd2E745T2Fk3hPxbfWq6x3g4ugFPQhGziF934+1VM7KEPUsPg/yplqI5OdY5U64mNn++H0iXWeldJ4Bzo9eUMzsIHzu58WhyVX95jXZ92td0LzrLTO7IXyCNe+7wGSXJBNOHbgD+FOS4/pEfFrl+hxk5jMzi6bPbG5mN+Mv/NWtIfgzcGNk++aY2fdI35ymOsdaVcdWJj4Hfpnkqd4P2LvD+8v48V7C5WqHHyQomoRgEVWkmAzU9z7wb/z5MVkGl0342r5M2pjXxravE8HT1/uBPyQ5x3wP34kc/E3epZH3c83s6uDP2qwV+ztwg4XShgbHxq/w4xalU9vHZkJ1rqlRrwSfrRRkm09POtRVzs6T7FiYAnw9eoNrZifhf59Umc/SCjoyO/N9H8LLHYDPhhWWZz6lfHi+Qvx4CSk7/ka+bycwN2jGFF3O1fhjKt3npwJDzWdAC3/+KPz4IGF3AROC98Lz9sGfg6Iew4+HsjNNQpCEN6nesbAxTaAIVLSYWG9m0eXm4celmFSXNXm1Ie4+CjjnpprZVuBO85ktSvAdnt5wzkUvUJvxI8PeGBzUZfh1eNOFUsk55+aaWQ8zuxvfbyEneL2EH6zou5HlPggUmNkf8Bdvh6+O/Fkt/IDPA783s7OBT9zeo+i+iB8kJ6MmHs65reZHFbzNzMrxnY4MfyP7amTeZWZ2LfBD8yNDJ0Z9nOUiOX0z/O5tZvaLYHnN8c0BcvDb6zvpmnkFfTLmmtmf8L/bT8JVrUFmkxwyD5jSeRD4kfmRELfj95Ht+FztP6rJAp1zzsx+APwguMglRvBeih+YJzzvDDM718yGAZ87536TYpnV2ffrhHPuafO5o/9sZjvx+30LfNO9KfVUhs/Mt4X9RbCfluJ/sw9cmhHK66gsE83sWjP7Iv4JYA7+nPA3fFvSavU5cM4tNbN/4s8BpezpJPgAaWrgqnOsVXVsZWgp/ibu1uDpbhn+JuVJFxrFPfi+N81sgJndh2+bnBPM/zP8TUF43tVm5oJz1m7nXNLjr773geAacQj+Ri3qZeDGTM79tbTt64xz7t9mtg4fLGzD78t5+BG3ZwbzvGZmXwvWYQf+9zR809i1+I6dtVWeHWZ2I/Dj4Hxfjh8X40GXJCVk5LO1emyGlpvxNTXJZ5ebH+3758H6JJrJbSGDY8E5tzY4P/w5uBYk7lWmOufuMLOMOqOncDs+iLkEKkY6XoofAC3cx+OX+N/DQvNtBn4CfNXMOqapyQmv331m9o3guuvw23Aj/rzwkwzK+zPgejNrG/ydg3/wcBc+A1nie8qCdfpe8FCtOPiu9fgxRqLl2hCc06ocl6kGx8Ik4FfBsTXF+XFfki33bjM7w8zuCZZbij+/TnQ1H6el3liSGiepR+abwxyRJIBocszsbufcVVXPKSIiIiJ1Le4+CuLb5T0RdyHiZmZHkD53uIiIiIjUIwUKMTLfqWt3PbUDb7CCKs8rgX/GXRYRERER8WLvo9AUmdlp+LR1zajcVrBJMZ+C8Af4QfD+mkFHIxERERGpJ6pRiDCzM81slpkVmx/Vb1I4Q0MGn88xs2lBp6U+kfeOMbMp+ExKl+FH+ivceylNg3NujXPueufcN51zM+Iuj4iIiIjsoUAhxHyu3afxKcHOw/fCPwt4qhqLuZIk4zIE6dJexfd2/zrwHeBo4N8WSV0nIiIiIhI3NT2q7Dr8EOTnJVLjBSnLHjSzQ5xz6UaNJMinexs+5V50oJIr8KncxjnntgbzL8Tn8T6VOh7YSkRERESkOlSjUNkI4LVI/uzXgn+HZ/D5u/F5f+9PsexpiSAh8A4+B3AmyxYRERERqTeqUahsf/ygGmGJv9P2JTCzc4Av4UdTTHbjvz9QqR2+c648GAwn5bKDkQ2jA8m0Bg4A5uAH6RIRERFpqPKBnsBbzrnNcRdGMqdAobI8fB+CCsEogJBmKHUzawf8CfiNc+4jM0sWKORHlx1IjICaypX4vhIiIiIi2ezLqKl1VlGgkLl0Q1j/BtiO759Q28u+B58lKexA4KnnnnuOAQMG1PArRaSpcM4xb/VW3pz3OW9+upb122pWEdmpdT7HHVDI8Qd25oAubQgeooiIpDV//nzOPPNM8P1AJYsoUKisFGgZnhDKSLQr2QeCTEmXA18EdptZM/b0/cg1s1znXBm+iVDLJIvIT7VsAOfcWiLNoRIX5wEDBjBkyJAqVklEmqrSsnJemLWSv7+9kHlrtgK50LIr+cGZqE3zZgzs0pruHVrRtV0L8nP9qaukrJxVm4tZsXEHn63ZxtZduwHYCry4DF5ctoFBXUr55hf6ccZh3cjLVXc3EcmImktnGQUKlX2O70sQ1iX4d1WKz9yEDwz+neS9+fisRsclW7aZ5QD7pVm2iEi1lZU7Hp+xjL/85zNWbi6umG4GR/ftyBeHdmV4344c0KUNuTnpawXKyh2frtnKu4s28MqcVUxftAHnYN6arXz/yQ+489V5XH3CQM4/qmeVyxIRkeyiQKGyacApZpbnnEv0Jzgl+LcoxWe+A7SJTDsd36/gDODT0LKvMbP9nHPrgmkn4PtFpFq2iEi1FC1cz89f/IiPV22pmNa5TXMmjOrDOcN6sH+7FtVaXm6OcVDXthzUtS3jR/Vh9eZinp65nElTF/P51l2s3FzMj5+dzcNFS7hl3GCO7teptldJRERiYs6lax7ftJjZMfgagFeBe4FuwO3AFOfcacE8Ofg+Ap+FgonociYAE4G+zrnFwbT++NSps4HfAq2AXwEbgcOD5kmZlnMIMGfOnDlqeiQiAGwtLuUXL37Ek+8vr5jWvX1Lrj5hAGcd3p0Webm1+n3FpWU8+78V/OU/81mxaWfF9POO6MFPxw2mTYu8Wv0+Eclec+fOZejQoQBDnXNz4y6PZE4NS0Occ+8AZwNd8R2Ib8WPyvy10GzDganAuGouewFwMrAbeBT4C74m4YvVCRJERKKmL1zPF/84uSJIaJGXw/dOOoA3vn8sXx3eq9aDBP8duXx1eC9e/96xfPfEA2iR5y8nT76/nC/+cTLvLtpQ698pIiL1SzUKWUg1CiICUF7uuPu/8/n965+SOJUfM2A/fnPuIXRvnyx3Qt1ZsWknP3rqQ96Z71tWmsH3TjyAq44fQI76Log0aapRyF6qURARyULbdu3m24++z52v+SChRV4OPz9jCA9dOrzegwTwzZweunQ4t4wbTPNmOTgHd772KVc+OpNtQdYkERHJLgoURESyzIpNOzn7nin8e+4aAPp0asULVx/D+FF9Yn16n5NjTBjdlxevOYY+nVoB8K+5qzn7nimV+jGIiEh2UKAgIpJFPl2zlXPumcqna7YBcNygzjx/1TEc0CWafC0+B3Rpw/NXHcOxB3QG4NM124Iyb425ZCIiUh0KFEREssT7SzZw3l+nsXqLHxvh0tF9uX/8UbRr1fAyDLVrlccDE47iktF9AFi9pZjz/jqN95eok7OISLZQoCAikgXeXbSBC+97l807fVbmH516IDefflCDHuQsN8f46emD+eGpgwDYvLOUC+97VxmRRESyhAIFEZEG7r3FG5gw8V12lpaRY/Dbcw7h28f1x6zhBgkJZsaVxw3gt+ccQo7BztIyLpn4rmoWRESygAIFEZEGbObSjUyYOIMdJWWYwe+/chhfOapn3MWqtq8c1ZM7v3IoZrC9pIzxD8zgf0s3xl0sERFJQ4GCiEgD9dmarVwycQbbdu3GDH537qGceXj3uItVY2cd3oPfneuDhW27djNh4gzmr1UHZxGRhkqBgohIA7R6czHjH9jTJ+H2sw/m3CN6xFyqfXfuET349VkHA77PwvgHZrAm6JwtIiINiwIFEZEGZvPOUiZMfJeVm/0N9PWnDOL8o3rFXKrac8HwXlx/iu/gvGLTTsY/8C5biktjLpWIiEQpUBARaUB2l5Vz9T9n8slq3yTnohG9ufK4/jGXqvZdeVx/LhrRG4BPVm/lqkdnsrvDcZ1mAAAgAElEQVSsPOZSiYhImAIFEZEG5Fcvf8Lkz9YBcPLgLtxyxpCsyG5UXWbGLWcM4eTBXQCY/Nk6fv3KJzGXSkREwhQoiIg0EE+8t4wHpiwC4KCubbnrgsMa9DgJ+yo3x7jrgsM4qGtbAO5/ZxFPvrcs5lKJiEiCAgURkQZg5tKN3PTsHAA6FeTzj4uPoFV+s5hLVfda5TfjHxcfQaeCfAB+8uwcZiptqohIg6BAQUQkZhu2l3DVozMpKSsnL9e498Ij6NGhVdzFqjc9OrTi3guPoFmOUVJWztWPzmTj9pK4iyUi0uQpUBARiVF5ueO7j89iVZDh6KbTBjO8b8eYS1X/hvftyM2nDwZg5eZivvvELMrLXcylEhFp2hQoiIjE6J435/PWp58DcNohXbl4ZO+YSxSfi0f25rSDuwLw5rzPufetBTGXSESkaVOgICISk3cXbeD3r30KQN/9Crj97IMbZYajTJkZt59zMH33KwDgzlfn8e6iDTGXSkSk6VKgICISg807S/nu47Mod9C8WQ73fH0YbVrkxV2s2LVpkcfdXxtGfrMcyh189/FZGoxNRCQmChRERGLw0+fnsGLTTgB+/KWDKlKECgzu1pYff/FAwI/c/NPn5sRcIhGRpkmBgohIPXvufyt4ftZKAI4f1LlJ90tIZfyoPhx7QGcAnpu1kudnrYi5RCIiTY8CBRGRerRy005ufm7PeAm/PffQJt0vIRUz43fnHULHYHyFm56dw6rNO2MulYhI06JAQUSknjjnuOGZ2WzdtRuA35xzCJ3bNI+5VA1XYZsW/PacQwDYums3Nzw9G+eUMlVEpL4oUBARqSdPvLeMt4NUqOce0YMTB3eJuUQN34mDu3DOsB4AvPXp5zz53vKYSyQi0nQoUBARqQcrN+3k1pc+BqBL2+YVg4tJ1X56+mC6tPU1L7986SM1QRIRqScKFERE6phzjh8/u6fJ0e1nH0K7lkqFmql2rfL49dkHA74J0o3PqAmSiEh9UKAgIlLHXvhgJW/O802OzhnWg+MPLIy5RNnnhAP3NEF6c97nvPjhqphLJCLS+ClQEBGpQ5t2lPCLFz8CYL/W+dx8+kExlyh73Xz6QXQKsiD94sW5bNpREnOJREQaNwUKIiJ16Ncvf8L67f6G9ubTB9O+VX7MJcpe7Vvl89Nxvm/Hum0l3P7KJzGXSESkcVOgICJSR4oWrufx95YBcOwBnTnj0G4xlyj7nXFoN74QDMT22IxlTF+4PuYSiYg0XgoURETqQGlZOTcFA6u1yMvh1jOHamC1WmBm3HbmUFrk+cvXTc/NobSsPOZSiYg0TgoURETqwMQpi5i/dhsA1409gJ4dW8VcosajZ8dWXDt2IACfrd3GpCmL4y2QiEgjpUAhwszONLNZZlZsZuvMbJKZdaziMwVmdquZzTezHWb2sZldb2bNIvP90MxcktdjdbtWIlKfVm8u5o+vfwZAv84FfOOYvjGXqPG57Jh+9OtcAMBdr3/Kmi3FMZdIRKTxUaAQYmajgaeBBcB5wM+As4Cnqvjoc8C3gN8DXwYeAm4DfhqZr3uw7DGR189qZw1EpCG47eWP2V5SBsAvzhhKfjOdamtbfrMcfn7GEAC2l5Rx2//7OOYSiYg0Ps2qnqVJuQ5YBpznnCsHMLOtwINmdohz7sPoB8xsMHAicJlz7v5g8mtmNggYT+VgoQfwmXPunbpcCRGJz7QF63nxg5UAfOng/Tlm4H4xl6jxGjOwM18cuj+vzFnNCx+s5KvDezGyf6e4iyUi0mjoMVdlI4DXEkFC4LXg3+EpPvMZ0BV4ODJ9G9AmMq07PhARkUaorNzxy5f8mAkt83K56bTBMZeo8bvp9MG0zMsF4JcvfURZuUZsFhGpLQoUKtsfWBuZlvg76VCqzrlS59xq51yJmeWaWVszOxeYgG+KFNYdGGVmS8xst5ktMLNLanMFRCQ+T89czkertgBwxbH96da+Zcwlavy6t2/Jt47tB8BHq7bwzMzlMZdIRKTxUKBQWR5QGp7gnCsL/pvJKEmbg9eTwGPOuVsj76/HBx5XAacBc4AHzOzUVAs0s0IzGxJ+Af0zWhsRqTc7SnZzx7/nAbB/2xZc/gV1YK4v3/xCP7q0bQ7AHa/OY0fJ7phLJCLSOKiPQuYyqc8eA7TD91m4wcw+cs5V1Co45w4Lz2xmrwOzgRuBf6VY5pWk6OxcVFTE6tWr6dGjB3379mXy5MkV740dO5bZs2ezdq2vEBk0aBAFBQXMnDkTgIKCAkaMGEFRURHbt28HYNiwYWzfvp158/zNTmFhIQcffDBvvPHGnhUcM4ZFixaxfLl/ate3b18KCwuZPn06APn5+YwZM4aZM2eyceNGAIYOHQrAnDk+p3yHDh0YNmwYkydPpqTEj1h79NFHs3btWhYtWgSgddI6Zd06TdvakbVbdwHwpZ5l7Nq+lQ1rs3udsul3+lLPMibOhTVbdnHDpDf4cv9mWb9OjfF30jo1zXUqKipCspM5p/acCWZWAvzBOfej0LR8YBdwo3Pu9mos66/A14EOzrmUj7fM7C7gIudc0h54ZlYIdI5M7g88P2fOHIYMGZJpkUSkjqzeXMzxd7zJztIyhnZvywtXHUNOjgZXq0/l5Y5xf3mHuSu30DIvlzevP44ubVvEXSwRAebOnZsIYIY65+bGXR7JnJoeVfY5vp9CWJfg31XJPmBmA83sCjNrH3nrQ6A1wU1+MN94M8uNzJdLmtoK59xa59zc8AufYlVEGog7Xp3HzlLfSvEnXxqsICEGOTnGT047CICdpWUVzcBERKTmFChUNg04xczyQtNOCf5NVW/WDbgXiPYzOBSf+Whd8Pd+wCTguMQMZmbACfi+CiKSheas2MzTQQfakwZ3UXrOGI3qvx8nHuSf7Tw1czlzV26OuUQiItlNgUJld+FrAF4wszPM7ArgDuBl59w8ADPLMbPBoWBiMj6I+IuZXW1mJ5vZL4BL8c2YEp2jpwNv4sdkuNTMTsN3eh4MRDs9i0gWcM5x2//7GOegWY5x4xcPjLtITd6NXzqQZjmGcwS/jZrXiojUlAKFkGAgtLPx4yI8ib+Bfwr4Wmi24cBUYFzwmXLgZGAi8D3gWfzozJU6IQfznY0fxfkW4BlgIHCBc+71OlwtEakjb3y8lmkL1wNw4Yje9OvcOuYSSf/OrblwRG8Api5Yz38+iWa8FhGRTKkzcxYKUqTOUWdmkfiUlpVzyl1vs/Dz7bRt0Yy3rj+eDgWZZFGWurZxewnH/u6/bCneTb/OBbz6nS/QLFfPxUTios7M2UtnThGRGnj6/eUs/NynBLx27EAFCQ1Ih4J8rjlhIAALP99e0YdERESqR4GCiEg1FZeW8ac3PgOgW7sWXDSyd8wlkqiLRvamazufHvWPr3/Grt1lVXxCRESiFCiIiFTTP6cvZeXmYgCuO3EgzZtFsx5L3Frk5XLdWF+rsHJzMf+cvjTmEomIZB8FCiIi1bCjZDf3vDkfgL77FXDOsB4xl0hSOeeIHvTp1AqAu/87nx0lKce+FBGRJBQoiIhUw8Qpi1m3rQSA75w4UJ1kG7C83By+e9IBAKzbVsKkqYvjLZCISJbRFU5EJEObd5byt7f8wOgH7t+GcYd0i7lEUpVxh3RjUJc2APztrYVs3llaxSdERCRBgYKISIb+8fZCthT75ivfP3kQOTkWc4mkKjk5xvdP9rUKm3eWct/khTGXSEQkeyhQEBHJwLptu3hgyiIADu3ZnhMPKoy5RJKpkwZ34dCe7QF44J1FrN+2K+YSiYhkBwUKIiIZuPfNBewo8Sk2rz95EGaqTcgWZsYPglqF7SVl3PvmgphLJCKSHRQoiIhUYdXmnTxctASAEf06MnpAp5hLJNV1zID9OLpvRwAeKlrCqs07Yy6RiEjDp0BBRKQKf/7PfEp2lwNw/SmqTchGZsb1pwwCoGR3OX/5z/yYSyQi0vApUBARSWPlpp08+d4yAI4f1JkjeneMuURSU0f26cjxgzoD8MR7y1i5SbUKIiLpKFAQEUnjr28toLTMAXDdiQfEXBrZV9cGozWXlrmKVLciIpKcAgURkRTWbCnmsRm+NuELB3TmsCBzjmSvw3t1YMzA/QD4vxnLWLulOOYSiYg0XAoURERS+NtbCyv6Jlx7woCYSyO1JVGrULK7nL+9rXEVRERSUaAgIpLEum27+Oe7PtPRqP6dOLKP+iY0Fkf16cjIfj5z1aPTl7BO4yqIiCSlQEFEJIl/TF5IcamvTbjmhIExl0Zq2zVjfQ1RcWk5901eFHNpREQaJgUKIiIRG7aX8PA0X5twVJ8OjOin2oTGZmS/ThzZuwMAD01bzMbtJfEWSESkAVKgICIS8cA7iypGYb527ECNm9AImVlFX4UdJWU8MEW1CiIiUQoURERCNu8oZdLUxQAc1rM9xwzYL94CSZ0ZM3A/Dg0yWU2aspjNO0tjLpGISMOiQEFEJGTi1EVs27UbgGvHDlBtQiNmZhXZrLbu2s2kKYvjLZCISAOjQEFEJLC1uJQH3vFNUIZ2b8vxgwpjLpHUtRMOLGRIt7YA3P/OQrYWq1ZBRCRBgYKISOCRoqVsKfa1CdecoL4JTYGZVWS12lK8m0enL425RCIiDYcCBRERoLh0T4fWA7q05qSDusRcIqkvJw/uwsDC1gDc/84iikvLYi6RiEjDoEBBRAR4ZuYKPt/qB9664tj+5OSoNqGpyMkxrji2PwCfb93Fs/9bEXOJREQaBgUKItLklZU7/v72AgC6t2/JuEO7xVwiqW9nHNaNbu1aAPC3txZQVu5iLpGISPwUKIhIk/evOatZvH4HAJeN6Uterk6NTU1ebg6XjekHwOL1O/j33NUxl0hEJH66GopIk+ac49635gPQoVUe5x/VM+YSSVwuGN6T9q3yALj3zQU4p1oFEWnaFCiISJM2Zf565qzYAsD4UX1old8s5hJJXFrlN2P8yD4AzF6xmakL1sdbIBGRmClQEJEmLVGb0DIvt+ImUZqu8aP60DIvF/C1CiIiTZkCBRFpsj5cvokp8/1T4wuG96RDQX7MJZK4dSzIr2h+9s78dcxevjnmEomIxEeBgog0WX99yz8xbpZjFR1ZRS4b05dmQXrcxD4iItIUKVCIMLMzzWyWmRWb2Tozm2RmHav4TIGZ3Wpm881sh5l9bGbXm1mzyHxDzOzfZrYteL1mZofU7RqJSDKL1m3nlTk+s80Zh3Wje/uWMZdIGooeHVpxRpAi95U5q1i0bnvMJRIRiYcChRAzGw08DSwAzgN+BpwFPFXFR58DvgX8Hvgy8BBwG/DT0LLbAW8A3YBLgcuC/79RVSAiIrXv728vIJHUJjHYlkjCt4J9otzB399eGHNpRETioUChsuuAZcB5zrkXnXN3A9cAx6d68m9mg4ETgRucc/c4515zzv0a+CcwPjTr14FC4Ezn3BPOuceA04BOwIV1t0oiErV2azFPv+9H3z3xoEIO6NIm5hJJQzNo/zaMPbAQgKffX14xareISFOiQKGyEcBrzrny0LTXgn+Hp/jMZ0BX4OHI9G1A+O5jBLDAOVfR4NU5tzj4fKpli0gdeHjaEkrK/GH+LdUmSAqJfaOkrJyHpy2OtSwiInFQoFDZ/sDayLTE34XJPuCcK3XOrXbOlZhZrpm1NbNzgQn4pkjplp1YftJli0jt21lSxiNFSwA4tGd7juzdIeYSSUN1VJ8OHNqjHQAPFy1hZ0lZzCUSEalfChQqywNKwxOcc4krQyZ5EzcHryeBx5xzt4bey48uO1CWbtlmVhh0gq54AXoEKlJDT81czsYd/lC8fExfzCzmEklDZbYnG9bGHaU8PXN5zCUSEalfGoI0cy6DecYA7Qj6LJjZR86531fxmaqWfSW+U/VeioqKWL16NT169KBv375Mnjy54r2xY8cye/Zs1q71lRiDBg2ioKCAmTNnAlBQUMCIESMoKipi+3af0WPYsGFs376defPmAVBYWMjBBx/MG2+8sWcFx4xh0aJFLF/uL5h9+/alsLCQ6dOnA5Cfn8+YMWOYOXMmGzduBGDo0KEAzJkzB4AOHTowbNgwJk+eTElJCQBHH300a9euZdGiRQBaJ61TnaxTu/bteeCdrQB0agH5az5i27a2Wb1OjfF3akjr1Lzc0b19C1ZsKubPr86ly7bP6N+vX1avU2P8nbRODXudioqKkOxkzmVy/9s0mFkJ8Afn3I9C0/KBXcCNzrnbq7Gsv+I7MHdwzu02s1eBds65oyPzTQW2OOdOTbGcQqBzZHJ/4Pk5c+YwZMiQTIsk0uS9Onc133z4fQBuPn0w3zimb8wlkmxw/zuL+OVLHwHwj4uP5KTBXWIukUh2mTt3biKAGeqcmxt3eSRzanpU2ef4vgRhiSvCqmQfMLOBZnaFmbWPvPUh0Jo9N/nJlp1YftJlAzjn1jrn5oZf+PStIlJN9032T8/atGhWMfquSFXOP6onbVr4Cvh/TFaqVBFpOhQoVDYNOMXM8kLTTgn+TVVv1g24F4jWCByKz3y0LrTsXuE0q2Y2EOiXZtkiUks+WLaJdxdvAOBrR/eidXO1vJTMtG7ejK8N7wXAu4s28MGyTTGXSESkfihQqOwufA3AC2Z2hpldAdwBvOycmwdgZjlmNjgUTEzG3+j/xcyuNrOTzewX+EHV/uCcS3RgfgRYAzxtZl81s/OAZ4DVwP/V2xqKNFGJJ8HNcowJo/rEWxjJOhNG96FZju/4rloFEWkqFCiEOOfeAc7Gj4vwJHArflTmr4VmGw5MBcYFnykHTgYmAt8DnsWPzlypE7JzbhMwFlgI3BfMvxIY65zbUpfrJdLULduwg5dn+xZ+4w7tRtd2LWMukWSbru1acvohXQF4Zc5qlm/cEXOJRETqnureI5xzzwPPp3m/CGgfmbYVuD54pVv2XPY0ZRKRejJxymLKg7wNl41RB2apmcvG9OO5WSspK3dMnLKYm08fHHeRRETqlGoURKRR27yzlMdnLAVgVP9ODOnWLuYSSbYa2r0do/p3AuDxGcvYUpxsaBwRkcZDgYKINGqPz1jK9mBE3cuDwbNEaiqxD23btZvH3l0ac2lEROqWAgURabTKyh0PTl0CwIDC1hx7QHRIEpHqOfaAzgwobA3Ag1OXsLusPOYSiYjUHQUKItJovf7xGlZs2gnAhFF9yAmy1ojUVE6OccnoPgCs2LST1z5aE2+BRETqkAIFEWm0Jk1ZDPgB1s46vHu8hZFG4+zDe9Cupc+Q/cCURTGXRkSk7ihQEJFG6ZPVW5i2cD0A5x/ZkwINsCa1pGV+Ll8NBmCbsXgjs5dvjrlEIiJ1Q4GCiDRKib4JZnDxyD7xFkYanYtH9iY3aMo2UbUKItJIKVAQkUZn044Snv3fcgDGHtiFXp1axVwiaWy6tW/JqUP3B+DFD1eydmtxzCUSEal9ChREpNF5fMYyikt9NpoJo/rEWxhptC4NOjWXljkeKVKqVBFpfBQoiEijUlbueGjanpSoowd0irlE0lgN69WBQ3v4Afz+OX0JxaVlMZdIRKR2KVAQkUYlmhLVTClRpW6YGZeM7gvAum0lvPjByphLJCJSuxQoiEijEk6JevYwpUSVuvWlg7tS2KY5ABOnLMY5F3OJRERqjwIFEWk0oilRW+UrJarUrfxmOVw0ojcAH63awvRFG2IukYhI7VGgICKNxoNTFwNKiSr162tH9yK/mb+cPvCOUqWKSOOhQEFEGgWfEnUFoJSoUr86tW7OmYd1A+C1j9ewdP2OmEskIlI7FCiISKMQTol6SZC2UqS+JDo1OwcPTlsca1lERGqLAgURyXq7y8orUqIOLGzNqP5KiSr166CubRnZz+93T8xYxrZdu2MukYjIvlOgICJZ7/WP11akRB2vlKgSk0uP8bUKW3ft5qn3lsVcGhGRfadAQUSyXqITs1KiSpxOOLCQXh1935hJUxdTXq5UqSKS3RQoiEhWC6dEveAopUSV+OTmGBNG9QFg8fod/Hfe2ngLJCKyjxQoiEhWC6dEvWhEn1jLInLekT1o3dwHqw9MUapUEcluChREJGspJao0NG1a5HHekT0AmDJ/PZ+u2RpziUREak6BgohkLaVElYZo/Mg+JPrTJ2q8RESykQIFEclKSokqDVWf/Qo4flAhAM/MXMHmHaUxl0hEpGYUKIhIVlJKVGnIEp2ad5aW8YRSpYpIllKgICJZadJU31FUKVGlIRozcD/6dy4A/EjNZUqVKiJZSIGCiGSdj1dtoWjhBkApUaVhMtuTKnX5xp288fGaeAskIlIDChREJOs8NG0x4FOiXjyyT5xFEUnp7GE9aBOkSp2kTs0ikoUUKIhIVtm4fU9K1BMP6kLPjkqJKg1TQfNmnHdkTwCmLljPvNVKlSoi2UWBgohklcff25MSNdG0Q6Shunhk7z2pUqctjrMoIiLVpkBBRLLG7rJyHlZKVMkiffYr4ISKVKnLlSpVRLKKAgURyRrhlKgTRislqmSHCcFggMWl5Tz+3tJ4CyMiUg0KFEQkayRSorZt0YyzDldKVMkOxwzYjwGFrQF4cOoSpUoVkayhQCHCzM40s1lmVmxm68xskpl1rOIzLczs12a2OPjcZ2Z2i5m1iMz3QzNzSV6P1e1aiWS/cErU85USVbKImTE+6E+zYtNOXleqVBHJEgoUQsxsNPA0sAA4D/gZcBbwVBUffQK4ErgbOAN4BPgJ8KfIfN2DZY+JvH5WO2sg0ng9GKSXVEpUyUZnH96dNi2CVKlTFsdbGBGRDOmRXGXXAcuA85xz5QBmthV40MwOcc59GP2AmQ0DxgGXO+fuCya/amadgcvN7CrnXKL3Wg/gM+fcO3W+JiKNyMbtJTw3SylRJXsVNG/G+Uf25L53FjFt4Xo+Wb2FA/dvG3exRETSUo1CZSOA1xJBQuC14N/hKT7THHgU+Fdk+kdAPtA6NK07PhARkWoIp0S9RClRJUtdPLLPnlSpGoBNRLKAAoXK9gfWRqYl/i5M9gHn3DTn3IXOueWRt44HZjnnNoamdQdGmdkSM9ttZgvM7JJ0BTKzQjMbEn4B/TNfJZHsFk6JekCX1oxUSlTJUr06tWLsgf5S8uz/VrBpR0nMJRIRSU9NjyrLAyoluXbOlQUpGPMzXYiZjQPOBU6PvLUe2ADcEHzPlcADZrbKORetkUi4khR9GIqKili9ejU9evSgb9++TJ48ueK9sWPHMnv2bNau9XHOoEGDKCgoYObMmQAUFBQwYsQIioqK2L59OwDDhg1j+/btzJs3D4DCwkIOPvhg3njjjYrljhkzhkWLFrF8uY+L+vbtS2FhIdOnTwcgPz+fMWPGMHPmTDZu9DHS0KFDAZgzZw4AHTp0YNiwYUyePJmSEn+hPProo1m7di2LFvmsNlonrVNind5fU8aKTf6wPLVfS/7zn/9k/To1xt9J65TZOh3Rdhuv41OlPjRlAaf3b57169QYfyetU+2uU1FREZKdzDmlaUswMwf83Dl3S5Lptzjnfp7BMg4FJgOTnHPXVjFvLjAb+Nw5d2yKeQqBzpHJ/YHn58yZw5AhQ6oqkkhWu+Dv0yhauIG2LZpR9OOxynYkWc05x0l/eJv5a7fRvX1L3rr+OJrlqnJfGre5c+cmApihzrm5cZdHMqezU2WlQMvwBDNL1CTsqurDZtYX31fhTeC7Vc3vnCsDXgWGpplnrXNubviFz5wk0uiFU6JeMLyXggTJembGhEqpUqOtXUVEGg4FCpV9ju+nENYl+HdVug8GT/5fBeYBXwmCgPD7A81sfFCLEJYLqFpHJIlwStSLRvSOtzAiteTsYaFUqcEggiIiDZEChcqmAaeYWV5o2inBvykb2JlZa+BlYDMwzjlXnGS2/YBJwHGhzxlwAjBnn0ot0ght3F7Cs/9TSlRpfFrlN+OCo3oCULRwAx+v2hJziUREklOgUNld+P4AL5jZGWZ2BXAH8LJzbh6AmeWY2eBEMBH8+ywwELgVONTMjgm9BgbLno5vkvSgmV1qZqcBTwKDg8+JSMjj7y1j126lRJXGSalSRSQbKFAICQZCOxvoir+JvxU/KvPXQrMNB6biB1kDn/L0RKAtPmCYHHn9JFh2ebDs54BbgGfwwcUFzrnX63C1RLKOUqJKY9ezYyvGHuhbtj43awUbtytVqog0POoZGOGcex54Ps37RUD70N+LActw2RuBq4OXiKTw+sdrWLFpJwDjR/XBLKNDTCSrXDK6D69/vIbi0nIef28ZVxyrIXJEpGFRjYKINDgTpywGoG2LZpx1ePd4CyNSR0b178TAwtYAPDxtCbvLymMukYhIZQoURKRB+XjVFqYvUkpUafzMjAmj+wCJVKlr4i2QiEiEAgURaVASHTtzlBJVmoCzDu9O2yBVaqImTUSkoVCgICINRjgl6lilRJUmoFV+My4Y3guA6Ys28NFKpUoVkYZDgYKINBiPzVBKVGl6LhrRmxylShWRBkiBgog0CD4l6mJAKVGlaenZsRUnHqRUqSLS8ChQEJEG4fWP17Bysx/UfMKovkqJKk3KhKAGbdfuch6bsSzewoiIBBQoiEiDEE6Jeubh3eItjEg9G9m/Ewd0SaRKXaxUqSLSIChQEJHYfbRSKVGlaTMzJozqC8DKzcW89pFSpYpI/BQoiEjslBJVBM48vBvtWuYBMFGdmkWkAVCgICKx2ri9hOdm+ZSoJyolqjRhrfKbccFRPQF4d9EG5q7cHHOJRKSpU6AgIrEKp0SdoJSo0sRdqFSpItKAKFAQkdgoJapIZT07tuKkwYlUqSvZoFSpIhIjBQoiEpvXPlJKVJGo8UHNWsnuch6bsTTewohIk6ZAQURiMyloWtGuZZ5SoooERvbrxKAubQB4eNoSpUoVkdgoUBCRWFRKiXpUT6VEFQmYGRNG9wFg1eZiXlWqVBGJiQIFEYlFOCXqhUqJKlLJmeFRC/QAACAASURBVId1r0iVOikYjFBEpL4pUBCRerdBKVFF0mqZn8sFw4NUqYs3MGeFUqWKSP1ToCAi9e6xGUv3pEQNmliISGUXKVWqiMRMgYKI1KvdZeU8Mm0JAIO6tGFkP6VEFUmmR4dWnDx4fwCe/2Al67ftirlEItLUKFAQkXoVTok6flQfpUQVSSNR4+ZTpS6LtzAi0uQoUBCRejUx6JjZtkUzpUQVqcLRfTty4P4+VeojRUsoVapUEalHChREpN7MWbGZdxf7lKhfHd5LKVFFqmBmTAgGYFu1uZhX5ypVqojUHwUKIlJvErUJOQYXjVRKVJFMfPmw7rRvFaRKnboo5tKISFOiQEFE6sW6bbt48YOVAJwyZH96dFBKVJFMtMzP5YKjegEwY/FGpUoVkXqjQEFE6sU/py+lJGhffcnovjGXRiS7XDRyT6rUSUqVKiL1RIGCiNS5kt3lPFzkU6IO6daWo/p0iLlEItmle/uWnDLEp0p9YdZK1ilVqojUAwUKIlLnXp69is+3+hubCUqJKlIjiU7NJWXlPPbu0ngLIyJNggIFEalTzjkmTvEdMDsV5DPuUKVEFamJ4X07clDXtgA8UrRUqVJFpM4pUBCROjVz6SY+WO47X3796F60yMuNuUQi2cmnSvXZwlZvKebfc1fHXCIRaewUKIhInUrUJuTlGheOUEpUkX1RKVVqkG5YRKSuKFAQkTqzavNOXpnjn3qednBXCtu2iLlEItmtRV4uXx3uU6W+t2Qjs5crVaqI1B0FChFmdqaZzTKzYjNbZ2aTzKxjFZ9pYWa/NrPFwec+M7NbzKxFZL4hZvZvM9sWvF4zs0Pqdo1E4vNI0RLKyh2glKgiteXCEb3JDXKlKlWqiNQlBQohZjYaeBpYAJwH/Aw4C3iqio8+AVwJ3A2cATwC/AT4U2jZ7YA3gG7ApcBlwf/fqCoQEclGxaVl/HO6z8wyrFd7Du3ZPuYSiTQOPlVqFwBe/ECpUkWk7jSLuwANzHXAMuA851w5gJltBR40s0Occx9GP2Bmw4BxwOXOufuCya+aWWfgcjO7yjlXCnwdKARGO+cWBJ8tAhYCFxIKKkQag+dnrWDjjlJAtQkitW3CqL68PHs1JWXl/N/0pVwzdmDcRRKRRihraxTM7Doz61LLix0BvJYIEgKvBf8OT/GZ5sCjwL8i0z8C8oHWoWUvSAQJAM65xcBnaZYtkpV8StTFAOzftgWnDt0/3gKJNDJH9enA4ESq1OlLlCpVROpE1gYKwB+AZWb2LzP7upm1qoVl7g+sjUxL/F2Y7APOuWnOuQudc8sjbx0PzHLObUyz7MTyky4bwMwKg74NFS+gf1UrIhKnaQvX88nqrQBcNLI3ebnZfKoRaXjMjAmj+wCwZssu/jVHqVJFpPZlc9OjgcCZ+D4Bk4C/mtnz+P4Br0ZqBTKVB5SGJzjnyoJRZPMzXYiZjQPOBU4PTc6PLjtQVsWyr8T3ldhLUVERq1evpkePHvTt25fJkydXvDd27Fhmz57N2rU+Nhk0aBAFBQXMnDkTgIKC/9/efcdXWd/9H399sggJEAghjLDClgSQiAw1LlSsEwfW2dvWUUfVX4ftrbdWbW3tbe1dta3bVlutAxcOHIhagxAUAkLCDmGEFQIBQiD7+/vjOokhYQRIcuWcvJ+Px3kccp3ruvxcJya53ue7Yhk/fjyZmZmUlJQAkJaWRklJCcuXLwcgMTGRESNGMHPmzNrzpqenk5eXR36+l4uSk5NJTExk7ty53kVGRZGenk5WVhZFRV5GSk1NBSA7OxuALl26kJaWRkZGBuXl5QCMGzeOgoIC8vK8qTR1TcF9TY996PXSiwyDU3pHUFRUFPTXFIrfJ11TcF9TxypHxyijuNzxl48XE1O4NOivKRS/T7qmEWRmZiLByZxzftdw1MwsAe+m/ALgTGA38CrwL+dc1mGcxwEPOOfu38/2+51zDzTiHKOADOAF59ztdbZ/AeCcO7Xe/l94m91pBzhfItCt3uaBwLTs7GxSUlIOVZJIi1q3bQ+nPPI5zsH3x/Thfy/VxF4izeXhj5bxxBdej9Z3f3IiI3tr0gBpfXJycmoCTKpzLsfveqTxQqU/wDZgDbAZ2AMk4AWGeWb2mZn1bOR5KoD2dTeYWc2n/YecVsLMkvHGKnwB/LTey+X1zx0QdbBzO+cKnHM5dR94szKJtEr/nLOGms8ffnhSfz9LEQl5mipVRJpT0AYF85xmZk8Am/CmHj0JeBjo7ZxLBY7F6///YiNPuxVvLEFdNQOmNx2inkTgE2A5cJlzrqoR5645/0HPLRIsdpdV8tq89QBMGNCVYT06+VyRSGjr1bk9Z6d4f1re/3YTW4s1VaqINJ2gDQp4rQef4o0FeA043jk30jn3J+fcFoDAdKaP4AWIxpgDTDKzyDrbJgWeD9jBzsw6ANOBncD5zrnSA5y7b90F1sxsMDDgYOcWCSavf7Oe4tJKAH4YGGgpIs2rZlBzeVU1r3y9zt9iRCSkBHNQmANcAvRyzt1xkLEInwHXNvKcj+KNB3jXzC4ws5vwgsZ059xyADMLM7PhNWEi8Pw23uDqB4FRZnZSnUfN5NYvAVuAN83sCjObAryFF3heObxLF2l9Kquq+ftX3sC4/l1jOOOYpp69WET2Z0y/LqT0CkyVmrmW8kpNlSoiTSNog4JzbrJz7h3nXOUh9lvnnHu9keecBVwM9ASm4t34vwFcWWe3scBsvEXWAJKAM4BOeIEho97jfwLn3gFMxFtg7TngH8BGYKJzbldj6hNpzT5ZsoX8or0AXHdSMmGBftMi0rzMjGtP6A9AQXEZH2arN6uINI1gnh61WTjnpgHTDvJ6JtC5ztdrgEbdEQUGIk865I4iQejZjNUAdI6J5JLjevtcjUjbcv6oXjz04TK2l5Tz4uw1XHhskt8liUgICNoWBRFpPeavLWLBuh0AXDWuLzFR+gxCpCVFR4Zzxdg+AGSt28G363f4XJGIhAIFBRE5as8FWhMiw43/mtDf32JE2ihNlSoiTU1BQUSOyrpte/g4ZzMAF4xKIrFTtM8VibRNPePac3aqN1Xqe99uZNPOvT5XJCLBTkFBRI7KP2bnUR1YYO369GR/ixFp464/yfsZrKx2alUQkaOmoCAiR2zn3gpe/8ZbYO2kQQkc01MLrIn4aXTfLoztHw/AvzPXUVxa4XNFIhLMFBRE5Ii9+vU6Ssq9RcjVmiDSOtxw8gAAissqeS0Q5EVEjoSCgogckYqq6tquDYMTO3DKkG7+FiQiAEwclsiAhFgA/vHVGiqqtACbiBwZBQUROSLTF29i085SwGtNMNMCayKtQViYcX2616qwYcdepi/WAmwicmQUFETksDnnahdYS+gQpcWdRFqZi9OS6BobBXiLITrnfK5IRIKRgoKIHLa5edvJ3rALgGvG9yc6MtznikSkrujIcH4QWNMke8Mu5qze5m9BIhKUFBRE5LDVLLDWLiKMq8f39bkaEdmfayb0o12E92f+2S9X+1yNiAQjBQUROSyrt+7m06UFAFyc1puuHdr5XJGI7E98bBRTxvQG4PPlW1mxpdjnikQk2CgoiMhheX5WXu2/rztJU6KKtGbXnTSAmnkGaloCRUQaS0FBRBpte0k5b2blA3D6sEQGJXbwuSIROZjkhFjOGt4dgHcWbKRgV6nPFYlIMFFQEJFGeylzLaUV3pzsWmBNJDjcGFiArbyqmhfnrPG1FhEJLgoKItIoe8urahdYS+nViQkDuvpbkIg0ynH94knr2xmAlzLXUVJW6XNFIhIsFBREpFGmzl/P9pJyAG46ZaAWWBMJIjWtCjv3VjB13nqfqxGRYKGgICKHVFlVXbvAWt/4GL6X2sPnikTkcJw5vAf9u8YA8PxXeVRWVftckYgEAwUFETmk6dmbWb99LwA3pCcTEa5fHSLBJDzMuC7da1VYv30vH+ds8bkiEQkG+msvIgflnOOpL3IB6BobxZQxfXyuSESOxKVpvekSEwnAM1/m4pzzuSIRae0UFETkoDJWFrJk0y4Arj2hP9GR4T5XJCJHon1UONdM6A/At/k7yVy93d+CRKTVU1AQkYN6+kuvNSEmKpxrJvTzuRoRORo/mNCPdhHen/4n/5PrczUi0topKIjIAS3O38lXq7YBcPnxfekcE+VzRSJyNBI6tOP7x3vdB79csZXsDTt9rkhEWjMFBRE5oKcCnzhGhBnXaYE1kZBwQ/oAwsO86Y2f/EKtCiJyYAoKIrJfawpL+DB7EwAXjOpFUuf2PlckIk2hT3wMF47qBcD07E3kFZb4XJGItFYKCiKyX89mrKY6MCnKj08Z6G8xItKkbjrV+5l2Dp7WWAUROQAFBRFpoKC4lKnz8wE4fVgiQ3t09LkiEWlKQ7p35IxjugPwZlY+m3eW+lyRiLRGCgoi0sDzGXmUV3ort96k1gSRkHTLad7PdkWV47nAyusiInUpKIjIPopKyvlX5loAxvaPZ2xyvM8ViUhzSOvbhfEDvJ/vf3+9jqKScp8rEpHWRkFBRPbxj9lr2FNeBcBPTh/kczUi0pxuOdX7Gd9TXsU/56z1uRoRaW0UFESkVnFpBS98lQfAyN5xpA9O8LkiEWlO6YMTSE3qBMALs/PYU17pc0Ui0pooKIhIrZcy17Gr1LtRuPW0QZiZzxWJSHMys9pWhaI9Fbzy9XqfKxKR1kRBoR4zm2xmC82s1MwKzewFM2tUJ20zSzOzx81s3gFe/6WZuf08Xm3aqxA5fKUVVTw/yxvQOKR7B84MzIgiIqFtUkoPBiTEAvDMl7mUVVb5XJGItBYKCnWY2YnAm0AuMAW4D7gIeOMQx11uZouA+cDNwIH6ayQFzp1e73FfU9QvcjRe/Xodhbu9wYy3njaIsDC1Joi0BeFhVruuwpZdZUydl+9zRSLSWigo7OsOYD0wxTn3nnPub8BtwGlmNvIgx50BZADjgZcPsl9vYKVzbla9x/KmugCRI1FeWc3TX3qtCf26xnDuiJ4+VyQiLemi0Un07uKtvv7kF7m10yOLSNumoLCv8cAM51zd35AzAs9jD3SQc+5659ytzrm5hzh/El4QEWlV3l6Qz6bAgks3nzKQiHD9ahBpSyLDw2rHKmzYsZe3F6hVQUQUFOrrARTU21bzdWITnD8JOMHM1ppZpZnlmtkPm+C8IkessqqaJ7/IBaBnXDQXp/X2uSIR8cMlxyXRKy4agL99nktllVoVRNo6BYV9RQIVdTc452pGdUU1wfm34QWPW4FzgWzg72Z29oEOMLNEM0up+wC0VK40mQ8Wb2LNtj0A3HjyAKIi9GtBpC1qFxFeO1Zh3fY9TFu40eeKRMRvEX4XEETcUZ/AuWPrfm1mnwKLgbuAjw5w2C0cYLBzZmYmmzdvpnfv3iQnJ5ORkVH72sSJE1m8eDEFBV6DyNChQ4mNjSUrKwuA2NhYxo8fT2ZmJiUlJQCkpaVRUlLC8uXekInExERGjBjBzJkza8+bnp5OXl4e+fles3RycjKJiYnMnev1uoqKiiI9PZ2srCyKiooASE1NBSA7OxuALl26kJaWRkZGBuXl3uDZcePGUVBQQF6eN4e/rqllrumbefP53w+3eNfQPoLT+7WrPU+wXlMofp90TbqmlrqmxOJVxEXBznL46+er6Ou2sGvnjqC+plD8PgXbNWVmZiLByZw76vvfkGFm5cCfnXO/qrMtCigD7nLO/aER53gBONU517+R/81HgWucc10P8Hoi0K3e5oHAtOzsbFJSUhrznxHZr/e+3chtrywA4O5zhnHjyWqsEmnrnp+Vx2/fXwLAY5cfy4XHJvlckQS7nJycmgCT6pzL8bseaTz1MdjXVrxxCnXVTCa/6WhObGaDzey/zCy83kvhHKS1wjlX4JzLqfvAm2JV5KhUVTsen7kSgK6xUVw9vp/PFYlIa3Dl2L4kdPB62/7ls1VUV+sDRZG2SkFhX3OASWYWWWfbpMDz0babJQAvAKfWbDBv2dvT8cYqiLSoDxZvYmXBbgB+fMoAYqLUE1FEoH1UODekDwBgVcFuPsze7HNFIuIXBYV9PYrXzeddM7vAzG4CHgGm16x1YGZhZja8XphojLnAF8CLZvYjMzsXmAoMBx5ssisQaYS6rQkJHdSaICL7unp8P7rEeH/m/vLZSrUqiLRRCgp1OOdmARcDPfFu4h/EW5X5yjq7jQVmA+cf5rmrA+d+B7gfeAsYDFzunPv0aGsXORwfLN7EqprWhJMHqjVBRPYR2y6C6wOtCss2F/PJki0+VyQiflBQqMc5N805d6xzrp1zLiGwmNrOOq9nOuc6O+feOsDx1x5oILNzrsg59xPnXN/A+Uc5515rpksR2a/6rQlXje/rc0Ui0hr9YEI/4tp7rQqPfrpCrQoibZCCgkgb8/6ijWpNEJFD6hgdyY0nf9eqMD37qOb0EJEgpKAg0oaoNUFEDse1J/QnPtabAenPM1ZQpVYFkTZFQUGkDXl/0UZyt3oL49x0iloTROTgYttFcPMp3voquVtLmLZwg88ViUhLUlAQaSMatCaM00xHInJoV4/vR7eO7QB4bOZKKqqqfa5IRFqKgoJIGzFt4YZ9WhPaR9Vf+09EpKH2UeHceqrXqrB22x7eysr3uSIRaSkKCiJtQHllNX/+dAUAiR3bqTVBRA7L5WP70jMuGoDHZ66irLLK54pEpCUoKIi0Aa99s4712/cCcNvEwWpNEJHDEh0Zzm2nDwZgw469vP7Nep8rEpGWoKAgEuL2llfx+GerAOgT357vj+njc0UiEoymjOlNn/j2APz181WUVqhVQSTUKSiIhLgXZq9ha3EZAD87cwhREfqxF5HDFxkexu2BVoUtu8p4ee46nysSkeamOwaRELZzbwVP/ScXgCHdO3DBqCSfKxKRYHbR6CQGJMQC8MTnqygpq/S5IhFpTgoKIiHs2S9Xs3NvBQC/OGso4WHmc0UiEswiwsO44wyvVWFbSTnPZeT5XJGINCcFBZEQtbW4jL9/5f0RP7ZPZ84c3t3nikQkFJw/shfDe3YC4JkvcyncXeZzRSLSXBQURELU3z5fxZ5yb7DhLycNxUytCSJy9MLCjP/+3jAASsqr+GtgsgQRCT0KCiIhKL9oD/8ODDQ8cVBXThiU4HNFIhJK0gcncOKgrgC8PHct67bt8bkiEWkOCgoiIeixT1dSXlUNwJ2ThvlcjYiEGjPjV2d7v1sqqhx/mrHc54pEpDkoKIiEmKWbdvFGVj4AZw7vzrF9OvtckYiEopG9O3PeyJ4ATFu4kewNO32uSESamoKCSIh56MNlOAfhYcavzh7qdzkiEsJ+cdZQIgKzqf3vR8t8rkZEmpqCgkgI+XLFVr5csRWAK8b2YVBiR58rEpFQ1j8hlivH9QUgY2Uhs1YW+lyRiDQlBQWREFFV7fj99KUAxEaFc8fEIT5XJCJtwW2nDyYmKhzwWhWqq53PFYlIU1FQEAkRb2Xls2xzMQA3nzqQbh3b+VyRiLQF3Tq244b0AQAs3rCT9xdv8rkiEWkqCgoiIWBveRWPfOLNOtKjUzTXnTTA54pEpC254eQBdI2NAuDhj5ZRWlHlc0Ui0hQUFERCwHMZq9myy1sd9ednDaF9oBuAiEhL6NAugv93ptfdMb9ob+2q8CIS3BQURILc1uIynvpPLgDH9OzExWm9fa5IRNqiK47vw5DuHQB44vNcCopLfa5IRI6WgoJIkHv00xWUlHvN/HefM4zwwFSFIiItKSI8jHvOHQ7A7rJK/u+TFT5XJCJHS0FBJIitKijm1W/WA3DKkG6kD+7mc0Ui0padPKQbpw9LBOC1eevJ2ahF2ESCmYKCSJByzvHAe0uoqnaEGdx9zjF+lyQiwt3nHENEmOEc/Pb9JTin6VJFgpWCgkiQ+nRpARmBxY0uH9uXoT20uJqI+G9QYgeumdAPgMzV2/lkyRafKxKRI6WgIBKESiuq+O37SwCIax/JL84a6nNFIiLfuWPiYOLaRwLw++lLKavUdKkiwUhBQSQIPT8rj3Xb9wDwszOHEB+Yv1xEpDXoHBPFT88YDMDabXv45+y1PlckIkdCQUEkyGzauZe/frYKgKHdO3LVuL4+VyQi0tBV4/sxsFssAI/PXEnh7jKfKxKRw6WgIBJkHpq+jL2BVU/vu2A4EeH6MRaR1ieyznSpxWWV/OHDZT5XJCKHS3cYIkHk67ztvPvtRgDOHdGTEwYm+FyRiMiBnTYskYmB6VLfmJ/PvDXbfa5IRA6HgoJIkKiqdtz3bg4A0ZFh3HXOMJ8rEhE5tPsvSKFdhHe7cc872VRWVftckYg0loJCPWY22cwWmlmpmRWa2QtmFt/IY9PM7HEzm3eA11PM7GMz2x14zDCzkU17BRKqXvl6HUs37QLg5lMG0btLjM8ViYgcWp/4GG49bRAAyzYX869MDWwWCRYKCnWY2YnAm0AuMAW4D7gIeOMQx11uZouA+cDNQIP+IGYWB8wEegE/Aq4P/HtmY4OItF3bS8p55JPlACR1bs+PTxngc0UiIo1348kD6NfV+3Dj/z5ZQcGuUp8rEpHGUFDY1x3AemCKc+4959zfgNuA0w7xyf8ZQAYwHnj5APtcBSQCk51zrzvnXgXOBboCVzfVBUho+t0HS9mxpwKAe887hujIcJ8rEhFpvOjIcB64IAXwBjb/fvpSnysSkcZQUNjXeGCGc65uB8oZgeexBzrIOXe9c+5W59zcQ5w71zmXW+e4NcDKg51bZE7uNt7Mygdg4rBEJqX08LkiEZHDd+rQRCaldAfgnYUbyVy9zeeKRORQFBT21QMoqLet5uvEZjh3zfkPeG4zSwyMbah9AAOPshYJEmWVVfzP24sBaB8ZzgMXpmBmPlclInJkfn1+CtGR3q3Hr6dlU6GBzSKtWoTfBbQykUBF3Q3OuarAjdnRLn0bVf/cAVWHOPcteGMlGsjMzGTz5s307t2b5ORkMjIyal+bOHEiixcvpqDAyyZDhw4lNjaWrKwsAGJjYxk/fjyZmZmUlJQAkJaWRklJCcuXe33hExMTGTFiBDNnzqw9b3p6Onl5eeTne59wJycnk5iYyNy5XmNKVFQU6enpZGVlUVRUBEBqaioA2dnZAHTp0oW0tDQyMjIoLy8HYNy4cRQUFJCXlwegawpc00frjdWFewE4r7+xPGsOsUF+TaH4fdI16Zp0TY2/pttOH8wfP17Oii27ufuFGdwxKSXorykUv09NeU2ZmZlIcDLnnN81tBpm5oAHnHP372f7/c65BxpxjheAU51z/ett/wLAOXfqfrY759xpBzhfItCt3uaBwLTs7GxSUlIOVZIEqdVbd3P2oxmUV1UzrEdH3rvtJCK1uJqIBLnyymrOfuxLVm8toX1kOJ/89GT6xGsWt1CWk5NTE2BSnXM5ftcjjae7jn1VAO3rbjCzmk/7j3bt+fL65w6IOti5nXMFzrmcug+8WZkkhDnnuOedbMqrqjGD3188QiFBREJCVEQYD072PvXeW1HF3W8vRh9airROuvPY11a8sQR1dQ88b2qGc9ec/2jPLSHm7QUbmJ3rDfS7alxf0vp28bkiEZGmc8LABL4/pg8AGSsLeTNrg88Vicj+KCjsaw4wycwi62ybFHg+2g52c4C+dadZNbPBwIAmOLeEkKKSch78wJs6sFvHdtw5SSswi0joufucY+jWsR0Av31/CVuLj7bhXkSamoLCvh7FGw/wrpldYGY3AY8A051zywHMLMzMhtcLE43xErAFeNPMrjCzKcBbwGbglaa7BAl2v5u+lO0l3oCz+84fTlz7w/1fTUSk9YuLieS3F3pdkHbureD+99R1XaS1UVCowzk3C7gY6AlMBR7EW5X5yjq7jQVmA+cf5rl3ABOB1cBzwD+AjcBE59yuoy5eQsJny7bwxnxvdotTh3bj3BE9fa5IRKT5nJ3ag++ler1yP1i0iRlLtvhckYjUpelR63HOTQOmHeT1TKDzQV6/9iCv5fBdVyaRfezcW8Fdb3lrJnRsF8HvLxqhNRNEJOQ9cGEKX60qZFdpJfe8s5hxA+LpFK2WVJHWQC0KIq3Eb99fwpZdXh/de847hl6d9zdJlohIaEnsGM095w4HYMuuMv7w4TKfKxKRGgoKIq1A3S5HpwzpxmWB2UBERNqCKWN6c+KgrgD8e+465gRmfRMRfykoiPhs5559uxz94RJ1ORKRtsXMeOiikURHerclv5j6LcWlFT5XJSIKCiI+++0H33U5uve84fSMU5cjEWl7+naN4a7vHQPAhh17+c17S3yuSEQUFER8VL/L0ZQxvX2uSETEP9eM70f64AQAps7P55OczT5XJNK2KSiI+ERdjkRE9hUWZjx86Ug6RXuTMt711mIKd2shNhG/KCiI+MA5x11vL1KXIxGRenrGtee3k72F2LaVlPPfby7GOedzVSJtk4KCiA+mzstn+mKvSf2MYxLV5UhEpI4LRvXi3JHegpOfLt3C1Hn5Plck0jYpKIi0sNVbd3P/ezkAJHZsx8OXjlKXIxGROsyM301OJbFjOwAeeC+H9dv3+FyVSNujoCDSgsorq7nj1YXsKa8C4E+XjSI+NsrnqkREWp/OMVE8fOlIAErKq/j5699SVa0uSCItSUFBpAX934wVLN6wE4AbTx5A+uBuPlckItJ6nTo0kavH9wXg6zXbeWzmSp8rEmlbFBREWsjsVYU8/WUuACm9OvGLs4b6XJGISOv3P+cMZ3BiBwD+8tlKZq8q9LkikbZDQUGkBRSVlPPT1xfiHLSPDOfxK0YTFaEfPxGRQ2kfFc7frkojOjIM5+CO1xaytVhTpoq0BN2piDQz5xx3vvHdVKi/Pn84A7t18LkqEZHgMaR7R35zgTdl6tbiMn762kKqNV5BpNkpKIg0s6e/XM2nS7cAMCmlO5cf38fnikREgs+UMb25aHQSALNWFfLEF6t8rkgkOeZZpgAAHq1JREFU9CkoiDSjObnbePijZQD0jY/h4Us0FaqIyJEwM347OZUBCbGANznE13nbfa5KJLQpKIg0ky27SrntlSyqHbSLCOPJq9OIi4n0uywRkaDVoV0Ef70yjaiIMKod3P7KAraXlPtdlkjIUlAQaQYVVdXc+nIWhbu9P2C/nZxKSq84n6sSEQl+w3t14tfnDQdg865S7nh1AZVV1T5XJRKaFBREmsEfPlzGvLVFAFx+fB8uG6NxCSIiTeWqcX05d2RPADJWFvLwx8t9rkgkNCkoiDSxDxZt4vlZeQCkJnXi/gtSfK5IRCS0mBkPXzKSYT06AvDMl6uZtnCDz1WJhB4FBZEmtKqgmF++8S0AnaIjePKq44iODPe5KhGR0BPbLoJnrhlDXHtv7Ncv31hEdmDlexFpGgoKIk1ke0k51704j5LyKgAevfxY+sTH+FyViEjo6ts1hr9dmUaYQVllNT/+13wKd2sxNpGmoqAg0gTKK6u56aX5rN22B4CfnjGE04d197kqEZHQd9LgBO4+5xgANuzYy60vZ1Ghwc0iTUJBQeQoOef4n7cX187nfcGoXtw+cZDPVYmItB3XnZTM5GN7ATA3bzsPvr/E54pEQoOCgshReubL1Uydnw/AsX068/ClI7WomohICzIz/nDJSFKTOgHw4py1/GvOGl9rEgkFCgoiR2HGki38IbDycq+4aJ75gQYvi4j4IToynKevGUNChygA7ns3hxlLtvhclUhwU1AQOUI5G3dyx6sLcA5iosJ5/trjSewY7XdZIiJtVlLn9vz92uNpHxlOtYPbXsni2/U7/C5LJGgpKIgcgY079nL9i/PYU16FGTx++WiO6dnJ77JERNq8kb0785crRhNmUFpRzXUvfsO6wEQTInJ4FBREDlNRSTk/+PvXbNpZCsBd3xvGGcM1w5GISGtxxvDuPHBhKgCFu8u59h9fU1RS7nNVIsFHQUHkMJSUVfLDF75hVcFuAK49oT83pA/wuSoREanvmvH9uOmUgQCsLizhhn/Oo7SiyueqRIKLgoJII9WslbAw0N/1glG9+PV5wzXDkYhIK/XLSUM5f5Q3beq8tUX8v1cXUqk1FkQaTUFBpBGqqx0/n/otGSsLATh5SDcemTKKsDCFBBGR1ioszHhkykjGJccD8FHOZn75xiKqq53PlYkEBwWFesxsspktNLNSMys0sxfMLL4Rx91gZkvNrMzMNprZo2YWXW+fy8zM7eeR2XxXJEfLOccD7+Xw3rcbARjdtzNPXZ1GVIR+fEREWrt2EeE884MxpPTyJpx4a8EG7p2WjXMKCyKHojudOszsROBNIBeYAtwHXAS8cYjjLgeeATIC+/8FuAV4ot6uSUAJkF7vcWOTXYQ0uT9/upIX56wFYHBiB/5x7fHEREX4XJWIiDRWXPtI/nXdOAYndgDg5bnrePCDpQoLIoegu5193QGsB6Y456oBzKwYeNHMRjrnFh3guJ8BGc65mhv+6WYWCdxrZnc657YFtvcG1jvnZjXjNUgTcc7x509X8vjMlYA3P/c/rxtL55gonysTEZHDFR8bxcvXj+Oyp+ewZtsenp+VR0xUOD8/a6jfpYm0WmpR2Nd4YEZNSAiYEXgeu78DzCwKGA18XO+lGXhBbHSdbUl4QURaufohoXundrx0/Th6xrX3uTIRETlSiZ2iefmG8SR19n6X/+WzVTzxxSqfqxJpvRQU9tUDKKi3rebrxAMck4AXCBpzXBIw0MyWm1mlmeWb2a8OVpCZJZpZSt0HMPCQVyJHzDnHn2esqA0JPTpF8+qNE0hOiPW5MhEROVpJndvz8vXjSOzYDoCHP1rO0//J9bkqkdZJXY/2FQlU1N3gnKsKTH95oP4mNdsr6m2vqvc6wA4gBngI2A5cBfzBzHY4554+wPlvwRsr0UBmZiabN2+md+/eJCcnk5GRUfvaxIkTWbx4MQUFXl4ZOnQosbGxZGVlARAbG8v48ePJzMykpKQEgLS0NEpKSli+fDkAiYmJjBgxgpkzZ9aeNz09nby8PPLz8wFITk4mMTGRuXPnehcbFUV6ejpZWVkUFRUBkJrqLXqTnZ0NQJcuXUhLSyMjI4Pycm8BnHHjxlFQUEBeXh6Ab9fknOOLbR14aYHXW6xLO/hFWgTJCbFBe02h+H3SNemadE26pqO9pttSHX9aGMGOvZU89OEyspev5NrjunHccccF7TW11u9TZqbmbAlWpoE83zEzBzzgnLt/P9vvd849sJ9j+gN5wA+dcy/sZ/u1zrkXD/Lf/AQY6pzrd4DXE4Fu9TYPBKZlZ2eTkpJyqMuSRnLO8adPVvDXz71maK8lYTz91ZIgIhKSlm8u5urn57K1uAyA605K5p5zj9H6OE0sJyenJsCkOudy/K5HGk9dj/ZVAezTCT0wBgGg7ADH1KwJX7/z+qGOqzEd6GtmHff3onOuwDmXU/eBNyuTNKHqascD7y2pDQk94xQSRERC3dAeHXn9xxPoFefNZv78rDzufjubKq2zIAIoKNS3FW+cQl3dA8+bDnDMNqD6UMeZWU8zu9bM4urtF36EtUoTKaus4vZXF/DC7DWAFxJeuUEhQUSkLUhOiGXqzSfQv2sMAK98vY6fv64VnEVAQaG+OcCkwNSmNSYFnvfbwc45VwYsAM6r99IkvNaGBXW2/QO4tN5+ZwBrnXPFR1q0HLndZZVc98I83l/k5cDBiR148+YTFBJERNqQpM7tef3HExjS3Vtn4Z2FG7n55SxKK6oOcaRIaFNQ2NejeOMB3jWzC8zsJuARYLpzbjmAmYWZ2fB6YeJPQJqZvWRm55rZncCdwHPOuV0AzrlNwD+BR8zsdjM728yeBM4GHmy5S5QahbvLuOKZTGatKgTguH5dmHrTBHp11hSoIiJtTWKnaF67cQIjkryG/xlLtnD5M5kU7j5UD2KR0KWgUEdgIbSLgZ7AVLwb+DeAK+vsNhaYDZxf57hX8FZXHgO8DfwUb3Xmn9b7T9wE/BVvYbdpwETgdufcc81wOXIQ67fv4dInZ7N4w04ATh+WyEvXjdNiaiIibViX2ChevmEcEwZ0BWDh+h1c9MRXrCpQo7+0TZr1KAgF1lLI1qxHR+abNdu56V/z2VbijUO/9LjePHTxCCLDlZtFRATKK6u5663FvJnlTUvaKTqCp68Zw4SBXX2uLDhp1qPgpTsjaVNe/XodVz6bWRsSbjplIH+8dKRCgoiI1IqKCOORKSP52ZlDANhVWskP/j6XN+fn+1yZSMvSgmvSJlRUVfO7D5bWzmwUGW78bvIILju+j7+FiYhIq2Rm3D5xMH3jY/jlG4sor6rm51O/JXfrbn5+1lDCw7TWgoQ+BQUJeUUl5dz67yxm53qrLSd0iOKpq49jTP94nysTEZHWbvLoJHrGRfPjl+azY08FT3yRy6L8nTx+xWjiYzWuTUKb+ltISFu6aReTn/iqNiQM79mJaT85SSFBREQabdyArrx9y4m106fOWlXIeY9nsHD9Dp8rE2leCgoSkpxzvJS5lgv/9hVrt+0B4NwRPXnj5gkkafpTERE5TMkJsbx9y4lcMKoXABt3lnLZU3N4KXMtmhhGQpWCgoScnXsruPXfWdzzTjblldWEGfzirCH89crRxESpt52IiByZ2HYRPHb5sdx3/nAiwozyqmrueSebX0xdRElZpd/liTQ53TVJSFm4fgc/+XcW+UV7AejeqR2PXz6acQM0pZ2IiBw9M+OHJyYzIimOW17OoqC4jDez8pm/djt//v6xjO7bxe8SRZqMWhQkJFRVO57+Ty6XPjm7NiScPiyRD+84WSFBRESa3Jj+8bx/+0mMH+CNeVuzbQ+XPjWHRz9dQWVVtc/ViTQNBQUJerlbd3PZ03N46MNlVFY7IsONe849huf/a4xmpBARkWaT2DGal68fz13fG0ZkuFFV7Xj005Vc+tQc1hSW+F2eyFFTUJCgVVXtePbL1ZzzWAbz1xYB0K9rDFNvOoHr0wdgpjmuRUSkeYWHGT8+ZSDv3PrdrEgL1+/gnMczeClzLdXVGugswUtBQYLSqoLdXPrUbH43fSllldWYwQ9P7M+Hd6RzbJ/OfpcnIiJtTEqvON79yUn88MT+AOwpr+Ked7K57Ok5rNhS7G9xIkdIQUGCSmlFFX/9bCXnPJ7BgnXe/NXJCbG8/uMJ3Hd+imY1EhER30RHhnPf+Sn867qx9IqLBmDe2iLOeSyDP368jNKKKp8rFDk8CgoSNGYu3cJZf/6SRz5ZQXmgFeH6k5KZfns6x2sBNRERaSXSB3djxs9O4bqTkgkzqKx2/O3zXCY9+iWzVhb6XZ5Io+njV2n18gpL+M17OXy+fGvttmE9OvK7i1I5rp8CgoiItD6x7SK497zhXDQ6ibveWsziDTtZu20PVz8/l7NTevDf3xtG/4RYv8sUOSgFBWm1dpVW8NQXuTyXkUd5YKq5jtER/PzMIVw9vh8R4WoQExGR1i01KY63bzmBF+es5U+fLGdPeRUf5Wxm5rIt/GBCf247fRCdYzRDn7ROCgrS6pRWVPHPOWt44otcduypqN3+/TF9uPPsoSR0aOdfcSIiIocpIjyM605K5pwRPfjjR8t5a8EGKqocz8/K4435+dw+cTDXjO9HVIQ+AJPWxZzTtF3BxsxSgOzs7GxSUlL8LqfJVFRV8/q89Tw+cyVbdpXVbh/VO44HLkzVbEYiIhISFufv5MEPljA3b3vttr7xMfzktEFclJZEZIi1mOfk5JCamgqQ6pzL8bseaTwFhSAUakGhoqqa9xdt5LFPV7Jm257a7QO7xfKLs4ZydmoPrYkgIiIhxTnHjCVbeOjDZeTVWZytT3x7bj11EBen9Q6ZFgYFheCloBCEQiUo7C2v4rVv1vFsRh4bduyt3Z7UuT13nDGYi0cnaRyCiIiEtIqqal79eh1/+zyXzbtKa7cndW7PLacN5NLjetMuItzHCo+egkLwUlAIQsEeFIpKyvnnnLW8OGcN20vKa7cndIjillMHcdX4vkH/S1FERORwlFVW8fq8fJ78fBUbd34XGBI6tOOa8f24anzfoB2jp6AQvBQUglCwBoUlG3fx0ty1vLNgA3vKv1t0pm98DDeePIBLj+tNdKQCgoiItF1llVW8MT+fJz7P3ae1PSoijMnH9uJHJyUzrEcnHys8fAoKwUuzHkmzKq2o4sPsTfxrzlqyAisp10jp1YmbThnI91J7qIuRiIgI0C4inKvG9WPKcX14f9FGnp+VR87GXZRXVvP6vHxen5fPCQO7cvnYvpw1vLs+YJNmpaAgzWLppl28lZXPG/PzKaozxakZnDY0kWtP6E/64AQNUhYREdmPqIgwLk7rzUWjk/g6bzt//yqPT5ZswTmYnbuN2bnbiGsfyUWjk5gypjcpveL8LllCkIKCNJlNO/fy7sKNvL1gA8s2F+/zWtfYKL5/fB+uGNuXPvExPlUoIiISXMyMcQO6Mm5AV9Zt28OLc9bwZlY+O/ZUsHNvBS/MXsMLs9eQmtSJS9J6c86InnTvFO132RIiNEYhCLWmMQoFxaV8uqSADxZvZHbuNur/7zQ2OZ6rxvXl7NQeGqAsIiLSBEorqpixZAuvz1vPrFWF+/ztNYPj+8dz/sienJ3ak24d/R8ArTEKwUtBIQj5HRTWbivh45zNfJyzhax1RQ3CwYBusVw8OokLj01S64GIiEgzyi/awxvzva6++UV793ktzGBcclfOHN6d04cl0j8h1pcaFRSCl4JCEGrpoFBaUcU3a7aTsbKQ/yzfyvItxQ32SejQjvNG9uTitCRGJMVp7IGIiEgLcs7xbf5O3v92Ix8s3sSmOlOs1hjQLZbThyZy+rBExvSPb7EF3RQUgpeCQhBq7qBQVe1YvrmY2bmFfLmykLmrt1FWWd1gv35dY5iU0oNJKd0Z3acLYWEKByIiIn6rrnYsWL+D9xdt5JOcLftMs1ojJiqcscnxnDCwKxMGJDC8VyfCm+nvuIJC8FJQCEJNHRTKK6tZvGEHX+cV8XXeNuatLaK4tLLBfmEGo/p05tQhiUxK7c7Q7h3VciAiItKKOedYsWU3M5dt4fNlBcxfW0T1fm79OkVHMH5AV8b070Ja3y6kJsU12dSrCgrBS7MetTHV1Y7VhSUsyt/BovydfJu/gyUbd+23xQC8JeRPHpJA+uBunDgwgbiYyBauWERERI6UmTG0R0eG9ujILacOoqiknC9XbuWrVYXMzt1WO65hV2klnyzZwidLtgAQGW4M7xXH6D6dGd23M6lJcSR3jVXvgTZGQSGE7S6rZPnmYpZvLmbZ5l0s21zMko272F3WsLWgRp/49hzfP56x/eMZmxxPckKsWg1ERERCRJfYKC481ptwBGD99j3Myd3G7NxCMldvZ/Mub2xDRZXj2/U7+Hb9Dl6Y7R0bExXOMT07kdLLewzt0YkB3WLpFK0PEUOVgkKQq6p2bNyxl7zCktrHmm0l5G7dzfrtDfsk1hUVEcbwnp0Y1TuOtH5dGJscT8+49i1UuYiIiPitT3wMfeJjuOz4PgBs3LGXrHVFLFi3g6x1ReRs2EV5ldfrYE95FfPXFjF/bdE+50js2I5BiR0Y2K3DPs/dO7XTh41BTkGhDjObDNwPDAN2A+8DP3PObT/EcTcAPwMGANuA14H/ds6V1tmnD/Bn4AygHbAwsM9/jrTe6178hu2R62p/gA+mfWQ4Q7p34JienRjZuzMje8cxtEdHIsNbZsYDERERaf16dW5Pr87tOW9kL8Cb+XDFlmJyNu4iZ+NOcjbuYummXZRWfHfvUVBcRkFxGbNzt+1zrtiocPrEx9Bx76YWvQZpOgoKAWZ2IvAm8A5wL9AX+H3g+fSDHHc58AzwLPBzYBTwANAJ+FFgnwjgQyAOuB0vhPwU+NjMRjrnVhxJzeu27SGq274hIToyjP5dY0lOiGVI944c07MjQ3t0om98TLPNZiAiIiKhKToyPPABY+fabVXVjrzC3awq8B65W0sCz7vZU15Vu19JeRXLNhdTvnXb/k4tQUBB4Tt3AOuBKc65agAzKwZeDNzMLzrAcT8DMpxzNwa+nm5mkcC9Znanc24bcBaQAkx0zn0WOPcMIB+4KXCOwzZhYDxpo5LpnxBLctdYkrvF0r1jtAYaiYiISLMJDzMGJXZkUGLHfbZXVzs27yold6sXIPIKS1i3fQ9LKwtQm0JwUlD4znjg45qQEDAj8DwWaBAUzCwKGI3XXamuGXitCqOBTwPnLgc+r9nBOVdsZpmBcx+R31w4gpSU4Ud6uIiIiEiTCQuz2q5L6YO71W7PyYkl9U8+FiZHTB3Uv9MDKKi3rebrxAMck4AXtg51XA+g0DVctKLgIOcWEREREfGNWhS+EwlU1N3gnKsKjNaPOsAxNdsr6m2vqvd61H72qdnvQOcGwMwSgW71Ng8DWLVq1cEOFREREfFdnfuVg97zSOujoNA4R7p8dWOOO9Q+twD37e+FyZMnH3ZBIiIiIj5JBRb4XYQ0noLCdyqAfRYRCIxBACg7wDHlgef6iw/UP658P/vU7Hegc9d4Aphab9sI4BXgUmDZIY6XpjEQmAZcCOT6XEtbofe85ek9b3l6z1ue3vOWNwx4AziiWR7FPwoK39mKN5agru6B5wMN1t8GVDfiuK1AVzOLdM5V1NvvoBMBOOcKqDcGos7iJcucczkHO16aRp33PFfvecvQe97y9J63PL3nLU/vecur857v9rMOOXwazPydOcCkwNSmNSYFnjP3d4BzrgyvCe28ei9NwmtFqGlemwOEA9+r2cHMOgETDnRuERERERE/KSh851G8QcPvmtkFZnYT8Agw3Tm3HMDMwsxseL0w8ScgzcxeMrNzzexO4E7gOefcrsA+HwE5wHNmdl1gBej38MYnPN0ylyciIiIi0ngKCgHOuVnAxUBPvDEBD+L1p7uyzm5jgdnA+XWOewW4ERgDvI234vJfAs81+1QCZwP/Af4Pb3xBBHCWc25Nc12TiIiIiMiR0hiFOpxz0/AGOB3o9Uyg8362Pws8e4hz5wNTjrbGgK14C7ptbaLzyaHpPW95es9bnt7zlqf3vOXpPW95es+DlDVcA0xERERERNo6dT0SEREREZEGFBRERERERKQBBQUREREREWlAQUFERERERBpQUAhiZpZiZtPMrMDMtpvZ52Z2kt91hbLAe/6ume0ws11mNtvMJh36SDkaZtbOzC4zs+lmdo/f9YQaM5tsZgvNrNTMCs3sBTOL97uutsDMhpjZ781svd+1hLrAWkh3mtkKMyszszVm9piZNZjNUJqGmUWY2c/MbKmZ7TGzXDN7yMxi/a5NGkdBIUiZ2VBgLtAJbx2HH+AtjT7DzFL8rC1UmdkA4Cu8tTZuwltjoxiYbmbpftYWqswsysyeADYBL+Otbq5pnZuQmZ0IvAnk4k3hfB9wEd46MtJMzOxMM5sNLMdbpLO3zyW1BY/irZH0GnAB8Bjwo8DX0jyewnvP/wFciPc9uBl40s+ipPH0Bzd4XYf3/TvPOVcCYGafAkXAZXh/7KVp/T+8cD3JObcdwMxmAuvw/thk+FhbqGoPnID3h+ZlYLO/5YSkO4D1wBTnXDWAmRUDL5rZSOfcIl+rC10nAyuBewP//rW/5YQ2M+sG3Ao85Jy7N7D5YzMLB/5oZv2cc2v9qzD0mFlH4Fq89/zhwOYZZtYVuNvMrnfOlftWoDSKgkLwehD4Y01ICCgHKoCO/pQU8vKBx2pCAoBzbq+Z5QHqptEMnHM7gWNrvjYzH6sJWeOBj2tCQsCMwPNYQEGhGdS5WUUtki0iDngFr/WsriWB53hAQaFpleC1lO2st3033v1ne7z7FmnFFBSClHNuF4B5d07heN1h7gGMQ6wSLUemzicitcwsAUjF+7RbJBj1AArqbav5OrGFaxFpFs65VcDV+3npNGALkNOyFYW+wIcPmwECLTfReC3EvwD+HvggSFo5BYVWysx64n0Csj+r6zTXXQJMDfy7DK8r0tLmri8UHcZ7XteTeH9knmu2wkLYEb7n0rQi8VoiaznnqgKtN1G+VCTSAszsOOB24Kf6XdPsvgVqxk9+hDdOQYKAgkLr9RDwXwd47RhgWeDfM4HjgV54fehfNbMJzrmVzV9iyGnsew6Amf0GOB84vV4XMGm8w3rPpcU5vwsQaQ5m1geYhnfTqoG1ze/7QBe8FoVfA3/CC2nSymnWo1bKOXetc84O8FhWZ78i59w859y7eDetBvzMt8KDWGPfcwAzuxn4H+Ba59xsfyoOfofznkuzqcDrK1zLzGpaEspavhyR5hWY+vdjYCNwpXNOgbiZOedynHOzAl14HwRuNbNeftclh6agEKTM7Dwzu6jutsCn2rnAAH+qahvM7BLgL8CtzrlX/a5H5ChtxRunUFf3wPOmFq5FpFmZWQzwAV5r2ffUGtx8zKyHmd0UaL2paxHe/Wc/H8qSw6SgELwuAZ6tu2hJ4N+DgDV+FRXqzOxUvIHL9zrnnvK5HJGmMAeYZGaRdbbVLCKY6UM9Is3CzCLwxvT1AM50zm3zuaRQ1w6vW9dV9baPCjxrlqkgoDEKwetx4ArgAzN7FO/TkZ/gzSrwuJ+FhSozG4XXp3UWkLGfVbAXaxYHCUKP4i2w9q6ZPYk33ukPwHTn3HJfKxNpWs8BZwI3AAMCi2jW2OmcW+xPWaHJObfWzF4H7g1MjjAfLyTcDfzLObfRz/qkcRQUgpRzboGZjQd+z3cz7nwNjHfOaZq35nER3krYEwOP+k4DvmjJgkSOlnNulpldDDyA92lrMd6qzD/3tTCRplczccIL+3ntP8CpLVZJ23EN3u+SH+INYl4H/A5vMLMEAdMYHhERERERqU9jFEREREREpAEFBRERERERaUBBQUREREREGlBQEBERERGRBhQURERERESkAQUFERERERFpQEFBREREREQaUFAQEREREZEGFBRERERERKQBBQUREREREWlAQUFERERERBpQUBARacPM7HQzc2b2qzrbTgtsu8fP2kRExF/mnPO7BhER8ZGZPQtcAQwFCoBvgQpgjHOuws/aRETEPwoKIiJtnJnFAUuADGAB8DtgnHNuvq+FiYiIrxQUREQEM7sAmAaUAY87537pc0kiIuIzBQUREQHAzJYDg4Ak59xmv+sRERF/aTCziIhgZtfihYQi4Df+ViMiIq2BWhRERNo4M+sJ5ACv4Y1TeBmY6Jz7zNfCRETEVwoKIiJtnJm9DZwIDHXOFZnZ50BfYIRzbo+/1YmIiF/U9UhEpA0zs8uBycCvnHNFgc23An2A3/tWmIiI+E4tCiIiIiIi0oBaFEREREREpAEFBRERERERaUBBQUREREREGlBQEBERERGRBhQURERERESkAQUFERERERFpQEFBREREREQaUFAQEREREZEGFBRERERERKQBBQUREREREWlAQUFERERERBpQUBARERERkQYUFEREREREpAEFBRERERERaeD/AxUOGkF+R7wSAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 720x480 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(dpi=120);\n",
"sns.lineplot(x=\"x\", y=\"y\", data=normal_df, ax=ax)\n",
"ax.set_title(\"Probability density function for normal distribution with \"\n",
" \"mean 0 and standard deviation 1\")\n",
"ax.set_xlabel(\"x\")\n",
"ax.set_ylabel(\"y\")\n",
"ax.set_xlim([-3, 3])\n",
"ax.set_ylim([0, 0.4]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Box-Muller-Gauss Method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* For normal distribution with mean µ and standard deviation σ\n",
"\n",
" * Compute \\\\(b \\sin(a) + \\mu\\\\) and \\\\(b \\cos(a) + \\mu\\\\) where \n",
" \n",
" a = uniform random number in [0, 2π) \n",
" rand = uniform random number in [0, 1) \n",
" b = \\\\(\\sigma \\sqrt{- \\ln(\\texttt{rand})}\\\\)\n",
"\n",
"Note that this means the *Box-Muller-Gauss method* outputs two random numbers, which is reasonable, since we need to provide two random numbers as inputs."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As an example, let's compute two samples from the standard normal distribution via the *Box-Muller-Gauss method*:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"random numbers = (-0.17825456071289894, 0.4784854000293813)\n"
]
}
],
"source": [
"a = 2 * np.pi * np.random.uniform()\n",
"b = 1 * np.sqrt(- np.log(np.random.uniform()))\n",
"rand1 = b * np.sin(a)\n",
"rand2 = b * np.cos(a)\n",
"print(\"random numbers =\", tuple((rand1, rand2)))"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.17825456071289894"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b * np.sin(a)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.4784854000293813"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b * np.cos(a)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### How to generate a distribution's random numbers from a probability distribution function"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You might be wondering how to derive the *Box-Muller-Gauss method* for sampling from a normal distribution.\n",
"Obtaining this formula, along with similar expressions for other functions, follows a standard recipe.\n",
"The requirement is that you are starting with a formula for the probability distribution function (PDF).\n",
"If you have that, then you can do the following:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<ol>\n",
" <li>\n",
" Assume that we have a continuous $\\texttt{PDF}(x)$ function defined for the range $[a, b]$.\n",
" Find the cumulative distribution function (CDF), which maps distribution values to their percentile rank, by integrating the PDF like so:\n",
"\n",
" $$\\texttt{CDF}(x) = \\int_{-a}^{x} \\texttt{PDF}(x') dx'$$\n",
"\n",
" where $a \\leq x \\leq b$.\n",
" </li>\n",
"</ol>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<ol>\n",
" <li value=\"2\">\n",
" The outputs of *any* CDF function, regardless of the initial PDF or the range of $x$, are uniformly distributed in the range $[0, 1]$:\n",
"\n",
" $$0 \\leq \\texttt{CDF}(x) \\leq 1$$\n",
" \n",
" This follows from 10% of the distribution existing below the 10th percentile, 20% of the distribution existing below the 20th percentile, and so on.\n",
" Thus, we can sample the $\\texttt{CDF}(x)$ outputs by generating a series of uniform random numbers `rand` between 0 and 1.\n",
" </li>\n",
"</ol>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<ol>\n",
" <li value=\"3\">\n",
" Once we have our uniform random number $\\texttt{CDF}(x)$, we then need to solve the following expression for $x$:\n",
"\n",
" $$\\texttt{rand} = \\texttt{CDF}(x)$$\n",
"\n",
" This involves inverting the $\\texttt{CDF}(x)$ expression, which may or may not have an analytic form.\n",
" If it does, then we will have an expression for $x$ that depends on constants and the random number `rand` only, like so:\n",
"\n",
" $$x_{\\text{random}} = \\texttt{CDF}^{-1}(\\texttt{rand})$$\n",
"\n",
" Note that $\\texttt{CDF}^{-1}$ indicates the inverse, not the reciprocal.\n",
" This allows us to draw directly from the distribution defined by the $\\texttt{PDF}(x)$.\n",
" </li>\n",
"</ol>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For example, imagine that we have $\\texttt{PDF}(x) = 3 x^{2}$ over the range $[0, 1]$.\n",
"We can sample from this distribution by first calculating the CDF:\n",
"\n",
"\\begin{equation}\n",
"\\texttt{CDF}(x)=\\int_{0}^{x}3x'^{2}dx'\n",
"\\end{equation}\n",
"\\begin{equation}\n",
"=x^{3}\n",
"\\end{equation}\n",
"\n",
"Now, we introduce `rand` and invert the $\\texttt{CDF}(x)$ expression:\n",
"\n",
"$$x_{\\text{random}} = \\left(\\texttt{rand}\\right)^{1/3}$$\n",
"\n",
"So, to sample from the distribution $\\texttt{PDF}(x) = 3 x^{2}$ defined in the range $[0, 1]$, we need to generate a uniform random number between 0 and 1, and then take the cube root."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"x_example = np.random.uniform(size=100000)**(1/3)\n",
"x_df = pd.DataFrame({\"x\": x_example})"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAApAAAAHQCAYAAAAWBVqhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAASdAAAEnQB3mYfeAAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XucHFWd9/HPLxNiICwIhCErt0SQ8JggS7wkKhE1CiqisLt4WbzgBR6WVRFWFlARVAS8g3eBddEFWW+LrguKGG+jMuyzBN0kyihsgkQIAxgQAiQk+T1/VE3odHpmuubWnZnP+/XqV89UnTp9qqtn5jun6pyKzESSJElq1qRWN0CSJEnbFgOkJEmSKjFASpIkqRIDpCRJkioxQEqSJKkSA6QkSZIqMUBKkiSpEgOkJEmSKjFASpIkqRIDpCRJkioxQEqSJKkSA6QkSZIqMUBKkiSpEgOkJEkalog4KyL+X0T8OSLuiYjvRsTcVrdLo8cAKamSiPg/EfGRiLggIp7X6vZIagvPBz4HPAd4IbAB+GFE7NrKRmn0RGa2ug3SmImIJwBvAZ5C8QtuI/AEoAe4JDM3tbB5TYuI5wOTM/OHY/y6uwLnA2/PzMfG8rVHQ0QcD/w8M29tYRv2Bw7NzMsj4lXAUzLzQ2P02lu83li/H2O9vwO0YyfgPOAVwJOA7YAzM/PDEfEV4CXArMxcO4S6nw78N/DWzPznEWx2W4uIHYEHgKMz87utbo9G3uRWN0AaKxGxPcV/yJ/MzM/VrXs28KmIeHv6X9VAXgFctK2Fx4g4LzPf2+p2NOGbQMdQNx7Cfg7r9arop21j9vqDuBJ4OXAtcAXFP5f/ERHPAF4HvGso4REgM2+KiG8D50XE1zLzoZFqdBUR8WHgGcABwHTgEeB24NvAZzLzvhF+yb+gOMv5pxGuV23CAKmJ5ATgQ416VzLzhojYDngV8LUxb9m2Y2/gqlY3Yrwqe8DHrBd8rF+v3V4fICIOpAiP12XmkXXrfgD8Gfj8MF/mAuBG4B0UPfitcCqwBLge6AWmAQuAc4ETI2JBZt4xgq93EfAroHsE61Qb8RpITSQzBjo1l5k/Aw4aw/ZsizooTvtL48ULy+dv1S6MiAOAFwFfz8xHhvMCmflfwC3A/42IVvW47pSZCzLzzZl5Zma+PTOfSRFonwScNVIvFBEfBZ4HHJuZ/r4Yp+yB1ESyLiL+MjPvGqDMFX1fRMTBwBspTmd1ADsA36q/7rDvujGKn6e+8juUqz+SmXeXvRyvL9c9Adge+Gxm/q5BPX8FHEoR1AJ4ELhwsD9i5Sn6twG7AlMoTiH9IDO/WVfuSODFwHqKfyJ3AC4v/8j1V/dB5b49B/iLiEjgj5n5yZp2HwS8AFgHfCUzl5bbTgVOA3aj6G16AnAP8NHMfHSk38e6dh9N8V4+t2wHwD9n5m/ryh0IvBVIHn/fL8zMNXXlmnqP+2nLDsAZ5TZJ8d7/DFhaU+b5NLi2dbBjNtB+DnR8Bni9Y4GFFO8zNPgMDnS9ZG29g7Rtq9cf4udl0GPXoI1/Q3EKvc8lEXFJ+fX/AY4v62p4RqLsnXwx8DeZ+e81ywP4F4rP8Icz88xy1b9R9Pa9CLhuoLaNhtr3rs7XgXdTXBe+hSHsIxHxcYrT/i9o5bXFGgOZ6cPHhHgAe1IExGlNlF1I8cu+o2bZJOBC4MV1ZY8HjgTOBibVLN8Z+AJwCPCeurqmUJwWe0JdPRcAL6urfyZwcd2y5wMvqvl+O4qehF3ryr0ZeFPN9y8E3lJXZjvgU8AzmnhfzqX4g1+//6cBr2lQfgpwObB/3fKnAv8KTBnp97Gfdp/Xz/LjgaPK9662LbsDH2/wPg36Hg/QhouBfeqWvQw4Bzi+0XGteswa7ecgx6f+c3Q88EHgqCY+g8fXH9f+6h2gbfWvX/Xz0tSx66eNzyk/z73AY+XX55bHo4Ni4MsG+vl9ARxMEVh/U/eZ/DhFmL2krvyLyuUfG6xtY/kA3lu2a6v3bAj7+CngbuCprd4vH6P/sAdSE0Zm/jEi3gNcGBErga9l5qp+ih+bme+o235TRLwX+AzFdUS1XkURArKm/AMRsQw4MTP/vq6u9RHxrxSjO79Ts+qhzLy2ruzKiFgaEYdm5s/7ae/xwKczc4sL1jPzSxHxqYi4IouBLy/IzLPryjwWEadTBIf/7qf+weybmZ9osPwUit7DLXoiMvM3ZU/FKcBHa1aN1PtYxauBN2bNqbbMvCciVkfErjXv6fE09x5vJSIOA36amX+o2/baiJg3SPtG4pj1d3waWZ91o2ab/AyOhKqfl2aP3VYy85cRcSPwT8BvM/PcvnURMY3iTMBvs5/BM5n56/Kz90aKXvHLI+LdFGH968BJdZv8v/J50KmvIuKdwBMHK1fjV5n57WYKRsS7gB0p/jF7BkXv8P9Q/HO8hSr7GBGfo+h5PBr4U0TMKFc9lC0aOKTRZYDUhJKZtwNvj4h9gL+NiCcDq4EvZ+Yfa4p+oZ/tN0TEnxusurY29NS4Dbi3n+b0APPrln2jn7L/StEz198f7+nZ/6n5vms7l1D+Ys/M1bUFMnNdGY6H6j/7WT4jM3/TaEVm/ioi3lK3eKTexyquzcbXad0K7Mvjo0ibfY8beTHwgX7WfQN49gDtG4lj1t/x6a89jQz2GRwJQ/m8NHPs+vNUissg6o/bnhS9kANd7gJF792rgXPLaWs+RHF6+vVZNyVY+Y/Qo8A+g9QJ8M6y/c36MsVo6ma8C9ij5vvvU/zTdk8/5Zvdx75/7hbXbf9+ip5djTMGSE1IZU/QRQARsSdwUkTcn5kfL9c3/CNWWt9gWX/B4RGK02ONPEDRC1Cr4YjUMixMa7SuXP6CiKivq8/OwIry6y8CH4uInwLfzcyHa16j0X41a0X9grJdvYNsd1dE7FjTQzFS72MV/V37+WBfvRXf40YmDfD+DjbIYCSO2UBtq7eh0cKBPoMjYYifl0GP3SD6en/rP3e7lc8DXkeZmasi4iLgTODTwC+Bvx7guPyJLcNbf/XOHKzMUGXmDICI2IPiNP6FwM0R8fLM3Ornr9l9zMwYrTarPRkgNeGVPY9nR8QREfGWLCf7LX/Bvh7YieJnpW/WgkY9CGMxd+TD/SzfDfj3zGzYa1qrDB8nRzG58dvLP9p3UQyqGNI8dwPYjUH+AFP0Ku4C9AWCdp2Ds+n3uB9D3q8xPmaD6e8zOBKG8nkZrr4AeXPd8r7BQlMZXG3P3VtqA34D29fU3VKZeTdwdUQsAX4HfAXo79aDVfZRE4QBUhNCObJzRmau7K9MZl4XER8py7+W4hqoT9afOoyI80azrQPo7z/8+yim4WhaZt4E3AQQEU8BvhARZ2TmncNr4hb+RDFaeSC7sG1MNFz5Pa6z/XAbMEbHbDC1n8GRDvut+LzMo9iPX9Ut7+sJ3Y0BlL8nPkZxGcwMims0/76fspMormsctDd4NK+BrJeZt0fEb4C/iojpmbnFpSJV9lETiwFSE8nrKG5XNpBHorjd4WGZWX8RfMtExBT6OYWbmWsjosr1UvXb/z4iTqX4w3D2YOUr1PtQ2Ys7kOkt6EWrbLjvMfBgRDwhM9eNUHtG5ZgNpPy5qP0M3k//IecJVJwvdKw/L+VUNAcDv8/MB+tW30XR6zZ7gO1fRnHt4XKKkfI/A94aERdn5i0NNplNEcDrw2ojo3kNZCN9/xxtccyGsI+aQJxIXBNCFnOg7R7F3WYaimKOv8kUc+w1nL+s/AP3tFFpZOHp/Sw/CvjeANv9T0S8qNGKsuej7+v397P9Gkagl6yB1RHx1H7aNZstT421u6be435cRzEQoZH+jnlf3WN9zBb0s/yVbPkZvJX+J95/2RBfeyw/LwdQzMlZf/qaciDXz4DpUdyrvL4th1LMIbkKOLwcgHI2xe+PrUYzl/re1x8P1rDMnJmZUeFx/ED1RcSBNaOia5dPiogPAZ3AL7Nm7swh7qMmEAOkJpJPU5z2a/QHYR/gk8AXy1+i+9cPmIiI/SguJF85im18URT336193b0p5vsb6JZgFwHHR8Qz67adSTHVS5/1EXFIXZlJFH8Ytpg+aIRcBJxRnnKtfc1ZFNOnXDwKr9nImnKw1HA0+x5vJTN/CcyNiL3qtn0mxZyjA6lyzEZiPw+NiC1CZPkZnFf7GczM5cCCiJheV/aVFIGkXjNtG8vPS38DaPr03ZnmiLq2HEwxqv0Bijlh7wLIYjL5/wZeGRGNjunhFD18Q51uajheAtwREYsj4pKIuCAivgT8nmIS8dUUt3oFhrWPmkA8ha0JIzNvjWIOtHdGxK4U1z5topiu40HgvTXX/5xLMW1F30jDAO6k+CP2zoiYVD9Nxwj5CPDCiDiOx0dkP0AxFUa/MnNjRLwJOK3c9tGyzffVbXs+RQh6Tfn9ZIp74n4jM380cruxuV3rIuIk4PSI2I3HR7D/CXjbSJ3SbcLlwPkR8RDwi6y5q0azKrzH/TmH4n3Yqfx+EsVdaC6imIuvP1WO2eUMcz8prnd7XkS8muIzmBSfwXMblP0AcF4UdyZaRzER+C8oRo7X/30ZtG1j/HlpJkDeDbwB+CxA+c/ndRTvyRGZeVvdNmdRzBH7UWp6cst/Ro8G/jNH9n7TzfohcAnwXIrT9k8E1lIMnvlX4FNZzpk51H3UxBONp1yTNNZigFvDSRp7EXEWRYCfl5lbnequUM/bKe7S8rzM7Bqp9kmt5ClsSZIa+yTwB/qfBH5Q5bXVZwHfMjxqPDFASpLUQDn47vXAf8fQJ1GfSXH6+F0j1S6pHXgKW2oTnsKWJG0rDJCSJEmqxFPYkiRJqsQAKUmSpEoMkJIkSarEAClJkqRKvBPNCCnvNHAYcAeP3z1BkiSpHU0B9gZ+mpkPVN3YADlyDqM19ziVJEkaqlcC/1F1IwPkyLkD4Nvf/jb7779/q9siSZLUr1tvvZWjjz4ayvxSlQFy5KwH2H///ZkzZ06r2yJJktSMIV125yAaSZIkVWKAlCRJUiUGSEmSJFVigJQkSVIlBkhJkiRVYoCUJElSJQZISZIkVWKAlCRJUiUGSEmSJFVigJQkSVIlBkhJkiRVYoCUJElSJQZISZIkVWKAlCRJUiUGSEmSJFUyudUNkCRJmihmnnlNU+VWXnjkKLdkeOyBlCRJUiUGSEmSJFVigJQkSVIlBkhJkiRVYoCUJElSJQZISZIkVWKAlCRJUiUGSEmSJFXiROKSJEnD1OwE4eOFPZCSJEmqxAApSZKkSgyQkiRJqsQAKUmSpEoMkJIkSarEAClJkqRKDJCSJEmqxAApSZKkSgyQkiRJqsQAKUmSpEoMkJIkSarEAClJkqRKDJCSJEmqZHKrGyBJktSOZp55Taub0LbsgZQkSVIlBkhJkiRVYoCUJElSJQZISZIkVWKAlCRJUiUGSEmSJFVigJQkSVIlBkhJkiRVYoCUJElSJQZISZIkVWKAlCRJUiUGSEmSJFXSdgEyIq6NiGzwOKlcv3dEfDMi7o+IRyLihog4rEE9J0TEbyNiXUTcGREXRcTUujJN1SVJkqTHTW51AxrYE/gmcHHd8lsjYjLwPWBn4B3AQ8CpwHUR8bTM/B1ARLwGuAS4FPhH4GDg/cBOwJvLMk3VJUmSpC21Y4DcC/hqZv68fkVEvAyYAyzKzB+Vy64HVgEnAaeVRU8DujLzxPL7ayNiO+DsiDg9M+8DDm+yLkmSJNVoq1PY5SnmXYE7+imyAFgP/LhvQWY+CHQDzyrrmAIcAlxXt+31FIH5kGbrkiRJ0tbaKkBSnL4GeFtE3BcRj0VEd0QsKJfPAO7NzKzbrhfoLL+eThEUexuUoaZcM3U1FBGdETGn9gHsN+jeSZIkjQPtdgp7KsUp5KXAe4DdgA9QXJf4FGAK8FiD7TaW66h5ri+3sW59M3X152TgnEYruru7Wb16NXvttRezZs2iq6tr87pFixaxdOlSenuLLDt79mymTZvGkiVLAJg2bRoLFiygu7ubtWvXAjBv3jzWrl1LT08PAJ2dnRx00EEsXrx4c70LFy5kxYoVrFq1CoBZs2bR2dnJjTfeWOzolCksXLiQJUuWsGbNGgDmzp0LwLJlywDYZZddmDdvHl1dXaxfvx6A+fPn09vby4oVKwDcJ/fJfXKf3Cf3aVzs0629D3Hvo8GlPR2cMHsj06cWfUlX3tbB7lOTw/fcBMAx+wZX397BWQdv2Fzvxcs7OHSP5OnTizJdqydxywPBCbOLmPHwBrh4+WSO228j++xY1Hv17ZPK+opt/vBQcOVtHZwyZwM7lEns0p4ODtw5WTijKNPT0zOqx6m7u5vhiK074NpLRMwCbgEuBPYFnp+ZM+vKXA4clpmzImImsAJ4U2ZeXlOmb/nxmfnlcpsB6xqgTZ3A7nWL9wO+s2zZMubMmVNtJyVJ0piZeeY1rW7CoFZeeOSo1r98+fK+sD43M5dX3b7deiC3kpkrIuK3wFzgPmD7BsWmAOvKr9eXz/Xl+noVa8sNVld/beql7hR5RAy0iSRJ0rjRVtdARsTCiGgUuTuABO4BditHVNfaA7ir/Po+YBPFNY71Zagp10xdkiRJqtNWARI4FPi3iHhi34KI2BuYDSwDbqAIky+tWb8T8GyK0dNk5jrgZuDldXUfQdHreHP5/aB1SZIkaWvtFiC/BDwMXBMRR0fEscB/AmuBzwLfB5YDl0XEWyLiaOC7FL2TX6yp5+PAvIi4IiKOjIjTgdOByzLzz2WZZuuSJElSjbYKkJl5N3AYcC9wGXA5xanmRZl5T2ZuAF4C/BT4BHAVxXWch2fmypp6rgJOBJ4BXE1xh5lPl899ZZqqS5IkSVtqu0E0mXkL8MoB1q8Cjm2inkspbmU4UJmm6pIkSdLj2qoHUpIkSe3PAClJkqRKDJCSJEmqpO2ugZQkSarV7J1jRvvuLXqcAVKSJI0L28ItCscLT2FLkiSpEgOkJEmSKjFASpIkqRIDpCRJkioxQEqSJKkSA6QkSZIqMUBKkiSpEgOkJEmSKjFASpIkqRIDpCRJkioxQEqSJKkSA6QkSZIqMUBKkiSpEgOkJEmSKjFASpIkqRIDpCRJkioxQEqSJKmSya1ugCRJmphmnnlNq5ugIbIHUpIkSZUYICVJklSJAVKSJEmVGCAlSZJUiQFSkiRJlRggJUmSVIkBUpIkSZU4D6QkSRpRzu84/tkDKUmSpEoMkJIkSarEAClJkqRKDJCSJEmqxAApSZKkSgyQkiRJqsQAKUmSpEoMkJIkSarEAClJkqRKDJCSJEmqxAApSZKkSgyQkiRJqmRyqxsgSZK2DTPPvKbVTVCbsAdSkiRJlRggJUmSVIkBUpIkSZW0dYCMiK9HREbE82uWzYmI6yLiofJxfUQ8rW67iIj3RMSKiFgfESsj4n0RMamu3KB1SZIkaUttGyAj4kjg2LplOwOLgScBbwbeWn69OCJ2rSl6BvBB4N+AVwJfA84F3juEuiRJklSjLUdhR8SOwOeAXwDPrVl1HNAJPDczbyvLdgP/C7wO+FREBPBO4IrMPKvc7nsRsTvw9og4LzM3NVPXKO+mJEnSNqldeyA/BDxWPtdaANzWF/gAMnMl8HvgWeWiWcAewHV1214PTAeeXKEuSZIk1Wm7ABkRzwT+Afh74JG61TOA3gab9VL0JvaVoUG5vu9ryw1WlyRJkuq01SnsiJgMXApclZnX1w6eKU2h6Jmst7FcR81zfbmNdeubqau/dnYCu9ct3m+gbSRJksaLtgqQwLuAvYEXD2HbHMFyg5U5GTin0Yru7m5Wr17NXnvtxaxZs+jq6tq8btGiRSxdupTe3qLjc/bs2UybNo0lS5YAMG3aNBYsWEB3dzdr164FYN68eaxdu5aenh4AOjs7Oeigg1i8ePHmehcuXMiKFStYtWoVALNmzaKzs5Mbb7wRgClTprBw4UKWLFnCmjVrAJg7dy4Ay5YtA2CXXXZh3rx5dHV1sX79egDmz59Pb28vK1asAHCf3Cf3yX1ynyb4Pp0weyPTpxZ/Iq+8rYPdpyaH77kJgFvuD66+vYOzDt6wud6Ll3dw6B7J06cXZbpWT+KWB4ITZhd9Og9vgIuXT+a4/Tayz45FvVffXpwcPWbfYps/PBRceVsHp8zZwA5larm0p4MDd04WzijK3HTvJH5+d3DKnI2bX/uCX0/mmH03cuATi3p/8MdJ3PNocNx+RZl7Hw0u7elo233q6ekZ1c9ed3c3wxGZzeau0RUR+wFLgVOAfykXHwb8EHgR8BPge8DOmTm/bttfAn/OzJdExHMoBt+8LDO/V1PmcIrrIp+dmd0R8YPB6hqgrf31QH5n2bJlzJkzp9rOS5K0DfBWhmNn5YVHjmr9y5cv7/sHZG5mLq+6fTv1QL4L2B64pHzU+iFwO0UwnN1g2z2AnvLre8rnGQ3KANxVU26wuhrKzF7qrp8sBn9LkrRtMRRqKNppEM2HgWfWPU4q150EHAXcAOxTO9l3RDyFYmR1X1/sbRTh8Ki6+o8A7gT+UH7fTF2SJEmq0zY9kOUUOitrl5XzQQL0ZObSiLiDYjLwb0XE+4ANwPuA1cBVZT2bIuKTwPkRcRHFaeuFwN8BZ+Tj5+yvGKwuSZIkba1tAmQzMvP+iFgEfAK4jGKwyy+AV2fmn2uKXggEcCLFgJc7gfcAHxtCXZIkSarR1gEyM39CEQRrly2nOB090HYJnF8+Bio3aF2SJEnaUjtdAylJkqRtgAFSkiRJlRggJUmSVIkBUpIkSZUYICVJklSJAVKSJEmVGCAlSZJUiQFSkiRJlRggJUmSVIkBUpIkSZUYICVJklSJAVKSJEmVGCAlSZJUiQFSkiRJlRggJUmSVMnkVjdAkiSNvJlnXtPqJmgcswdSkiRJlRggJUmSVIkBUpIkSZUYICVJklSJAVKSJEmVOApbkqRtiKOr1Q7sgZQkSVIlBkhJkiRVYoCUJElSJQZISZIkVWKAlCRJUiUGSEmSJFVigJQkSVIlBkhJkiRVYoCUJElSJd6JRpKkNuAdZrQtsQdSkiRJlRggJUmSVIkBUpIkSZUYICVJklTJsAJkRLwhIqb3s27/iHjdcOqXJElS+xluD+S/AE/pZ9084PPDrF+SJEltZkjT+ETEj/q+BD4TEQ/UFekADgFWDKNtkiRJakND7YGM8pE1X9c+NgDfAl49Am2UJElSGxlSD2RmvgAgIjYBb8vMX45oqyRJktS2hnsN5JuA341EQyRJkrRtGNatDDPzyxGxW0S8BNid4vR1fZmvDOc1JEnalnmLQo1HwwqQEfEK4EpgBxqER4prJA2QkiRJ48iwAiRwIfAr4B+Bu4ffHEmSJLW74QbImcBpmflfI9AWSZIkbQOGO4jmFmC3kWiIJEmStg3DDZBnAB+IiINGojGSJElqf8MNkF+nOI39q4jY2OCxoUplETE5Ik6LiN9GxMMRcVtEXBAR02rK7B0R34yI+yPikYi4ISIOa1DXCWU96yLizoi4KCKm1pVpqi5JkiQ9brjXQF5MMdJ6pHwB+DvgXOBm4EDgg8CewBsiYjLwPWBn4B3AQ8CpwHUR8bTM/B1ARLwGuAS4lGKAz8HA+4GdgDeXZZqqS5IkSVsa7jyQ545QO4iIvwCOBy7IzI+Ui6+PiN2Ad0fEW4EXAXOARZn5o3K764FVwEnAaeV2pwFdmXli+f21EbEdcHZEnJ6Z9wGHN1mXJEmSagx3HsjnDVYmM3/WZHVrgb2AB+qWP0TRzu2BBcB64Mc19T8YEd3As8o2TQEOoejFrHU9RS/kIcAPm6lLkiRJWxvuKeyfUJzCrp9EvPa0dkczFWXmJmA1QER0AFOB5wDvAr6UmQ9ExAzg3sysP23eC8wvv55OsV+9DcoAdJbPzdQlSZKkOsMNkC/oZ/kewGeA9w6x3l9TnF4G+D7w9+XXU4DHGpTfWK6j5rm+3Ma69c3U1VBEdFLcurHWfgNtI0kaX7xFoSay4V4D+dP+1kXEdOB1FINZqno1sAtFD+T7gI9TDHQZsDlN1t1MucHKnAyc02hFd3c3q1evZq+99mLWrFl0dXVtXrdo0SKWLl1Kb2/RGTp79mymTZvGkiVLAJg2bRoLFiygu7ubtWvXAjBv3jzWrl1LT08PAJ2dnRx00EEsXrx4c70LFy5kxYoVrFq1CoBZs2bR2dnJjTfeCMCUKVNYuHAhS5YsYc2aNQDMnTsXgGXLlgGwyy67MG/ePLq6uli/fj0A8+fPp7e3lxUrVgC4T+6T++Q+uU81+3TWwcVEIz/44yTueTQ4br+in+LeR4NLezo4YfZGpk8t/pxceVsHu09NDt9zEwC33B9cfXvH5joALl7ewaF7JE+fXpTpWj2JWx4ITphd1PvwBrh4+WSO228j++xY1Hv17cVkKsfsW2zzh4eCK2/r4JQ5G9ih/At/aU8HB+6cLJxRlLnp3kn8/O7glDkbN7/2Bb+ezDH7buTAJ6b71Cb71NPTM6o/T93d3QxHbH0Gd2RExCLgW5n5xGHWcybwIWBviusaX5mZe9SV+SowLzMPjIgnAX8ETs7Mz9eUOQDoAV6bmf8WEZcMVtcAbeqvB/I7y5YtY86cOQ22kiSNJ/ZAajStvPDIUa1/+fLlff9Uzc3M5VW3H+4p7IbKaxhfA6ypsM0M4Gjgmsy8o2bV/1DMV7kvcA+wW0Rsl5m1p5/3AO4qv74P2ERxjSN1Zagp10xdDWVmL3XXWEbUXwYqSZI0Pg1rIvGIWBER/1v/AO6nmG/x4grVPQH4PHBc3fKDy+fbgRsoBuW8tKYNOwHPBroBMnMdxRySL6+r5wiKUdc3l98PWpckSZK2NtweyJ+y9fWCSTGa+geZ+ZNmK8rM2yPi6xRzNQLcRBEe3w38a2beGRG9wHLgsog4i6K38dTyNb9YU93Hga9GxBXAVcBTgdOByzLzz2WZ7zdZlyRJkmoMdxDN8SPUjj6vp7hzzJsoBs/8geL6x4+Xr7chIl4CfBL4BMVo6SXA4Zm5sqZdV0XEjmVdrwLuBT5NEUapUpckSZK2NCLXQEbEk4FFFHMw3gP8KDP/t2o9mbkeuKB89FdmFXBsE3VdSnErw4HKNFWXJEmSHjfsABkRFwFvY8uiA2BRAAAW60lEQVTrKTMiPpuZg029I0mSpG3McAfRnE4RHj8IzKK43eB+FNPtnBQRpw63gZIkSWovw+2BPBH4cGa+v2bZCuCD5T2pT6K4xlCSJEnjxLB6IIG9gK5+1v0C2GeY9UuSJKnNDLcH8naK2w1+v8G6Z1PcEUaSpG2Gd5iRBjfcAPl54KMRsR64giIw7kExGfiZwNnDrF+SJEltZrjzQF4cETOB95ePWl/IzI8Mp35JkiS1n2FP45OZp0bEZ4EXALsDfwJ+npnLhlu3JEmS2s+wAmQ50voiijvHHJyZvyuXr4qIbwPvyMxNw2+mJEmS2sVwR2GfBbwBOAe4o2b52cDrqLl1oCRJksaH4QbI1wNnZ+ZHMvORvoWZ+S8U97J+yzDrlyRJUpsZboD8S+A3/az7HTBjmPVLkiSpzQw3QPYAL+9n3ZHArcOsX5IkSW1muKOwPwpcGRG7At8F7gY6gaOA1wJvHmb9kiRJajPDnQfyqoiYCnyIIjAmEMC9wD9k5peH30RJkobHu8tII2sk5oH8l4i4HDgQ2I1iHshbnL5HkiRpfBp2gATIzAR+OxJ1SZIkqb0NdxCNJEmSJpgR6YGUJKkVvLZRag17ICVJklSJAVKSJEmVGCAlSZJUiQFSkiRJlRggJUmSVIkBUpIkSZUYICVJklSJAVKSJEmVGCAlSZJUiQFSkiRJlRggJUmSVIkBUpIkSZUYICVJklTJ5FY3QJI0ccw885pWN0HSCLAHUpIkSZUYICVJklSJAVKSJEmVGCAlSZJUiQFSkiRJlRggJUmSVIkBUpIkSZUYICVJklSJAVKSJEmVGCAlSZJUiQFSkiRJlRggJUmSVIkBUpIkSZUYICVJklTJ5FY3QJK07Zt55jWtboKkMdRWPZARMSkiTo+I30XEuohYGREXR8QTa8rsHRHfjIj7I+KRiLghIg5rUNcJEfHbsp47I+KiiJhaV6apuiRJkvS4tgqQwEXAecDXgFcAFwNvLr8nIiYD3wPmA+8AjgM2ANdFxAF9lUTEa4BLgC7gGODTwMnA52rKNFWXJEmSttQ2p7AjYnfgH4ALMvPscvF1EdEBfDQi9gXmlI9FmfmjcrvrgVXAScBp5XanAV2ZeWL5/bURsR1wdkScnpn3AYc3WZckSZJqtFMP5M7AVcC36pb/pnzeFVgArAd+3LcyMx8EuoFnAUTEFOAQ4Lq6eq6nCMyHlN8PWpckSZK21jY9kJl5K/C6BqteANwNLAf+Hrg3M7OuTC/FqWiA6RT71dugDEBn+TyjiboaiohOYPe6xfsNtI0kSdJ40TYBspGIeDrF9YmnZub6snfxsQZFNwJTyq/7nuvLbaxb30xd/TkZOKfRiu7ublavXs1ee+3FrFmz6Orq2rxu0aJFLF26lN7eIsvOnj2badOmsWTJEgCmTZvGggUL6O7uZu3atQDMmzePtWvX0tPTA0BnZycHHXQQixcv3lzvwoULWbFiBatWrQJg1qxZdHZ2cuONNxY7OmUKCxcuZMmSJaxZswaAuXPnArBs2TIAdtllF+bNm0dXVxfr168HYP78+fT29rJixQoA98l9cp8m4D59+0fdAPzhoeDK2zo4Zc4Gdij/clza08GBOycLZ2zirIPhpnsn8fO7g1PmbNz82hf8ejLH7LuRA59Y/K/+gz9O4p5Hg+P2K8rc+2hwaU8HJ8zeyPSpRZkrb+tg96nJ4XtuAuCW+4Orb+/grIM3bK734uUdHLpH8vTpRZmu1ZO45YHghNlFvQ9vgIuXT+a4/Tayz45FvVffXpx0O2bfTU3tE7hP7lPr9qmnp2dUf0d0d3czHLF1B1x7iIi9gRuAm4CjMzMj4nLg+Zk5s67s5cBhmTkrImYCK4A3ZeblNWX6lh+fmV9upq4B2tZfD+R3li1bxpw5cyrtqyS1K6fnkVpj5YVHjmr9y5cv7/tHcW5mLq+6fVv2QEbErhTXMN4J/F3Naeb1wPYNNpkCrKspQ4Nyfb2KteUGq6uhzOyl7hR5RAy0iSRJ0rjRToNoAIiIHYBrgARemplra1bfA+xWjqiutQdwV/n1fcAmimsc68tQU66ZuiRJklSnrQJkOTfjNyjC34vL6XZq3QB0AC+t2WYn4NkUo6fJzHXAzcDL67Y9gqLX8eZm65IkSdLW2u0U9mXAi4ETgCdHxJNr1j0AfJ9iNPZlEXEWRW/jqRS9lV+sKftx4KsRcQXF1EBPBU4HLsvMP5dlmq1LkiRJNdotQL6xfL68wbqfZubzI+IlwCeBT1Bcr7gEODwzV/YVzMyrImJH4B+BVwH3UtyN5t01ZTY0U5ckSZK21FYBMjMHHYmSmauAY5sodylw6UjUJUmSpMe11TWQkiRJan8GSEmSJFVigJQkSVIlbXUNpCRp9Hl3GUnDZQ+kJEmSKjFASpIkqRIDpCRJkioxQEqSJKkSB9FI0jjh4BhJY8UeSEmSJFVigJQkSVIlBkhJkiRVYoCUJElSJQZISZIkVWKAlCRJUiVO4yNJbc7peSS1G3sgJUmSVIkBUpIkSZUYICVJklSJAVKSJEmVGCAlSZJUiQFSkiRJlRggJUmSVIkBUpIkSZUYICVJklSJd6KRpBbxDjOStlX2QEqSJKkSA6QkSZIqMUBKkiSpEgOkJEmSKnEQjSSNMAfHSBrv7IGUJElSJQZISZIkVWKAlCRJUiUGSEmSJFVigJQkSVIlBkhJkiRVYoCUJElSJc4DKUlNcn5HSSrYAylJkqRKDJCSJEmqxAApSZKkSgyQkiRJqsQAKUmSpEoMkJIkSarEAClJkqRKnAdS0oTn/I6SVE1b9kBGxAERcX5E3NFg3d4R8c2IuD8iHomIGyLisAblToiI30bEuoi4MyIuioipQ6lLkiRJj2urABkRL46IXwI9wOnAXnXrJwPfA+YD7wCOAzYA10XEATXlXgNcAnQBxwCfBk4GPle1LkmSJG2prQIk8Dzg98CLgPMbrD8cmAO8MTO/kpn/DrwMWAecVFPuNKArM0/MzGsz8wLgPOD1EbFbxbokSZJUo60CZGaenZlvzMzFQDYosgBYD/y4ZpsHgW7gWQARMQU4BLiubtvrKa75PKTZuiRJkrS1bW0QzQzg3sysD5e9FKeiAaZT7FdvgzIAnRXqkrSNcmCMJI2ebS1ATgEea7B8Y7mOmuf6chvr1jdTV0MR0QnsXrd4v4G2kSRJGi+2tQA5kEanvIdabrAyJwPnNFrR3d3N6tWr2WuvvZg1axZdXV2b1y1atIilS5fS21t0hs6ePZtp06axZMkSAKZNm8aCBQvo7u5m7dq1AMybN4+1a9fS09MDQGdnJwcddBCLFy/eXO/ChQtZsWIFq1atAmDWrFl0dnZy4403AjBlyhQWLlzIkiVLWLNmDQBz584FYNmyZQDssssuzJs3j66uLtavXw/A/Pnz6e3tZcWKFQDuk/u0Te3TCbM3cmlPByfM3sj0qcWP9JW3dbD71OTwPTcBcMv9wdW3d3DWwRs213vx8g4O3SN5+vSiTNfqSdzyQHDC7OJ/0Ic3wMXLJ3PcfhvZZ8ei3qtvL64GOmbfYps/PBRceVsHp8zZwA7lb9lLezo4cOdk4YyizE33TuLndwenzNm4+bUv+PVkjtl3Iwc+saj3B3+cxD2PBsftV5S599Fwn9wn92mC7FNPT8+o/i7v7u5mOGLrM7jtISLOBc7JzKhZdgnwyszco67sV4F5mXlgRDwJ+CNwcmZ+vqbMARSju1+bmf/WTF0DtK2/HsjvLFu2jDlz5gxhjyWNJE9hS9qWrbzwyFGtf/ny5X3/0M/NzOVVt9/WeiDvAXaLiO0ys/b08x7AXeXX9wGbKK5xpK4MNeWaqauhzOyl7hrLiOintCRJ0viyrQXIG4AO4KXAfwBExE7As4GLATJzXUTcDLycLU8zH0Ex6vrmZuuS1H7sWZSk1tvWAuT3geXAZRFxFkVv46kU1yx+sabcx4GvRsQVwFXAUykmJr8sM/9csS5JkiTV2KYCZGZuiIiXAJ8EPkExWnoJcHhmrqwpd1VE7Aj8I/Aq4F6Ku9G8u2pdkiRJ2lLbBsjMPBc4t8HyVcCxTWx/KXDpIGWaqkuSJEmPa9sAKWli8dpGSdp2tNWtDCVJktT+DJCSJEmqxAApSZKkSgyQkiRJqsQAKUmSpEoMkJIkSarEAClJkqRKnAdS0qhyfkdJGn/sgZQkSVIlBkhJkiRV4ilsSUPiqWlJmrjsgZQkSVIlBkhJkiRVYoCUJElSJQZISZIkVeIgGklbcHCMJGkw9kBKkiSpEgOkJEmSKjFASpIkqRIDpCRJkipxEI00QTg4RpI0UuyBlCRJUiUGSEmSJFVigJQkSVIlBkhJkiRVYoCUJElSJY7ClrZhjqyWJLWCPZCSJEmqxAApSZKkSjyFLbUhT01LktqZPZCSJEmqxB5IaQzZsyhJGg/sgZQkSVIlBkhJkiRVYoCUJElSJV4DKY0Ar22UJE0k9kBKkiSpEgOkJEmSKjFASpIkqRKvgZQG4LWNkiRtzR5ISZIkVWIPpCYkexYlSRo6eyAlSZJUiQFSkiRJlXgKW+OKp6YlSRp99kBKkiSpEnsg1fbsVZQkqb1M+AAZEUcD5wIHAg8B/wmclpl/amW7JgKDoSRJ26YJfQo7Ip4LfAu4DTgWOAc4BvhmK9slSZLUziZ6D+QpwB3AsZm5CSAiHgS+HBFPy8z/aWnrtlH2LEqSNL5N9AC5ALiuLzyWri+fnwUYIGsYDCVJEhggZwC9dcv6vu8c47a0jMFQkiRVMdED5HbAY7ULMnNjRABM6W+jiOgEdq9bfCDArbfeOsJN3NqLP/HTUX8NSZLUOsuXLx/V+mvySr95ZyATPUAOJAdYdzLFgJutHH300aPTGkmSNGHM/dKYvdTewM1VN5roAfIxYPvaBRHRl8TXDbDd54Bv1C3bETgAWAasH0ab9gO+A7ySYnS42ofHpn15bNqXx6Z9eWza11gcmykU4XFIpzUneoC8h+I6yFp7lM939bdRZvay9bWTADcOt0Hl6XOA2zJzdPuvVYnHpn15bNqXx6Z9eWza1xgem8o9j30m9DyQwA3AERGxXc2yI8rn7ha0R5Ikqe1N9AB5EcVgmP+IiFdExEnAx4BrM7OntU2TJElqTxM6QGbmz4G/Bv6S4prG8yjuQvN3rWyXJElSO5vo10CSmd+huFC1XdwDvL98Vnvx2LQvj0378ti0L49N+2r7YxOZA81WI0mSJG1pQp/CliRJUnUGSEmSJFVigJQkSVIlBkhJkiRVYoAcQxFxdET8KiIejYh7I+LyiNi1ie1OiIjfRsS6iLgzIi6KiKlj0eaJYijHJiKmRsQFEbGy3O73EXGux2ZkDfXnpmb7SRFxQ0RkRMwcvZZOPMM5NhHxrIj4RkTcFRH/PdptnWiG+Dttu4g4KyJ+FxEPl89nRMSEn7FlNETEARFxfkTc0UTZiIj3RMSKiFhf/t15X0S0Lsdlpo8xeADPBTYC3wKOAv4BeAD40SDbvQZI4BLgZcBZFPfa/lKr92m8PIZxbP6jLHc6cDhwDsX91S9p9T6Nl8dQj01dHW8rf4YSmNnqfRovj+EcG+B15bY/A04EntXq/RlPj2H8TrsA2AC8r/yd9kFgE3Beq/dpPD2AFwO/LH8nPVZEsUG3ObM8FhcALwU+XH7/vpbtR6vfyInyAL4OrAQm1Sx7Q/kBetoA2/0X8LO6Ze8rP3S7tXq/xsNjKMcGmFeuf2vd8s8A64DtWr1f4+Ex1J+bmrJ7ln84f26AbI9jA+wDrAU+RTmVnI+2OTZ/BL5St+yrwB2t3qfx9CiD+ZeBRRRzPeYg5QNY3eDYfIlinshJo9HOwR6ewh47C4DrM3NTzbLry+dnNdogIqYAhwDX1a26nmIS+ENGupETVOVjAzwBuBL4ft3y3wBTgB1HtIUT11COTa3PAkuBfx7phmnIx+ZE4EHgn7L8K6gRN9RjswPFP1y17sPfZyMqM8/OzDdm5mKKUD+YWcAeNM4C04Enj3ATm2KAHDszgN66ZX3fd/azzXSKoFh1O1VT+dhk5g2Z+brMXFW36gXArzJzzQi3caIays8NABHxNxSXfZxIc7+kVc1Qj82LgOXA5RHxUHmN3nci4kmj0cgJaqjH5ovAGyLiiIjYKSKOoOi5/MIotFHNm1E+t1UWMECOne0oTjtvlpkbyy+n9LNN3/LH6pYPtp2qGcqx2UpEHAX8LfDekWvahDekYxMRO1OcIv1wZv5m9Jo3oQ315+YpwAuBDuCVwP+luGbv3yMiRqGdE9FQj827gZ9QnFl5oHy+DnjPyDdRFbRlFnBkVXsYau+IvSqjr6n3OCIOpjil/enMvGZ0m6TSQMfmwxTX2X1ojNqiLQ10bHYGvp+Zx/YtiIhHgK9R9OD/aJTbNtENdGzeRTHA458orr9/DsXgwHcAF41+0zRELckCBsix8xiwfe2C8hpHKAZdNLK+fN6+bvlg26maoRyb2rKzKP5T/wlw6kg3boKrfGwi4rnACRQjFTeUU5D0nW3piIiOmt4YDd1wfm5uqvu+7/q8gzBAjoSh/NzsCnwAOD8zP1ou/mk5LdkFEfFlL81pmbbMAp7CHjv38Ph1DH32KJ/v6meb+yiG6VfdTtUM5dgAEBGdwA+AHuBVBpMRN5Rj816K323XUfwhfYzHB9HcCiwe4TZOVEP9ubmTYiR2I/Wn6DQ0Qzk2B1AMDqyfk/O/gKnlerXGPeVzW2UBA+TYuQE4IiK2q1l2RPnc3WiDzFwH3Ay8vG7VERT/kdw80o2coCofG4CI2BG4luJaoaMy89HRa+KENZRj807gmXWP95frXkFxzZ2Gb0g/NxTB/qjyOtX67X41gu2byIZybPpCyMF1y/+qfF49Qm1TdbdRhMij6pYfQfEP2R/GvEXgPJBj9QAOpbjg9XsUf8ROAu4HrqkpMwl4KjVzCAKvpbi+4QrgSIpJq9cBn231Po2Xx1CODcVF6tdThMejyzpqH09p9X6Nh8dQf24a1HM8zgPZFscG2BdYA9wIHEMR6P8EXNfqfRovj2Ecm28BDwFnUEwkfgbFtcTfbPU+jdcHcC5180CWf1+eypbzeJ5V/g67iOLynPMpzlCe3rK2t/rNm0gPihGHvyoD4L3AZcDONesXlD/kf1233QnALRS9jncCHwOmtHp/xtOj6rEBZvL43U0aPS5v9T6Nl8dQf27q6jBAttGxAZ4G/BB4pAyP/ww8sdX7M54eQzk2FKeq30dxSc7D5fP7gKmt3p/x+ugnQP5teWyeVbMsKEbJryyzwMoyVLZsMv4oGyZJkiQ1xWsgJUmSVIkBUpIkSZUYICVJklSJAVKSJEmVGCAlSZJUiQFSkiRJlRggJUmSVIkBUpIkSZUYICVJklSJAVKSJEmVGCAlSZJUiQFSktpMRLwwIjIizqhZ9oJy2Xtb2TZJAojMbHUbJEl1IuJS4LXAbKAX+DXwGPCMzHyslW2TJAOkJLWhiNgZ+A3QBdwMfAiYn5k3tbRhkoQBUpLaVkS8AvgOsA74VGb+U4ubJEmAAVKS2lpE9AD7A3tm5upWt0eSwEE0ktS2IuJ4ivC4BvhAa1sjSY+zB1KS2lBE/CWwHPgaxXWQVwKLMvNHLW2YJGGAlKS2FBFXA88FZmfmmoj4MbAPcFBmPtza1kma6DyFLUltJiJeAxwNnJGZa8rF/wDsDZzfsoZJUskeSEmSJFViD6QkSZIqMUBKkiSpEgOkJEmSKjFASpIkqRIDpCRJkioxQEqSJKkSA6QkSZIqMUBKkiSpEgOkJEmSKjFASpIkqRIDpCRJkioxQEqSJKkSA6QkSZIqMUBKkiSpkv8PICZtCWyDDdYAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 720x480 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(dpi=120)\n",
"sns.distplot(a=x_df[\"x\"], bins=50, kde=False, ax=ax, hist_kws={\"alpha\": 1})\n",
"ax.set_xlabel(\"x\")\n",
"ax.set_ylabel(\"count\")\n",
"ax.set_title(r\"Samples from the distribution $f(x)=3x^2$\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Random Walk"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Random walk** refers to the apparently random movement of an entity.\n",
"\n",
"* An entity lives on a cell in a rectangular grid.\n",
"\n",
"* At each time step, the entity can step onto a neighboring cell\n",
"\n",
"* The step direction is chosen using a random number generator\n",
"\n",
"* There can be constraints on how probable each step direction is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Cellular automata"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* **Ceullular automata** are one style of computer simulation that involves movement on a grid.\n",
"\n",
"* These are computational models that are discrete in space, state, and time.\n",
"\n",
"* Structure of a cellular automata simulation\n",
"\n",
" * Space is built as a one-, two-, or three-dimensional **grid**.\n",
" \n",
" * A **site**, or **cell**, of the grid has a state, and the number of states is finite\n",
" \n",
" * There are **rules**, or **transition rules**, specifying local relationships and indicating how cells are to change state and regulate the behavior of the system"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"initial_state = np.random.randint(low=0, high=2, size=100)\n",
"nine_states = [np.random.choice(a=initial_state, size=(10, 10), replace=False)\n",
" for _ in np.arange(9)]"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"def draw_heatmap(*args, **kwargs):\n",
" data = kwargs.pop('data')\n",
" sns.heatmap(data[\"state\"].values[0], **kwargs)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAArQAAALICAYAAACO3mNgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAASdAAAEnQB3mYfeAAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3X2wZXtZJ/bvj9sXuJfBiKAERoVSSTOAZWkRw4sYxtJBhYoOOqNpRKU1RDRSVCVllzDEJXFgGp0SGdQI0mKqbCUq4qhlUY5FqARQfEEsUJsJEkcQFHlJUO+FC/fJH3v37X0Pp2+/rHXuWc85n09VV52z9upn/9bLs/d3//ba+4yqCgAAdHWPwx4AAADMIdACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsnDnsAR8UY4+Yk35lkJHl9Vf3xIQ8J2EOfwvrpU67HkQy0Y4x7JvmuJPdMUjv/XlFVf3dAd/vsJC+vqg9tx/B1Sf6hqn7rgO7vmo0xHlNVv7Pz+4OTfFtVvegQh3WXxhifl+SfJbnfmsfJtdOn++vWp2OMf5LkXyb5+2yO372S/ERVfeRQB8Yi9On+mvbpv0hyW5KPJbk5ya9U1TsOdWALOpKBNsnzk7ykqj54ccEY435Jvn+M8fyquv0A7vNjF5svSarqVw/gPuZ6TJI7GrCq/irJWpvvM5J8R5L/mOSVSb7ncEfEAdCn++vUp5+Z5L+pqh/cWXZzkjNJfuDQBsaS9On+OvXpA5N8XVW9YM/y7x9jvKeq/t9DGtqijtw1tGOMr09yfrf5kqSqPpzkfJKnHNBd1wHVPZaq6kNV9cNV9dqquu2wx8Oy9OmR8Y1JXrq7oKr+Icn7tzN7NKZPj4ynJvnxfZb/cpLH381jOTBHcYb2IVX12v1uqKp3jDEefvH3McYNSZ6Z5N7ZNNDNSd5aVb+5s85zkvxikqcn+fh2vduq6mXb278oyVclecIY40SSj1bVT40xnpjkI1X1R9v17pvNLOPHsrkuqJL8SpInVtWrLt5XVb1k77h3l2/H85YkX5HkliQ/U1UfGmM8KckXZvN2wo3b+j9aVbePMR6X5HE7Y/yDqnr9fvc5xvjWJPdPcnuSm5L8RVX9/J79cS6bt6A+uV3vxu19CZ5cLX16NPr0j6rqln2W/12S+yb54D630Yc+PRp9+vKq+uQ+y/9RklsXvJ9DdRQD7X4H7Q5V9cs7vz4vybmqes/FBWOMbxhjPKWqfn276F5JnpHkRRdPiDHGl44xnlRVr6uqtyV52xjjE/s1z45nJfmxiw/+20b4gSTvusbte2CSR1bVD+2M+cuS3FBVP7Kz7OHZNMlPVNWbkrzpSmMcY3xvNhfgv31n2ePHGN9RVa/cWfU5SV5cVbdu1/mcJE9L8qpr3BaOL32a/n1aVW++zE0P3DurR0v6NEeiTz/lOI4xHprkq7PSyySux1EMtDdczUpjjJNJ3r3bfMmmQccYU5KLDfjYJE/fPSGq6i1jjGcned1V3teDk7xvdyajqj4xxnh1kkdfTY0dD03y3D3LHrD3VXRV/dn2QvqrMsb4tCT33m2+bZ03jjG+fIxxr6r6WJJHZvNgdOvOOn+5/f+Xq/3kJA+7zM0/XQf3wQLWS5/maPbpGOOxSd56Neuyevo0R6dPxxiflc1nU74wm2P7HVV1ZC7vOIqBdlzlel+R5NWXue2dY4wHVdX7kryxqj46436S5BHZ/wH+Pbn2BnzL3hPwcm8JZfN2zNV6fJLfvsxtb0jyJUnenOSdVfXn+6xz2f1RVb9xDePgeNCnlxyZPh1jPCDJf1lVL73iynSgTy9p36dV9TfZzsiOMf7zbD7Y98Kq+vvrrbkmxznQ3lQ7n6Lc4wNJHpDkfdlcQzPnfpLks5Is9T16+76aGmM8LJsL9O+RS28Tff411H1Akt+7zG0fSHLxWqkl9gfo0yPWp2OMeyf575P8m4O6D+52+vSI9elFVfX+McYPJ/n27P+BsXaOYqC98a5uHGM8tapek+SWMcZNtf8HGh6S5O37LL9ef53NCf43V1jvuqb+xxjfnM0F+D9eVR/fWf6cayjzwWwuXv/bfW57yGWWX+34XHLAXvr00vL2fTrGGEn+p2yua7zL6y5pRZ9eWt62T8cYX1M7H867qKo+Mo7Qt5EcxUD77jHGI6rqT/besH3V9Yntr29I8jVJXrNPjYdu3x5ZyjuyeUtm75hu2vP7pzTgGOPTk9znCvU/r6peuM/yK/2/XW9K8i1JLuxz22OTvPgaat2JSw7Yhz695Cj06XOS/Oxl3k6mL316Sec+/aIxxn+oPd+cMDbfTHEtl1Ks2pH7Htokv5TkaWOM++8uHGPcJ8k3ZXtx+rZBHzXG+Ow96z05ye8uOaCqen+SB23fkrt4PyObrzi506pj8wcFdj0jlx40LmdsT8zdBU9N8gV71rvsC5ja/FWfe4wxHrWnzqOTvH97ATssRZ/maPTpGOPbkvwfVfWXd/d9c+D0aY5En/5CkueNMfZmvmdm/xchLR25GdqqqjHGC5J89/akvC2bV2q3Z/PXTnb/qskLkzxrjHGvbL4T78Ykf1pVv3YAQ3t5kmePMT6ezQuJE0lem80F7he9crvOLdux3JzNyfaVV6j9s0meO8b4cDbX3tyY5PVJ/nDPeu8aY3xfkvdW1c/tLVJVLx1jPGOM8bXZfDfdDdk03yuubVPhrunTo9GnY4z/Kpsvbb/fGOO/3nPze6vqF+/uMbEcfXo0+rSq/p8xxs9l8yGwT2Tznbv3TPK67QuEI2EcoW9saGeM8cVJTlTV5S4eBw6ZPoX106ccxUsOVmmM8dW70/3bt0iekOT3D29UwC59CuunT9mPGdq7ybj01z9uzeZrQE4k+YWFL5YHZtCnsH76lP0ItAAAtOaSAwAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhNoAQBobfafvp2myfd+caxN0zQOewxXok857vQprN+cPp0daJPk7LkLs2ucOX0ySTJNp2bXmqbzOz/Pq7dkrd16S+6zTd3ltvO4jG3pc23t1rzv9OnVOajtTOZv65LbualxMMdg7Zbcd/r06h3HPl3zc8L1cMkBAACtCbQAALQm0AIA0JpACwBAawItAACtCbQAALQm0AIA0JpACwBAawItAACtCbQAALQm0AIA0JpACwBAawItAACtCbQAALQm0AIA0JpACwBAa6OqZhWYpmleAWhumqZx2GO4En3KcadPYf3m9KkZWgAAWjuxRJFpOrVAjfOL10qSs+cuzKp15vTJnbrLjW3uuJJlx7bkPkvuPLYlj8GSY1t6O9duye1dc5+ueTuX7NM1j02fXr/jcv7azqt3UM/1a84018MMLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGujqmYVmKZpXgFobpqmcdhjuBJ9ynGnT2H95vSpGVoAAFoTaAEAaO3EEkWm6dQCNc4vXmuJeru1zp67MKtWkpw5fXJbd73beVzGtvR2rt2S27tkL2zqrvMcWXo759Zbstbeems9Bvr0emqs//xd8zmy1tywRL0ux+B6mKEFAKA1gRYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhtVNavANE3zCkBz0zSNwx7DlehTjjt9Cus3p0/N0AIA0JpACwBAayeWKDJNpxaocX7xWkvU26119tyFWbWS5Mzpk4vXWqLebq01H4Mlx7b0MVi7JbfXOXL1luytg3o82tRe5zFY+lxbu+PSp2t+Pl3rdi5R76Ce69fwfGqGFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoLVRVbMKTNM0rwA0N03TOOwxXIk+5bjTp7B+c/rUDC0AAK0JtAAAtHZiiSJnz12YXePM6ZNJkmk6NbvWNJ2/4+e5Y7s4riVq7dY7Ltu5RL3dWkvut6W3c+3WfP7OrbdkLySXtnPNvbDmsenT67fkvltzn675HFlrLyTrfQxZQ5+aoQUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaG1U1q8A0TfMKQHPTNI3DHsOV6FOOO30K6zenT83QAgDQmkALAEBrJ5YoMk2nFqhxfvFaS9TbrXX23IVZtZLkzOmTi9daot6StfbWW+sxWPpcW7s19+mS56/tvHodxqZPr6fGZnvX/Fi+5nNkrb2QrPe5fg19aoYWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgtVFVswpM0zSvADQ3TdM47DFciT7luNOnsH5z+tQMLQAArZ1Yosg0nVqgxvkkydlzF2bXOnP65E7deWO7OK4lau3WW7JWMn+/LbnPNjUOZmxLnh9Lb+farblPj8s5ok+vvZ4+vZ4a638+XXJsa348WnOfrvlcux5maAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFobVTWrwDRN8wpAc9M0jcMew5XoU447fQrrN6dPzdACANCaQAsAQGsnlihy9tyF2TXOnD6ZJJmmU7NrTdP5nZ/n1dutZTuv3sXtXKLebq0l99vSx2DtljyuzpFrr5WsdzuXqNflsXLt1nz+rvUcWfrxaMntXPNjyFHrUzO0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGujqmYVmKZpXgFobpqmcdhjuBJ9ynGnT2H95vSpGVoAAFoTaAEAaO3EEkWm6dQCNc4vXitJzp67MKvWmdMnF6u1W++4bOcS9XZrrflcW7slj+vS58jcY7F7HJY8rselF5L1jk2fXrsOzzNrPkc8n169g3qsvB5maAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFobVTWrwDRN8wpAc9M0jcMew5XoU447fQrrN6dPzdACANDaiSWKnD13YXaNM6dPLl4rSabp1Kxa03R+sVq79da8nUuPbW69JWvt1lvyeHaw5vN3yXNkye1c+hxZ63YuUe+gHkP06fXU0KfXWytZ73NWst5Ms4Y+NUMLAEBrAi0AAK0JtAAAtCbQAgDQmkALAEBrAi0AAK0JtAAAtCbQAgDQmkALAEBrAi0AAK0JtAAAtCbQAgDQmkALAEBrAi0AAK0JtAAAtCbQAgDQ2qiqWQWmaZpXAJqbpmkc9hiuRJ9y3OlTWL85fWqGFgCA1mbP0AIAwGEyQwsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANDaicMewFExxrg5yXcmGUleX1V/fMhDAvbQp7B++pTrcSQD7Rjjnkm+K8k9k9TOv1dU1d8d0N0+O8nLq+pD2zF8XZJ/qKrfOqD7u2ZjjMdU1e/s/P7gJN9WVS86xGFdtTHGTyd5WVX90WGPhfn06f669ekY458m+cok79tz0x9V1f91CENiQfp0f936NEnGGDcl+fYk903yySQ/c3EfHwVHMtAmeX6Sl1TVBy8uGGPcL8n3jzGeX1W3H8B9fmz3xKiqXz2A+5jrMUnuaMCq+qskq22+XWOMb0jiVfrRok/3161PH5Tkx6rqbw57IBwIfbq/Vn06xrh/kudk06t/e9jjOQhH7hraMcbXJzm/23xJUlUfTnI+yVMO6K7rgOoee2OMByb57Ai0R4Y+PVIekOQDhz0IlqdPj5TvTPK/HNUwmxzNGdqHVNVr97uhqt4xxnj4xd/HGDckeWaSe2fTQDcneWtV/ebOOs9J8otJnp7k49v1bquql21v/6IkX5XkCWOME0k+WlU/NcZ4YpKPXHx7fIxx3yTfk+Rj2VwXVEl+JckTq+pVF++rql6yd9y7y7fjeUuSr0hyS7ZvGYwxnpTkC5PcluTGbf0frarbxxiPS/K4nTH+QVW9fr/7HGN8a5L7J7k9yU1J/qKqfn7P/jiXzVtQn9yud+P2vm7bb78v4LuT/OvtNnA06NMj1KdVJYAcTfr0CPTpGOPzk7yjqj6+VM01OoqB9pN3dWNV/fLOr89Lcq6q3nNxwRjjG8YYT6mqX98uuleSZyR5UVV9crvOl44xnlRVr6uqtyV52xjjE/s1z45nZTPVf8u2xokkP5DkXde4fQ9M8siq+qGdMX9Zkhuq6kd2lj08myb5iap6U5I3XWmMY4zvzeYC/LfvLHv8GOM7quqVO6s+J8mLq+rW7Tqfk+RpSV51jdtyRWOMf5Hk16rq42OMpctzePRpjkyfjm2f/uNsjsPFYPDRhe+Hu58+zZHo03+a5PwY46lJPjfJfZJ8OJvjdeuC93OojmKgveFqVhpjnEzy7t3mSzYNOsaYklxswMcmefrF5tuu85YxxrOTvO4q7+vBSd53sfm2NT4xxnh1kkdfTY0dD03y3D3LHrD3VXRV/dn2QvqrMsb4tCT33m2+bZ03jjG+fIxxr6r6WJJHZvNgdOvOOn+5/f+Xq/3kJA+7zM0/XZf5YMH2UoPPrapfvNrtoA19mqPRp9nMkv3+xT4dYzwgyfPHGM+tqk9c7baxSvo0R6JPPyOb2fN/X1Wv2db53CT/avvvSDiKgfZqp/G+IsmrL3PbO8cYD6qq9yV542VmGq5luvARSd66z/L35Nob8C1739673FtC2TzRXK3HJ/nty9z2hiRfkuTNSd5ZVX++zzqX3R9V9RvXMI5dz07yguv8v6ybPr2kdZ9W1f+65/e/HWP8WJJTSf6366nJaujTSzr36T9O8oN15w/a/acxxmvHGF9ZVf/hOmquznEOtDfV5b+u4gPZfNDhfdlcQzPnfpLks7LcB5r2vVZtjPGwbC7Qv0cuvU30+ddQ9wFJfu8yt30gycVrpZbYH1c0xvimbF5NXsuDCH3o0yPQp5dTVe8dY3zG3XFfHCh9ejT69N37HZ+q+v3t7LhAu1I33tWNY4ynbqfcbxlj3LT7tsWOhyR5+z7Lr9dfZ3OCX+lrba7rgxVjjG/O5gL8H6+di763F5xfrQ9mc/H6fp+AfMhlll/t+K7nLZIvT/KgMcZjd5Z9QTbH7YlJ/vfafE0KPenTS8vb9unYfK/l/S7Ti3d5/SUt6NNLy9v2aTbXLl/OPa93LGtzFAPtu8cYj6iqP9l7w/ZV18Vrut6Q5GuSvGafGg/dvj2ylHdk85bM3jHdtOf3T2nAMcanZ3MB9135vKp64T7Lr/T/dr0pybckubDPbY9N8uJrqHUn1/MWSVV9z95lez/pSmv69JK2fZrNJ9X/ZZI7fThmjHFjNp/Ypjd9eknnPv1PY4zP23t5wxjjnyT5v693LGtz5L6HNskvJXna2HyJ8B3GGPdJ8k3ZXpy+bdBHjTE+e896T07yu0sOqKren81s47137mdkc5H2nVbd5226Z+TSg8bljLH5ypTdBU/NZkZz12VfwFTVR5LcY4zxqD11Hp3k/d76Z2H6NP37dPvhnrH9AOeub8/ma5ToTZ+mf58m+eUkz9wet4tjuUc2X592vZ9xWZ0jN0NbVTXGeEGS796elLdl80rt9mz+2snurMELkzxrjHGvbGYabkzyp1X1awcwtJcnefYY4+PZvJA4keS12VzgftErt+vcsh3Lzdm84v3KK9T+2STPHWN8OJtrb25M8vokf7hnvXeNMb4vyXur6uf2Fqmql44xnjHG+Nokt2bzCdf3V9Urrm1T4a7p0yPVpy9L8t9tZ2VvzGam7HdcEtSfPj0afVqbr708m+R7xxi3J/m7JP8oyb+rg/v++LvdKN+HfWjGGF+c5ERVXe7iceCQ6VNYP33KUbzkYJXGGF+9neK/+PtI8oQkv394owJ26VNYP33KfszQ3k3Gpb/+cWs2n/49keQXFr5YHphBn8L66VP2I9ACANCaSw4AAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoLXZf/p2mibf+8WxNk3TOOwx3BU9Cuvv00Svwpw+nR1ok+TsuQuza5w5fTJJMk2nZteapvM7P8+rt2St3XpL7rMl6u3WOi7HYOntXLs177u1nr9r7oU1P4YsObalt3Pt1rzv1nr+rvnx6LgcgzX0qUsOAABoTaAFAKA1gRYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhtVNavANE3zCkBz0zSNwx7DXdGjsP4+TfQqzOlTM7QAALR2Yoki03RqgRrnF6+VJGfPXZhV68zpkzt117udc+stWWtvvbUeg7njSu48trVb8/m75Dmy5HE9Tn265NiWPAbHrU/XvO/W+li+5sejNY9tzY+V18MMLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGujqmYVmKZpXgFobpqmcdhjuCt6FNbfp4lehTl9aoYWAIDWBFoAAFo7sUSRs+cuzK5x5vTJxWslyTSdmlVrms4vVmu33pq3c+mxza23ZK3dektv59otef4uWStZ9hw5Ln269DFY62PI0tu5dmved2s9f9f8eHRcHkPW8HxqhhYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhNoAQBobVTVrALTNM0rAM1N0zQOewx3RY/C+vs00aswp0/N0AIA0JpACwBAayeWKDJNpxaocT5Jcvbchdm1zpw+ecfPc+vt1lpyO5estUS93VpLH4O1jm3pY7B2a+7TtZ4jx2U7l6i3ZK3dektv59qted8dl3Nkrdu5RL0umeZ6mKEFAKA1gRYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhNoAQBobVTVrALTNM0rAM1N0zQOewx3RY/C+vs00aswp0/N0AIA0JpACwBAayeWKDJNpxaocX7xWkvUW7LWbr2z5y7MrnXm9Mk7fp5bb7fWmo/Bkvtt6WOwdmved8flHFmyT9d8DNb8nLB2a953a+3TNW/ncRnbGp5PzdACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArY2qmlVgmqZ5BaC5aZrGYY/hruhRWH+fJnoV5vSpGVoAAFoTaAEAaO3EEkXOnrswu8aZ0ycXr5Uk03RqVq1pOr9Yrd16x2U7k/nburudaz7X1m7J7V36HFny/F3zObLW7dzUXufYlj7X1u64nL9LHtfj8ni0qX08jsH1MEMLAEBrAi0AAK0JtAAAtCbQAgDQmkALAEBrAi0AAK0JtAAAtCbQAgDQmkALAEBrAi0AAK0JtAAAtCbQAgDQmkALAEBrAi0AAK0JtAAAtCbQAgDQ2qiqWQWmaZpXAJqbpmkc9hjuih6F9fdpoldhTp+aoQUAoLUTSxSZplML1DifJDl77sLsWmdOn9ypO29sF8e1RK3demvezqXHNrfektu5qXF+8VodrHnfrfX8PS69sES9gzoG+vR6aujT662VrHc7N7XXmWnW0KdmaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFobVTWrwDRN8wpAc9M0jcMew13Ro7D+Pk30KszpUzO0AAC0JtACANDaiSWKnD13YXaNM6dPJkmm6dTsWtN0fufnefV2ay25nUvWWqLebq2lj8Fax7b0MVi7Ne+7tZ4jx6UXkmXHtubHyrVbc58u+Xy65u1cay8sUa/LY+X1MEMLAEBrAi0AAK0JtAAAtCbQAgDQmkALAEBrAi0AAK0JtAAAtCbQAgDQmkALAEBrAi0AAK0JtAAAtCbQAgDQmkALAEBrAi0AAK0JtAAAtDaqalaBaZrmFYDmpmkahz2Gu6JHYf19muhVmNOnZmgBAGhNoAUAoLUTSxQ5e+7C7BpnTp9MkkzTqdm1pun8zs/z6i1Za7fekvtsiXq7tRyDq7e739ZuzftuyXPE49Hh1OtyDNZuzX261ueZNT8eHZc+XcPzqRlaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1kZVzSowTdO8AtDcNE3jsMdwV/QorL9PE70Kc/rUDC0AAK2dWKLINJ1aoMb5xWstUW+31tlzF2bVSpIzp08uXmuJeru1lj4GS45tyf229Hau3XHp0yXHtnSfrvXxaFP7eByDtVvzvjsu58hatzNZ73P9Gp5PzdACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGvi3gmAAAAQ7klEQVQCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0NqpqVoFpmuYVgOamaRqHPYa7okdh/X2a6FWY06dmaAEAaG32DC0AABwmM7QAALQm0AIA0JpACwBAawItAACtCbQAALQm0AIA0JpACwBAawItAACtCbQAALQm0AIA0JpACwBAawItAACtCbQAALQm0AIA0JpACwBAawItAACtnTjsARwVY4ybk3xnkpHk9VX1x4c8JGAPfQrrp0+5Hkcy0I4x7pnku5LcM0nt/HtFVf3dAd3ts5O8vKo+tB3D1yX5h6r6rQO6v2s2xnhMVf3Ozu8PTvJtVfWiQxzWvsYYz0pyr31uOpHk5qp6wd08JBamT/fXqU+TZIzxRUn+eZKLx+w/S/ILVfWOwxsVS9Gn+2vYp/9Fkm9Kcks2LxRuSPITVfX/HerAFnQkA22S5yd5SVV98OKCMcb9knz/GOP5VXX7Adznxy42X5JU1a8ewH3M9ZgkdzRgVf1VklU2X1X95H7LxxjfmORNd/NwOBj6dH9t+nSM8cAkX1tV086ykWQaY7ysqj5waINjKfp0f5369DOT/LdV9YM7y25O8rztvyPhyF1DO8b4+iTnd5svSarqw0nOJ3nKAd11HVBd7uzB2wcOGtOnR8ZXJnnl7oKqqiQ/ub2NxvTpkfEtSX54d0FV/UOSN4wxvvhwhrS8ozhD+5Cqeu1+N1TVO8YYD7/4+xjjhiTPTHLvbBro5iRvrarf3FnnOUl+McnTk3x8u95tVfWy7e1flOSrkjxhjHEiyUer6qfGGE9M8pGq+qPtevdN8j1JPpbNdH8l+ZUkT6yqV128r6p6yd5x7y7fjuctSb4im7cOfqaqPjTGeFKSL0xyW5Ibt/V/tKpuH2M8Lsnjdsb4B1X1+v3uc4zxrUnun+T2JDcl+Yuq+vk9++NcNm9BfXK73o3b+7ptv/2+lDHGE5L8nwd5H9xt9OnR6NOPJbnPPstv3t5Gb/r0aPTpxQC715uzCbtvXfK+DstRDLSfvKsbq+qXd359XpJzVfWeiwvGGN8wxnhKVf36dtG9kjwjyYuq6pPbdb50jPGkqnpdVb0tydvGGJ/Yr3l2PCvJj1XVLdsaJ5L8QJJ3XeP2PTDJI6vqh3bG/GVJbqiqH9lZ9vBsmuQnqupNSd50pTGOMb43mwvw376z7PFjjO+oqt1ZmOckeXFV3bpd53OSPC3Jq65xW67Vl1TVjx3wfXD30Kc5En3675P8qzHGS6vqb7f385nZBJZ/veD9cDj0aY5En+47411VHx1j7PdZlZaOYqC94WpWGmOcTPLu3eZLNg06xpiSXGzAxyZ5+sXm267zljHGs5O87irv68FJ3nex+bY1PjHGeHWSR19NjR0PTfLcPcsesPdVdFX92fZC+qsyxvi0JPfebb5tnTeOMb58jHGvqvpYkkdm82B06846f7n9/5er/eQkD7vMzT9dV/HBgjHGFyT5j1ezLbSgT9O/T6vq42OMlyT5mTHGe7OZYXpgku+sqk9c7XaxWvo0/fs0ya1jjM+oneuSt7Uen+S+V7lZq3cUA+24yvW+IsmrL3PbO8cYD6qq9yV5Y1V9dMb9JMkjsv+U/nty7Q34lqq606uty70llGt7y+/xSX77Mre9IcmXZPP2xDur6s/3Weey+6OqfuMaxnE5X5vk3y1Qh3XQp5e07dOx+XDQ/5jkmVX119tl90/yP4wxfnLvEyjt6NNL2vZpkp9N8j+PMV54cf9vL514RC59O0l7xznQ3nQXD7YfSPKAJO/L5hqaOfeTJJ+VZKnv0dv3rYMxxsOyuUD/Hrn0NtHnX0PdByT5vcvc9oEkF6+VWmJ/XJPtk+aH9z7w0Jo+PRp9+u1J/s1uSKmqD44x/m02bwv/6ML3x91Lnx6BPq2qW8YYL0rybdvLM05k86Lg15L8syXv6zAdxUB7413dOMZ4alW9JsktY4ybdt+22PGQJG/fZ/n1+utsTvC/ucJ61xXYxhjfnM0F+D9eVR/fWf6cayjzwWwuXv/bfW57yGWWX+345l5y8I1Jfu56759V0qeXlnfu0xv2m3Grqlu3HxKiN316aXnnPk1tvm/2ZXtqfVN2vnqsu6MYaN89xnhEVf3J3hu2r7ouXtf1hiRfk+Q1+9R46PbtkaW8I5u3ZPaO6aY9v39KA44xPj37f4p41+dV1Qv3WX6l/7frTdl82vHCPrc9NsmLr6HWncy55GCMcWOSG2v/T2jSlz69pHOf3vsubtu73+hHn17SuU8/xRjjHtl80Ppyl4q0c+S+hzbJLyV52vY6rjuMMe6TzV/J+PUk2Tboo8YYn71nvScn+d0lB1RV70/yoDHGHQ/+Y4yRzVec3GnVMcZn7Fn2jFx60LicsXc2ZIzx1CRfsGe9y76AqaqPJLnHGONRe+o8Osn7txewH4Z/ns3XsXC06NMciT797THGd+1dODZfV/Tmu3ksLE+f5kj06Z1sv9lgyp7vkO7uyM3QVlWNMV6Q5Lu3J+Vt2bxSuz2bv3ay+1dNXpjkWduD+/Fs3l7506r6tQMY2suTPHuM8fFsXkicSPLabC7KvuiV23Vu2Y7l5mxe8V7pC8p/NslzxxgfzubamxuTvD7JH+5Z711jjO9L8t6q+pS38KvqpWOMZ4wxvjbJrdl8wvX9VfWKa9vURT1o4Vf3rIA+PRp9WlVvHmNkjPG8JH+fS98/+qaqesPdPR6WpU+PRp8md3wW5WnZbM+9kvxUVb33MMZyUIbP2RyesfkLHSeq6nIXjwOHTJ/C+ulTjuIlB6s0xvjq7TUrF38fSZ6Q5PcPb1TALn0K66dP2Y8Z2rvJuPTXP27N5mtATiT5BW+nw3roU1g/fcp+BFoAAFpzyQEAAK0JtAAAtCbQAgDQmkALAEBrAi0AAK0JtAAAtDb7T99O0+R7vzjWpmkahz2GK9GnHHf6FNZvTp/ODrRJcvbchdk1zpw+uXitJJmmU7NqTdP5xWrt1lt6O+fWW3KfbWost992a635XFu7JY/rcTlHjst2bmqv87Fy6WOwdmt+jFvr88ya+3TpsS15DNZ8rl0PlxwAANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0NqpqVoFpmuYVgOamaRqHPYYr0accd/oU1m9On5qhBQCgtRNLFJmmUwvUOL94rSQ5e+7CrFpnTp/cqbvc2OaOK7nz2Na6nUvUW/J4Jpe2deljsHZLbu+a+3TN58hx6dM1Pyes3Zr33Vofy9e8ncdlbGvoUzO0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArY2qmlVgmqZ5BaC5aZrGYY/hSvQpx50+hfWb06dmaAEAaE2gBQCgtRNLFDl77sLsGmdOn1y8VpJM06lZtabp/GK1dustWWuJege1nUvUO6ixLX2urZ0+vfZ6a97Opcc2t96S27mpcTCPlWu35vN3refImvv0uDyfrqFPzdACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArY2qmlVgmqZ5BaC5aZrGYY/hSvQpx50+hfWb06dmaAEAaE2gBQCgtRNLFJmmUwvUOL94rSQ5e+7CrFpnTp9crNZuvaW3c269JfdZcuf9tuTYltxvS2/n2q35/F2yT9d8jqx1O5P1PlYetz49Lufvmh+Pjsvz6ZqPwfUwQwsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANCaQAsAQGsCLQAArQm0AAC0JtACANDaqKpZBaZpmlcAmpumaRz2GK5En3Lc6VNYvzl9aoYWAIDWBFoAAFo7sUSRs+cuzK5x5vTJxWslyTSdmlVrms7f8fOSY5s7rk2N8zs/L7edS49t7n5b8nhuamzGtvS5tnZr3nfHpU+X7AXH4Ortjm3tltze4/I8s+ZeWPMxOGp9aoYWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaG1U1awC0zTNKwDNTdM0DnsMV6JPOe70KazfnD41QwsAQGsCLQAArZ1YosjZcxdm1zhz+uTitZJkmk7NqjVN5+/4ecmxzR3Xpsb5nZ+X286lxzZ3v+0ezzUfg7VbcnvX3Ke283DqdTkGa7fmfbfWx/I1b+dx6dM1PJ+aoQUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhtVNWsAtM0zSsAzU3TNA57DFeiTznu9Cms35w+NUMLAEBrJ5YoMk2nFqhxfvFaS9TbrXX23IVZtZLkzOmTi9daot6StfbWW/IYLHl+LL2da7fk9h6XPj0u27lEvYN6DFn6GKzdmh/j1nr+HpdeWKLeksdzU2M9z6dmaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFoTaAEAaE2gBQCgNYEWAIDWBFoAAFobVTWrwDRN8wpAc9M0jcMew5XoU447fQrrN6dPzdACANCaQAsAQGsnligyTacWqHF+8VpJcvbchVm1zpw+uVN3ubHNHVey7Nh295ljcPV2x7Z2S27v0ufIWs/fpc+RtfZCsuzY1nyurd1xOX+XHNtx2c5kvY+Va+hTM7QAALQm0AIA0JpACwBAawItAACtCbQAALQm0AIA0JpACwBAawItAACtCbQAALQm0AIA0JpACwBAawItAACtCbQAALQm0AIA0JpACwBAa6OqZhWYpmleAWhumqZx2GO4En3KcadPYf3m9KkZWgAAWhNoAQBo7cQSRc6euzC7xpnTJxevlSTTdGpWrWk6v1it3XpL1krm77cl99mmxnL7bcntTC5t69LbuXZrPn/Xeo4cl8ejJep1OQZrt+Y+XfJ5Zs3nyJqfT9d6DNbwfGqGFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoLVRVbMKTNM0rwA0N03TOOwxXIk+5bjTp7B+c/rUDC0AAK2dWKLINJ1aoMb5xWslydlzF2bVOnP65GK1dustWWuJeru1lj4Gc+stWWu33tLHYO2W3N419+maz5Ele+G4PIYsfa6tnfP36h3U8+lat3NTe51jW0OfmqEFAKA1gRYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhNoAQBoTaAFAKA1gRYAgNYEWgAAWhNoAQBobVTVrALTNM0rAM1N0zQOewxXok857vQprN+cPjVDCwBAa7NnaAEA4DCZoQUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1gRaAABaE2gBAGhNoAUAoDWBFgCA1v5/jz7lZOoD+1cAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 720x720 with 9 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"nine_states_df = pd.DataFrame({\n",
" \"Configuration\": [f\"{n}\" for n in range(1, 10)],\n",
" \"state\": nine_states,\n",
"})\n",
"g = sns.FacetGrid(\n",
" nine_states_df,\n",
" col=\"Configuration\",\n",
" col_wrap=3,\n",
" col_order=[f\"{n}\" for n in range(1, 10)],\n",
")\n",
"g.map_dataframe(draw_heatmap, xticklabels=False, yticklabels=False,\n",
" cmap=\"YlGnBu\", cbar=False, square=True,\n",
" linewidths=1, linecolor='tab:gray')\n",
"\n",
"g.fig.set_dpi(120)\n",
"g.fig.set_size_inches(6, 6)\n",
"plt.tight_layout();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Defining the Random Walk\n",
"\n",
"* We will consider a random walker that can move in one of four directions, NE, NW, SE, SW\n",
"* For each step, we need to generate two random numbers, each of which is either 1 or -1\n",
"* The pair of numbers corresponds to the four directions:\n",
" * NE step: $(1, 1)$\n",
" * NW step: $(-1, 1)$\n",
" * SE step: $(1, -1)$\n",
" * SW step: $(-1, -1)$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Algorithm for the Random Walk\n",
"\n",
"1. Select the number of steps `n` for the simulation\n",
"2. Initialize the walker's position at the origin $(x, y) = (0, 0)$\n",
"3. Initialize a list to store the walker's history, put in the origin as the first entry\n",
"4. Do the following `n` times:\n",
"\n",
" rand <- a random 0 or 1\n",
" if rand is 0\n",
" increment x by 1\n",
" else\n",
" decrement x by 1\n",
" rand <- a random 0 or 1\n",
" if rand is 0\n",
" increment y by 1\n",
" else\n",
" decrement y by 1\n",
" append point (x, y) onto end of walker's history\n",
" \n",
"5. Return the walker's history at the end of the walk"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Code for Random Walk (based on pseudocode)\n",
"\n",
"First, we implement the Random Walk as specificed in the above pseudocode.\n",
"The code is as follows:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Choose the number of steps, initialize walker_x and walker_y, and initialize\n",
"the walker history as an array of 3 columns and 101 rows (origin start + 100\n",
"steps)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"nsteps = 100\n",
"walker_x = 0\n",
"walker_y = 0\n",
"walker_history = np.zeros((101, 3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first column in walker_history is the step number, which runs from 0 to 100."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"walker_history[:, 0] = np.arange(start=0, stop=101, step=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Simulate the walker's steps with a for loop"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"for step_number in np.arange(start=1, stop=nsteps + 1, step=1):\n",
" # Generate the first random number rand for x coordinate\n",
" rand = np.random.randint(low=0, high=2)\n",
" # Conditional test\n",
" if rand == 0:\n",
" walker_x += 1\n",
" else:\n",
" walker_x -= 1\n",
"\n",
" # Generate the second random number rand for y coordinate\n",
" rand = np.random.randint(low=0, high=2)\n",
" # Conditional test\n",
" if rand == 0:\n",
" walker_y += 1\n",
" else:\n",
" walker_y -= 1\n",
" \n",
" # Store the step in the history\n",
" walker_history[step_number, [1, 2]] = (walker_x, walker_y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Convert walker_history into a pandas DataFrame"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"walker_history = pd.DataFrame(data=walker_history, columns=[\"step\", \"x\", \"y\"])\n",
"walker_history[\"step\"] = pd.to_numeric(walker_history[\"step\"], downcast=\"signed\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Check the first and last 5 rows to verify that everything is worked as expected."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>step</th>\n",
" <th>x</th>\n",
" <th>y</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>1.0</td>\n",
" <td>-1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2</td>\n",
" <td>2.0</td>\n",
" <td>-2.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>3</td>\n",
" <td>3.0</td>\n",
" <td>-3.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>4</td>\n",
" <td>2.0</td>\n",
" <td>-4.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" step x y\n",
"0 0 0.0 0.0\n",
"1 1 1.0 -1.0\n",
"2 2 2.0 -2.0\n",
"3 3 3.0 -3.0\n",
"4 4 2.0 -4.0"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"walker_history.head()"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>step</th>\n",
" <th>x</th>\n",
" <th>y</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>96</th>\n",
" <td>96</td>\n",
" <td>-4.0</td>\n",
" <td>-16.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>97</th>\n",
" <td>97</td>\n",
" <td>-5.0</td>\n",
" <td>-15.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>98</th>\n",
" <td>98</td>\n",
" <td>-6.0</td>\n",
" <td>-16.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>99</th>\n",
" <td>99</td>\n",
" <td>-5.0</td>\n",
" <td>-15.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>100</th>\n",
" <td>100</td>\n",
" <td>-4.0</td>\n",
" <td>-16.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" step x y\n",
"96 96 -4.0 -16.0\n",
"97 97 -5.0 -15.0\n",
"98 98 -6.0 -16.0\n",
"99 99 -5.0 -15.0\n",
"100 100 -4.0 -16.0"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"walker_history.tail()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"While this approach works fine, it is not the most efficient way to run the simulation.\n",
"Next we show how to make full use of `numpy` to obtain dramatic speedups in the calculation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Using `numpy`\n",
"\n",
"Like in many other situations, we can speed up our random walk by making use of the vectorization features in `numpy`.\n",
"This allows us to write code that does not contain any `for` loops.\n",
"In particular, we generate all quantities related to random number generation up front, and then perform a cumulative summation operation to obtain the walker's history.\n",
"The code is as follows:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Choose the number of steps and initialize the walker history as an array of 3 columns and 101 rows (origin start + 100 steps)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"nsteps = 100\n",
"walker_history = np.zeros((101, 3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first column in walker_history is the step number, which runs from 0 to 100."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"walker_history[:, 0] = np.arange(start=0, stop=101, step=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Generate the sequence of random numbers for the x and y coordinates for all 100 steps"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"random_walk_steps = np.random.choice(a=[1, -1], size=(nsteps, 2), replace=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Apply the cumulative sum function to random_walk_steps, summing down each column separately"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"random_walk_cumsum = np.cumsum(a=random_walk_steps, axis=0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Insert random_walk_cumsum into columns 2 and 3 of walker_history, starting at the row with index 1"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"walker_history[1:, [1, 2]] = random_walk_cumsum"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Convert walker_history into a pandas DataFrame"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"walker_history = pd.DataFrame(data=walker_history, columns=[\"step\", \"x\", \"y\"])\n",
"walker_history[\"step\"] = pd.to_numeric(walker_history[\"step\"], downcast=\"signed\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Check the first and last 5 rows to verify that everything is worked as expected."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>step</th>\n",
" <th>x</th>\n",
" <th>y</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2</td>\n",
" <td>2.0</td>\n",
" <td>0.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>3</td>\n",
" <td>1.0</td>\n",
" <td>-1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>4</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" step x y\n",
"0 0 0.0 0.0\n",
"1 1 1.0 1.0\n",
"2 2 2.0 0.0\n",
"3 3 1.0 -1.0\n",
"4 4 0.0 0.0"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"walker_history.head()"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>step</th>\n",
" <th>x</th>\n",
" <th>y</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>96</th>\n",
" <td>96</td>\n",
" <td>-12.0</td>\n",
" <td>20.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>97</th>\n",
" <td>97</td>\n",
" <td>-13.0</td>\n",
" <td>19.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>98</th>\n",
" <td>98</td>\n",
" <td>-12.0</td>\n",
" <td>20.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>99</th>\n",
" <td>99</td>\n",
" <td>-13.0</td>\n",
" <td>19.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>100</th>\n",
" <td>100</td>\n",
" <td>-12.0</td>\n",
" <td>20.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" step x y\n",
"96 96 -12.0 20.0\n",
"97 97 -13.0 19.0\n",
"98 98 -12.0 20.0\n",
"99 99 -13.0 19.0\n",
"100 100 -12.0 20.0"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"walker_history.tail()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you compare these with the previous outputs, you'll see that they're identical, so the implementation is working as expected."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Benchmarking the `numpy` version\n",
"\n",
"As always, don't just take my word that the `numpy` version is more efficient than the `for` loop version, check out the following benchmarks to verify for yourself.\n",
"\n",
"The main part of the algorithm that we're interested in is the `for` loop part, the other lines are just for setup and teardown purposes.\n",
"So, when doing the benchmarks, let's isolate and only time that part.\n",
"We do that by putting the initialization code on the same line as the `%%timeit` magic.\n",
"\n",
"We will increase the number of steps from 100 to 1000 for this simulation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### `for` loop version"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Set up the calculation with the important parameters on the first line and then time the simulation:"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"6.2 ms ± 129 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
]
}
],
"source": [
"%%timeit nsteps = 1000; walker_x = 0; walker_y = 0; walker_history = np.zeros((1001, 3)); walker_history[:, 0] = np.arange(start=0, stop=1001, step=1)\n",
"\n",
"for step_number in np.arange(start=1, stop=nsteps + 1, step=1):\n",
" rand = np.random.randint(low=0, high=2)\n",
" if rand == 0:\n",
" walker_x += 1\n",
" else:\n",
" walker_x -= 1\n",
"\n",
" rand = np.random.randint(low=0, high=2)\n",
" if rand == 0:\n",
" walker_y += 1\n",
" else:\n",
" walker_y -= 1\n",
" \n",
" walker_history[step_number, [1, 2]] = (walker_x, walker_y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### `numpy` vectorized version"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Set up the calculation with the important parameters on the first line and then time the simulation:"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"40 µs ± 187 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
]
}
],
"source": [
"%%timeit nsteps = 1000; walker_history = np.zeros((1001, 3)); walker_history[:, 0] = np.arange(start=0, stop=1001, step=1)\n",
"\n",
"random_walk_steps = np.random.choice(a=[1, -1], size=(nsteps, 2), replace=True)\n",
"random_walk_cumsum = np.cumsum(a=random_walk_steps, axis=0)\n",
"walker_history[1:, [1, 2]] = random_walk_cumsum"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Benchmark conclusions\n",
"\n",
"If we take the ratio, then we are able to determine how many times faster the numpy version is when compared with the `for` loop version."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The numpy version is 174 times faster than the \"for\" loop version\n"
]
}
],
"source": [
"print(f\"The numpy version is {round(5.13E-3 / 29.5E-6, 0):.0f} times \"\n",
" \"faster than the \\\"for\\\" loop version\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For this simulation, the vectorized implementation results in a 174x speedup, which is an excellent improvement."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Animating the Random Walk\n",
"\n",
"The random walk is an excellent example of a simulation that you want to animate.\n",
"We can use matplotlib to create a series of still plots that are combined into a single animated movie.\n",
"The following code will generate the movie:"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"# Initialize matplotlib viewing window\n",
"fig, ax = plt.subplots(dpi=150, figsize=(4, 4))\n",
"\n",
"# Label the horizontal and vertical axes as x and y, then give animation a\n",
"# title\n",
"ax.set_xlabel(\"x\")\n",
"ax.set_ylabel(\"y\")\n",
"ax.set_title(\"2D Random Walk: 100 steps\")\n",
"\n",
"# Initialize the walker as a point and its history as a line\n",
"point, = ax.plot([], [], 'o', color=\"tab:blue\")\n",
"line, = ax.plot([], [], '-', color=\"tab:orange\", lw=2)\n",
"\n",
"# Set the viewing window to a square region centered about the origin.\n",
"# The distance from the center to the window edge is set by the absolute\n",
"# maximum value in walker_history.\n",
"max_xy = np.max(np.abs(walker_history[[\"x\", \"y\"]].values))\n",
"ax.set_xlim([-max_xy, max_xy])\n",
"ax.set_ylim([-max_xy, max_xy])\n",
"\n",
"# Show the grid after by placing a tick at every integer value along the x and\n",
"# y directions\n",
"xticks_start, xticks_end = ax.get_xlim()\n",
"yticks_start, yticks_end = ax.get_ylim()\n",
"ax.xaxis.set_ticks(\n",
" np.linspace(xticks_start, xticks_end, np.int((xticks_end - xticks_start) + 1)),\n",
" minor=False,\n",
")\n",
"ax.yaxis.set_ticks(\n",
" np.linspace(yticks_start, yticks_end, np.int((xticks_end - xticks_start) + 1)),\n",
" minor=False,\n",
")\n",
"ax.axes.grid(True, linestyle=\"-\", linewidth=1, color=\"tab:gray\", which=\"major\")\n",
"\n",
"# After generating the grid, remove the numerical and text labels\n",
"ax.tick_params(\n",
" labelbottom = False, labelleft=False, bottom=False, left=False\n",
")\n",
"\n",
"# Running the above commands creates a blank matplotlib window, so let's get\n",
"# rid of that\n",
"plt.close()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The initialization function, clears out the data in the point and line objects."
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"def init():\n",
" point.set_data([], [])\n",
" line.set_data([], [])\n",
" return (point, line,)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The animation function. Input i is the frame number of the animation, and is to be used for referencing how the data changes over time."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"def animate(i):\n",
" # Get walker's position at step i and set as the underlying data for the\n",
" # point object\n",
" x = walker_history.loc[i, \"x\"]\n",
" y = walker_history.loc[i, \"y\"]\n",
" point.set_data(x, y)\n",
"\n",
" # Get the walker's position up through step i and set as the underlying\n",
" # data for the line object\n",
" xline = walker_history.loc[:i, \"x\"]\n",
" yline = walker_history.loc[:i, \"y\"]\n",
" line.set_data(xline, yline)\n",
"\n",
" return (point, line,)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use animation.FuncAnimation to put the animation together.\n",
"\n",
"* `frames` controls the number of frames in the movie.\n",
"* `interval` controls the delay in milliseconds inbetween each frame\n",
"* `blit` optimizes the animation size by only storing the changes between frames instead of as a series of full plots"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"anim = animation.FuncAnimation(\n",
" fig=fig,\n",
" func=animate,\n",
" frames=101,\n",
" init_func=init,\n",
" interval=100,\n",
" blit=True\n",
");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we've obtained the animation object `anim`, use the method `to_html5_video()` to render the movie in a format that can be run natively in your web-browser.\n",
"We wrap the output in the `HTML()` function (imported from `IPython.display`) to tell the Jupyter notebook to render the output HTML instead of just printing it as text."
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"#HTML(anim.to_html5_video())"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}