Notebooks/RandomWalk.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PYTHON VERSION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "print(sys.version)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Escape times \n",
    "\n",
    "Random walks on $\\mathbb{Z}$ and $\\mathbb{Z}^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.random import choice\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams['figure.figsize'] = (10, 6)  # default figure size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Symmetric walk on $\\mathbb{Z}$\n",
    "\n",
    "For the simple symmetric walk on $\\mathbb{Z}$ the well-known `simple_markov_chain_lib` library will not be necessary here, a chain from scratch will be implemented.\n",
    "\n",
    "To implement a Markov chain (at least) 2 things are necessary:\n",
    "\n",
    "1. A variable to store the current state\n",
    "2. A \"function\" to give us the next state, which can depend on the present state and time.\n",
    "\n",
    "In the following code cell the simple symmetric random walk is implemented for iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1\n",
      "0\n",
      "-1\n",
      "-2\n",
      "-1\n",
      "0\n",
      "-1\n",
      "-2\n",
      "-3\n",
      "-4\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(2017)  # for reproducibility\n",
    "\n",
    "x = 0  # current state of MC\n",
    "steps = [-1, 1]\n",
    "probs = [0.5, 0.5]\n",
    "\n",
    "for n in range(10):\n",
    "    x += choice(steps, p=probs)  # select 1 step at random\n",
    "    print(x)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable `x` stores the existing state starting from 0. In each iteration one of the 2 possible steps (left or right) is added with equal probability.\n",
    "\n",
    "The function `choice` picks one of the 2 possible steps (`steps`) with probability (`p=probs`). For the case where all possibilities are equally likely, `p` can be omitted but is added here for completeness. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# help on 'choice' function\n",
    "choice?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python - Generators\n",
    "\n",
    "For convenience (and less writing) the walk is implemented as a `generator` function. \n",
    "A generator function works like a regular function except that it \"remembers\" its previous state.\n",
    "\n",
    "The previous example using a generator is written likewise:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1\n",
      "0\n",
      "-1\n",
      "-2\n",
      "-1\n",
      "0\n",
      "-1\n",
      "-2\n",
      "-3\n",
      "-4\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(2017)  # same results are previous run\n",
    "\n",
    "# random walk generator:\n",
    "def rand_walk_Z(start = 0, probs = (0.5, 0.5)):\n",
    "    x = start\n",
    "    steps = [-1, 1]  # probs order should be (left, right)\n",
    "    while True:  # continue for ever in needed\n",
    "        x += choice(steps, p=probs)\n",
    "        yield x  # it's a generator!\n",
    "\n",
    "walker = rand_walk_Z()  # start a new markov chain\n",
    "for n in range(10):\n",
    "    x = next(walker)  # take 1 step!\n",
    "    print(x)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The command\n",
    "\n",
    "`walker = rand_walk_Z()` \n",
    "\n",
    "creates a new chain. The chain is in the initial state and waits to move on.\n",
    "When we give the command\n",
    "\n",
    "`x = next(walker)`\n",
    "\n",
    "Python runs the code until it finds the command **`yield`** and returns the variable x.\n",
    "\n",
    "A normal function would stop and \"disappear\" at this point, the next time you call it it would start over. But because of the `yield` keyword Python understands that the `rand_walk_Z` function is not a regular function but a *function generator*!\n",
    "\n",
    "Therefore, the next time the function is called with `next`, the execution will resume from `yield` (the point where it had stopped). Due to the `yield` definition inside of an *infinite loop* (`while true`) Python will continue to return values no matter how many times the chain is advanced (`next`) (that's why it's called a generator)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Application: Escape Time from Interval\n",
    "\n",
    "In this application simulations will be done for the simple symmetric walk to exit an interval $[a,b]$.\n",
    "\n",
    "\n",
    "The following code cell implements this case and visualizes the result with a time vs position graph for 3 different walks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXl8Y2d97/9+JFmLF8lje7zb8tg6yjKZJGQjC52QlSTlxXKBFnp/KVDa0L4ufd3ee39cLr+2QG9Ly4VcCiUQGkgIBVqWsIRCyMISIIWQjWQySWYk2WPZHo9nPJ6xZFu21uf3h3TOyLJsS7I2e5736+WXdM55znO+Oj7nfM6zfL9fIaVEoVAoFIpCMdXaAIVCoVBsL5RwKBQKhaIolHAoFAqFoiiUcCgUCoWiKJRwKBQKhaIolHAoFAqFoiiUcCgUCoWiKJRwKBQKhaIolHAoFAqFoigstTagEnR0dMihoaGS9l1aWqKpqam8BpUBZVdxKLuKQ9lVHDvRrmefffaklHJ3QYWllDvu79JLL5Wl8rOf/azkfSuJsqs4lF3Foewqjp1oF/CMLPAZq7qqFAqFQlEUSjgUCoVCURRKOBQKhUJRFEo4FAqFQlEUSjgUCoVCURRVEQ4hxH1CiBNCiINZ69qEEI8JIfyZz13r7PvOTBm/EOKd1bBXoVAoFOtTrRbH/cAtOev+F/ATKaUG/CSzvAohRBvwYeDVwBXAh9cTGIVCoVBUh6oIh5TyF8CpnNVvBL6c+f5l4E15dn0d8JiU8pSU8jTwGGsFSKFQbBO+euCrzK/MAxBNRPnic18kJVM1tkpRLEJWKee4EGII+IGU8oLM8ryUsjVr+2kp5a6cff5fwC6l/LvM8l8Dy1LKO/PUfwdwB0BnZ+el999/f0l2xuNxGhoaStq3kii7ikPZVRzVsGtqeYo7XryD9wy8h7f0vIXH5x7n46Mf5+/P/Xsudl5cM7tKYSfaddtttz0rpbyskLL1HnJE5FmXV+mklPcA9wDs3bu35JAjwWAQt9td0r6VRNlVHMqu4qiGXYcnDwNw0nSSoaEhTodPA7BgX2C9+/VsPl+lUC27ajmr6rgQogcg83kiT5kpYCBruR+YroJtCoWizPhD/vRnOP0ZCAdWrVdsH2opHN8H9FlS7wQezFPmEeBmIcSuzKD4zZl1CoVim6ELxGh4lGQqiS/kAzA+FduHak3H/Tfg18A5QogpIcR7gI8BNwkh/MBNmWWEEJcJIb4IIKU8Bfwt8HTm739n1ikUim2GP+zHIixEk1EOzR/i6NJRLMJCIBxQA+TbjKqMcUgp37HOphvylH0G+OOs5fuA+ypkmkKhqAIriRUmFifY372fx489zo8mfwTAa7pfw+PHHudY5Bh9TX01tlJRKMpzXKFQVJyxhTFSMsXrBl6HQPDQ5EMA3Dp4K6C6q7YbSjgUCkXF0cc39rXtY6B5gJnIDA6Lg2t7rl21XbE9UMKhUCgqjj/kx2qyMtg8iNflBcDj9OC0Oult7FXCsc1QwqFQKCqOL+Rj2DmMxWRBc2kAqz6VcGwvlHAoFIqKEwgF1giG3vLQXBpjC2PEU/Ga2acoDiUcCoWiooRiIWaWZwyhuKTjEgabB7mq6yogLRyJVILgQrCWZiqKoN5DjigUim1OIJT2EPe4PAB0Ojp5+LaHje16C8Qf8htlFPWNanEoFIqKoo9faE4t7/bhlmHMwqzGObYRSjgUCkVF8Yf8NDc009PYk3e71WzF3eJWvhzbCCUcCoWiovjDfjSnhhD5gl2n0ZyaEfxQUf8o4VAoFBVDSlnQ2IXm0phanCKSiFTJMsVWUMKhUCgqxuzKLKFYyJhRtR5elxeJZDQ8WiXLFFtBCYdCoagYxsC4K//AuI7eIlED5NsDJRwKhaJiFCocA00D2M12JRzbBCUcCoWiYvhCPjrsHeyy7dqwnNlkZsQ5ooRjm6CEQ6FQVIzsUCOboWJWbR+UcCgUioqQkikC4cCmA+M6mktjdmWW+eh8hS1TbBUlHAqFoiJMLk6yklzB4ywsjEh26BFFfVMz4RBCnCOEeD7rLyyE+IucMq8VQoSyynyoVvYqFIri0B36Cu6qcirh2C7ULMihlPIwcDGAEMIMHAW+m6foL6WUr6+mbQqFYuvkBjfcjE5HJ06rU3mQbwPqpavqBmBUSqniKisUOwRfyMdA0wCNlsaCygsh0qFHVIuj7qkX4Xg78G/rbLtKCPGCEOJHQoi91TRKoVCUTilh0vWZVVLKClmlKAc1z8chhLACbwA+mGfzc4BbSrkohLgN+B6Qt8NUCHEHcAdAZ2cnwWBpjZdYLFbyvpVE2VUcyq7iKLddsVSMI+EjXN58eVH1tifbWYgv8Kz/WXbbdp8156tcVMuumgsHcCvwnJTyeO4GKWU46/tDQojPCSE6pJQn85S9B7gHYO/evdLtdpdkTDAYpNR9K4myqziUXcVRbrsOzx8mRYrLBi/DPVh4va9ufDWfHf8skZYI7h73WXO+ykW17KqHrqp3sE43lRCiW2RiMQshriBt71wVbVMoFCVQaKiRXNSU3O1BTVscQohG4CbgvVnr/hRASvl54K3AnwkhEsAy8HapOj8VirrHF/JhMVlwNxf39uuyuuhydBkzshT1SU1bHFLKiJSyXUoZylr3+YxoIKW8S0q5V0p5kZTySinlr2pnrUJR/zww9gAzkZm822aXZ/nm6DerMvAcCAXY07IHq9la9L6aS1PZAOuceuiqUigUZeBY5BgfeuZDfH3063m3f33063zk2Y8wHZmuuC161r9S8Dg9jIZHSaaSZbZKUS6UcCgUOwR9XGC98QG9+6fS4wdL8SWOLh0tenxDx9vqJZaKMbE4UWbLFOVCCYdCsUPYTDj07p9KdwMFwmmBKlU4jNAjyoO8blHCoVDsEHTBmFqaYim+tGrbSmLFeIOvdIuj1BlVOsPOYQRCzayqY5RwKBQ7BH/IT4OpAWBN7u7R8CgSSYOpoeIzlnwhHw6Lg76mvpL2d1gcDDYPKuGoY5RwKBQ7gEQqwWh4lP09+4G1rQq922d/z35GF0aJp+IVs8Uf8uNxejCJ0h8vKqlTfaOEQ6HYAUwuThJLxbiu97q8ubv9IT9Wk5Ub+m4gkUoQXKhcWIpisv6th+bSCC4GiaViZbJKUU6UcCgUOwB9wPvc1nPxOD1rBpb9IT8jzhHObT3XWK4EcytzzEXnCs76tx6aSyMlU0wsq5lV9YgSDoViB+AP+REIhluG8bg8eVscmktjT8sezMJcMeHQ6y02Km4ueoslGKm/QIIKJRwKxY4gEA4w2DyI3WJHc2mcXDnJ6ehpAEKxEMeXj6O5NGxmG+5md8Wmuhozqkp0/tNxN7tpMDUwvjxeBqsU5UYJh0KxA/CFfEb3kP6Z69ehtwKKCemRkinuP3y/IUK5jC+M850j3zGW/SE/u2y76LB3lPZDMlhMFoadw4xHxrdUj6IyKOFQKLY5uo9GtjDAGcHQp9/qgqK5NKYWp4gkIpvWfWj+EB9/4eM8OP5g3u1f8X2Fv3r6r1iILaSPmQk1kglqvSU0p0ZwWXVV1SNKOBSKbc7YwhgpmTIEY7d9Ny6ry2hV+EN+mhua6XZ0A2kBkUjGwmOb1r1ZmBK9yysQDiClNMZSyoHX5WU2NmuIkqJ+UMKhUGxz9Ie63qIQQqC5NOOh7wv5VrUC9JZJId1V2eKTi5RyVRiT6cg0kUSkbMJhtJxU6JG6QwmHQrHN0T3GB5sHjXWaU8MfTufu9odXtwIGmgby+nqsVzekWxQpmVq1bXZllnAsbJTbaqiRXHSBU46A9YcSDoVim+MP+Rl2DmMxncnLprk0FuOLHDh1gHAsvOphbjaZGXGOFBR6xB9OOw6uJFeYWppac1wAq8lKIBQwWh/lEo7exl4cJocSjjpECYdCsc3xh/xrHO70h/cPJ364ajl7+2ZdVeFYmJnIDNf2XAus7drSl/f37McX8uEP+elp7KG5obn0H5OFEIKhxiElHHWIEg6FYhsTjoWZWZ5Z4zehd/M8PPkwsFY4PE4PsyuzzEfn161bD49+y8At6eWcFoo/5KfD3sHluy9nPjbPb078pmytDR23w40/5K9K1kJF4dRcOIQQ40KIF4UQzwshnsmzXQgh/kkIERBCHBBCXFILOxWKemS93Bcuq4tuRzcnV07SYe9gl23Xqu3eVu+q/fOhv+lf2H4h/U39a9789ZhUel0nV05uOdRILkONQ8zH5jm5crKs9Sq2Rs2FI8N1UsqLpZSX5dl2K6Bl/u4A7q6qZQpFHeObT3cX5QvxkevXkY3eQtH3z4c/5KfJ0kRvY++arq1kKkkgHMDr8q5q7Ww11EguQ44hwxZF/VAvwrERbwT+RaZ5EmgVQvTU2iiFoh4IhAPGwz0XXTDyhf/odHTibHBuONVV98nQp/cGF4LEkulotVNLU6wkV/A4PbTZ22i3ta97rK3gbnQbtijqh3oQDgk8KoR4VghxR57tfcBk1vJUZp1Cse2RUvIvvn9hJjKTd/vRpaN81f9Vo4//tyd/y6NTjxrbsx/uuWR7iueii0H2A/nHUz/m2dlnDbuynfm8Li8JmWB8YTx93Izg6N1UmkvDLMwMO4eL+v2b0drQSru9fZXAfWv0W4bzYiKV4J5X7jGmBSuqg2XzIhXnGinltBCiE3hMCHFISvmLrO35YhesGSnLiM4dAJ2dnQSDpYUqiMViJe9bSZRdxbFd7Dq6cpSPvfAxJmYnuL3/9jXl7524lweOPcC58lx223Zz5yt3MhoZxZtIP7APnT7ENbuuyftb+2J9eJu8DCQG8m7vNnXz+NzjjI+PE4/H+fDBDzPoGORj532MU7FTzMfmaU+2EwwGaYo0AfCr0V9h67Dx9NGnAbDOWwkuBHl106tpTjUzM5VfAEslFosxYB3g4ImDBINBlhJLfPjZD3Nb5238+Z4/50D4AJ965VMkF5Pc1nlbWY+9mV3b4fqqFDUXDinldObzhBDiu8AVQLZwTAEDWcv9wHSeeu4B7gHYu3evdLvdJdkTDAYpdd9Kouwqju1il38q/SY9y2xee09MnAAg0hLB3eNm8oVJwokwzd3NSClZSCxwSf8lefd14+Z73u+ta8sl8Uv44YkfYu+0E54Oczp+GmESuN1ujs4cBeDVe16Nu9NNT7IHy0ELpxtO43a7mZ2eZaBpgHOGzwHgve73ln5SNiAYDLKvax8PjD3AwOAAB+YOADCTmsHtdvMf/v8A4JT5VFX/39vl+qoUNe2qEkI0CSFa9O/AzcDBnGLfB/4wM7vqSiAkpTxWZVMVioqgdxWt51ORHeH2dPS0MbvIH/Jv2eFO388X8hnBBE9FT3Fy5eSZujNjFlazlT3OPavsLffU2/XwurwsJ5eZWppaFQIlO+SJGgOpLrVucXQB3830z1qAf5VSPiyE+FMAKeXngYeA24AAEAHeXSNbFYqyoz/4JhYnWEmsYLfYjW0LsQWORdLvSIFQYNXDMRAKkJRJoPQBaY/zTEiPSORMpFw9fEi7vZ02e5uxXnNqvHDqBWLJGMGFIDf13VTScYu2Myv0iH4OFuILzCzPGNOJfSEfUsqyROVVbE5NhUNKOQZclGf957O+S+C/VNMuhaJa+MN+bGYb0WSUsYUxzt91vrFNfyjazDb84TMPTZvZhj/kJymTax7uxdBqa6XT0UkgFCC+HDfsCIQCBMKBvE6FD00+xMHTB0nKZNVaHNkCFwgHDDt1b3Wb2cZ8bJ656NyW84AoCqMeZlUpFGcl+pv7eiE9dKG4tudaRsOjHJo/hMvq4qK2i4yH5lanv3pdXnwhH+ORcfa17aPN1sah+UMEQoE1znz68o8mfgSULybVZjQ1NBkOiL6QzzhfTxx7goX4wrrnT1E5lHAoFDVibGGMpExyfd/1WE3WNf30vpCPRksj+3v2E01GefzY42guLT2NNuxnNDy65Ye3x+lhNDzK+PI4mjNd9y9mfsFycjlvfCuAh6cexmKy4G6u3uCw5tJ4avYpTkdPc2nHpXQ5unh4Kh1O5bbB9GyqQoI2KsqDEg6FokboQnFe63kMO4fXD+mRedOfW5kzHu7LiWWWk8tbDvGhuTRiqRiRZDqPhsfpYW5lDljrBd7X1IfD4mBuZY49LXuwmq1bOnYx5NqluTRj+dWdr6bN1qYGyKuIEg6Fokb4Q37jzT3XGU+fMaQ5NUacI4iMO1O2kMDWQ3xk17Wmbufquk3CZKwrd0yqzdAdDfVj662fLkcXLqurqDzqiq2jhEOhqBH+kN94c/e6vBxfPk4oFgJgLjrHfGwezaXhsDgYaE67MnlcnlVikftwL5Zh57AhSvqbPEB/Uz9NDU1ryuvbt3rcYtHHctpsbbTb28+EU8n6zJdsSlEZlHAoFDUie3BbF4PsdK/AqpAf+nJzQzO9jb3rPtyLwWFxMNg8SIe1A5fVtWFgxGw7slsA1WCoZQiLsOQ9H5AWluXEMtNLa3yDFRWg1n4cCsVZyWJ8kenING8bfhtw5o3aH/Jz6e5LDQHRH4y3DtyKzWzDZXUB8OY9bzZaClvlTUNv4sTptId6c0Mzb3C/gdd0vyZv2f09+3l06lFe1f6qshy7UKxmK2/e82b2te0D0i2eq7qu4oa+G4AzQuYL+ehv7q+qbWcjSjgUihqQKwx65rxsT/F2Wzvt9nTU2VsHb+XWwVuN/f/L3vK5Nr33/Peuim/0sVd/bN2yQy1DfPX6r5bt2MXwN5f9jfHdarZy77X3GssjzhEgfd6u77u+6radbaiuKoWiBugCoQuHEALNqRlRYH0hX9lzW+xk9O67jRJTKcqHEg6Fogb4w34cFgd9TWcyBHhcHnwhHymZSntuV8nBbqegOzMqKo8SDoWiBvhDfjxODyZx5hb0uryEY2Gen3ue5cTWfTTONjwuD0fCR4xkU4rKoYRDoagB2UmSdPTlhyYeWrWsKAzNpZGQCYKL9ZcnY6ehhEOhqDJzK3Ocip5a06LQfSMemXwEODPgqygM/XwqD/LKo4RDoagyej987uB3mz3t3DYXnaO3sZfmhuZamLdt2dOyB7MwK+GoAko4FIoqY0zFzRPZ1nCwU+MbRWM1W3G3uJVwVAElHApFlfGH/Oyy7cqbOyLXk1xRHF6XVwlHFVDCoVBUGT14Yb5sdbkxmBTFoTk1JpcmiSQimxdWlIwSDoWiimzmo3F199Vcvvtyruy8ssqW7Qy01vR5HQ2P1tiSnY0SDoWiiszGZokkIusKR09jD1++7svsduyusmU7A72rzzevHAErSc2EQwgxIIT4mRDiFSHES0KI/5qnzGuFECEhxPOZvw/VwlaFolwciRwBVFdUpehv6sdutqvQIxWmlkEOE8D/kFI+J4RoAZ4VQjwmpXw5p9wvpZSvr4F9CkXZCUbSzmlKOCqD2WTG4/So0CMVpmYtDinlMSnlc5nvC8ArQN/GeykU25vx5XEjEq6iMnhcHjWzqsIULBxCiKuFEH8ghPhD/a9cRgghhoBXAb/Js/kqIcQLQogfCSH2luuYCkUtGI+Mq9ZGhdFcGidXTnI6errWpuxYCuqqEkJ8BRgBngeSmdUS+JetGiCEaAa+DfyFlDKcs/k5wC2lXBRC3AZ8D8h71wkh7gDuAOjs7FyVX6AYYrFYyftWEmVXcdSjXYlUgqmVKS7n8rqzrR7PF5Rml2slnezqCd8TXOi8sBJm7ajzVQqFjnFcBpwvpZTlPLgQooG0aHxNSvmd3O3ZQiKlfEgI8TkhRIeU8mSesvcA9wDs3btXut3ukmwKBoOUum8lUXYVRz3aFQgFSMgElw5eWne21eP5gtLscnQ64DCE7eGK/aaddL5KodCuqoNAdzkPLNLeT/cCr0gpP7lOme5MOYQQV5C2d66cdigU1cJI3pQn1IiifOy278ZpdapxjgpSaIujA3hZCPEUENVXSinfsIVjXwPcDrwohHg+s+7/AwYzdX8eeCvwZ0KIBLAMvL3crR6Folr4Q35MmBh2DtfalB2NEEKFHqkwhQrHR8p9YCnlE8DamAury9wF3FXuYxdCPBXnroN3cbv39rwxhXL5ydGfkJIpbuq/adOyp6On+dLhL/G+ve/DarYWZVcyleSzL32Wt428jZ7GHuaj89x76F7ed8H7sJltRdWlqC7+sJ8+e5/6P1UBzanx7xP/jpQyb2gXxdYoqKtKSvlz4BDQkvl7JbNux3Jg7gBfOPQFHp58uKDynzn4GT794qcLKvvI5CN88dAX+e3cb4u263DoMJ9/5fP8e/DfgbRg3Xv4Xp6efbrouhTVxRfy4W6sv37xnYjm0liMLzKzPFNrU3YkBQmHEOL3gKeAtwG/B/xGCPHWShpWa/RmbiHN3XgqztjCGMHFINFkdNPyep2lhEXQQ3Ln2qdCLNQ3kUSEqcUphhxDtTblrECf8qwcAStDoYPjfwlcLqV8p5TyD4ErgL+unFm1pxjhGF8YJ5FKkJRJxsJjm9cdTtdZSlgE/UYw7NtCXYrqMRYeQyIZahyqtSlnBXpYev1FS1FeChUOk5TyRNbyXBH7bkv0B7I/5Gez8fjsi1Pfbz2klEWJ0hq7MvuMLYwRT8XXCImiPtH/T26H6qqqBi6riy5Hl7ovKkShD/+HhRCPCCHeJYR4F/BD4KHKmVVbpJT4Qj4cZgdLiSWmI9MblveFfJiFmQZTA/75jS/U2ZVZQrEQDrOjIFHKxR/24zA7SKQSPH/yeeZW5nCYHQTCAZKp5OYVKGqCP+THZrbRY++ptSlnDZpLU11VFaLQwfH3k3auuxC4CLhHSvmBShpWS2ZXZgnHwlzXex2weXPXH/LjbnEz3DK8aYtDfwO6rve6gkQpm3AszExkxrDrocmHjLqiySiTS5MF16WoLoFQgBHnCGZhrrUpZw2aS2MsPEYilai1KTuOgrubpJTfllL+dynlf5NSfreSRtUa/eF+6+Ctq5bXLR/2ozm1goKr5dZdTB+sPo5x88DNmIWZRyYfAeC2wduKrktRXfSsf4rqobk0YqkYk4vqharcbCgcQognMp8LQohw1t+CECI3rtSOQW/eXtJxCd2N3Rs2d/XZMl6XF6/Ly7HIMRZiCxvW3WHv4IrdV6w6ViHoorN3117czW7mY/O4rC6u7LoSgVDN8jplPjrP7MqsCm5YZbwuL6BmVlWCDYVDSvmazGeLlNKZ9dcipXRWx8Tq4w/56bB3sMu2C82pbdj9NBoeRSLRXJrxYNhohlMglE4b2mJtobuxu6jBO3/IT6Olkd7G3lW5qRstjfQ396uBwDpFvx6UcFSX4ZZhBELdFxWgUD+OrxSybqegP9xh835S/aL0uDzGPutdqEa+6UyXxWailO9YmktDCGEcS3+rKrYuRfXQfWyUcFQXu8XOYPOgmqpeAQod41iVB0MIYQEuLb85tScpkwTCAeOB7HV5iafiBBfToYrHF8b5zMHPGLOh/CE/drOdgaYBeht7abQ0Gk1jKSWfe+lzxtjD5OIkK8mVVQ/9sXB6Wm0+Hhx/kJ9P/9yoyx/ynxGKTB0ep8dYDi4EiSVjADwz+wxf83+tvCenBGaXZ/nEC59Y9zduR5KpJP944B+ZWpzKu90f8nP3S3efuUbCfpwNTrocXdU0U0H6HsvuqvrR5I94dOrRGlq0M9hsjOODQogF4MLs8Q3gOPBgVSysMsejx1lJrhgPZN2RSG9FfOfId7j75buNGUz+kD89W8ZkNloCulAcixzjrpfu4ttHvp0um2kRZLdm4qk4E4sTeW355IFPcs8r9wBwOn6a+di8Ydfluy/nNd2v4Xd6fgdI3yBJmWRsIe2A+GXfl/n4Cx+v+QP7oYmH+NLhL3Hw1MGa2lFORsOjfOHQF/h+8Pt5t39z9Jt85qXPMLsyC6SvEY/Lo2Im1QCPy8PE4gQriRUgHRrosy99tsZWbX82G+P4ByllC/CJnPGNdinlB6tkY1U5EjkCgLc1/WY/3DKMSZjWOO1lf2Z3QWjO9NxxKeUZr+6MkOifI84RYK0oZaMPqPrDaV+P8eXxdP2ZY7XaWrln/z30NfXlrcsf8m8oStUi25Fyp7DZb9K7RnQ/neyWoqK6eF1eUjLF2MIYK4kVJhYnOBI+YrTMFaWxWYvj3MzXbwkhLsn9q4J9VSe4nO6S0h/udosdd7P7TEyoLE/tfLNlNJfGfGyekysnDWdAfR9fyMdA0wBNDU3GMbJFKRt93WJ8kWORYwQjabvWewANtQxhMVnwh/wsJ5aNKYi1nlGSe952AvqYxXrCkf2bjy8fZyG+oMY3akT2uOPYwhgpmSIhE0bXs6I0Ngur/t9Jp2P9v3m2SeD6sltUY4KRIANNAzRaGo11mkvj8Pxh4yEO6QvRSMyTIxyQfivV30xnV2aZj84bXRY6NrNtlShlk73OH/YzvjxOu62dNntbXrsbTA1pB8SQ34iLBJlWzkBJp2LLpGRqTWtrJ6C3KPSgltlh0k+unORU9BSQ6aLKGoNSVJ/B5kGsJiv+kB+RlcUht6dAURwbCoeU8o7M53XVMaf2HIkcQWtbfUFpLo3Hph7jxVMvAuCwpMOF6G/R2a0A/bt/Pi0sDouD5cQyL59+mfGFcW7su3FN3YfmD62xQw8tspxcJhAKMB4Z3/RC11waz518zrBLt7NWHF06ynJyGYfFYXTf7YR+/uz/65GFI5zbeu6qbZA+94FQAJ8z/b/QBURRXSwmC8POYeOeaDA1rIoXpyiNQqfjvk0I0ZL5/ldCiO8IIV5VWdOqTywZ4+jK0TUPaM2lIZFGbo7re69nfGGcl0+/jNPqZLd9t1G2zd5Gu62dQ/OHGAuPcX1vulH2yNQjJGUyb92Ti5MsJ5ZXrfeH/Jy/63y6Hd0cmj9EcDm4qXB4nB6ORY7x/NzzWE1Wrum6pqZdRPqxr++93ui+2+4sxZeYWpoy/q+551dvWV3fez2BcIDD84fpdHTSamutuq2KNPqEFX/Iz7BzGHeLe0d1ndaCQqfj/rWUckEI8RrgdcCXgc9XzqzacGThCClSax/uGb+LR6cepdHSyLU915KQCX46/VM0p7bmLdrj8vD4sceJpWK8pvs1OBucxhTA9UQpOxy7/kakudJhTJ6YeYJoKrqpcOgD+o/BziF/AAAgAElEQVROPcqIc4RzWs/JK0rVwgivMpAJ3bID/Ez0bqqb+m/CYrKs6YLzhXzssu3iys4rWUmu8MuZX6pQIzVGc2rMLM9w4NQBNKe2auajojQKFQ497OrvAndLKR8Eist5mgchxC1CiMNCiIAQ4n/l2W4TQnwjs/03QoihrR5zI/S3kNwHtN5PGoqF0Fwa57SeA0AoFso7WO11eQnFQkZdmksjFAthMVlwN7vXlM0+NmAMqOphTLLr2gj9AaXb5XV514hSNQmEAvQ39XNx+8XG8nZH/w3ntZ7HSMvImi4PfQaV/r/SrxlF7dBa194Xk0uTRBKRGlu2fSlUOI4KIf6ZdPa/h4QQtiL2zYsQwgx8FrgVOB94hxDi/Jxi7wFOSyk9wD8C/2crx9wMf8iPRVgYahlatd5sMhuzrDSnhrvZjUWkh4fyPRT0dSZhYrhl2Fje07JnTY7xgaYBY/BORxeRbG902LyfvLepF4fFsWbfWjXLfSEfmkszuu92QveAHm6/r6kPj8uz6jfpkQE8To9xvYAaGK812S0+vRUPO+NFplYU+vD/PeAR4BYp5TzQBrx/i8e+AghIKceklDHg68Abc8q8kXS3GMADwA2igqOrgVCAfns/DaaGNduynfasZqshLtmzpHLLDjYPYrfYz+ybp8tCF6XssAjZs7X0fbtsXcY03vUwCdOqWTz5RKlaxJIxxhfGVzlS7oQBSX1mnEmY0FwaxyLHWIwvAjC9NE0kEUFzaTQ1NDHQlJ7OpoSjtvQ09tDc0AxkXqicG4cGUmxOofk4IsAo8DohxPuATinlVv32+4DseMdTmXV5y0gpE0AIaN/icfMSj8R56lN/hvOJ30+viMXo+od/wHIinfgwWziANaE/stEflrll1nuA5IZF8If8dDm6cFldhgNiobmqs0Ol5BOlbJKpJHe+cCfBha3PaU/JFJ888EmOhNMOlOML4yRkwhh38bq8BMIBUjK15WPVEn/Yv+YayE3jq6/3uDwIxKrWh6L6CCHQnBpNliZ6G3vpb+rHbrave18sxBb4u+f+zohyPRYe45MHPmlcu0/MPMEPjv+gavbXI5v5cQAghPivwJ8A38ms+qoQ4h4p5We2cOx8LYfcdHiFlEkXFOIO0j4ndHZ2EgwW9zBMpBLEg9dA1EIwGGTXiy9y3le+wvHGRoJvehPnifO4atdVtC61EgwGubLxSkSnYP7YPPPMr6nvDV1v4FWNryIYDOJMOblm1zWcbzo/r10dqQ5OLJ/g4OhBWiwtvDT7Ev3WfqPsG7veyDmOcwr6TZfbLifcEWblxApBEaTX0svzc8/n3fdI5Aj3Hb6PlcUVbu+/vajzpROLxQgGg0wtT/HFQ19kIbzAuwbexa9P/hqApqUmgsEgbYk2lhPLPOV7qipZ8HS7ysl8fJ65lTk6kh0Eg0Eao2lfnyePPEnbUhtPH30aAGvISnApyGuaXoOz28mJo2eyLlfCrnKw0+26znUd59vPZ2IiHUlhwD7AgZkDeev++dzP+dfAvzLIINe2X8v9E/fzwLEHuNJ2JX32Pu5+5W4OLR7idzt/t+6ml1fr/1iQcJAea3i1lHIJQAjxf4BfA1sRjilWu6b1A7np8PQyU5nAii7gVL7KpJT3kM5SyN69e6XbXXxu5/3WKU4e9eJ2m9n1q18B0D03B243btxcfc7VRlm3281beeu6dX3M/bFVy1/Y84V1y15uvZx7J+9lxbnCee3nMfn0JL/T/zvov+Gj7o8SDAYp5De5cfMG3mAsX7x8MT85+RNc3a41U0IPTqTjR80yW1Dd+dDtOjx1GIATnMDtdnM6fBqLsHC192qsZitXNl/Jp498mkhzBHdf5fNuF3q+imHmxAwAVwxdgbvbzYAcoPFgI6fMp3C73cwem6WnsYfzR9JDdW63m3fwjorbVQ52ul1/7P7jVct7j+/liZkn8tY9H06/CM43zON2uzkxkRb+SFMEd7+byRcmWU4t07C7wQj5Uy9U6/9Y6BiH4MzMKjLftyq1TwOaEGKPEMIKvB3IjRr3feCdme9vBX4qi03SXQTn9J3mcHSEZCyJzZfuOrL5K98Pmt3lMbE4QSwVM7p4tspGOUJyu1i2gh6GIzsky1DLkDEZQO++284D5LmRAvRxjuz4YCom1fZAc2mcXDnJ6ejpNdv0+yH7WtaXT0dPG/5IZ/PgeqHC8SXgN0KIjwghPgI8Cdy7lQNnxizeR3rQ/RXgm1LKl4QQ/1sIob8y3wu0CyECpMOfrJmyW04858ZZwcHMU3OGYNj8fkhVtl++y9FFS0NLOoyJHkG3THP/c/vhs9HXZUcPLRVdmI4uHWUpvrQmpENTQxN9TX3b+mbzh/y0WlvpsHcY6/SglvFUnLGFMTUQvk3Y6L4wBCPsXxNmKPv63c4vQVul0MHxTwLvJt1NdBp4t5TyU1s9uJTyISmlV0o5IqX8aGbdh6SU3898X5FSvk1K6ZFSXiGlrKhDwvBl6amsR369gC0QINnYiDkSoeHYsUoe1gjHrse/EgiGncNlqTtblHLRMwrq0UO3gl4XwIFTB5hamsrrSLmdnQCzE2np6EEtnzv5HIlUIu8sO0X9sV7SNT0VdKOlkanFKQ6cOgBAo6URf/hMmCGHqbbhfGrNZtFx7UKIvxBC3AVcDnxOSvlpKeVvq2NedRn8nXYEKcaeWsIcDrN4ww0ARrdVJckWjsHmQcMfY6voopT7dqSHzrihL/0bt3ITrCRWCC4Gjbr00Cz5vOS3a0jrbG/+bPTlhyYeAtaPXqyoL3bbd+O0OtfcF3qA0Bv6bkiHGZpIX8s39N1AcCHIS6dfwml1srdl77Z+Cdoqm7U4vgxcBrxI2lHvzopbVEPsbXaGzUH8h9PRTsO3pkNlVGOcQ3NqhONhnjzxZNm7O/QQC9nDQ3rX0g19N9BgathSs1sPV/3antfiMDt4ZOoR47i5dmzXkNbTkWmWEkvrCscjU49gFmb2tOyphXmKIhFCpKeI5wkZA3DbwG0APHo0HWZof89+kjJphBkaahzaMHvnTmcz4ThfSvn/SCn/mfTg9P4q2FRTvM4JXkqkQ4pELr6YeE9PdYQj8wAKx8LlF46MKJ1YPjMtVL9hzm09l2Hn8JZaHHpd3lYvHpeHcCyMw+ygv6l/tR2b5GSvZ/TfmPu/abe302ZrIxwL4252rwqxrqhv9K7T7Bcqf8iPzWzjqq6rsJlthGNhPE6P0ZLU788hx1BdJEqrFZsJhyGnmcHsHc9w7xx+NBZ2D5BqbWXF662qcOR+L2fduU6G+sM935tXMfhCPhpMDbib3WdyoWe8q7PZ07IHi7Bsy0FFI45ZnkkLmzl4KuoTzaWtGvyG9AvCiHMEq9lqOG56Xd50mCGTxVgeahwCtudLUDnYTDguysk1ruceXxBChKthYLUZ1KIksfBiJmx2VNOwjY1BvLJN0lZbqxGevVLCkZscygidkYkeGo6V9i/Vw1VbTBbjwZovrpbVbMXdkj9xVb3jD/npbuymxdqyZpsSju1JdtI1HX/Yb1zDxmcmzJDeDelxeRhwDKybvfNsYLOc4+acXOOWrO/OahlZTXouNgPwQkva2S/q8SASCazj4wBYAwE677zTmKLb9POfs+trXyvLsb0ur/HmXk50Ucr25fCFfGuy0+k3wdGlo3z0uY8a/bdPn3iaL7yyvgPjqpstJxxHLrmtmx8Ef8D3g7nuO/VH9m/MJTvMi2L7kBvscD46z4nlE8b67NYzrA4hZDVZ183eWQgH5g5w18G7jOWHJx/mO0e+s8Ee9cWWItzuRNpe1YzVFOf5tvRwTtSbvlhsgfTF1fq979F+3300ZEIXtH/5y+z+1KegDH6Jbx1+K+8+591Gk7icZM+sOrVyirmVuTVvyvpN8MOJH/K1wNd4+fTLAHwt8DU+ffDTeX09FhOLzERmjDouar+Im/tv5vq+/FmFNZfG5NIkS/ElAD7/8uf555f/uYy/tPzEU3HGwmPrCsP+nv3c2Hcjl3deXmXLFFvBZXXR5egy7gv9xUr/P9/YfyM399/MRe0XAfCmoTfxB54/wGV1AWdmQpbCN8e+yede/hzz0bSX+j2v3MPnXvrcln5PNVHCkYPZamKPluTQfHpgNzY8jDSbsWc7BAL2jJDY/H7MS0tYpnOjpRTP6wZex1/s+4st15MPzaUxGh4lmUquuUF6GntosjQZ6/WbITtf+Hq+HsHl4Kq6mhqa+NTVn2KgOX+ic/3tbTQ8SjQZJbgYNHJ31ysTixPEU/F1u6I6HZ380zX/ZDxQFNuHXM9/fR2ko1t/6upPGf5J13Rfw19d8ler9p1YnCgpUVr2MROpBGPhMaYj00ak5XpHCUceNC2K35+eHSOtVmJu92pP8syn+dQpLHNzwBkhqVc0l0Y0GWVyaXJNwqpcXw9dMPwhv/Fwz16fTTCS3lao41u2x+6RhSMkZTItSjVKNlUI+k2unPt2HppLYyw8RiKVwBfy0dLQQpejq+B9S0mUpudtgfS1Nbk4SSyV9m0aDY8W9wNqhBKOPGhalGPHGlhYSJ+eqKZh8/kwhcM0zKQD3dl8vlWzrarhJLgVsruj/CE/LqtrdeiMzJtXPBVndCF98frDfsbCYyRl0tg3l/HlcSNcdSHoIa39If+aWV71ij/kxyRMKjz6DkRzacRSMSYWJ/JGBtgI4yWoSEfAo0tHjVZKtjc6nIn5Vu8o4ciDpqW7TQKBdKsjqmk0TE3hePFFAJKNjdj8fkM49OV6ZsQ5gkAY8Xa8Lu+a0BmhWIhnZp8hkUqkQyxkxc5qtDTmnUY7Hhkv6mbTk035w2k7LCYLDaaGuhcO5aOxM8lO3RwIBYqaGacnSiv2Ya9f68Y9lgkzZDfbt403uhKOPGhaehBY766Ker0IKWl5JO0RvXDjjViDQewvv0zC5WL50kvrXjgaLY30N/fjC/lWJSPS0W+gH038CIAb+27k5MpJnj7xNBaThf09+9dc1FJKgsvBoqehZodXGW4ZZrhluK5vmHyhRhQ7g+GWYQSCX838inA8XNTMuM0Spa2HLhzX916/KszQVgbbq40Sjjz09iZobEzi86WFY8WT7tt2PvooyaYmlvbvRySTtPz0p0S9Xla8XqxV8PXYKppT48kTT7IYX1zjZ6EvP3r0UczCzM39NwPw2NRjjLSMcF7recxEZoysaAAnV04SToQ3zYW+xo5MSOvfnvytkR63Xp0ClxPLTCxOKOHYodgtdgabB3n0aDqhaakvQcXgD/npa+rjVR2vYiG+wG9O/Aavy4vX5VXCsZ0RIt1dpXdVxQcGSNntmMNhoprGipa+uPTlqKZhisexTtR3+AHNpRlOfrn5PtrsbbTb243QGRe0XQBAOB5O52nO5yyVuciLzR1ihFeJhw3hyBWlekEPeqeEY+eSfV+U8hJ0fPk4oVio4H10n6Ds+0C/x05F01Pl6x0lHOvg9Ubx+Wxp9wyzmehIemA0qmnEhoaQFouxHM0ISb13V2U3w/PdINm5snfbd6+ar57X+1yfvlhk7pBsO/Q3LcifbKrW5OYRV+w89Gu709G5JktmofsW2lKIJWMcCR/B4/KsugfXu8fqFSUc6+DxRJmft3DyZNqTXBeHqKZBQwPR4WFjOTY8jDSZ6l449Omk3Y5unNa1jv/Z3t/6FF19ubex1xjM0/GFfLRaWmmztxVlR4e9Y5Uo6XbVY3eVb96H1WRloCm/X4pi+2Nc5yUkTzMybBYY6y24GCQhE3hdXlptrXQ6Oo16lHDsAPSZVX6/HUiHHoEcAcmslzYbMbcbe9aU3I6776bxqaeqafKmDLUMYTFZ1vVHyA0Xot9I+gys3P7cQDhgBHsrBr2uRksjPY099Db20mRpqssbJhBOB70zm8y1NkVRIXJDixRDt6ObloaWgl96cp0MNadmhBlqt7Wzy7arLu+DXMof22KH4PXqU3KtXH31Egs33YT9lVdYvvBCAObf/GYSHR2knOk396imYTt8GACxvEzHXXcRvuUWIldcUZsfkIcGUwN/cu6fcN6u8/Ju39+zn1sGbuGKzrTNbxh6AwmZoKexB0gLyGNTjyGlRCIJhALcvPvmkmx5h+cdBBeCRgRdPWdIveEP+bmy88pam6GoIEPNQ7xlz1t4/eDri9433wvVRvhCvlV5W35/5Pe5sP3CVZF367HlnYsSjnVoa0vS3p7A50u3OOIDA0x/4hPG9sjVVxO5+mpjOapptDz2GGJ5GdvoKELKVS2QeuHPL/jzdbftduzmk1d90li+qP0iI04PpMdFvhX7FidXTrKcXGY5ucyQY6gkO24duHXVsubSDFEq1Cek0oRiIY4vH1ce4zscs8nM317+tyXv73F6eHjq4YKuXX/Iz1DLEFazFUjHw7qx/8ZVdX13/LukZGpNWoJ6oiaWCSE+IYQ4JIQ4IIT4rhAi74iUEGJcCPGiEOJ5IcQz1bbT4zkTemQzopqGkBLb2Jgx1mENBhGx7ZcmdT2yZ1bpb1ildFXlw+P0MB+b5+TKybLUVw5yuxUUinx4XV7CsdWJ0tZjMydDzaURSURW5QipR2olaY8BF0gpLwR8wAc3KHudlPJiKeVl1THtDF5vekpuJoL6hhhRdH0+I/xIdjj2nYARYmH+jHAMOgbLU3frmRhW9YIx3VjNqFJsQL6p6vlYii8xuTS54SB8tid7PVMT4ZBSPpqVUfBJoH+j8rVC06IsL5s4erRh07KxgQFSVms6FEkgQLKpCaj/GFbF0GZvo93WbrQ4+pr6aDQ3lqVu/WaqN+EoJuid4uyk0NlQegBDrXV94dC7RevpPshHPXSi/RHwo3W2SeBRIcSzQog7qmgTsDb0yIaYzcRGRrAFAtj8fhZf+1qkxVL3U3SLxePyGGESyvkmni1K9UKxQe8UZyd6orTNHvaF+D01NzTT09hT98JRscFxIcSPge48m/5SSvlgpsxfAglgvRR610gpp4UQncBjQohDUspfrHO8O4A7ADo7OwkGgyXZHYvFjH2tVjOwh6efXmZkZPP6nL29dD75JA1LSxzr7sbU30/qwIGSbVnPrlrSberm4fmHScgElzZfWla7BmwDHDxxsC7Ol5SSw6cPs799f1nPe738H3NRdhVHrl391n5emn1pQ1ufnXoWm8lGci5J8NT65QasA5vWVahdlaJiwiGlvHGj7UKIdwKvB26QMn/6PCnldObzhBDiu8AVQF7hkFLeA9wDsHfvXul2l5Z+NRgMkr1vX1+M2dku3O7kpvtaLr6Yhp/8BIDGK64gdewYrS+8QKm2bGRXrbg0eSkPzjyY/j5wKVasZbNr36l9PDD2AAODA1ueUbLV83U8cpzF5CKX9F1S1vNeL//HXJRdxZFr14WnL+Qbo9+gf6B/XZ+f4+PH0Vwaw0PDG9Z9YehC7vfdT+9ALw2mzbvJN7KrUtRqVtUtwAeAN0gpI+uUaRJCtOjfgZuBg9WzMk12UqfN0J0C9e9Rjwfr0aOYlpYqZV7VyZ4RUu7ZRppLYzm5zNGlo2WttxT0LjM1o0pRCJpLYyW5wuTS5Lpl/CF/QVO7PS4PiVSC4EL9tbR0ajXGcRfQQrr76XkhxOcBhBC9QoiHMmW6gCeEEC8ATwE/lFI+XG1DNS3KkSM2CplVq8+sSjqdJLq6jGVrnWcHLAY9vo5FWAwnpnKhP6TrYUaJnmNBCYeiEDYLPXI6epqTKycLup6ys2TWK7WaVeWRUg5kptleLKX808z6aSnlbZnvY1LKizJ/e6WUH62FrZoWJZEQjI9v3upIdHWRbGlJhycRwmiB1KMjYKk0NTTR19SHu8VtODGVC12USr1hHhx/kG+MfiPvtlgyxoee/hDjC+MF1RUIB9ht31100DvF2Ul2ojRIX29//fRfM7GYjphdzNTuPS17MAuzsU88FefDz3y44Gu3GtTDrKq6Jjcb4IYIwcn3vpdT//k/AxDv6yPlcGDbQS0OgPec8x7e5X1X2evVRanU0CNfOvwlvnjoi3m3HZo/xANHHuDhycIarSp5k6IYshOlAbx0+iW+feTbPDKZTv5WjDOpzWzD3ew2uksPzx/mW2PfMpKs1QMq5Mgm7NkTxWKR+Hw2brtt8/Kn/uiPziyYTEQ9nh3lywHwds/bK1Z3qbF64qk4YwtjJFIJluJrx5T0G7eQ1kwylSQQDvD2kcr9TsXOQ3NqxsNev4b1680X8uG0Otlt311YXS6Nl+dfXlVXPXTh6qgWxyZYreB2xwoeIM8lqmk7rsVRSTxOD+ML48SSxYVqGV8YJ5FK+5Tmy+tRjHBMLk0STUZVi0NRFJpLI7gQJJaMrbneAqEAmrNwnyDNpTG1OEUkETlTVx35OCnhKIDsbIDFEtU0LHNzmOfqP6tXPaC5NBIyQXCxuBkl2d1b+cRBv+kKESW9LiUcimLwurwkZZKxhTHjGhpbGCOeiuMPF+cwq7k0JJKx8Jm6dFGqB5RwFIDXu8LkpJWlpeI9iFe2SXbAeqHUWD16uGqH2ZF3X1/IR3NDMwmZ2HSQ0RfyIRCMOEeKskFxdpOdkEy/3uKpOL858RsW44tFvYhkzzDU69JFqR5QwlEAHk96gHx0tPhWx3ZJK1svDLUMYRGWomdW+UN+3C1uIyRKNqdW0nmcb+xL+6Ru1uT3h/z0N/fTaClPHC7F2YGeKO3J408yH5s3rreHJtIeBsWE5x9oGsButvPM7DPMrsyeuXbrZIquEo4C0JM6lTLOkezoILFrlxrnKBCr2cpQy1DxwhH2ozm1vAmh9DGPm/tvLkiU9LoUimJoMDUw3DLMj4/+GIDX9b8OkzAZy8W0OMwmMyPOEWPfm/pvwmIq/oWqUijhKID+/jh2e6q0AfKMP8dOm1lVSYrJqAYQSUSYWpzC6/LidXmZi84xH583tutdV3t37d1UlGLJGMGFoAqlrigJzaWxGF8E4IK2C3A3u1mML9Ll6MJldZVc1/mt5zPcMqyEYzthMhWX1CmXqMeT7qrKH5JLkYPm0phamso7rTYfo+FRJBLNpRndAeORcWO7P+THZXXRYe/YVJTGFsZIyqQaGFeUhO7E2m5rp93efia3eAnXk16Xs8FJp6Oz6BeqSqKEo0DSMavsJe0b9XoxRyJYpqfLbNXORL/J9PwFm6HfTB6Xx+hiCi6fmZWlZ13T80NvJEoq659iK+gJyfQXmK0Ih16Xfu16nB6mI9NGK6SWKOEoEE2LcvKkhdOn05Evn3yykTvv7DS2P/BAK1//ev7wFEboETVAXhC5iXHmo/O8/8n3c2rlVN7y/pAfu9nOQNMAHfYOdtl2GS0OKeWqqZBGTKE8vh56XRaTBXdL/UVkVdQ/+otLrmCUIhy5delCUmpkhXKihKNA9NAjenfVt7/dyn33tbOwkD6F99/fxpe+1J5336gn/fahZlYVRn9TPw6zwxCOX5/4NT+c+CFPzDyRt7w/5GfEOYLZZE63Kpwa48vjAByLHGMxvmg0+zcLRucP+RluGS46nLVCAdDX1Mft2u280f1GAK7qvIo3Db2J/d37i66r09HJO73v5E1DbwLqK0umEo4CyZ1ZpTsEBgI2olFBMGhlaqqBSGStr0eqpYV4T48SjgIxCVN6Wm34jNctrH/D5MaV0lwawUjQaG3Ambe1XFHarC6FohiEEHzwVR9kb9teAFqsLfz9FX9Pm72tpLo+cPEHuLD9QgB6m3pxWBx14UGuhKNAOjoSuFwJfD4b8fgZnw6fz8bYmJVkUiClWNfXI6ppSjiKwOP0rInRk8+xbz46z+zK7BrhWE4tMx2Zxj/vN+qDM6KUr67F+CLTkWklHIq6xCRM6XhYqsWxfRDiTOiRYNBKIpFuWfj9tlWzrdabebWiaVjHxiAer4q92x3NpTG3MseplVMbxurJN5idPUbiD/vpdnTjtDqN7R6nJ29dKtSIot6pl5lVSjiKwOtNT8n1+dKzq1pakoZwWCxyQ1+PqKZhisexTkxU0+Rtiz6Y/eKpF5lcnKSloYWZyAwLsYVV5fSWQ7bfRW64hlwh8Lq8hihtVpdCUU94XV5ORdOREGqJEo4i8HiiLC6a+eUvmzCbJdddt2AIychIlJGR9X09jNAjyhGwIPSH/cOTDyOR3NR/E7C21eEP+deEq25uaKbT2smh+UOMhcfWCIfRIslTV6OlkZ7GnrL/HoWiHGTHw6olSjiKQJ9Z9eMft+B2x9i7d4X5eQvPPefA44lumJ88NjyMNJtV6JEC6bB34LK6jJALtw2mk6HkNtMD4fzhqt2Nbn557JfEU/H1hSNPXR6nB5NQt4WiPtlsVmC1qMkdIoT4iBDiaCbf+PNCiLwpkoQQtwghDgshAkKI/1VtO3PRhWNpyYymRY2ZVktLZrzetHDMzjYwP29es6+02Yi53WqAvECEEHhdXpYSS9jMNq7YfQVNlqZVD3spJf5Q/nDVQ44hlhJpJ79c4eiwd9BqbV1Tly/kU91Uirqmw95Bm63trG5x/GNWzvGHcjcKIczAZ4FbgfOBdwghzq+2kdk4nSm6u9OD25oWNYQkd3nd7iqPZ0flH680+gN/xDmCxWRZMzA4szzDQnwh72D2UOMQkJ6JMtwyvGqb7kGeXddcdI7T0dNqYFxR99TDAHk9t8mvAAJSyjEpZQz4OvDGGttkhFjXtChtbUna2xPG8qbCoWk0TE4ilperY+w2x3Day/Kg9Yf8yEzMr+xQI7m4HWnP78HmQeyWtaFiiqlLoagnPE4PgXCAlEzVzIZaCsf7hBAHhBD3CSF25dneB0xmLU9l1tUUvXtK01Yyn1EaG5P09sbp7EzgdCbx+dYRDq8XISW20cJiMJ3tZMfq0T/nY/OcXDkJbBxXasAxgFmY1+160lwaS4klpiPTq+pSXVWKesfr8hJJRJheql3sO0ulKhZC/BjozrPpL4G7gb8FZObz/wJ/lFtFnn3XDS8rhLgDuAOgs7OTYLC41KM6sVhsw1iCyv8AAA7zSURBVH0vv3yO06f7kTJAMAi33BLhwgsbmZhI7zM42MXBgyJvHY2NjfQDi08+ydGWlrLaVSsqaVdLqoXbOm/jAtMFBINBnMtpX4wn/E9wiesSnp9+nvaGduaPzTPP/OqdE3B73+2c03xOXvtaltPn/1f+X3HFriv47dHf4rK4WDy+yCKVCyJ3Nv4ft4Kyay36tfsfgf/gyl1XrtpWLbsqJhxSyhsLKSeE+ALwgzybpoCBrOV+YF2JlVLeA9wDsHfvXul2lxakLhgMstG+bjdce+0y4DaWM1sA2LcPfvADJ4ODbtbkpe/vJ2Wz0XvqFJYi7dvMrlpRabvu3HOn8b1lpQUOQdgWxu12M314mvPaz8t7/GAwyP+8+n+uW++u2C54GUK2EG63m2P+Y5zbdm7Fz/HZ+n8sFWXXWtrj7auu3VrYVatZVdkT5d8MHMxT7GlAE0LsEUJYgbcD36+GfVtB06IsLJg5fjyPJpvNREdG1JTcEmmzt9Fub8cX8pFIJRgNj5Y8mO20Oulu7MYX8pGSqXTWPzUwrtgGNDc009vYW9MB8lqNcXxcCPGiEOIAcB3w3wCEEL1CiIcApJQJ4H3AI8ArwDellC/VyN6C0QfP1x3nUNkAt4Tm1AiEA0wuThJLxbb0sNecGv6wn+mlaZYTyypdrGLbUOuZVTURDinl7VLKfVLKC6WUb5BSHsusn5ZS3pZV7iEppVdKOSKl/GgtbC2WgmZWnTiBaX4+73bFxug5xQ+HDqeXt/Cw11waY+ExDs0fMpYViu2A5tIYWxgjnqpN7Lt6no67LXG5UnR1xTcPPaK6q0rC6/KynFzm8enHEQhGnCNbqiueihve6WoqrmK7oLk0EqkE4wvjNTm+Eo4KsFGaWZUNcGvoD/cfH/3xuj4apdTV29hLc0NzWWxUKCpNrUOPKOGoAJoWZXTUSjK5dluiq4uk06lCj5SI7hQYSUS23LU03DKMSZjKUpdCUU2GW4YxC3PNxjmUcFQATYsSi5mYmLCu3SgEUY9HDZCXSFNDE/1N/cDWxyTsFjvuZndZ6lIoqonVbMXd4q5ZzColHBWgkAFyWyAAcl1/RsUGZHuS11NdCkU10WcF1gIlHBVgZCSKEHLD0CPmcBjL8eMANP7613T/zd9U08Rtjf6QL0d4kHLWpVBUE6/Ly9TiFJFEpOrHrpjn+NmM3S4ZHIwRCKyfRhbA5veT6O7G9f3v0/rgg8y+730k29uraeq25I3uNxJPxRlqGdpyXa8ffD2L8UVj7ESh2C5oLg2JZDQ8yr62fVU9tmpxVIiNkjpFPemHlD5Anvup2Jg9zj28/6L3lyXhkrvFzQcu/gBm09ocKgpFPaPPCqzFALkSjgqhaVGCQSvR6NpYjanWVuKdnelxjmTSiJarhEOhUBTKQNMAdrNdCcdOQtOipFKC0dE8M6s4E3rEOjmJKZoeTFfCoVAoCsVsMjPiHFHCsZPQZ1atN84R9XiwjY5iO5QOd6F8OxQKRbHUKmaVEo4K4XbHaGhIrT/O4fViikZp+elPkUKwcOONaeFQU3QVCkWBaC6N2ZVZ5qPVjX2nhKNCWCwwPBzD59s49EjLT35CfGCA5QsvxByJ0DBdu6xeCoVie6FPJ692q0MJRwXZcGbVyAhSCEzLy0Q17UzwQ9VdpVAoCkSPDq2EYwfh9UaZmWlgYWHtaZYOB/GBdILDFa9XCYdCoSiaTkcnTquz6qFHlHBUkEJCj0B6oDzV3Ey8p0fFsFIoFAUjhDCSm1UTJRwVxONZATbOBpj7qVocCoWiGPSZVbKKE2t2ZMgRIQSJRKKkfaWUJe+bS2dngqamJD6fNW+dc69/PalYjMjgICQSRDweOn71KxLLy9DQUDG7yomyqziUXcWh7NqckZYRFuILHF04WjW7dqRw2O129u0rLXbL3Nxcyfvm48IL4dixdvbtyxODat8++N3fpVNfvu46uO8+9tntcP75FbWrXCi7ikPZVRzKrs2Zd87z0ec/SqItgWPZURW7atJVJYT4hhDi+czfuBDi+XXKjQshXsyUe6badpaDCy6AgwcLdM+44IL058GDFbVJoVDsHC7oTD83Dp6o3nOjJsIhpfx9KeXFUsqLgW8D39mg+HWZspdVybyysm8fzM3BzEwBhc89F8xmePHFitulUCh2Brscu+hr6ePFE9V7btR0cFwIIYDfA/6tlnZUkqIaEXY7aJpqcSgUiqK4oPOCqrY4aj3G8TvAcSnlelOJJPCoEEIC/yylvGe9ioQQdwB3AHR1dfH444+XZNDi4mLJ++YjFGoAruF73wvQ0DC1afnzu7pofvppnsqxodx2lQtlV3Eou4pD2VUYrdFWfnr8p4R6Q9WxS0pZkT/gx8DBPH9vzCpzN/A/NqijN/PZCbwA7C/k2JdeeqkslZ/97Gcl77senZ1SvvvdBRb+yEekFELKxcWK21UOlF3FoewqDmVXYdz/2/slH0F++YdfLrkO4BlZ4PO9Yi0OKeWNG20XQliA/wRcukEd05nPE0KI7wJXAL8op53VYN++Inqf9u1Lj6S/8gpcti2HdRQKRZXRB8iPLB2pyvFqOcZxI3BISpm3/0YI0SSEaNG/AzeTbrFsOy64AF56CVKpAguDGudQKBQFc/7u8xEIxpfGq3K8WgrH28kZFBdC9AohHsosdgFPCCFeAJ4CfiilfLjKNpaFCy6ASASOFPIyMDICNpuaWaVQKArG0eDA0+ZhbGmsKser2eC4lPJdedZNA7dlvo8BF1XZrIqQ3YgYGdmksNmcdv5TLQ6FQlEEF3RewDPB6ri71XpW1VnB3r3w+teD01ngDm95CywtVdQmhUKxs7jVcyv/f3t3GyNXWYZx/H/ZdatgSamlslpiW9KaNDFABdOCGkQUMKTG6IeaRmtUEk1E1PjSZhOj3wSJookRE9QYrWgtFUsTA4rEj20B7QsvW0qosgi2/QCCGAF7++G5Nx3W2cCjnfNM2+uXTGbOc86mV+/ZM/fMc87Oee7J54gIyl86DI4bRwfmzIHbbqv4gfHxgWUxsxPTVW+5iqVPLx140wB/O66ZmVVy4zAzsypuHGZmVsWNw8zMqrhxmJlZFTcOMzOr4sZhZmZV3DjMzKyK4mVd0/T4IukQ8Of/8cfnA4ePYZxjxbnqOFcd56pzIuZ6Y0Sc8XI2PCEbx/9D0t0xhJepda46zlXHueqc7Lk8VWVmZlXcOMzMrIobx3+b8brmjTlXHeeq41x1TupcPsZhZmZV/InDzMyquHEkSZdLmpC0X9L6hjnOknSXpAck3SfpmhyfJ+m3kh7K+9Mb5Zsl6Y+StuXyYknbM9cvJI02yDRX0mZJD2bdVg1DvSR9Lp/DvZJulvSqVvWS9ENJByXt7RnrWyMV38l9YbekFR3n+kY+l7sl/UrS3J51GzLXhKTLuszVs+4LkkLS/FxuWq8cvzprcp+k63rGB1OviDjpb8As4GFgCTAK7AKWN8oyBqzIx3OAfcBy4DpgfY6vB65tlO/zwM+Abbm8CViTj28EPtUg04+BT+TjUWBu63oBbwAeAV7dU6ePtqoX8A5gBbC3Z6xvjSiXb/4NIGAlsL3jXO8BRvLxtT25lue+ORtYnPvsrK5y5fhZwO2UvxObPyT1eifwO2B2Li8YdL0G/gt7PNyAVcDtPcsbgA2tc2WWXwPvBiaAsRwbAyYaZFkI3AlcAmzLHeVwz07+ojp2lOm0fIHWtPGm9crG8Sgwj3KlzW3AZS3rBSya9oLTt0bA94EP9duui1zT1r0f2JiPX7Rf5gv4qi5zAZuBc4ADPY2jab0ob0Yu7bPdwOrlqapiaiefMpljTUlaBJwHbAdeFxGPA+T9ggaRbgC+BBzJ5dcCT0bEC7ncom5LgEPAj3IK7SZJp9K4XhHxGHA98BfgceAp4B7a16vXTDUapv3hY5R389A4l6TVwGMRsWvaqtb1Wga8PadA/yDpgkHncuMo+l2kt+npZpJeA9wCfDYi/t4yS+a5EjgYEff0DvfZtOu6jVA+un8vIs4D/kGZdmkqjxe8jzJF8HrgVOCKPpsO42mNw/C8ImkceAHYODXUZ7NOckk6BRgHvtJvdZ+xLus1ApxOmSb7IrBJ5cLjA8vlxlFMUuYupywE/tooC5JeSWkaGyNiSw7/TdJYrh8DDnYc6yJgtaQDwM8p01U3AHMljeQ2Leo2CUxGxPZc3kxpJK3rdSnwSEQciojngS3AhbSvV6+ZatR8f5C0DrgSWBs5z9I419mUNwG7ch9YCNwr6czGuch/f0sUOygzAvMHmcuNo9gJLM0zXkaBNcDWFkHyncIPgAci4ps9q7YC6/LxOsqxj85ExIaIWBgRiyj1+X1ErAXuAj7YMNcTwKOS3pRD7wLup3G9KFNUKyWdks/pVK6m9ZpmphptBT6SZwutBJ6amtLqgqTLgS8DqyPi2Wl510iaLWkxsBTY0UWmiNgTEQsiYlHuA5OUk1ieoHG9gFspb+SQtIxygshhBlmvQR3AOd5ulDMj9lHOPBhvmONtlI+Tu4E/5e29lOMJdwIP5f28hhkv5uhZVUvyl3E/8EvyzI6O85wL3J01u5Xysb15vYCvAQ8Ce4GfUM5uaVIv4GbKsZbnKS96H5+pRpQpju/mvrAHOL/jXPspc/NTv/839mw/nrkmgCu6zDVt/QGOHhxvXa9R4Kf5e3YvcMmg6+W/HDczsyqeqjIzsypuHGZmVsWNw8zMqrhxmJlZFTcOMzOrMvLSm5jZTCRNndIKcCbwb8pXoAA8GxEXNglmNkA+HdfsGJH0VeCZiLi+dRazQfJUldmASHom7y/OL5/bJGmfpK9LWitph6Q9ks7O7c6QdIuknXm7qO3/wKw/Nw6zbpwDXAO8GfgwsCwi3grcBFyd23wb+FZEXAB8INeZDR0f4zDrxs7I7y+S9DBwR47voVyIB8oXIy4vX20FwGmS5kTE050mNXsJbhxm3fhXz+MjPctHOLofvoJyoZ1/dhnMrJanqsyGxx3Ap6cWJJ3bMIvZjNw4zIbHZ4DzJe2WdD/wydaBzPrx6bhmZlbFnzjMzKyKG4eZmVVx4zAzsypuHGZmVsWNw8zMqrhxmJlZFTcOMzOr4sZhZmZV/gOrTj7OSB5FOAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()  # initiate plots\n",
    "\n",
    "a = -7  # lower boundary\n",
    "b = 10  # upper boundary\n",
    "\n",
    "for col in ['red', 'green', 'blue']:  # 3 random walks r-ed, g-reen, b-lue\n",
    "    walker = rand_walk_Z()  # start new symmetric walk from 0\n",
    "    x = 0  # start position 0\n",
    "    chain = [x]  # save positions to plot\n",
    "    while a <= x <= b:\n",
    "        x = next(walker)\n",
    "        chain.append(x)\n",
    "    ax.plot(chain, color=col)\n",
    "\n",
    "# To look nicer\n",
    "plt.axhspan(a, b, color = '#cccccc', alpha = 0.8)  # box between a & b\n",
    "ax.grid()  # adds grid lines\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Position')\n",
    "\n",
    "# Show image!\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Further Information\n",
    "[Wikipedia - Donsker's invariance principle](https://en.wikipedia.org/wiki/Donsker%27s_theorem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random walk on $\\mathbb{Z}^2$\n",
    "\n",
    "To implement the random walk on $\\mathbb{Z}^2$ a new generator function is necessary. In the cell below we give you `rand_walk_Z2`.\n",
    "\n",
    "The logic is the same as `rand_walk_Z` just this time:\n",
    "\n",
    "* 2 variables `x, y` are used to determine the position\n",
    "* `steps` is a list of tables representing the possible shifts\n",
    "* In each iteration one element of `steps` is picked and added to `x` and `y` respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# random walk in Z^2 generator:\n",
    "def rand_walk_Z2(start=(0, 0), probs=(0.25, 0.25, 0.25, 0.25)):\n",
    "    # probs order (left, right, bottom, up)\n",
    "    x, y = start\n",
    "    steps = [(-1, 0), (1, 0), (0, -1), (0, 1)]\n",
    "    while True:\n",
    "        dx, dy = steps[choice(4, p=probs)]  # choice(4) = choose a number in [0,3]\n",
    "        x, y = x + dx, y + dy\n",
    "        yield x, y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application: Escape Time from Disk\n",
    "\n",
    "For familiarization purposes, the simple symmetric walk on $\\mathbb{Z}^2$ is implimented until the disk of radius $r = 20$ is abandoned.\n",
    "\n",
    "The code cell below simulates this chain 3 times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAEKCAYAAAA8bsGsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnX2QZFd12H+nP6Znpud7Z2clrZB2V8KxJBYvxUaAjWAoOza7xsG4oirshJBA2bgKYyeFs7AocahQWETYuBJiOyiYAju2MTGSwfbKWKYYhIKwvRhFHwisZbQrtDvamd2ene/+Pvmj+/X09Lye/nrd73X3+VVNTfft+947t7vf6XPPPfccUVUMwzC8IuS3AIZh9BamVAzD8BRTKoZheIopFcMwPMWUimEYnmJKxTAMTzGlYhiGp5hSMQzDU0ypGIbhKRG/BfCC6elpPXTo0K72jY0N4vF45wVqA700FrDxBB238Xzzm9+8oqr7ax3bE0rl0KFDnD17dlf73Nwcs7OznReoDfTSWMDGE3TcxiMiF+o51rfpj4i8RES+IiLPiMjTIvIrxfYpEXlYRJ4t/p/0S0bDMBrHT59KFnivqt4GvBp4t4jcDrwf+LKqvhT4cvG5YRhdgm9KRVUXVPUfio/XgGeAg8Cbgc8Uu30G+Gl/JDQMoxkC4VMRkUPAK4C/BQ6o6gIUFI+IzPgomlGDfD5PPp8nl8uV/udyOTKZDNlsttTupNhw/qfTaS5cuICIACAihEIhwuEwkUik9Oe0hcNhQqEQIlI6xggm4nc+FREZAb4KfFhVHxCRa6o6Ufb6sqru8quIyC8AvwBw4MCBV372s5/dde719XVGRkbaJ3wHCcJYVLX0V64o9qKaAkgmkwwODrpeo55zOkooKEomCJ+Pl7iN5w1veMM3VfV4rWN9tVREJAp8HvhDVX2g2HxZRK4vWinXA4tux6rq/cD9AMePH1c3z3sveeQ7OZZ8Pk82myWdTpNKpUgmk2QyGaBw04dCodIfVFcce/H0009zxx13NHyco9AqlVo0GiUWizE4OEg0GiUajZbk6wS99F2D1sbjm1KRwjfx94BnVPVjZS99EXg78JHi/y/4IF5foapkMhmSySTr6+uk0+lSu6M8otFoICwCESlNhxwcRbOxscH6+nqpPRqNEo/HGR4eDoz8/YCflsqPAG8DnhSRx4ttH6CgTD4nIu8Engfu9km+nkZVSaVSbG5usrGxQS6XQ1WJRCJddwPupWiuXbvG8vIy4XC4pGBisVhHrZh+wzeloqqPAtW+uT/aSVn6hVwuRyqVYmNjg83NzdL0IRKJMDAw4LN03lKpaPL5PGtra6ytrQEwPDxMPB5ncHBwhzIyWicQqz9G+3AsktXVVTY3N4HCDReJRLrKGmmVUChUUpyqytbWVun9GBwcZHx8nMHBwb56T9qFKZUeJZfLsbGxwcrKCtlsNlB+Eb8REaLRKLCtdF988UUikQhjY2OMjIyY9dICplR6CFUlnU6ztrbG+vp6yUcSi8X8Fi2wlCuYXC7H8vIyy8vLxONxRkdHicVipogbxJRKD5DP59nc3GR1dbW0cmNWSeM4PhhVZXNzk/X1daLRKOPj48TjcXPu1okplS5GVVlfX2d5eZl8Pm9THI9wrBdnBenq1askEgkmJycZGRkx5VIDUypdiPNLury8TCaT6cnVmyBQvoKUz+dJJBJcu3aNqakp4vG4Ke8qmFLpIlSVZDJJIpEgnU6bv6SDOKtHuVyOpaWlknIZGhoy5VKBKZUuQVW5fPkyW1tbhMNhBgYG7MvsA47lks1muXz5MrFYjKmpKb/FChSmVAJOJpMpTXNSqZQpk4AQiURKymVhYaG0V8qmoaZUAouqsra2RiKRAHYufRrBwAkidHwuly5dYmJigvHx8b5W/KZUAkgmk2FpaYlUKtXx3bZG4zjpF6LRKMvLy2xsbDA9Pd23/i77tgYIVWV1dZWLFy+SyWQYGBgwhdJFiAixWIxcLsfCwgLXrl2rKz9Mr2GWSkAw66R3iEQiqGrfWi32zfUZs056k362WsxS8ZFsNsvS0hLJZNKskx6l0mqZmZnpeYe7fYt9IpVKcenSpdIysSmU3sWxWrLZLJcuXSqlXOhV7JvcYZz9OgsLCwAWd9JHONbo5cuXWVlZ6dnpkK9KRUQ+JSKLIvJUWdsHReSiiDxe/Dvpp4xeoqokEgmWlpZKJSiM/iIcDhONRkkkEly5coV8Pu+3SJ7jt6XyaeCNLu2/parHin9nOixTW8jlciwuLrK6umrTnT7H2Ue0vr7Oiy++SDab9VskT/H1m62qjwAJP2XoBOl0moWFBba2tmy6YwAFP8vAwACZTKbkW+sVglBM7BDwF6r6suLzDwL/BlgFzlKot7zsclxXFBNzauhAc/VxHKoV3+pWbDzblCcgD4oF20oxsSAqlQPAFUCBDwHXq+o79jrH8ePH9ezZs7va/S7wtL6+XvKftJrztNniW0HFxrOTfD5PJpNhamqKsbEx361Zt3tHROpSKsFQi2Wo6mVVzalqHvhfwJ1+y9QMa2trLC0tEY1GLYmyURMna18ikej6laHAKZViqVOHtwBPVesbVFZXV7ly5YoFtBkN4Thwl5eXuzoC1+9ayn8MzALTIvIC8J+BWRE5RmH6cx54l28CNsHKygqJRMIUitEUjgP32rVrAExMTPg+FWoUX5WKqv6sS/PvdVwQj1hdXTWFYrRMtysW++Z7xNraGlevXjWFYnhCuWLpNh+Lffs9YH193Xwohuc4imV5eZnV1VW/xakbuwNaZHNzs7TKYwrF8Bono1wikWB9fd1vcerC7oIWSKfTLC4uBipoyeg9nOXmK1eudEXkrd0JTeLs5XEKThlGOwmFQqUdzkHfK2RKpQlUlaWlJbLZbM8n3DGCQyQSIZ/Ps7i4GOjdzaZUGsRJX7C1tWUKxeg4AwMDpNNprl69GtgVIVMqDbK+vl5KX9BNsQNG7xCNRllfX2dlZcVvUVwxpdIAyWSyFItiCsXwi/Kl5iCmpjSlUifZbJbFxcWSw8ww/MSpjri4uEg6nfZbnB3Y3VEHqlpyjlkKSCMohMNhRCRwjltTKnWwurpaKvJlGEEiGo2SyWRK+4SCgCmVGqTTaZaXl82PYgSWgYEBVlZWAhMYZ0plD1SVK1euICLmRzECixOAubS0FIhpkN0pe+BMe8yPYgSdSCQSmGmQKZUq2LTH6DaCMg0ypeKCTXuMbiQo06AgViicEpGHReTZ4v/JTstl0x6jWwnCNMjvn+FPs7tC4fuBL6vqS4EvF593jEwmY9Meo6vxexrkd47aR4p1f8p5M4Vk2ACfAeaA93VKpuXlQt0ym/YYXnPmzCBzc7Fd7bOzKU6eTDZ1ztOnxwG4997tfUDOtD2RSHDdddd1/McxiPb9AVVdAFDVBRGZcetUUaGQubm5XX3W19dd26uhqmQymUBaKMlkkqefftpvMTyjH8fzwAO3cfEiHDy4vV/n4sVhEolNbr75maau+9hjrwJwvbaq8swzzzT1A9novVNOEJVKXajq/cD9UKhQ6FaJsJEKharK5cuXAxs5axX9gk0945maGmdqCu67b9sqOXUqBsSafi+GhgqlVt2Oz2azhEIhDh482PAPZSvVPYNo4192CooV/y924qLJZJKtrS1zzho9g+O03djY6Oh1g6hUvgi8vfj47cAX2n1BJ/GSs0HLMHqFSCTC8vJyR5eYg1ih8CPA50TkncDzwN3tlmNjY4N0Ok0sttuJZhjN8J73THDunPvtderUeOnxk08WptonTkzv6nfsWGaHA3Yvys9ZSS6XIxbLcvToAO99b12nawm/V3/cKhQC/GgHZWB5edmmPYanVFMofhAKhUinU+TzEToxOQnOyH1ifX2dbDZrVorRFh566ErpsWON3Hefd2kgjx7N1HXOdDrNxEQOmPDs2tUIok+lY+Tz+VKgm2H0MpFIhGvXrpHL5dp+rb5WKpubm+RyOQt0M3oe5zveiSqHfXs3qSorKyvmSzH6hkgkwurqattLe/TtHZVOp0mn0wwMDPgtitHlnD49zmOPvaoUiFaLypWaVsL0GyEUCpFKpUgmkwwNDbXvOm07c8BZXV1FRCwuxfCV+fmI636gdhEKhdpeL6gvLZVcLsfGxoY5aA1PuPfelV1h+m5xJ24rNXvFl7SDSCTC1tYWmUymbd//vrRUnLBls1KMfsOxztvpsO07peI4aMPhsN+iGIYvOA7bdoXu9930J5lMWrBbH9OunCbVHLVuIfn1hunD9pSpGvPzEY4cyTYkbygUIpPJsLm5ycjISEPH1nV+z88YcFZXVy0upY+Zm4sxP7/zt7TTzlIvOXIky+xs4xnewuEwKysrbVle7itLJZfLsbW1ZQ7aPufIkaynzlI3R61zTrfrlLc5Fkp5OH95u5ch/eWEw2EymQzZbNbz+6GvfrKTyYJ5aw5ao98REVSVra0tz8/dV0plY2PDFIphFAmHw21J4NQ30598Ps/m5qZNfXoAN2er4/Asd2w6bW6UT3macXaWs5ej1s0BW2/brus8ehqAe197b6ntPV95D+eunXPt/9BbHtrzfOFwmFQqRS6X83Q1tG+UilOuwCyV7sdxtraiCMpp1tnZaR5fenxXWzWFUg/OvZBMJonH402fp5K+USqbm5u1OxldQ6Wz1c2xuZdj1EsH6F4RtZUO2Erq7VeLcqvkxIMn6j5ORNjY2OgPpSIi54E1IAdkVfV4s+dSVTY2NmxHsmFUEIlE2NzcJJ/PexZqEfS77A2q2poKp1B1MJfLmVIxjAqcVaBUKuXZzuW+uMs2NzfNl9JDuEWmOjQSrergVuWvkb4/95t/ynN/9wPcdfP2tZ9fe55UNsXRu3cWCds3NM2d1/3TmtfZC7fpjVvbqa+d2vF89sZZTh4+uaufiLC5udkXSkWBvxYRBT5RLB5WopEKhZnM3qHO3UA/VvSrxtZWoSpfIrFa1jZWtc2N8mvvVeWvEre+3/q/B0gvTZIYXd6+dqYQ/5GtSN+YSiVJJMr6FWWs59qNxpQklhOlxxe3LpJYTnDz5s2ufVV1R26hXq1Q+COqeqlY9vRhEfmOqj7ivFhvhcLXve51PP/8811fcL0fK/pVw1m6/eQny3faOo+3/QInThT6uSWfLr/2XlX+ql27vG8k/HXYf55PfvIl29d58F8Vrl1jWXdb7trXvjNxJwD33XVf2XVO7LqOY6GU93Paqo0xnU5z4403llwEvVahEABVvVT8vwg8CNzZzHkcK6WbFYphdAKvLPpAKhURiYvIqPMY+HHgqWbOlU6n256T0zB6ASeWq1WCOv05ADxYtC4iwB+p6l81c6JUKmW7kvscN4dus1Gtbn23lm8jr/m6HahHp4/ueF7NgVrJk1eeBHY7YN2YX5nf0W9+ZZ4j40eq9nfy13pBIJWKqs4DP+TFuZLJpCmVHuPYse53vDvMr8wD1KVU6mX2xtldbUfGj7i2O4TDYZLJJKrasqsgkErFS9qxtdvwl3rrC9dbvc8tqrVWSoLy9qN3P0M2l3ONaq3lqK3H6qjEzVFbzsnDJxtWUk68ihfxXD39E+74UsxJaxj14YWztueVijlpDaN+vPCr9LRS8XI/g2H0Ol45a3vap6KqplSqcOa5M8y9MNfSOepdtWgWt/whjfDo2U3WMmv87dvO72i/Pn4DN429ZFf/elaJVlIrrGXWOPW1Xytr/enC8U34R2D3So1Du9/fSpz8Kq1a9z2rVJypjykVd+ZemKu5zLgX7Vi1qMQtf0gjrGXWSGV3/vKupteAS65KpdlztkK1FRm39/fY/mOeXdcNESGfz5tSqYZT08SctNU5Mn5kx0pCIzT7q9xpYpEYX/2DHyw93w5hL2trICH10bsK59zxvv35OInEMvfddVfD8lVbqXF7f5u12BolV7FfqVF69me8XYWSDKPXafXe6Vml0qq2NYx+pdV7p+enP4Y7Tsh3vakHa++27Rwvf+vXAHjju/+81Pbo38yw9u0f3tFvK3MbALe+ujxlQMH5+9BLvr7d7/uFfpw65L2wLVDpwHU+Mzfc+tXaDuDmCFfVlu+dnlUqZqn0LpvnX1Z8tK1U1r79w6QWbyY2c8EfoTxmr5D6ZnBz/Lo5wkWEbLa1hOI9q1RafWN6HedXzC3nhltbECmX86H//nViMxd48v9sWyu3vrqQDOncNyZ3tZX3O3p3wWq5777tfn7j5sCtN3eK06+Zz1FEWo6q7VmfiikVw2gcEbHVn2pks1lbTjaMBgmFQjb9qYb5VGpTb84Nt4jPRh29lWxtbTF3x1zp+b989zUublziznd+elff3df4I2B72gLbztbyNqje9vq3fafU4vhi6pkipLL3EIvEdkXaFnLNrrsf1ASnHz1dNfivniTXlW17BTqW93OCRieyE8wy26DUBXpWqdjqz97Um3PDa4dhNZ54YpBU9mBHrlVJbOYCo7d/vXZHYP9NV9Glf9JmibynVj4VByeqthUCq1RE5I3AfwPCwCdV9SONHG+7k/em3pwb1frVmy/ErV8166YyUrXaNV7/QMHK+Oof7Ha2ujlq3drco2zriC6+CyADbEfanjo1Xszi7503wS16tl453Ry1brg566GwU/nCheZX0aq+CyJyRkQONX3mFhCRMPDbwAngduBnReR2P2QxjH6jVV/kXqr10xTq7twjIp1OnXYncE5V51U1DXwWeHMjJzBLxTCao20bClX1cyLyl8CvAWdF5A+AfNnrH2vpyntzEPh+2fMXgFe18XpGk9QTyQk7pzxXvvAz5NNpxh/fPlbW3wfRKCdO3Fhqe3b5jxg+9BS8Zee5EskEqcWbXR21zvSmnKN3JfeUx40zZwaZm4vt2cdhfj7CxERdXVumWpoEN2r1q+a8bXeO2gywAcSAUcqUSptxG9UO9VmrQmE6nSaVSvVMVb+gVSh0quWVV8Hbq23HsUuHIZ9nefhb2425HKrK1ta2AhiQAUZkdNe4Y7fOkc3dRbZsgS+v1b+a2bKOTj83ucqv88ADt3HxIhw8uFn1vA4TE/BDP7TA009fq9m3FW7RW0jkEzveXzfcPgc3JpjgFr1l1/urqqTT6aYrFEo1U6foKP0Y8EXgv6hq7XfXI0TkNcAHVfUnis9PA6iq697v48eP69mzZ3e0nT9/nmeffbZnqvp1Q4XCvZyy5W2vf9t3SKVTfONPtgsmvPwVVwF44lv7mrqOY408+bXBUptbRK1bP9dIYpd0CHvRDZ9PvaTTac6fP7+rQqGIfFNVj9c6fi9L5R7gblX14+fx74GXishh4CLwVuDnGjmBBb4Zhj/s5VNpPOOMR6hqVkR+CfgShSXlT/mk3Ayj7+jZzG+qegY40+zxZqn4R63ozrW1n0Iv3sQ9r9+OaiU/RSq806m6ldli+NBTO1IcADz/mQ+TevEwRz+23TeVTblEuhamP27O29eX5a1dOFeIkn1ypP4KhUZ1enbvjymV4HL7gUe4LfTcjrYxBolFBqscsZPUi4d3tcUiMUajo03J08qxvYZVKNyDcDjstwh9h1uEplt05/ifn2L56P8m9MlPltpc88Q++JbisTsjap+cdM8fWxnpevSuoUL/MqesswxdHo27fe0ru9r6DVVt+d7pWUslEolYAJxhNIiqWtnTarT6xhhGP5LP51u+d3r2zjOl0nnc0iHIygqhtTXe9PHtKFuZSXFocpyPVxw/Px/ZMe149pFCioMT97s7UGtPUXZH0jqRt5WpCyrPNz8f4ciR/kv05YWl0rN3nvlUgkFobQ1JpdDYdsi7xmJ8Z2zn1HR2dneBruHosOfy1HvOI0eyrjL1A6ZUqmBKxT/KI13f9PGjaCzGX7xnOxP8iQdPkKkIkz95MsnJkxWWxakZYKcDFbatjFrRrtvlTrfTHLzxNbvPWe/5+gVz1FbBlIphNEerpYJ7VqlYDWXDaA6zVKrgKBVbVjaMxmj1B7lnfSoigoh4EiHYz0yfqD+xdeQNzwLwnz5QttIzUfCdlK/+hIForP5cJdVWeSpXcI4dy3Dvvdt+kcT5Q679DHece8WUyh44NUxsKhQsJJUiXEcZiH5dffGLXC7HwMCAhenvRSgUsulPi2SOFiyMlft2JkceP3VqV/vLf+P1AHzo179aanMslPLVnzd9/CjUoVRcV4SaoDyc36yW6uTzeQYH69t/tRc9/RNu0x7DaIxYndPSvegLpWLWimHURzTaeo77np7+QGF5rN+dtRPveQ+Rc+d2tWeOHWPlXtcMnTWJPlmYzjjTIIAnpi8CO52ysrXbUSupFITDu5JP3zpxKx9/Q2Xwfn2cee4Mcy/M7WirVk0QbBpUiXOPeLG9pactFSiYc/1eAtVNofiJxmLk4vFd7eeuNS/n3AtzzK/M72hrJE/Krbf23z6fcrxy0kIALRUR+SDw88BSsekDxSxwTTE4OMjmZsdydgeaKw/tXU2wGcodtcc+UDj/h05tO2V/6r5bAVzD9B+uo2phIxwZP7Ijb8t9X4PKHCuGO145aSGASqXIb6nqb3hxooGBAS9OYxg9jxdOWuiD6Y/jeDJnrWHsjRdOWgiupfJLIvKvgbPAe1V1d+biOgmHw4TDYfL5vG0yrMHgmTPEKgpIuTlkAf7shhW+dGANyqJnHy9GzzpTnnbh5pStVm3PqE0+n/fMSQs+KRUR+RvgOpeX7gF+F/gQhYqEHwJ+E3iHyzn2rFAIsL6+ztzcHLlcjlwu19UrQK1UKHxVcQWm1vG3PfAAXLzI5sGDpbax4rGriZ3V7s7sW+Z7wxluWd3+dcvnq1cJLL/2TeGbyA3ldrXVIyPAA997gItbFzk4tC1ntWp7nSJoFSQbQVUJhUI899x2MnLn3mmGqhUKg4CIHAL+QlVftlc/twqFAHNzc8zOzpJMJnnxxRe72r/SSgU8Z/9OLUetW5RstWOd/T0f+vVtB6zb8W7VBKG18bhVFPSbbq5QmE6nmZmZYXh4O4GVc++UU2+FwsD5VETk+rKnbwGeavWcsVgMEdnzl9Qw+hHHqPBq5QeC6VO5T0SOUZj+nAfe1eoJRYR4PM7GxoZtLjSMMnK5HIODg57eF4FTKqr6tnacd3h4mPX19Xac2hfcnKq1KHe2Og5YZ8MgQGR+nuwRd2fnrhQIr4BzIynuKW4iBHhiZqEQ2Pbgk1RSGYdyU/gmPnHHJ/aU180hC+aU9ZJcLkfcJRCxFfrmZ9tZgw+yD6kRYnNzRObna3dsgOyRI6Qq5tHV+InLo9y6vjOuQWMx8qPeVfpzi5KFQpDb7I2znl2nX3FC84eGhjw9b+AslXYRDocZHBwknU73TPmO7JEju1ISuOFYGW4O2FrHV0t98JOnTvGTWVj51b2dsns5auuhMkrW8I58Pk80GvX8fugbSwUgHo/3/T4gw3Box9QH+kypOGZer0yBDKMVVNXzqQ/00fQHCkWSBgcHyWQyPTEFiszP74p0dRyw7b7OXk7dWpx+9DSJRKKmo9ZoH7lcjkgk0pbYre6/sxpkbGyMpaWl2h0DTr0OVYfsrTtD5zPHjrV0nUacupU8vvQ4WxXFxIzOks1m2bdvX1uizPtOqQwNDZUC4bo5ZiV58iTJkyd3tbtFtbpRb3Kmatcxuhdn1acd/hToM58KFJJhj42Nka0j8bJh9CKZTIaRkZG2bbDtO6UCMDIyApjD1ug/nO/8qIfxRJX03fQHCnkjhoaGSCaTnuWQ8Irx06eJPv74rvZXbW0RvvPOmse3w1HrxulHT/P40m45of4sbs7GwGo8eeXJqv1mb5zl5GGbljWKkzaynZtr+9JSgYLD1jYYdifzK/Ou4ftGbXK5HGNjY21NA9KXlgoUdmVGIhFyuVygkje5OVDrjX4t79tu7n3tbjmrRc9WcnT6KInlRM1I2WopDmpZOIY7zuJEeYqDdtC3loqIMD4+bhG2Rt+QzWYZGxtr+6pn3yoVoLSkZtMgo9dxHLTOIkU76WulEg6HmZiYsOVlo+fJZDKMjo52ZGGib30qDmNjY6ysrLQ9GG6vpNLlOU06tXrjNxe3Lu7wjTgrPW5U+lAsn0pjOJb4+Ph4R67X15YKFILhJicnyWQybb1OO/KfuJE5dqzuEHy/mL1xdkfS6kaxfCqNkc1mGR8f79h+N7+y6d8NfBC4DbhTVc+WvXYaeCeQA35ZVb/UbnlGRkZYWVlp+0pQZf6TvfKclCeabmRFp9nayF5wdPpo7U7AycMnuXnz5h2JooOYzLoXcMpvjI2NdeyaflkqTwE/AzxS3igitwNvBe4A3gj8joi0fb3XsVbMt2L0GtlslsnJyY6GTfiiVFT1GVX9rstLbwY+q6opVX0OOAfUDiP1gHg8TjQaNcVi9Ay5XI5QKNSRFZ9yguaoPQh8o+z5C8W2tiMiTE1NcfnyZcLhcHu2hLvkJYHOBax5yelHTwPuQXDNMr8yv8spa+H4zZPNZpmenu74bvy2KZW9qhCq6heqHebS5rrrr5EKhY2QyWRKW8O9ZOaWW9iXSEBZpb+xPXKKlOdw/cGbbiKXy/FsgCrgPfb8Y8BOORPLiV1t1ais6HeL3kIinyidAworRInlBDdv3uyV2G0jqBUKL1y40NRxrVQobJtSUdUfa+KwF4CXlD2/EbhU5fz3A/dDoUJhZTU1cK+yVotUKsWlS5cYGBjwVrHccQe8a2cJI7eCIY4ls6Pa3Sc+wXcCVgFv6B8LaQjLZZpKTO1qq0ZlRb87uIN3VZR4cqyWII27GkGqUKiqpaqDzeZMaebecQjakvIXgbeKSExEDgMvBf6ukwLEYjHGxsbavsRsGO3CyZfS7j0+1fBFqYjIW0TkBeA1wF+KyJcAVPVp4HPAt4G/At6tqh3fnDM5OUkkEjGnrdF1OM7Zqamptu5E3gtfHLWq+iDwYJXXPgx8uLMS7SQUCjE9Pc3CwkLJadtMRcBKUrOzO1Izup2zlYTSfuAWFVuZT+XY/mOeOnTrdRK3w5kcZFSVbDbLzMyMrzvvg7b6ExgGBwcZGxtjbW2NgYGBUkRssze8E01brlTcztlKQul+oVpyqGb79Qp+T3scTKnsweTkJFtbW6VpUL0VAd1wW0pu9Zx+4kTPlkfh+GdnAAAMcklEQVTA1ptPxfCeIEx7HILmqA0UzjQol8uB5bM1Aooz7dm3b18gEo6ZUqmBMw2ynCtGUAnKtMfBpj91MDk5yWYohLaoWCojap00B5VTo0qHrp+cee6Maz7Yak5ZL6l1jXZeu1sI0rTHwSyVOgiFQgy+/OWkDx9u2mJJzc7W5eSNzM+3vMrkJXMvzDG/0v6UDe3i1olba3fqUvL5fCBWeyoxS6VOIu97H7GNDVYWF5uKtnWr9OdWTbCaQ9dPjowf8SUlgZszGNzTJPSbk1hVyWQy7Nu3j8HBQb/F2YFZKg0Qj8eZmJggnU5bITLDN5ww/LGxsbYWBWsWUyoNMjExwfDwMOl02m9RjD4lm80yODgYKD9KOTb9aRARYf/+/SwsLJDJZFpOJFzpvPUzonavqoOt1NppNn1BtQqFnXASB5VsNksoFGJmZiaQCgXMUmkK50MFWqob5Oa87bWIWqsm6B35fJ58Ps+BAwcC5ZitxCyVJolGoxw4cICFhQVEpKlEOG7OWz9pxx4ZL6oJ1uMk7vWqhY5jdmZmpq11kL3ALJUWGBwcZHp6upTYyTDageOYnZiYaDo/SicxpdIiIyMjtiJktI3ylZ6JiQm/xakLUyotIiJMTEyYYjE8x1Eoo6OjgV3pccN8Kh7gKJZ8Ps/q6qr3qSi7hL1WYypf2y/7+f07fr/0fK/cJ9XO28uBbo5CGRkZYd++fV31fTKl4hFONn6grxVLvTy/+fyO5/2W+2QvHIUSj8eZnp7uuu9RoCoUisgh4BnAqQn0DVX9RR9EbApTLAXKLQi38PlG40sqLZJejk8pt1C6UaGAf5aKU6HwEy6vfU9Vg10MeA8cxSIirKys9K1iMRqn3IfSbVOecvzKUfsM0LVvWi1EhMnJSUSEa9euEY1GO17Qyeguyld5uskp60YQfSqHReRbwCrwH1X1a34L1AyO8zYcDnP16lUikUigoyC9pN7piVu/ettgZ8Db/Mo8R8a7J2F4Ofl8nkwmU1pF7GaFAsGrULgA3KSqV0XklcCficgdqrrqcv62VChsB07eC2jeOgtqBbxytuqsuLhf9nNBL+zZvx7KqxlOMMEteotv71Gzn48TghCJRAJlzbZy74ifcRUiMgf8quOobfR1h+PHj+vZs7u7tFJlzWsymQyLi4ulTYiNKpcgVcBrBLfcJ7B7PN2eI6XRz8fJKxsKhThw4EDgQu/d7h0R+aaqHq91bHBUIyAi+0UkXHx8hEKFwu5NO1ZGNBrluuuuY2hoyILk+hzHfzIwMMANN9wQOIXSKoGqUAi8DnhCRP4f8KfAL6pqotp5uo1wOMzMzEwp+taSafcf+Xy+5JAN+m7jZglUhUJV/Tzw+c5L1DkcB240GuXKlSuEQiEikSD6yw2vyeVypVIao6OjXe+QrYZ9m31ARBgZGSEajbK4uEgqlbJ4FheO7e/acKUdOGkLRITrr78+cDllvcaUio/EYjFuuOEGrl27xurqal8tO9dDL9RAdqyTkZERpqam+uLzNaXiM+FwmH379jE8PMyVK1fMaukRyq2TmZkZhoeH++YzNaUSEIaGhsxq6RH60Topx5RKgOhVq2V+ZX5XusfEcoKpxFTpuZPMupvpZ+ukHFMqAaTSaunmX7rZG2f9FqFjOLuL+9E6KceUSkApt1oSiQSqSi6X67ov68nDJ13Lc1RGoHZrOgNnqgP0tXVSTqAiao3dOFaLE8uSSqUsaC4A5PN5UqkUqsr09DTRaJR4PN73CgXMUukKnBIgBw8eZGNjg0QiUdpDFKRNaP2AE2IfCoWYmppidHTUPoMKTKl0EaFQiNHRUeLxOKurq6ysrJDNZpvaoOgnpx89zWPPP8bQPw75LUrdOE5YKJS+HRsb67qpaKcwpdKFhEIhJiYmGB0dZWVlhdXVQmaIoG2fb5RbJ271W4RdlKetGB0dZXx83LZV1MDenS4mHA4zNTXF2NgYa2trrK6ulrbTh8PhwFov97723kCncnCc4rlcjlAoxNjYGKOjoy3Xze4XTKn0AJFIhMnJScbHx9na2mJlZYVUKoWIdN3UyE+cKY6qMjAwwNTUFMPDw11t/fmBKZUeIhQKEY/HGR4eJpPJsLa2xtraGqpqEbp74FglAPF4nLGxsZ4IOvQLUyo9iIgwMDDAvn37mJiYYHNzc4f1Eg6HCYVCfXvTqGrJV+K8H1NTU8TjcVO8HmBKpccJh8OMjo4yMjJCOp1ma2uLjY2NkpkfDocD7X/xCsdP4sT4RCIRJiYmGB4eNqvEY0yp9AkiQiwWIxaLMTExQTabLSmYZDJZ6uNYMb2Akwc2n8+Xxh+PxxkaGiISiZgiaRN+VSj8KPBTQBr4HvBvVfVa8bXTwDuBHPDLqvqlqicymiYSiTA6Osro6Cj5fJ5kMsnm5iYbGxulJVSga6ZKzpTG8Y1AQUkODQ0xMjJCLBazqU2H8MtSeRg4rapZEfmvwGngfSJyO/BW4A7gBuBvROQHVDW3x7mMFgmFQgwPDzM8PMy+ffvIZrNkMhlSqRRbW1uk0+kd/f1WNG4KBGBgYIB4PE4sFiMajdrKl0/4laP2r8uefgP4F8XHbwY+q6op4DkROQfcCTzWYRH7FmcZOhqNMjw8zOTkZGka4SiaZDJJKpXadayqlrYUiMiOv3pQ1V1/jg+k8hymQIJLEHwq7wD+pPj4IAUl4/BCsc3wkUpFAzutBee/s2O3/K/a5kdnD00lTiLwcDhMNBotPXasI7+tJKM2vlYoFJF7gCzwh85hLv1dC+R0U4VCL+ilsUAh98j58+f9FsMzeu3zaWk8biZnJ/6At1OY1gyXtZ2m4Gtxnn8JeE2tc73yla9UN77yla+4tncjvTQWVRtP0HEbD3BW67i3/Som9kbgfcA/V9XNspe+CLxVRGIicphChcK/80NGwzCawy+fyv8AYsDDxbnxN1T1F1X1aRH5HPBtCtOid6ut/BhGV+HX6k/VPe6q+mHgwx0UxzAMDxHtgULhIrIEXHB5aRq40mFx2kUvjQVsPEHHbTw3q+r+Wgf2hFKphoicVdXjfsvhBb00FrDxBJ1WxtMbmzwMwwgMplQMw/CUXlcq9/stgIf00ljAxhN0mh5PT/tUDMPoPL1uqRiG0WF6TqmIyEdF5Dsi8oSIPCgiE2WvnRaRcyLyXRH5CT/lrBcRuVtEnhaRvIgcr3it68YDhYjqosznROT9fsvTKCLyKRFZFJGnytqmRORhEXm2+H/STxnrRUReIiJfEZFnit+zXym2Nz2enlMqFHK1vExVXw78I4X9RFTkankj8Dsi0g1Ze54CfgZ4pLyxW8dTlPG3gRPA7cDPFsfSTXyawntezvuBL6vqS4EvF593A1ngvap6G/Bq4N3Fz6Pp8fScUlHVv1ZVJ3XZN4Abi49LuVpU9TnAydUSaFT1GVX9rstLXTkeCjKeU9V5VU0Dn6Uwlq5BVR8BEhXNbwY+U3z8GeCnOypUk6jqgqr+Q/HxGvAMhXQjTY+n55RKBe8AHio+Pgh8v+y1bs/V0q3j6Va5a3FAVRegcKMCMz7L0zAicgh4BfC3tDCeICRpaph252rpNPWMx+0wl7ZAjKcG3Sp3TyMiI8DngX+nqqutJMHqSqWiqj+21+si8nbgTcCP6vaa+QvAS8q63Qhcao+EjVFrPFUI7Hhq0K1y1+KyiFyvqgsicj2w6LdA9SIiUQoK5Q9V9YFic9Pj6bnpTx/launW8fw98FIROSwiAxSczV/0WSYv+CKFxGMU/1ezMAOFFEyS3wOeUdWPlb3U/HjqyeTUTX8UHJbfBx4v/v3PstfuoVAS5LvACb9lrXM8b6Hw654CLgNf6ubxFOU+SWFl7nsUpni+y9Sg/H8MLACZ4mfzTmAfhVWSZ4v/p/yWs86xvJbC9POJsnvmZCvjsYhawzA8peemP4Zh+IspFcMwPMWUimEYnmJKxTAMTzGlYhiGp5hSMTpGcUfscyIyVXw+WXx+s9+yGd5hSsXoGKr6feB3gY8Umz4C3K+qbpUQjC7F4lSMjlIMCf8m8Cng54FXaGG3stEjdOXeH6N7UdWMiPwH4K+AHzeF0nvY9MfwgxMUwtxf5rcghveYUjE6iogcA/4ZhSxj/764A9boIUypGB2juCP2dynk7Hge+CjwG/5KZXiNKRWjk/w88LyqPlx8/jvAD4rI632UyfAYW/0xDMNTzFIxDMNTTKkYhuEpplQMw/AUUyqGYXiKKRXDMDzFlIphGJ5iSsUwDE8xpWIYhqf8f6HULqG79IvQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()  # initiate plots\n",
    "\n",
    "r = 20  # radius of circle boundary\n",
    "r2 = r**2  # square of radius\n",
    "\n",
    "for c in 'rgb':  # 3 random walks r-ed, g-reen, b-lue, 1st loop c=r, 2nd loop c=g, 3rd loop c=b \n",
    "    x, y = 0, 0  # start position 0\n",
    "    X, Y = [0], [0]  # store positions\n",
    "    walker = rand_walk_Z2()  # start new symmetric walk from 0, 0\n",
    "    while x**2 + y**2 < r2:\n",
    "        x, y = next(walker)\n",
    "        X.append(x)\n",
    "        Y.append(y)\n",
    "    #X = np.array(X) + np.random.randn(len(X))/3\n",
    "    #Y = np.array(Y) + np.random.randn(len(Y))/3\n",
    "    ax.plot(X, Y, color=c, alpha = 0.8)\n",
    "\n",
    "# To look nicer\n",
    "circle = plt.Circle((0, 0), r, color=\"#cccccc\", alpha=0.5)\n",
    "ax.add_artist(circle)  # draw a circle\n",
    "ax.set_xlim([-r-1, r+1])\n",
    "ax.set_ylim([-r-1, r+1])\n",
    "ax.set_aspect('equal')  # set aspect ratio 1:1\n",
    "ax.grid()  # adds grid lines\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "\n",
    "# To Show image!\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `alpha` parameter makes the lines more transparent to make the areas that the chain has visited multiple times a bit more obvious."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application: Output time distribution\n",
    "\n",
    "Next a simulation of the chain for `N = 1,000` times follows.\n",
    "For each simulation the time it took for the chain to exit the disk for the first time is stored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = 20  # radius of circle boundary\n",
    "r2 = r**2  # square of radius\n",
    "\n",
    "N = 1000  # sample size\n",
    "escape_time = np.zeros(N, dtype=int)  # to store the exit time\n",
    "for n in range(N):\n",
    "    walker = rand_walk_Z2()\n",
    "    x, y = (0, 0)\n",
    "    t = 0\n",
    "    while x**2 + y**2 < r2:\n",
    "        x, y = next(walker)\n",
    "        t += 1\n",
    "    escape_time[n] = t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that N=1000 random samples of the output time have been created an estimation of its distribution can be given, by plotting the histogram of its values in our sample.\n",
    "\n",
    "In the code cell below this histogram is plotted by dividing the range of values of `escape_time` into 30 regions (`bin=30`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEKCAYAAAAvlUMdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFZtJREFUeJzt3X+w5XV93/HnK8tC0qiUH6tuWXA3ssYskwnSlZiaSdMQ5UdpVycyWVINaUlJppBo02QA06GUCTOuk8jYqdoSoSI6LhRt3Di0iODodMYAqyHKQlZWwLCwyiIUTBt+7PLuH+ez4eRy7j1n1/vZe8/d52PmzPl+P9/P93M/nz2X++L743y+qSokSZpvP7TQHZAkLU0GjCSpCwNGktSFASNJ6sKAkSR1YcBIkrowYCRJXRgwkqQuDBhJUheHLXQHFtKxxx5bq1evXuhuSNJU+epXv/p4Va0YV++QDpjVq1ezdevWhe6GJE2VJN+epJ6nyCRJXRgwkqQuDBhJUhcGjCSpCwNGktSFASNJ6sKAkSR1YcBIkrowYCRJXRgwna1cdQJJxr5WrjphobsqSfPqkJ4q5mD4ziMP85qLPze23rc3nX0QeiNJB49HMJKkLgwYSVIXBowkqQsDRpLUhQEjSerCgJEkdWHASJK6MGAkSV0YMJKkLgwYSVIXBowkqQsDRpLUhQGzWCxbPtGsy868LGlaOJvyYrH3+YlmXQZnXpY0HTyCkSR1YcBIkrowYCRJXRgwkqQuugZMkjOSbE+yI8klI7YfkeSGtv2OJKuHtl3ayrcnOb2VHZ/ki0nuS7ItybuH6h+d5NYk97f3o3qOTZI0t24Bk2QZ8CHgTGAdcG6SdTOqnQ88WVUnAlcBm9q+64CNwEnAGcCHW3t7gH9XVT8BvAm4cKjNS4DbqmotcFtblyQtkJ5HMKcCO6rqgap6DtgMbJhRZwNwXVu+CTgtSVr55qp6tqoeBHYAp1bVrqr6GkBVfR+4DzhuRFvXAW/rNC5J0gR6BsxxwMND6zt5MQxeUqeq9gBPAcdMsm87nfYG4I5W9Kqq2tXa2gW8ch7GIEk6QD0DJiPKasI6c+6b5GXAp4H3VNXT+9Wp5IIkW5Ns3b179/7sKknaDz0DZidw/ND6KuDR2eokOQw4Enhirn2TLGcQLp+sqs8M1flukpWtzkrgsVGdqqqrq2p9Va1fsWLFAQ5NkjROz4C5C1ibZE2SwxlctN8yo84W4Ly2/A7g9qqqVr6x3WW2BlgL3Nmuz1wD3FdVH5ijrfOAz877iCRJE+s2F1lV7UlyEXALsAy4tqq2JbkC2FpVWxiExfVJdjA4ctnY9t2W5EbgXgZ3jl1YVXuT/CzwLuAbSe5uP+q9VXUz8D7gxiTnA38FnNNrbJKk8bpOdtn+8N88o+yyoeVnmCUIqupK4MoZZf+b0ddnqKrvAaf9gF2WJM0Tv8kvSerCgJEkdWHASJK6MGAkSV0YMJKkLgwYSVIXBowkqQsDRpLUhQEjSerCgJEkdWHASJK6MGAkSV0YMJKkLgwYSVIXBowkqQsDRpLUhQEjSerCgJEkdWHASJK6MGAkSV0YMJKkLgwYSVIXBowkqQsDRpLUhQEjSerCgJEkdWHASJK6MGAkSV0YMJKkLgwYSVIXBowkqQsDRpLUhQEjSerCgJEkdWHASJK6MGCm0bLlJBn7WrnqhIXuqaRD2GEL3QEdgL3P85qLPze22rc3nX0QOiNJo3U9gklyRpLtSXYkuWTE9iOS3NC235Fk9dC2S1v59iSnD5Vfm+SxJPfMaOvyJI8kubu9zuo5NknS3LoFTJJlwIeAM4F1wLlJ1s2odj7wZFWdCFwFbGr7rgM2AicBZwAfbu0BfKyVjXJVVZ3cXjfP53gkSfun5xHMqcCOqnqgqp4DNgMbZtTZAFzXlm8CTkuSVr65qp6tqgeBHa09qurLwBMd+y1Jmgc9A+Y44OGh9Z2tbGSdqtoDPAUcM+G+o1yU5OvtNNpRoyokuSDJ1iRbd+/ePdlIJEn7rWfAZERZTVhnkn1n+gjwWuBkYBfwR6MqVdXVVbW+qtavWLFiTJOSpAPVM2B2AscPra8CHp2tTpLDgCMZnP6aZN+/o6q+W1V7q+oF4I9pp9QkSQujZ8DcBaxNsibJ4Qwu2m+ZUWcLcF5bfgdwe1VVK9/Y7jJbA6wF7pzrhyVZObT6duCe2epKkvrr9j2YqtqT5CLgFmAZcG1VbUtyBbC1qrYA1wDXJ9nB4MhlY9t3W5IbgXuBPcCFVbUXIMmngJ8Hjk2yE/gPVXUN8P4kJzM4lfYQ8Bu9xiZJGq/rFy3brcI3zyi7bGj5GeCcWfa9ErhyRPm5s9R/1w/UWUnSvHKqGElSFwaMJKkLA0aS1IUBI0nqwoCRJHVhwEiSupgoYJJ8Osk/TWIgSZImMmlgfAT4FeD+JO9L8vqOfdJ88cmXkhbQRF+0rKovAF9IciRwLnBrkocZzPn1iap6vmMfdaB88qWkBTTxKa8kxwC/Bvw68OfAB4FTgFu79EySNNUmOoJJ8hng9cD1wD+rql1t0w1JtvbqnCRpek06F9lHZz6COMkR7YmT6zv0S5I05SY9RfYHI8q+Mp8dkSQtLXMewSR5NYNHFf9Ikjfw4pMmXwH8vc59kyRNsXGnyE5ncGF/FfCBofLvA+/t1CdJ0hIwZ8BU1XXAdUl+qao+fZD6JElaAsadIntnVX0CWJ3kd2Zur6oPjNhNkqSxp8h+tL2/rHdHJElLy7hTZP+1vf/Hg9MdSdJSMelkl+9P8ooky5PcluTxJO/s3TlJ0vSa9Hswb62qp4GzgZ3A64Df69YrSdLUmzRglrf3s4BPVdUTnfojSVoiJp0q5k+T/CXwN8C/SbICeKZftyRJ026iI5iqugT4GWB9m5r//wIbenZMkjTdJj2CAfgJBt+HGd7n4/PcH0nSEjHpdP3XA68F7gb2tuLCgJEkzWLSI5j1wLqqqp6dkSQtHZPeRXYP8OqeHZEkLS2THsEcC9yb5E7g2X2FVfXPu/RKkjT1Jg2Yy3t2QpK09Ex6m/KXgIeA5W35LuBrHfu16K1cdQJJxr4k6VA16V1k/xq4ADiawd1kxwH/BTitX9cWt+888jCvufhzY+t9e9PZB6E3krT4THqR/0LgzcDTAFV1P/DKXp2SJE2/SQPm2ap6bt9K+7KltyxLkmY1acB8Kcl7gR9J8hbgvwN/2q9bkqRpN2nAXALsBr4B/AZwM/Dve3VKkjT9JrrIX1UvJPkT4E+qanfnPkmSloA5j2AycHmSx4G/BLYn2Z3kskkaT3JGku1JdiS5ZMT2I5Lc0LbfkWT10LZLW/n2JKcPlV+b5LEk98xo6+gktya5v70fNUkfJUl9jDtF9h4Gd4+9saqOqaqjgZ8G3pzk3861Y5JlwIeAM4F1wLlJ1s2odj7wZFWdCFwFbGr7rgM2AicBZwAfbu0BfKyVzXQJcFtVrQVua+uSpAUyLmB+FTi3qh7cV1BVDwDvbNvmciqwo6oeaHegbealz5DZAFzXlm8CTsvg24kbgM1V9Wz72Ttae1TVl4FRT9Qcbus64G1j+idJ6mhcwCyvqsdnFrbrMMtH1B92HPDw0PrOVjayTlXtAZ4Cjplw35leVVW7Wlu7mOV7OkkuSLI1ydbdu72cJEm9jAuY5w5wG8CoeVJmfndmtjqT7HtAqurqqlpfVetXrFgxH01KkkYYdxfZTyV5ekR5gB8es+9O4Pih9VXAo7PU2dm+vHkkg9Nfk+w703eTrKyqXUlWAo+NqS9J6mjOI5iqWlZVrxjxenlVjTtFdhewNsmaJIczuGi/ZUadLcB5bfkdwO3toWZbgI3tLrM1wFrgzjE/b7it84DPjqkvSepo0i9a7rd2TeUi4BbgPuDGqtqW5Iok+54jcw1wTJIdwO/Q7vyqqm3AjcC9wP8CLqyqvQBJPgV8BfjxJDuTnN/aeh/wliT3A29p65KkBTLp82AOSFXdzOBb/8Nllw0tPwOcM8u+VwJXjig/d5b63+MQnt1ZkhabbkcwkqRDmwEjSerCgJEkdWHACJYtn+jxzytXnbDQPZU0Rbpe5NeU2Pu8j3+WNO88gpEkdWHASJK6MGAkSV0YMJKkLgwYSVIXBowkqQsDRpLUhQEjSerCgJEkdWHASJK6MGAkSV0YMJKkLgwYSVIXBowkqQsDRpLUhQEjSerCgJEkdWHASJK6MGAkSV0YMJKkLgwYSVIXBowkqQsDRpLUhQEjSerCgJEkdWHASJK6MGAkSV0YMJKkLgwYSVIXBowkqQsDRpNbtpwkY18rV52w0D2VtAgc1rPxJGcAHwSWAR+tqvfN2H4E8HHgHwLfA365qh5q2y4Fzgf2Ar9dVbfM1WaSjwH/GHiqNf9rVXV3z/EdcvY+z2su/tzYat/edPZB6Iykxa5bwCRZBnwIeAuwE7gryZaquneo2vnAk1V1YpKNwCbgl5OsAzYCJwH/APhCkte1feZq8/eq6qZeY5IkTa7nKbJTgR1V9UBVPQdsBjbMqLMBuK4t3wScliStfHNVPVtVDwI7WnuTtClJWgR6BsxxwMND6ztb2cg6VbWHwemtY+bYd1ybVyb5epKr2uk3SdIC6RkwGVFWE9bZ33KAS4HXA28EjgYuHtmp5IIkW5Ns3b1796gq+kF5M4Ak+l7k3wkcP7S+Cnh0ljo7kxwGHAk8MWbfkeVVtauVPZvkvwG/O6pTVXU1cDXA+vXrZwae5oM3A0ii7xHMXcDaJGuSHM7gov2WGXW2AOe15XcAt1dVtfKNSY5IsgZYC9w5V5tJVrb3AG8D7uk4NknSGN2OYKpqT5KLgFsY3FJ8bVVtS3IFsLWqtgDXANcn2cHgyGVj23dbkhuBe4E9wIVVtRdgVJvtR34yyQoGp9HuBn6z19gkSeN1/R5MVd0M3Dyj7LKh5WeAc2bZ90rgyknabOW/8IP2V5I0f/wmvySpCwNGktSFASNJ6sKAkSR1YcBIkrowYCRJXRgwkqQuDBhJUhcGjCSpCwNGktSFASNJ6sKAkSR1YcBIkrowYCRJXRgwWjg+Wlla0ro+D0aak49WlpY0j2C0+E14pOPRjrS4eASjxW/CIx3waEdaTDyCkSR1YcBIkrowYCRJXRgwkqQuDBhJUhcGjJYWv7wpLRrepqylxS9vSouGRzCSpC4MGElSFwaMDk1eq5G68xqMDk1eq5G68whGktSFASPNxVNp0gHzFJk0F0+lSQfMIxhJUhcGjCSpCwNGmg9eq5Fewmsw0nyY9FrNH76dJGPrvfq449m186/mo2fSgjFgpIPJmwZ0COl6iizJGUm2J9mR5JIR249IckPbfkeS1UPbLm3l25OcPq7NJGtaG/e3Ng/vOTZJ0ty6BUySZcCHgDOBdcC5SdbNqHY+8GRVnQhcBWxq+64DNgInAWcAH06ybEybm4Crqmot8GRrW9J+WrnqBK8naV70PEV2KrCjqh4ASLIZ2ADcO1RnA3B5W74J+M8ZnKDeAGyuqmeBB5PsaO0xqs0k9wG/APxKq3Nda/cjfYYmddZuGhinx7Wa7zzy8KI+jbdy1Ql855GHJ6rrtayF1TNgjgOGfwt2Aj89W52q2pPkKeCYVv5nM/Y9ri2PavMY4P9U1Z4R9aXp0+Fazf78YZ7IhCG47PAfZu9zz4ytN2kYTBqAcOhdy1ps4Zuq6tNwcg5welX9elt/F3BqVf3WUJ1trc7Otv4tBkcqVwBfqapPtPJrgJsZnNJ7SZtD9U9s5ccDN1fVT47o1wXABW31x4HtcwzjWODxA/sXmDqOdWlyrEvXQo73NVW1YlylnkcwO4Hjh9ZXAY/OUmdnksOAI4Enxuw7qvxx4O8nOawdxYz6WQBU1dXA1ZMMIMnWqlo/Sd1p51iXJse6dE3DeHveRXYXsLbd3XU4g4v2W2bU2QKc15bfAdxeg0OqLcDGdpfZGmAtcOdsbbZ9vtjaoLX52Y5jkySN0e0Ipl1TuQi4BVgGXFtV25JcAWytqi3ANcD17SL+EwwCg1bvRgY3BOwBLqyqvQCj2mw/8mJgc5I/AP68tS1JWiDdrsEsBUkuaKfUljzHujQ51qVrGsZrwEiSunCyS0lSFwbMCOOmuJlGSR5K8o0kdyfZ2sqOTnJrm17n1iRHtfIk+U9t/F9PcsrC9n68JNcmeSzJPUNl+z2+JOe1+vcnOW/Uz1pos4z18iSPtM/37iRnDW3br2mXFpMkxyf5YpL7kmxL8u5WvuQ+2znGOr2fbVX5GnoxuHngW8CPAYcDfwGsW+h+zcO4HgKOnVH2fuCStnwJsKktnwX8TyDAm4A7Frr/E4zv54BTgHsOdHzA0cAD7f2otnzUQo9twrFeDvzuiLrr2u/wEcCa9ru9bFp+z4GVwClt+eXAN9uYltxnO8dYp/az9Qjmpf52ipuqeg7YN8XNUrSBwbQ6tPe3DZV/vAb+jMF3jFYuRAcnVVVfZnAn4rD9Hd/pwK1V9URVPQncymAuvEVllrHO5m+nXaqqB4F90y5Nxe95Ve2qqq+15e8D9zGYpWPJfbZzjHU2i/6zNWBeatQUN0th2pkCPp/kqxnMZgDwqqraBYNfbuCVrXyp/Bvs7/imfdwXtdNC1+47ZcQSGmsGs62/AbiDJf7ZzhgrTOlna8C81KjJlZbCrXZvrqpTGMxEfWGSn5uj7lL9N9hntvFN87g/ArwWOBnYBfxRK18SY03yMuDTwHuq6um5qo4om6rxjhjr1H62BsxLTTLFzdSpqkfb+2PA/2BwGP3dfae+2vtjrfpS+TfY3/FN7bir6rtVtbeqXgD+mBdnH5/6sSZZzuAP7ier6jOteEl+tqPGOs2frQHzUpNMcTNVkvxokpfvWwbeCtzD352qZ3h6nS3Ar7Y7ct4EPLXvdMSU2d/x3QK8NclR7TTEW1vZojfjGtnbGXy+sJ/TLh3MPk8iSRjMynFfVX1gaNOS+2xnG+tUf7YLfefEYnwxuBPlmwzuxPj9he7PPIznxxjcSfIXwLZ9Y2LwmIPbgPvb+9GtPAwe7PYt4BvA+oUewwRj/BSD0wfPM/g/uPMPZHzAv2JwsXQH8C8Xelz7Mdbr21i+zuCPycqh+r/fxrodOHOofNH/ngM/y+D0zteBu9vrrKX42c4x1qn9bP0mvySpC0+RSZK6MGAkSV0YMJKkLgwYSVIXBowkqYtuT7SU9KIk+26rBXg1sBfY3db/X1X9owXpmNSRtylLB1mSy4G/rqo/XOi+SD15ikxaYEn+ur3/fJIvJbkxyTeTvC/Jv0hyZwbP8nltq7ciyaeT3NVeb17YEUijGTDS4vJTwLuBnwTeBbyuqk4FPgr8VqvzQeCqqnoj8Ettm7ToeA1GWlzuqjbvW5JvAZ9v5d8A/klb/kVg3WDqKgBekeTlNXiGiLRoGDDS4vLs0PILQ+sv8OJ/rz8E/ExV/c3B7Ji0vzxFJk2fzwMX7VtJcvIC9kWalQEjTZ/fBta3JxzeC/zmQndIGsXblCVJXXgEI0nqwoCRJHVhwEiSujBgJEldGDCSpC4MGElSFwaMJKkLA0aS1MX/BzwmzdOhzZ3rAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(escape_time, bins=30, edgecolor='black', density=True)\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Density')\n",
    "plt.show()"
   ]
  }
 ],
 "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": 4
}