Call center data modeling.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3 (Anaconda 5)",
      "env": {
        "LD_LIBRARY_PATH": "/ext/anaconda5/lib",
        "PROJ_LIB": "/ext/anaconda-2019.03/share/proj",
        "PYTHONHOME": "/ext/anaconda5/lib/python3.5",
        "PYTHONPATH": "/ext/anaconda5/lib/python3.5:/ext/anaconda5/lib/python3.5/site-packages"
      },
      "language": "python",
      "metadata": {
        "cocalc": {
          "description": "Python/R distribution for data science",
          "priority": -1,
          "url": "https://www.anaconda.com/distribution/"
        }
      },
      "name": "anaconda5"
    },
    "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.7"
    },
    "colab": {
      "name": "[Official] CS146 Assignment 2.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Hu41CcFa8c9A",
        "colab_type": "text"
      },
      "source": [
        "Link if people prefer Google Colab Notebook: https://colab.research.google.com/drive/1_Z1RehijDeqJ6fZzOaHGQCeSsamxWxhX "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "md6p3e6FPanf",
        "colab_type": "text"
      },
      "source": [
        "# Section 1: Call center data modeling\n",
        "\n",
        "You are advising a client on the number of call center agents they need to handle their customer support load. We start by modeling the number of phone calls we expect to come into the call center during each hour of the day.\n",
        "\n",
        "You are provided with a data set of the intervals between phone calls arriving during one day. Each value is a time in minutes indicating the amount of time that passed between receiving two consecutive phone calls."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ehihE-OtqYvo",
        "colab_type": "text"
      },
      "source": [
        "## Loading and initializing the data"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dRjqzkP0Pang",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# importing packages\n",
        "\n",
        "import numpy as np\n",
        "import scipy as sp\n",
        "import scipy.stats as sts\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "%matplotlib inline"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FSkSBmk_Pank",
        "colab_type": "code",
        "outputId": "2a75d1c4-2710-44e9-d5c8-512504bb1749",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 70
        }
      },
      "source": [
        "# Load the data set containing durations between calls arriving at the call\n",
        "# center during 1 day. All values are in minutes.\n",
        "waiting_times_day = np.loadtxt('https://docs.google.com/spreadsheets/d/e/2PACX-1vTDuQCIYts48eH8mNrw60_c5H0cjRzhwxD2OLveOaN5DnelQ99s22oU2rzRyLS0Oh1iNlEatp68fdmv/pub?gid=485723753&single=true&output=csv')\n",
        "\n",
        "# Display some basic information about the data set.\n",
        "print('Size of data set:', len(waiting_times_day))\n",
        "print('First 3 values in data set:', waiting_times_day[:3])\n",
        "print('Sum of data set:', sum(waiting_times_day))\n",
        "\n",
        "# Make 24 empty lists, one per hour.\n",
        "waiting_times_per_hour = [[] for _ in range(24)]\n",
        "\n",
        "# Split the data into 24 separate series, one for each hour of the day.\n",
        "current_time = 0\n",
        "for t in waiting_times_day:\n",
        "    current_hour = int(current_time // 60)\n",
        "    current_time += t\n",
        "    waiting_times_per_hour[current_hour].append(t)"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Size of data set: 5856\n",
            "First 3 values in data set: [30.   3.4  3.2]\n",
            "Sum of data set: 1441.6838153800093\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UWRNCcZuqz-s",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "outputId": "22288cc4-c832-4058-a708-acc7006776cc"
      },
      "source": [
        "number_call = [len(waiting_time_hour) for waiting_time_hour in waiting_times_per_hour]\n",
        "print(\"Percentage of calls between 8 AM and 8 PM: \",\n",
        "      np.sum(number_call[8:21])/np.sum(number_call))"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Percentage of calls between 8 AM and 8 PM:  0.922301912568306\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "l8cj6UTtPan7",
        "colab_type": "text"
      },
      "source": [
        "## 1. Building important functions\n",
        "\n",
        "1. The **likelihood** function, which is P(datum | lambda_), and this likelihood follow an exponential distribution: the exponential distribution of the datum given the parameter value\n",
        "2. The **prior** function, which is P(lambda_), and this likelihood follow a Gamma distribution: the distribution Gamma(λ | α, β).\n",
        "3. The **compute_posterior** function, which uses Bayes' Theorem.\n",
        "4. The **compute_percentile** function, which compute the parameter value at a particular percentile of the given probability distribution values."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Vm8QCkIjPan3",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# THE LIKELIHOOD FUNCTION\n",
        "#\n",
        "# This function takes 2 inputs, namely the parameter (λ) value at which to\n",
        "# compute the likelihood and a value from the data set. You need to evaluate\n",
        "# the exponential distribution of the datum given the parameter value.\n",
        "#np.array([np.sum(np.log(likelihood(param, data))) for param in parameter_values])\n",
        "def likelihood(lambda_, datum):\n",
        "    return sts.expon(scale=1/lambda_).pdf(datum)\n",
        "\n",
        "#####################################################\n",
        "# THE PRIOR DISTRIBUTION\n",
        "#\n",
        "# This function takes 1 input, namely the parameter value (λ) at which to\n",
        "# compute the prior probability density. You need to evaluate the distribution\n",
        "# Gamma(λ | α, β).\n",
        "\n",
        "def prior(lambda_):\n",
        "    return sts.gamma(1, scale=1/0.25).pdf(lambda_)\n",
        "\n",
        "#####################################################\n",
        "# THE POSTERIOR DISTRIBUTION: T\n",
        "#\n",
        "# The function below is provided to help make computing the posterior easier. It\n",
        "# follows the same pattern as in the previous class, where we multiple the prior\n",
        "# and the likelihood evaluated at various values of the parameter (λ).\n",
        "#\n",
        "# You will see the function uses the logarithms of various distributions. This\n",
        "# is for numerical reasons. When you multiply lots of likelihood values, the\n",
        "# result can become very large or very small, causing numerical overflow or\n",
        "# underflow in Python. Taking logarithms avoids this problem.\n",
        "def compute_posterior(parameter_values, data):\n",
        "    # Notes: The detail of this function will be presented in Section 2 \n",
        "    # (stretch goal: problem 2)\n",
        "    log_prior = np.log(prior(parameter_values))\n",
        "    log_likelihood = np.array([np.sum(np.log(likelihood(param, data))) for param in parameter_values])\n",
        "    unnormalized_log_posterior = log_prior + log_likelihood\n",
        "    unnormalized_log_posterior -= max(unnormalized_log_posterior)\n",
        "    unnormalized_posterior = np.exp(unnormalized_log_posterior)\n",
        "    area = sp.integrate.trapz(unnormalized_posterior, parameter_values)\n",
        "    posterior = unnormalized_posterior / area\n",
        "    return posterior\n",
        "\n",
        "#####################################################\n",
        "# THE PERCENTILE FUNCTION\n",
        " \n",
        "# Compute the parameter value at a particular percentile of the given\n",
        "# probability distribution values. This function uses the cumulative trapezoid\n",
        "# integrator in SciPy.\n",
        "#\n",
        "# Inputs:\n",
        "#\n",
        "# parameter_values (array of float) This is the list of parameter values\n",
        "# at which the probability distribution has been evaluated.\n",
        "#\n",
        "# distribution_values (array of float) This is the list of values of the\n",
        "# probability density function evaluated at the parameter values above.\n",
        "# \n",
        "# percentile (float) This is the value between 0 and 1 of the percentile \n",
        "# to compute.\n",
        "#\n",
        "# Returns: (float) The parameter value at the given percentile.\n",
        "def compute_percentile(parameter_values, distribution_values, percentile):\n",
        "    cumulative_distribution = sp.integrate.cumtrapz(\n",
        "        distribution_values, parameter_values)\n",
        "    percentile_index = np.searchsorted(cumulative_distribution, percentile)\n",
        "    return lambdas[percentile_index]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "u0876HGRPan8",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "heF2IEVSPan_",
        "colab_type": "text"
      },
      "source": [
        "## 2. Compute and Visualize the results\n",
        "\n",
        "1. Compute a 98% posterior confidence interval over the number of calls per minute for every hour of the day: 24 intervals total\n",
        "2. Compute the posterior mean of λ for each hour of the day\n",
        "3. Visualize the results above using an errorbar() plot: Each dot is at the posterior mean and each line shows a 98% confidence interval for a λ"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6K_gOYBgPaoH",
        "colab_type": "code",
        "outputId": "06f3b1e4-5af7-4881-cb09-a19bc23209b4",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 404
        }
      },
      "source": [
        "# Create empty lists for storing the data\n",
        "#\n",
        "# upper_percentiles_diff to store the difference between the 99% \n",
        "# and the posterior mean\n",
        "#\n",
        "# lower_percentiles_diff to store the difference between the 01% \n",
        "# and the posterior mean\n",
        "#\n",
        "# lambda_ to store the posterior mean\n",
        "upper_percentiles_diff = []\n",
        "lower_percentiles_diff = []\n",
        "lambda_ = []\n",
        "\n",
        "# iterate through every hour in the 24 hours series data\n",
        "for waiting_times_hour in  waiting_times_per_hour:\n",
        "    # initialize the lambdas values\n",
        "    lambdas = np.linspace(0, 20, 501)[1:]\n",
        "    #\n",
        "    # compute the posterior by using the compute_posterior function\n",
        "    # with the initialized lambdas and the data (of each hour)\n",
        "    posterior = compute_posterior(lambdas, waiting_times_hour)\n",
        "    #\n",
        "    # Add the newly calculated posterior mean to the lambda_ list\n",
        "    lambda_.append(sp.integrate.trapz(posterior * lambdas, lambdas))\n",
        "    #\n",
        "    # Calculate the difference between the posterior mean with the 99% percentile\n",
        "    # and append to the upper_percentiles_diff list\n",
        "    upper_percentiles_diff.append(compute_percentile(lambdas, posterior, 0.99) - lambda_[-1])\n",
        "    #\n",
        "    # Calculate the difference between the posterior mean with the 1% percentile\n",
        "    # and append to the upper_percentiles_diff list\n",
        "    lower_percentiles_diff.append(lambda_[-1] - compute_percentile(lambdas, posterior, 0.01))\n",
        "\n",
        "# initialize the graph\n",
        "plt.figure(figsize=(8, 6))\n",
        "# plot the error bar with the mean at lambda_ and the range from the mean derives\n",
        "# from the upper_percentiles_diff and lower_percentiles_diff\n",
        "plt.errorbar(list(range(24)), lambda_, yerr = [upper_percentiles_diff,lower_percentiles_diff]\n",
        "             , fmt=\"o\", capsize=3)\n",
        "# label the axes and plot\n",
        "plt.ylabel('Call rate (Calls per minute)')\n",
        "plt.xlabel('Hour of the day')\n",
        "plt.title('Mean and 98% confidence interval over call rate')\n",
        "plt.show()"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfsAAAGDCAYAAAAs+rl+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucHFWZ8PHfQwg6YiQoEU0gBATi\noijREURcF9QY9FWJrC4iKnhD31ddzLpxiesq4rrEzbsCivfLiyAiq4aIgmIWULwgmBAgyE0k3Abk\nFqIR44Ykz/tH1bidTs9Mz0xfZmp+38+nP+k6VV31dHVnnj6nTp0TmYkkSaqu7bodgCRJai+TvSRJ\nFWeylySp4kz2kiRVnMlekqSKM9lLklRxJntVWkTcHhEv7XYcnRARPRHxvYj4fUR8KyKOiYgfDbL9\njyPi7Z2McTAR8fmI+Jdux9GMqn2vIiIjYu/y+ZkR8a/djkmtZbKf4Mo/WhsjYpe68lXlH4BZ3Yms\n/SJiakR8LSLuLx8n1a0/ICJ+WibPu2sTUUTsHhG/jIi1EfEfda/7QUT0duht1HotsCvwpMx8XWae\nk5kv60IcI5KZ78rMjzWzrQlp/KjaD6PxymQvgDXA0f0LEbE/8LjuhdMxp1K8z1nAgcCbIuItNeu/\nAVwOPBH4G+D/RMSry3WLgK8BewLz+5N7RBwFrMnMFR15B1vbA7glMzd14djjSkRs3+0YRmOsxD9W\n4tDQTPYCOBt4c83yscBZtRtExGMi4v9GxJ0RcV/Z5NpTrts5Ir4fEQ9ExMPl891qXvvjiPhYRPw8\nItZHxI/qWxJqth3VviLiTRFxR0Q8FBH/PMT7fhXw75n5p8y8HfgK8Naa9bOAczJzc2b+FvgZ8Ixy\n3Z7ApZn5e+BXwF4R8QTgROCDQxyXiHhGRCwvWwbui4gP1pzn0yLinvJxWkQ8plx3aNnC8P6yJeLe\n/h8nEfFR4MPAURHxx4h4W0QcFxE/qznm3Ii4qWypOAOIupjeGhE3luf94ojYo2ZdRsS7IuI3EbEu\nIj4TEVGz/h3la9dHxA0R8ZyyfHpEfKf8PNdExN8Pck7+Ulsf4r0eDxwDfKB8r98b6lgRcVJEfDsi\nvh4RfwA+GBEbIuKJNdvMiYgHI2JyRDwtIi4tv0cPRsQ5ETF1qM+13M9OEXFWGccdEfGhiNiu/GzX\nRcQza7adVsbx5HL5lRFxTbndLyLiWTXb3h4R/xQR1wGPRINEO8j36sCIuKLc770RcUZE7NDM+6nb\n/3Hl/71TI+Ih4KTBzlVEnA3MBL5XflYfKMufX76/dRFxbUQcOtxYNEyZ6WMCP4DbgZcCNwN/BUwC\n7qaoJSYwq9zuVOACilruFOB7wCnluicBf0tRS54CfAtYVnOMHwO/BfYFesrlxQPEM+J9AfsBfwRe\nBDwG+CSwCXjpAMd6EDiwZvmfgYdrlv8NWAxMBmaX5+V55bolwHuAqcBvKH4EnA4c28Q5nwLcC7wf\neGy5fFC57mTgl8CTgWnAL4CPlesOLd/PyWVMrwD+BOxcrj8J+HrNcY4DflY+3wVYT9HUPxlYUO7r\n7eX6I4Bby+/A9sCHgF/U7CuB75fvdybwAHB4ue51QB/wPIofEHtTfH+2A1ZS/AjZAdgLuA2YN8B5\nORP41ybf61+2LZcHPVZ5bh4F5pfb9gCXAu+o2ccS4PPl872BuRTfo2kULTyn1f+/GeB9nAV8t/xc\nZwG3AG8r130V+HjNtu8Gflg+nwPcDxxE8f/w2PI4j6k55jXA7kDPML9XzwWeX362s4AbgffVfb57\nNzq3dcc4rvxc3lvuq2e45wqYATxUfqbbla99CJjW7b+HVX50PQAfXf4C/E+y/xBwCnA4sLz8j5zl\nH4YAHgGeVvO6gymaqxvt8wC2Tpo/Bj5Us/x/+v/ANRFf0/ui+EP/zZp1OwIbB/mj/HVgaflHcW+K\nHxH/XbP+BRQJcFN5Lj5as+6JwHnAtRSJcw5wWVne3/z/ngGOezSwaoB1vwVeUbM8D7i9fH4osAHY\nvmb9/cDzy+cnMXCyfzPwy5p1QfHjpT/Z/4AyIZXL21Ek1z3K5QReWLP+P4ETy+cXAyc0eC8HAXfW\nlS0C/t8A7/1Mtk72g73Xv2zbzLHKc3N53fq3U7TO9J+Pu4AXDRDb/NrPjAGSPUWS3gjsV1P2TuDH\n5fOXAr+tWfdz4M3l889R/rCrWX8z8Dc1x3zrIP9XBvxeNdj2fcD5NcvDSfZ3DrHvQc8V8E/A2XWv\nuZgmfij7GPnD6y3qdzZFgtqTuiZ8il/rjwNW1rbcUvxhIyIeR1HzPxzYuVw/JSImZebmcvl3Nfv7\nE/D4RkGMcl/TKf5gA5CZj5RNjQP5e+DTFDXzh4BzKfsulM27P6SovX8DeArw7Yi4LzM/m5lrgaPK\nbbejOHfvomjGv57ij+LVEXFJZt5Yd9zdKZJ6I9OBO2qW7yjL+j2UW1+TH/BcNthv7bnJiLirZv0e\nwOmxdWfDoKiF9ccz0Hkf6P3sAUyPiHU1ZZOAnzYRLwzvvTZzrLu2fgnfAT4dEU+laCna0r99ROxK\n0VLz1xQ/BrcDHm4i5l0oWiLqP8MZ5fPLgMdFxEHAfRQ/Zs+veQ/HRsR7a167A1t//vXvodaA36uI\n2JeipauX4v/y9hQtISOxVQwjOFd7AK+LiFfVlE2mODdqE6/ZC4DMvIOio94rKGq7tR6kqGU9IzOn\nlo+dMrP/D+/7KZq5D8rMJ1A0o0PdNeEmjWZf91L8wSteUPxweNJAG2fm2sw8JjOfkpnPoPj/cFW5\nei9gc2aelZmbMvNu4JsU56fe8RS15uuB/YEVmbkRWF0u17ur3H8j91D8Mew3sywbrfpzE7XLZUzv\nrPl8p2ZmT2b+ool93wU8bYDyNXX7nJKZjc7hcNVP19nMsbZ6TWY+DPyI4kfbGyhahfq3+bdy+/3L\n7+Ebae47+CDF5YL6z7CvPOZmilaRo8vH9zNzfc17+Hjde3hcZp47yPuuPwcDfa8+B9wE7FO+nw82\n+X4aqY9hqHPV6LM6u+597piZi0cYj5pgslettwEvzsxHagszcwvwJeDUmo5EMyJiXrnJFIofA+vK\nGvFHRhHDaPb1beCVEfHCsvPRyQzyHS87Fj0pIiZFxMspknb/7Vy3FJvEG8rOVU+hSArX1e3jyRTX\nXU8qi9YAh0XE4ylqUbc1OPT3gadGxPvKTltTypoeFK0LHyo7bu1CcWni68M4BwO5EHhGRBxZduz6\ne4rWin6fBxZFxDPK97VTRLyuyX1/GfjHiHhuFPaOonPfVcD6slNZT3menxkRz2vB+7mPrRPbSI/1\nDYpLHK8tn/ebQtH/4/cRMQNY2ExQNcn84+XnugfwD2z9GX6D4rt0TN0xvwS8KyIOKs/jjhHxvyJi\nSjPHZvDv1RTgD8AfI+LpwP9ucp/NGOpc1X9WXwdeFRHzys/psVF0yNwNtY3JXn+Rmb/NgW8Z+yeK\n69e/jKI3839R1MABTqPoqPMgReeyH44ijBHvKzN/TZF4v0FRk32Y4rr0QJ5LUfteT9Ff4ZhyH2Tm\nH4AjKa7HP0zRMep6/ufHQL//C5ycmX8sl08BXkxRe/leo/NZ1uTmUtwN8DuKywiHlav/FVhB8aNi\nNXB1g2MOW2Y+SNGRbjHFJYt9KK4X968/H/gE8M3y870eeHmT+/4W8HGK874eWAY8sUx8r6Roql5D\n8Zl+GdhptO+H4s6J/cre3MtGcawLKM7F7zLz2pryjwLPAX5P8UOpvrVrMO+l6ONyG8UdHN+g6JgH\nQGZeWa6fTtFXor98BfAO4AyK79ytFJeDmjLE9+ofKVov1lP8qDhvGO9nKEOdq1MofsCui4h/zMy7\nKDqEfpCio+ddFD8QzEdtFP/TaiVJkqrIX1KSJFWcyV6SpIoz2UuSVHEme0mSKs5kL0lSxVVqBL1d\ndtklZ82a1e0wJEnqiJUrVz6YmdOG2q5SyX7WrFmsWNGNmUUlSeq8iLhj6K1sxpckqfJM9pIkVZzJ\nXpKkijPZS5JUcSZ7SZIqzmQvSVLFmewlSao4k70kSRVnspckqeJM9pIkVZzJXpKkijPZS5JUcZWa\nCEdSa5y6/BZOv+Q325Sf8JJ9WDB33y5EJGk0IjO7HUPL9Pb2prPeSa1z1BeuAOC8dx7c5UgkNRIR\nKzOzd6jtbMaXJKniTPaSJFWcyV6SpIoz2UuSVHEme0mSKs5kL0lSxZnsJUmqOJO9JEkVZ7KXJKni\nTPaSJFWcyV6SpIoz2UuSVHEme0mSKs5kL0lSxZnsJUmqOJO9JEkVZ7KXJKniTPaSGlq2qo9Vd67j\nyjVrOWTxpSxb1dftkCSNkMle0jaWrepj0dLVbNy8BYC+dRtYtHS1CV8ap7bvdgCSxp4lF9/Mhkc3\nb1W24dHNLLn4ZubPmdH0fk5dfgunX/KbbcpPeMk+LJi776jjlNQck72kbdyzbsOwygeyYO6+LJi7\nL0d94QoAznvnwaOOTdLw2YwvaRvTp/YMq1zS2Na2ZB8RX42I+yPi+pqykyKiLyKuKR+vGOC1h0fE\nzRFxa0Sc2K4YJTW2cN5seiZP2qqsZ/IkFs6b3aWIJI1GO2v2ZwKHNyg/NTMPKB8X1a+MiEnAZ4CX\nA/sBR0fEfm2MU1Kd+XNmcMqR+7PDpOJPxIypPZxy5P7Dul4vaexo2zX7zLw8ImaN4KUHArdm5m0A\nEfFN4AjghtZFJ2ko8+fM4Nyr7gS81i6Nd924Zv+eiLiubObfucH6GcBdNct3l2WSJGkEOp3sPwc8\nDTgAuBf4j9HuMCKOj4gVEbHigQceGO3uJEmqnI4m+8y8LzM3Z+YW4EsUTfb1+oDda5Z3K8sG2ucX\nM7M3M3unTZvW2oAlSaqAjib7iHhqzeJrgOsbbPYrYJ+I2DMidgBeD1zQifgkSaqitnXQi4hzgUOB\nXSLibuAjwKERcQCQwO3AO8ttpwNfzsxXZOamiHgPcDEwCfhqZv66XXFKklR17eyNf3SD4q8MsO09\nwCtqli8CtrktT5IkDZ8j6EmSVHEme0mSKs5kL0lSxZnsJUmqOJO9JEkVZ7KX1FbLVvWx6s51XLlm\nLYcsvpRlqwYcI0tSm5jsJbXNslV9LFq6mo2btwDQt24Di5auNuFLHWayl9Q2Sy6+mQ2Pbt6qbMOj\nm1ly8c1dikiamEz2ktrmnnUbhlUuqT1M9pLaZvrUnmGVS2oPk72kbZy6/BZmnXghV65Zy5Vr1jLr\nxAuZdeKFnLr8lmHtZ+G82fRMnrRVWc/kSSycN7uV4UoaQmRmt2Nomd7e3lyxYkW3w5BUY9mqPj7w\n7evYuHkLM6b2sHDebObPmdHtsKRKiIiVmdk71HZtmwhHkgDmz5nBuVfdCcB57zy4y9FIE5PN+JIk\nVZzJXpKkijPZS5JUcSZ7SZIqzmQvSVLFmewlSao4k70kSRVnspckqeJM9pIkVZzJXpKkijPZS5JU\ncSZ7SZIqzmQvSVLFmewlSao4k70kSRVnspckqeK273YAkqrr1OW3cPolv/nL8qwTLwTghJfsw4K5\n+3YrLGnCiczsdgwt09vbmytWrOh2GJIkdURErMzM3qG2sxlfkqSKM9lLklRxJntJkirOZC9JUsWZ\n7CVJqri2JfuI+GpE3B8R19eULYmImyLiuog4PyKmDvDa2yNidURcExF2r5ckaRTaWbM/Ezi8rmw5\n8MzMfBZwC7BokNcflpkHNHNLgSRJGljbkn1mXg6srSv7UWZuKhd/CezWruNLkqRCN6/ZvxX4wQDr\nEvhRRKyMiOM7GJMkSZXTleFyI+KfgU3AOQNs8sLM7IuIJwPLI+KmsqWg0b6OB44HmDlzZlvilSRp\nPOt4zT4ijgNeCRyTA4zVm5l95b/3A+cDBw60v8z8Ymb2ZmbvtGnT2hCxJEnjW0eTfUQcDnwAeHVm\n/mmAbXaMiCn9z4GXAdc32laSJA2tnbfenQtcAcyOiLsj4m3AGcAUiqb5ayLi8+W20yPiovKluwI/\ni4hrgauACzPzh+2KU5KkqmvbNfvMPLpB8VcG2PYe4BXl89uAZ7crLkmSJhpH0JMkqeJM9pIkVZzJ\nXpKkijPZS5JUcSZ7SZIqrqne+BGxMzAd2ADcnplb2hqVJElqmQGTfUTsBLwbOBrYAXgAeCywa0T8\nEvhsZl7WkSglSdKIDVaz/zZwFvDXmbmudkVEPBd4U0TslZkN752XJEljw4DJPjPnDrJuJbCyLRFJ\nkqSWGrKDXhTeGBEfLpdnRsSAE9NIkqSxpZne+J8FDqa4dg+wHvhM2yKSJEkt1Uxv/IMy8zkRsQog\nMx+OiB3aHJckSWqRZmr2j0bEJCABImIa4K13kiSNE80k+08B5wNPjoiPAz8DTmlrVJIkqWWGbMbP\nzHMiYiXwEiCA+Zl5Y9sjkyRJLTFkso+IszPzTcBNDcokSdIY10wz/jNqF8rr989tTziSJKnVBkz2\nEbEoItYDz4qIP0TE+nL5fuC7HYtQkiSNyoDJPjNPycwpwJLMfEJmTikfT8rMRR2MUZIkjUIz99n/\nICJeVF+YmZe3IR5JktRizST7hTXPHwscSDEu/ovbEpEkSWqpZm69e1XtckTsDpzWtogkSVJLNdMb\nv97dwF+1OhBJktQezdxn/2nKoXIpfhwcAFzdzqAkSVLrNHPNfkXN803AuZn58zbFI0mSWqyZa/Zf\n60QgkiSpPZppxj8EOAnYo9w+gMzMvdobmiRJaoVmmvG/AiyguN1uc3vDkSRJrdZMsv99Zv6g7ZFI\nkqS2aCbZXxYRS4ClwH/3F2amPfIlSRoHmkn2B5X/9taUJY6gJ0nSuNBMb/zDOhGIJElqjwGTfUS8\nMTO/HhH/0Gh9Zn6yfWFJkqRWGaxmv2P575ROBCJJktpjwGSfmV8o//1o58KRJEmt1sygOnsC7wVm\n1W6fma9uX1iSJKlVmumNv4xiYJ3vAVuGs/OI+CrwSuD+zHxmWfZE4DyKHw+3A3+XmQ83eO2xwIfK\nxX912F5JkkammSlu/5yZn8rMyzLzJ/2PJvd/JnB4XdmJwCWZuQ9wSbm8lfIHwUcobvs7EPhIROzc\n5DElSVKNZpL96RHxkYg4OCKe0/9oZueZeTmwtq74CKC/lv41YH6Dl84Dlmfm2rLWv5xtfzRIkqQm\nNNOMvz/wJopBdPqb8UczqM6umXlv+fx3wK4NtpkB3FWzfHdZJkmShqmZZP86YK/M3Njqg2dmRkSO\nZh8RcTxwPMDMmTNbEpckSVXSTDP+9cDUFh7zvoh4KkD57/0NtukDdq9Z3q0s20ZmfjEzezOzd9q0\naS0MU5Kkamgm2U8FboqIiyPigv7HKI55AXBs+fxY4LsNtrkYeFlE7Fx2zHtZWSZJkoapmWb8j4x0\n5xFxLnAosEtE3F3uazHwnxHxNuAO4O/KbXuBd2Xm2zNzbUR8DPhVuauTM7O+o58kSWpCZI7qkvmY\n0tvbmytWrOh2GJIkdURErMzM3qG2a6ZmL0mqqFOX38Lpl/xmm/ITXrIPC+bu24WI1A7W7CVJHPWF\nKwA4750HdzkSDUezNftBO+hFxKSIOKd1YUmSpE4bNNln5mZgj4jYoUPxSJKkFmvmmv1twM/L2+0e\n6S/MzE+2LSpJktQyzST735aP7YAp7Q1HkiS12pDJPjM/ChARj8vMP7U/JEmS1EpDjqBXznZ3A3BT\nufzsiPhs2yOTJEkt0cxwuadRTDn7EEBmXgu8qJ1BSZKk1mkm2ZOZd9UVbW5DLJIkqQ2a6aB3V0S8\nAMiImAycANzY3rAkjYSjoUlqpJlk/y7gdGAGcA/F7HPvbmdQkkZmwdx9WTB3X0dDk7SVZnrjPwgc\n04FYJElSGzTTG3+viPheRDwQEfdHxHcjYq9OBCdJar9lq/pYdec6rlyzlkMWX8qyVX3dDkkt1kwH\nvW8A/wk8FZgOfAs4t51BSZI6Y9mqPhYtXc3GzVsA6Fu3gUVLV5vwK6aZZP+4zDw7MzeVj68Dj213\nYJKk9lty8c1seHTrG6w2PLqZJRff3KWI1A7NdND7QUScCHwTSOAo4KKIeCJAZq5tY3ySpDa6Z92G\nYZVrfGom2f9d+e8768pfT5H8vX4vSePU9Kk99DVI7NOn9nQhGrVLM73x9+xEIJKkzls4bzaLlq7e\nqim/Z/IkFs6b3cWo1GrN1OwlSRU1f84MAD7w7evYuHkLM6b2sHDe7L+UqxpM9pI0wc2fM4Nzr7oT\ncCCmqmpqbHxJkjR+NTOoziERsWP5/I0R8cmI2KP9oUkaCQdIkVSvmZr954A/RcSzgfcDvwXOamtU\nkkbEAVIkNdJMst+UmQkcAZyRmZ8BprQ3LEkj4QApkhpppoPe+ohYBLwReFFEbAdMbm9YkkbCAVIk\nNdJMzf4o4L+Bt2Xm74DdgCVtjUrSiAw0EIoDpEgT25DJPjN/l5mfzMyflst3ZqbX7KUxaOG82fRM\nnrRVmQOkSBqwGT8i1lMMh7vNKiAz8wlti0rSiDhAiqRGBkz2mWknPGkccoAUSfUGq9k/cbAXOtud\nJEnjw2C98VdSNONHg3XOdidJ0jgxWDO+s91JUsWduvwWTr/kN39ZnnXihQCc8JJ9WDB3326FpRZr\naiKciNgZ2Ad4bH9ZZl7erqAkSZ2xYO6+JvUJYMhkHxFvB06guL/+GuD5wBXAi9sbmiRJaoVmavYn\nAM8DfpmZh0XE04F/a29YkqSh1DfB97MJXvWaSfZ/zsw/RwQR8ZjMvCkiRjxCR/na82qK9gI+nJmn\n1WxzKPBdYE1ZtDQzTx7pMSWpivqb4I/6whWAt1pqYM0k+7sjYiqwDFgeEQ8Dd4z0gJl5M3AAQERM\nAvqA8xts+tPMfOVIjyNJkgpDJvvMfE359KSIuAzYCfhhi47/EuC3mTniHw+SJGlwA46NHxHPi4iX\n15Zl5k+ATcD+LTr+64FzB1h3cERcGxE/iIhnDBLn8RGxIiJWPPDAAy0KS5Kk6hhsIpxPADc0KP81\nLZj1LiJ2AF4NfKvB6quBPTLz2cCnKS4hNJSZX8zM3szsnTZt2mjDkiSpcgZrxp/SqHk9M++IiF1a\ncOyXA1dn5n0NjvGHmucXRcRnI2KXzHywBceVKssBUtQt3hkwtg2W7HceZN3jWnDsoxmgCT8ingLc\nl5kZEQdStEA81IJjSpXmACnqFu8MGNsGa8b/r4j4eET8ZWz8KJwMXDqag0bEjsBcYGlN2bsi4l3l\n4muB6yPiWuBTwOszs9F0u5IkaQiD1ezfD3wZuDUirinLng2sAN4+moNm5iPAk+rKPl/z/AzgjNEc\nQ5IkFQabCOcR4OiI2Avo7w3/68y8rSORSZKklhjs1rtZAJl5W2Z+r3zcVrM+ImK39ocoSRrIslV9\nrLpzHVeuWcshiy9l2aq+boekMWiwZvwlEbEdxbC1K4EHKGa92xs4jGJAnI8Ad7c7SEkTl728B7Zs\nVR+Llq5m4+YtAPSt28CipasBmD9nRjdD0xgzWDP+6yJiP+AY4K3AU4E/ATcCFwEfz8w/dyRKSROW\nvbwHtuTim9nw6OatyjY8upklF99sstdWBh0uNzNvAP65Q7FIkobhnnUbhlWuiWuwW+8kSWPY9Kk9\nwyrXxGWyl6RxauG82fRMnrRVWc/kSSycN+JZyFVRJntJGqfmz5nBKUfuzw6Tij/lM6b2cMqR+3ft\ner13BoxdQ05xW46gdwywV2aeHBEzgadk5lVtj06SNKj5c2Zw7lV3At3tvOidAWNbMzX7zwIHU4xl\nD7Ae+EzbIpIkjTuD3Rmg7huyZg8clJnPiYhVAJn5cDk9rSRJgHcGjHXN1OwfjYhJQAJExDRgS1uj\nkiSNK94ZMLY1k+w/BZwPPDkiPg78DDilrVFJUg07fo193hkwtg3ZjJ+Z50TESorhcQOYn5k3tj0y\naYJxWNjG7Pg1PvR/Fh/49nVs3LyFGVN7WDhvtp/RGNFMb/yzM/NNwE0NyiS1iMPCNuaQsOPHWLkz\nQNtqphn/GbUL5fX757YnHEnamh2/pNEbbIrbRRGxHnhWRPwhItaXy/dTzIQnSW1nxy9p9AZM9pl5\nSmZOAZZk5hMyc0r5eFJmLupgjJImMDt+SaPXTAe9RRGxM7APxXz2/eWXtzMwSQI7fg2mvlPnrBMv\nBOzUqW0100Hv7cAJwG7ANcDzgSuAF7c3NEkq2PGrsf5OndJQmumgdwLwPOCOzDwMmAOsa2tUkiSp\nZZpJ9n/OzD8DRMRjMvMmwItlkiSNE82MjX93REwFlgHLI+Jh4I72hiVJklqlmQ56rymfnhQRlwE7\nAT9sa1TSBNU/LOzGzVs4ZPGldkST1BKDJvtyAJ1fZ+bTATLzJx2JSpqAHBZW45l3Boxtgyb7zNwc\nETdHxMzMvLNTQUkTkcPCajzzzoCxrZlr9jsDv46Iq4BH+gsz89Vti0qagBwWVlK7NJPs/6XtUUhi\n+tQe+hokdoeFlTRazXTQ8zq91AEL581m0dLVWzXlOyyspFZopmYvqQMcFrYxO35JoxeZ2e0YWqa3\ntzdXrFjR7TCkUXE+e0nNioiVmdk71HbNjKAnSZLGsQGb8SNiNdCo2h9AZuaz2haVJElqmcGu2b+y\nY1FIkqS2GTDZZ6bj30uSVAGDNeOvZ/Bm/Ce0LSpJktQyg9Xsp7TzwBFxO7Ae2Axsqu9NGBEBnA68\nAvgTcFxmXt3OmCRJqqKm77OPiCcDj+1fbtFY+Ydl5oMDrHs5sE/5OAj4XPmvJEkahiFvvYuIV0fE\nb4A1wE+A24EftDkugCOAs7LwS2BqRDy1A8eVJKlSmqnZfwx4PvBfmTknIg4D3tiCYyfwo4hI4AuZ\n+cW69TOAu2qW7y7L7q3dKCKOB44HmDlzZgvCkrrDkeImjvrPup+ftdplyBH0ImJFZvZGxLXAnMzc\nEhHXZuazR3XgiBmZ2VdeHlgOvDczL69Z/31gcWb+rFy+BPinzBxwiDxH0JM0njhaokar2RH0mqnZ\nr4uIxwOXA+dExP3UTHU7UpnZV/57f0ScDxxYHqNfH7B7zfJuZZkkSRqGZpL9EcAGYAFwDLATcPJo\nDhoROwLbZeb68vnLGuzzAuA9EfFNio55v8/Me5HGGJtkJY11g91nvzewa2b+vCzaAnwtIl4ITAUe\nGsVxdwXOL+6uY3vgG5n5w4jfnFNDAAAO/UlEQVR4F0Bmfh64iOK2u1spbr17yyiOJ7XNgrn7smDu\nvjbJShqzBqvZnwYsalD++3Ldq0Z60My8Ddjmmn+Z5PufJ/DukR5DkiQVBrv1btfMXF1fWJbNaltE\nkiSppQZL9lMHWdfT6kAkSVJ7DJbsV0TEO+oLI+LtwMr2hSRJ1bdsVR+r7lzHlWvWcsjiS1m2ypuN\n1D6DXbN/H0UnumP4n+TeC+wAvKbdgUlSVS1b1ceipavZuHkLAH3rNrBoaXHVdP6cGd0MTRU1YM0+\nM+/LzBcAH6UYIvd24KOZeXBm/q4z4UlS9Sy5+GY2PLp5q7INj25mycU3dykiVd2Q99ln5mXAZR2I\nRRq3+ptkN27ewiGLL2XhvNnW0DSge9ZtGFa5NFpDToQjaXADNcl6DVYDmT61cR/ngcql0TLZS6Nk\nk6yGa+G82fRMnrRVWc/kSSycN7tLEanqmp7PXlJjNslquPov8Xzg29excfMWZkzt8dKP2spkL43S\n9Kk99DVI7DbJajDz58zg3KvuBBxiWe1nM740SjbJShrrrNlLo2STrKSxzmQvtYBNshOHUxprPDLZ\nS9IwOKWxxiOv2UuSVHEme0mSKs5kL0lSxXnNXpI6rL6T36wTLwTs5Kf2MdlLUof1d/KTOsVkL42S\ntbTxoZW3zDnLocabyMxux9Ayvb29uWLFim6HIWkMG+0tc/2zHNZOftQzeRKnHLm/CV8dFxErM7N3\nqO3soCdJw+AshxqPTPaSNAzOcqjxyGv2mrAc9lQj4SyHGo9M9pqwHPZUI7Fw3uyG1+yd5VBjmcle\nkobBWQ41HpnsJU0YrbplzlkONd7YQU/ShNB/y9zGzVsA6Fu3gUVLV7NsVV+XI5Paz2QvaULwljlN\nZCZ7TWj9zbpXrlnLIYsvtZZXYd4yp4nMZK8Jy2bdiWWgW+O8ZU4TgcleE5bNuhPLwnmz6Zk8aasy\nb5nTRGFvfE1YNutOLK26Zc6JjzQemew1YTkS2sTTilvmnJ5W45HN+JqwbNaVNFFYs9eE5Uho0tjj\nnBXt0fFkHxG7A2cBuwIJfDEzT6/b5lDgu8CasmhpZp7cyTg1MTgSmjS2OGdFe3SjZr8JeH9mXh0R\nU4CVEbE8M2+o2+6nmfnKLsQnSVKldDzZZ+a9wL3l8/URcSMwA6hP9pLUMvai10TW1Wv2ETELmANc\n2WD1wRFxLXAP8I+Z+esB9nE8cDzAzJkz2xOopHHPXvSayLrWGz8iHg98B3hfZv6hbvXVwB6Z+Wzg\n08CygfaTmV/MzN7M7J02bVr7ApYkdYTDWLdeV5J9REymSPTnZObS+vWZ+YfM/GP5/CJgckTs0uEw\nJUkd5jDW7dGN3vgBfAW4MTM/OcA2TwHuy8yMiAMpfpQ81MEwNQF4DVcaewYbxtrbYkeuG9fsDwHe\nBKyOiGvKsg8CMwEy8/PAa4H/HRGbgA3A6zMzuxCrKsxruNLY4zDW7dGN3vg/A2KIbc4AzuhMRJKk\nscJhrNvD4XIlSWOGw1i3h8PlSpLGDIexbg+TvSRpTHEY69Yz2WtccrIMSWqeyV7jkpNlSFLz7KAn\nSVLFmewlSao4m/E1bvWPn71x8xYOWXypPXalCnBky/aIKg1M19vbmytWrOh2GOqA/vGza4fV7Jk8\niVOO3N+EL2nCiIiVmdk71HY242tcGmz8bEnS1kz2GpccP1uSmmey17g00DjZjp8tSdsy2Wtccvxs\nSWqevfE1Ljl+tiQ1z2SvccvxsyWpOTbjS5JUcSZ7SZIqzmZ8dVSrZqtzlC1Jap4j6KkrnK1Okkav\n2RH0rNlLkiqlVS2IVWLNXh23bFWft8xJaruJ0ILo2Pgak/onsNm4eQsAfes2sGjpapat6utyZJKq\npH9WzCvXrOWQxZdO+L8xJnt1lBPYSGo3KxXbMtmro5zARlK7WanYlsleHeUENpLazUrFtkz26ign\nsJHUblYqtmWyV0fNnzODU47cnx0mFV+9GVN7OOXI/e2NL6llrFRsy/vs1ZRW3rfqBDaS2slZMbfl\nffZqWivuj3ewC0md4n32NduZ7NWM/ltZanu49kyeZBO8pDFnIlUqTPZqqUMWX0pfg56sM6b28PMT\nX9yFiCSpvcbDjwbHxldLeSuLpIlmz112pGfypG1aNPfcZcdh7Wcs/Ggw2U8ArfiiTZ/a07BmP5Fv\nZZFUbYMNzjOcy5et+tEwGt56NwH0f9FqDfeL5q0skiaaVrVojoUR/azZN7BsVR9LLr6Ze9ZtYHqX\nb9loRSyt+HXqrSySJppWtWiOhcug1uzr9Pc671u3gWR0EygsW9XHIYsvZc8TLxzRrEutiqUVX7RT\nl9/C+867ZquJJd533jWcuvyWYcUiSeNFq1o0x8KIfl2p2UfE4cDpwCTgy5m5uG79Y4CzgOcCDwFH\nZebtnYjtX757fcNa8L989/ph1WLrb1XrT9RA0/tpVSyt+HW6YO6+Y6b3qSR1wpoHH2n4N3jNg48M\naz8L581ueOtyJy+DdjzZR8Qk4DPAXOBu4FcRcUFm3lCz2duAhzNz74h4PfAJ4KhOxPfHP28aVvlA\nWpGoWxXLWPiiSdJ406pKTqt+NIxGN2r2BwK3ZuZtABHxTeAIoDbZHwGcVD7/NnBGRER2YFCAVl2j\naUWifvxjt2d9g+0f/9jhfWxj4YsmSRPVWGgZ7UaynwHcVbN8N3DQQNtk5qaI+D3wJODB+p1FxPHA\n8QAzZ84cdXDPnP6Ehsn+mdOfMKz9tCJRf+yIZzaskX/siGcOK5ax8EWTJHXPuO+gl5lfzMzezOyd\nNm3aqPf3hTf3ctpRBzBjag9BMULcaUcdwBfePOQARVv52BHPbNixYziJun+GuNpYHJ5WkjRc3ajZ\n9wG71yzvVpY12ubuiNge2Imio15HzJ8zY9QJtf/1o71trhWxSJImtm4k+18B+0TEnhRJ/fXAG+q2\nuQA4FrgCeC1waSeu17eaiVqSNBZ0PNmX1+DfA1xMcevdVzPz1xFxMrAiMy8AvgKcHRG3AmspfhBI\nkqQR6Mp99pl5EXBRXdmHa57/GXhdp+OSJKmKxn0HPUmSNDiTvSRJFWeylySp4kz2kiRVnMlekqSK\nM9lLklRxJntJkirOZC9JUsWZ7CVJqrgYh0PODygiHgDuaOEud6HBtLpqGc9ve3l+289z3F6e36Ht\nkZlDTvlaqWTfahGxIjOHN7etmub5bS/Pb/t5jtvL89s6NuNLklRxJntJkirOZD+4L3Y7gIrz/LaX\n57f9PMft5fltEa/ZS5JUcdbsJUmqOJN9AxFxeETcHBG3RsSJ3Y6niiLi9ohYHRHXRMSKbscz3kXE\nVyPi/oi4vqbsiRGxPCJ+U/67czdjHM8GOL8nRURf+R2+JiJe0c0Yx7OI2D0iLouIGyLi1xFxQlnu\nd7hFTPZ1ImIS8Bng5cB+wNERsV93o6qswzLzAG+taYkzgcPryk4ELsnMfYBLymWNzJlse34BTi2/\nwwdk5kUdjqlKNgHvz8z9gOcD7y7/7vodbhGT/bYOBG7NzNsycyPwTeCILsckDSozLwfW1hUfAXyt\nfP41YH5Hg6qQAc6vWiQz783Mq8vn64EbgRn4HW4Zk/22ZgB31SzfXZaptRL4UUSsjIjjux1MRe2a\nmfeWz38H7NrNYCrqPRFxXdnMbxNzC0TELGAOcCV+h1vGZK9ueWFmPoficsm7I+JF3Q6oyrK47cZb\nb1rrc8DTgAOAe4H/6G44419EPB74DvC+zPxD7Tq/w6Njst9WH7B7zfJuZZlaKDP7yn/vB86nuHyi\n1rovIp4KUP57f5fjqZTMvC8zN2fmFuBL+B0elYiYTJHoz8nMpWWx3+EWMdlv61fAPhGxZ0TsALwe\nuKDLMVVKROwYEVP6nwMvA64f/FUagQuAY8vnxwLf7WIsldOfhEqvwe/wiEVEAF8BbszMT9as8jvc\nIg6q00B5C81pwCTgq5n58S6HVCkRsRdFbR5ge+AbnuPRiYhzgUMpZgm7D/gIsAz4T2AmxWyQf5eZ\ndjIbgQHO76EUTfgJ3A68s+b6soYhIl4I/BRYDWwpiz9Icd3e73ALmOwlSao4m/ElSao4k70kSRVn\nspckqeJM9pIkVZzJXpKkijPZS+NIRPyxbvm4iDijwzG8LiJujIjL6spnRcQb2hFbuW/vY5dGyGQv\niYjYfhibvw14R2YeVlc+C3jDtptL6jaTvVQRZe330nJilksiYmZZfmZEvLZmuz+W/x4aET+NiAuA\nGxrs7+iIWB0R10fEJ8qyDwMvBL4SEUvqXrIY+OtybvcFZdn0iPhhOR/5v9fs+2URcUVEXB0R3yrH\nRK8//nMj4tqIuBZ4d937/Gn52qsj4gVl+VkRMb9mu3MiwhkrJUz20njTUybTayLiGuDkmnWfBr6W\nmc8CzgE+1cT+ngOckJn71hZGxHTgE8CLKUaJe15EzM/Mk4EVwDGZubBuXycCPy3ndj+1LDsAOArY\nHzgqInaPiF2ADwEvLSdDWgH8Q4PY/h/w3sx8dl35/cDc8rVH1bzPrwDHlfHvBLwAuLCJcyBV3nCa\n7iR134bMPKB/ISKOA3rLxYOBI8vnZwP/ztCuysw1DcqfB/w4Mx8oj3MO8CKKIXiH45LM/H25jxuA\nPYCpwH7Az4sh0dkBuKL2RRExFZhaziPf/35eXj6fDJwREQcAm4F9ATLzJxHx2YiYBvwt8J3M3DTM\neKVKMtlL1beJshUvIrajSK79Hmnzsf+75vlmir85ASzPzKNHuM8FFOPTP5viff25Zt1ZwBspJrB6\nywj3L1WOzfhSdfyCIskBHEMxsQgUk7Q8t3z+aoqa8VCuAv4mInaJiEnA0cBPhnjNemBKE/v+JXBI\nROwNf5kFcavLCJm5DlhXTpACxfvptxNwbzm17JsoJqzqdybwvnIf2/RDkCYqk71UHe8F3hIR11Ek\nwRPK8i9RJO5rKZr6h6zNl7O3nQhcBlwLrMzMoaYXvQ7YXHaqWzDQRuWlgeOAc8tYrwCe3mDTtwCf\nKfsmRE35Z4Fjy/fz9Nr3k5n3ATdSXO+XVHLWO0mVERGPo5gm9Tn9fQUkWbOXVBER8VKKWv2nTfTS\n1qzZS5JUcdbsJUmqOJO9JEkVZ7KXJKniTPaSJFWcyV6SpIoz2UuSVHH/H7eBYGHaEHKJAAAAAElF\nTkSuQmCC\n",
            "text/plain": [
              "<Figure size 576x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Hlc0hUWGPaoK",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ghfiZ1MJsTt0",
        "colab_type": "text"
      },
      "source": [
        "## 3. Summarized Findings"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DwT7aHoVse6h",
        "colab_type": "text"
      },
      "source": [
        "The number of calls per minute in each hour increases exponentially from midnight and peaks at 11 AM. The calls rates decrease until 2 PM and start increases again, peaks at 6 PM. After 6 PM, the calls rates decrease exponentially. \n",
        "\n",
        "\n",
        "Starting from 11 AM to noon, the average calls rate per minute is 15.4 calls, with 98% in $\\pm$ 1.2 calls interval. At noon, the average calls rate per minute is also high, at 14.3 calls, with 98% in $\\pm$ 1.1 calls interval\n",
        "Another local mode is at 6 PM, where the average calls rate is 8.8 calls, with 98% in $\\pm$ 0.9 calls interval. Based on the data, 92.2% of the number of calls happened between 8 AM to 8 PM. \n",
        "\n",
        "\n",
        "Business action items: The company should arrange more employees at noon (from 11 AM till 1 PM), with the average 12.8 calls, and at 6 PM, with an average of 8.8 calls. Busy hours are between 8 AM and 8 PM, which might require additional staffs than other periods. More data needs to be collected (call duration, etc.) if we want to have an accurate estimate of how much time each employee respond to a phone call.  "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lbxQ4SURPaoM",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tLon4mKM0cG3",
        "colab_type": "text"
      },
      "source": [
        "$\\pm$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vSkpV9O1shff",
        "colab_type": "text"
      },
      "source": [
        "# Section 2: Stretch goal\n",
        "Explained the **compute_posterior** function."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WK01DTIvtG0v",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def compute_posterior(parameter_values, data):\n",
        "    log_prior = np.log(prior(parameter_values))\n",
        "    log_likelihood = np.array([\n",
        "                               np.sum(np.log(likelihood(param, data))) \n",
        "                               for param in parameter_values])\n",
        "    unnormalized_log_posterior = log_prior + log_likelihood\n",
        "    unnormalized_log_posterior -= max(unnormalized_log_posterior)\n",
        "    unnormalized_posterior = np.exp(unnormalized_log_posterior)\n",
        "    area = sp.integrate.trapz(unnormalized_posterior, parameter_values)\n",
        "    posterior = unnormalized_posterior / area\n",
        "    return posterior"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oBu2rNkis07D",
        "colab_type": "text"
      },
      "source": [
        "## 1. How is the function multiply the prior and likelihood using logarithms"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Qgbc5_OyuPp3",
        "colab_type": "text"
      },
      "source": [
        "As we know: $$posterior = prior * likelihood$$\n",
        "\n",
        "We have: \n",
        "1. prior is $Gamma(λ | α, β)$\n",
        "2. likelihood is $P(datum | λ)$\n",
        "\n",
        "If given more data (assuming all the data points are independent: which is a valid claim because a call from 1 customer will not likely to affect other customers' calls).\n",
        "\n",
        "Then, with $n$ data points\n",
        "$${posterior = prior * P(datum_1 | λ) * P(datum_2 | λ) * ... * P(datum_n | λ)}$$\n",
        "$$$$\n",
        "$$\\Rightarrow \\ln(posterior) = \\ln(prior * P(datum_1 | λ) * P(datum_2 | λ) * ... * P(datum_n | λ))$$\n",
        "$$$$\n",
        "$$\\Rightarrow \\ln(posterior) = \\ln(prior) + \\ln(P(datum_1 | λ)) + \\ln(P(datum_2 | λ)) + ... + \\ln(P(datum_n | λ))$$\n",
        "$$$$\n",
        "$$\\Rightarrow \\ln(posterior) = \\ln(prior) + \\sum_{i=1}^{n} \\ln(P(datum_i | λ))$$\n",
        "\n",
        "\n",
        "Hence, the log of the posterior equals to sum of the log of the prior and the series of likelihood. This is demonstrated in the code where **unnormalized_log_posterior = log_prior + log_likelihood**. \n",
        "\n",
        "Where the **log_likelihood** variable is already a sum of all log of each datum likelihood by using **np.sum** above. \n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "VYXVzf6nsghJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def compute_posterior(parameter_values, data):\n",
        "    # take the log of the prior\n",
        "    log_prior = np.log(prior(parameter_values))\n",
        "    \n",
        "    # for each lambda (parameter_values), we compute the sum of all of log of\n",
        "    # each datum likelihood by using np.sum\n",
        "    log_likelihood = np.array([np.sum(np.log(likelihood(param, data))) for param in parameter_values])\n",
        "    \n",
        "    # the log of posterior is the sum of log of prior and log of likelihood\n",
        "    unnormalized_log_posterior = log_prior + log_likelihood\n",
        "    unnormalized_log_posterior -= max(unnormalized_log_posterior)\n",
        "    unnormalized_posterior = np.exp(unnormalized_log_posterior)\n",
        "    area = sp.integrate.trapz(unnormalized_posterior, parameter_values)\n",
        "    posterior = unnormalized_posterior / area\n",
        "    return posterior"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "O4666lmlyt-p",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6MZ8maiztBem",
        "colab_type": "text"
      },
      "source": [
        "## 2. The purpose of np.sum() in line 4"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GqrEO_OgyvAd",
        "colab_type": "text"
      },
      "source": [
        "The np.sum takes the sum of all values within an array. In line 4, it is computing $\\sum_{i=1}^{n} \\ln(P(datum_i | λ))$ with $n$ is the number of data points. For each lambda, we compute the sum of the likelihood of each datum given that lambda value."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "N3B_JBcGPaoP",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jNP90Z-BtQ-k",
        "colab_type": "text"
      },
      "source": [
        "## 3. The purpose of the maximum of the unnormalized log posterior is subtracted in line 7"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MZU-oG7W03nt",
        "colab_type": "text"
      },
      "source": [
        "If we don't subtract the max, in some real-life cases, the maximum will be too big, which makes line 8 occurs Warning like **overflow encountered in exp** \n",
        "\n",
        "(as Python cannot fully represents np.exp(x) as x is too big and occur numerical overflow - if x > 710, this warning will occur).\n",
        "\n",
        "Subtracting will make the largest value in the posterior to be 0 and others values got smaller -> avoid numerical overflow for extremely big numbers. \n",
        "\n",
        "The act of subtracting the max in line 7 does not affect the visualization/intervals at all as we are just simply rescaling an additional step by dividing by a constant factor (subtract for log functions is dividing for after ward functions).\n",
        "\n",
        "This is a pratical step to avoid Python numerical overflow error.\n",
        "\n",
        "Example: Given the same dataset, we modify the code to see the warning. The warning occurs, which makes some values will be represented as **inf**. Then, these values are unprocessable, which lead to another error **invalid value encountered in true_divide**. This will create missing values at some intervals.  \n",
        "\n",
        "As we can see from the code below, at hour[11], the value got too big, lead to overflow. Then, if we look at the graph, the values for hour = 11 does not exist. This cause the figure to be incorrect. Hence, subtracting the max is the practical move. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4aNP2Q7WTjwM",
        "colab_type": "code",
        "outputId": "6da32921-4b86-4469-e49c-7a2716021b16",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 922
        }
      },
      "source": [
        "def compute_posterior(parameter_values, data):\n",
        "    log_prior = np.log(prior(parameter_values))\n",
        "    log_likelihood = np.array([\n",
        "                               np.sum(np.log(likelihood(param, data))) \n",
        "                               for param in parameter_values])\n",
        "    unnormalized_log_posterior = log_prior + log_likelihood\n",
        "    #unnormalized_log_posterior -= max(unnormalized_log_posterior)\n",
        "    unnormalized_posterior = np.exp(unnormalized_log_posterior)\n",
        "    area = sp.integrate.trapz(unnormalized_posterior, parameter_values)\n",
        "    posterior = unnormalized_posterior / area\n",
        "    return posterior\n",
        "\n",
        "# Create empty lists for storing the data\n",
        "#\n",
        "# upper_percentiles_diff to store the difference between the 99% \n",
        "# and the posterior mean\n",
        "#\n",
        "# lower_percentiles_diff to store the difference between the 01% \n",
        "# and the posterior mean\n",
        "#\n",
        "# lambda_ to store the posterior mean\n",
        "upper_percentiles_diff = []\n",
        "lower_percentiles_diff = []\n",
        "lambda_ = []\n",
        "\n",
        "for i, waiting_times_hour in  enumerate(waiting_times_per_hour):\n",
        "    print(i)\n",
        "    # initialize the lambdas values\n",
        "    lambdas = np.linspace(0, 20, 501)[1:]\n",
        "    #\n",
        "    # compute the posterior by using the compute_posterior function\n",
        "    # with the initialized lambdas and the data (of each hour)\n",
        "    posterior = compute_posterior(lambdas, waiting_times_hour)\n",
        "    #\n",
        "    # Add the newly calculated posterior mean to the lambda_ list\n",
        "    lambda_.append(sp.integrate.trapz(posterior * lambdas, lambdas))\n",
        "    #\n",
        "    # Calculate the difference between the posterior mean with the 99% percentile\n",
        "    # and append to the upper_percentiles_diff list\n",
        "    upper_percentiles_diff.append(compute_percentile(lambdas, posterior, 0.99) - lambda_[-1])\n",
        "    #\n",
        "    # Calculate the difference between the posterior mean with the 1% percentile\n",
        "    # and append to the upper_percentiles_diff list\n",
        "    lower_percentiles_diff.append(lambda_[-1] - compute_percentile(lambdas, posterior, 0.01))\n",
        "\n",
        "# initialize the graph\n",
        "plt.figure(figsize=(8, 6))\n",
        "# plot the error bar with the mean at lambda_ and the range from the mean derives\n",
        "# from the upper_percentiles_diff and lower_percentiles_diff\n",
        "plt.errorbar(list(range(24)), lambda_, yerr = [upper_percentiles_diff,lower_percentiles_diff]\n",
        "             , fmt=\"o\", capsize=3)\n",
        "# label the axes and plot\n",
        "plt.ylabel('Call rate (Calls per minute)')\n",
        "plt.xlabel('Hour of the day')\n",
        "plt.title('Mean and 98% confidence interval over call rate')\n",
        "plt.show()"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "0\n",
            "1\n",
            "2\n",
            "3\n",
            "4\n",
            "5\n",
            "6\n",
            "7\n",
            "8\n",
            "9\n",
            "10\n",
            "11\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:8: RuntimeWarning: overflow encountered in exp\n",
            "  \n",
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:10: RuntimeWarning: invalid value encountered in true_divide\n",
            "  # Remove the CWD from sys.path while we load stuff.\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "stream",
          "text": [
            "12\n",
            "13\n",
            "14\n",
            "15\n",
            "16\n",
            "17\n",
            "18\n",
            "19\n",
            "20\n",
            "21\n",
            "22\n",
            "23\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAAGDCAYAAAA72Cm3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucXGV9+PHPlxh00cCioDUrIVBI\nWhEluooYa1VM492UnxYR7xfw91OLqQ1NrL96af2RNq1AvRWr1huiFWPEokZK8I5gwiJBMKAEAgty\nEVcjxiZsvr8/zlmcXfYyuzuzM2f283695pU5zzlzznfOTPY753me8zyRmUiSpOrZp9UBSJKkqTGJ\nS5JUUSZxSZIqyiQuSVJFmcQlSaook7gkSRVlElclRcSNEfGsVscxEyKiKyK+EhG/iogvRMTJEfGN\ncbb/ZkS8fiZjHE9E/FtE/N9Wx1GPTvteRURGxBHl809ExD+0OiY1lkm8Q5V/jHZHxEEjyvvK/9gL\nWxNZ80VEd0R8MiLuKB/vGrH+mIj4TpkUb6lNMBFxSET8ICLujoh/GfG6r0VE7wy9jVovBh4BPCwz\nX5KZ52bmn7UgjinJzDdm5t/Xs62Jpjo67QdPVZnEO9t24KShhYg4GtivdeHMmDMp3udC4EnAKyLi\nNTXrPwt8G3go8KfA/4mIF5br1gCfBA4DVgwl7Yg4EdiemZtn5B0MdyhwXWbe24JjV0pEPKDVMUxH\nu8TfLnFoYibxzvZp4JU1y68CPlW7QUQ8MCL+OSJ2RMTtZdVnV7nuwIj4r4i4MyJ+WT5/VM1rvxkR\nfx8R34uInRHxjZFX/jXbTmtfEfGKiLgpIn4REX87wft+AfBPmfnbzLwR+Bjw2pr1C4FzM3MwM38G\nfBc4qlx3GLApM38F/BA4PCL2B1YDb5/guETEURFxUXklf3tEvL3mPJ8VEbeWj7Mi4oHluqeXNQJv\nK2sObhv60RER7wb+DjgxIn4TEa+LiFdHxHdrjrksIn5S1ix8AIgRMb02Iq4tz/vGiDi0Zl1GxBsj\n4vqIGIiID0ZE1Kx/Q/nanRFxTUQ8viyfHxFfLD/P7RHxl+Ock/uurid4r6cAJwOnl+/1KxMdKyLe\nFRHnR8RnIuLXwNsjYldEPLRmmyURcVdEzI2IP4yITeX36K6IODciuif6XMv9HBARnyrjuCki3hER\n+5Sf7UBEPKZm24PLOB5eLj8/Iq4st/t+RDy2ZtsbI+JvIuIq4J4YJYGO8716UkRcWu73toj4QETs\nW8/7GbH/V5f/986MiF8A7xrvXEXEp4EFwFfKz+r0svzJ5fsbiIgfRcTTJxuLJikzfXTgA7gReBaw\nDfhjYA5wC8VVXQILy+3OBC6guCqdB3wFOKNc9zDgf1Fc1c4DvgBsqDnGN4GfAYuArnJ57RjxTHlf\nwKOB3wBPAx4IvA+4F3jWGMe6C3hSzfLfAr+sWf5/wFpgLrC4PC9PLNetA94MdAPXUyT3s4FX1XHO\n5wG3AW8DHlQuH1uuew/wA+DhwMHA94G/L9c9vXw/7yljei7wW+DAcv27gM/UHOfVwHfL5wcBOymq\n3OcCK8t9vb5c/yLgp+V34AHAO4Dv1+wrgf8q3+8C4E7g2eW6lwD9wBMpfhgcQfH92QfYQvHjYl/g\ncOAGYPkY5+UTwD/U+V7v27ZcHvdY5bnZA6wot+0CNgFvqNnHOuDfyudHAMsovkcHU9TInDXy/80Y\n7+NTwJfLz3UhcB3wunLdx4H31mz7JuDr5fMlwB3AsRT/D19VHueBNce8EjgE6Jrk9+oJwJPLz3Yh\ncC3w1hGf7xGjndsRx3h1+bm8pdxX12TPFdAD/KL8TPcpX/sL4OBW/z3s5EfLA/DRpA/290n8HcAZ\nwLOBi8r/oFn+hw/gHuAPa153HEW18Wj7PIbhyfCbwDtqlv/P0B+uOuKre18Uf8A/V7PuwcDucf7Y\nfgZYX/6xO4Lix8H/1Kx/CkViu7c8F++uWfdQ4PPAjygS4hLgkrJ8qBr+zWMc9ySgb4x1PwOeW7O8\nHLixfP50YBfwgJr1dwBPLp+/i7GT+CuBH9SsC4ofJUNJ/GuUiaZc3ociaR5aLifw1Jr1/wmsLp9v\nBE4b5b0cC+wYUbYG+I8x3vsnGJ7Ex3uv921bz7HKc/PtEetfT1GbMnQ+bgaeNkZsK2o/M8ZI4hTJ\ndzfw6JqyU4Fvls+fBfysZt33gFeWzz9M+YOtZv024E9rjvnacf6vjPm9GmXbtwJfqlmeTBLfMcG+\nxz1XwN8Anx7xmo3U8QPYx9Qftnt0vk9TJJ7DGFGVTvHrej9gS20NKsUfLCJiP4or9WcDB5br50XE\nnMwcLJd/XrO/3wIPGS2Iae5rPsUfYgAy856yym8sfwm8n+JK+hfAeZR9A8pq1q9TXG1/FvgD4PyI\nuD0zP5SZdwMnltvuQ3Hu3khRnX41xR+7KyLi4sy8dsRxD6FI1qOZD9xUs3xTWTbkFzm8zXvMcznK\nfmvPTUbEzTXrDwXOjuGd9ILiqmkonrHO+1jv51BgfkQM1JTNAb5TR7wwufdaz7FuHv4Svgi8PyIe\nSVGzs3do+4h4BEXNyp9Q/MjbB/hlHTEfRFFzMPIz7CmfXwLsFxHHArdT/Ej9Us17eFVEvKXmtfsy\n/PMf+R5qjfm9iohFFDVTvRT/lx9AUXMxFcNimMK5OhR4SUS8oKZsLsW5UZPYJt7hMvMmig5uz6W4\nOq11F8VV0VGZ2V0+DsjMoT+ob6Oobj42M/enqM6GEW2udZrOvm6j+ENWvKD4QfCwsTbOzLsz8+TM\n/IPMPIrie355ufpwYDAzP5WZ92bmLcDnKM7PSKdQXOVeDRwNbM7M3cDWcnmkm8v9j+ZWij9yQxaU\nZdM18txE7XIZ06k1n293ZnZl5vfr2PfNwB+OUb59xD7nZeZo53CyRk6rWM+xhr0mM38JfIPix9jL\nKGpxhrb5f+X2R5ffw5dT33fwLopq+5GfYX95zEGKWoyTysd/ZebOmvfw3hHvYb/MPG+c9z3yHIz1\nvfow8BPgyPL9vL3O9zOakTFMdK5G+6w+PeJ9Pjgz104xHtXBJD47vA54ZmbeU1uYmXuBfwfOrOmA\n0xMRy8tN5lEk+YHyCvad04hhOvs6H3h+RDy17LTzHsb57pYdch4WEXMi4jkUyXjotqXrik3iZWWn\npD+g+GN/1Yh9PJyiXfNdZdF24BkR8RCKq54bRjn0fwGPjIi3lp2d5pVXZlDUBryj7PB0EEUTwWcm\ncQ7GciFwVEScUHaI+kuK2oUh/wasiYijyvd1QES8pM59fxT464h4QhSOiKJT3OXAzrIzVld5nh8T\nEU9swPu5neEJa6rH+ixFU8OLy+dD5lH0r/hVRPQAq+oJqiZJv7f8XA8F/orhn+FnKb5LJ4845r8D\nb4yIY8vz+OCIeF5EzKvn2Iz/vZoH/Br4TUT8EfC/69xnPSY6VyM/q88AL4iI5eXn9KAoOjI+CjWN\nSXwWyMyf5di3Rv0NRfvwD6Lo3fvfFFfMAGdRdHC5i6JT1tenEcaU95WZP6ZIqJ+luPL8JUW771ie\nQHG1vJOiP8DJ5T7IzF8DJ1C0d/+SokPR1fw+yQ/5Z+A9mfmbcvkM4JkUVxtfGe18lldeyyh6x/+c\nojr/GeXqfwA2U/xY2ApcMcoxJy0z76LogLaWoungSIr22KH1XwL+Efhc+fleDTynzn1/AXgvxXnf\nCWwAHlomtOdTVBlvp/hMPwocMN33Q3EnwaPL3s0bpnGsCyjOxc8z80c15e8GHg/8iuIH0MjaqfG8\nhaIPyQ0UdzR8lqJDGwCZeVm5fj5FX4Sh8s3AG4APUHznfkrRLFOXCb5Xf01R27CT4sfC5yfxfiYy\n0bk6g+KH6UBE/HVm3kzRkfLtFB0kb6ZI/OaZJorf1zJJkqQq8ReSJEkVZRKXJKmiTOKSJFWUSVyS\npIoyiUuSVFGVGLHtoIMOyoULF7Y6DEmSZsSWLVvuysyDJ9quaUk8Ij5OcX/nHZn5mLJsaFzqhRTj\n7v5FObrSuBYuXMjmza2YAVKSpJkXETdNvFVzq9M/QTFOdq3VwMWZeSRwcbksSZKmoGlJPDO/Ddw9\novhFwCfL55+kmBVHkiRNwUx3bHtEZt5WPv858IixNoyIUyJic0RsvvPOO2cmOkmSKqRlvdPLWYXG\nHPM1Mz+Smb2Z2XvwwRO27UuSNOvMdBK/vZzjl/LfO2b4+JIkdYyZTuIXAK8qn78K+PIMH1+SpI7R\ntCQeEecBlwKLI+KWiHgdxXSJyyLieuBZ5bIkSZqCpt0nnpknjbHq+GYdU5Kk2cRhVyVJqiiTuCRJ\nFWUSlySpokzikiRVVCVmMZMkTd6ZF13H2Rdff7/y044/kpXLFrUgIjVaFAOntbfe3t50FjNJmpoT\nz7kUgM+felyLI1G9ImJLZvZOtJ3V6ZIkVZRJXJKkijKJS5JUUSZxSZIqyiQuSVJFmcQlSaook7gk\nSRVlEpckqaJM4pIkVZRJXJKkijKJS5JUUSZxSZIqyiQuSVJFmcQlSaoo5xOXKsK5oSWNZBKXKmLl\nskWsXLbIuaEl3cfqdEmSKsokLklSRZnEJUmqKJO4JEkVZRKXJKmiTOKSJFWUSVySOtiGvn76dgxw\n2fa7Wbp2Exv6+lsdkhrIJC5JHWpDXz9r1m9l9+BeAPoHdrFm/VYTeQcxiUtSh1q3cRu79gwOK9u1\nZ5B1G7e1KCI1mklckjrUrQO7JlWu6jGJS1KHmt/dNalyVY9JXJI61Krli+maO2dYWdfcOaxavrhF\nEanRnABFkjrUiiU9AJx+/lXsHtxLT3cXq5Yvvq9c1WcSl6QOtmJJD+ddvgNw5rtOZHW6JEkVZRKX\nJKmiTOKSJFWUSVySpIoyiUuSVFEmcUmSKsokLklSRZnEJUmqKJO4VCHODS2plklcqgjnhpY0kklc\nqgjnhpY0kklcqgjnhpY0kklcqgjnhpY0kklcqgjnhpY0klORShXh3NCSRjKJSxXi3NCSarWkOj0i\nVkbEjyPi6og4LyIe1Io4JEmqshlP4hHRA/wl0JuZjwHmAC+d6TgkSaq6VnVsewDQFREPAPYDbm1R\nHJIkVdaMJ/HM7Af+GdgB3Ab8KjO/MdNxSJJUda2oTj8QeBFwGDAfeHBEvHyU7U6JiM0RsfnOO++c\n6TAlSWp7rahOfxawPTPvzMw9wHrgKSM3ysyPZGZvZvYefPDBMx6kJEntrhW3mO0AnhwR+wG7gOOB\nzS2IQ5I62pkXXcfZF19/3/LC1RcCcNrxR7Jy2aJWhaUGmvEknpmXRcT5wBXAvUAf8JGZjkOSOt3K\nZYtM1h2uJYO9ZOY7gXe24tiSJHUKx06XJKmiTOKSJFWUSVySpIpyAhRJajMje5UPsVe5RjKJS1Kb\nGepVfuI5lwLOWKexWZ0uSVJFmcQlSaook7gkSRVlEpckqaJM4pIkVZRJXJKkijKJS5JUUSZxSZIq\nysFepIpwbmi1gqPHtTeTuFQRzg2tVnD0uPZmdbokSRVlEpckqaJM4pIkVZRJXJKkijKJS5JUUSZx\nSZIqyiQuSVJFeZ+4pJZwEBFp+kziklrCQUTGt6Gvn74dA+we3MvStZtYtXwxK5b0tDostRmr0yWp\nzWzo62fN+q3sHtwLQP/ALtas38qGvv4WR6Z2YxKXpDazbuM2du0ZHFa2a88g6zZua1FEald1VadH\nxIHAfGAXcGNm7m1qVJI0i906sGtS5Zq9xkziEXEA8CbgJGBf4E7gQcAjIuIHwIcy85IZiVKSZpH5\n3V30j5Kw53d3tSAatbPxqtPPB24G/iQzF2fmUzOzNzMPAdYCL4qI181IlJI0i6xavpiuuXOGlXXN\nncOq5YtbFJHa1ZhX4pm5bJx1W4AtTYlIkma5oV7op59/FbsH99LT3WXvdI1qwjbxiAjgZODwzHxP\nRCwA/iAzL296dJI0S61Y0sN5l+8AvP1OY6und/qHgOMo2sYBdgIfbFpEkiSpLvX0Tj82Mx8fEX0A\nmfnLiNi3yXFJkqQJ1HMlvici5gAJEBEHA95iJkmzxNDocZdtv5ulazc56EwbqSeJ/yvwJeDhEfFe\n4LvAGU2NSpLUFhw9rr1NWJ2emedGxBbgeCCAFZl5bdMjkyS13Hijx9lbvvXq6Z3+6cx8BfCTUcok\nSR3M0ePaWz3V6UfVLpTt409oTjiSpHYy1ihxjh7XHsZM4hGxJiJ2Ao+NiF9HxM5y+Q7gyzMWoSSp\nZRw9rr2NmcQz84zMnAesy8z9M3Ne+XhYZq6ZwRglSS2yYkkPZ5xwNPvOKdJFT3cXZ5xwtO3hbaKe\n+8S/FhFPG1mYmd9uQjySpDbj6HHtq54kvqrm+YOAJ1GMm/7MpkQkSZLqUs8tZi+oXY6IQ4CzmhaR\n1GHOvOg6zr74+vuVn3b8kaxctqgFEbWPoUFEdg/uZenaTU7yIU1SPVfiI90C/HGjA5E61cpli1i5\nbBEnnnMpYHXkkLEGEQFM5FKd6rlP/P2UQ65SdIQ7BriimUFJ6nwOIiJNXz1X4ptrnt8LnJeZ32tS\nPJJmCQcRkaavnjbxT85EIJJml/ndXfSPkrAdRESq34QjtkXE0oi4KCKui4gbImJ7RNwwE8FJ6lwO\nIiJNXz3V6R8DVlLcVjY4wbaSVJehdu/Tz7+K3YN76enusne6NEn1JPFfZebXmh6JpFnHQUSk6akn\niV8SEeuA9cD/DBVmpj3UJakJRo4tsHD1hYBjC+j+6knix5b/9taUJY7YJklNMTS2gDSRenqnP2Mm\nApEkSZMzZhKPiJdn5mci4q9GW5+Z75vqQSOiG/go8BiKq/rXZualU92fJEmz0XhX4g8u/53XhOOe\nDXw9M18cEfsC+zXhGJIkdbQxk3hmnlP+++5GHjAiDgCeBry63P9uYHcjjyFJ0mxQz9jphwFvARbW\nbp+ZL5ziMQ8D7gT+IyIeR3H/+WmZec+I454CnAKwYMGCKR5KkqTONeGIbcAG4Ebg/cC/1Dym6gHA\n44EPZ+YS4B5g9ciNMvMjmdmbmb0HH3zwNA4nSVJnqucWs99l5r828Ji3ALdk5mXl8vmMksQlSdL4\n6kniZ0fEO4Fv0IDBXjLz5xFxc0QszsxtwPHANVPZlyRJs1k9Sfxo4BUUg7vsLcumO9jLW4Bzy57p\nNwCvmca+JEmalepJ4i8BDi97kTdEZl7J8BHgpI62oa+fvh0D7B7cy9K1m5zoQ5XhELDtrZ4kfjXQ\nDdzR5FikjrShr58167eye7CoyOof2MWa9VsBTORqew4B297qSeLdwE8i4ocMbxOf6i1m0qyybuM2\ndu0ZPovvrj2DrNu4zSQuaVrqSeLvbHoUUge7dWDXpMolqV71TIDyrZkIROpU87u76B8lYc/v7mpB\nNJI6ST2DvUiahlXLF9M1d86wsq65c1i1fHGLIpLUKeqpTpc0DUPt3qeffxW7B/fS091l73RJDTFu\nEo+IOcCnMvPkGYpH6kgrlvRw3uU7APj8qce1OBpJnWLcJJ6ZgxFxaETs28j7xCXJ+4+l6aunOv0G\n4HsRcQHFZCUAZOb7mhaVpI7n/cfS9NWTxH9WPvYB5jU3HEmSVK96bjF7N0BE7JeZv21+SJIkqR4T\n3mIWEcdFxDXAT8rlx0XEh5oemSRJGlc994mfBSwHfgGQmT8CntbMoCRJ0sTqGuwlM28eUTQ46oaS\nJGnG1NOx7eaIeAqQETEXOA24trlhSZKkidRzJf5G4E1AD3ArcEy5LEmSWqie3ul3AY7YJklSm6mn\nd/rhEfGViLgzIu6IiC9HxOEzEZwkSRpbPdXpnwX+E3gkMB/4AnBeM4OSJEkTqyeJ75eZn87Me8vH\nZ4AHNTswSZI0vnp6p38tIlYDnwMSOBH4akQ8FCAz725ifJIkaQz1JPG/KP89dUT5SymSuu3jkiS1\nQD290w+biUAkSdLk1DVimyRJaj/1VKdLmoYzL7qOsy++/r7lhasvBOC04490Pm1J0xKZ2eoYJtTb\n25ubN29udRiSNK6RP9iG+INNkxURWzKzd8LtJkriEbEUuDIz74mIlwOPB87OzJsaE+rETOKSquTE\ncy4F4POnHtfiSFRV9SbxetrEPwz8NiIeB7wN+BnwqWnGJ0mSpqmeJH5vFpfrLwI+kJkfBOY1NyxJ\nkjSRejq27YyINcDLgadFxD7A3OaGJbWe7ZuS2l09SfxE4GXA6zLz5xGxAFjX3LCk1lu5bBErly2y\nfVNS26pnsJefA++rWd6BbeKSJLXcmEk8InZSDKt6v1VAZub+TYtKkiRNaMwknpl2XpMkqY2NdyX+\n0PFe6OxlkiS11nht4lsoqtNjlHXOXiZJUouNV53u7GWSJLWxuiZAiYgDgSOBBw2VZea3mxWUJEma\n2IRJPCJeD5wGPAq4EngycCnwzOaGJkmSxlPPsKunAU8EbsrMZwBLgIGmRiVJkiZUTxL/XWb+DiAi\nHpiZPwEWNzcsSZI0kXqS+C0R0Q1sAC6KiC8DMzYNqSRVyYa+fvp2DHDZ9rtZunYTG/r6Wx2SOlg9\nw67+efn0XRFxCXAA8PWmRiVJFbShr58167eye3AvAP0Du1izfisAK5b0tDI0dagxr8Qj4okR8Zza\nssz8FnAvcHSzA5Okqlm3cRu79gwOK9u1Z5B1G7e1KCJ1uvGq0/8RuGaU8h/jLGaSdD+3DuyaVLk0\nXeMl8XmZeb+277LsoOaFJLUP2zc1GfO7uyZVLk3XeEn8wHHW7dfoQKR2M1b7polcY1m1fDFdc+cM\nK+uaO4dVy72hR80xXhL/74h4b0TcN3Z6FN4DbGp+aFJr2b6pyVqxpIczTjiafecUf1p7urs444Sj\n7dSmphmvd/rbgI8CP42IK8uyxwGbgdc3OzCp1Wzf1FSsWNLDeZfvAODzpx7X4mjU6cabAOUe4KSI\nOBw4qiz+cWbeMCORSS02v7uL/lEStu2bktrFeLeYLQTIzBsy8yvl44aa9RERj2p+iFJr2L4pqd2N\nV52+LiL2Ab5MMbf4nRSzmB0BPAM4HngncEuzg5RaYagd8/Tzr2L34F56urtYtXyx7ZuS2sZ41ekv\niYhHAycDrwUeCfwWuBb4KvDeoTHVpyIi5lC0r/dn5vOnuh+pmWzflNTOxh12NTOvAf62Scc+jeIH\nwf5N2r8kSR2tnglQGq5sS38eRe93SZI0BRNOgNIkZwGnA/PG2iAiTgFOAViwYMEMhSVptjrzous4\n++Lr71d+2vFHsnLZohZEJE1sxpN4RDwfuCMzt0TE08faLjM/AnwEoLe3N2coPEmz1Mpli1i5bBEn\nnnMpYB8IVcOE1enlrWQvj4i/K5cXRMSTpnHMpcALI+JG4HPAMyPiM9PYnyRJs1I9beIfAo4DTiqX\ndwIfnOoBM3NNZj4qMxcCLwU2ZebLp7o/SZJmq3qq04/NzMdHRB9AZv4yIvZtclySJGkC9STxPeU9\n3QkQEQcDextx8Mz8JvDNRuxLkqTZpp7q9H8FvgQ8PCLeC3wXOKOpUUmSpAlNeCWemedGxBaKYVYD\nWJGZ1zY9MkmSNK4Jk3hEfDozXwH8ZJQySVJp5L3mC1dfCHivuZqnnjbxo2oXyvbxJzQnHEmqrqF7\nzaWZMt5UpGsiYifw2Ij4dUTsLJfvoJjZTJIktdCYSTwzz8jMecC6zNw/M+eVj4dl5poZjFGSJI2i\nno5tayLiQOBIivnEh8q/3czApFazfbMaHPNcs1k9HdteTzFt6KOAK4EnA5cCz2xuaFJr2b5ZDY55\nrtmsnvvETwOeCNyUmc8AlgADTY1KklpgQ18/fTsGuGz73Sxdu4kNff2tDkkaVz1J/HeZ+TuAiHhg\nZv4EWNzcsCRpZm3o62fN+q3sHiwGpOwf2MWa9VtN5Gpr9STxWyKiG9gAXBQRXwZuam5YkjSz1m3c\nxq49g8PKdu0ZZN3GbS2KSJpYPR3b/rx8+q6IuAQ4APh6U6OSpBl268CuSZVL7WDcK/GImBMR943U\nlpnfyswLMnN380OTpJkzv7trUuVSOxj3SjwzByNiW0QsyMwdMxWUNF3edqTJWrV8MWvWbx1Wpd41\ndw6rltsFSO2rnmFXDwR+HBGXA/cMFWbmC5sWlTRN3nakyVqxpAeA08+/it2De+np7mLV8sX3lUvt\nqJ4k/n+bHoUktYEVS3o47/Ki0tEffqqCejq2fWsmApEkSZNTzy1mkiSpDZnEJUmqKJO4pMpzuFTN\nVmO2iUfEViBHWwVkZj62aVFJUp3GGi4VsGe5Ot54HdueP2NRSNIUjTdcqklcnW7MJJ6Zjo8uqe05\nXKpms/Gq03cyfnX6/k2LSpLqNL+7i/5RErbDpWo2GLNjW2bOy8z9R3nMM4FLaherli+ma+6cYWUO\nl6rZop4R2wCIiIcDDxpadix1Se3A4VI1m02YxCPihcC/APOBO4BDgWuBo5obmjQ9Q7cd7R7cy9K1\nm/zD3sEcLlWzVT33if898GTgusw8DDge+EFTo5Kmaazbjrx/WFInqac6fU9m/iIi9omIfTLzkog4\nq+mRSdPgbUearJHT1y5cfSHg9LVqb/Uk8YGIeAjwbeDciLiDmilJpXbkbUearKHpa6Uqqac6/UXA\nb4GVwNeBnwEvaGZQ0nSNdXuRtx1J6iRjJvGIOCIilmbmPZm5NzPvzcxPAlcA3TMXojR53nYkaTYY\n70r8LODXo5T/qlwnta0VS3o444Sj2XdO8RXv6e7ijBOOtj1cUkcZr038EZm5dWRhZm6NiIVNi0hq\nEG87ktrHyI6DQ+w4OD3jJfHxqsxtWJQk1W2o4+CJ51wK+MO6UcarTt8cEW8YWRgRrwe2NC8kSZJU\nj/GuxN8KfCkiTub3SbsX2Bf482YHJkmSxjfeVKS3A0+JiGcAjymLL8zMTTMSmSRJGteEg71k5iXA\nJTMQiyRJmoS6ZzGTpHbkcKmazUzikirN4VI1m9Uz7KokSWpDJnFJkirKJC5JUkWZxCVJM2JDXz99\nOwa4bPvdLF27iQ19/a0OqfLs2KaOZI9lqb1s6Otnzfqt7B7cC0D/wC7WrC+m53BioqmLzGx1DBPq\n7e3NzZs3tzoMSdIULV27if6BXfcr7+nu4nurn9mCiNpbRGzJzN6JtrM6XZLUdLeOksDHK1d9TOKS\npKab3z365Jdjlas+JnFJUtNrZu4hAAAMmElEQVStWr6YrrlzhpV1zZ3DquWLWxRRZ7BjmySp6YY6\nr51+/lXsHtxLT3cXq5YvtlPbNJnEJUkzYsWSHs67fAcAnz/1uBZH0xlmvDo9Ig6JiEsi4pqI+HFE\nnDbTMUiS1AlacSV+L/C2zLwiIuYBWyLiosy8pgWxSJJUWTN+JZ6Zt2XmFeXzncC1gI0ikiRNUkt7\np0fEQmAJcFkr45AkqYpa1rEtIh4CfBF4a2b+epT1pwCnACxYsGCGo1OrjBwudYjDpUrS/bUkiUfE\nXIoEfm5mrh9tm8z8CPARKIZdncHw1EIrly1i5bJFnHjOpYA9WCVpPK3onR7Ax4BrM/N9M318SZI6\nRSvaxJcCrwCeGRFXlo/ntiAOSZIqbcar0zPzu0DM9HElSeo0jp0uSVJFOeyqJKnpRt55snD1hYB3\nnkyXSVyS1HRDd56osaxOlySpokziajsb+vrp2zHAZdvvZunaTWzo6291SJLUlkziaisb+vpZs34r\nuwf3AtA/sIs167eayCVpFCZxtZV1G7exa8/gsLJdewZZt3FbiyKSpPZlEldbuXVg16TKJWk2M4mr\nrczv7ppUuSTNZiZxtZVVyxfTNXfOsLKuuXNYtXxxiyKSpPblfeJqKyuW9ABw+vlXsXtwLz3dXaxa\nvvi+cknS75nE1XZWLOnhvMt3AE5FKknjsTpdkqSKMolLklRRJnFJkirKJC5JUkWZxCVJqiiTuCRJ\nFWUSlySporxPXA1x5kXXcfbF19+v/LTjj2TlskVT3s/C1RdOaT+SNBtEZrY6hgn19vbm5s2bWx2G\n6nDiOZcCDtIiSdMREVsys3ei7axOlySpoqxOlyRVRqOa7jqFSVySVBkrly1i5bJFNt2VrE6XJKmi\nTOJqmA19/fTtGOCy7XezdO0mNvT1tzokSepoJnE1xIa+ftas38ruwb0A9A/sYs36rSZySWoik7ga\nYt3GbezaMzisbNeeQdZt3NaiiCSp85nE1RC3DuyaVLkkTZVNd79nEldDzO/umlS5JE2FTXfDmcTV\nEKuWL6Zr7pxhZV1z57Bq+eIWRSSpE9l0N5z3iashVizpAeD0869i9+Beerq7WLV88X3lktQINt0N\nZxJXw6xY0sN5l+8AHIBBUnPM7+6if5SEPVub7qxOlyRVhk13w3klLkmqDJvuhjOJS5Iqxaa73zOJ\ny1mBJKmiTOJqyKxAI38ILFx9IeAPAUlqJpO4GmLoh4AkaeaYxCVJlWGt33CRma2OYUK9vb25efPm\nVofR0Tb09dvbU9KsUIV+QBGxJTN7J9zOJK6hsYhrhzLsmjuHM0442kQuqWNNpx9Qs9WbxK1O17hj\nEZvEJWl07XBFbxKXYxFLmnWGpjPdPbiXpWs3TakJsRF39kyXSbzCGvUr0LGIJc0mY01nCkw6kTfi\nx8B0OHZ6ha1ctogb1z6PYw97KMce9lBuXPs8blz7vElX4zgWsaTZpFHTmbbD3OYm8Yob+hV42fa7\nWbp205S+PCuW9HDGCUez75zi69DT3WWnNkkdq1FNiO0wt7nV6RXWyCohxyKWNFs0qgmxHfoTeSVe\nYe3wK1CSqqZRTYhjJf2Z7E/klXiFNepXoCMgSZpNGjWd6arli0cdY2Mm+xPNqiS+oa+fdRu3cevA\nLua3eFSyRsTSqCohxz2XNNs0ogmxHeY2nzUjtp36qc1svOb2+5Uvf/QjOOeVEw6KM8x0E3CjYnGk\nNUmavEYP0tKM+8QdsW2Eq2/99aTKxzIyAfcP7OKtn7+Sr229re4E3KhY2uFXoCRVTSfVPrYkiUfE\ns4GzgTnARzNzbbOP2aj240Yk4Eb2aLRXuSS1Rjv0J5rxJB4Rc4APAsuAW4AfRsQFmXlNM4/bTrcU\nOEKaJFVfO1zRt+IWsycBP83MGzJzN/A54EXNPmg73VLwmPn7T6p8LGdedB0LV1/IZdvv5rLtd7Nw\n9YUsXH0hZ1503aT2I0mqplZUp/cAN9cs3wIcO3KjiDgFOAVgwYIF0z7o9rvuGfWe6u133TOp/Txm\n/v6jXkVPJgGf88rehvROb4dfgZKk1pnx3ukR8WLg2Zn5+nL5FcCxmfnmsV7TbvOJt9OtapKkztPO\nvdP7gUNqlh9VllXGiiU9Jm1JUsu1ok38h8CREXFYROwLvBS4oAVxSJJUaTN+JZ6Z90bEm4GNFLeY\nfTwzfzzTcUiSVHUtuU88M78KfLUVx5YkqVM4i5kkSRVlEpckqaJM4pIkVZRJXJKkijKJS5JUUSZx\nSZIqyiQuSVJFmcQlSaook7gkSRU147OYTUVE3Anc1MBdHgTc1cD9aTjPb/N5jpvL89tcnt+JHZqZ\nB0+0USWSeKNFxOZ6pnjT1Hh+m89z3Fye3+by/DaO1emSJFWUSVySpIqarUn8I60OoMN5fpvPc9xc\nnt/m8vw2yKxsE5ckqRPM1itxSZIqb9Yl8Yh4dkRsi4ifRsTqVsfTaSLixojYGhFXRsTmVsdTdRHx\n8Yi4IyKuril7aERcFBHXl/8e2MoYq26Mc/yuiOgvv8dXRsRzWxljVUXEIRFxSURcExE/jojTynK/\nww0yq5J4RMwBPgg8B3g0cFJEPLq1UXWkZ2TmMd5C0hCfAJ49omw1cHFmHglcXC5r6j7B/c8xwJnl\n9/iYzPzqDMfUKe4F3paZjwaeDLyp/Jvrd7hBZlUSB54E/DQzb8jM3cDngBe1OCZpTJn5beDuEcUv\nAj5ZPv8ksGJGg+owY5xjNUBm3paZV5TPdwLXAj34HW6Y2ZbEe4Cba5ZvKcvUOAl8IyK2RMQprQ6m\nQz0iM28rn/8ceEQrg+lgb46Iq8rqdqt7pykiFgJLgMvwO9wwsy2Jq/mempmPp2iyeFNEPK3VAXWy\nLG4v8RaTxvsw8IfAMcBtwL+0Npxqi4iHAF8E3pqZv65d53d4emZbEu8HDqlZflRZpgbJzP7y3zuA\nL1E0Yaixbo+IRwKU/97R4ng6TmbenpmDmbkX+Hf8Hk9ZRMylSODnZub6stjvcIPMtiT+Q+DIiDgs\nIvYFXgpc0OKYOkZEPDgi5g09B/4MuHr8V2kKLgBeVT5/FfDlFsbSkYYSTOnP8Xs8JRERwMeAazPz\nfTWr/A43yKwb7KW8VeQsYA7w8cx8b4tD6hgRcTjF1TfAA4DPen6nJyLOA55OMevT7cA7gQ3AfwIL\nKGb3+4vMtGPWFI1xjp9OUZWewI3AqTVtuKpTRDwV+A6wFdhbFr+dol3c73ADzLokLklSp5ht1emS\nJHUMk7gkSRVlEpckqaJM4pIkVZRJXJKkijKJS20gIn4zYvnVEfGBGY7hJRFxbURcMqJ8YUS8rBmx\nlfv2HmxpikziUgeLiAdMYvPXAW/IzGeMKF8IvOz+m0tqNZO41ObKq9VN5WQcF0fEgrL8ExHx4prt\nflP++/SI+E5EXABcM8r+TirnfL86Iv6xLPs74KnAxyJi3YiXrAX+pJxXe2VZNj8ivl7OB/1PNfv+\ns4i4NCKuiIgvlGNmjzz+EyLiRxHxI+BNI97nd8rXXhERTynLPxURK2q2OzcinH1QwiQutYuuMkle\nGRFXAu+pWfd+4JOZ+VjgXOBf69jf44HTMnNRbWFEzAf+EXgmxYhkT4yIFZn5HmAzcHJmrhqxr9XA\nd8p5tc8sy44BTgSOBk6MiEMi4iDgHcCzyklwNgN/NUps/wG8JTMfN6L8DmBZ+doTa97nx4BXl/Ef\nADwFuLCOcyB1vMlUtUlqnl2ZeczQQkS8GugtF48DTiiffxr4JyZ2eWZuH6X8icA3M/PO8jjnAk+j\nGMp1Mi7OzF+V+7gGOBToBh4NfK8YMpt9gUtrXxQR3UB3OYf30Pt5Tvl8LvCBiDgGGAQWAWTmtyLi\nQxFxMPC/gC9m5r2TjFfqSCZxqbrupaxNi4h9KJLmkHuafOz/qXk+SPG3JICLMvOkKe5zJcXY5Y+j\neF+/q1n3KeDlFJMWvWaK+5c6jtXpUvv7PkXyAjiZYkIJKCbmeEL5/IUUV7ITuRz404g4KCLmACcB\n35rgNTuBeXXs+wfA0og4Au6b1W5YdX5mDgAD5cQYULyfIQcAt5XTf76CYpKiIZ8A3lru437t/NJs\nZRKX2t9bgNdExFUUye20svzfKRLyjyiq3Ce8+i5n4loNXAL8CNiSmRNNA3kVMFh2Rls51kZlFf2r\ngfPKWC8F/miUTV8DfLBs+4+a8g8Bryrfzx/Vvp/MvB24lqI9XVLJWcwktb2I2I9iOsvHD7XFS/JK\nXFKbi4hnUVyFv98ELg3nlbgkSRXllbgkSRVlEpckqaJM4pIkVZRJXJKkijKJS5JUUSZxSZIq6v8D\nRlBR2YkXhUEAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 576x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gJmdGvx250rV",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xk7WPdiwtWKD",
        "colab_type": "text"
      },
      "source": [
        "## 4. The purpose of dividing by the area in line 10 even after having subtracted the maximum of the unnormalized log posterior in line 7"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i4UELqOu52a6",
        "colab_type": "text"
      },
      "source": [
        "The area under the curve has to be 1. Hence, we divide by the area. \n",
        "\n",
        "Subtracting the maximum of the unnormalized log posterior is simply a rescaling method to ensure computation efficiency and accuracy, but does not link anything with the condition $area = 1$. Therefore, we still have to divide the area at the end. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "itNW-GA4TrYI",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def compute_posterior(parameter_values, data):\n",
        "    # take the log of the prior\n",
        "    log_prior = np.log(prior(parameter_values))\n",
        "    \n",
        "    # for each lambda (parameter_values), we compute the sum of all of log of\n",
        "    # each datum likelihood by using np.sum\n",
        "    log_likelihood = np.array([np.sum(np.log(likelihood(param, data))) for param in parameter_values])\n",
        "    \n",
        "    # the log of posterior is the sum of log of prior and log of likelihood\n",
        "    unnormalized_log_posterior = log_prior + log_likelihood\n",
        "\n",
        "    # rescaling to ensure no numerical overflow occurs. \n",
        "    unnormalized_log_posterior -= max(unnormalized_log_posterior)\n",
        "\n",
        "    # as we are using log value the whole time, we now \n",
        "    # convert back to the original values.\n",
        "    unnormalized_posterior = np.exp(unnormalized_log_posterior)\n",
        "\n",
        "    # compute the area under the curve and divide the unnormalized_posterior\n",
        "    # by the area to ensure area under the curve = 1 \n",
        "    # \n",
        "    # (the -= max calculation above is an additional rescaling step to ensure\n",
        "    # computation efficiency and accuracy and avoid numerical overflow, and \n",
        "    # does not have anything relate to area under the curve = 1)\n",
        "    area = sp.integrate.trapz(unnormalized_posterior, parameter_values)\n",
        "    posterior = unnormalized_posterior / area\n",
        "    return posterior"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "X-KNJBtl6sDS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5zQz1ABmtdJk",
        "colab_type": "text"
      },
      "source": [
        "## 5. Create an example where not taking logarithms would cause a problem. \n",
        "Create a prior, likelihood, and data set that fails to produce the correct posterior when we don’t take logs. Show all your code and visualize your results on one or more plots.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VUbD6G_f7sv8",
        "colab_type": "text"
      },
      "source": [
        "Not taking the log would possibly face the same problem: numerical overflow.\n",
        "\n",
        "Instead of taking log and take the sum, we then have to compute: $\\prod_{i=1}^{n} P(datum_i|λ)$. \n",
        "\n",
        "For 5856 data points given, the product of all the likelihood of the datum given the parameter can exceed the numpy numeric memory (numerical overflow). Hence, causing warnings and inaccuracies. \n",
        "\n",
        "Therefore, it's a better practice to take the log and take the sum as it will less likely to be numerical overflow. \n",
        "\n",
        "Example: Given the same dataset, we modify the code (not using log) to see the warning.  The warning occurs, which makes some values will be represented as **inf**. Then, these values are unprocessable, which lead to another error **invalid value encountered in true_divide**. This will create missing values at some intervals.  \n",
        "\n",
        "As we can see from the code below, at hour[11], the value got too big, lead to overflow. Then, if we look at the graph, the values for hour = 11 does not exist. This cause the figure to be incorrect. Hence, taking the log is the practical move. Also, we can see the graph is similar to Section 2, part 3: both graphs got missing values at hour[11] due to numerical overflow"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "THgmtlrc7rdN",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 922
        },
        "outputId": "09bf9911-9635-46cf-9cd4-c25eeaa17a5e"
      },
      "source": [
        "def compute_posterior(parameter_values, data):\n",
        "    priors = prior(parameter_values)\n",
        "    likelihoods = np.array([np.prod(likelihood(param, data)) for param in parameter_values])\n",
        "    unnormalized_posterior = priors * likelihoods\n",
        "    area = sp.integrate.trapz(unnormalized_posterior, parameter_values)\n",
        "    posterior = unnormalized_posterior / area\n",
        "    return posterior\n",
        "\n",
        "# Create empty lists for storing the data\n",
        "#\n",
        "# upper_percentiles_diff to store the difference between the 99% \n",
        "# and the posterior mean\n",
        "#\n",
        "# lower_percentiles_diff to store the difference between the 01% \n",
        "# and the posterior mean\n",
        "#\n",
        "# lambda_ to store the posterior mean\n",
        "upper_percentiles_diff = []\n",
        "lower_percentiles_diff = []\n",
        "lambda_ = []\n",
        "\n",
        "for i, waiting_times_hour in  enumerate(waiting_times_per_hour):\n",
        "    print(i)\n",
        "    # initialize the lambdas values\n",
        "    lambdas = np.linspace(0, 20, 501)[1:]\n",
        "    #\n",
        "    # compute the posterior by using the compute_posterior function\n",
        "    # with the initialized lambdas and the data (of each hour)\n",
        "    posterior = compute_posterior(lambdas, waiting_times_hour)\n",
        "    #\n",
        "    # Add the newly calculated posterior mean to the lambda_ list\n",
        "    lambda_.append(sp.integrate.trapz(posterior * lambdas, lambdas))\n",
        "    #\n",
        "    # Calculate the difference between the posterior mean with the 99% percentile\n",
        "    # and append to the upper_percentiles_diff list\n",
        "    upper_percentiles_diff.append(compute_percentile(lambdas, posterior, 0.99) - lambda_[-1])\n",
        "    #\n",
        "    # Calculate the difference between the posterior mean with the 1% percentile\n",
        "    # and append to the upper_percentiles_diff list\n",
        "    lower_percentiles_diff.append(lambda_[-1] - compute_percentile(lambdas, posterior, 0.01))\n",
        "\n",
        "# initialize the graph\n",
        "plt.figure(figsize=(8, 6))\n",
        "# plot the error bar with the mean at lambda_ and the range from the mean derives\n",
        "# from the upper_percentiles_diff and lower_percentiles_diff\n",
        "plt.errorbar(list(range(24)), lambda_, yerr = [upper_percentiles_diff,lower_percentiles_diff]\n",
        "             , fmt=\"o\", capsize=3)\n",
        "# label the axes and plot\n",
        "plt.ylabel('Call rate (Calls per minute)')\n",
        "plt.xlabel('Hour of the day')\n",
        "plt.title('Mean and 98% confidence interval over call rate')\n",
        "plt.show()"
      ],
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "0\n",
            "1\n",
            "2\n",
            "3\n",
            "4\n",
            "5\n",
            "6\n",
            "7\n",
            "8\n",
            "9\n",
            "10\n",
            "11\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/numpy/core/fromnumeric.py:86: RuntimeWarning: overflow encountered in reduce\n",
            "  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)\n",
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in true_divide\n",
            "  \n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "stream",
          "text": [
            "12\n",
            "13\n",
            "14\n",
            "15\n",
            "16\n",
            "17\n",
            "18\n",
            "19\n",
            "20\n",
            "21\n",
            "22\n",
            "23\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAAGDCAYAAAA72Cm3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucXGV9+PHPlxh00cCioDUrIVBI\nWhEluooYa1VM492UnxYR7xfw91OLqQ1NrL96af2RNq1AvRWr1huiFWPEokZK8I5gwiJBMKAEAgty\nEVcjxiZsvr8/zlmcXfYyuzuzM2f283695pU5zzlzznfOTPY753me8zyRmUiSpOrZp9UBSJKkqTGJ\nS5JUUSZxSZIqyiQuSVJFmcQlSaook7gkSRVlElclRcSNEfGsVscxEyKiKyK+EhG/iogvRMTJEfGN\ncbb/ZkS8fiZjHE9E/FtE/N9Wx1GPTvteRURGxBHl809ExD+0OiY1lkm8Q5V/jHZHxEEjyvvK/9gL\nWxNZ80VEd0R8MiLuKB/vGrH+mIj4TpkUb6lNMBFxSET8ICLujoh/GfG6r0VE7wy9jVovBh4BPCwz\nX5KZ52bmn7UgjinJzDdm5t/Xs62Jpjo67QdPVZnEO9t24KShhYg4GtivdeHMmDMp3udC4EnAKyLi\nNTXrPwt8G3go8KfA/4mIF5br1gCfBA4DVgwl7Yg4EdiemZtn5B0MdyhwXWbe24JjV0pEPKDVMUxH\nu8TfLnFoYibxzvZp4JU1y68CPlW7QUQ8MCL+OSJ2RMTtZdVnV7nuwIj4r4i4MyJ+WT5/VM1rvxkR\nfx8R34uInRHxjZFX/jXbTmtfEfGKiLgpIn4REX87wft+AfBPmfnbzLwR+Bjw2pr1C4FzM3MwM38G\nfBc4qlx3GLApM38F/BA4PCL2B1YDb5/guETEURFxUXklf3tEvL3mPJ8VEbeWj7Mi4oHluqeXNQJv\nK2sObhv60RER7wb+DjgxIn4TEa+LiFdHxHdrjrksIn5S1ix8AIgRMb02Iq4tz/vGiDi0Zl1GxBsj\n4vqIGIiID0ZE1Kx/Q/nanRFxTUQ8viyfHxFfLD/P7RHxl+Ock/uurid4r6cAJwOnl+/1KxMdKyLe\nFRHnR8RnIuLXwNsjYldEPLRmmyURcVdEzI2IP4yITeX36K6IODciuif6XMv9HBARnyrjuCki3hER\n+5Sf7UBEPKZm24PLOB5eLj8/Iq4st/t+RDy2ZtsbI+JvIuIq4J4YJYGO8716UkRcWu73toj4QETs\nW8/7GbH/V5f/986MiF8A7xrvXEXEp4EFwFfKz+r0svzJ5fsbiIgfRcTTJxuLJikzfXTgA7gReBaw\nDfhjYA5wC8VVXQILy+3OBC6guCqdB3wFOKNc9zDgf1Fc1c4DvgBsqDnGN4GfAYuArnJ57RjxTHlf\nwKOB3wBPAx4IvA+4F3jWGMe6C3hSzfLfAr+sWf5/wFpgLrC4PC9PLNetA94MdAPXUyT3s4FX1XHO\n5wG3AW8DHlQuH1uuew/wA+DhwMHA94G/L9c9vXw/7yljei7wW+DAcv27gM/UHOfVwHfL5wcBOymq\n3OcCK8t9vb5c/yLgp+V34AHAO4Dv1+wrgf8q3+8C4E7g2eW6lwD9wBMpfhgcQfH92QfYQvHjYl/g\ncOAGYPkY5+UTwD/U+V7v27ZcHvdY5bnZA6wot+0CNgFvqNnHOuDfyudHAMsovkcHU9TInDXy/80Y\n7+NTwJfLz3UhcB3wunLdx4H31mz7JuDr5fMlwB3AsRT/D19VHueBNce8EjgE6Jrk9+oJwJPLz3Yh\ncC3w1hGf7xGjndsRx3h1+bm8pdxX12TPFdAD/KL8TPcpX/sL4OBW/z3s5EfLA/DRpA/290n8HcAZ\nwLOBi8r/oFn+hw/gHuAPa153HEW18Wj7PIbhyfCbwDtqlv/P0B+uOuKre18Uf8A/V7PuwcDucf7Y\nfgZYX/6xO4Lix8H/1Kx/CkViu7c8F++uWfdQ4PPAjygS4hLgkrJ8qBr+zWMc9ySgb4x1PwOeW7O8\nHLixfP50YBfwgJr1dwBPLp+/i7GT+CuBH9SsC4ofJUNJ/GuUiaZc3ociaR5aLifw1Jr1/wmsLp9v\nBE4b5b0cC+wYUbYG+I8x3vsnGJ7Ex3uv921bz7HKc/PtEetfT1GbMnQ+bgaeNkZsK2o/M8ZI4hTJ\ndzfw6JqyU4Fvls+fBfysZt33gFeWzz9M+YOtZv024E9rjvnacf6vjPm9GmXbtwJfqlmeTBLfMcG+\nxz1XwN8Anx7xmo3U8QPYx9Qftnt0vk9TJJ7DGFGVTvHrej9gS20NKsUfLCJiP4or9WcDB5br50XE\nnMwcLJd/XrO/3wIPGS2Iae5rPsUfYgAy856yym8sfwm8n+JK+hfAeZR9A8pq1q9TXG1/FvgD4PyI\nuD0zP5SZdwMnltvuQ3Hu3khRnX41xR+7KyLi4sy8dsRxD6FI1qOZD9xUs3xTWTbkFzm8zXvMcznK\nfmvPTUbEzTXrDwXOjuGd9ILiqmkonrHO+1jv51BgfkQM1JTNAb5TR7wwufdaz7FuHv4Svgi8PyIe\nSVGzs3do+4h4BEXNyp9Q/MjbB/hlHTEfRFFzMPIz7CmfXwLsFxHHArdT/Ej9Us17eFVEvKXmtfsy\n/PMf+R5qjfm9iohFFDVTvRT/lx9AUXMxFcNimMK5OhR4SUS8oKZsLsW5UZPYJt7hMvMmig5uz6W4\nOq11F8VV0VGZ2V0+DsjMoT+ob6Oobj42M/enqM6GEW2udZrOvm6j+ENWvKD4QfCwsTbOzLsz8+TM\n/IPMPIrie355ufpwYDAzP5WZ92bmLcDnKM7PSKdQXOVeDRwNbM7M3cDWcnmkm8v9j+ZWij9yQxaU\nZdM18txE7XIZ06k1n293ZnZl5vfr2PfNwB+OUb59xD7nZeZo53CyRk6rWM+xhr0mM38JfIPix9jL\nKGpxhrb5f+X2R5ffw5dT33fwLopq+5GfYX95zEGKWoyTysd/ZebOmvfw3hHvYb/MPG+c9z3yHIz1\nvfow8BPgyPL9vL3O9zOakTFMdK5G+6w+PeJ9Pjgz104xHtXBJD47vA54ZmbeU1uYmXuBfwfOrOmA\n0xMRy8tN5lEk+YHyCvad04hhOvs6H3h+RDy17LTzHsb57pYdch4WEXMi4jkUyXjotqXrik3iZWWn\npD+g+GN/1Yh9PJyiXfNdZdF24BkR8RCKq54bRjn0fwGPjIi3lp2d5pVXZlDUBryj7PB0EEUTwWcm\ncQ7GciFwVEScUHaI+kuK2oUh/wasiYijyvd1QES8pM59fxT464h4QhSOiKJT3OXAzrIzVld5nh8T\nEU9swPu5neEJa6rH+ixFU8OLy+dD5lH0r/hVRPQAq+oJqiZJv7f8XA8F/orhn+FnKb5LJ4845r8D\nb4yIY8vz+OCIeF5EzKvn2Iz/vZoH/Br4TUT8EfC/69xnPSY6VyM/q88AL4iI5eXn9KAoOjI+CjWN\nSXwWyMyf5di3Rv0NRfvwD6Lo3fvfFFfMAGdRdHC5i6JT1tenEcaU95WZP6ZIqJ+luPL8JUW771ie\nQHG1vJOiP8DJ5T7IzF8DJ1C0d/+SokPR1fw+yQ/5Z+A9mfmbcvkM4JkUVxtfGe18lldeyyh6x/+c\nojr/GeXqfwA2U/xY2ApcMcoxJy0z76LogLaWoungSIr22KH1XwL+Efhc+fleDTynzn1/AXgvxXnf\nCWwAHlomtOdTVBlvp/hMPwocMN33Q3EnwaPL3s0bpnGsCyjOxc8z80c15e8GHg/8iuIH0MjaqfG8\nhaIPyQ0UdzR8lqJDGwCZeVm5fj5FX4Sh8s3AG4APUHznfkrRLFOXCb5Xf01R27CT4sfC5yfxfiYy\n0bk6g+KH6UBE/HVm3kzRkfLtFB0kb6ZI/OaZJorf1zJJkqQq8ReSJEkVZRKXJKmiTOKSJFWUSVyS\npIoyiUuSVFGVGLHtoIMOyoULF7Y6DEmSZsSWLVvuysyDJ9quaUk8Ij5OcX/nHZn5mLJsaFzqhRTj\n7v5FObrSuBYuXMjmza2YAVKSpJkXETdNvFVzq9M/QTFOdq3VwMWZeSRwcbksSZKmoGlJPDO/Ddw9\novhFwCfL55+kmBVHkiRNwUx3bHtEZt5WPv858IixNoyIUyJic0RsvvPOO2cmOkmSKqRlvdPLWYXG\nHPM1Mz+Smb2Z2XvwwRO27UuSNOvMdBK/vZzjl/LfO2b4+JIkdYyZTuIXAK8qn78K+PIMH1+SpI7R\ntCQeEecBlwKLI+KWiHgdxXSJyyLieuBZ5bIkSZqCpt0nnpknjbHq+GYdU5Kk2cRhVyVJqiiTuCRJ\nFWUSlySpokzikiRVVCVmMZMkTd6ZF13H2Rdff7/y044/kpXLFrUgIjVaFAOntbfe3t50FjNJmpoT\nz7kUgM+felyLI1G9ImJLZvZOtJ3V6ZIkVZRJXJKkijKJS5JUUSZxSZIqyiQuSVJFmcQlSaook7gk\nSRVlEpckqaJM4pIkVZRJXJKkijKJS5JUUSZxSZIqyiQuSVJFmcQlSaoo5xOXKsK5oSWNZBKXKmLl\nskWsXLbIuaEl3cfqdEmSKsokLklSRZnEJUmqKJO4JEkVZRKXJKmiTOKSJFWUSVySOtiGvn76dgxw\n2fa7Wbp2Exv6+lsdkhrIJC5JHWpDXz9r1m9l9+BeAPoHdrFm/VYTeQcxiUtSh1q3cRu79gwOK9u1\nZ5B1G7e1KCI1mklckjrUrQO7JlWu6jGJS1KHmt/dNalyVY9JXJI61Krli+maO2dYWdfcOaxavrhF\nEanRnABFkjrUiiU9AJx+/lXsHtxLT3cXq5Yvvq9c1WcSl6QOtmJJD+ddvgNw5rtOZHW6JEkVZRKX\nJKmiTOKSJFWUSVySpIoyiUuSVFEmcUmSKsokLklSRZnEJUmqKJO4VCHODS2plklcqgjnhpY0kklc\nqgjnhpY0kklcqgjnhpY0kklcqgjnhpY0kklcqgjnhpY0klORShXh3NCSRjKJSxXi3NCSarWkOj0i\nVkbEjyPi6og4LyIe1Io4JEmqshlP4hHRA/wl0JuZjwHmAC+d6TgkSaq6VnVsewDQFREPAPYDbm1R\nHJIkVdaMJ/HM7Af+GdgB3Ab8KjO/MdNxSJJUda2oTj8QeBFwGDAfeHBEvHyU7U6JiM0RsfnOO++c\n6TAlSWp7rahOfxawPTPvzMw9wHrgKSM3ysyPZGZvZvYefPDBMx6kJEntrhW3mO0AnhwR+wG7gOOB\nzS2IQ5I62pkXXcfZF19/3/LC1RcCcNrxR7Jy2aJWhaUGmvEknpmXRcT5wBXAvUAf8JGZjkOSOt3K\nZYtM1h2uJYO9ZOY7gXe24tiSJHUKx06XJKmiTOKSJFWUSVySpIpyAhRJajMje5UPsVe5RjKJS1Kb\nGepVfuI5lwLOWKexWZ0uSVJFmcQlSaook7gkSRVlEpckqaJM4pIkVZRJXJKkijKJS5JUUSZxSZIq\nysFepIpwbmi1gqPHtTeTuFQRzg2tVnD0uPZmdbokSRVlEpckqaJM4pIkVZRJXJKkijKJS5JUUSZx\nSZIqyiQuSVJFeZ+4pJZwEBFp+kziklrCQUTGt6Gvn74dA+we3MvStZtYtXwxK5b0tDostRmr0yWp\nzWzo62fN+q3sHtwLQP/ALtas38qGvv4WR6Z2YxKXpDazbuM2du0ZHFa2a88g6zZua1FEald1VadH\nxIHAfGAXcGNm7m1qVJI0i906sGtS5Zq9xkziEXEA8CbgJGBf4E7gQcAjIuIHwIcy85IZiVKSZpH5\n3V30j5Kw53d3tSAatbPxqtPPB24G/iQzF2fmUzOzNzMPAdYCL4qI181IlJI0i6xavpiuuXOGlXXN\nncOq5YtbFJHa1ZhX4pm5bJx1W4AtTYlIkma5oV7op59/FbsH99LT3WXvdI1qwjbxiAjgZODwzHxP\nRCwA/iAzL296dJI0S61Y0sN5l+8AvP1OY6und/qHgOMo2sYBdgIfbFpEkiSpLvX0Tj82Mx8fEX0A\nmfnLiNi3yXFJkqQJ1HMlvici5gAJEBEHA95iJkmzxNDocZdtv5ulazc56EwbqSeJ/yvwJeDhEfFe\n4LvAGU2NSpLUFhw9rr1NWJ2emedGxBbgeCCAFZl5bdMjkyS13Hijx9lbvvXq6Z3+6cx8BfCTUcok\nSR3M0ePaWz3V6UfVLpTt409oTjiSpHYy1ihxjh7XHsZM4hGxJiJ2Ao+NiF9HxM5y+Q7gyzMWoSSp\nZRw9rr2NmcQz84zMnAesy8z9M3Ne+XhYZq6ZwRglSS2yYkkPZ5xwNPvOKdJFT3cXZ5xwtO3hbaKe\n+8S/FhFPG1mYmd9uQjySpDbj6HHtq54kvqrm+YOAJ1GMm/7MpkQkSZLqUs8tZi+oXY6IQ4CzmhaR\n1GHOvOg6zr74+vuVn3b8kaxctqgFEbWPoUFEdg/uZenaTU7yIU1SPVfiI90C/HGjA5E61cpli1i5\nbBEnnnMpYHXkkLEGEQFM5FKd6rlP/P2UQ65SdIQ7BriimUFJ6nwOIiJNXz1X4ptrnt8LnJeZ32tS\nPJJmCQcRkaavnjbxT85EIJJml/ndXfSPkrAdRESq34QjtkXE0oi4KCKui4gbImJ7RNwwE8FJ6lwO\nIiJNXz3V6R8DVlLcVjY4wbaSVJehdu/Tz7+K3YN76enusne6NEn1JPFfZebXmh6JpFnHQUSk6akn\niV8SEeuA9cD/DBVmpj3UJakJRo4tsHD1hYBjC+j+6knix5b/9taUJY7YJklNMTS2gDSRenqnP2Mm\nApEkSZMzZhKPiJdn5mci4q9GW5+Z75vqQSOiG/go8BiKq/rXZualU92fJEmz0XhX4g8u/53XhOOe\nDXw9M18cEfsC+zXhGJIkdbQxk3hmnlP+++5GHjAiDgCeBry63P9uYHcjjyFJ0mxQz9jphwFvARbW\nbp+ZL5ziMQ8D7gT+IyIeR3H/+WmZec+I454CnAKwYMGCKR5KkqTONeGIbcAG4Ebg/cC/1Dym6gHA\n44EPZ+YS4B5g9ciNMvMjmdmbmb0HH3zwNA4nSVJnqucWs99l5r828Ji3ALdk5mXl8vmMksQlSdL4\n6kniZ0fEO4Fv0IDBXjLz5xFxc0QszsxtwPHANVPZlyRJs1k9Sfxo4BUUg7vsLcumO9jLW4Bzy57p\nNwCvmca+JEmalepJ4i8BDi97kTdEZl7J8BHgpI62oa+fvh0D7B7cy9K1m5zoQ5XhELDtrZ4kfjXQ\nDdzR5FikjrShr58167eye7CoyOof2MWa9VsBTORqew4B297qSeLdwE8i4ocMbxOf6i1m0qyybuM2\ndu0ZPovvrj2DrNu4zSQuaVrqSeLvbHoUUge7dWDXpMolqV71TIDyrZkIROpU87u76B8lYc/v7mpB\nNJI6ST2DvUiahlXLF9M1d86wsq65c1i1fHGLIpLUKeqpTpc0DUPt3qeffxW7B/fS091l73RJDTFu\nEo+IOcCnMvPkGYpH6kgrlvRw3uU7APj8qce1OBpJnWLcJJ6ZgxFxaETs28j7xCXJ+4+l6aunOv0G\n4HsRcQHFZCUAZOb7mhaVpI7n/cfS9NWTxH9WPvYB5jU3HEmSVK96bjF7N0BE7JeZv21+SJIkqR4T\n3mIWEcdFxDXAT8rlx0XEh5oemSRJGlc994mfBSwHfgGQmT8CntbMoCRJ0sTqGuwlM28eUTQ46oaS\nJGnG1NOx7eaIeAqQETEXOA24trlhSZKkidRzJf5G4E1AD3ArcEy5LEmSWqie3ul3AY7YJklSm6mn\nd/rhEfGViLgzIu6IiC9HxOEzEZwkSRpbPdXpnwX+E3gkMB/4AnBeM4OSJEkTqyeJ75eZn87Me8vH\nZ4AHNTswSZI0vnp6p38tIlYDnwMSOBH4akQ8FCAz725ifJIkaQz1JPG/KP89dUT5SymSuu3jkiS1\nQD290w+biUAkSdLk1DVimyRJaj/1VKdLmoYzL7qOsy++/r7lhasvBOC04490Pm1J0xKZ2eoYJtTb\n25ubN29udRiSNK6RP9iG+INNkxURWzKzd8LtJkriEbEUuDIz74mIlwOPB87OzJsaE+rETOKSquTE\ncy4F4POnHtfiSFRV9SbxetrEPwz8NiIeB7wN+BnwqWnGJ0mSpqmeJH5vFpfrLwI+kJkfBOY1NyxJ\nkjSRejq27YyINcDLgadFxD7A3OaGJbWe7ZuS2l09SfxE4GXA6zLz5xGxAFjX3LCk1lu5bBErly2y\nfVNS26pnsJefA++rWd6BbeKSJLXcmEk8InZSDKt6v1VAZub+TYtKkiRNaMwknpl2XpMkqY2NdyX+\n0PFe6OxlkiS11nht4lsoqtNjlHXOXiZJUouNV53u7GWSJLWxuiZAiYgDgSOBBw2VZea3mxWUJEma\n2IRJPCJeD5wGPAq4EngycCnwzOaGJkmSxlPPsKunAU8EbsrMZwBLgIGmRiVJkiZUTxL/XWb+DiAi\nHpiZPwEWNzcsSZI0kXqS+C0R0Q1sAC6KiC8DMzYNqSRVyYa+fvp2DHDZ9rtZunYTG/r6Wx2SOlg9\nw67+efn0XRFxCXAA8PWmRiVJFbShr58167eye3AvAP0Du1izfisAK5b0tDI0dagxr8Qj4okR8Zza\nssz8FnAvcHSzA5Okqlm3cRu79gwOK9u1Z5B1G7e1KCJ1uvGq0/8RuGaU8h/jLGaSdD+3DuyaVLk0\nXeMl8XmZeb+277LsoOaFJLUP2zc1GfO7uyZVLk3XeEn8wHHW7dfoQKR2M1b7polcY1m1fDFdc+cM\nK+uaO4dVy72hR80xXhL/74h4b0TcN3Z6FN4DbGp+aFJr2b6pyVqxpIczTjiafecUf1p7urs444Sj\n7dSmphmvd/rbgI8CP42IK8uyxwGbgdc3OzCp1Wzf1FSsWNLDeZfvAODzpx7X4mjU6cabAOUe4KSI\nOBw4qiz+cWbeMCORSS02v7uL/lEStu2bktrFeLeYLQTIzBsy8yvl44aa9RERj2p+iFJr2L4pqd2N\nV52+LiL2Ab5MMbf4nRSzmB0BPAM4HngncEuzg5RaYagd8/Tzr2L34F56urtYtXyx7ZuS2sZ41ekv\niYhHAycDrwUeCfwWuBb4KvDeoTHVpyIi5lC0r/dn5vOnuh+pmWzflNTOxh12NTOvAf62Scc+jeIH\nwf5N2r8kSR2tnglQGq5sS38eRe93SZI0BRNOgNIkZwGnA/PG2iAiTgFOAViwYMEMhSVptjrzous4\n++Lr71d+2vFHsnLZohZEJE1sxpN4RDwfuCMzt0TE08faLjM/AnwEoLe3N2coPEmz1Mpli1i5bBEn\nnnMpYB8IVcOE1enlrWQvj4i/K5cXRMSTpnHMpcALI+JG4HPAMyPiM9PYnyRJs1I9beIfAo4DTiqX\ndwIfnOoBM3NNZj4qMxcCLwU2ZebLp7o/SZJmq3qq04/NzMdHRB9AZv4yIvZtclySJGkC9STxPeU9\n3QkQEQcDextx8Mz8JvDNRuxLkqTZpp7q9H8FvgQ8PCLeC3wXOKOpUUmSpAlNeCWemedGxBaKYVYD\nWJGZ1zY9MkmSNK4Jk3hEfDozXwH8ZJQySVJp5L3mC1dfCHivuZqnnjbxo2oXyvbxJzQnHEmqrqF7\nzaWZMt5UpGsiYifw2Ij4dUTsLJfvoJjZTJIktdCYSTwzz8jMecC6zNw/M+eVj4dl5poZjFGSJI2i\nno5tayLiQOBIivnEh8q/3czApFazfbMaHPNcs1k9HdteTzFt6KOAK4EnA5cCz2xuaFJr2b5ZDY55\nrtmsnvvETwOeCNyUmc8AlgADTY1KklpgQ18/fTsGuGz73Sxdu4kNff2tDkkaVz1J/HeZ+TuAiHhg\nZv4EWNzcsCRpZm3o62fN+q3sHiwGpOwf2MWa9VtN5Gpr9STxWyKiG9gAXBQRXwZuam5YkjSz1m3c\nxq49g8PKdu0ZZN3GbS2KSJpYPR3b/rx8+q6IuAQ4APh6U6OSpBl268CuSZVL7WDcK/GImBMR943U\nlpnfyswLMnN380OTpJkzv7trUuVSOxj3SjwzByNiW0QsyMwdMxWUNF3edqTJWrV8MWvWbx1Wpd41\ndw6rltsFSO2rnmFXDwR+HBGXA/cMFWbmC5sWlTRN3nakyVqxpAeA08+/it2De+np7mLV8sX3lUvt\nqJ4k/n+bHoUktYEVS3o47/Ki0tEffqqCejq2fWsmApEkSZNTzy1mkiSpDZnEJUmqKJO4pMpzuFTN\nVmO2iUfEViBHWwVkZj62aVFJUp3GGi4VsGe5Ot54HdueP2NRSNIUjTdcqklcnW7MJJ6Zjo8uqe05\nXKpms/Gq03cyfnX6/k2LSpLqNL+7i/5RErbDpWo2GLNjW2bOy8z9R3nMM4FLaherli+ma+6cYWUO\nl6rZop4R2wCIiIcDDxpadix1Se3A4VI1m02YxCPihcC/APOBO4BDgWuBo5obmjQ9Q7cd7R7cy9K1\nm/zD3sEcLlWzVT33if898GTgusw8DDge+EFTo5Kmaazbjrx/WFInqac6fU9m/iIi9omIfTLzkog4\nq+mRSdPgbUearJHT1y5cfSHg9LVqb/Uk8YGIeAjwbeDciLiDmilJpXbkbUearKHpa6Uqqac6/UXA\nb4GVwNeBnwEvaGZQ0nSNdXuRtx1J6iRjJvGIOCIilmbmPZm5NzPvzcxPAlcA3TMXojR53nYkaTYY\n70r8LODXo5T/qlwnta0VS3o444Sj2XdO8RXv6e7ijBOOtj1cUkcZr038EZm5dWRhZm6NiIVNi0hq\nEG87ktrHyI6DQ+w4OD3jJfHxqsxtWJQk1W2o4+CJ51wK+MO6UcarTt8cEW8YWRgRrwe2NC8kSZJU\nj/GuxN8KfCkiTub3SbsX2Bf482YHJkmSxjfeVKS3A0+JiGcAjymLL8zMTTMSmSRJGteEg71k5iXA\nJTMQiyRJmoS6ZzGTpHbkcKmazUzikirN4VI1m9Uz7KokSWpDJnFJkirKJC5JUkWZxCVJM2JDXz99\nOwa4bPvdLF27iQ19/a0OqfLs2KaOZI9lqb1s6Otnzfqt7B7cC0D/wC7WrC+m53BioqmLzGx1DBPq\n7e3NzZs3tzoMSdIULV27if6BXfcr7+nu4nurn9mCiNpbRGzJzN6JtrM6XZLUdLeOksDHK1d9TOKS\npKab3z365Jdjlas+JnFJUtNrZu4hAAAMmElEQVStWr6YrrlzhpV1zZ3DquWLWxRRZ7BjmySp6YY6\nr51+/lXsHtxLT3cXq5YvtlPbNJnEJUkzYsWSHs67fAcAnz/1uBZH0xlmvDo9Ig6JiEsi4pqI+HFE\nnDbTMUiS1AlacSV+L/C2zLwiIuYBWyLiosy8pgWxSJJUWTN+JZ6Zt2XmFeXzncC1gI0ikiRNUkt7\np0fEQmAJcFkr45AkqYpa1rEtIh4CfBF4a2b+epT1pwCnACxYsGCGo1OrjBwudYjDpUrS/bUkiUfE\nXIoEfm5mrh9tm8z8CPARKIZdncHw1EIrly1i5bJFnHjOpYA9WCVpPK3onR7Ax4BrM/N9M318SZI6\nRSvaxJcCrwCeGRFXlo/ntiAOSZIqbcar0zPzu0DM9HElSeo0jp0uSVJFOeyqJKnpRt55snD1hYB3\nnkyXSVyS1HRDd56osaxOlySpokziajsb+vrp2zHAZdvvZunaTWzo6291SJLUlkziaisb+vpZs34r\nuwf3AtA/sIs167eayCVpFCZxtZV1G7exa8/gsLJdewZZt3FbiyKSpPZlEldbuXVg16TKJWk2M4mr\nrczv7ppUuSTNZiZxtZVVyxfTNXfOsLKuuXNYtXxxiyKSpPblfeJqKyuW9ABw+vlXsXtwLz3dXaxa\nvvi+cknS75nE1XZWLOnhvMt3AE5FKknjsTpdkqSKMolLklRRJnFJkirKJC5JUkWZxCVJqiiTuCRJ\nFWUSlySporxPXA1x5kXXcfbF19+v/LTjj2TlskVT3s/C1RdOaT+SNBtEZrY6hgn19vbm5s2bWx2G\n6nDiOZcCDtIiSdMREVsys3ei7axOlySpoqxOlyRVRqOa7jqFSVySVBkrly1i5bJFNt2VrE6XJKmi\nTOJqmA19/fTtGOCy7XezdO0mNvT1tzokSepoJnE1xIa+ftas38ruwb0A9A/sYs36rSZySWoik7ga\nYt3GbezaMzisbNeeQdZt3NaiiCSp85nE1RC3DuyaVLkkTZVNd79nEldDzO/umlS5JE2FTXfDmcTV\nEKuWL6Zr7pxhZV1z57Bq+eIWRSSpE9l0N5z3iashVizpAeD0869i9+Beerq7WLV88X3lktQINt0N\nZxJXw6xY0sN5l+8AHIBBUnPM7+6if5SEPVub7qxOlyRVhk13w3klLkmqDJvuhjOJS5Iqxaa73zOJ\ny1mBJKmiTOJqyKxAI38ILFx9IeAPAUlqJpO4GmLoh4AkaeaYxCVJlWGt33CRma2OYUK9vb25efPm\nVofR0Tb09dvbU9KsUIV+QBGxJTN7J9zOJK6hsYhrhzLsmjuHM0442kQuqWNNpx9Qs9WbxK1O17hj\nEZvEJWl07XBFbxKXYxFLmnWGpjPdPbiXpWs3TakJsRF39kyXSbzCGvUr0LGIJc0mY01nCkw6kTfi\nx8B0OHZ6ha1ctogb1z6PYw97KMce9lBuXPs8blz7vElX4zgWsaTZpFHTmbbD3OYm8Yob+hV42fa7\nWbp205S+PCuW9HDGCUez75zi69DT3WWnNkkdq1FNiO0wt7nV6RXWyCohxyKWNFs0qgmxHfoTeSVe\nYe3wK1CSqqZRTYhjJf2Z7E/klXiFNepXoCMgSZpNGjWd6arli0cdY2Mm+xPNqiS+oa+fdRu3cevA\nLua3eFSyRsTSqCohxz2XNNs0ogmxHeY2nzUjtp36qc1svOb2+5Uvf/QjOOeVEw6KM8x0E3CjYnGk\nNUmavEYP0tKM+8QdsW2Eq2/99aTKxzIyAfcP7OKtn7+Sr229re4E3KhY2uFXoCRVTSfVPrYkiUfE\ns4GzgTnARzNzbbOP2aj240Yk4Eb2aLRXuSS1Rjv0J5rxJB4Rc4APAsuAW4AfRsQFmXlNM4/bTrcU\nOEKaJFVfO1zRt+IWsycBP83MGzJzN/A54EXNPmg73VLwmPn7T6p8LGdedB0LV1/IZdvv5rLtd7Nw\n9YUsXH0hZ1503aT2I0mqplZUp/cAN9cs3wIcO3KjiDgFOAVgwYIF0z7o9rvuGfWe6u133TOp/Txm\n/v6jXkVPJgGf88rehvROb4dfgZKk1pnx3ukR8WLg2Zn5+nL5FcCxmfnmsV7TbvOJt9OtapKkztPO\nvdP7gUNqlh9VllXGiiU9Jm1JUsu1ok38h8CREXFYROwLvBS4oAVxSJJUaTN+JZ6Z90bEm4GNFLeY\nfTwzfzzTcUiSVHUtuU88M78KfLUVx5YkqVM4i5kkSRVlEpckqaJM4pIkVZRJXJKkijKJS5JUUSZx\nSZIqyiQuSVJFmcQlSaook7gkSRU147OYTUVE3Anc1MBdHgTc1cD9aTjPb/N5jpvL89tcnt+JHZqZ\nB0+0USWSeKNFxOZ6pnjT1Hh+m89z3Fye3+by/DaO1emSJFWUSVySpIqarUn8I60OoMN5fpvPc9xc\nnt/m8vw2yKxsE5ckqRPM1itxSZIqb9Yl8Yh4dkRsi4ifRsTqVsfTaSLixojYGhFXRsTmVsdTdRHx\n8Yi4IyKuril7aERcFBHXl/8e2MoYq26Mc/yuiOgvv8dXRsRzWxljVUXEIRFxSURcExE/jojTynK/\nww0yq5J4RMwBPgg8B3g0cFJEPLq1UXWkZ2TmMd5C0hCfAJ49omw1cHFmHglcXC5r6j7B/c8xwJnl\n9/iYzPzqDMfUKe4F3paZjwaeDLyp/Jvrd7hBZlUSB54E/DQzb8jM3cDngBe1OCZpTJn5beDuEcUv\nAj5ZPv8ksGJGg+owY5xjNUBm3paZV5TPdwLXAj34HW6Y2ZbEe4Cba5ZvKcvUOAl8IyK2RMQprQ6m\nQz0iM28rn/8ceEQrg+lgb46Iq8rqdqt7pykiFgJLgMvwO9wwsy2Jq/mempmPp2iyeFNEPK3VAXWy\nLG4v8RaTxvsw8IfAMcBtwL+0Npxqi4iHAF8E3pqZv65d53d4emZbEu8HDqlZflRZpgbJzP7y3zuA\nL1E0Yaixbo+IRwKU/97R4ng6TmbenpmDmbkX+Hf8Hk9ZRMylSODnZub6stjvcIPMtiT+Q+DIiDgs\nIvYFXgpc0OKYOkZEPDgi5g09B/4MuHr8V2kKLgBeVT5/FfDlFsbSkYYSTOnP8Xs8JRERwMeAazPz\nfTWr/A43yKwb7KW8VeQsYA7w8cx8b4tD6hgRcTjF1TfAA4DPen6nJyLOA55OMevT7cA7gQ3AfwIL\nKGb3+4vMtGPWFI1xjp9OUZWewI3AqTVtuKpTRDwV+A6wFdhbFr+dol3c73ADzLokLklSp5ht1emS\nJHUMk7gkSRVlEpckqaJM4pIkVZRJXJKkijKJS20gIn4zYvnVEfGBGY7hJRFxbURcMqJ8YUS8rBmx\nlfv2HmxpikziUgeLiAdMYvPXAW/IzGeMKF8IvOz+m0tqNZO41ObKq9VN5WQcF0fEgrL8ExHx4prt\nflP++/SI+E5EXABcM8r+TirnfL86Iv6xLPs74KnAxyJi3YiXrAX+pJxXe2VZNj8ivl7OB/1PNfv+\ns4i4NCKuiIgvlGNmjzz+EyLiRxHxI+BNI97nd8rXXhERTynLPxURK2q2OzcinH1QwiQutYuuMkle\nGRFXAu+pWfd+4JOZ+VjgXOBf69jf44HTMnNRbWFEzAf+EXgmxYhkT4yIFZn5HmAzcHJmrhqxr9XA\nd8p5tc8sy44BTgSOBk6MiEMi4iDgHcCzyklwNgN/NUps/wG8JTMfN6L8DmBZ+doTa97nx4BXl/Ef\nADwFuLCOcyB1vMlUtUlqnl2ZeczQQkS8GugtF48DTiiffxr4JyZ2eWZuH6X8icA3M/PO8jjnAk+j\nGMp1Mi7OzF+V+7gGOBToBh4NfK8YMpt9gUtrXxQR3UB3OYf30Pt5Tvl8LvCBiDgGGAQWAWTmtyLi\nQxFxMPC/gC9m5r2TjFfqSCZxqbrupaxNi4h9KJLmkHuafOz/qXk+SPG3JICLMvOkKe5zJcXY5Y+j\neF+/q1n3KeDlFJMWvWaK+5c6jtXpUvv7PkXyAjiZYkIJKCbmeEL5/IUUV7ITuRz404g4KCLmACcB\n35rgNTuBeXXs+wfA0og4Au6b1W5YdX5mDgAD5cQYULyfIQcAt5XTf76CYpKiIZ8A3lru437t/NJs\nZRKX2t9bgNdExFUUye20svzfKRLyjyiq3Ce8+i5n4loNXAL8CNiSmRNNA3kVMFh2Rls51kZlFf2r\ngfPKWC8F/miUTV8DfLBs+4+a8g8Bryrfzx/Vvp/MvB24lqI9XVLJWcwktb2I2I9iOsvHD7XFS/JK\nXFKbi4hnUVyFv98ELg3nlbgkSRXllbgkSRVlEpckqaJM4pIkVZRJXJKkijKJS5JUUSZxSZIq6v8D\nRlBR2YkXhUEAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 576x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "T3Gkgs7UBCTm",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JLGtP7W2b5nf",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}