15-HMMs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c88549bb",
   "metadata": {},
   "source": [
    "# Hidden Markov Models\n",
    "\n",
    "In this notebook, we'll implement and assess a simple HMM. To do so, we'll use the `hmmlearn` package. There is more information in [the hmmlearn tutorial](https://hmmlearn.readthedocs.io/en/latest/tutorial.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5b6b8f28",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the necessary libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from hmmlearn import hmm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f977d2f",
   "metadata": {},
   "source": [
    "## Example #1: Weather prediction\n",
    "\n",
    "### Step 1: Define our model\n",
    "\n",
    "In this example, The state space is defined as a state which is a list of two possible weather conditions: “Sunny” and “Rainy”. The observation space is defined as observations which is a list of two possible observations: “Dry” and “Wet”. The number of hidden states and the number of observations are defined as constants.\n",
    "\n",
    "The start probabilities, transition probabilities, and emission probabilities are defined as arrays. **The start probabilities** represent the probabilities of starting in each of the hidden states, the **transition probabilities** represent the probabilities of transitioning from one hidden state to another, and the **emission probabilities** represent the probabilities of observing each of the outputs given a hidden state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "59e62046",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of hidden states : 2\n",
      "Number of observations : 2\n"
     ]
    }
   ],
   "source": [
    "# Define states & observations\n",
    "states = ['Sunny', 'Rainy']\n",
    "\n",
    "# Define the observation space\n",
    "observations = ['Dry', 'Wet']\n",
    "\n",
    "# Define the initial state distribution\n",
    "state_probability = np.array([0.6, 0.4])\n",
    "\n",
    "# Define the state transition probabilities\n",
    "transition_probability = np.array([[0.7, 0.3],[0.3, 0.7]])\n",
    "\n",
    "# Define the observation likelihoods\n",
    "emission_probability= np.array([[0.9, 0.1],[0.2, 0.8]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63fe7ae6",
   "metadata": {},
   "source": [
    "### Step 2: Train the model\n",
    "The HMM model is defined using the `hmm.CategoricalHMM` class from the hmmlearn library. An instance of the CategoricalHMM class is created with the number of hidden states set to `n_states` and the parameters of the model are set using the `startprob_`, `transmat_`, and `emissionprob_` attributes to the state probabilities, transition probabilities, and emission probabilities respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "732eddf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = hmm.CategoricalHMM(n_components=len(states))\n",
    "model.startprob_ = state_probability\n",
    "model.transmat_ = transition_probability\n",
    "model.emissionprob_ = emission_probability"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38a6ba40",
   "metadata": {},
   "source": [
    "### Step 3: Decode\n",
    "\n",
    "With our model now setup, we can create a list of observations and decode it.\n",
    "\n",
    "A sequence of observations is defined as a one-dimensional NumPy array. The observed data is defined as `observations_sequence` which is a sequence of integers, representing the corresponding observation in the observations list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "568c83f8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0],\n",
       "       [1],\n",
       "       [0],\n",
       "       [1],\n",
       "       [0],\n",
       "       [0]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Define the sequence of observations\n",
    "observations_sequence = np.array([0, 1, 0, 1, 0, 0]).reshape(-1, 1)\n",
    "observations_sequence"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5291caad",
   "metadata": {},
   "source": [
    "Below, the Viterbi algorithm is used to calculate the most likely sequence of hidden states that generated the observations using the decode method of the model. The method returns the log probability of the most likely sequence of hidden states and the sequence of hidden states itself. Log probability isn't particularly useful in itself but can be used to compare models -- like non-log probability, higher is better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "84c82493",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Most likely hidden states: [0 1 1 1 0 0]\n",
      "Log Probability : -6.360602626270058\n"
     ]
    }
   ],
   "source": [
    "# Predict the most likely sequence of hidden states\n",
    "hidden_states = model.predict(observations_sequence)\n",
    "print(\"Most likely hidden states:\", hidden_states)\n",
    "\n",
    "# Decode\n",
    "log_probability, hidden_states = model.decode(observations_sequence, \n",
    "                                              lengths = len(observations_sequence), \n",
    "                                              algorithm ='viterbi' )\n",
    "\n",
    "print('Log Probability :',log_probability)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9c1d7a8",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\"><b>Tasks</b>: Modify the code above to match the transition probabilities described in the lecture slides.    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56bec70a",
   "metadata": {},
   "source": [
    "## Example #2: Dishonest casino\n",
    "\n",
    "We suspect that a casino is trading out a fair die (singular or dice) for a loaded die. We want to figure out 1) when the loaded die was used (i.e. the most likely path) 2) how often the loaded die is used (i.e. the transition probabilities) and 3) the probabilities for each outcome of a roll for the loaded die (i.e. the emission probabilities).\n",
    "\n",
    "First, we'll act as the casino and exchange a fair die for a loaded one and generate a series of rolls that someone at the casino would observe. If we have a fair die, we'll have a 95% chance of continuing to use it, but a 5% chance to switch. In the loaded die state, we have a 90% chance of staying and a 10% chance of leaving."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "8a60fd97",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3gAAAFNCAYAAABSRs15AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAB6+klEQVR4nO2defwdZXX/PycbCVsQCAgkMVEBRRDFiFAVwa2gFmrVCnWt1rihtuKutdTlZxVb64JVXFERUHGJSkFlETcoQQEFBAOyBATCFtaELOf3x8zcO3fuLM82c+/MfN68wvfeO89znjPLnfucec4iqgpCCCGEEEIIIe1nxqQVIIQQQgghhBASBhp4hBBCCCGEENIRaOARQgghhBBCSEeggUcIIYQQQgghHYEGHiGEEEIIIYR0BBp4hBBCCCGEENIRaOARQsYQkZeIyE8mrQchhBA7RORaEXlm/PpYEflGQ+N+VUQ+5Ni3VE8RuUxEDs62FZHFInKviMws6XuviDzcRa8SmSoijwwpk5CQ0MAjnUdEniIivxaRtSJyh4j8SkSeGG97pYj80kLWkvjGPqs+jZslb59U9SRVffYk9SKEkL4SG2kPxMbJzbHxtPWk9ZoUqvoYVT035/PrVXVrVd0EACJyroj8U6bN1qp6TUOqBsPWiMzbd9JfaOCRTiMi2wL4EYBPA9gewG4A/h3A+knqNSm6ZJgSQkjH+RtV3RrA4wA8HsC7J6lM2SoZIWS6oIFHus4eAKCqJ6vqJlV9QFV/oqqXisijAXwOwIHxU9K7AEBEnisivxORu0XkBhE5NiXvvPjvXXGfA+M+rxKRK0TkThE5U0QeFn8uIvIJEbk1lvd7Edk7T1ER2VVEVsSrjKtE5DWpzx8Qke1TbR8vIreJyOyy8eNtKiJvFJE/AfhTztBj+5Rd2YxlvEFE/iQi94jIB0XkEfHK6N0i8i0RmZNq/zwRuVhE7orbPNbgXBFCCMmgqjcDOBORoQcAEJHDY7fFu+KVm0fbyhWRg0VktYi8J/49uVZEXpLa/lUR+R8ROV1E7gNwiIg8Oh7vrnj8wzNidxSRn8a/Ez/P/BZ9Mv5NvVtELhKRp2b6zhWRU+O+vxWRfVN9B26nmX0YeKCIyIcBPBXAZ+Lfss/EbQYrYSKyhYh8XESuF5FbRORzIjIv3rajiPwo3rc7ROQXIlI2T36OiFwTH7vjRGSGiMyJ++6T0nEnEblfRBbk6P/I+DitjeWcGn+e/C5fEu/Li0XkIbF+a+Lf+h+JyMK4fdG+Pyo+H3eIyJUi8vcl+0M6BA080nWuArBJRE4UkcNE5CHJBlW9AsDrAPwmduHYLt50H4CXA9gOwHMBvF5E/jbedlD8d7u4z29E5AgA7wHwdwAWAPgFgJPjds+O++wBYD6Avwdwe4GupwBYDWBXAC8E8P9E5OmqehOA3wB4QartPwD4jqpuqBg/4W8BPAnAXjnjju1TgX5/DeAJAA4A8A4AJwB4KYBFAPYGcBQQGZ8AvgzgtQB2APB5ACtEZIsCuYQQQgqIJ/GHAVgVv98D0T3+nxHd808H8MP0QzYLHgpgR0TeLa8AcIKI7Jna/g8APgxgGwAXAPghgJ8A2AnAmwCclGn/EgAfjGVeDOCk1LYLERmp2wP4JoBvi8jc1PYjAHw7tf37Ej/ENEFV34vo9+/o+Lfs6Jxm/4Ho9/hxAB4Z7/f7423HIPoNXgBgZ0S/q1oy5PMBLAOwX6z7q1T1QUS/5S9NtTsKwFmquiZHxgcRHc+HAFiIyNsIqpr8Lu8b78upiObsXwHwMACLATwA4DNF+y4iWwH4KaJjuROAIwF8VkTy5gGkY9DAI51GVe8G8BREN+kvAFgj0SrZziV9zlXV36vqZlW9FNEP6dNKhnkdgI+o6hWquhHA/wPwuPjJ5QZEP4yPAiBxm79kBYjIIgBPBvBOVV2nqhcD+CIiQxOIbtCJASWIbtTfNBg/4SOqeoeqPlCyH1V8TFXvVtXLAPwBwE9U9RpVXQvgfxG5EAHAcgCfV9UL4lXTExG5xB7gMTYhhPSN74vIPQBuAHArgH+LP38xgB+r6k9VdQOAjwOYB+CvHMf5V1Vdr6o/B/BjRA8iE36gqr9S1c2IjKKtAfyHqj6oqmcjCoE4KtX+x6p6nqquB/BeRB4yiwBAVb+hqrer6kZV/U8AWwBIG4cXqep34n36LwBzEfB3I/7tXA7gX+Lfw3sQ/V4eGTfZAGAXAA9T1Q2q+gtVLTPwPhrLuR7Af2N4HE4EcFQ8HgC8DMDXC2RsQGSw7Rr/9hfmBIiP3Wmqen+s+4dRPjd5HoBrVfUr8TH/HYDTALyopA/pCDTwSOeJDZ9XqupCRCtNuyK6GeciIk8SkXNiN4i1iAyoHUuGeBiAT8ZuHXcBuAOAANgt/gH8DIDjAdwqIidIFBeYZVcAyQ9OwnWIni4C0U35QBHZBdGK22ZET+tKx0/JuqFEf1NuSb1+IOd9kgDgYQCOSfSJdVqEaB8JIYSY8bequg2AgxE9JEx+h3ZF9PsAAIiNrxswes835U5VvS/1/jqM3qvTvx27ArghHi/dPve3RlXvRfR7tCsAiMjbJAolWBv/LszH6G9ruu9mDD1aQrEAwJYALkr9Np0Rfw4AxyFaJf1J7Hr5rgp56WMzOG6qegGA+wEcLCKPQrRSuKJAxjsQ/V7/n0Qur68qGkxEthSRz4vIdSJyN6Lwiu2kODbyYQCelPktfgmiVVvScWjgkV6hqn8E8FVEhh6Q737xTUQ340WqOh9RnJ6UtL8BwGtVdbvUv3mq+ut4zE+p6hMQuUfuAeDtOTJuArC9iGyT+mwxgBtjGXcicuN4MSKXmVNSTxZLxy/R22SbCzcA+HBGny1VNes2SgghpIJ4Ze2riFbqgOj3Ih3bJogeot3oIP4hsStfwuJY/mD41OubACzKxKUNfqdiFqX02hqRu+VNEsXbvQPR6uBD4pCItRj+tmb7zkDkspjWxYSy37PbED2MfEzqt2l+nMgGqnqPqh6jqg8HcDiAt4rIM0rkLUq9zh63ExG5ab4MUTjFulxlVW9W1deo6q6Iwho+K8WZM49BtOL5JFXdFsPwiqL5yQ0Afp75Ld5aVV9fsk+kI9DAI50mDjA+JhWIvAiRG8X5cZNbACzMxC5sg2g1bZ2I7I/IoEpYg2j1LF1T53MA3i0ij4nHmC8iL4pfPzFeEZyNKLZvXdx/BFW9AcCvAXxEROZKlJTk1QDSdYG+ichl84UYumeWjm9I3j758AUAr4v3W0RkK4kS12xT2ZMQQkge/w3gWRIlHvkWgOeKyDPi35ZjELnB/7qkfxn/LlFykKcicuv7dkG7ZGXqHSIyW6K6dH+DKOYs4TkSlSaagyi+7Pz4920bABsR/d7MEpH3A8h6szxBRP5OomzP/xzv0/mw4xYU/JbFq4JfAPAJEdkJAERkNxH56/j18yRKeiKIjM9NyPm9TvF2iRKfLALwFgCnprZ9A1GM3ksBfK1IgIi8KJmfALgTkZGWjJndl20QGah3SZR07d8wSrb9jwDsISIvi8/X7HhOYp2Qh7QPGnik69yDKLnIBRJlATsfUfzYMfH2swFcBuBmEbkt/uwNAD4Qxz68H9GPKQBAVe9H5Pf+q9jl4QBV/R6AjwI4JXab+AOigHgg+gH7AqIb93WIEqwcV6DrUQCWIHoK+D0A/6aqP0ttXwFgdwA3q+olKZ3Kxq8kb59M+xbIWwngNYhcU+9E5PLySh+ZhBDSZzRK0PE1AO9X1SsRGQ6fRrQq9TeISio86CD6ZkT36ZsQJUR5XezpkqfDg/FYh8XjfhbAyzPtv4nI8LgDUVKuJNnImYjcIa9C9Fu4DuOhAz9A5KVyJ6KVr7+L4/Fs+CSAF0qUZfJTOdvfieg36fz49/JnGMYB7h6/vxdRYrPPquo5JWP9AMBFiJLJ/BjAl5INsVH7W0QG2y/yOsc8EdH85F5Ev/Fv0WHNvmMBnBj/Lv89IiN/HqJjfz6i41m473HIx7MRxRjehOhcfxRR7CPpOFIeP0oIIYQQQrpGvAL3jTg+nQRGRL4M4CZVfd+kdSH9g0WPCSGEEEIICYSILEFUuujxFU0JqQW6aBJCCCGEEBIAEfkgolCJ41T1z5PWh/QTumgSQgghhBBCSEfgCh4hhBBCCCGEdAQaeIQQQgghhBDSEVqXZGXHHXfUJUuWTFoNQgghDXDRRRfdpqoLJq1HW+BvJCGE9IOy38fWGXhLlizBypUrJ60GIYSQBhCR6yatQ5vgbyQhhPSDst9HumgSQgghhBBCSEeggUcIIYQQQgghHYEGHiGEEEIIIYR0hNbF4BFCCCHEnA0bNmD16tVYt27dpFWZCubOnYuFCxdi9uzZk1aFEEJqgQYeIYQQ0jAiciiATwKYCeCLqvofme2vBHAcgBvjjz6jql90GWv16tXYZpttsGTJEoiIh9btR1Vx++23Y/Xq1Vi6dOmk1SGEkFqgiyYhhBDSICIyE8DxAA4DsBeAo0Rkr5ymp6rq4+J/TsYdAKxbtw477LBD7407ABAR7LDDDlzNJIR0Ghp4hBBCSLPsD2CVql6jqg8COAXAEXUOSONuCI8FIaTr1GbgiciXReRWEflDwXYRkU+JyCoRuVRE9qtLF0IIIWSK2A3ADan3q+PPsrwg/n38jogsaka1evjUpz6FRz/60XjJS16Su33lypV485vf3LBWhBDSTeqMwfsqgM8A+FrB9sMA7B7/exKA/4n/EkIIIX3nhwBOVtX1IvJaACcCeHpeQxFZDmA5ACxevLhS8Cd+elVANYF/edYelW0++9nP4mc/+xkWLlyYu33ZsmVYtmzZ2OcbN27ErFlMF0AIITbUdtdU1fNEZElJkyMAfE1VFcD5IrKdiOyiqn+pS6c8LrruDlx1y72lbQTAIY/aCTtvO7dS3uU33Y1LVt8FAJg5Q/DXj3ko5s+rN1PXDXfcjwc3bcYjFmxt3ffXq27DE5duj9kzp99b9677H8RPLrsFm1St+wqApz9qJ+xkcA4B4O51G3DGH27Gps3lY+0yfy5EBDfd9QD22W0+9t5t/sj2n15+C267d31u33mzZ+KwfR6KLWbNLB3jjzffjd9df5eR3mn22mVb7LtoOwDA2gc24MzLqvenjCU7bIUDH7GDcft1Gzbh9N//Bes3bi5s84SHPQR77LyNs055rLlnPc7+4y3w2NWJM3OGYO9d5w/uJWlmz5yBw/Z+KLbaYhZuvOsBnHfVGr+xRPCsvXbGQ7aaU9jm11ffhutuv3/kszkzZ+A5++yCeXPKr1+Sy40A0ityCzFMpgIAUNXbU2+/COBjRcJU9QQAJwDAsmXLpu7Kf93rXodrrrkGhx12GF760pfi+9//PtatW4d58+bhK1/5Cvbcc0+ce+65+PjHP44f/ehHOPbYY3H11VfjmmuuweLFi3HyySdPehcIIaRVTPKxWJGLypiBZ/t00oY3nvQ73Hx3dbD1K/9qCY49/DGV7d713Utx6eq1g/f3rtuIVz2l3kxdT/3YOQCAa//juVb9fnf9nfiHL16A5Qc9HO95zqPrUC0op154Az7yv3907v9PT1mK9z0vL4/BON+9aDWO/eHlVvL32mVbnP6Wpw7e33L3OrzmaytL+8zf8ok4ZM+dStu893t/wEXX3WmlCwAs3XErnPO2gwEA3155Az704yusZaTZcs5MXP6BQ43bn3XFrXjrty4pbXPAw7fHKcsP9NIry5d/9Wf8z7lXB5U5bagqXrRsET75s6vwrZWrveW9/d498cZDHlm4/R+/cmGuoT5n1gz8zb67eo/fQy4EsLuILEVk2B0J4B/SDTIPPA8H4PcFniCf+9zncMYZZ+Ccc87BnDlzcMwxx2DWrFn42c9+hve85z047bTTxvpcfvnl+OUvf4l58+ZNQGNCCGk3rfB7qPPp5PqNm/CC/Rbi7X+9Z2Gb537qF6WrECPyNmzGQXsswAcOfwwO/vi5eHCTWb9JcPu9DwIArr61fAVzWkjOwS/feQhmzbBbcXzWJ35udS6Ssc4+5mnYck7+1+Tk/7senzzrT4P3WfkPxjLe99xH43mPHZ0EX3XLPXj5l/9v0KZcl0148iN3wH++6HHG+n/wR5fj4hvuSsmIxjnv7Ydgziz71dr/OXcVTrrgeqs+6zduAgB8+3UHYtFDthzb/qaTf2u0/7as37AZW86ZibOPOTi47CZYv3ETnnbcuQCARdvPw7df+1eDbbfdux7P+/QvB9fa+o2bsdt283Da6/8qT5QRB3zkrMr72/qNm/GqJy/F8oMeDgBYfef9eOHnflPL+esDqrpRRI4GcCaiMglfVtXLROQDAFaq6goAbxaRwwFsBHAHgFdOTOGArF27Fq94xSvwpz/9CSKCDRs25LY7/PDDadwR0kfO+cikNWiGQ95dq/hJGniVLipNsdUWM/HQ+cWuezNm2GXcmjd7hpE7J3Fjl/nzMNPynNi2T3jo/LmFBt62hq632205Z+z6uv2+fLfNIubOKr9GsxS5ze08f4tKl9A8ttrC/VaxYOstcnXfYtZM3P/gRme5ZcwUsTpe08S6DZsGr2fPmDGyH3mX8ayZzezr1nNnDcZJjHfijqqeDuD0zGfvT71+N4B6ZwAT4F//9V9xyCGH4Hvf+x6uvfZaHHzwwbntttpqq2YVI4SQDjHJwKsVAF4eZ9M8AMDapuPvAEARxWeVIYOWJvIUAkGShdkhXIwUkBxLF1NNYHcukqZSMlp2i2YGKNM3kWuikypgm9VbcvRJj2uLiOk3YMhg/wuGdJFpNC7U7SKZEkaOV3Y/MvcVDbWrBhfiqFo2d0VChqxduxa77RYlDP3qV786WWUIIaSj1Fkm4WQAvwGwp4isFpFXi8jrROR1cZPTAVwDYBWALwB4Q126lBFNnsunSCLmxkF2Mq6cAgUjOZYuJYxExOpcVBkneduy0sv0HX5WrVOki91OZ42nxNhzLf8kkFyDsQwTI7mOByDBjJ4JkT5e4/bdqGGl8K/pVWVo5107wwdYvL8RO97xjnfg3e9+Nx7/+Mdj48Z6VvAJIaTv1JlF86iK7QrgjXWNb4rpBMXYwHPsR6rxPZZ2K3gOg2W6mK7OWYo1Ji07xHVov4JnYLy6qULiYxvKwCoTU7otyOikaUzKGoTm2muvBQDsuOOOuOqqYZmGD33oQwCAgw8+eOCueeyxxzasHSGEdIvpz40/Bdi6tYm4r5SQalxWLOo4HaYyy1fw3GWUti/QzvU4+FzPxS6a9X1J6pRdN6MrZVK4bfCZ73jG7dp7TAkhhJA+0XsDL3JxKm8TuTCZrvTFMXicDAXHZ7XANt7LzEVzdOO4i2bcLudasIlhiq4pO7LX7EAXR8PHNoYxPWaZzDqWuFW11Q9YDELwRlw0fW81Ve7Lw2sn3SezkRBCCCFTQ+8NPNMJiruLJmdAwfA8lnWfivEkK6bxdfUQ2kXTXoEgTeoYtvXoiIUXUF7utjLjrw9HmxBCCGkXNPBQ7Xpk/YCcLpq14X5cXVeuSrJoGoqciItmkVuknRh3BQy61vkdafPXL73Kmj1GeSuwzblolutBCCGEkOmg9waemYummD+njjP4DTyY+IA7GCYlLYqIzrFNFs3qrJNjZRIq3uf1NVkBibJCWjtpjmbR9MhAGklLdLE4hoMxi+MBa8ui2WIDJK8cQXZbch4UWn8WzVS7cT28hiaEEEJIDdDAM4xvsnHRFJHBpIvzn3D4TNxtjQmjmnvZGLyCLJq5qy6ZemalusA+pixb2qNMF1N5aTkmVB1D29IVxuPCPmZxmsiNdcu8H3hoBigJEZXAKN6ed+1k9SCEEELI9EADz7ideZKV0feWCpFCfI0Bl0LnVvKt1vDMx3G/hsaTrPhilajGpE1NK3hdJ13oPIi80iQrJdt6cKxJGLbeeusgco499lh8/OMfn8jYhBDSFmqrg9cmTLJoWslDu2OAppmm0/yXrXi5xC7Z9070sGoe/PrzyQobPB7QY8w2UH7NFa8Guw/o0qXFB5gA53wkrLxD3h1WHiGEEC+4gmfg9icC82ybGK2Dxyxz4YjOlVtfgZ07oImLZlaXYhfN4r6mxcBdajGOiPYsHWCjb2rIqG+RiZvVMRBx5F94wRNg7N6UddGEehtbVeGpedcx72/EFVXF29/+duy9997YZ599cOqppwIA7r33XjzjGc/Afvvth3322Qc/+MEPBn0+/OEPY4899sBTnvIUXHnllYPPr776ahx66KF4whOegKc+9an44x//CAD485//jAMPPBD77LMP3ve+9zW7g4QQMgX0fgXP2PXSVF7BJJ/443so7Vw07UcbO/fWEorkuknSgtc+2LlohnFRtaUP37lBkpVgLpqueoQZn/SH7373u7j44otxySWX4LbbbsMTn/hEHHTQQViwYAG+973vYdttt8Vtt92GAw44AIcffjh++9vf4pRTTsHFF1+MjRs3Yr/99sMTnvAEAMDy5cvxuc99DrvvvjsuuOACvOENb8DZZ5+Nt7zlLXj961+Pl7/85Tj++OMnvMeEENI8vTfwgOpn/darJ2h3Fr9pxnW1wt1F01+X/ELntorYNg+bTt+rb2Hn+r4jbf/6Jaub2d0IUXLDtX/6mmr54SUT5Je//CWOOuoozJw5EzvvvDOe9rSn4cILL8Rhhx2G97znPTjvvPMwY8YM3Hjjjbjlllvwi1/8As9//vOx5ZZbAgAOP/xwANGK369//Wu86EUvGshev349AOBXv/oVTjvtNADAy172Mrzzne9seC8JIWSy9N7AU4Pc+9Fky3SlbzRtOR9wh8PkXBVh4WU7HAsV8VAVupS7aMpIm1I5sN/t7DXrWzqgniya7quTFSO33gBJrtcCD80BIY5elEWzOpHKiC7MokkCc9JJJ2HNmjW46KKLMHv2bCxZsgTr1q0rbL9582Zst912uPjii3O38yErIaTPMAYP1SsxNsZB+qm71FXoq6f4pL8XKU8FPz6Wgcxsn2wG1ZJoMMm0qVLGdrKSvWZ9SwcMy37Y1MFLKZMn00Of0nE9YjWnheR4FxU6T2fRDFIHrywGL+c6HtwzeX8jljz1qU/Fqaeeik2bNmHNmjU477zzsP/++2Pt2rXYaaedMHv2bJxzzjm47rrrAAAHHXQQvv/97+OBBx7APffcgx/+8IcAgG233RZLly7Ft7/9bQDR/feSSy4BADz5yU/GKaecAiAyHAkhpG/03sAzTp7i2I7Tn4B4HkyruDqHietYkQQT281wBc+FvDp4vljJMUkgwzIJTgyv5TA7Wyal3PgjxI7nP//5eOxjH4t9990XT3/60/Gxj30MD33oQ/GSl7wEK1euxD777IOvfe1reNSjHgUA2G+//fDiF78Y++67Lw477DA88YlPHMg66aST8KUvfQn77rsvHvOYxwwSs3zyk5/E8ccfj3322Qc33njjRPaTEEImSe9dNAGTMgluOepbvogwlTS5MhOqfEaIuCl7F80a0uk7UpZFs+kx28LACyCzH2Wrwb5jVbbLyaJJWsoEyhrce++9AKJ703HHHYfjjjtuZPuOO+6I3/zmN7l93/ve9+K9733v2OdLly7FGWeckft5WtaHPvQhH9UJIaR19H4Fz8R1zTZ+K5mU2boFknJcygUk2JS6GI5VITOddCLHzW34Pi/hiUUMnmOJg5EYPPgZPU7jV/S1LV1hPq5fSYhpINF/3EUz+jvqouk7Vvl9anAec5Ks8P5GCCGETB808EwmSDZJVjKTcdaJCoeroQPEBpjVWAYrtzL6MnuuB7FLJSt4Jjq55pbJlMHzWuqxMUhHxiwZtrY6eB34yg0eEhV8nnbQ9DbwUH6fSu59oyt4MrKNEEIIIdND7w28Omn5IgKxYNoytk2ZOo3T890nhBBCSI/pvYEXPIsmRrNo8gF3OPLqgplSlQp+bCxD192EGSUummVyTHRyyZQoGF2y9M+iOZRjynDlpzgGr54iCdNncFszkoo39fGYi6b6xxtWZtEsVo+3t/bA1dYhPBaEkK5DA8/A7U8sZqJpl88oxoiEwmfi7uaiWSUzHZNUfK7LEp6YuWjaG2dj++sZq+USczWM3SqSaWd0G4/bgS+dZP5mSQztUC6aZeTVhHSpi0gmx9y5c3H77bfTsEH0m3/77bdj7ty5k1aFEEJqo/dZNM1X5gxj8LJxWPw9DYbvsQxdB2+E0iQrfgO57nc2yUoIAldJqGkFr/tfunSSlTDy3GohdP9Id4OFCxdi9erVWLNmzaRVmQrmzp2LhQsXTloNQgipjd4beIBJtkRbeayTUBfuLpoufapddxNmlDTNTW1fc52E/HT6zWbRrOxbZ5mEln/3qrJojnzmPZaZBBl53fID3DNmz56NpUuXTloNQgghDUEXTQMfJ5tYulEXzX6sJjSFwj0IT8TOXdYk62T6spmRU4ehNItmpk2VLrYT6qyLpk8GUiCdRdMmjnG077hM1LWE134DryKLZoIC3jtb5b6cex0PXDR5fyOEEEKmjd4beIBZvTNjAw+ZiRDnP8HwS7Jia5wYJFnJlkkoSrJSVibBUCW3Onijr8MkWbEZvzzLjK3RbTxuDTKbZriCV5RkRQd/vVfwUJFkJec0tt2AJoQQQrpMrw082wm/mczh67qyBIZimnWrA6v9tTw4ecZKoBA8p1WSSJ8aYvACXzT1JFkJkFlyygl+Hsrq4DWoByGEEEL86bWBl1CdLdFaYvz/bk8yJ4Fz+nuHbtV1zsezCla1K/usfCw7gsT9BaJo2Dq1afsKU1EWzfzz6jmWaQzeSNZYQgghhEwrvTbwhq5H1dMV8yfVmnKvmu4YlbZN0nziyGzjvYzqI465aGYzqBa7KNq4aDqlws/Ejfq7aNr7aFbVAay1Dl4NcpskOd7jSVbiWMj4ve95BUxcNMdj8IZ6TO/9jRBCCOkr/TbwamjrlCqfGOF7KF2KdJvi7qJplmTFhZEkK4Em4lbHsMZ9a1rmtDEokxDsvDpu68GxJoQQQtpGrw28BJuC1kbyMn9JOJyTrDgs/dl0sS6TYKuLtUtn2Hz6QVb/Asp0HbMtDO8hkvv5yGcBsmja6FSkByGEEEKmg14beAPXo4p2VS5MIzIxmgGPD7jDEZWgcJta2pzDwVhVMlO6zJDxTKvDLJrFxpaZi6a9a2q2akOwLJpBXTSlFhe/LrhoFj0lGnrKJlk0Q+xreZbgwbZ0DJ5DVlVCCCGENEO/Dbz4r1mSFdMsmqMZ/OjCFA6T0gVF2NQyjMaqNiazaePH3TqLHyAM6sqZ6OJY1y1rPPms9CQ9Xdyay5Ks1OOi2f4vXWGSlST2LXHRDFDzr+r+pjnX8bAuot/YhBBCCAlPvw08h9WIynap19EiyvTOgKZXs3x8J5N2MXi20nNi8Eytt6omtqpg3HgKZfRYlRYxWZ2sw8ADOrCEV84gyUqoGDzHILxpvr8RQgghfaXXBl5C5UqNi3sc0PlJ5iRwz6LpEINnoUtpDF5JFs1w2tQxZrqvz+pf83US2v7VK8qimdvWc29Ne49m0fQakhBCCCE10msDz/Tps+SszhTKTMXE1OWCFoq2zdHMIibzsXfRrA5uGquDl43By2k37DvaplQXBze87DXrG5fmEnM1cO0rtO9qugI9YjWnhUEcb961k7qxaICAw6rvRtl1PM33N0IIIaSv9NvAs3LRNI/BI/Xg76JZ31h5DwHCuSi67Xj6Wgx1WYZ2a67j+9IHt0HN/PWXVxKD1/3DSQghhHSKXht4CWZJVmzkychfEo4mD2m9LppupTeM2wcY02d8I5l00Sxk4AVgUO3Cd19NV1LpokkIIYS0Axp4qJ7g5HjfFTKSZCU3syJxxyeLZnkq+KI+pdvH5I8OUFaGQzJtynBz0RwvdO41J0+yNzqsGZVm0fRQqYgQmSUnTVkM3si1HCiLZrmLZlkWTd7fCCGEkGmj1wbesE5ZRUMb4yAz4eL0Jxw+E3fzogTJWNW157IrGvlFEpBr4dnEtEVhVvY7PppF0zPJSloZ4/ETwyB/YNu4SPNxw8tsmmEcb3785qAOHtTp2hiXV0zefdKlLiIhhBBCmqHfBp7NhN+h3bQnWZli1XLxjsGzrINnQ56xEioGz2mVJGPNhYvVsmhrsm+1FDr3N3qmnXQdvJDycreV9QszPCGEEEIC0msDL6FyAc9aHmPw6sJ14u5yKqq7DFvMKBkgfxXGv+xBaXuLT+sY36RvnUZY2796gyyaBuUu/F00DWPw0llj/YYkhBBCSI3UauCJyKEicqWIrBKRd+VsXywi54jI70TkUhF5Tp36ZDF10bSJpUu79qVdqaaRtk3SFNVuk0XkuVCWjmWQan/EZQ3j57q0TIC1i6YdQxe6YTp9PxfNJObKvM8wvX6BTLpollD+wCCdRTOEMVueRTO5UaZ0EPvrgYxS9fuYavcCEVERWdakfoQQQtpLbQaeiMwEcDyAwwDsBeAoEdkr0+x9AL6lqo8HcCSAz9alTx61JHjIvucEKBj+Lpo2Lrl2g+UmcQnkfumz38O+YS5EK7dmIxfN8PThKzd00Qzlo1k9Vn63Phzt8Bj+PkJEtgHwFgAXNKshIYSQNlPnCt7+AFap6jWq+iCAUwAckWmjALaNX88HcFON+hRikkXTTl78t21LZC3AOYumQ08b192yc527gOdYesO4fUlxdRf8XDSLk6zURdvdo4cumtV1EryTrBh2l4LXxAmT30cA+CCAjwJY16RyhBBC2k2dBt5uAG5IvV8df5bmWAAvFZHVAE4H8KYa9RljkOWv0kXTPIvmqCvcePFr4k7kjuYeg2fvolklczQGr2gBL09nm71wWaXJZukMlUXTzkWzuExEIrUuF822GyCS+ZvdNsyiGahMglG7VAwes2j6Uvn7KCL7AVikqj9uUjFCCCHtZ9JJVo4C8FVVXQjgOQC+LiJjOonIchFZKSIr16xZE2xw07mJTSxdFCc2nAhxAhQOn2Npm9G0yjTJbpVhp6GMkhhPmxgmn0n8SAxegCQrLlk0C5Os5B20ILT/S1eZZCVgFk3BeA3HNDkheMPr1394kkP8O/hfAI4xbF/LbyQhhJB2UqeBdyOARan3C+PP0rwawLcAQFV/A2AugB2zglT1BFVdpqrLFixYEExBqwm/xQpeQn0T2DBMr2b5+Mb7hE7xnyZaBclPslKuUz2BaoMVN5txDLCLYzSR565LmcyWe2hWopm/oeTlb3MM0CNlVP0+bgNgbwDnisi1AA4AsKIo0Updv5GEEELaSZ0G3oUAdheRpSIyB1ESlRWZNtcDeAYAiMijERl4jT9+tMmWaCQv85eEw3ni7tDRptB5eZkEs89KxwpQViFEFk2nvkUxeM4STcasUXgDDEqtlGwbvPfcWeMYvJYf0ymj9PdRVdeq6o6qukRVlwA4H8DhqrpyMuoSQghpE7UZeKq6EcDRAM4EcAWibJmXicgHROTwuNkxAF4jIpcAOBnAKzVYWjgTJaM/1ck0LGLwUgLrSgMfitbN1zxWZiIXTas1vOrrIpM2vqjQeZmxVZeLZtYF1DsuzSXmqqKxbVyk8bDwTzwyacqSrIyUbdHq67RyLFRkyiy4jus6f33A8PeREEIIcWJWncJV9XREyVPSn70/9fpyAE+uU4cySuuUpclxvysROjK5nGYDb4pVy8Vn4m5rIBklWckUfi6ukpCXZMU8hkk9JvGjyTg8VuGcxq7KLloe++VKk8+I6qIyyUpi3yFEkpXyZFCDZEHZlUNM9/1t2qn6fcx8fnATOhFCCOkGk06y0mnavopAzKH72nTB80EIIYSQvtJrAy8vO1weNk+qoyyacT+blb8J0LaVDlX1dNG0GcvAQM/E4GWPZ1kZjvpdNEfl+55qm6yfCVVuofW6aLab5HgXZWAdKX/hOxbK7wVF13GkR7vuIYQQQkgf6LeBF/81SbJiOo3JTrim2YaaYtVy8Zm4205G04Z6ocwR+cUummVijDJtqn8iDZP9KWOYldP2GJYnn6kri2Z3yHPvTbtolh9j0yFMXDTL9CCEEELI9NBvA8++MJpVs7y4rGmibZMzX33rKIuRkDvZNVmdM1rBs9/xQYxf4Hz6oY9hLTF4QOd9NAexlQ3USShNwBJoeEIIIYSEo9cGXoJNMg0bed5P1skYrsfUpVd1Fk3JfV3WbvhZWF1M5Hut4Hm4xhbLrO/70fZvXlmh8+zOhXDRNGqXUYa3N0IIIWQ66bWBZ+JCB9jF0kUZD9uRRbNtz9/9XDQtV59g4Lqbej0jp+kgS2tuX/M9UYcdz7pU+pYOGMT0WfQxiR2sJQbPI1ZzWhgYeHnbkCl/ESSLZrWTZnYYm/IxhBBCCGmOfht4wyC80nY2xkF2UjvNSQjaNjlTj4wSAssYPIOm2Tp4Ub9hR7M6eAYxeHA3zoaGgG8M3vj+mYxdprdMuw/zBBkUOi9IspJQ/AjBZqwKN8yi63jKk0gRQgghfaXfBp5VwgjDdqmG0z6BnWLVcvHV124Fz260wYpZqptZDJqRMtZkV9waCNXKaWtmvIYmRGbJaScxtEPFMDrH2bXtJkIIIYT0gF4beAnVZRIsY/CSv12fZU4A50Pq0NEmNnNGabbI/EyIIXWpY0yf8U0GrbNOZNvjX4cumtXxm/4umobtcgqdE0IIIWT66LeBV+JClyZy0bR4VD3mQjWdtM1FEx7lAqwXUw1im0a258SolXkAD1w6jVRR58n0cKUnjNFjdc2Y1MGrJYtm2y7sccoeEmWvZf8kK+Xuy0UumnXVMSSEEEKIH7028GwTRlS2yUxWoyQE0zsFattE2Fvf8FUxBgxdNNMxeGY17kK0GdOnJhdNG0k1eZ9Wy+yFi+bo31DycreVGn/tuocQQgghfaDXBl5CaFcxumjWh08WTfuxbLJournxGrcPcC35uWg6rpyWumjWR9u/e8nxrkqyUtTGbizDdmPvW36QCSGEkI7SawOvLMthGhGzdOBZeVOeY6V1Lpo+KeHts2gaZJ1MbZ8Rf5PMXTSTNvVk0Rxm9YxleC5r5SWRqSJbMmRMpmXpCvNx2298DB4SFcRSDstflB9jU0qLJJS5aLbsHkIIIYT0gX4beIZJxk0NtWFdvVQM3hRPgKZYtVx8Ju5OdfCqZKZaDMsIpGToeLuhPuPtC3XxKXGgwz8hkqxYuTVXGOTV9dfcaJvrcS6S+ZveJBhx0QxSB8+gTEJWmWl/gEUIIYT0lX4beLaz1commRg8kameALUtfsZ34h74dI+QvyIXpkyAy14PVtzSFl4AbI3kkPKMx/W1ZltAk7GVJglYCCGEEDI99NrAS7DKlmghr+NzzIng7KLpVCahIgZP0q/LXRF9CRGz55NF03nl1HGbL23/7pUs4I196r2CZ9puzEWz7UeZEEII6Sa9NvDyXCrzsHfRTH3GR9zB8DmUthlNzVw0x1/numgWCBKB4cpwiRCTvvArtYDU8HZxjBVGQE1p9rvwjRsmWclz7x2NrfSNwatyXx66Gmf6oSPusIQQQkjH6LeBZzPhd0lnP+V1otpme/qqa+eiaTdanh0zqTIBQNpFM/47ERdNky+Nuy5lMru/uKSp/4eS5tCvZfcQQgghpA/02sAbYJBF00ocXTRro450/cWdzHUpK5NQtMJio5K9i2b+yo8rzuUpSrfV9w3pThbN4m3FH1iOZVomYcxH029cQgghhNRDrw28ItejLKauSIOsnGn3qgafcNuuOrXNvcq3gHXwLJrpMgmJC6OVi2Z1+Y3knLrGgSb9fUsHSM7+VVGV7CT6etSTRbPtK3iDh0QVWTTh+Z0Aqt2Xy1w0CSGEEDJ99NrASzBJpuHkoolmjShbd6n2uVe5T9ytM5pWxY8hL2n86PmuKsNh8uCgrNSCCUM7wNfoGd8/817F22rLotlykvOdv4I3fDAQhWcGiMEr2T58cJXtZxfXSgghhJBm6LWBZ70aYUldE9gibIdq29zMW1+rJCtuMXh5K3g+Krnucjbpy0Ri8EwSyHjoUiaz7St4VQwKnQc6sWZ18PL0IIQQQsi00WsDL6F6LugYg9fxSeYkcF7Bq6FPnotmVTuTz33bunWoR1zd5SMKZbfcgXB4D6mOpfR30TRsN7aC5zkwIYQQQmqh1wZeketRFtNQujx3ukZX8Kxj8NpFFNLlnmTFLoumyQRWUq8SF8aUjJx2o72r3UYHMXhVqhT1T2Vb9KuD5zJ2+TG0LV1hPG7blqZLKEqyMuqi6TtI+XVYVE6maQ8FQgghhJjRbwMvuLvZaJtoAt9gDJ5t+5bNznz1rdMlN5vUJPvadRxfF81EwGRcNA3auKtSKrPrq0vhy1+UJVkp2da6x0SEEEJI9+m1gZdQuYLnKK/rk8xJ0KyLZnXynYTSMgllWVZMdXHMouk4XI48x5XTUpluuvSBwfHOzaKZl97HYyzHhr7JXQghhBBSD7028Ipcj7KYzmPy0ok366Jp2b4eNWrDpHRBEVEWTbskKzaG/2AFr6Jd9nPjLJq2tRjHXEb9smgOkrZYH8MSwxf1ZdFsu/Ehg7/5+zF00Sw/xkZjVXQf3icz/dDs/Y0QQgghZvTbwDOsMZZOS14qL2kv4581gbW7VMsmZ+rhe2c7GTVpm55YJ6/z6+AVxOBFFl65Hp4nKZ1F08vAy8kSajJ26ZA1GWEtu6xzqayDN8iiGWYFzySL5tjKoWVcKyGEEEKaod8GnlVbk3iqTAyeQSHrkNiv4LVretZglQTrsbIxb9FL/zIBrtfPcEVxmGQlBHbfGcN2ob8kqi3PoWlA6NhKozQrOVvadQshhBBCekGvDTxTrOOfSooUEz/cXTRd+lS47qZel5ZJKPzcXCn7a9BvPN/xTfrV+f1ouYfmcAWvZNvwA9+xzATkFGzwG5gQQgghtdBrA880vknELiPgqLjpfcTdtqfvqj5xZNUlCUbHqp6+pnUZuGimo950vF22v+nqlatxNnTR9I3BS1xQLWLwNOpZKNPB7dNo3LDiJsLgIVFuspzRa9m35p+5i2amnwDdONqEEEJIt+i1gWczOTFpmZ0kmRqGoeh6khVf7FwB7Y7OIAlJTgxe6ShVMXjeLpqjf32xk2PWOvR1GCIubdpJruVQ7q2lBp5jP0IIIYRMhp4beBGVKzUexbVJWJp10azYni50XpotsiDJSkBdTMb0uhxrcdGs7wvS+iyaAxfNnPOYu5LmP1Zlu5xC54QQQgiZPnpt4FW50A0wdaXLuHxmXanqxjZpStuevvukv7ftZZJ1crQOXtwvLSOn3Wj/6usjOafWez1wfxxmW/SxBPJWKKuoWknLKw4fAkX7k6wMyiTkumiOrsz6Z9EsLyFS5qLZtnsIIYQQ0gf6beDFfyvr4MHMjSxvMh48Q2DZ+J3Pouk+cbedjEYTZ/skKOnzXXXuTUo3GD+EqOjvawgMDWu7GDyTJCt1uGi2nuQhUW6ZhGF2Xt/yF5HAqhi8/IcMVYYhIYQQQiZDvw28wCFZuTF4Vhr5YTtW2ybCvvpaFem2HGy8sLjdQ4Hi7W5kJ+NNxGqNtTWNwQtfJaHz7tHDFbxA59V1W8vuIYQQQkgf6LWBl1Dtimc3WyxLcU78cE7X73A2rFw0S75JhXIsVLLVP++a9cui6dqvOotmPbT72yc5r4o+CZFF06Vh141oQgghpK302sAzjW8yd9Ecth981mgWzTpyHE4PUUyXe8IbexfNCpk5SVY0Zwmv6AGBmYtmfI16umgm47mSzcppOnapi2ZOaYkQtO26zmPwkCjv+KVigkOsVkrFDW7gJpyTZKULx5oQQgjpGv028Azjm0zrlQ0n48PZ2TS7aLbNv0rhnv/e1l3WJCnJSB28pF+6Dl7FAwSTlWFfF81EB58ENZG8HAO2AtO4v/Aumn41/6YByfzNbhtJsuJr4FUlWUH+Q4Z0LCAhhBBCpgcaeKZtHdpEKzTTnGSlXXjH4NkU6baUnZeDxKwOXkUMnuM+Z4uIB4vVsi50Phlabt9Vk0qyEkScYxAek6wQQggh00etBp6IHCoiV4rIKhF5V0GbvxeRy0XkMhH5Zp36FFOdRdNKWpl7FfGijlgw17FGyyTYx5rZXB+ucaAjn1lJqJZn1m8yMXht/+4NSq0UZNHMfOI5lmE7r1EIIYQQ0hSz6hIsIjMBHA/gWQBWA7hQRFao6uWpNrsDeDeAJ6vqnSKyU1365FHkepTF1BVpGKsyIWxX8Fr28F0xrDdnjbWLZrWb30gMXtIvLSOnXUalap08r6mBK59nrJZLSYOq1R0Xt0+jcVt2XecxdNHMT7KiKSfNEDF4Jgt4WcOyKnaPEEIIIZOhzhW8/QGsUtVrVPVBAKcAOCLT5jUAjlfVOwFAVW+tUZ8xQqd8z7YxSaIREvtC5y2bnfmqW+Pu5iVZsXkoULjdUemh8TSMwQtB6NIiUbPQSVbUO7PktJOugxdGXkkMXqmLJiGEEEKmjToNvN0A3JB6vzr+LM0eAPYQkV+JyPkicmiN+hRikkXTTmDiXtXtSeYkcM6iWUOfkSQrpdkiiz4318r2Usp30Wy+ToLLcQlB2796ZW7eY8lOfMcylDCeZMVzYEIIIYTUwqSTrMwCsDuAgwEcBeALIrJdtpGILBeRlSKycs2aNcGVqJxom6bYz3GnazIJQeeTrHi4o4llRlOTrJPprUkMnlUWTVRfH75uvyMFsb3sO/uSBlUZHgdun3TRHCM53vlZNCWVPCdMmYRyF8386zjSowMHe0JUxaiLyOtE5PcicrGI/FJE9pqEnoQQQtpHnQbejQAWpd4vjD9LsxrAClXdoKp/BnAVIoNvBFU9QVWXqeqyBQsWBFPQdPKcnlCVykvaJ0/f0exk03aots3NfOLIbDOaRm5+FTLzyiTkuWiWJFmpdtFM2rrteNqVL0iSFRu3Zi13lXSprWc0LvyNnokzWMHLicGTdPmLMO6opW6YyX0yZwWvZbeQqSEVo34YgL0AHJVjwH1TVfdR1ccB+BiA/2pWS0IIIW2lTgPvQgC7i8hSEZkD4EgAKzJtvo9o9Q4isiMil81ratRphOCxPxlxtsW1/ce3jMGrSY+6aDIEz/a85RkrLqU1xvVwjMHLWGShznUNIXjBV4FCGT3TTHoFL4g8121tu4lMD5Ux6qp6d+rtVmjfLZsQQsiEqC2LpqpuFJGjAZwJYCaAL6vqZSLyAQArVXVFvO3ZInI5gE0A3q6qt9elUxEmhc6t5A3cq7o9yZwEzjF4Dt2q+wwblJZJKHXSDKVLteQQWTSt+5W6aNYZhFef6CawUd/fRdNUgJS8I5bkxag/KdtIRN4I4K0A5gB4ejOqEUIIaTu1GXgAoKqnAzg989n7U68V0Y/XW+vUo4gi16Mspu59eWUXGo3Bs23fssfvJqULirB1l41cGiti8NKbkxW8HB/Nsjp4xi6a5c2K+w+D8LyS/uRlCTUZu2zEOl002055kpVhPKmv6y0Q9y/Nopl/HdvGtRJ7VPV4AMeLyD8AeB+AV+S1E5HlAJYDwOLFi5tTkBBCyFQy6SQrEyW0u9nYHKlxF83mxpoE/i6adjF4NuTVdDOTYJZkxZas8RTqQYNtkhWjdjVYeF1fXRrGVoY6r47bun7TqQ+TGPU0pwD426KNdcWpE0IIaSe9NvASrFZqjOSN/iWTx2n1ymBlN6GsALu/gyasL8LcAtk+Lpoe2UsnQdtLlJS5eY9ls/TcV9Pu41k0iQeVMeoikk449lwAf2pQP0IIIS2mVhfNaWfw9LlyIu+WRTP9WV2kn6DbFzoPrU29mJQuKMLNRbNCZkqXvBi8oQtwviQzF83yUgtVjGTRDBCDZ3cMyzORSk0+mi27rHOproMXLnlO5Xej6DpmFk1nDGPUjxaRZwLYAOBOFLhnEkIIIVn6beDFf03S4ZsYT4NYleTpuwC62UNBA0YmZpazrSbjA0Og8Fg1sHSXVZSvyiGjy8BWGQnBKzfOjB4cGMaJjske2E5DQ8AnqYmLLaZA6Qkb5vmsI4tmu0mXWsnbNhpb6TuWlJ6D4jp4oIXngUGM+lsaV4oQQkgn6LWLpu2Kjm0bQfnEKQQe9l3rVvCmuahgnrFiViahIgbPTo1xfULHatlayUYy3XQpG7blHpqVaOavt7wSQaXbaOERQgghU0evDbyEKrc/68liiXsV8cM9i6Z9R5vYzLJrqCyLZihdTGT7XY/urrGF21gloZC0F0DRtqL39mMZtsvJokkIIYSQ6aPnBp5pfJNdOvC0vLoXnUZi8KxdNNuFj4tm5NZml0WzunzGsEG+i+Z4u9H+BjF4ji6a2f6+59rNRVPLDd9Bu7C0bmU6h2EMXn6ynORa9ikdMiqveHvRdWwb10oIIYSQZui1gVe7i2YDSQhGXTS7n2TFq3+NYw3KJFiOVx2C5+ekOYjBC+bLZ9HU2EUzcAxehWHZBYK7aLpua9k9hBBCCOkDvTbwEipXaqwTXBSnOCd++GTRtB/LfHtpQpZCF01zray9hHNdND2SrNTQr04jrCvfvNwkKwZt/EfJaTXmouk9MCGEEEJqwMjAE5EtReRfReQL8fvdReR59apWP8MsmhWxViOty+SNu3zWXQg4zyXQuG/LnDQV7tkRTUoSjI5lZ6DPiC28UZdZHYxdOE6DLpo+8/HEGLMqdF6R4bGmKgmdWFUaGL8FaTRDlb+IxZXep4q2NZFEihBCCCH2mK7gfQXAegAHxu9vBPChWjRqENPJs6lxkJXXjIumXdbGkb4tm5v5TGZtJ6MmsU152/NGKCyTYFB+w/QhROWYnrFa2aycJkQGuUEMXugsmr7W7BQgg785MXhIu2j6u6NW9R5cgzkreG27hxBCCCF9wNTAe4SqfgxRwVWo6v1o/RTKNumGW5v6k6ykX3d7tuUdgxc+w/+AsiQrPgO5ntPBilugJCtDfcK3rWMVqOvu0cMkK6HkOW4LMzwhhBBCAmJq4D0oIvMQ/56LyCMQreh1gqqpoGuK+q4nepgMjjF4NZyK9PmdUVomoSCLpo1OtnGg/iJG+7qunE6oTkLbv3ppL4CibYP3gcaqbNdxo5kQQgjpCrMM2x0L4AwAi0TkJABPBvCPdSnVFIOnz0YumgYxeEnMVUpgk0+4rWPwWrbi51PA2tZdNnIHNYnNzL5Ou8yWl+EQg/IbwxT1bqSzaPolWTGNRE2PXSXTsKElbbuu85DM37xtQPLEzXesKvfl/FhSScUC9hUR2a9su6r+tildCCGEkAQjA09VfyIiFwE4ANF04i2qelutmjVAaJe9bBuxHcQBH/Ftm5z5Ttyn0UWzLmMkq08oN0grt2ZjF82wdCAEr5LBsZ2wiyadNPGfJdsUwNObUoQQQghJMDLwROQsVX0GgB/nfNZ6zLJoWsgrca8ifjSZ6L/SdXekTEKZi6bd5/lt3dyERz6zklAtz6yf/XEJQdu/e4NSKwblLnxdJ41dNAO7hnYBVT1k0joQQgghWUoNPBGZC2BLADuKyEMw/E3fFsBuNetWO1rgepTF1BUpr02jWTStyyS0D79YMLsMIZXXRWqKOywjkBKR0260f7VG/i6aQzkhjB67a6a8rMXA7bOOLJotpzqLZux6C7/sqEC1+3LRdcwsmoCI/F3ZdlX9blO6EEIIIQlVK3ivBfDPAHYFcBGG8467AXymPrUawmLybOaalhiMMpDbaBZN2+l3yyZnqu6rFbbnwsTNLz2xzkv5X1WGw+TBgelDiDHZA+MpnW3RY/0zxwW1CvM6eGEvxKryDG2gKsnKSB0837EgFXXw8nVpogxMC/ibkm0KgAYeIYSQxik18FT1kwA+KSJvUtVPN6RTY4RMGAGMT35F6i8EPLJiZL2C167pma++tklWbMgzVkz0NV3Bs2agj9k45ljE4Jm2q2EFr+0umlUEL39Ruq3M+GvXPSQ0qtr6ZGOEEEK6h2mSlU+LyN4A9gIwN/X51+pSrEkqsyU6pqjv+BxzItSSrr+wj3lsZlkMnkn/al38ZfsVOndcOS1bwXPUxXfcdiCp/2e3jLtKBhjKulnbV0lDIiLzAfwbgIPij34O4AOqunZyWhFCCOkrRnXwROTfAHw6/ncIgI8BOLxGvRqhyoUuQWDmizSIVUnJq99FM71iZNs3rC5147MyU+WGNjZWRfxYLHTAjJIsmoU6S/UKSFUcXxVDVz6D/SnBzUWz3FVyuOoZlpZd1rkMrpmCbDnD8hf+7qhV4al00TTiywDuAfD38b+7AXxlohoRQgjpLaZ18F4IYF8Av1PVfxSRnQF8oz61mqGqTlmC6URmmBBjmAGvdgNvZHzbFBjjMqaZKC7OfSXJ1kXTJslK0jjvFPglWXGMwRt00NRndjJG5WWlVaMVY2bjBEPhG284DZR5AaTjSauOsdFYUl6PMXUFFepB8AhVfUHq/b+LyMWTUoYQQki/MVrBA/CAqm4GsFFEtgVwK4BF9anVFdo9ySTmTNuZnjZ9mqb9LpqkZTwgIk9J3ojIkwE8MEF9CCGE9BjTFbyVIrIdgC8gyqZ5L4Df1KVUU5i7aJqtMuRlPKy9TMJIFk23zm1JlKAeKQOts2garALl1cEbLVtRvvomBsuKvmdmNNuiRxZNh5IGlafLwe3TcOTQAidGURbNzUGzaJbfAwqv44qVv57xegAnxrF4AuAOAK+YrEqEEEL6immSlTfELz8nImcA2FZVL61PrWYocj3KYu+imerXoI+max28tkzSTEoXFBGVJLDLAGli+A/lx/3yYvBK+ldl2hw+hLDb86xLpW+9NJeSBprumCfTXZ3ycQMYPdNCfh28jO+353Kl6f1tPMlKex4O1Y2qXgxg39jDBQDuA3AkgNb/ThJCCGkfpklWzkpeq+q1qnpp+rO2YjXhN0myki2TYKmPC6MTbssYPB39O/V46mkXg2c3WN65tnkoUNLCSo+E7IpbqHNsm2QltEwjeei+i6ZCgxpXZaJac3+YACKyrYi8W0Q+IyLPQpRo5eUAViFKtkIIIYQ0TukKnojMBbAlgB1F5CEYzmO3BbBbzbo1RuVKjePqSdcnmZPA9lx4jVW1PaVLWZmEIp1tdsV2r4tc+1xxXjkt21bjuexKCn+T89hQlYSx88X7GwDg6wDuRBSy8BoA70V0SJ8fr+oRQgghjVPlovlaAP8MYFdEsXcJ9wD4TE06NcYwBX05Jq50GGkzlFh/mYT810Z9k1TrAfWpEz8XTVjvaP0umlJ5zkzjRIv7D8+xVwyeY0kDkzp4vgXss3TJbTA3i2bsoel7bQzlSXkx84Jsw8yiCQB4uKruAwAi8kUAfwGwWFXXTVYtQgghfabKRfPXAP4KwNtU9eEA/h3AHxAVcf1mzbrVj2l8k2G5g+yES1A+cQqBu4Nm2n2vHbM0Vfc4suhc2IxVbRCldUmuoZEkKzlJd7L9K2PwkraWxlm2SIJPgpq0RFu35vIVvGG7kHTJRTO/DF50LbteG+PyzFw0x+vg1X9/awEbkhequgnAahp3hBBCJk2Vgfd5AOtV9dMichCAjwA4EcBaACfUrVzdWCeMsKSROnjpQueOSVbagn9GSZvz7RaDZ7uiarqCZ03GeJrEuTY9hqF161KSlSJUpyMGryXPhupkXxG5O/53D4DHJq9F5O5JK0cIIaSfVLlozlTVO+LXLwZwgqqeBuC0LhVxrXbRdFs9IeHxctEM3Cd9XcwoeVQSItYsxIqUjwjnldPS2ERHZXpE3vEbX0nzHcS0WSYGz3PYLqCqMyetAyGEEJKlagVvpogkRuAzAJyd2mZaQ29qMa6DZxi/lU1pb5p+3IdRF81uZ9FUdTeWbEPworEMhA5eJi6aKRlV3Q3qiBXFP1Ux1Ge4hOdjaOatUFZR6aLp4PZpNq42moxnEoy6aPpR9QBrME6OYdmWewchhBDSJ6qMtJMB/FxEbgPwAIBfAICIPBKRm2arsZqsWiVZiYiSaNQcgxckyUo7Zmm+etqdbzuG8WR2J6QuF83BZDywi6aVW7Nh0+AumoHlTSPpJCth5JUkWSnb1oujTQghhLSLUgNPVT8c17vbBcBPdPhLPwPAm+pWrikqk2lYy3PsSCpp1EXTIcmKzbg2Ktnqn5t90U5EZnzHldOy/a/x+9GVBbz8MglhyxX4JC4ihBBCyPRR6WapqufnfHZVPeo0S5HrURZTV6Q8l8/6XTTdk6wMVnda8hBePeok2Gb8M8k6md48I6eMQJW60XVSkUVzsNltx9NZNEPUwbNd9TYxAoJffy25nk3IO36Ju/EwS2uALJol2wvvk3TRJIQQQqaSqhi8TmPqPmmaYj87ETKYv/uT9gi0jcFL/rZkkmZSuqAI25pdJrZkemI9jCdLyaiIezN5cFBVaqFKt3QWzRBJVmxj8MoGlawfaSCiYbuxulRY6Fw12Pe28jrMxBYP+qFTtjQhhBDSGfpt4Nm0NYqnysTgGSTR8GVkxcg2Bm9QBLsf0zSrvbQ8JMO5b2pF1Shus2K7ZwzeIM4y0Cm2+s6Ytgt8+fmuVraB8IuejnF2/bh1EEIIIa2i1wZeglEWTRt5HVk9mEacJ+5OZRLMYzNnFDQ1ySJppItxyxIZXlk03VdOXbaRCJNYSu8YPMMzEXpcQgghhNRDrw28QcycQZIVGxfNRFzkSdVcFk3Xvu1x0VT3JCuwq1lh5qI5OgKQ56JZ3t/YTdhyNp2NmYvi4dzJzRJahcH+x82C4uuOOk0UJVlJZ9H0faBU5aJZVE4mcl1vyc2DEEII6RG9NvBsppY2SVYSmqmD555kRTN/px1fPa2TrFiQZ6xYPRQo1MNKjUJ9JuOiadY6vItm91eXNP4vnDzHbW25eRBCCCE9oucGXkR1QWu31ZOOzzEngnNKdycXzYrtqTM8o6hMQslVYFUmwaJtYY8JXJDl+1+fQl0pdF6URXPkfUNlErK6dOQQE0IIIZ2jVgNPRA4VkStFZJWIvKuk3QtEREVkWZ36ZClyPcpiPo8ZT1te9xPuEZdA2yyaLVvCm74smsPXgzIJGRfNqjoJdWXRHOowTLLiY1C5JLw0cVGNRIbOotmSC9qAoiyaoy6a3qO4uWg24KFACCGEEHtqM/BEZCaA4wEcBmAvAEeJyF457bYB8BYAF9SlSxGDsgaGBa2r3PayEy7b2msueGXRTDIstmSapnDPjmg7Ga0qcZAnP+o3mkWzKsmIqYumfZmERIfxz1wYlF2wcXOtGDMbJxgKgxKGLSeOwUveBVnBq86iObZyCKk9xpgQQggh9tS5grc/gFWqeo2qPgjgFABH5LT7IICPAlhXoy652JcVqNieeW+7auTCqEFh23f077Tjq6fNZNTW6B3UwRsV4q2T6y5nqzaEmojb1cGbUAwe0HULLyp0HvDAmazgFelBCCGEkOmiTgNvNwA3pN6vjj8bICL7AVikqj+uUY9KbGKtQsgj7jiv4IWJYivUpUgvExdFM11COOI137c0Bo/fk0ryVpHzsll6jeHYkOePEEIImU4mlmRFRGYA+C8Axxi0XS4iK0Vk5Zo1a4LpUOR6ND5+0r5CXjZteUX68RCMxny5DdaWp/BmZysfFxfNqqHSE+uBC+NITKR/khXN+v0aknWprHKXrJaHgRxTqse0d/s0HbgrtSjz9iIpvxLSRdMki2beMW3L6j8hhBDSJ+o08G4EsCj1fmH8WcI2APYGcK6IXAvgAAAr8hKtqOoJqrpMVZctWLAgmIL2LppVMXij25ueZNq7aOrI32nH30WzvrGGOUhSLrMGQmzdfq31CeyGW8cxDO+i6R6r2SZCHrfSa7VkWzvuHIQQQki/qNPAuxDA7iKyVETmADgSwIpko6quVdUdVXWJqi4BcD6Aw1V1ZY065RIui+aovD5MMpvG3UXTpY9Z8h2gpExCqYumRRIX45bF43pl0azhYUWd34+ufPWKsmgGHcPwaI1n0ezKUZ4MVVmmReStInK5iFwqImeJyMMmoSchhJD2UZuBp6obARwN4EwAVwD4lqpeJiIfEJHD6xrXhuHTZ8Msmoby0tLqXh0bS8tv0zfzd/opz0pZisNk1MbwnzFcwhtQlc0xyqJpmJnVcTI9WMHzXNVyKWkQuWj6uai60JIFaSPyXTTjsgae18ZAnrGLZo5uXTrYDWKYZfp3AJap6mMBfAfAx5rVkhBCSFuZVadwVT0dwOmZz95f0PbgOnUpGBOAwUQ+J74qX17SYfin7unP6ITbbrQ2ZtH0XcFTVaMJsWq1MTkiJ+chQKWBYxSj6RZ5mM3q6XPs0uNbXSsVx9D0e2WLb7zhNFGUZEXj/wB/Q7kq22/RQwbWwfNikGUaAEQkyTJ9edJAVc9JtT8fwEsb1ZAQQkhrmViSlTZSudqS2W42gffDbwVvmICjDYTQ0zguzFLuwKCyPB/GDw1s9UkMTg17jm2TrJi1C3sFRsZ5Ryy8AtKFzsPIK4mzK90WToeeUZllOsOrAfxvrRoRQgjpDLWu4LWF0FPBrk8uJ4nrsXVZ0XFy0SxpM77NIgYvyCUVYAnPtltZDKKbyF5hEoPnn0XTMAav4j2pBxF5KYBlAJ5W0mY5gOUAsHjx4oY0I4QQMq30egXPNL5puBpSJXC0vUDCp4DPH3LstVFfzb6YbiL3Sre+uYXIS8eyTLIyYzzlv6I6CK96VThp6hiDl/z1dtE0/RKkxq7afXuRZuMilEE8eQpj8FAcG+eCUQxeTpKVuu9vHaYqyzQAQESeCeC9iBKQrS8SVlemaUIIIe2k3waeaR08wylUdsLVjItmOi2/Zd/M32lH4T6ZzbosVo9VbUymHwzkxagZJVkxdNG0NVjGjSePBDVpeRZ9omNYHoNoK9No3Irj3iry06FGdfAGMcT+SVbKTsJYfc+hGm15NjSNlGaZBgAReTyAzyMy7m6dgI6EEEJaSq8NvLrpyioCqWbaznXv3YSn7YQQksIwy/RxALYG8G0RuVhEVhSII4QQQkbodQye6eqIqStZnstn/Vk00+P3IYumYwxeIsNqLPsRsvKrVrA2V15TvpkS4yQrE8iiWb2CmSSmackFOAHyXTQx6qLpG4MXu3wWocOGo/0a8FDoMlVZplX1mY0rRQghpBP0egWvyPUoy9A4MMuiOZxwSbNZNO17x/9vxyzNxxCwjfeyNavyXECrYgYF1XnmiybX5voM5XgVOhd7Y6zSqKzBRdPfIJ4u8pOsRNfN8P7lP4ZJFs1xXcoNQ0IIIYRMhn4beLbtLVPaRxOi5tbwrGPw2raCF0SGYQye42C2SW9MC53bko1vC7VKVkuZhKDp/qO/XffQjFbwwh04V0lcfSWEEEKmj14beAmmLprG8txVIRU4Z9F06GjTpaipSRZJM/n+V1UIF82QPfk9qSbvvI8vpHkmWXFs13UjmhBCCGkrvTbwTJ8+m6bYz8bENJFlbtRFs9sxeDAoXVApwnhf7bJOSnbJDNUxgyYxTONuv4b6ZAqv+2QgTY9vF4NX4aI62Bh+JaorSWaK6uCp6mBnw7hoFm8vKifTjSNMCCGEdI9+G3jxX/MkK1XudKOFEqQ6xMobLXxj0jeJwWsHCp8VPMuxLJOS5MVpaoWRWJXcIhYyIt9Yn4G9mU6y4rMOZ/aQY7xf9bawLppuBvG00kiSlYp6dkXlZJhkhRBCCJlOem3g2RtEdtsFUnuMik+SleEKXjtmaSH0tEuy4iffZKzKhwaOemSNp2AxeJZJVozaOepSJqsj9l0hqqENY8dtrXk8RAghhPSHfht4Mb6FgsflBRVHUrgeWheXPZs+Ree83EXRQpdJX1Q1rJxOepfaQGEWzfR7T3PWOAYvWyah82Y0IYQQ0k56beAVuR5lGaaIrxQ4Iq8ZF02PLJqBdambEC6a5lk0y+PHxuRj/BoxuboqLynHrJBjWTQdZOTJsyFyUS1LsmL4vbIZs2NZNIuSrGj8HxBgXyvuU0VxjXTRJIQQQqaTfht41j6NVZtHGzDJSlhC6Fm7i6blWKZJVuzJ1K0LdI5tC52btQuZZCUxejpi4RUQ2kWz7Pood9EkhBBCyLTRawMvoTLJirW8bk8uJ4nrsXXpZZdF0+7zaJuFC6iFLsUy6ih2UNGPLpp+FGVZqWhiN4SZhDEXTZ4/QgghZCrptYFnmlLd1L1v4B426NdwkhXnLJrteA5flZWyjJwqBuVjOWadHD3fVVk063fRHGri66KZZNG0SLJSMWb4IgntWY02Jd++E6SqJATIoll+Xou2NZFEihBCCCH2zJq0ApPEdPJsms49T16jMXj2naM/LZmjqcI92UfWZbFyLLuDkmdAVpVaEAMf3uFWux3Py6LpVQcvI8+EaMxqC6+O668rq0t5DxmyZVtCJFkxq4M33rEltw5CcvnET6+atAq18y/P2mPSKhBCJkC/DTzbmDXL7WKyROPJ6Aqe2/60ZZIWQk/jFTxX+bZlEirluWmSXXELdY5t5Jgf64AxeINV9I5YeAWoatjkNI13JNNMHwwf0i36cM3SWCc29NrASzDNomkur9uTy0niemxdVnRcsmiafh5tq0eXOmQ49y1dwOP3pAqDEDzvIDzTczuWRdNvWEIIIRb0wYgFgAOuv33SKjTCgYfUK7/XBt7gCXiVi2bGJapY3mjacoE0+oDbOilorG9r4mgqXB6NRFgs4bnVwUu7zJaXWohiNCvViNoaazLaPp0p1Sf5j1NJAy3Xe6BO0JWoQKUDpoT8OnijhyxEkpWye0D2vjbUo9n7GyHEnr4YBYSQUZhkJWD7PBfNJpOs2O5Q+1w0Q9RJCNpsvJ+1i6ahhWfJ8KGE4TiG2CZZCdnOSFYm0VFXCV0moUxUeXxeW+4ehBBCSH/otYGXUJlFsyE9SDXuWTTte1q5aBaVSSjrY6WL/1XolWTFNblNSUd+r6rJXcEbKzjumWTF2EWz/D0hhBBCpoN+G3gFrkdjiJl7WjbbXAM5VjJZNC2TrLRsCa8qK2UZQ29A8yyaVgZY4sKYloEKA0cMrqnE5dBClzx9fDKQjuhjnUWzGDH8XtkQqnTAtJD38Ckpa+B6bYyPUbFKNxg3a1i25tZBCCGE9IpeG3im8U3mxsFo2nKTCbwvfnXwRv9OO5GN4pdkxfQYKeyMhDz51WMZxOAVpag31meYRdOr0LlD16pjaFpf0mrMQKUDpoXCGDx1vzaMBklR5PZqUOWDEEIIIROg3wZe4CC8rLwoCUHNMXgl41f2bVmSlRB6GseFOQ5lu6JaXSbBTY8x+YFOsc05MDamuYJnjSL8g5mic+tSBJ0QQgghk6PXBl5CVQyL6+oJCY+vi6bdWBZZNAtllMm30CXAipRXmQTXlVPHbaSMcVfJcNJK2uVk0SSEEELI9NFrA2/ozlVOXnxVrrykfToGr3YXzfSKkWVfx36TwieMbBjvZRiDB8sYvDwX0IqYQTHQx3VFKluCwHZ/iuRZxeBBK2MQo3bhaMlitDF5x2/oohnGHbXq3A5dQcfr4HXteBNCCCFdoN8Gnm17SxdNNJCEYNRF0y0Iry2TtBB61u+iaT9WuR5uUgYGLRI33ADKwG6fzF00w/todn91KbxzZJE81xIKhBBCCJkMvTbwEqrmgtarJ3Q+qw/HibtTN6s++Y3LroWmbRA/F806+vF7UkXeERpLduLtotmg3zMhhBBCaqfXBt4wO5zZNLTqmfkgbfnARbP+JbyRrI22fZPVndY4afobGlZZNC1Gy2atTF6Xu2gaZNHMyDfWZ6DDUI5fFk37kgZVpRlqcdEMVDpgWjDJohlqjMKV1ILrmA+yCCGEkOmk3wZe8sJwBc84pX2qX/3GU9qgsOzZIhfNQbyRexBeJMf0fFTEz42JzxdRUQfO4KGB4UOIPNmJDpGccmOzUl4ix/J6Ll3BTF6E9NAMVTpgSsitgwcZPJoB/N1RM+GaYxRdx5Gh2YKbByGEENIz+m3gWU5OTJOsJDSTZKVMA7O+bZijBdPR2L5zG9C2LmF1G8cYPIyuuAU7fFYreOYJbUIxWPEMJnE6UdXg39uqJCu528KqQAghhJAA9NrAS6iMwXNcPSHhqSNdf4g+RasoJlkkzeRbKFMkw6eva5hWqYsmvyhVFLlojryvYYzxNnkriYQQQgiZRmjgwWCiUhWjguz2qIM0nkXTtm97nr+7xqIl2MZ7qaOL5mih86o6cFK9Kuzocjjc31QWTQ+DapiV05zq/Y/bBXXRDOO2OC0UJVlR+H8nBvIq3JeLSmwksYCEEEIImS56beAV1XfKYjsRTSdZqTtGxSvJysBFc/pnaaY1C4vIuixWjpfqYyQ/J05TVSstHOM6eLaryAMdxj/zwuJaqTKSXWrrVY6Zkd128lfwZCTJSqhkJ2UumkVJVtr0kIgQQgjpC7028AghhBBCCCGkS/TawDNNqW7q7pWfRbNeRtPyW/bN/J1mwrloGib+sMw6mbdCaOKiWJeL5sCtGAEykCYiLa/nyLWvLIumvdtn5ZiZ72D7yd8TVR0ry+I8QkX/otVsumgSQggh00m/DTzDybOpi+ZwwiWDfrVn0cwZ37hvC7NousZWudXBs5CfF+OnVUlWqi0m17pugwl5Ktuiryuf7fVs7qIZMotmx+ok5FwgyWUTypitcl8uqmfYxAMsQgghhNjTbwPPun1VkpXMBw1MMm3T8md6x/+f/mlaKB1tkqw4yU+vqBqMZrqCZ0va4Ay6QmYhzfhYu6lSKqwr5l0h6n5tFIssTrJS2Gf6bx2EEEJI76jVwBORQ0XkShFZJSLvytn+VhG5XEQuFZGzRORhdepTqGfFdNA6g6GHLqQenGztAAa6SRZJIzkhyiQEyrbY5Jgkv9D5yPuaXTTztSj+lBBCCCGTpTYDT0RmAjgewGEA9gJwlIjslWn2OwDLVPWxAL4D4GN16ZOHsYumYba/rLyhW2B9j7mzafmt+rbSRdOt/9ANzXxnfV00TVwUTbNo2k6m0y6pvhlI0zJtXTTLRmUWTTcGLpqByrpXui8XZdGUeCMhhBBCpoo6V/D2B7BKVa9R1QcBnALgiHQDVT1HVe+P354PYGGN+oxh6m5mmgwim9K+jgls4aCwNyTblGQlwTmOzOJcuCQlkRwLL6qSUJZkxMRF0y2RxqC2mWowo8c+5qo8Uc3w2ASMwQtcOmAaiQxt9X7oMZCXFz+aojDJCtrxcIgQQgjpG3UaeLsBuCH1fnX8WRGvBvC/NeozhnXWSeNC55nP7Yaxwkd2om8bJmlN6ugz1uiKqkEMXgNJeMIWEg/fNuwKXgsu5gDUsZeF96+SE9SPo00IIYS0i1mTVgAAROSlAJYBeFrB9uUAlgPA4sWLaxjfb3tR+y6vIkwKdxdNlz4Whc6LPi910fSXb9PeNQPpUKZDDF7Zthq/Hl130Rx57yvPQEKxiyYhhBBCpo06V/BuBLAo9X5h/NkIIvJMAO8FcLiqrs8TpKonqOoyVV22YMGC4IqaTlxNXTQHcmtIAz82pkcWTffoveZxLReQMHRZNBkr6WMjH2Pyq8aKXDTNMrPaGmdpfYKtaol9Fk2jMgl+Wo2O2bksmnllEiSOrRy+96HSRVMLHhhAar23EUIIIcSNOg28CwHsLiJLRWQOgCMBrEg3EJHHA/g8IuPu1hp1ycU6Zs3SwgsfYZQ3pJ1L4EjfFiZZ8ZZj5DbpPpgWvC4eq0qemy7puNGg59fKRdOsMZOs2KPQcIZ7IrNAXNkoLbh1EEIIIb2jNgNPVTcCOBrAmQCuAPAtVb1MRD4gIofHzY4DsDWAb4vIxSKyokBcrVS7aLqtnpDwNOuiadO2IJF8eZaRWnQplOGbjMOpT1mSGX5RqslPblLeogYtci4e3ucIIYSQ6aTWGDxVPR3A6ZnP3p96/cw6x6/C1J1ruL0qpX2S8bC5LJohXDTb8BQ+m6HUFptzEdJFs7wOnliX3jBXKOmfLrTg78pnc62Yu2iGzKKZuPJ22PoQZFw0PcVV+GgWu2i2Y/WfEEII6Ru1FjqfdoYT+fIZknUdvIzcOjP7jbgEOmYFbUMcjWu5gASbeC+X2KahO+6ok6aJgWOki20dvFTzYIaAZcxVtYE7bBeKgaxO23cS18GL3wdamS26TynyT2QUCzj99w5CCCGkb/TbwHNOSmK3vd4VPP9UKW2YooXS0WRC6mOQ266omhc6d9cn7AqZTVvDGDxHXXqNhn8wUxiDVzIMz507InKoiFwpIqtE5F052w8Skd+KyEYReeEkdCSEENJOpqJMwqSpdtFkDF7bqd1lr0B8oBA8+1IdnuOF0CHqUxKDV2eZhPpET5zxMglhsmiWtvEagWQRkZkAjgfwLEQ1Yi8UkRWqenmq2fUAXgngbc1rSEi7OOD6EyatQu2cv3j5pFUgLaLXBt4wZq68nb2L5mgMXp2Mumh2P4uma0p4NxdNC/mprJUDOSiffJvEtLmu0qTdg8O5aNrH4FVLDLsSFap0wPSQUyYBSRbN1AcemGT7LUyy0oJ7x5SyP4BVqnoNAIjIKQCOADAw8FT12njb5kkoSAiZLvpgxJJw9NvAs3bRNEuykjCY9NdbJyHvpWHX9sTghZpI1r6raZdZE3fQyjIJbqTj24JWSQhr4QXHt15iW0gnWQknsyAGr2SgFtw5ppXdANyQer8awJNchYnIcgDLAWDx4sV+mqE/E8k+rIj05VwSQkbptYGXUJlkxVqeuy6knCYPrY3rW9E5L3fR9Jdvg78rn31/3yQzpLo8gb/rbbWEvCadzlTaMlT1BAAnAMCyZctodxtC44cQ0lX6nWTFsJ17Fs1knDqzaLov4Q1cNMOpUxum7rRFDCex5qtqdi6a49IV5ZPvyEXTbAnPtRajwj8D6UAm7K7lShfVpF0NWTS7bDwmWTQx+E6EicErTBKFgphOacnq/3RyI4BFqfcL488IIYQQb3pt4EHL09gPMXO1zMbE1DGBHRtzxL5zi8Frg4VnWrOwCJtz4eLmN4h5y2bRrEgyUn1Nubkcpt2Dg51eA33TaMX3q44yIqFKB0wziWHl+50YyIv/lmXRzI3BQytuHdPKhQB2F5GlIjIHwJEAVkxYJ0IIIR2h1waea8xacYOiOlL14VfoPL0OMN0EK5Ng0sZjMNuyFdVJVtx1ieRrrQ8Yysc2bBd0Ba8NV7M/ivDf29I6eEV9+nG4g6OqGwEcDeBMAFcA+JaqXiYiHxCRwwFARJ4oIqsBvAjA50XksslpTAghpE30PgbP5Om3fYr65rJo9g3fLJp19SlqGiq60/oaLChM7YNLb5NC53XQ5fiw7J5532dMYvDcupESVPV0AKdnPnt/6vWFiFw3CSGEECv6vYJX4HqUxdS9L+seNnSTqzMGL/+1Ud9BmYTpfwzvG0dmk9F0cB4dDLCR81HpomhxTTkaLJpa6gmRjMPmWom+X2Xyhu1C0T0XzZyDIzKSRdM7eU7JUEDxeRSRWuOLCSGEEOJGvw28wJOT7ETVpvaa+5jptPyWfTN/p5lwLpomSVbcR7PvWt7BVZX0hLyOGDeztmatJ1bGoeWEfjBTlmSlWIegKhBCCCEkAL028ABTF03LDIZuqhADnJOs1O6imd84lItikDIJvlk0nY5hWRZNflOqyU9uMvK+kfNarQchhBBCpoNeG3hVLmQJxi6aAzdCyXzuoJwhoy6abkt4bXgKP9DRNQYvK6dsLI8BRl00KwycRrJohs22aJdFs7pMRNQu6BpeLLu75kfWM8D/vJa7LxfeJ6Udq/+EEEJI3+i3gQezVQTTenbZCZc04KMZJovm9E/TXA2dBJt4r2EtNYcYvJEsmlqxgieVR14dZ/Hp/R3axr710uxirhQw0rsOF83umnfDe9bmQDtbfX/Lv46FFh4hhBAylfTbwLM1iNyqJNRsQNml5R/p2aIVvFCH0OhchBrLyJg0e2jgrIPBGFby7ILwwrZrVtTUE/p7W7aCV9inV0ecEEIIaQe9NvAAGD39dk1R3+VVhEnhvghl39EqRq7oc4Mskmby/eNAvWO1nPqUu6jWRYc9NMdWo4Nl0TQYs+ozQgghhEyeXht4VS50CYMYlUp5o+3rSAM/NmY26Mumb+bvNONbLsDKRXMQx2UjfzyOKYpBK0syYnDsHctDpPUJFqtl6ZGnqCoTIYN2oQhl9EwP48dmEE/qcJ3mUeVJXnQd28ZkEkIIIaQZem3gWeckqXKny2wfTsTqQwteG/Ud+GiG0qY+mpxI+oyVNlbMSjJUyXMjbQQErTNnmWQltMxKWYGMnjawObiLZr7Asuu4BbcOQgghpHf028CD4UTQPoUhqQn3Quf19ilsa7CCFVqXkgF9BQTtwa+JCTkrZ5mEPqGyaJa2oYsmIYQQ0hp6beAZZ9FMtTdhEIM3cJOr7zm3XxbN5O/0P4f3z6JZngp+dKzRPmbyMSa/qkxA1KZqVdhel6w+vsduVKZdFk2TGERm0bRj6No6+t5dYPSntExCbjep9d5GCCGEEDf6beBpeYxQgqlxkJ1cNlAlYdQl0DYGr0VZNIeGjlv/bNxS+VgOMXg5cZomBk51mQQ342zECPA8dgOZsHXR1IoYRHOj23zMWHaXLbwYl+s0j8qHEMg3Im1jMgkhhBDSDD038Kx7VGwtiGFpKMmKdQyeY79JEErH2gqd58k3EVRTDN5QQNj1WSsDz7hdwCQrrbiaw7B5c1h5TmUS+nO4CSGEkNbQawMPMEwRbiszftrdg0WExvHNomk3lr98kxUsH/k2eKfTdzmGgcpEOIxcp/CJkvUM8He9dfxOeY5LCCGEkHrotYFX5HqUxTTF/lj8Tw1p4MfGzBnfvLPGf6b/MbyOHVw7bNxlB4fDYeKbzaLp76IZt3VJcBLLD+W2KBCra7kqBlFS7ULRBxfNoWtrWBfNQg+Eouu4yweZEEIIaTH9NvAMkmAADnXwkiQr2Q01kDbOOu2i6Zk8IzspLh3LISlJYZKVMgPPQJ/BVscVyNBJVqyNsaoDgHquvy6bHskh3exzceTIKzy3Bdfx0EBvwx2EEEII6Q/9NvAsp5amK3jj4zRDl5OshMJoVwMdjwAheN6T5+B18Ezb1VEwr1lRU09ow6rEvqvQI6gahBBCCPGk1wYeAKOH37aeSMlqET2YwuOcEj5w/Ni4+PzGpTF4NvKd9JfMe3sZI/0D9/GNCSwdt8NfvmzZFu/zanIPLKnHRwghhJDpotcGnrmLZtLeLIvm0EUzfBr4sTE9ZLcp46C/i+aonNKxBn1c6uCNlq0oj8ET4/PnalxFLprJe19XPnN9TWLhaqmDF8gddZrJulR6u95WuC8XXcemruuEEEIIaZZeG3i2uLoq1ZtkJW1QWPYdiReb7mlauGNoEIPnYzRblq2o2i/f0xIlWWm+DIHNiHUkWekDm6fGRbNHB50QQghpAb038IxcuTwzN5JwOHtoOnS0c9G0xy6JSwD9J3A9lpeJqHHcDn/3sitnvu6orrfALh9jQgghpM302sCrcqFLsHW1zGbRbK7QuXvSmGl/CO+b/t7ORdMli+a4u1qURbM8i2Rl4h4Pl0NB5FIZzJVPYLw0Z5LCf3DMgq4wJrKDiZw6su7AoXa1rNB53nWcjQUkhBBCyHTQbwMPhjF4g1ihKne6ZMI1mmSlzgmQj5GmBa+nEd84Mpc6eC4reFkXzaokI6Zuv04Gi4xesyFWeoyzaA5VKJaXaRuC7HewywSrbyjlZ6LoOjatD0oIIYSQZum3gWc7MXEtk1DjDMjHSMsmBJlmQulnk2TFST4sLe7KFTxPNOwE3PQ8TKhKwtQ/qAhJ8Bg8x/tXm5I1EUIIIX2g1wYeYLaqYfuAPJtFk4TD3UXTIYbNIYum6edV23x0GfYpf28vL2wcY61ulB3+6g1dW+P3vtlRHRt1uRQFIYQQ0mZ6beAp1NBF0ywd+JhbWgMuTKOrcB5yAuhSJ7765ZUxKBzLIWgtO+kGTFw0zbNoutXBi7NoZsp3uGLnopmMWZZkJXya/VDxhtNMXXXwyrJolg0x5Yv/hBBCSO/ot4HnUVYghLwQjLpodj/JirecmsfKJlnxHcvX/U1VJ1KGwM5FM6iJF1DWdNOUi2aPDikhhBDSCXpt4AGGKcItn5AnqxZdXkWYFK5uYXWXMSiUUbaCVfMFknXd817p8evegMCU6A67D1bkRLGXZ3AicrNodvcQE0IIIa2m1waeaZLx4XzKLNlA0j7PbS84Xlk0Nff1dOKZEt7BXdbWSIjKCIwe00lm0ZS4DINvBtKhvGp983Rw2eZKv1w0w7nepuVlKbqObcvHEEIIIaQZ+m3gqd0KnqmL5lgdvBqNJx/ZbXTR9E2yYnK8XI2EbJm4qusrMsDMDrxrkhVFqnRAgBU82yyaJoXO68ii2eXVpeTBw+bASVbK6+Dl6RFvn/qHQ4QQQki/6LWBZ+vjZJpkZezzWpOspF+H3Z9pIpiuRpULmjOaq1fwfGPwAicxMW5nPmpIA2HaH1SEpLkyCW79CCGEEDIZajXwRORQEblSRFaJyLtytm8hIqfG2y8QkSV16pOrY8BWg9ZJDF6HVxEmhW+h8zr75Lt0Vq9g1aFL1CfwBRj4GNYZJ9flEiVZj4JQLpqlbXKOZ3ePMCGEENJuajPwRGQmgOMBHAZgLwBHichemWavBnCnqj4SwCcAfLQuffKwd9G0e5RdlX48BLZZG0f6Biqx0AT+LpqxnBrHypY90CoZIuZuv3aqjOgz3J8ArnyWWTSrykSk24YglDvqNBO6TEIisTQGr9RFkxBCCCHTxKwaZe8PYJWqXgMAInIKgCMAXJ5qcwSAY+PX3wHwGRERDZs3fYQb73oAl9909+C1Db+9/i5s2FSs2tW33Zf7+S9X3YZVt95rNZYpyb4AwKo19+Knl99i3Pe+9ZsGr8/64y3YYtbMoLqF5MY77w8i56Lr7sT9D24qbXPL3euc5f/5tvsG5+CedRsAzKvsU3bOVq3xu26uv/1+/Obq27xkpLn57nVG19i6DeXHOM1Vt9hdt+Wy7gkipw2kv/shOP+aO3DTXePX/q33rC/td/Yfb8W82eP3DgHwzL12DqUeIYQQQgyp08DbDcANqferATypqI2qbhSRtQB2ADAyIxWR5QCWA8DixYu9lPrVn27DO067dPB+j523ruwzf95sAMCnzvqTcVsA2G7eHADAv37/D7ZqWjFDooQLP7j4Jvzg4puc+h79zd/VpF1Y0sfXhu22jM7Ff/30qqBjJccvGmM2Tv/9zTj99zcPth/w8B2KdYrlv+ZrK0vH2HqLWZjhsEwzf8vZOOuPt+KsP94avXc8dgnbbTkHF113Z6W+IzqUjDl75gzMmz0Tp/12NU777Wov3WzGbQMH7bEAP7n8FjxiwVZj2+bH1/LJ/3c9AGDbub7nNer/wR9dXthm/6Xbj/eL729vPjn/3jFzhuDq//ccL90IIYQQYo/UtVgmIi8EcKiq/lP8/mUAnqSqR6fa/CFuszp+f3XcpnDJYdmyZbpypfkEM8sd9z2Im1Ird7ttNw8P2WpOZb/rbr8P96zbWNlup223wE7bzAUQuYutuvVerN+42VlfEx6y1RzMniGVT9rzeMSCrXHt7fdh0+bpd7SaO3sGHrFga2dXwz/fdh/uW199DgFgi1kz8Midqsd6IF4NnDdnJm69Zx1uvXv0HDxyp60xN2d1AwA2bVZcefM9lckydtpmC+y07VwjvdPcfu96/GVttCLje+wAYO0DG3DDHeYrqTNnCPbceRvMmFE85s1r1+G2e+2v2zK2nDMTD19Q/eBmmlFV3Hn/Bmyfc29SVVx1y73YsGkz5s+bjUXbb+k91tVr7itddV20/ZZjRvPmzYorb7mn9N6x927zvXQDABG5SFWXeQvqCb6/kQDwmy+9LZA2hBBC8jjw1R/3llH2+1jnCt6NABal3i+MP8trs1pEZgGYD+D2GnXC9lvNyZ00VfGwHcafpFchIth9522s+7niYgQAwKN32TawJtPJ0h3tz2EV8+YMjbedtpk7MO5NmDlDsNeu9R37HbbeAjtsvUUwefPnzcb8ABP2NA+dPxcPne923XYZESm8T4kI9nxouPuKiOCRO9kbxDNmSG/uHYQQQkibqDOL5oUAdheRpSIyB8CRAFZk2qwA8Ir49QsBnF1n/B0hhBBCCCGEdJnaVvDimLqjAZwJYCaAL6vqZSLyAQArVXUFgC8B+LqIrAJwByIjkBBCCCGEEEKIA3W6aEJVTwdweuaz96derwPwojp1IIQQQgghhJC+UGuhc0IIIYQQQgghzUEDjxBCCCGEEEI6Ag08QgghhBBCCOkINPAIIYSQhhGRQ0XkShFZJSLvytm+hYicGm+/QESWTEBNQgghLYQGHiGEENIgIjITwPEADgOwF4CjRGSvTLNXA7hTVR8J4BMAPtqsloQQQtoKDTxCCCGkWfYHsEpVr1HVBwGcAuCITJsjAJwYv/4OgGeIiDSoIyGEkJZCA48QQghplt0A3JB6vzr+LLeNqm4EsBbADo1oRwghpNXUWgevDi666KLbROQ6TzE7ArgthD4dg8dlHB6TcXhM8uFxGSfEMXlYCEW6jIgsB7A8fnuviFw5SX0c6cP3pw/7CHA/u0Qf9hGYxH7+03+GkFL4+9g6A09VF/jKEJGVqroshD5dgsdlHB6TcXhM8uFxGYfHpJAbASxKvV8Yf5bXZrWIzAIwH8DtecJU9QQAJ9SgZ2P04Vrpwz4C3M8u0Yd9BLq5n3TRJIQQQprlQgC7i8hSEZkD4EgAKzJtVgB4Rfz6hQDOVlVtUEdCCCEtpXUreIQQQkibUdWNInI0gDMBzATwZVW9TEQ+AGClqq4A8CUAXxeRVQDuQGQEEkIIIZX01cBrtStLjfC4jMNjMg6PST48LuPwmBSgqqcDOD3z2ftTr9cBeFHTek2QPlwrfdhHgPvZJfqwj0AH91Po8UEIIYQQQggh3YAxeIQQQgghhBDSEXpn4InIoSJypYisEpF3TVqfphCRL4vIrSLyh9Rn24vIT0XkT/Hfh8Sfi4h8Kj5Gl4rIfpPTvD5EZJGInCMil4vIZSLylvjzvh+XuSLyfyJySXxc/j3+fKmIXBDv/6lxcgiIyBbx+1Xx9iUT3YEaEZGZIvI7EflR/L7Xx0RErhWR34vIxSKyMv6s198fYkfeb1PXKPqt6RpFvx1dJPtb0EXy7u9dRES2E5HviMgfReQKETlw0jqFoFcGnojMBHA8gMMA7AXgKBHZa7JaNcZXARya+exdAM5S1d0BnBW/B6Ljs3v8bzmA/2lIx6bZCOAYVd0LwAEA3hhfD30/LusBPF1V9wXwOACHisgBAD4K4BOq+kgAdwJ4ddz+1QDujD//RNyuq7wFwBWp9zwmwCGq+rhUium+f3+IHV/F+G9T1yj6rekaRb8dXST7W9BVsvf3LvJJAGeo6qMA7IuOnNdeGXgA9gewSlWvUdUHAZwC4IgJ69QIqnoeokxsaY4AcGL8+kQAf5v6/GsacT6A7URkl0YUbRBV/Yuq/jZ+fQ+iL/Vu4HFRVb03fjs7/qcAng7gO/Hn2eOSHK/vAHiGiEgz2jaHiCwE8FwAX4zfC3p+TAro9feH2FHw29QpSn5rOkXJb0enyP4WkPYiIvMBHIQoazFU9UFVvWuiSgWibwbebgBuSL1fjQ7eZC3YWVX/Er++GcDO8eveHafYhe7xAC4Aj0vifnIxgFsB/BTA1QDuUtWNcZP0vg+OS7x9LYAdGlW4Gf4bwDsAbI7f7wAeEwXwExG5SESWx5/1/vtDSBGZ35rOkf3tUNUu7ud/Y/S3oKvk3d+7xlIAawB8JXa5/aKIbDVppULQNwOPFBAX0O3ckzYTRGRrAKcB+GdVvTu9ra/HRVU3qerjACxEtPL9qMlqNFlE5HkAblXViyaty5TxFFXdD5H75RtF5KD0xr5+fwjJo+y3pitkfztEZO8JqxSUnv0WlN7fO8IsAPsB+B9VfTyA+zAMK2g1fTPwbgSwKPV+YfxZX7klcZGK/94af96b4yQisxH94J6kqt+NP+79cUmIXRXOAXAgIpe6pHZmet8HxyXePh/A7c1qWjtPBnC4iFyLyLX76Yj89vt8TKCqN8Z/bwXwPUQPA/j9ISRDwW9NZ0n9dnQtvnLst0BEvjFZleqh4P7eNVYDWJ1aaf4OIoOv9fTNwLsQwO5x5rs5AI4EsGLCOk2SFQBeEb9+BYAfpD5/eZz17gAAa1MuV50hjon6EoArVPW/Upv6flwWiMh28et5AJ6FKGbkHAAvjJtlj0tyvF4I4GztWIFNVX23qi5U1SWI7htnq+pL0ONjIiJbicg2yWsAzwbwB/T8+0NIlpLfmk5R8Nvxx4kqFZiC34KXTlit4JTc3zuFqt4M4AYR2TP+6BkALp+gSsGYVd2kO6jqRhE5GsCZAGYC+LKqXjZhtRpBRE4GcDCAHUVkNYB/A/AfAL4lIq8GcB2Av4+bnw7gOQBWAbgfwD82rnAzPBnAywD8Po4ZAID3gMdlFwAnxllnZwD4lqr+SEQuB3CKiHwIwO8QByXHf78uIqsQJUs4chJKT4h3or/HZGcA34tzx8wC8E1VPUNELkS/vz/EgrzfJlX9Unmv1pH7W6Oqp09OpVrI/e2YsE7Ejdz7+2RVqo03ATgpXvi5Bh35bZKOPVQmhBBCCCGEkN7SNxdNQgghhBBCCOksNPAIIYQQQgghpCPQwCOEEEIIIYSQjkADjxBCCCGEEEI6Ag08QgghhBBCCOkINPAIIYQQQsjUIiKbRORiEblMRC4RkWNEZEa8bZmIfGrSOhIyTdDAIwSAiOwQ/3hcLCI3i8iN8et7ReSzNY35zyLy8sxne4rIiSIyQ0R+4yDzuPgH8DhHnY4VkbfFr78qIi8saXu0iLzKZRxCCCHEggdU9XGq+hhEBdQPQ1TPF6q6UlXfPFHtCJkyelXonJAiVPV2AI8DIiMHwL2q+vG6xhORWQBeBWC/zKanAjgPwD4A/uAgejmA7VV1U9nYqrrRQXaWLwP4VfyXEEIIqR1VvVVElgO4MP69fhqAt6nq80RkKwCfBrA3gNkAjlXVH0xOW0ImA1fwCClBRA4WkR/Fr4+NV9d+ISLXicjficjHROT3InKGiMyO2z1BRH4uIheJyJkiskuO6KcD+G1iaInIU0XkYgAfA/A2AD8G8NcisjJHJ4lX6v4Qj/3i+PMVALYGcFHyWarPsSLydRH5FYCvi8gSETlbRC4VkbNEZHHFcfgPEbk8bv9xAFDV+wFcKyL7WxxSQgghxAtVvQbATAA7ZTa9F8DZqro/gEMAHBcbfYT0Chp4hNjxCETG2eEAvgHgHFXdB8ADAJ4bG3mfBvBCVX0CotWtD+fIeTKAi5I3qvoLVX0cgCsB7AXgpwAOU9VlOX3/DtFq474AnonoB2wXVT0cQzeWU3P67QXgmap6VKzjiar6WAAnASiMXxCRHQA8H8Bj4vYfSm1eiWjVkRBCCJk0zwbwrviB6bkA5gIofYBJSBehiyYhdvyvqm4Qkd8jenp4Rvz57wEsAbAnIteQn4oI4jZ/yZGzC4Ar0h+IyJYA1quqisjuiIy9PJ4C4OTYDfMWEfk5gCcCWFGh+wpVfSB+fSAiQxEAvo5o5bCItQDWAfhSvJr5o9S2WwE8qmJcQgghJBgi8nAAmxD9Bj06vQnAC1S16PeTkF7AFTxC7FgPAKq6GcAGVdX4882IHpgIgMviVbTHqeo+qvrsHDkPIHqyCGDgXnkxgL1E5FIAjwWwMutq6cl9Lp1iN9L9AXwHwPMwNGqBaB8eyOtHCCGEhEZEFgD4HIDPpH6DE84E8CaJn7CKyOOb1o+QaYAGHiFhuRLAAhE5EABEZLaIPCan3RUAHpm8id0rvwDg9QDeDOBzJa6WvwDwYhGZGf/QHQTg/yz1/DWAI+PXL4ll5iIiWwOYr6qnA/gXRK6hCXvALRkMIYQQYsq8pEwCgJ8B+AmAf89p90FEyVUujdt+sEEdCZka6KJJSEBU9cG4tMCnRGQ+ou/YfwO4LNP0fxG5RqY5CMDXEGXC/HnJMN9D5GJ5CQAF8A5VvdlS1TcB+IqIvB3AGgD/WNJ2GwA/EJG5iFYo35ra9mQAx1qOTQghhBijqjNLtp2LKN4OcRjCa5vRipDpRcZXtwkhTSAi30NknP1p0rq4ELu+vFVVXzZpXQghhBBCSAQNPEImhIjsCWBnVT1v0rq4ICLPAvAnVb120roQQgghhJAIGniEEEIIIYQQ0hGYZIUQQgghhBBCOgINPEIIIYQQQgjpCDTwCCGEEEIIIaQj0MAjhBBCCCGEkI5AA48QQgghhBBCOsL/B1fR71oBqpe2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# make our generative model with two components, a fair die and a loaded die\n",
    "gen_model = hmm.CategoricalHMM(n_components=2)\n",
    "\n",
    "gen_model.startprob_ = np.array([1.0, 0.0])\n",
    "\n",
    "gen_model.transmat_ = np.array([[0.95, 0.05],\n",
    "                                [0.1, 0.9]])\n",
    "\n",
    "# the first state is a fair die with equal probabilities and the\n",
    "# second is loaded by being biased toward rolling a six\n",
    "gen_model.emissionprob_ = \\\n",
    "    np.array([[1 / 6, 1 / 6, 1 / 6, 1 / 6, 1 / 6, 1 / 6],\n",
    "              [1 / 10, 1 / 10, 1 / 10, 1 / 10, 1 / 10, 1 / 2]])\n",
    "\n",
    "# simulate the loaded dice rolls\n",
    "rolls, gen_states = gen_model.sample(30000)\n",
    "\n",
    "# plot states over time, let's just look at the first 500 rolls\n",
    "fig, ax = plt.subplots(1,2,figsize=(15,5))\n",
    "ax[0].plot(gen_states[:500])\n",
    "ax[0].set_title('States over time')\n",
    "ax[0].set_xlabel('Time (# of rolls)')\n",
    "ax[0].set_ylabel('State')\n",
    "\n",
    "# plot rolls for the fair and loaded states\n",
    "ax[1].hist(rolls[gen_states == 0], label='fair', alpha=0.5,\n",
    "        bins=np.arange(7) - 0.5, density=True)\n",
    "ax[1].hist(rolls[gen_states == 1], label='loaded', alpha=0.5,\n",
    "        bins=np.arange(7) - 0.5, density=True)\n",
    "ax[1].set_title('Roll probabilities by state')\n",
    "ax[1].set_xlabel('Die')\n",
    "ax[1].set_xticks(ticks=[0,1,2,3,4,5],labels=['1','2','3','4','5','6'])\n",
    "ax[1].set_ylabel('Roll')\n",
    "ax[1].legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bac1856b",
   "metadata": {},
   "source": [
    "Now we'll train and decode, seeing if we can recover our hidden states, transmission matrix and emission probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "fc78a160",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model #0\tScore: -26192.686320731424\n",
      "Model #1\tScore: -26325.519153251287\n",
      "Model #2\tScore: -26247.37450144936\n",
      "Model #3\tScore: -26186.69051567507\n",
      "Model #4\tScore: -26198.92880039691\n",
      "Model #5\tScore: -26180.251318807237\n",
      "Model #6\tScore: -26219.376180428026\n",
      "Model #7\tScore: -26182.606841632427\n",
      "Model #8\tScore: -26228.56673432129\n",
      "Model #9\tScore: -26219.247899254307\n",
      "Generated score: -26175.4414621326\n",
      "Best score:      -26180.251318807237\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAEWCAYAAABYGk2QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7cElEQVR4nO2debwV1ZXvf+vcC1zmUVHCHEEUAWVQEI2CL2oSQ9JpacW0QseERLuj/dKJiSbdSVqStB1j1LSJoZ9GE4kDakykk0gU1BYUkSgOgAKCMs/TBS7cYb8/quqcOnVq2HvXeE6t7+dzP/ecU7v2VLtWrb322qtICAGGYRgmeQppV4BhGCavsABmGIZJCRbADMMwKcECmGEYJiVYADMMw6QEC2CGYZiUYAHM5AIiep6Ivph2Paod7sdoYQEcI0R0HhEtJaIDRLSXiJYQ0QTz2Cwiekkhr8FEJIioPr4a5xMi+h4RPZR2PeKGx1D24AsRE0TUDcACANcBeAxAewDnAziWZr1qASKqF0K0pF2PrEFEdUKI1rTrwcjDGnB8DAcAIcTDQohWIcRRIcRCIcSbRHQagHsBTCKiRiLaDwBE9Ckiep2IDhLRJiL6ni2/F83/+81zJpnnfIGIVhPRPiJ6hogGmb8TEf2UiHaa+b1FRGe4VZSIehHRr4hoq5nPU7ZjXyKidaYG/wci6mc7JojoeiJaS0SHiOhWIvqoqfUfJKLHiKi9mfZCItpMRLcQ0W4i2khEn7fl5dl2m+Z2LRF9CGCRX9vNYx8nojXm7OO/AJBH2y8FcAuAK8x+XWn+3s9s716z/V/yutBE1JuInjbrvpyI5thnN0Q0goj+Yub1LhH9ne3YA0R0DxH9j9mHy4joowrn/oKI/khEhwFMiXIMqfQjo4kQgv9i+APQDcAeAA8C+ASAno7jswC85PjtQgCjYDwYRwPYAeCz5rHBAASAelv6zwBYB+A0GLOZ7wBYah67BMAKAD1g3DSnATjZo67/A+BRAD0BtANwgfn7VAC7AYwF0AHAzwC8aDtPAPi92daRMLT75wAMBdAdwCoAM21tawFwh5nXBQAOAzhVoe2/BtAZQMeAtvcBcAjA5WZ7/q9Z9hc92v89AA85fnsRwM8BNAA4E8AuAFM9zn/E/OsE4HQAm6xra9Z3E4B/MOt5ltmnp5vHHzDHydnm8XkAHlE49wCAyWa/NUj2o+wYUupH/tOQE2lXoJb/zEH9AIDN5sD9A4C+5rFZcAhgl/PvBPBT87PbzfMnANfavhcAHAEwCIbwfA/ARAAFnzJOBtAGxwPCPHYfgP+0fe8CoBnAYPO7ADDZdnwFgG/avv8EwJ3m5wvNPuhsO/4YgH9VaPtQybZfA+AV2zEyr4GUAAYwAEArgK62334E4AGXc+vMPjnV9tsclATwFQD+13HOLwF81/z8AID/Zzv2SQBrFM79dYxjSKkf+U/9j00QMSKEWC2EmCWE6A/gDAD9YNwQrhDROUS0mIh2EdEBAF+BoYV4MQjAXUS03zRj7IVxk3xECLEIwH8BuAfATiKaS4Zd2skAAHuFEPtcjvUD8IGtPY0wtLWP2NLssH0+6vK9i+37PiHEYdv3D8wyZNu+SabtZp7FtMKQHvZzg+gHo08OOer6EZe0J8DQHO35O+t5jlVPs66fB3CSLc122+cjKPWZzLll7YpyDCF8PzIBsABOCCHEGhgai2WHdQtD91sYWvIAIUR3GHZi8km/CcCXhRA9bH8dhRBLzTLvFkKMgzEtHg7gGx559CKiHi7HtsK4QQEARNQZQG8AW/za6kNPMw+LgWYZgH/bLex94Nf2bTAeLFa9yf7dBWffboXRJ10ddXVr9y4Ymn1/22/2sjYBeMFRzy5CiOt86qNyrrPuUY4h1X5kFGEBHBPm4sm/EFF/8/sAADMAvGIm2QGgv7VIZdIVhubVRERnA7jKdmwXDFPBUNtv9wK4mYhGmmV0J6Lp5ucJpjbUDoattck8vwwhxDYY09CfE1FPImpHRB8zDz8M4B+I6Ewi6gDghwCWCSE26vYLgO8TUXsiOh/AZQDmS7TdDc+2w7BpjySiz5HhcnUDyrVGJzsADCaiAgAIITYBWArgR0TUQESjAVwLoMJVTRheB08C+B4RdSKiETCm7hYLAAwnoqvNvm1nXpvTAtqne25kYwjq/cgowgI4Pg4BOAfAMnOF+hUAbwP4F/P4IgDvANhORLvN364H8O9EdAjAv8GwkQIAhBBHAPwAwBJzujhRCPE7ALcBeISIDpr5f8I8pRuA/wawD8b0eQ+AH3vU9WoYdsw1AHYC+GezzGcB/CuAJ2BoQx8FcKVmfwDGVHsfDA1zHoCvmDMD37a74dd2IcRuANMB/AeMdg8DsMQnO+shsIeI/mp+ngHDZroVwO9g2F2f9Tj/n2AsOm4H8BsYD65jZl0OAbgYRr9tNdPcBmMh0hfNcyMbQxr9yChChlmHYeKFiC6EsdDVPyBp1UNEtwE4SQgxM+26MNmGNWCGCYlpbhpNBmfDMFf8Lu16MdmHd8IxTHi6wjA79INhT/4JDP9ohvGFTRAMwzApwSYIhmGYlFAyQfTp00cMHjw4pqowDMPUJitWrNgthDjB+buSAB48eDBee+216GrFMAyTA4joA7ff2QTBMAyTEiyAGYZhUoIFMMMwTEqwAGYYhkkJFsAMwzApwQKYYRgmJVgAMwzDpETysSA2vgS8/0LixbrSfzww/JK0a5Esh3cDr/0KaD2udl6/s4ARn/Q+fmg7sOJBoC2GlxX3PgUYc4X38aYDwKv/DbTE/MLpTr2As78MFBT1li0rgHf/HJyu7+nAyL/Rq5sbrS3AsnuN/glDfXtg/LVG+1VpbjLqcPxwcNq0adcAnD0b6NA1OG1EJC+An/0+sPlVpP9yVQH0GJQ/Abz6D8DiOeYX2WsggK79/AXwW48Dz/9QMV/JskHAqOnegm/dc8CiW2Mo21kPAB+9CDhhuNqpL94OvPtH+NdNAB26RSuAd74DLPy2+UW3X8x2d/sIcGZQjHwXNi0Dnv1uyDokgdnOPsOB0z6dWKnJC+C2ZmDYxcDn5wenjZOnrs+OJp4kraaG+o33gc695c55+kZgzR/907Q1G/9v2Qa076RfPycv/Cew+Adwf5uOVbbZpq/+Fej9Ue90YXjnd8D8WXoafmuzMYOY/bx3moXfAZbfp1s7j3LNul41Hxh+sV4e+zcBd56hP7OxxsUXFgIDz9HLIwl2rAJ+MSmeGZwPyduARRtAGTA9Exl1yRtWm0lBG6FCcF8V84342lr19Ctfp03K9SgE18MLmTEv08c65Vp56xKm3QBgRVvMwj3vR9h2apJjARzDgK8GdG7KVAWwxI0RV9mq9fAidQEc4sEUWgAncG2ioNjOZMPzpiCARTYuBgtg+XOkBHBMmg4LYH0yoQEnMDuJgnxpwBm4GCyA5c+hQrBmwBqwNzUhgDU1w6rRgCVMXTHAJoi8EbsJIuKHa+YEsIYgkpn1ZV4A58UEwQI4GVgAy58jLYApBwI4Rg0YiNYGGYkADqkZsgD2JecCOIfvw9MSwBIeI3Fd18wI4BCCSMbsFocAYA1YHhbACZNbDVhjsUxWA45VAPs8LJNwdYpdA47BBskCWB4WwAmTWwEcowmipjXgpEwQUQrgCB5M7AccKzkWwLwRQ5pUBXCONmLo5u9Xrj1vHVgDjhX2A84bOotlsn7ArAG7k7oAzsJGjAy4nvqRn40Y7AecKjqaKhUAiAA7LJsgPEldALMGHEh+NOCsmCBYAEsjtRDGAtgTWT9g3fw9y+WNGNLwRoyEYQEsj6wQjGNmkzkBrLMRgzXgTNzzfrAGnDAy0+paREsAyyyE5cQG7BcW0xOh4Aec1Y0YuhpwtXlB5MIGnIGLkVKHp06sGnAtC+CwGzGqVQMmACE8hlgD9oUFcN7MEDqaaiYEMG/E0CoXCN8vYcx1LIB9ybEATsfonjqsAcdXDy+qeSOGdT4L4FjItx8wkFMBrLhYlqoA5o0Y2kTVL6EEcLXZgGteAGdFA86zAI5DA055ES7uMVXVAjgLGnAGfP/9yNciXAYuBgtgearBBMEC2L1ce966sAkiNlgDZgEcjPRGjBT9gDMtgKt4I4Z1Pm/EiAUWwCyAg2ENOIKNGFUaD9g6v+Y14JDudpqwAGYBHIzsQlguBHDcJoiMbcQAECpyYLUIYCDcg0YT9oLIHSH8gP12gaXtBxy7AE5qI0aUi0DshqZEPgQwa8CpUqsbMapeA+aNGKlT8wJYCGhpYHGQ640YNegHHLdnTVV6QWRBA64SP2AgLwIY2bgYudWAeREutnp4kbobWpobMVgD9iNhAWw1jv2AUyM2DTjtaGisAbuWa89blzxsxADCudtpko4AzsLFyLUAjksDZj9gV3LvB6z4Cqy0yI0GnIXpCAtgeaQ3YrAJwpW8a8BZuN9lSOFFvSyAWQAHwzbgnG/ECOkHnIX7XQbWgBOEBbA8vBEj5xsx8qIBswBOjlgc36uAqvUDTmEBsKwecW/EqGE/4Czc7zKwAE6QXGvA1eQHnKL92bUebANWrkMW7ncZWAAnSK43YlSTBswbMbTJykaMLNzvMtS+AOaNGKkTmwDmgOyepK4Bp7wRIwv3uwzsB5wgLIDlYT9g/fEiq3TEJoAj8MENLYAzcL/LUPsacJZMECyApcmEDbhaBbDkmI9LAEfRL2E3YmThfpeB/YAThAWwPLwRIwIBnJIfcCQCmP2A44IFMAvgYNgPWN9tsWY0YBbAccACOHcCuBb9gKtBA05pIwYLYHlYACdIbjdi1KINOMMbMaQFcEwbMVgAy8MCOEFy7QdcixsxElhp17lBUzVBRPRgYj/g2MixAM6rCSImDTiuN52kaX+uqIuGN0DabmhRCWC/9wH6VyIb97sModqpB2/EYAEcjLQZoIb9gK26VJUGHNHMIFd+wDUtgHkjRurUpA2YBbBn2ewHLA/7ASdJnm3Aqv2fZjyGDJkgoHGDFgVXUN/EtQgXhQYc1g84AwqXDCyAE4Q1YHl4I0apLrwRI506JAEvwiUIC2B52ARRqgtvxEinDknAAjhBciuAdTZiZH0nXEIr7YnYgHkjRmqwAE6QXG/EqEY/4DTsz8666NiAeSNGJu53GVgAJ0huNeBqNUGk/EYMqy5V5QXBGzGUqH0BnCU/YPaCkCZNM0DmNmJUkwBmP2Al8uMHnAUBzBqwNGmaATK3CFdtAjhtDZhNEH7wRgwWwMFkwgTBAlgZ3oihBvsBJwgLYHnYD7hUF92NGFUtgNkPOC5YALMADoY14FJdtP2Aq3kjBpsg4oJNELkTwAII3BbroCr8gDkcpWvZLIDlyY8AzsAFybUfMGvAenWJ0w+YN2KkDgvgBMmtBlytGzFqPSB7HBsx2A9YCRbACcJ+wPKk6gecJQ24Gk0QWfADzsD9LkPtC+AsbcTIswZcTSYI3oihTWZMEBlY85GBN2IkCAtgeVLdiEFAUBxeFsDeZbMfsDz58QPOwBORBbA8aZsBggRfpgUw+wFXjwCueRMEa8CpU20bMazyq1YA14oGzAI4DlgA504A68QDzroGnGQ8YN6IkUodkoAFcILkVgDraMApL4RlRgNmP+DU6pAELIATJNcbMarID9gqP40FQNV6uJGqCSILfsAsgP1gAcwacDBBfRW3e2HQ1L8mbMC1+kYM3ojhR44FMG/EkCZQAMd8XTNjgohbAEfsBpWZjRgZ8HqSoeb9gC0yIYDzaoIIswjn0VdxuxcSAfC7TjXwUs5imihtwFGaINgPOA5yrAGzCUIa1oDl6uGGqgCuORNENQng3JggMjAlYQEsTyZswBkRwMoaqkLfZFYAUwgNuNpswLVsgmANOF2EgNZ0XdoEkQMBrOwHXAsCmDXguGABnDsBjBo0QSRlAw7jBywx62MBnC4pCOD6REp542Hg6D5g83LjexYuiFWHD18GXv55unVJCl1BaQmPTcvc+6qlSS9f6fILwPa3va9TS1NyGvD+D9XGy85VpXNl8t/yenTjsXEn0K1f+HyoALQ169WraX827ncZqAA0H/Fu59mzgbpoRWYyAnjJncCuNcbnDt2Bhu6JFOtLoR7o2g9Yu9D4yxM9BqilJwJ6DATWP2f8edFdMV9ZegwENr0CbHnNP03c9BwEbHgBeOZmtfPqG4AuJwan6zEQ+OAl4y8qekyLII+BxsNbtd3286uBHgMNAezVzvFfiFwAk1CwaY0fP1689prPTeBF08GS9tWuI1DfQT2POGg5bnR4nijUAR26qp8X1Fe6+crQ2gIcb/Q+TgWgoVs8ZdsRAmg6oH5efQdj3AfR2gwcP6yevx8N3aNZ9Lbfw2nVIQmaDnjb+UO0g4hWCCHGO39PRgNO4ubQob698ccEk2Zf1dUDHXukU7YdonjrUdcuG+10I6v3cNQkPDuvEuMMwzBM7ZGMBswwTGZobm7G5s2b0dTUlHZVao6Ghgb0798f7dq1k0rPAphhcsbmzZvRtWtXDB48GFQtttkqQAiBPXv2YPPmzRgyZIjUOWyCYJic0dTUhN69e7PwjRgiQu/evZVmFiyAGSaHsPCNB9V+ZQHMMExuufPOO3HkiJor6vPPP4/LLrsskvJZADMMU7MIIdDW5u2/rCOAo4QFMMMwiXPrrbfi1FNPxXnnnYcZM2bg9ttvx/r163HppZdi3LhxOP/887FmjbF7dtasWbjhhhtw7rnnYujQoXj88ceL+fz4xz/GhAkTMHr0aHz3u98FAGzcuBGnnnoqrrnmGpxxxhnYtGkTrrvuOowfPx4jR44sprv77ruxdetWTJkyBVOmTAEALFy4EJMmTcLYsWMxffp0NDYaG4D+/Oc/Y8SIERg7diyefPLJyPqBvSAYJsd8/+l3sGrrwUjzPL1fN3z30yM9jy9fvhxPPPEEVq5ciebmZowdOxbjxo3D7Nmzce+992LYsGFYtmwZrr/+eixatAgAsG3bNrz00ktYs2YNpk2bhssvvxwLFy7E2rVr8eqrr0IIgWnTpuHFF1/EwIEDsXbtWjz44IOYOHEiAOAHP/gBevXqhdbWVlx00UV48803ccMNN+COO+7A4sWL0adPH+zevRtz5szBs88+i86dO+O2227DHXfcgZtuuglf+tKXsGjRIpxyyim44oorIusrFsAMwyTKkiVL8JnPfAYNDQ1oaGjApz/9aTQ1NWHp0qWYPn16Md2xY8eKnz/72c+iUCjg9NNPx44dOwAY2urChQtx1llnAQAaGxuxdu1aDBw4EIMGDSoKXwB47LHHMHfuXLS0tGDbtm1YtWoVRo8eXVavV155BatWrcLkyZMBAMePH8ekSZOwZs0aDBkyBMOGDQMA/P3f/z3mzp0bSV+wAGaYHOOnqSZJW1sbevTogTfeeMP1eIcOpfgxVvwaIQRuvvlmfPnLXy5Lu3HjRnTu3Ln4fcOGDbj99tuxfPly9OzZE7NmzXJ1FRNC4OMf/zgefvjhst+96hQFbANmGCZRJk+ejKeffhpNTU1obGzEggUL0KlTJwwZMgTz588HYAjDlStX+uZzySWX4P777y/aabds2YKdO3dWpDt48CA6d+6M7t27Y8eOHfjTn/5UPNa1a1ccOnQIADBx4kQsWbIE69atAwAcPnwY7733HkaMGIGNGzdi/fr1AFAhoMPAGjDDMIkyYcIETJs2DaNHj0bfvn0xatQodO/eHfPmzcN1112HOXPmoLm5GVdeeSXGjBnjmc/FF1+M1atXY9KkSQCALl264KGHHkJdXV1ZujFjxuCss87CiBEjMGDAgKKJAQBmz56NSy+9FP369cPixYvxwAMPYMaMGUXzx5w5czB8+HDMnTsXn/rUp9CpUyecf/75RaEdlmTCUTIMkxlWr16N0047LdU6NDY2okuXLjhy5Ag+9rGPYe7cuRg7dmyqdYoKt/5NNxwlwzCMjdmzZ2PVqlVoamrCzJkza0b4qsICmGGYxPntb3+bdhUyAS/CMQzDpAQLYIZhmJRgAcwwDJMSLIAZhmFSggUwwzCMAl26dIksLxbADMOkSlDIyCRoaWlJpVwWwAzDJI4zZOStt95aEVYSAH79619j9OjRGDNmDK6++uriuVOnTsXo0aNx0UUX4cMPP8SBAwcwaNCgoiA/fPgwBgwYgObmZt8wl1/5yldwzjnn4KabbvJMt2HDBkyaNAmjRo3Cd77znUj7gf2AGSbP/OlbwPa3os3zpFHAJ/4jMJkVMvLgwYN4/PHHK8JK9u7dG3PmzMHSpUvRp08f7N27FwDw1a9+FTNnzsTMmTNx//3344YbbsBTTz2FM888Ey+88AKmTJmCBQsW4JJLLkG7du18w1xu3rwZS5cuRV1dHS666CLXdDfeeCOuu+46XHPNNbjnnnsi7SoWwAzDpIIVMvLrX/+6a1jJlStXYvr06ejTpw8AoFevXgCAl19+uRgU/eqrr8ZNN90EALjiiivw6KOPYsqUKXjkkUdw/fXXo7Gx0TfM5fTp01FXV+ebbsmSJXjiiSeK5X3zm9+MrA9YADNMnpHQVOPCChnpFVbyZz/7mVJ+06ZNwy233IK9e/dixYoVmDp1Kg4fPuwb5tKqQ1A4zLheYso2YIZhUsUrrOTUqVMxf/587NmzBwCKJohzzz0XjzzyCABg3rx5OP/88wEY3gkTJkzAjTfeiMsuuwx1dXXo1q2bVJhLv3STJ08uKy9KWAAzDJMqF198Ma666qriQtfll1+OQ4cOYeTIkfj2t7+NCy64AGPGjMHXvvY1AIZm/Ktf/QqjR4/Gb37zG9x1113FvK644go89NBDZa8NmjdvHu677z6MGTMGI0eOxO9//3vXenilu+uuu3DPPfdg1KhR2LJlS6Rt53CUDJMzshCOspZRCUfJGjDDMExKsABmGIZJCRbADMMwKcECmGFyiMraDyOPar+yAGaYnNHQ0IA9e/awEI4YIQT27NmDhoYG6XN4IwbD5Iz+/ftj8+bN2LVrV9pVqTkaGhrQv39/6fQsgBkmZ7Rr1w5DhgxJuxoMUhDAP/zjavzudX9n5n+88KOYNdl7gBxracXf/mIpdhw09mr/n9NOxI8+NzrSegLAgje34o6F7+HZr12AQkFuK+Ltz7yLdTsbce/V4yKvTxi2HTiKq/57GRqPqYXd69ejIx7/yiS0q/O3Vr20dje+8fhKtLR5T2s71BdQIEJLaxvuvPIsnD3E2Nu/p/EYpv/yZRxqqqxb+7oC7p5xFsYN6umZ77UPLMebWw5Itgjo3bk9nrz+XHRqX497Fq/DA0s3Sp9r59KRJ+HWz54RmO5YSysu/8XL2H6wyfV4ny4d8OR156Jj+zqteti56fGVWPxudWi21ng42txa9nu7AuEnf3cmJg7thSvnvoL3dx/Wyt+vX5es242vzy8frw3tCrhv5gQM79tVqzwdEhfAy97fg/oC4cJTT3Q9vmDlVqz4cD9mTfbOY/+RZry95SAmDO6JnYeO4ZX398ZS12/MfxNHm1vR1NKKTu3luuq/Fq+LpS5h2bj7CDbsPoyLRpyIE7vJ2aje23EIKz7Yh8amFvTs3N437dtbD2DbgSZMH9cf9S7CeuWm/Vi17WDx++ptB4sCeNO+o3h/12FcMPwE9OvRsZjm6PEWPPXGVry7/ZCvAF6yfjcG9+6MswZ6p7HYsLsRr7y/F7sOHcOg3vV4dcNetLYJXDLypMBz7fzv2l1YtmGPVNoDR5rx1pYDmDC4J045sfzmfn9XI5Zt2IvdjccwoFcnpTq48fL7e9ClQz0mDu0dOq84OXi0Gf/z1jYAwIiTuhav3fGWNjzx181Yte0gxg3qiWUb9mJ0/+4Y2a+7Uv5B/fqOY7weOHocf3xrO9btbKxtAdwmgNNO7oYffW6U6/FlG/agLWBxwDr+ubH98cr7e/DGpv1RV7PmsBZcvvSxodI354NLN2LFB/sCrwdQuia3fvYMNLSr1Dh+svDdMgFsz9P6/A+TB5c9mHceasJTb2yVGA/AhaeeiG99YkRgPZ96fQteeX8vLMWnTQgM7NXJczx6cf28FXhvR6NUWqusvzmrP646Z2DZscdXbMayDXsR1XpYWxswdnBP5fYkzbqdjUUB/LHhJ+CWTxo7xw41NeOJv242grSbnXLJyJPwj1NOUcr/CbNfvcaOdU3+/TNnoGP7Ory34xD++NZ2qbEeJYl7QbQJAb/ZfIEocHXW6rwCGenj7rRaWCwu9Zl8VCfrOvlYFYpYfeSVvTOalD1P63o762Z9DxoPImBMldfDKl+Y50L63PJ85MddW7F9lccKjvqERaUv0sReR/tlt655mxDFMaUyZov5mJLNa+xa/W1lrTLWoyQFAewf2q1AxlPcN482q/MIJJE+LDUgf32FgBckKQCB0jXxulmc5YoyDdhK4y6Ag26KNiF/kzqFuqEQaNzgRNIP5jaPB4z9t6gEsEpfpIm9jm6f24TemK3Mx71fncJdZaxHSeICOOgJLaPR2jtPRmMOS9LTkjgoPfFVNGA5AWhP4y2AnRqwTQC3ud9ostph0KzKrR52E4SeAJYfF36zg5JGrlwFV9qEKGp/WaZc6JZ+t89Q/B5cQQQJVOeYi/pBKEtKJgjvDjWmdsF5AJYJIv5pg0j3fYGRUHpoyZ+jMj0O0lacv9uvWVtRQJUnIokHgDCnqrIPFmebjBmZ1KmOfHRMEN4acFRKRNAMMyu4mR3sn4Wwjwv1/INMCk6FoZg+4Xs9FROEnwAuUPBgtA/oJGzAtaQBq9mA5bUCEaBhV9qAS3kKD+FtffcbD6p2wqJQbyvlradhyd+sxZvd5W5TmWXIUDU2YFslyUUbbmsTnmsDUvkHjN1KG3CONGC//pQRqPYno4zGHJbaEMDGf5XBTEUBKJd/kGnJjj3PkoDysgF7V0DVTuimAetM2VVMX/4acHmasOiaVJLGfr3KP9ttwJXHVfP3ekgKUw5Zwl9lrEdJCjbgYA04SKCKMg04fsN50iujceB84sugohUE3fgVJghbp3oJURnt0Et4e9ejNMW1yta3Acul9ZsdUMSaV7UvwrnagDUkcFC/OvspVxqwX3/KuPfYtblk3NCqXwLrTOeCXHnsBJuWnCYI+2d3AeV0GXMvV+3BUmqT3QasN8WVtwFb57jlY/yPzA84YIaZFcptwPbfDc8mux+w7vUBvPvVKYeiNgXJkrlFOJmHnXMRLu4+qwUNuGgrVbjiKgtEAv43vvOYsF01Lztu0E3kd653PRyajqbNlIikx51fHaPWvIJmmFmhXOutvO5tQm/huJSH8V94XCXngzcofVwkL4DbgvyAZTTg0pORiMqms3GQXxuwvFYQbFoK1oA93dB8KqBuAy5vk+6UXcX05bsRQ2GWIVtWVSzCeZggjO/h3dCCNFrnYqXKWI+SbPoBB6wu25+MKg7xutSGAFYTVPa0shsx/K9r+feoNmKo3qTONukKLEtLk8FvKh29Dbi6F+GA0sJ6mEW4IPOVs59UxnqUZM4NjUjO8R4oLcLxVuRgktiI4asBO+4it1gQztPlbMBWWrVFOLsGrGdjlB93MiaI/PkB+2vAQoiyHa+qBPWr5yJcwipwOjZgn1JlNFq7X2WhwG5oMqjaSo20xn9ZLwh/G7C3CcJrgdC+IOOFlw+xdz1K9bXO144FITnw5GJBqNfBjarxAy6b/juPGWZInTFrzwPw24hRPl5zswgnEGADLgTf8Ha3HhmNOSw1IH+1TBAq02MhhK+7UOVOuOBFOOs3v9JVb9Lo3NBUYkGUl+1Xn7BU4yKcm+lJiNKCWDg/YPeOdfYTObxjkiKbNmBFNzS2AQejFw1NXjiouqHZ8wxy04p2I0b51DTMIpzqVmS3YmTMLCpU5yJc+TFDqdIbs6U81Bbhon4QypJBG3CwScE+7UzCBlwLbmh6GzHKzw3KX2URzm0jhtcilcxGDPVYEKWydRRGFdOXnw921G5o1WMDLn122wFp94LQuj4Bi2qVNmDr9xrXgGX8gOXjAfNGDFm0NmIoLsL53fj2Y5aGUzrX30bqd31V2+U0q+hO2VVMX0mZIMLETkgafz9gK8iSfnssoe5vA3ZZhKt5DbhNJhZEQB62J2MysSDizT8JwsSCkLYBB1xXi/pC+UMz0AYsoQHrx4LQd0OTtgF7hNt0q08YwrhtJY2fG5olA8KYIIL61Rm7JGpTkCwZjQURpAGXnoxJ+O/Vhg1YfUFDxUVKZYdjXaE8kE1QuEa5jRhJL8Kpa8D+fsDKVXApx+yLKpDAfotwVjgCnTFrzwPwvkbOKHhRuwPKktFYEP552DWmJKYOtSGAjf+p+QGXacAFhwnCqlvleU5zRWW5anZCp801aGemXz7yfsBJacD6NtOkKXcBKz9mBTqyNmRxLIgIiTYecDLGc52ss2Y3VvWXtaeV8XcN9gMufa5zmCD8NLcgQafqhlaaalrnh/ADljVBWHX0aJ9Vj7CE8ZtNmvI1gcpFOBFSA5YxQbjGIc6DBhw+FoTx34oFYeUbFzp5Z81uHPTONjfiigVR7/AgCNIQ/YRTSXgH19FeD3s0tDA2RlnzjP0c1/pE8CaGMAIrTbxiQYR5oMhsxLCPmRzFgvCfIqkswlmxIKx840LnomTNbKHnB2z8l7cBex93asBusSAI7jZSqXjAihpwWSwIzYDs9vL98HezK08ThjCLVmniHQtC7eFanofx39sG7C742QYssbjhDMhu5RsXehpw1gSwKQRUwlEGuPKU5y+vAXuaIDw0YJmNGGnEgrCX70eQl4e9PmHQifeRBSoEobkbNkx7gmNBVC6+JuHS6iSDfsASsSDaSmmTMJ7rPBUzJn9TjwXhXPBwX4RT14BVbdvONoWxAdvz8cP3AWPegZHYgIv3ReisEsV52S0ZEM4NLcgE4RWDQrmoUKSyCOdvA5Z3QzP8gMt/i4PaMEGo2weVY0H4XNfKjRilPMPZgK10siYIpwasHwsCkN+m7VXHODTg6jNBuGuiOgvHpTyM/96LcJXXXWb2HTWJCmCZDtWJBQEg1lfH64Soy9winJYGrCBk2uRNEM5Zjt8CYdB4UI8FYfyPIhaEvXw//NzD4nBDqzYN2F0QhtOAVWNBWOUkrTclKoBlOtQI8+efT1GQF5JxH6klDVhlLMdngijPM0hD9F2EU/QVrfSC0IwFoeQh4v2AidKLR8fXOwu474SLORaEi8JQID1lKwwJC2AZDVgxFkQhugHshZYNOEaNXIckYkGoaMBusSDcFgiDpoW6O+EswR3kPueFiumLY0H44xULIoxJxfmgdeLmDlvzNmCZVU1VN7Qk/PdqQwM2/scaC8JnNNmLrbQBe9ctaFpYOjewimX1iCIWBCD3oA3y8rCnCUM1xYKw46oBt4XbWBLkX+2MBQHkwgZsFuongCUCstsFeVZjQWRPAKvbB6ONBeG0AdtNEP4CKlINuFCucYaJBWEv3w8/0wAvwrnZgKOIBWH8l40FARhjo6b9gGU6tFZiQWR1ES7OWBBBMxv7Z7tmEsoGrGgnrIyGpulnqmD68lt8jtKLp5piQdhxzpyKsSBC2LSdD1onbjOfHJggzEID3NA4FkT06Pi7qi7CBdn2LZxTPT/BEWwDtvJXXYQzvof3Aw5O66eZ5jUWhJ0sxIKwzqlpE4TME1rVDY1jQcihM9WOMhaEM/i1eywILxuwdwXUA7Ib/xONBWHbOFSZTwwmiMS9+8PhHQsiikU49+NuCkMSscWdJGsD9hmIFiqLcETRriJ7l6dzTrYksI6giTIWRJkbWgEOG7CVxsME4bPQpR6QvVzjDLsIpxYLwi2f8jRhqLlYEBLywotgG7C74GcbsMQ0gGNBqKPj7xrkylOZv9/iqlMDlluEk3VD04kFYbz2JsFYEC4NjCMge/X5AbtrwOH8gP1NO7mMBVGaIvlMVSERC6JonK+M7xoHtRILQnUgF/tWytUKLrHMSpTZgFF+vYK8BOJahBO2caSKFblNaSecWz4Ks4wghE85WcZZX6LyWBB618fAzwThLNg5LpMglUW4oFgQ8otwFPikC0MY4Z41DTgoVoMbxb6VKkDeBkxEZXn6LYQ5zRUVxSpOu+0ap5VrmCmuzGX2K6fkr5pfE4TXYlgUNmD/tyJX1iPp2zZ7sSAKMjZg47+hAUc3hfMurxYW4dQHsoqLVJvw34hR7obmtAF7PxzijAURbpVdfu1BbiOGeh1Uysky7jZgEeqBEmSjd/cDrnkbsFmor6akZwOOdyOGzjnZksBhbMBRv5Sz0gbsPSZkA7LrxIIIYzO1HjZhN2JEGgsixDvU0sRpkiz5Aes/UCjg+riNufzYgH1Xy2XiAVeaILJnA86WAA5yE3MjuY0Y3g8H+Z1wwXW018NYhKusmywqC5QyL+WMYrjUigbs9AOO66WczmxzsBEjuEPl4gGX0mbJC8LNtSor6Lhbqa30y7+SyDnL8Xs4BMeCULMT2s0q4ba6KjycAsJtWvUJS5gHSpp4BcURtvtcFZmNGJVboGt8I4bMAJHbiFES5IlsxJCMbOYW4SsrxL0RQzUamr172tp8FuGCNGBFX1G7ZhTOxggzH3kTBG/EcMdLEIaJbSFnA648p6YX4aKPBaG2GKKLrDAts2tmLBylTsyDaDdi2ASwI+BSsA042ASh7IbWFo2fqdJGDI9wm/Y0YaheP2Dnd+ulnKXvqgT1q7sfcI1rwLKxIADJV5FTMhsxZLN221yQFfRiQci7SAXbgMvzdfaVvw3Yv1x7XYMoswGH2GmltRHDRwOOYs2gWt3QvHakxboRo83L9FHTAji4Q2U0i7KA7DEuwlnXQt4G7P45CwS9MsgNlb4NtgFT2WdnLAivzTnSsSAkR3L0NmC52YFXOVG6oYV5h1qauC+G2fyANRrEsSBckFkwkdEs7II8yimcd3my6bKrAesswgW58lTmL6sBVy5Y+vsB+5VbSieDsW7g8AMOcYPLXGY5GzBrwBbWa8nCBJgPkiNuC785iAVhFuozQGQ0i3I/4OimcF7I24DVz0kKPRuwmpBRWYRz9pXXTRa0Mq2jxTptjHHHgvCb+ZUUCOUqKJWTZYJiQejZgCU04Io4xDWuAcv6AQP+N325G1p8JggLWeHuXFjKEkGvDHJDVcjImJasfJ195SUE5WNByN+klVtdpU8tq5dRfnDaoJdyWhp5WMIIrDTx8kYIEwvCytcvGE/+4gFL7NRRMUGULcLFKPFks7a/HyxrGzF03NDUhEzQzMb+2cUG7LMIF2UsCHv5ca6y2wkqJ6rFH52+yAIVgrDgfEDqtcevX93Ga83bgFU0YLntp9mKBZFlDTiRWBB+17UsHGV6sSDs5UcTCyLcIpz1O8eCKFGMBeGzgUUuX2+B6jZea94GLIqCM1hTCrIBW8I3S7Egsr4IpzqOo4wFYT/iFgvC60wif59qHS2WUP7SxzBbXWU3qfiVE+TrLEsYm3aauMVksJsgdB8ofusHHAvCg+JN73vTlW720ssRo6mjG/I2YPvnbAngIBOBG0pCxsWv0i0v67Ozr3SFk87Ck6Vxhpmyq5i+glz0DM1LuQoVVK8G7L8Ip/tA8dvZ5hkLIuENVCkJ4LA24FK6rMaCyJj81XJDizIWRFlAdodmIoS3H2+QcNKxE1qaTmKxIAJd9CiSNYywNtO0cPcDDhcLwjrPq1/dbcC1vggnsaop87pvu8aUzEs51dNlTQNOOxaE86WcwtFX+jbgUjpZiJKPBREc/0S5CpXlhNjZlyZOP2xLEIb16lC3Add4LAiZJ7TMTS/KNGAq/hYXtbIIp+vvKm0D9g3IXv5ZNhZELItwBacNWP5cex5G+cFpg+zvUWle1esH7PzutAHrNcjfBuwekD0XGnDoWBBtNhtwAiYI2ayzvAinEwsCkF+YUI4H7Ogrr1MtbdWvXCOdugki3CtvrPJlzDPBGjDHgrB/d9qANfMtePer23jlRTjIu6GVBHD8i3Cy9jn7tcueH7CuoJGbHgfbgEsHra2m5efqacA6mylKb1yorJssSrEgfMJt2usTFtW4GFnB3QYclR+w+zG38ZobP+AoNmJYWaj4quqi5YaWuXCUehqw7PQ40A3NNtIq/IDb9H1kdXxFicrfuBDODzg4rZwNmDVgi8qNMrr5qrqh5cQPWCZqVpDW49SAORaEP/o2YLmFifCxIJJbhDNWxxHqHWpRbtOOSvOqJTe08gekrg1YfREuFxqwX/Qp2VgQzkW4WMJRmi8Ur5lYEFqanpwJJjgWhO2zS0B2Pz9g/7FgzqoURnIUbmiqsSD8x3y0sSCqbyOG87szWJJ+vt7xgHMYC0JGA5bRLASSWYSzBoBszuUXO1sSWMcNDVCxActrwJWaiX4sCLf8ZdIKhNuIoWL6kpkdRDl8q80E4SYIhRCAKO141cE3FgS8TR9JkkEbsNwiXKUfcESVtGHd+HomiOjrEwadgOyAqg3YPx+LaOMBq2uxFX6mGneB+jZtv7yidUOrPhNE+Xe7DTjMw8R/Ec5N886NDdhnsUZGA7YN6DhjQVgXL6+xIAB/Vx5n/mo24PK+8tSAA3wz9WzAlp+p/pRdNRaEf/yTnG/E8HAH0104tgj0Ay64l5skmXVD8/cDdnNDi0MAK2rANs+HrGnAQSYCL2RNEOp+wHLnBgkn/VgQIpTXQKTbtAtR24BDZ5Uo7otw+gvH9ny9utUzFkRtmyDMQkObIOwasJk+YrcvIUTx4sneG+XxDbIlgYN2qnkhOz0OFwsi2XjApVgD+lP2yGNBRDBewti006RCEBZKfsBhNGBVNzRZc1uUZPClnOVp3fOx24CD0+tgz05nI0YWTRC6Gw4yHQtCyw84Og04U7EgQrptpYWXINQds/Z8VTZicCwIlG5Uv46wv17HsuNE3W/2/GRvDmE7K2sbMQTC+LuG14CdeVZuxNBdhLPXU7Z8Mmc4+lN2lY0YQviXEZXmFXbjQlr4uaGFEcDkM3vK5UYMlVgQQTvhnG5oUXeczoJa1jdiaMeCkHiYqNjr3DZieJ0aJJx0FtKsNoXTgOXXHvy2Wlt5RTFcqtcP2F0Q6i4c2/P1twG7LcLpl6dDhhfh/PJxW4SLpo6lMtTtueXnRFufsAQJAS9k7JOqe/adQdb9FgiDhJOOndAS6mFiDZRMX8Fpk3JDC2PTThOvxTDdhePyfCr71VrfcXd/y4EG7LtaXrDS+ms9VhaJ2IBlTRAZd0PTjwURlLfxX80E4aibT0D2oPUA1ZvUudNKa2YgEbfagmNB+OOUB8V3wsXkhua1WBkU/D8OUrIBe6eRjwecnAZcKxsx9GNB+DdGZou5M89KP2DNRTgNzd5y+wrnB2z8l/WR9rcB53sRzoklCFvb4lmE85qJ58gPWMYNLeimK08fvQ3Y/bPvOW3qQjspdLUJmemxqv+pW0B2fT9g9UW05GNBBK97RBkPWCUuRhaxy4BQfsAe/tVeM/HaD8gusVNHKiC7yyJcFO/UcpZhofNSzqz5Aeva02QWJrymdF44hWo4P2C9Vy1ZNkajjDA2YDkNOAk3tGp9J5wTayy0tIb1Aw7SgN1MH/rl6ZBBP+BgzSKJWBCibFebXObO+AZZQt8PWF4DlrcBl89awsaCUL1JnW9cCOcFEZw2yASR91gQTqx7OqwJwmtRrfTgLf+95t3Qig33GSFFzcJnZLvFgoh66uCcIsud435+FtCOBSHhIqW6+FO6ZqW6+Qdkj34RTtgW4cL5AYdfhIvOBlxet2rFqn9LWxQ74Sp/93rwpuGGVp9EIV98cDk+2HMEB442AwD8+tTqlK89thKd2te5ptmy/yiGntC5LP29L6zHo8s3RVbnVtuVeOad7Vi5aX/gOUeOtxY/37N4PX677MPI6hOWD/cewen9uiufVyDCC+/twsfveMEzTaviw8Z6AF9y54sgAB/sPYKTu3d0TWvFA/Yqf+ehY/4Dyq18AlZ8sA/rdzWaZaidb+UBAHc/tw6/efkD37Rb9x/F4D6dffNavmGvbx/LsO/I8VDnZwWrbxet2YnuHduFyIdc+7XVYyZeIML+I8c9r8OCG85Dh3p3maRLIgJ4YK/OaF9vKNsndm3ASd0aPNOO7NcN08f1x+HjLZ5phvXtgqkj+gIAOrWvw5cvGIpNe49EW2kAo/p3R58uHbDtwFHpcyaf0hvt6wvYezhbN8Owvl0wfVx/5fO+cN5gvPDersB0Z/TrjqkjTvRN88O/GYXTTu6KTu3rsWb7IbSaiwLD+nbB3451r9slI0/C+l2NnlrwsL5dMLp/j8D62bl60mD8+e1tAIALO7bDKSd2UTofALp3bIcvnjcEWyXGhn28ujHz3MF45p3tynVwY3DvzmhoF62QiIun/+k8vL5pX8XvF512It7eehCtbW04Z0hv7fyvmTTIs1/P6NcdUxzj9dNjTsauQ8fKdrTaIdUnvQSkYvMYP368eO211yKvBMMwTC1DRCuEEOOdv1e5wwrDMEz1wgKYYRgmJVgAMwzDpAQLYIZhmJRgAcwwDJMSLIAZhmFSggUwwzBMSrAAZhiGSQmljRhEtAuA/75Lb/oA2K15brXCbc4H3OZ8EKbNg4QQJzh/VBLAYSCi19x2gtQy3OZ8wG3OB3G0mU0QDMMwKcECmGEYJiWSFMBzEywrK3Cb8wG3OR9E3ubEbMAMwzBMOWyCYBiGSQkWwAzDMCkRuwAmokuJ6F0iWkdE34q7vCQhovuJaCcRvW37rRcR/YWI1pr/e5q/ExHdbfbDm0Q0Nr2a60FEA4hoMRGtIqJ3iOhG8/dabnMDEb1KRCvNNn/f/H0IES0z2/YoEbU3f+9gfl9nHh+cagNCQER1RPQ6ES0wv9d0m4loIxG9RURvENFr5m+xju1YBTAR1QG4B8AnAJwOYAYRnR5nmQnzAIBLHb99C8BzQohhAJ4zvwNGHwwz/2YD+EVCdYySFgD/IoQ4HcBEAP9oXs9abvMxAFOFEGMAnAngUiKaCOA2AD8VQpwCYB+Aa8301wLYZ/7+UzNdtXIjgNW273lo8xQhxJk2f994x7YQIrY/AJMAPGP7fjOAm+MsM+k/AIMBvG37/i6Ak83PJwN41/z8SwAz3NJV6x+A3wP4eF7aDKATgL8COAfGjqh68/fiOAfwDIBJ5ud6Mx2lXXeNtvY3Bc5UAAtgvPq01tu8EUAfx2+xju24TRAfAWB/VfFm87dapq8QYpv5eTsA622MNdUX5jTzLADLUONtNqfibwDYCeAvANYD2C+EsN4ca29Xsc3m8QMA9N8smR53ArgJQJv5vTdqv80CwEIiWkFEs83fYh3bibwVOa8IIQQR1ZyfHxF1AfAEgH8WQhwk2/u9a7HNQohWAGcSUQ8AvwMwIt0axQsRXQZgpxBiBRFdmHJ1kuQ8IcQWIjoRwF+IaI39YBxjO24NeAuAAbbv/c3fapkdRHQyAJj/d5q/10RfEFE7GMJ3nhDiSfPnmm6zhRBiP4DFMKbfPYjIUmDs7Sq22TzeHcCeZGsamskAphHRRgCPwDBD3IXabjOEEFvM/zthPGjPRsxjO24BvBzAMHP1tD2AKwH8IeYy0+YPAGaan2fCsJNav19jrp5OBHDANrWpCshQde8DsFoIcYftUC23+QRT8wURdYRh814NQxBfbiZzttnqi8sBLBKmkbBaEELcLIToL4QYDOOeXSSE+DxquM1E1JmIulqfAVwM4G3EPbYTMGx/EsB7MOxm307b0B5x2x4GsA1AMwwb0LUwbF/PAVgL4FkAvcy0BMMjZD2AtwCMT7v+Gu09D4ad7E0Ab5h/n6zxNo8G8LrZ5rcB/Jv5+1AArwJYB2A+gA7m7w3m93Xm8aFptyFk+y8EsKDW22y2baX5944lq+Ie27wVmWEYJiV4JxzDMExKsABmGIZJCRbADMMwKcECmGEYJiVYADMMw6QEC2AmsxDRt80IZG+aEarOIaJ/JqJOEudKpWOYNGE3NCaTENEkAHcAuFAIcYyI+gBoD2ApDJ9L39eDm7u4AtMxTJqwBsxklZMB7BZCHAMAU5BeDqAfgMVEtBgAiOgXRPSaI1bvDS7pLiail4nor0Q034xnwTCpwhowk0lMAfkSjBCQzwJ4VAjxglOzJaJeQoi9Zuzp5wDcIIR4057O1J6fBPAJIcRhIvomjF1c/55C0ximCEdDYzKJEKKRiMYBOB/AFACPkvsbVf7ODB1YD0NrPh3GtmE7E83fl5iR29oDeDmuujOMLCyAmcwijDCQzwN4nojeQikoCgDjFTkAvg5gghBiHxE9ACMugRMC8BchxIx4a8wwarANmMkkRHQqEQ2z/XQmgA8AHALQ1fytG4DDAA4QUV8Yr4mxsKd7BcBkIjrFzLszEQ2PsfoMIwVrwExW6QLgZ2YoyBYYkbZmA5gB4M9EtFUIMYWIXgewBsbbCZbYzp/rSDcLwMNE1ME8/h0YUfoYJjV4EY5hGCYl2ATBMAyTEiyAGYZhUoIFMMMwTEqwAGYYhkkJFsAMwzApwQKYYRgmJVgAMwzDpMT/BwzvoWBPYFGEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# split our data into training and validation sets (50/50 split)\n",
    "X_train = rolls[:rolls.shape[0] // 2]\n",
    "X_validate = rolls[rolls.shape[0] // 2:]\n",
    "\n",
    "# check optimal score\n",
    "gen_score = gen_model.score(X_validate)\n",
    "\n",
    "best_score = best_model = None\n",
    "n_fits = 10\n",
    "\n",
    "for idx in range(n_fits):\n",
    "    model = hmm.CategoricalHMM(n_components=2, random_state=idx, init_params='se')\n",
    "\n",
    "    model.transmat_ = np.array([[0.95, 0.05],\n",
    "                                [0.1, 0.9]])\n",
    "    \n",
    "    model.fit(X_train)\n",
    "    score = model.score(X_validate)\n",
    "    \n",
    "    print(f'Model #{idx}\\tScore: {score}')\n",
    "    \n",
    "    if best_score is None or score > best_score:\n",
    "        best_model = model\n",
    "        best_score = score\n",
    "\n",
    "print(f'Generated score: {gen_score}\\nBest score:      {best_score}')\n",
    "\n",
    "# predict the most likely sequence of states given the model\n",
    "states = best_model.predict(rolls)\n",
    "\n",
    "# plot our recovered states compared to generated\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(gen_states[:500], label='generated')\n",
    "ax.plot(states[:500] + 1.5, label='recovered')\n",
    "ax.set_yticks([])\n",
    "ax.set_title('States compared to generated')\n",
    "ax.set_xlabel('Time (# rolls)')\n",
    "ax.set_xlabel('State')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c070fc60",
   "metadata": {},
   "source": [
    "Let’s check our learned transition probabilities and see if they match."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "969e7223",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Transmission Matrix Generated:\n",
      "[[0.95 0.05]\n",
      " [0.1  0.9 ]]\n",
      "\n",
      "Transmission Matrix Recovered:\n",
      "[[0.939 0.061]\n",
      " [0.101 0.899]]\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(f'Transmission Matrix Generated:\\n{gen_model.transmat_.round(3)}\\n\\n'\n",
    "      f'Transmission Matrix Recovered:\\n{best_model.transmat_.round(3)}\\n\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4393e7c8",
   "metadata": {},
   "source": [
    "Finally, let’s see if we can tell how the die is loaded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "760cd376",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Emission Matrix Generated:\n",
      "[[0.167 0.167 0.167 0.167 0.167 0.167]\n",
      " [0.1   0.1   0.1   0.1   0.1   0.5  ]]\n",
      "\n",
      "Emission Matrix Recovered:\n",
      "[[0.16  0.17  0.168 0.173 0.171 0.157]\n",
      " [0.116 0.106 0.112 0.092 0.1   0.474]]\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(f'Emission Matrix Generated:\\n{gen_model.emissionprob_.round(3)}\\n\\n'\n",
    "      f'Emission Matrix Recovered:\\n{best_model.emissionprob_.round(3)}\\n\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e44560c",
   "metadata": {},
   "source": [
    "In this case, we were able to get very good estimates of the transition and emission matrices, but decoding the states was imperfect. That’s because the decoding algorithm is greedy and picks the most likely series of states which isn’t necessarily what happens in real life. Even so, our model could tell us when to watch for the loaded die and we’d have a better chance at catching them red-handed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db632a17",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\"><b>Task</b>: Modify the code above for a hypothetical coffee --> water --> alcohol example.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77739873",
   "metadata": {},
   "source": [
    "## About this notebook\n",
    "This notebook is largely based on the examples in [Geeks for Geeks](https://www.geeksforgeeks.org/hidden-markov-model-in-machine-learning/) and the [hmmlearn tutorial](https://hmmlearn.readthedocs.io/en/latest/auto_examples/plot_casino.html)."
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}