Homework about Hidden Markov Models

Timur Khuzin HW HMM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework Program HMM Timur Khuzin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function that read parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All lines started with `#` are ignored\n",
    "\n",
    "First line of file: M,K,L tab separated<br>\n",
    "Than M lines of M numbers (MxM transition matrix)<br>\n",
    "Row index is old state, column index is next state<br>\n",
    "Than M lines of K numbers (MxK emisstion matrix)<br>\n",
    "Each row of emission matrix is probabilities of K emissions.<br>\n",
    "Than M numbers in line (Vector of M — initial distribution)<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "class HMM_Params:\n",
    "    def __init__(self, state_num, emission_num, to_observe):\n",
    "        self.state_num = state_num\n",
    "        self.emission_num = emission_num\n",
    "        self.observation_length = to_observe\n",
    "        self.transition_matrix = np.zeros((state_num, state_num), dtype='float64')\n",
    "        self.emission_matrix = np.zeros((state_num, emission_num), dtype='float64')\n",
    "        self.beginning_distribution = np.zeros((state_num), dtype='float64')\n",
    "\n",
    "def read_params(filename):\n",
    "    def read_non_comment(f):\n",
    "        s = f.readline()\n",
    "        while(s[0]=='#'):\n",
    "            s = f.readline()\n",
    "        return s\n",
    "    with open(filename, 'r') as f:\n",
    "        M, K, L = map(int, read_non_comment(f).strip().split())\n",
    "        params = HMM_Params(M, K, L)\n",
    "        \n",
    "        for i in range(M):\n",
    "            params.transition_matrix[i,:] =  [float(x) for x in read_non_comment(f).strip().split()]\n",
    "            params.transition_matrix[i,:] /= np.sum(params.transition_matrix[i,:])\n",
    "            \n",
    "        for i in range(M):\n",
    "            params.emission_matrix[i,:] = [float(x) for x in read_non_comment(f).strip().split()]\n",
    "            params.emission_matrix[i,:]/= np.sum(params.emission_matrix[i,:])\n",
    "        \n",
    "        params.beginning_distribution[:] = np.array([float(x) for x in read_non_comment(f).strip().split()])\n",
    "        params.beginning_distribution[:] /= np.sum(params.beginning_distribution[:])\n",
    "        \n",
    "        return params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HMM generator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.random import choice as npchoice\n",
    "\n",
    "# params is HMM_Params\n",
    "# returns sequence of tuples: (pi, x)\n",
    "def generate_random_HMM(params):\n",
    "    current_pos = npchoice(params.state_num,1, p=params.beginning_distribution)[0]\n",
    "    current_emission = npchoice(params.emission_num, 1, p=params.emission_matrix[current_pos,:])[0]\n",
    "    result = []\n",
    "    while len(result)<params.observation_length:\n",
    "        result.append((current_pos, current_emission))\n",
    "        current_pos = npchoice(params.state_num, 1, p=params.transition_matrix[current_pos,:])[0]\n",
    "        current_emission = npchoice(params.emission_num, 1, p=params.emission_matrix[current_pos,:])[0]\n",
    "    return result\n",
    "\n",
    "def beautify_HMM_sequence(hmm_seq):\n",
    "    def beautify_pair(p):\n",
    "        return 'p%d'%p[0], 'e%d'%p[1]\n",
    "    return [beautify_pair(x) for x in hmm_seq]\n",
    "\n",
    "def print_HMM_to_file(hmm, filename):\n",
    "    hmm = beautify_HMM_sequence(hmm)\n",
    "    with open(filename, 'w') as f:\n",
    "        f.write('\\n'.join(map(lambda p:'%s\\t%s'%p, hmm)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = read_params('test.txt')\n",
    "hmm = generate_random_HMM(params)\n",
    "print_HMM_to_file(hmm, 'generated_path.txt')\n",
    "#beautify_HMM_sequence(hmm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# input is HMM_Params and observation list (with ints inside)\n",
    "# output is hidden states in integer list\n",
    "def viterbi_solver(params, observations):\n",
    "    # matrix LxM where stored previous states for each state in row\n",
    "    parents = np.zeros((params.observation_length, params.state_num), dtype='float64')\n",
    "    # current probabilities\n",
    "    current_iter = np.log(params.beginning_distribution) + np.log(params.emission_matrix[:,observations[0]])\n",
    "    # probabilities of next step\n",
    "    next_generation = np.zeros(params.state_num,dtype='float64')\n",
    "    for i in range(1, params.observation_length):\n",
    "        for state in range(params.state_num):\n",
    "            probs = current_iter + np.log(params.transition_matrix[:,state])\n",
    "            prev = probs.argmax()\n",
    "            next_generation[state] = probs[prev] + np.log(params.emission_matrix[state, observations[i]])\n",
    "            parents[i,state] = prev\n",
    "        # avoid garbage collection\n",
    "        current_iter, next_generation = next_generation, current_iter\n",
    "    \n",
    "    path = np.ones(params.observation_length, dtype='int64')\n",
    "    path[-1] = current_iter.argmax()\n",
    "    for i in range(params.observation_length - 2, -1, -1):\n",
    "        p = parents[i+1, path[i+1]]\n",
    "        path[i] = p\n",
    "    return path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy: 0.755000\n"
     ]
    }
   ],
   "source": [
    "# checking viterbi\n",
    "params = read_params('test.txt')\n",
    "data = [x.strip() for x in open('generated_path.txt', 'r')]\n",
    "observations = [int(x.split()[1][1:]) for x in data]\n",
    "actual_states = [int(x.split()[0][1:]) for x in data]\n",
    "most_probable = viterbi_solver(params, observations)\n",
    "with open('viterbi.txt', 'w') as f:\n",
    "    f.write(\"RPath\\tObserv\\tViterbi\\n\")\n",
    "    p = map(lambda x: '%s\\tp%d'%x, zip(data, most_probable))\n",
    "    p = list(p)\n",
    "    f.write('\\n'.join(p))\n",
    "same_indices = [i for i in range(len(data)) if actual_states[i]==most_probable[i]]\n",
    "print(\"Accuracy: %f\"%(len(same_indices)/len(data)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import logsumexp\n",
    "\n",
    "# returns the matrix for posterior\n",
    "def forward_algorithm(params, observations):\n",
    "    observe_len = len(observations)\n",
    "    prob_matrix = np.zeros((observe_len+1, params.state_num), dtype='float64')\n",
    "    # initial state\n",
    "    # {1, 0,...,0} leads to NANs, but \n",
    "    prob_matrix[-1,:] = np.log(1/params.state_num) \n",
    "    for i in range(observe_len):\n",
    "        for nxt in range(params.state_num):\n",
    "            probs = np.log(params.transition_matrix[:,nxt]) + prob_matrix[i-1,:]\n",
    "            prob_matrix[i, nxt] = np.log(params.emission_matrix[nxt, observations[i]]) + logsumexp(probs)\n",
    "    return prob_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "matr = forward_algorithm(params, observations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# returns the matrix for posterior\n",
    "def backward_algorithm(params, observations):\n",
    "    observe_len = len(observations)\n",
    "    matrix = np.zeros((observe_len, params.state_num), dtype='float64')\n",
    "    matrix[-1,:] = 1\n",
    "    for i in range(observe_len-2, -1, -1):\n",
    "        for prev in range(params.state_num):\n",
    "            probs = np.log(params.transition_matrix[prev,:])+\\\n",
    "                    np.log(params.emission_matrix[:,observations[i+1]])+\\\n",
    "                    matrix[i+1,:]\n",
    "            matrix[i, prev] = logsumexp(probs)\n",
    "    return matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-135.82776627, -135.30817911],\n",
       "       [-135.33918632, -134.81965626],\n",
       "       [-134.85045917, -134.33115293],\n",
       "       [-134.36115517, -133.8427262 ],\n",
       "       [-133.86959477, -133.35459986],\n",
       "       [-133.36926906, -132.86765287],\n",
       "       [-132.83578423, -132.3853528 ],\n",
       "       [-132.18823033, -131.92163045],\n",
       "       [-131.24450407, -131.53678022],\n",
       "       [-130.71381007, -130.22118283]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "backward_algorithm(params, observations)[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# returns probabilities for each state in each iteration\n",
    "# iteration is row, state is column\n",
    "def posterior_decoding(params, observations):\n",
    "    forward_matr = forward_algorithm(params, observations)\n",
    "    backward_matr = backward_algorithm(params, observations)\n",
    "    probabilites = forward_matr[:backward_matr.shape[0],:] + backward_matr\n",
    "    total_probability = logsumexp(probabilites, axis=1)\n",
    "    out_probs = np.zeros((len(observations), params.state_num), dtype='float64')\n",
    "    for i in range(len(observations)):\n",
    "        out_probs[i, :] = probabilites[i,:] - total_probability[i]\n",
    "    out_probs = np.exp(out_probs)\n",
    "    return out_probs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1.])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "posterior_decoding(params, np.array(observations))\n",
    "np.sum(posterior_decoding(params, np.array(observations)), axis=1)[:20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unfair casino"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "casino_params = read_params('casino_in.txt')\n",
    "casino_hmm = generate_random_HMM(casino_params)\n",
    "print_HMM_to_file(casino_hmm, 'generated_casino_path.txt')\n",
    "hidden_states = np.array([x[0] for x in casino_hmm], dtype=int)\n",
    "casino_observes = np.array([x[1] for x in casino_hmm], dtype=int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### viterbi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy: 0.790000\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:7: RuntimeWarning: divide by zero encountered in log\n",
      "  import sys\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x178c2356320>]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXv0XlV557/PeRMCyi2QYGJCCJdUAbERf0WRjlYERDpjaEtbWKs17dKVLkfstNZWGGfUWu3CdqZ2tYvWxpqK2CUUWoe0pkVu6rQK8qOGS2AgIWAJCRDuAnLJe57549z2e959bvvc9n7P81nrt37ve27v3mfvs5+zn9smZoYgCIIgRHh9F0AQBEGwCxEMgiAIwgQiGARBEIQJRDAIgiAIE4hgEARBECYQwSAIgiBMIIJBEARBmEAEgyAIgjCBCAZBEARhggV9F8CEJUuW8OrVq/suhiAIglPcdtttjzPz0qLjnBQMq1evxvz8fN/FEARBcAoi+mGZ40SVJAiCIEwggkEQBEGYQASDIAiCMIEIBkEQBGECEQyCIAjCBI0IBiLaRESPEdFdGfuJiP6MiHYQ0R1EdLKybz0RbQ//1jdRHkEQBMGcpmYMXwZwds7+9wBYE/5tAPCXAEBEhwH4JIC3ADgFwCeJaHFDZRIEQRAMaCSOgZm/Q0Srcw5ZB+ArHKwjejMRHUpEywH8DIDrmPlJACCi6xAImK81US4dd+56Btfd/Yh23wmvPQRnv2FZWz/tPC/v8/F/tj6M805eCc8jo2tcd/ej+MmVh+CIg/dvuHR2cO22R/CmVYfiiIOard99j/4I/3T77tLHLz/0AFxwyqrKv8PM+Id/fxjnnLQcB+w3qnx+Uzzw+PP4+g8eBsKlh4894kCsW7sCALD3Ry/ha9//D+wb+9pzDz5gIX79tKMxMuyjVfm3HY/jlp1PaPe99djD8bZjl3RSjibpKsBtBYCHlO+7wm1Z26cgog0IZhtYtap6h4/48xu345t3PwpK9RlmYMmB+4lgyOG79z+O37v6Dhy/7GCctPKQyuf7PuM3Lp/Hb5/xE/jwu9a0UMJ+eWXs44NfvQ2/c9br8KF3Htfotb/4nZ246rZdU/1WR7SM+zlvWI5DXrWw0u888Pjz+J2rbsf+C0f42TcuNyhpM1z+vR9i0789AKKgPiOPYsHwj7fvxp9cdx8AaJ9jADjtuCU4fvnBnZT1D7fcg227n9WW5dv37cU1F/50J+Vokq4Eg647c8726Y3MGwFsBIC5uTntMWXY5zNOWnEI/vHDk431qc3bgjcUIZNXxsFtf8XXv6kVMWaGz8EAOouM/fbqt89nrDrsVfjO772z8NjLv/cg/uc124zaKW7jntton+9j8asW4gefOAt/ev19+NPrt4OZQUTYF9Zr2++/G69eNDmEXX/3o/jAV+axb2w8RFQv65hx9onL8IVfffPE9g9cNo/dT/+4s3I0SVdeSbsAHKl8Xwlgd8721hj7DN0Mkyh4oxWyGYf3hw1vU3T+2PQCluOH9WqjH2X1Wx0UvrqalCNuo56fhaC+QT2i/0nZMLFdxQtHtC772Jg5/t2JslDSJ1yjK8GwGcD7Qu+ktwJ4hpn3ALgWwFlEtDg0Op8VbmsNn1mrHx8ROduIXcHx/TG7T9Hpsyp//Rbrl9VvdUS6dZNyxMKt52fBZ8T1TdcnKpt+MKaJY7rAZ9YKqZHn7pjSiCqJiL6GwJC8hIh2IfA0WggAzPwFAFsAnANgB4AXAPx6uO9JIvoDALeGl/p0ZIhui6xG9Dya2TfZpojuj/GMocU3ahtoc0aU1W91RPLDpBzWCAZlhhRVOz0j0z7HNWZLpvh+xphC1PvMy5SmvJIuKNjPAD6UsW8TgE1NlKMMY58xymhEQ9X5YIg6uWlft0VN0RbRYNSWKknXb3XUGRzT6pq+GHNS31FalRQKCN39iGYXXfaxMbPWA8rzyNnZ8eAin4Mp6vT2kdf/W5LtRLeHDe8Tcz3BYjttvm2rqpUiEtWLyYwh+t+3Kok1qqTJ/qO7H4kqqYNChvi+fvYyEhuDO+RO+xxtxK6Ijc81z3f1YSki6j9tvG2rqpUi0sbaSr9jlSopqEdiTFf36c/zUmqnLgjUfLqyuKtKGpxgyJz2EYHZ/G14CIxrDhrJwDmb9zgeuFroQ1n9VodXY8Zgi7pvzMlMYZSymeTdi15USX6OKsnRvj44weBzxrSvhifHUIiFpuE9SrySZvMmt65KKmljGNVQp6TVNX2hvoVPq5KyDfF1hKIpWWq+wNOxs2I0yvAEQ8Y0NPbkcLUlOyBSkYgqSU+bb9vVVEmT5an2O8nv9YlelcRT+9L0566qKYvnbszO4ARD3rQPmN1BqwnqqpL8WVcltThjyOq3Orwa6hRfUdf0iVrf9Gx+7CNblRTbV9ovY0S+p6ObfX1wgiEzGKWHNw3XiFRJprco0cE3VCDLiD162jA+V4hjiPqySTvVFf5NoarO0nEZWW/oQOJx2PmMQadKcjjATQRDSB1PjqFQ2ytJAtyMqRTgViMtRJuxGJXKoaSZSMdl5EWBS4BbMwxQMOinoZ4Ynwup6+NuiytkW3CL9cvqtzrq6NnbTOtRrRxKgJvG+JwV7NeHE0nmmBJ6OrrI8ASDz9rUxbH/c99PhMXE98ZYlRS9UTdUIMtoc0Y0zui3OpqJfO7fxkCxKmk6iR5lGp/D8ztOopc1pvRtqzFlcIIhywc69n92tCG7IM6VZCgZhqNKav7afoU4hjq+/LbM6tT6ph1DfJ8xyhi5+lIlZaXn6FvAmjI4wZA1De3Dzc01kgfT8PwWA8BsoM04jTz1SZo6aSGsEQw+pnIlqarMYlVSt8bnLPW0q119eIIhYxqavGl0XSJ3iN7CTPv6rLurxnEaraiSstUnaeqkhbApiV5U3XRcRrAv3/jcVR9jDhZn0o8p7mogBicYxhnT0FEPC3y4RhzgJsZnLW2m/MhTn6Spo0qyJTrd94tUSXZEPke3WKtKEq8kd8ic9vWgm3SN6GGrH/ncUIEso83ssVVsDE3kSur7OZiwMaST6OV4aNVJB2JCdI91QjtqBxfzrzUiGIjobCK6l4h2ENFFmv2fJ6Kt4d99RPS0sm+s7NvcRHny8DOmoWJjKCYWDBL5rGXcchK98qqkGoLBlshnRT2Tns1neQEB3ae2iX4nb0xxsb/XXqiHiEYALgVwJoI1nG8los3MfHd0DDP/tnL8hwG8SbnEj5l5bd1ylCUrfL2PrIyuUXfNZ1ty/bdF27mSyhqf66SFYLZkxuBznFU1Lejy7kX3qqRoxpDv6djIimgd0sSM4RQAO5h5JzO/DOAKAOtyjr8AwNca+F0jJMDNnJphDDOfRM+aALcaaSHGlqQt0auSeGpfmlHHKuE8G0NUbhe7exOCYQWAh5Tvu8JtUxDRUQCOBnCjsnl/IponopuJ6NwGypNLYYCbi63YEXWNx36Lb9Q20Gaeoc4C3GxRJSkBbunZfH6AW/SW3kEhoaqSdGWZPMYlmpjh6Foo606cD+BqZh4r21Yx824iOgbAjUR0JzPfP/UjRBsAbACAVatWGRdWXUtWJb2urDBNY6qkGXUJblWVVCGOoU6wpjWqJKW+6biMYMagPy+eLXU1Y/DLqZJco4kZwy4ARyrfVwLYnXHs+UipkZh5d/h/J4BvYdL+oB63kZnnmHlu6dKlxoXNC0aJ9gt6anslzbi7auLq2fy1K3kl1fDMsUXdp6rO0rN5mwLc8mwMsSrJwRehJgTDrQDWENHRRLQfgsF/yruIiF4HYDGA7ynbFhPRovDzEgCnAbg7fW6TSICbOXGAW11V0owKBusC3GrlSqp8aqOoKt9pVVKJALeO+lj0O1kBbl2WpUlqq5KYeR8RXQjgWgAjAJuYeRsRfRrAPDNHQuICAFfw5KhyPIC/IiIfgZC6RPVmaoOsaeiohzzurhEbn41VSbMdx9BmAB/nqE/S1Hlr5riN+54xZAe4cY4hvmuDb/Q7eZ6OLo4pjXhRMfMWAFtS2z6R+v4pzXnfBXBSE2Uoy5j1udOp4zcNF6mdRM+S4Km2SFY/a/7aWf1WRx3/eWuMz0p90+6qYz9noZ6e4hh05UkvSeoSg4p8ZmYw6xdV79rNzUWSRVwMzw9v7awa+Mctrpc8zlgMRkcde5ktQYi+nzyno5SaN09Idh2PFAsGMT67S+xznNOIMzpmNUJd4/Os50pqV5VUPo6hTlqIxI5U/dwmUVW+lNLVc44hnohA1J0qLFeV1HF6jiYZlGDIn/ZNHiNMUzeJni0eL23R5tt2nvokTR11StTGfT8H6gwpfmlTjM95syePqHPjs6cZSaMiuqiFGJRg8Dln2ie5kgpRjX91zu970GmLNgXfmLPXOU7TiCrJAuNzVI/0bH7M+mc4Ishq2noRg7LEL5uzlWZnkIJh1jwIuiJRJRm6q9YULLaTBGA1f22uEuBWw15WN1FiU/gMJcAt2DahSsq5FZ7XpSqpOMDNxTFlUIIhT7pTDU+OoVA38nmsGA9nkcQ437PxuUZaCJvWfI7G2nSKj1KqpK6Mz1w8pohgsJzIqyHPg8DFRuyKunEIM69KakkNE60SVl6VFPw3mzEE/7vKNZRZDl+nSlIEQ5EqqSsbQ54qqUaW274ZlmCIVUnT+9IuccI00b0xViVZ4vHSFm15JeV5vuio85JjjyppOleSOmPNuxddrrUct432ZTP47+LL5qAEQ+JBoJv2TR4jTJPoeM3On/U4hrpxHlkk6opyx9dJC2FLBlzV2J6shKbuyz7XI7sC3Pq+lyYMSjD4edM+h5fh64ro3pjeoVlPotfWQkR53nQ66qSFsKWNfE7qkTY++wVR4COPOk+il+fp6GJ3H5ZgyJn2eQ7rA7simcrXUyX1Pei0RVuG29g2VjGJnkk52pr1mJQjqkc6Jb5fYHwm6kEw6JLopZYkdYlBCYa8KXl6XVlhmsggKXEMetqyMUR9smoSPSPBoKhr+mSsSaIXvZCo+3SMuvRK8pPfTOPyms+DEgx5qqRk+u1eI3ZFrEqqHfncWJGsoq3ssXlvpTrqpIWwQZWUzmmWHmDVPEo6AlVS++UEVFWSvhyAm2PKsARD3rTPYeneFXUHdnUx91mkrXQSeS80WZimhWgzFqN0GcKfTifRi2asgY0h+3yi7spf5mXTxTFlUIIhaqDcZfgcbMSuiG0MhufboqZoi7YEX16/zcI0LYQNKTGS+gbfo7dxVlSRuaokr/tcSbl2Swf7+6AEQ/wmkrO0p4Nt2BnRvZEkenraMq7n9dssTNNCJKnDK5/aGGlPnylVEhfnSupOlRT8z/d07KYsTdKIYCCis4noXiLaQUQXafb/GhHtJaKt4d8HlH3riWh7+Le+ifJkkaiSpve5vAxfV9SNY0gWmm+qRHbRVuRzXr/NwjQtRNxGPT4HaZVvel0DO1VJ0/u6XjSoSWqv4EZEIwCXAjgTwC4AtxLRZs0SnVcy84Wpcw8D8EkAcwg0FLeF5z5Vt1w64ilqbvi6e43YFYkqyXTGEP6fUeHbWoBbTr/NwjQtxFhR1/RFur7puIyxn59QcOR16ZWUo0oa+EI9pwDYwcw7mfllAFcAWFfy3HcDuI6ZnwyFwXUAzm6gTFryglHSLnHCNFxzxmCDx0ub2BLgFh1rFOBmgbovrTpLv3mrKbl1eJbEMSQBbu719yYEwwoADynfd4Xb0vwCEd1BRFcT0ZEVz22EvEAhlz0IuiIZ2M3OVwWLiw9LEa2pkioGuAXHmvVljoVb5VMbI62eSed+Kgpws0UwuBw024Rg0LVQulX+EcBqZn4jgOsBXFbh3OBAog1ENE9E83v37jUqaF6gUNolTphmXDOJnjpQzaIAVpMENin4qga4BccaqpIsyJWU9vSJ01f7ieAt9ErqOsBNq4WIjnGvrzchGHYBOFL5vhLAbvUAZn6CmV8Kv34RwJvLnqtcYyMzzzHz3NKlS40KWiZ8fRbfZJuiKVUSMJtBbn5L9asa4BYd62qAm66+qqBT8yjp8HoIcNO+bDqsnm5CMNwKYA0RHU1E+wE4H8Bm9QAiWq58fS+Ae8LP1wI4i4gWE9FiAGeF21phVoNRuqJuriT1tFm0M6gqgyb7kXGAWx2vpD4D3DSqM9UFVc2jpMOj7vpX9DuUN6Y42NdreyUx8z4iuhDBgD4CsImZtxHRpwHMM/NmAL9JRO8FsA/AkwB+LTz3SSL6AwTCBQA+zcxP1i1TFnlJ9NLrygrT1F2aUx2oZlEwTM4YGhQMOf02C9O0EDakLdG9hasuqH6JXEldC4a8XEkujim1BQMAMPMWAFtS2z6hfL4YwMUZ524CsKmJchQRdXrdi1e0bRYHrKaIOrjpHZp5G0NLgiEv538Wpr78kY3NhshnSquSfI5Xs9O9oUd0urRnrkNL8N/FFDADi3zOlu4Sx1BMXVdGVQU1i0FubQk+ExuDqfHZClWS5jmNVEnx7CnXxtBd/4rVfDk2BhfHlGEKBlnz2YjaqqSW3qhtQX3+2zA+V82VVE+V1KdgCP6r9aXQbpBn7I3oY6GevFxJLvb1QQkG3RQ1Iu0SJ0wTu2Oanq+c6KJBrgi17zTZj/IWnM/CVJVkw5oZOpVvNNjnPcMRppllTRjnzOY8h182ByUYit68uszK6CLxvTG8R20NnLag9p0m+5FJ5LOpL3+kgrHD+JzUN7IblJk9eURWpN1O1NOdFKVRhiUYclZbirbP4HjVGEnKB7PzJ72SGiiQZbTulVTRXdWkDDbFMaj19cIZQ5l70e1CPclvponsDjJjsJyo02f1qS6zMrpIokoynDEMSpXU3HVNvJJMBYNdqqSU8dnP9yyMME0HYkJe24iNwRGiBzdXlSSCIZO6abcn3qhn8D6PWxJ8iedLF6okC2YMmjQTHgX3tOgZDo7tIVeSzqHFYU/HYQmGnGkfIKqkIpJBw+z8QQW4NequGvyv4pVkmhaibqLEJtB5HnleYDcoY2PoxSsp1/jcSVEaZVCCIfEg0O+PXOIEPUmAm6kqSTHOuvi0FOC3JPiK+q0O07QQfkvrVlchUflOBrj5zNp9aSTArT6DEgxFOWdElZTPOJEMRrRlnLWF1gLcDNxVR4aDow3qvlhdpBqfiTDmYgcSoJ8kerkBbg729WEJhhLuqrM4YDVFnA/f2F1VvVYTJbILmwLcPMO+bIPw1qnOohlQqQC3LpPoaYRYhBifHaEoUIg6NFq5SJORz7M4M2s/V1K1hXpMPKMmZj09PQs6z6MoNqF0gFtXqqS8ADeHg2YHJRiifp7l3RG5xAl66ibRa0sHbwttxzFUzZVkNmNIPvfVRFG+plHKxhAk0cPUvjSmy5qa4OeMKS5nbB6UYCgy4kUucYKeunl0JvXXjRTJKlpPolfhaTVNC2GDg0D8nKYin31W92Wf32UcQ97aEOm1ql1iWIIhRx8IJC5xwjTqfTFXJamfZ+8+tyX4ivqtDtO0EDapktQZkudhIldS3uypy9Q2ecuMEpGzno6DEgyseRNREeNzNk3cl5lXJU0Y19uYMbS/UM/EC0BPs7pYXaTUN/KyitVMBQFuXS2n6TPnCylH7ZaNCAYiOpuI7iWiHUR0kWb/R4jobiK6g4huIKKjlH1jItoa/m1On9skRW8bkUucME0TKbNtcIVsk9aT6FU0Ppu5qyaf+58xJNsiL6s8Y298bIfG50CVVFSWTorSKLVXcCOiEYBLAZwJYBeAW4loMzPfrRz2AwBzzPwCEX0QwB8B+OVw34+ZeW3dcpQhGvQzVUmSKykT9W3YWJXUkg7eFtpLux387zKJXvC7PdsYKG1jqKBK6jDALXf24rk5O25ixnAKgB3MvJOZXwZwBYB16gHMfBMzvxB+vRnAygZ+tzKJKkm/X1RJ2aj3pU7kc7KEahOlsou26mdifDb2SvKTOnSljkmjUxclqiRM7UsTqJJaLWJMoErK3j/qMAV4kzQhGFYAeEj5vivclsX7Afyz8n1/IponopuJ6Nysk4hoQ3jc/N69e40KWkqV5GAjdsGkKsnsGj4DC8PIpFkUwGOlfn1HPpv2ZZ85qUNvqqTg/7TxuVym2S69C33mXNtPl4sGNUltVRIA3V3R3gki+hUAcwDeoWxexcy7iegYADcS0Z3MfP/UBZk3AtgIAHNzc0Z3eqx5E1HpMiujazTileQzFnqElzG7qqSofm3kSuokiZ4FbZTUN9nmEWHf2Ne6sqbpVpXExek5HOzrTcwYdgE4Uvm+EsDu9EFEdAaAjwN4LzO/FG1n5t3h/50AvgXgTQ2USUsc4JaTK8nBNuyEyftidpOYGQtmeMbgt1Q/owA3QzdJZsR16DvATa1vpBrTBb+l6TrArUhIuTimNCEYbgWwhoiOJqL9AJwPYMK7iIjeBOCvEAiFx5Tti4loUfh5CYDTAKhG60YpmoZ2GRjjGuMmZgzMWDhyN39MEWM/qV87qqTy55iqktQ26m3GoFGdUegxmJfNNKJTVVJOgFvXZWmS2qokZt5HRBcCuBbACMAmZt5GRJ8GMM/MmwH8MYADAVwV5jj5D2Z+L4DjAfwVEfkIhNQlKW+mRokDhbJUSWJ8zkQ1RJreo7EPLAgtqLMY+cyc1K/JbmSaRM+kDGOfkzayKIneKPQYLGOI7zJ2wOcCVVKHMRVN0oSNAcy8BcCW1LZPKJ/PyDjvuwBOaqIMZYgaKCsBl9gYslHfekxvETNj4QJ3UxEXMVbr16i7an6/1WE6+2VGXIfeBIMmiV6kSipjiKfQK4mZK90zE8YFv+GqQ8ugIp/zwtcB8xz2Q2BClVTjGgvjGcPs3WdfqV8bAW5VZgymaSHGrNShp1mdzthO4bNZxhAfr4PQQR/z/YIxxXMzwG1QgsHnIqPVbPrXN4E6xhirkpixILYxNFEqu/CV+jWpPohVKxXdVU3KMPbVNupLlTRtYI7UQ2UM8V1mNfW5OMDNRVXSsASDEryjwzTx2BCYePsyViUlOvhZVSUtaOFtW7c+QREmKoxoAEvq0LcqKe2VVM4QnwQZtl/+MRePKS729UEJhnGZaZ+DjdgFEzYG02soXjuzKIB9H4lHT5MzhgKnCR0mvvzR8X17jumcRCg0Phc5kADJTKMzVVJBEj0X1dODEgyFqiRy0+e4C5rwSlKjamfRyK/Wzw5VktnvJG1U7fym0NU3Nj6XSKKXqJI6EAwFdssuYyqaZGCCoWjaN5tvsk0wbiCJnq/or118iypi3FL9otlHNVVS9VlLNJD23UZROSgV+TwuKRgiFVQXLtFjv2iZUTf7+qAEQylVkoON2AWNeCUNaMbQdICbR9XcVeupkvptI93CRF647G6caTZXlRRep7MZQ/Z+sTE4QOG0T+IYMmliPePJJHqNFMsq1Po1HeBWxb4AmKkwonaN26i3GUPwfyLALUxfnbjuZp/fuSqpQK0lXkmW45cIRhHBoGfivtRRJXXoY941Y7V+DSfRqxqoZaRKCt/G26hDpXJoVGeRl1Wyr4wqqQOvJF8C3JynKBOiqJKymVQlmccx9K2maBPfb0+VVMXwDJh5w4ynZgyVTm8MrSrJm1yop+g5BrpUJeUbn11cFXJQgqE4GMVND4IuUMcY0wFDNT7PopG/zQA3E1VS1XKkjc+9B7ilFurxWa9mShMJjU4C3PwCIUUS4GY9xQFusxl41QTNrOCmLGQzg7dZnRE1nSupasqfyGunSjkiYd33Qj26ALfIu0eXRylNHODWhSqpTICbgy9BwxIMJXIlzaKKowmaWKjHV1I6u/gWVUQg+Jp/W2UD47NJWogkjqHfNtLNkCJVUpm8UV0an4vaxtWMzYMSDGMu9n+exXTQTdDE0p5qSmcX36KK8FtKWT1mrrRID2CWFiJq4zbSelQhqO/ktihdTbyCW4HBF+imj439/LYJYqNaL0bjDEowFC2qMfJmc8BqgsnOXT/yeRZVdu2pkqqt3gaYpYWYUiX1mCspXd8oXU2ZtNtehzOGMZdYZtTBvt6IYCCis4noXiLaQUQXafYvIqIrw/23ENFqZd/F4fZ7iejdTZQni0JVkqPTvi6YsDEYq5JUNUUTpbIHZg7WMmhNlVTtHBN1ShLH0Lcqafo59Sw1PjNzHFCnw1UX+NqCgYhGAC4F8B4AJwC4gIhOSB32fgBPMfNxAD4P4HPhuScgWAr0RABnA/iL8HqtUDTtI0cbsQuaSqLXd7qFtoiq04ZHT1G/1WGSFiJqkwUtJAKsgm6GFKWrKVqeV91nhyrJzYzNTcwYTgGwg5l3MvPLAK4AsC51zDoAl4WfrwbwLgp67joAVzDzS8z8AIAd4fVawS/Q1bqaCbEL1M5tMuhF58+qjSGdTqLpXEnVVUnJuWVJRz73mSspPfDHqqTIxlBg8AU6FAwzqEpqYmnPFQAeUr7vAvCWrGPCNaKfAXB4uP3m1LkrGiiTliJ/8JFHeO6lffjgV29rqwjO8uizLwKIQvyB517ah09t3obnX9pX6vzo2VjgEYiAf7pjN+579EeNlO1txy3Br771qKntd+56Bl/49v2lBdkhByzE7687EYsWTE5aX3xljE9esw3PvvhK5rlpwbD59t24Z8+zZauQy+0PPY0FFXVJUT//6FW3Y9GCcuc+F7ZlVIdN//oA/uWuRyr9LgD86qlH4W3HLom/XzX/EG78f4+VPn/b7me1qiRm4OrbdgEoih0I9v3hlntwyAELqxQ9+C2P8Junr8Hrlh0EALj0ph246+FntMc+/NSPsfrwV2dfiwg/fOKFRseUT/6XE7HskP0bu56OJgSDroXST2LWMWXODS5AtAHABgBYtWpVlfLFvPbQ/fGq/bI1Vaceezi+d/8TuH/vc0bXn3XmjlqMp154GQxg28PP4OrbdmHFoQfg1YvKaf9ev+wgzK0+DP/5ja/FvY8828h93vP0i7h7z7NawbDlrj34xp178BOvObDwOs+9uA+7n3kR7zt1NU547cET++579Ee4cv6hwroev/xgzB21GOectAw7HnuusX504P4L8PY1Syuds/bIxTjxtQdj11MvVDrvjSsPwemvPwL/d/tePPn8y5XrsHPv8zhw0YIJwfDl7z6IBx9/HisWH1DqGvsv9HDaccsntp1y9GFt9fA4AAAb0ElEQVR4/bKD4DPjjOOPwAELs9vhdcsOwk+uPASPP/cSHn/upUrlZwa2P/Ycjl92UCwY/vzG7Thg4QhLD1o0dfzyQ/fHO16X3TY/87ql+I8nn290THl5X/tuTk0Ihl0AjlS+rwSwO+OYXUS0AMAhAJ4seS4AgJk3AtgIAHNzc0Zzs8+ce1Lu/nefuAzvPnGZyaUHw8/9xb+BOXEb/F+/+JM49djDK12j6vF5fOTKrfj+g09q9/k+Y9ECD9/87XcUXueb2x7Bhstv084uotnAZ859A975+iMKr/WWY5qrnyknrTwE3/jN/2R8/hUbTjU677RLbpxSnYx9xtuOW4Ivvm/OuDynHbcE//Jbby917JGHvQrXXPjTRr8z9hnH/vctE666vg/88k+twkXveX3l6/3KW4/Cr2heWmynCRvDrQDWENHRRLQfAmPy5tQxmwGsDz+fB+BGDlweNgM4P/RaOhrAGgDfb6BMQktEU/ro2a8aeNV4eXLSmFTJSprnxROpqvN0yULASNMezNUWGeqTqInT2YSreoW5Tu0ZQ2gzuBDAtQBGADYx8zYi+jSAeWbeDOBLAC4noh0IZgrnh+duI6K/A3A3gH0APsTM47plEtqDEKTEKOMd0gV5C6FU8f/PC4pKFocxLOSA0LXHmBmeIwMrUWADUwWDifHfdZpQJYGZtwDYktr2CeXziwB+MePczwL4bBPlENqHwkjOcQnvkC7I8/rQebdkkQRFTe8rk9FTCPA07VHkDWgbamqcKD7FpfI3gSNyXLAFIgKD4+CnvgfLQLWVIxhKSgadCkG9DtC/EHSBkaY9dJHMNhMkvgs+lwmom0VEMAiVIAQ64+jB6fuBz8teWbT+hkpeCgnfkrq6gK49xgZJAPvE85KXAVtUpl0jgkGoBFEgGJK36H7LE6Qx0e/zC/LYqOTl1ymznKQQ4GnawzfI9dQnIyVaeaizRenqQiUIgSqpTDKzTspD2Xn3i5ImqkT10KWQiHTmVZfXHCKepj2q2HpswKPETpI4HjhUgQYQwSBUwvNCVVKJvPhdMKJs4/OYK6iSvOScNL4Yn0ujcwYY+66pkpIZw1AdD0QwCJUgRAumBN/7fpPKy4hbzfhcHMfg0uDWF1EWVJUqKj0bUNWTQ41hEcEgVIIoyFniW2KUy1tcqYo3TKJKyo58HthLoxGzoUpKZo629POuEcEgVILCyOd4it23KsnLziBalDRx8jrZcQxllpMUAkbetFeSX0GlZwOqC7QtKtOuEcEgVCKIfLbHKJe3TrduicgsykQ+uzS49YVuYZqi1NS2oQo3W/p514hgECoRuKuWy4vfTXmi3E16o3FpVVL4JOQl0ROvpGJ0gsHpALeBxrCIYBAq4cWqpOB732/R0RQ/602/tCop1/g8THWCCXpVklv3LvC8S6uS+ixR9wysukJdoiR6NgW4AVk5jsq/5VOeKskSIegCugC3MbNThvuJOAZL4nW6RgSDUIkoiZ4tutfo57Pe9Mu+6eWl3U4C3MzKOCQ8mr6HfoXUJDagLvFrSz/vGhEMQiWCJHr2BH0VqYCq5krSub6yqJJKo3MGqKLSswF1jY+hxrCIYBAqESTRY4wtCXDL8yYa+1xBlRSeozU+T/6WkA0phlsgclRwy3Cvrikx1BgWEQxCJeIker4dNoY4+Z3mTd9oBTedgLHEnuICI2/yHsZv3A6NrLpcSTJjqAARHUZE1xHR9vD/Ys0xa4noe0S0jYjuIKJfVvZ9mYgeIKKt4d/aOuUR2scL12Ow5YEZ5dkY/PIDUp4R25a1J1wgnaLExcy0wfKkKcEwsLav21wXAbiBmdcAuCH8nuYFAO9j5hMBnA3gT4noUGX/7zLz2vBva83yCC0TzRjGlhjlohmDVgVUwRsmX5VkR11dgFJJDV2MAVHXlHCx/E1QVzCsA3BZ+PkyAOemD2Dm+5h5e/h5N4DHACyt+btCT8RJ9CwZLPNyHLFBHIMuUC4WDANTJ5igrmUAIDbi9j2zrILqchu7KjtU/iaoKxhew8x7ACD8f0TewUR0CoD9ANyvbP5sqGL6PBEtqlkeoW2iJHqWPPD5cQzVk+jpjNguDm59kV44KZlZ9lQgA1SXWxdVYU2woOgAIroewDLNro9X+SEiWg7gcgDrmTkyFV4M4BEEwmIjgI8B+HTG+RsAbACAVatWVflpoUG8ML2qLUseenkqoArpnmOVVJ7x2aHBrS+IJu+hi2o4NY5hqIs0FQoGZj4jax8RPUpEy5l5TzjwP5Zx3MEAvgHgfzDzzcq194QfXyKivwHw0ZxybEQgPDA3N6fPmia0DgHhegyB/r7vB6ZQlVSyeNFsQJePz8XBrS/ScQwuxoB4igHdlnidrqk7QdoMYH34eT2Aa9IHENF+AL4O4CvMfFVq3/LwPyGwT9xVszxCy8TrMViSSjk3YrmSKik8R3MdFwe3vkh7JbkoVD1KbAu2qEy7pq5guATAmUS0HcCZ4XcQ0RwR/XV4zC8BeDuAX9O4pf4tEd0J4E4ASwB8pmZ5hJZRk+jZ8LAXBbhVXcFNf53JY4Rs0gFuSQyIO/dOXZ50qAFuhaqkPJj5CQDv0myfB/CB8PNXAXw14/zT6/y+0D2qKsmGgC8vZ8bAXD2OQeuVJDaG0ow8pFRJ4XaHRlY1dbjEMQhCGZTIZxseliRX0vS+cQXhlcwYpvdxuOBP3/YUF0jbGGxxUqiCp7jc2hLI2TUiGIRKxAMo27H4SmwbyFiPoaqNoa6tYugQTa7HYMuCTlVQXW4lwE0QShAl0eMKrqBtkqdKqrJyGBEFKcW1WVrtqKsLjJTMpICbK6CpSfSGGsMigkGoRDB4Rm/RfZdGdVed3jeumO55RNOrjwHRzMO4iINCHVQBN1dAU20MLqrCmsCh5hJsIEqiV3XQbYtowNG5mfoVPad0q48BweBggz3FBTxvMleSiwvdqC63tuQE6xoRDEIloiR6bI2NIUeVVPFNX7f6WHwdC4SgC4yIJjy7bMmpVQU1id5QY1hEMAgVIUWV1P/Dkhf5PPYbUiVZUlcX8FL3cOzgwOpNGJ/DbQNrfxEMQiWC55sx9u142KMy6G0D1YzGaioEFVvUZi4QDarxegYODqwjZebooo2kCQZWXaEuE6okC3pPokqa3lddlUTamYfPbg1sfZKkLw++JzaGvkpUHZ0qaWjtb8GjLbhEtB6DbXEMWfEHVYzGo5ThNMK3xAPLBdI5p6IB1qUZl+clLwgu5npqAhEMQiWiJHq2eOrkq5KqGY0DN8Xp7VVtFUMmnb7cyQA3mg5wG1r7i2AQKhEl0Zu1ADcgyqopqqQ6pNOXu+iu6nnJjCeqhw19vUtEMAiV8ZntC3DLiFiu5JWUYXy2JWGgC6RVSXHaapcEg+JyO9QEitLdhUoQAeBqeYjaJE6ilxH5XKWIXipldIQta0+4QFpQuxg5rBqfJbuqIJQgiHy2RzBQ6g1VpWoGWM/LNmLbUFcXSMeVxAFuDkmGkacIBkmiJwjFROsx2GKQjVdwyzA+Vw1wk8jneqSdAVxcAS2yowFifDaCiA4jouuIaHv4f3HGcWNl9bbNyvajieiW8Pwrw2VABYuJ4hhsyTiaLO05uZ2Z4XO1N7101G6E7w9PlWCKl2oPF3X0HrltI2mCujOGiwDcwMxrANwQftfxY2ZeG/69V9n+OQCfD89/CsD7a5ZHaJkoiZ4tGUez1mo2eaDzIp8HNi4Yk44rcTFX0oQqKawHDUy3Ure66wBcFn6+DMC5ZU+k4FXudABXm5wv9ISSdtuGt6isXEnJylvlrzUi0hqxfUvUZi4Q9Ykp461D98/zNKokC/p6l9QVDK9h5j0AEP4/IuO4/YlonohuJqJo8D8cwNPMvC/8vgvAiqwfIqIN4TXm9+7dW7PYgikESrySLHjYE1XSpGAwWXmLKMOILbmSSpOOK3ExclirShpY+y8oOoCIrgewTLPr4xV+ZxUz7yaiYwDcSER3AnhWc5wm7jTcwbwRwEYAmJubyzxOaBePEKiSfFjh2++l3lAjTFbeClYf06mShueVYkp64aSoWVwSDKoTQqxKcqf4jVAoGJj5jKx9RPQoES1n5j1EtBzAYxnX2B3+30lE3wLwJgB/D+BQIloQzhpWAthtUAehQyLj85gZCy2QDFmRzyZGz2zjM2M0sIHBlPTCSUlKjL5KVB2Ko/tZVEmGbAawPvy8HsA16QOIaDERLQo/LwFwGoC7OXg1uwnAeXnnC3YRJdGzJY4hDnBLjecmKoxg9bHp7aJKKk9WgJtLA6vqcuuijaQJ6gqGSwCcSUTbAZwZfgcRzRHRX4fHHA9gnohuRyAILmHmu8N9HwPwESLagcDm8KWa5RFaxguT6NmyeE3slTSlSqr+QI8IelWSz6JKKslUgJuLSfQUl9uhBrgVqpLyYOYnALxLs30ewAfCz98FcFLG+TsBnFKnDELHhNNsWxavyVQlmcwYslRJzFjgki6kR+K37bQqyaGBNSpqlF7ehn7eNdLbhUpEj8jYt+NhH2W4q44N3lQ9L0swDE+VYEra+BzlnnJKlUSqKsmtsjeFCAahEvHblGXZVdO2gdgrqcpCPUoqBJVAlWRawmExFeDmoPFZdYH2B9r2DjWXYAPRQLzP9614i44GnLRtwCSrp5qHX4UHqk4wIR1X4mLkMymznqE6HohgECqRqJJsMT7r4xhMjJ5eRhI9W5YxdYF0e7gYIDZSZj22qEy7RgSDUInoGdnn2xX5PJUrKdRtVzU+67K0DnVwMGEq8tnBADFP6VO25ATrGhEMQiVIeSO0IegrGrDTL/rj2F21/LVGHuWokoyLOCjScSW+g3EMqsutqJIEoQQTMwYLHvasOAYTN0kvI4meLWozF0i3h4sBYmocw1DbXgSDUAlCMmOwSpWUFgxGcQwZK7hZkjDQBWJVkj8Z+exSgJiayt2WZJFdI4JBqET0jOwb+1aoB4gozN+kj2OomkRPF8fAA/VlNyG9cJKLM4YJVdJAF2kSwSBUghRVgS2+6R5N2waMjM9ZC/VYErPhAumFk1xcAU11uR2L8VkQiolUSbbYGIBg0BmnbAOJjaH8dQJ31enttqjNXCCdKylRJfVWpMqoLrf+QNteBINQCXXGYIt6wPOyA9yqJtHTq5LsWK3OBbIC3GzpK2VQXW7FK0kQShAZEW2aMeiS3xkFuOWs+WxLXW0nM8DNofunutyOeZgxLCIYhEqoj4gtD8xIowIyd1fVJ9EbojrBBE8ZVAFHA9yUWbEEuAlCCdSHxJYHhjRupiZvqjoBA9iTMNAFIoeEqD04HFidcldNJdGz5QWoS0QwCJVQH3BbdK86N9Mmk+gNNSe/CaOUKsnFADE1dbhNtrQuqSUYiOgwIrqOiLaH/xdrjnknEW1V/l4konPDfV8mogeUfWvrlEdoH/UZt0W9MtLYBuIAt6pJ9DLWfHZtcOsLXa4kW/pJWdR1q32xMRhxEYAbmHkNgBvC7xMw803MvJaZ1wI4HcALAL6pHPK70X5m3lqzPELLTNoYeivGBKTJimqS1VMnYKJrDXFwMCG95rOLwYFqHYLI554L1AN1q7wOwGXh58sAnFtw/HkA/pmZX6j5u0JPTKiSLHngRxqvpLFhHIPOXTVQJ9Qq4mBIVEkI/7tnn1FjMYJkkY5VoAHqdvfXMPMeAAj/H1Fw/PkAvpba9lkiuoOIPk9Ei7JOJKINRDRPRPN79+6tV2rBGHtVSZPbzHIlZRifHVSH9EXa+OxicKCaf2uobV8oGIjoeiK6S/O3rsoPEdFyACcBuFbZfDGA1wP4KQCHAfhY1vnMvJGZ55h5bunSpVV+WmgQUpRJtqhXiKbXfDbJ0TPy9En0fIljKE068tnF1e9Ul9uhtv2CogOY+YysfUT0KBEtZ+Y94cD/WM6lfgnA15n5FeXae8KPLxHR3wD4aMlyCz2hPuO2PPC6dRTGhjOGTFXSAAcHE9ILJ7kYHKiuWz3Utq+rStoMYH34eT2Aa3KOvQApNVIoTECB4vpcAHfVLI/QMhOqJEseGH2AW/C/ahI9zYRBAtwqkA5wc9Fwr6b1CNq+5wL1QN0qXwLgTCLaDuDM8DuIaI6I/jo6iIhWAzgSwLdT5/8tEd0J4E4ASwB8pmZ5hJaZVCX1WBCFPFVSlYfaI90SodWN2EMmftuOUmI4aHxWVykcqqtyoSopD2Z+AsC7NNvnAXxA+f4ggBWa406v8/tC95CtqqSMALeqkc9Z3k1DVCeYkF44ycUAsXTabdfK3wQDnCQJdVDdVW15k/K0cQxmSfSAyUytJtcZMukANydVSUrks4vlbwIRDEIlbAxwyxUMFY3PwGTqbZMFf4ZMOsDNxQCxOLV8nCup3/L0gWNNJvSN+pDbMsXWxTFEAVaVVEkpjxpAdXutV8ahMEoZn1306olVSb6sxyAIpVCNz7ZkzPRoeoEdM+NzpEpKto0NZh5DRl3ICXAzOFBdt3rsszX9vEtEMAiVsNH4rFtgxyzyOfg/qUoSwVAF9W0bcDNATF232h/o6n0iGIRK2JorqakkesG5qiqp+nWGTFqV5Pv29JOyJDPHII5hiG0vgkGohPqI2PK86yKWTVRAah7++DoSx1AJ1XAb/XdVlRTFMdjSz7tEBINQCTtVSZODOWAWmOalBjVA3FWrQkTwyO0AN9U7TeIYBKEEqvHZlgdGu1CPURI9nSpJAtyqoraHi149STwLxMYgCGVQn3FbvDU8yk6iV6WMFKuSkmuZJOMbOqS0x5jt6SdlmTA+++6VvwlEMAiVmFAlWfLA6JbkrDNjmFAlRQFujr319slIaQ/fZ4wcu3XqutVDXaRpgFUW6qGqknoshoJ2oZ7Im6hiriT13OCzBLhVRW0Pt1VJEuAmCKWwU5U0HeCWqJLKXyc6dkKVJAFulSGaTKJnSz8pi2p89tm98jeBCAahEjbGMWhzJfk1VEmKYGARDJUZeRTfNxeNt7EqiYPZo2vlbwIRDEIlJpLoWdJ79En0gv8muZLUa40liV5lVGcAFwPEon7N4QpujhW/EWo92kT0i0S0jYh8IprLOe5sIrqXiHYQ0UXK9qOJ6BYi2k5EVxLRfnXKI7SPKgxsGSy16zGwiSpJJxjExlCVIOAw+Dx2MEBsQpXkuxeg1wR1u/tdAH4ewHeyDiCiEYBLAbwHwAkALiCiE8LdnwPweWZeA+ApAO+vWR6hZWyMY9AtyckcvOlV0Q/nGZ9tEYIuMPISFRw7aLxVk+i5qAprglqCgZnvYeZ7Cw47BcAOZt7JzC8DuALAunCd59MBXB0edxmCdZ8Fm1GeEVsGS92SnGODJRm1SfREMFRGTVEydjCJXuyEEK7gNsQZQ62lPUuyAsBDyvddAN4C4HAATzPzPmX71PKfgl2oD7ktD/yICA8/9WOc+Sffjrc9/txLlR/o6PjfuPw2LFoQvDO9uG8c/MYABwdTPCL8y7ZHsPWhp/HDJ1/AsoMP6LtIlYhmCJv+9QG8+IpvTT/vkkLBQETXA1im2fVxZr6mxG/o7irnbM8qxwYAGwBg1apVJX5WaIO1Kw/FeW9eCQLwhhUH910cAMB5b14ZD+ARa15zIE5YXq18c0ctxs+fvAIvvjJ5rVNWH461Rx5au5xD4TfecQxu3vkEgKAdfuHklT2XqBoLRh4+fPpxuH/vcyAi/OxJy/suUucQp5WzJhch+haAjzLzvGbfqQA+xczvDr9fHO66BMBeAMuYeV/6uDzm5uZ4fn7qpwRBEIQciOg2Zs50FIrowtfiVgBrQg+k/QCcD2AzBxLpJgDnhcetB1BmBiIIgiC0SF131Z8jol0ATgXwDSK6Ntz+WiLaAgChDeFCANcCuAfA3zHztvASHwPwESLagcDm8KU65REEQRDq04gqqWtElSQIglAdm1RJgiAIgkOIYBAEQRAmEMEgCIIgTCCCQRAEQZhABIMgCIIwgZNeSUS0F8APDU9fAuDxBovTJ1IXO5G62Mms1KVOPY5i5qVFBzkpGOpARPNl3LVcQOpiJ1IXO5mVunRRD1ElCYIgCBOIYBAEQRAmGKJg2Nh3ARpE6mInUhc7mZW6tF6PwdkYBEEQhHyGOGMQBEEQchiUYCCis4noXiLaQUQX9V2eKhDRg0R0JxFtJaL5cNthRHQdEW0P/y/uu5xZENEmInqMiO5StmnLTwF/FrbTHUR0cn8lnySjHp8ioofDttlKROco+y4O63EvERWuNdIlRHQkEd1ERPcQ0TYi+m/hdhfbJasuzrUNEe1PRN8notvDuvx+uP1oIrolbJcrw2UMQESLwu87wv2raxeCmQfxB2AE4H4AxwDYD8DtAE7ou1wVyv8ggCWpbX8E4KLw80UAPtd3OXPK/3YAJwO4q6j8AM4B8M8IVvl7K4Bb+i5/QT0+hWChqvSxJ4T9bBGAo8P+N+q7Dkr5lgM4Ofx8EID7wjK72C5ZdXGubcL7e2D4eSGAW8L7/XcAzg+3fwHAB8PP/xXAF8LP5wO4sm4ZhjRjOAXADmbeycwvA7gCwLqey1SXdQAuCz9fBuDcHsuSCzN/B8CTqc1Z5V8H4CsccDOAQ4nIivUVM+qRxToAVzDzS8z8AIAdCPqhFTDzHmb+9/DzjxCsl7ICbrZLVl2ysLZtwvv7XPh1YfjHAE4HcHW4Pd0uUXtdDeBdRPUWqh6SYFgB4CHl+y7kdxzbYADfJKLbwvWvAeA1zLwHCB4MAEf0VjozssrvYltdGKpXNikqPWfqEaof3oTg7dTpdknVBXCwbYhoRERbATwG4DoEM5qnOVj4DJgsb1yXcP8zCBY+M2ZIgkEnQV1yyTqNmU8G8B4AHyKit/ddoBZxra3+EsCxANYC2APgf4fbnagHER0I4O8B/BYzP5t3qGabVfXR1MXJtmHmMTOvBbASwUzmeN1h4f/G6zIkwbALwJHK95UAdvdUlsow8+7w/2MAvo6gszwaTeXD/4/1V0IjssrvVFsx86Phg+wD+CISlYT19SCihQgG0r9l5n8INzvZLrq6uNw2AMDMTwP4FgIbw6FEtCDcpZY3rku4/xCUV3dqGZJguBXAmtCyvx8CI83mnstUCiJ6NREdFH0GcBaAuxCUf3142HoA1/RTQmOyyr8ZwPtCL5i3AngmUm3YSErP/nMI2gYI6nF+6DVyNIA1AL7fdfmyCPXQXwJwDzP/ibLLuXbJqouLbUNES4no0PDzAQDOQGAzuQnAeeFh6XaJ2us8ADdyaIk2pm8LfJd/CLwq7kOgr/t43+WpUO5jEHhQ3A5gW1R2BHrEGwBsD/8f1ndZc+rwNQRT+VcQvOG8P6v8CKbGl4btdCeAub7LX1CPy8Ny3hE+pMuV4z8e1uNeAO/pu/ypuvw0ApXDHQC2hn/nONouWXVxrm0AvBHAD8Iy3wXgE+H2YxAIrx0ArgKwKNy+f/h9R7j/mLplkMhnQRAEYYIhqZIEQRCEEohgEARBECYQwSAIgiBMIIJBEARBmEAEgyAIgjCBCAZBEARhAhEMgiAIwgQiGARBEIQJ/j/my2/mD0VmzwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "most_probable = viterbi_solver(casino_params, casino_observes)\n",
    "with open('casino_viterbi.txt', 'w') as f:\n",
    "    f.write(\"Hidden\\tObserved\\tViterbi\\n\")\n",
    "    p = map(lambda x: 'p%d\\te%d\\tp%d'%x, map(lambda x:(x[0][0], x[0][1], x[1]),zip(casino_hmm, most_probable)))\n",
    "    p = list(p)\n",
    "    f.write('\\n'.join(p))\n",
    "same_indices = [i for i in range(len(casino_hmm)) if hidden_states[i]==most_probable[i]]\n",
    "print(\"Accuracy: %f\"%(len(same_indices)/len(casino_hmm)))\n",
    "# plot difference between actual and estimate\n",
    "# if upper zero, real is loaded dice, estimate is fair\n",
    "# if lower zero, real is fair, estimate is loaded\n",
    "# if same, there is correct estimate\n",
    "plt.plot(np.arange(len(casino_hmm)), hidden_states-most_probable)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### posterior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot of hidden values is reflected on Ox to better visualization.<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy: 0.780000\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x178c22dc080>]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXe8JFd1Lvqt6nTyOXMmR81IGoFGoABjIRBBGIwB+4IDDyMbc+H6GWM/OYcLF4N58HP2w77G2M9cX4ONAwhsgwzCIhgJIQlQAEkojUZpcp45uWPt+8fau2pXd6W9q7pnaujv95vpPt0VdldYtfa31voWCSEwxBBDDDHE+QXnbA9giCGGGGKI/DE07kMMMcQQ5yGGxn2IIYYY4jzE0LgPMcQQQ5yHGBr3IYYYYojzEEPjPsQQQwxxHmJo3IcYYoghzkMMjfsQQwwxxHmIoXEfYoghhjgPUT5bO16zZo3Yvn372dr9EEMMMUQhce+9954QQqxNWu6sGfft27fjnnvuOVu7H2KIIYYoJIjomTTLDWmZIYYYYojzEEPjPsQQQwxxHmJo3IcYYoghzkMMjfsQQwwxxHmIoXEfYoghhjgPkWjciehviegYEX034nsioj8nor1E9AARPS//YQ4xxBBDDGGCNJ77xwC8Oub71wDYKf+9HcBfZR/WEEMMMcQQWZCY5y6E+BoRbY9Z5PUA/l5wv75vENEMEW0UQhzOaYwB3P30KXxtz3H84vfvRLXc+2xaarTxd3c9jXqz04/dFx4lx8H1L9iK2/ecwP7Ty3jj7q3YNDNqvJ17nzmN2x471ocRBjE5UsHbrt2Ocik9g/jJu/fh4OkV632+/Nnr0HEFvrbnuPU2wrB2soY3X3MBiCjV8q4r8LE7n8aZ5abZjojwI1duwoVrJyxGydhzdAGfu/+Q9foAMD1WxdtetB2Ow7+30e7gY3c8jaVGGyPVEt76ou0Yq/av1Obup0/hdnkOX3TxGlxz4Wqr7RxfaOCfv7UP7Y6b29hecel6XLF1JrfthSGPI7sZwH7t7wPysx7jTkRvB3v32LZtm9XO7n3mND70n3vx89ddhGrIxOPOJ07ij/7jMbk/q12ct1DtckerDn7v5kcBAA4RfukVO4239Wdf3oPbHz/R12Osxnv1jtnUN8Jio43//i8PArA7/0IADx6cQ73l4q4nT+b2+9RvedVlG7B+aiTVOk+fXML7P/cwALPfIgQwv9LC+153mekwPfzv25/CJ+/Zb/371e992SVrcPG6SQDAAwfm8PtfeNRb5tKNU3j5s9ZZjzEJf3LLY/jmU6cAAHc8cRL/8vMvstrO5x84hA9+aQ+A/GzKuqmRQhj3sJ8b2nVbCPERAB8BgN27d1t15pZOAKL6equn6xd++SW4dOOUzS7OWyzUW3ju+76IZtv3QNquXYP0VsfF1dtnceM7XpjX8Hpw62PH8NaP3m00xk6Hl33PD+/Cz7x4h/E+f+TDd6DtCnRcgWsunMUn3p7P7/vk3fvw3//lQaPfopb9i5+8Cj98+abU6z3vA19C283mZbZcF5tnRnHHO7/fav2bHzyMX/jH+wK/tyXvzXe/9lL87s2PeOeqX2i7AtdevBolx8HcSivTdgDggfe9ClMjlbyG13fkkS1zAMBW7e8tALLN52JA8lniRlh39enQa++FogMC9iXqKZkAIRD+WM8RPn2RfoxCLms7NLVLVwjvWssDalvC4HirRU3HQbA+rdrOs91D3pkLXGvyOwr82TcIeQ6pZyCm2+HXopmUPDz3mwDcQESfAPACAHP94tuB5AvD9oY453DyCeC+vwcOfRtYOQ2068DEemDz84GrfhpYc7HxJtUR0R+Mtpf8AGx7uIFIgOgyIDb7FEL+vjx/YMKMMwzeg8pwHETZDWfW308hv1e9dcj8QWcD9RuyHg//PBTLpiQadyL6ZwDXAVhDRAcA/A6ACgAIIf5/ADcDeC2AvQCWAbytX4OV4wHvO/x7dSKcqPPgusATXwFOPQlsfwmwflcfRhmB5jLwyL8Dc/uBC64FLoiY8n/9z4CvvJ+vyg2XA1ObgHINmD8E3PVh4M4/B17xO8C1v2x0B6qbytVcd+v7S8ALlPULnhEwWKfbgNjsU0BACGG9jajtmkKdG9PDTESZPfesv9+7T9F7rZUc8/NqAyF4HE7G42F7Hs420mTLXJ/wvQDw/+Q2ogT43lz42XLjPLfmMnDjW4C9X/I/+/73AC/9jVzHGIqlk8DHXw8cedD/7MW/CrzyfcHlvvW/gC//DrDrR4DX/jEw0RVwWjwGfOG3eJnRVcDz/2vqIahj0tGOXRS9lQRXCJT77Ml4FIkBT61+j7XnTvz8dzPSEj3bla8mx9tf1oaWyWY6XZFtZhY261K/RzkFfffchQDJsdhe54BmUwrGBhSuQjVsuqdDxN0QX/1dNuyv+SPgVx4EnvMG4D8/ADx8U1/GGsBn3gGceBz4iX8A3rmfqZWv/ynw6M3+MovH2GO/8OXAGz7aa9gB/uzH/xbY8VLgP97F3rwhdFuZiZbpN+eu7SstsvKjBOm5W64fuV0bWsaSYiLKzrnz+c3Bcw+lZdDzXT8QoGWyeO6W9NjZRuGMe9qpes8U6sTjwDf+Etj934AX/Bwwsw340b8G1lwC3PoH/b3S9n0TePyLwHXvBC79L8DIFPBDHwTWXcZeuCtz8u/6MNBaBl77J4ATc2ocB/gv/xNoLQH3fTz1MPKkZfKmLcKQRMGFISs/6hmCPtEyNofbdByKWsoCIUQmY+YZ8AAtw+9LGY6FCYTgY0FE2Tj3jHGcs4XCGXc9myEM/rS860w8cCO/Xvcu/7NSGXjxrwHHHgL2fjnvofq460PA2Brg6rf7n5WrTAfN7QeeuZOvoIc/A1x4Xbpg6eyFwI6XAd/5B+YRUsCjZdzeG84UlhmURvC9XYsMkwy0jBB9oGUSrtsweNey6b6Q/fyIrLSM93uD2wR8WiYLVZIGrkbLZKGAhHceimXdi2fc5Ws0LRNczvvwu/8CbH9xL9Xx3DcA1Qlgz3/kPFKJTgt44lbg0h8GquPB7y75QaAyzmM7fD9w+mlg1+vTb/uqnwbO7AMO3J1qcZ/39T/LRstoR/ngvcC/vh34ygcst9iLbLSMpefu0TKiL7fyYGiZHAKqEPa0jJZGqhtVP9nBfEZmOYx8aJmk89BcAr72J8AX3wM0Fux3lDPOWps9W4RF4XX4kW3tTBx7GDj1BPCiX+xdoVQBLngR8NTX8h4q48DdQHMBuCikGKQ6Djzr1ZxBM7EeIAd49g+n3/aF18l9fAvY9oLExT1aRr/h7HkZn/pqN4Ab/yuwdJxTNtddyg/NjPADbwbDSsqWStwnIDr+lD4v+NsymIX0rJsOnPqXlZaxOIannwH+7R3AkQcx8f3/yNvp2iYAKCWJvtMy0LJlMuwt8Tzc+vvAnR/i950m8Jo/tN5Xniie554QjAnNljh4H7/ueGn4SjteCpzYA8znkJ7/bz8P/P424JZ38997v8JGe8fLwpe/4Fpg+QTw6OeB1TuBsdn0+5pYy7GDA+l60YbRMrbT90A2xb1/x/TST/wjsHk3ezA5uGV2GSZyXes8d4IrRB9pGe3DhaPAv/4cz3hWTves413Los1xoVv/EFg5k2pfWQ+/VRHXrX8AHP4OUBnBrrt+DVW0As6D2+V4DSxbhrLRVLH02Jl9wDf/GrjyzcD3/d/Atz7C8b1zAMUz7oj35kI/PvIg0x+rdoSvpIz+07dnG9zBe4H7/wkYmQbu+gvg2CPsuW+4HBiN0JHYcDm/Hn0Q2PAc831u3s37TQEK9dzNdwl0Tdsf/BSw8Urg4lcAz38rsHCIZ0sZYVPJmJUfVQUvIucS3NDr9vb/j4/dAzcCd/5Fzzpq2Q1Pf5a9w1t/j7OpUuwrq+EUpg+3xePAdz8NXPlTwGv/BONze/EC55FgtoxKhRw0LZPxeMTSMo9+nr31l/4Gx++EC+y5xXpfeaJwxj0sCh9AV9AGABv3Dc+JzkBZdxnglNkYZ8HtH+Tc87f+O1AZA+74n/wUX/vs6HXW74JnRNZbGPctu9lrXjiaanH2Ynp5UFN40/ZWnb21C1/GG1f00xP/abXd4FjNPbzsAVXSqD27bYSh57pdOsEVyFdcD+x6HXt8PV65QAkdXPDQX7ETcNWbge/8I6+bsK+sdtM4FfLhz7CRu/pngZ2vgutU8RLnwSAtI1/9IqY+e+4QMlsm2/FQ64Yej2fu4Nnz7A5gejOw+uLsTmJOKJxxD53eauiZQgkBHP0usOG50RstlYGZC7hq1Rauyyd11+uBVds55XHPf7AXG5f9Uh3nCwKwM+4br+TXow/GLydB6M6WMd8loI4/sTxCpwlslZz/9GZ+mCnj3mkDy6es9pEUPA+Db9xtA6r8MDH2XJO2q65bldj02BeA9gpwzTuAa34BaMwDT90WWMcVwHPoKYwuPsPxohfewDGN7/5rwr4oh2wZw4DywfuA8XWcWlwdw/za5+ElzoNdabdBzz2jtlkiFLWWNcAcmRYqBPDMXcA2TW1yx0uBp+/g6/4so3jGPUGAyX/KyjdnnuEbJ864A5xaeOoJ+4GdehKozzFNArAGjOJRVydI6qqx2dAyM1I6+cz++OUkum982+mqd8Hv/wZ/sFUL6F5wLccBhABuugH44C42ZoZICp6HjstGOKxVB/7pTcAt7/a8vLyFw9SIvN9y6NtAbYpnjZuuApyKHxuSEAK4wpHX5AUv4kD1xHpeN2FPA6dlDt8PbLrSW2lu44txqbMP5RVfE797RtT3gKrwhcOy0jKhh+LkXo6X6TIi21/CCRRH7rfeX14onnFPCKj2ZMuc2Muva54Vv+HVFwGnnrJ3ZQ/KoObm5/PrJq3b4JoE4/7cN7DHP7nRfL+TGwEqAXMHer9rN4Ev/jbTQ+0GAL6xAt6U+R49OAQ24rMXAuNr/C/W7OQH6p5bgPv/GXBKwKd/ho2o6fZh57nH1YD14Au/Bez5AnDXX+DKlW96wmFG20hAz2859G1g4xW8k3KN6blD3cZd4ArnSTRH1gBTm/nDjVewIY1BXsJhqbN0WivA8Ud5bBKLa/j9yGk/uNhDy/Q7oAo+FllpKkXv9EA9ZLdc7X+mYmgnHme5kzs/xNIjZwEFNO7xwZgeWmZBludPJehhz14INBdZAsAGB+/loO1a+RDZ8Bzm8UG87Tg8+4eAN/69HQ9QKvONP7cf2P8t4JNv5nx5gA3WnR8CvvReDvCCZz55acsQiPfVPTNRv/dL7wFGZvxq2gPfMtqHL+9sNi593UR0WhzQvOJ6YPYivHruRrhCeF5fXghct+0mcPQh9nQVNj0POHR/gKtwBXA5PYnF1Zf718bGK9iQtqI7TTENkc1wuiYVqkcf4vzRjf7vaU5dAACoLTwT2CZgl+JqA559kJyt2u8sMnPqjPxts1qixsxWQN0Xt/wPdq5uuqH/PzYExTPu8jUyz717wYUj/Dq5IX7Dsxfxqy3vfug7fLM6Jf67MsrT6Jmt/L6fmN4CPP114G9lzvzX/4w/v/8TwMQGNgiPfp4/60oLy5Dmzhf83H55QWtQxv3EHmDr1cDOV/Hs4snberYTB6sK1a51E3H4Aea+L/lB4OJXYHtzD8jtIOdkmeB1e/wRoNNgOkZh8/OAxlzg+qPWAi6iQ1hac7m/3MYr2JAejc5GYhoi23iNKlSPPCDH5o+zObYRTVHCyOK+wDYBXYphQKmQyHY8+FjIo/HEfwIffgHf72f2cZxBv7/LNXa2HrsZuPejwLpd/D6HBANTFM64q6ly5MnqCtpg/hAwtpoPehzU09eWdz+zL/gEB4Dr/gfw8t+2254JZrYC8wf5pr/4B9ion3qS9WwufyNw6et4ZjF/ODdaRgAYc5c5zjDdZdxntnFuP8De3MgUG6+nzIy7zRiNb+J9d/Hr1muAzc/HiKhji7vfjJZIgcB1e0y2mluvxYHWyZZ4J/Z4H43OPQmHBFZWadlWivqI4XQHLhw2d4Af3lNbvI+cchkHxFqMLGjGXX1nQbfZQNEymQOqkLMYIYAvv49nTh//Uc7CUzEvHau2+9TZ//UxPjbP3Gk/AEsUzrgndWLy5TklFo4AkylalE3LC3PBopCp0wIWj/q8qMKzXwtc8RPm2zOFMq6zFwEvfxd7ol//U8Btc2ris17L3+/9MtMyuWjLCMx2ZPplt+dervnHU1EPF7yIA4ad9O3O7PKhux7uSdh3F9+MUxu9YPglrT3m2SIJCFy38wf5w2ntelG04YKv8llZ4uPbGtdiMVNbomMsEgMXDps/zDNjLUhBAJ4R6wO0TLdwWL/1iVSVMT/sMgZUCeycHL4feO4bgZVTzLmHGnempDC2mrOH1u1KDIL3A8Uz7gmRdq+IRS24cCiZkgHYIFUn7dL2Fo7wiJJ4/X5BGdftL2ZvkEp+dsqanRwHcCrAyb29ee4ZXPfVLWncuz13wKdmFA87exHPLJREcXM5cRc2tIxxherhB/wMp9kLseRM4Fmdx3KvUIV+3S4cAWrTQa2hiXV83rQq6eoKU4pN3bg7Dl/Pim6MwECFwxYO9yQDELFxH13c511kfrB7MLRMQDgsw3a8+MvTd/CM9DV/yOcKiPbcAU6uIGIH59B9A+fdC2jc47257qkfFo6wV5YGY7PAskVkWxmsbs99UFAX2PYXs9rk7IWs81Ie5VmLU2Jv4vRTcLqCS/YVqsBsWwafw4z7xiukRywfeOoBdPop4FNvAz54aeKD1K5CVa6bxjS5HfailaflODhU3YFtnf3RGRKWCMxCwhwOp8RpjtrMsbp0BC1Rgju2Orjs5IbYGWZewmGpf//C4Z7fQ0TYJ9aj3Fr0zvPZEA4DIZdOTA6BYyWzF7KdUGnLccZdZcxtuorTos8807tsH1E84y5fEzsxgSRdcix9iuHY6sTqv1CoafbZ8ty3vxT44T/j7k2An7Gz+mJ/qrxqB3DqqZ4ipizZMrPto0CpykapGy9/N/D223wLPS1vgi++B3joX4H6GVbDjIGNETASDls8yrMJ7aF8prwGs+4pM881BQLX7fzhcIdjahNfS6eeAk7sRXX5KI5ile8lKkxujPXc+bdnzJZxU8xcmst8vywc7rn2CcABIdNj55h3V+fREw4bgCeraJns2TIEHH/Mrzbf8n38OnNB7wrrZfxk+4v5VVGTCSmseaN4xj0lLQOCTGsUZsY9k+d+lox7qQzsfht77YBv3PXK2NkdnJ5FIh/JXwHMto6wYQxLCC/Xgno6ioM/8gB7Pxuey6X0MVC2xUg4TGYSpnI65xT37QcC58prsUacgsjaZ64LlcYp/Eb5k9jyzfdzhlFYHGhqI/DkrcCfXwn8xfOx4el/w1Gxqve3JHruOdAySEgF/e6/Ar+/BfjjizmoHuK5nxKT/If03HuEw7INMREeLZNDnnsVLW5ar4z7jpcyRbP2kt4V1l8G/OpDwI6X8N9K08qia1oWFM+4JwiHKTgE/wZIa9zH19hx7vOHWEtmJEIcbNBYo3nuCqt2AI15zGIxN+Gwqc7p9Me2MuJ7+JuuAp7z4xxkipkpWdEyfjJk8sLzMiipee5zlTWooYkJsZArLbPxiU/hhvJnseGRj8rge8hxCzH4R8Sq3nFMbuBp/kdeDnzip/yHlMRAhMMeuYlnPep4d43dIeA0pHGXldpnTTgsB1pmBx3m37vuUv7w0tcBv/xAOC0DBBwGjK7iGa5NskYGFM64q+l2qk5MygvXqyfjYO25H2SvPdcIXAYoPlBdiIDHA26lo7noubsuMNZZiFa7DIPi5jdeCayVY4upK8giHJaKlgmZcc2V1wIANohjGO3k13hh6vg9OCQ0Oeewh6Iax/RWT6n0iFjd+5hShvTQfcCjn+MKZA25CIfFGXchgH3fCDaW6fbcQTgjJvgPj3NnlAbUickTDkPWbBmBiyEdAeW5E/VmiUWBiOtNUor75YXCGfe08gMEAPV5/qM2lW7jY7NcSRlT/ReK+UNnj5IJw/rLgP92i8/BA14O/lYcDaZCZtjNaGfBbLaiboZNV2p1BU9FLu7z1Ol34Z3/NA/auYM84xpd5X00X+bg5d+03oUPPPpDwGd+If3Oo+C6mDp+L27rXAHXqfBnYcZ9nB8smL3Qy2dfQTWcllFY/xxf+kIhD+GwOFrm9NPsher9Ebo5dwLOQBr3FWncB+y5+8JhWWkZYCPJGaYKlpoigUrrBwpn3HsEmLrgTcoJrG8CcBFNGqisBFPvffkk90g9l7DtGr9aFvCmietxKqDGl0U4bNRdYO36tFi1HQCx/sbMBfz+dIxx77dw2PwBpmQ06zlXYQNbgVT1e+Tfs8sXHn8U5eY87nGfheVpKdXQnQED+DPMna/ioioAZ8REr5HVHwwXv4LTOTXdnqyeKpDgue+TYnHbXghc+yv8PiRTrIMSWhU/vbhb96ffnDvvj3LRc58Vc5x9Vpuw28jkeqbjBojCGfek6rZAiXPD1HOXN5epcW/Mmxm5s4HqBOCUMdXNuVtuzhEdjLgrZrTMNb8AvPnTvE5lhL29GM+978Jhcwd7PM6Fiv+QfnjyWj63WtWoFaRnfa/YiT0v+XPma3VdGYWdrwLe8lk+Ts/+Idz34o/gY50f7DWyiq/f9kLO2nBbXC0pkYcOvUBMIdjxR5hDXnsp8Mr3Ae/c12P01Lqt6ozvucvvShZ0mx24FWQewmGzNJ+e3g3D5EbOtvnLFwE3/WJqFdcsKJxxTyscBoBpGSql13ax9dzrc+lnB2cLRMDoKkxjqUs4zG5zY2KR35jQMhPrgItf6f+9ake8595v4bCFIz3G3XWq3vvb117Pb1I2II/E3AEIEPaLdVie3AH8xMfDr0ki7ovrOAARjm14Gdoo9xr30VXAT/8b8JM3+gVY2hizCmUBCcJhC0f8ilSiUMdGrduszmjZMpKWGZBwmK7n7mbgqVwBrEJG4z6xHoAAjj3ETVoe73+3puIZd/maNFX3aJmRqfSBTs+4G2TMtOrcrCLt7OBsYnQWU1gIbaBgikll3E08927Mbo/n3C0qVL0l05zy+lzvw0lbb9/E5fx9VuO+cBit0TXooGRYlRnzoLro+/nantrIldVagQzTMtmGHLt+SEVqNwLGfaWLlhmocFj2tEshgFViPhv1qh+vX7wPuOotGUaUDuW+7yFnJAmHBYI2jQUzo6uezCaFTB6vf47TMoD03BdzyXOfEEv8Jkv656odwNIxoLkULMWXyFKhmpjG6Ha4qULXeXOI8K7yb6HuAuOOw2mbSvXQFgtH0BpbD5zuE8XUxedmrcgEEmiZhSPBTKwQBGgZKR7WS8tkG2MS+DfIsWTal8AqzGWkZWQacKnKvSMGgAJ67vFpVAHhsPq8GV1Sk3m5DYMUuPocvxbIuOchHJaL567y3iNoMLtUyJQB1YhgOwG4tXQNbsX38bU2vTl7Cpsy7jAsyPKu5YRfM7E+OMaMFZkA4oXDFo4ke+7ytVGd6c1zl1an38JhriuknnvGCtWOyE7LKFZAb+LTZxTOc0eCNxdIhWvMm3nupQoLbLWW0q9jmm55NjG6ClNiIZciJt9zz/BQUymIK6dDi0GsUiHVuknWPeK8KblcTzhsYj3PLlzXrDXTwhHgc7/KHvXcfrS3XhYYXxp4mT9Jv2VifaC0nQz3E7pvEfGAbCzyfZUgxufTMtO8fKc1eFpGjQXZjkfFXcYImtlomfXPBV70S8AL3pFhJGYonOeelCMb0BapGxp3gOmBFIqFHhrF89zzMO6TsAiohowHgN9rtgs2ZeqpaZkIOk0JqwkhhbMm1rN08oph5fI3/oqbNBy8F1g5jfa4nKXY0DJpjHsXLZPVbkYKh6n9JHLuvG6zIq+PldO+5G9YQLXTAm7+rdgmJMYQfCwcJxtNNd4+I99kMO6lMvCqDwRlnvuMVMadiF5NRI8R0V4iemfI99uI6KtE9G0ieoCIXpv/UOW+5GuicJjy3E2zWKrjzAFH4dHP+80WAN8DPNezZQBgdBXGUAd1mt5HttPVSchjlIWWUeuunAn9WtkWMyojJS3j0Wkhnjs02zixjl9Nc5SPPhQojrKjZXTfMwYT67hFZIMfuFlpCCBGOMyT9Ejw3OVrQxn35VO92jL6GPfdBXzrr4G/eqGdBEgIVMYPIdvxGGsp4742l3ENConGnYhKAD4M4DUAdgG4noh2dS322wBuFEJcBeBNAP4y74H64+HXyFOln8TGvM+jp0VlLJqWObEX+ORPA198t/9ZoTh3vtHGXf/32V7yk2IRLapmayGY4LlnqVBNtO6RdBp7eV4Rj4oLmBr3Yw9zV6wJNoLKc7exMYm0jDK0cow5OO7RpWBe28q0nrvMf2/Mx3dieup2//0DNxqONhwC8hdkPB7jHXl9nmuFiglI47lfDWCvEOJJIUQTwCcAvL5rGQFA3SXTAPomf5ZMy8iLRwhLWmYsmpa57Q9YPOipr3leknGh1NmENKYTYt77yNahmcIi6iXDB2fEeKKMe1J8JQx+hWqCRYx4KDvEHmWAlgHMGqcvn2K9ofWXccMGAB1l3NNvJT3F5M0ueIz5CIeJcDoopeeu1m2XxvhNa7lXfkBf4enbOdhYmwZOPm4/cA3qAZ2Vpppoy+tzPKSq+BxGGuO+GYBeTnVAfqbjfQDeTEQHANwM4BdzGV0I1PUWJxxGREBrmQ2xKV1SiaBlhAAe+Rz3uuw0gSe/yp/X53lUVcuy5EFijIWrJlw/G8g+W2YJ9XJG414ZBUo11nYPgWfULCKqiTx1BOeu0zIE2HnuxyRvvP45wJbdAMhrldcXimlCee7sVefiuYuIGcPSCU46SHBm1MO15YzwB5rD5HQLhzWXgQP3sETu6otYWjcH6MJhWWiZ0ba8VjSarQhIY9zDTnH3kboewMeEEFsAvBbAx4moZ9tE9HYiuoeI7jl+/Lj5aLXRxMkPEOCnM9oEVMNomaXj3Jv0qp9i72LvV/hzVZ1qkklxtqA8d3fR+8ialsEyGqUcHmijqxJpGbMKVblukrcblS0jvV7PuNUm+IFv4rmroOD6XcALfg746X+DK6f0diJoCQt2zS6on8JhjYVUhYHkee7SuLfibbI0AAAgAElEQVSWPQPbE1A99SRLKGy8Mlfj7sondNaHXVVI3Z4iOHAa0likAwB0bcst6KVdfgbAjQAghLgLwAiAHoJKCPERIcRuIcTutWvtghNJaVRe8YVtimIULXNaVgDOXsSKhqpBcWOejX0RII37pKvTMnaX/QgaaDoZ+HZ9TFHG3SbPPW36YP0MC0GVq4GPHWKj4AotW2RiXWLP0gDmD3KxyuRGdhYuerlGrZj8FjWmhB8ztpobRyjOHchcISREhL/SXExl5HzjLq+R1nJQsRXakdAzcGYv4mYmmhCaNQQ/rLmoy/54VDsraKAaFOIrANIY97sB7CSiHURUBQdMb+paZh+AVwAAEV0KNu6Wrnk8koJsrnLdTRUhFSrjTOl0Q5V3z2xjPRLFPZoWSp1NyLTFCeF77rYeXg1NtEu17GManYnMlvECbwabS50+GJFJRdIQCMC/2CYMFf2WTwKjs4EnjJ/5k34zqakEx2EJAjlbzUPP3RUxnnsKh0k9mD1aprUSeFgRaQ9tdWwn1skGM4JlhTNCSOEwpmXst1NxV7CCkczjGTQSjbsQog3gBgC3AHgEnBXzEBG9n4heJxf7dQA/S0T3A/hnAG8VfZJ8SxIOg0fLKM/dkBeujrF30o0zXEKNmW2+NnOnxXRNETJlAK/Ev4aG95HtSRpFA20nD+O+KjoVMoNwWCJTXY9W8hTyP8+4mXboWj7VI+mbtoNY70BSSiPVJrxYUR7CYYGHm47GQqp7Sq3a9jj3Ja2RTpf+jWfc1wOrL+T3J/faDVyDV4iWsYlOzV1BnYpn3FNVqAohbgYHSvXP3qu9fxjAtfkOLRy+N5dAy6hpXWXMbAdRRUxnnuEbtjbB08flk8DfvBI4/B1g5w+a7eNsoVRBGyVUXX/Ka/sMrqGJBSeHC350FeuRh8GLr5hTGcm0zFyoB6p0WVw9W6Q6wTo0abF80gte+9tV4zOnmFK1+6tOeJ57wHDaQkTst7HgZ+fEQK3LDgCx507B771jsXiMx1+b8PuNKmcqA5RwmKNdR6mauHSh6taxUkDjXoAoYBBJ01vWkwAHPwGgbHhSKuNApwF02sHPz+zzS+RVju/h7/Dr1T9rto+ziAZqqMEvYrI1AiNo+l5ZFozMxGTLWGwvbfpghAa/ogsEtAdEUmFbN5ZP9nruVrRMcN1Y1Ca8GWfWnqG8bxEhP5DSc1e/F5C1I1q2jKf3Ij9YOOIHhdU5ibgmTKDSom1mgDqK6rkXzrjDm95Ge+4EAG1JPZQNqYOqysvtuplPPyO7ByFYwPGGvwV2/oDZPs4i6lTDiNBpGduAal6c+yo2Su1mz1fKy+pbhWoY5w6ZCqlzzrUJv64hDVZ6aZmk6zYMfgAyreeeX4Vq4OGmI21AVW1HgFNeW8ue1DSRykqSCy0e8427U+IEhQiqzgRC/ggb6WgdVVFH/Xzk3M81JAXZhJpOtiX1YOy5S+PeTc0sHPFbiemd69ddZrb9s4wG1YKcu831LoT03HMKqAJ+UZGGpOB5xNB43cSA6mKoB+o4zFe7AkFaptPgGEsS3A5n/3QZd5tZSEAnKQlVjXM331XvvuNomVSeu1aoVB0LBFQVDe7TMkd8SVyAr4mMnruuDmoTmNfBnnsOmWEDRuGMe1J6nJctozz3iqFxV16JnjHTabEnr6aMynN3KgPTZs4LDdQwIo17ybEMvHWacEigkwctU/ErGLthJRzWtW4k2iucCtmFAF+ttqGuibBAezfqc4BwQ2gZm1mINqgk1Py4QB567qG0jNvh82RAywgh+Bw3lwIVt0p9E0DQcwdiM6jSIrgv82Ovg417Do7MgFE84y5f484T0zKWnruiZXSOVRVEKeM+uoorK9c+i2WCC4QmVVETTIGUbI2ANMS50DLq4dvuzWvOIhyWiHYjnLIjnQ6RUP1Bk6iZdsOX3u3JlmGYZcso7/Ms0DIev6nBKwxMny3DtAx77sFsGZl73lzm+IcepB2ZiZakSAl9Xwq2h6Qq6kPPfRBI1JYRgsubVbaMqQGqhBh3NUVUHC0RsGanpxtSJDQ1Wsa6ck8eW69AJQtiPHeF3GkZIfhhEiJ65hB5PWY97191iUoKqn71d4GP/wi/78mWMU+F9GchKRYOBFSz57nzfrt2rIx7Cs49oPwoA6r+hIikhg9YKx/o9dyz0jJqX0gxi0vAiLuCRgEDqoVr1pHkzblC89xLVXNZAHUj6wHVeogOyVs+az4rOAfQoBGMgnO2S45d5Z5orYAAdPLg3NUxDKlItLspU6QPxgTbueCly+urSk81jpZx3aCaYWS2jMEsxAtApvHcJ1nzqN30CrGyIJSWUb/fJFtGgGfDS8cB4Xd38iQSPAE3TTp6JEdaxiGrY69vqCrqqOdBQQ4YhTPuCtE1TDKXtd2wM76el6Z5kuoC1POiswj3n0U0qIZJZKNlRHOZjXsutEy05+7dlAY5bKnSB7002RDOPYyW8a6JGOO+7y6/ahkIyZZh2MQPUj3itDHmkeceKhxmoNcUaExdGQWay77jBZWVJPzZkN5DV3nukeplydANuRUlptBagQMxpGUGgSShQKF77jbGPYyWKVIT7AQ0qYbRjLSMkPx4Ow9tmYpfnt4N76Y02Fyq9MFYz91fz7MraTj3x7/IAXaFCM/dphNT6jx3gI07UQ567iHyA55xTyGgFQiojstsGa2ISD1E1X2mzwZGZngWEnJNmIJIsxk2G5BOx5CWGQCSBJg4tTUHzz1Ay4R37SkiGlTDKEnPPQMtAwBunp57SEDVjqdOkT6ojEYo5+6/9wxRmmyZ+UOcIvuGjwGP3OQH5r3tKk/WpEK1axxxqPoPIKVJnwWhwmEGAdXA8a+MAq0lmV6pvqfgNgOeu6bzXzWsMJfQs2VCOz+lhTzn9TwcmQGjcMY9qdLP69rerpsXMAG+sdG9tDDOvaBoUo2b/UKlQppvQ0jKKhdapuxLwnbDLlsmuG4o4grcQsS+Uhn3xSOsq77l+fwvYrMmx1t0c/9xUAa3uQRCObPkL6/ftWP1+1OpQmrph1WVLePPjLyMnihaBmBqxrLnaFgxm9UxaRbXcy8eLZMgwJSZlqlNAtPbgKdu8z8rUrelBDQ8WkbYT98lZ51vtkxYKqRFnnsa4bA4zj3wXv6VhpZZOBqruWIjHNbD/cfBewAtyPOanZiJ5txtUiGXIUTH+8KLC3jGXXtgqOBqhqCqXqpAWXgZOb6hcR8AEoXDVGWdredOBFx5PfDEV33N9vocZyMUTM85DC1nBGVyUUHHb0doCCENcT4B1WjP3d+h+RhtPXc9y8ajFipjACg+FXLxSGzruf4Lh0nPt5FfQLWH2jKiZboCqgDKnaZn9D3hsGZIeqXuuVvCr1Alq2PvQdKzufQuGDAKZ9wThcM8WqZh37z5ijcBEMDDn+W/i6TZnoCmrLQbRQMlsuzYo2iZPNLDlPccwrkD5trkOtcaiRjOXV/Ne08ky/sjPPfWCjsAEzF9RS1oGWPhMMAPqOZSoRpCy5SqqQr3ApRahR88FXfFOy+ecFhziQPRetOUvD13qLFYbEg+0IfZMgNBfHDEK6yz9dwBYPZCFi9SDQPqZ84LSgbwjfsImnAcy+m7pDU6eeT5Ow4XmkV47qba5KmEwxLy3P332l/V8Wjj7nUSWh/+vb6tvgmHaZx7HgFVhDxU2o1QKit2OwLeQ7TUXtG2ST4t05194ylD9uoNpd6vK/eiyQ/YBVT5umwWMM+9cMY9qVqPU2MzZMsoTK73W6tFyMMWEZ7nTg17DRJJy+SSLQPIbIrwtDdTiiFV+qAnTRGSLaNdYIFtxClDLqhmE2lomfRI3TJQjQ8AGguesmUWhAqHtRs9bQmjEFhXZryU3bpGy4BH2QhRmUxRtZwE79jB7th7KHC2TOGMe5IIkJct01qx99yBYGu184mWkR7ICJr22TIqFTIvb0Zqj4TBMQz6qmXjK1SVcY+/PgIpiLGeu3QC4jx3dd0aHHCjPPfyCEAlYOFwbsJhPeg0U8t5BArQpLGudLpoGRfhEsJl2eAjgqpLA5+e8wdjVaG6wtXcjSEt038kVZt500mLKWQAkxvYc//8r3NTjvPGc2eDPIqmfT50ewV1UQGZSjtEoTISXbBCpqmQKZZVRiOJc9e/qE5GB1RTeO52BVl+UDARRMCzXgN8639h18LXMwuHIaw41IDqDPxeuY7jNrVsGRVQXQqmQQK843LMNZECvoSEdvRMD4nrAt/5JzxR3ol6aTx5+XMMhTPuqYTDsmTLKEys52yZu/+G/86hWu5cQKsUpGXsNlJHHVWrlmWhqIxFemlq9p4aabzdVrRiaGi2DMAGqBHRam/xKEBOrCSFVUGW7n2mwY//DTC2GpfNfy2HRMgoWibdPRX4vSWmckqiFaBlmHNf7DXuAD94s3ju8pXITjoaALD3S8DJvfj38R/LLD52NlA4455eOCwj5z6xHhAd/+/LftR+W+cQWtRNy9ikh62wcc9rUOWRSH7VnJZJkT4YIwcdCKjq24gzNvUznOERlyqbcN2Gwc+WSXmkK6PA7IVY1TrWHz13A+MeuE+lcXfclhfT8ITDmkvhqZWV0dDah7QQ2rGzFg7b/02ASriz9uKsPbbPCgpXoaoQdZqYlsnBc9dzln/lQb9/asHRdPxUSIcIHQvSndorWBHV/C74mBuZudm8K1RjjHtYKiTA11K70bM8AEktxFdt2hwrPSiYGtNbMH3sW7nID/SMudMw4Nw1b1mmTpLbDvwWpmUiPPfyiF9sZgG9E1MSlRsJ2du1jVKerszAUDjP3UkIfQshUEIHcFvZPXeAc3Sntthv5xyDCqiOognHsRVTYlomt6mq7LEZBtPMj3TCYdFy0PpvCnjMpSoHFMMQZaBCtpt75k83pjdjunUsh1RIEULLNFNnywBy3LrnLlreMXUc8IltxNAyWTx3bQzWtMwCF6aFpoUWAIXz3P2ChOgK1Rpkr0vTFns6lOe+9hJzTfhzGG0tFbJkmEPub2QFDeTtucdky5gYxFTCYemkKQKbKFUTPPd445503YZBaEHB1JjeiopoYgb2OeKAnAGFBVQNssZIbUdx7m7LOy8E8rVlwmY9mT13fnW0KqbUM8DbPwg89gV+aK/aDveU4Tk4R1A8456QsyogfOOeh+e+9tn22zgHEaBlTMs/Fdr1nDn30egb2Thbxl8vEjG6QxTw3PUx1qIbZKcx7lZ57hYe4zTPMteLE4Yr9u68Z/bTaXqGOg08jRuHzYzjtqFODHv1HT7vYcY95oGfBm4ILZMKj34e+Mr/y++dMrDtGuBkSPyhACicS5qcLQPUSBn3DJz7yDRw8Q8Az/4h+22cg3AdvjkraBsHKxWo00RTlHPMlokvYjJCmvTBGOOue/xODy1jz7nb0jLG1Jc07htw3Gy97n1D9M5+DMX4vIwYL1um6T2sHCJUO/Kch+nDZ0yFDKVl0hz7vV/237ttYGKDzByyHspZQ/E8d/kanS0jUBM5eO5EwJs/bb/+OQpXnvIKOva0jNtGG6WBBFQdQ815v4gpZqF2PZKyo4j37Lk3wyONKTh3BdOcfeNDPL0VALAB2Tx3NzTPvWnkMDH1Ao1z72i0DFDpyDhLKOc+krGISctzN8mWmTsYLGCcXC/1qopn3QvnuXsSHRFfs+cuA19ZPPfzFY6DtnBQoTYHVG1c906LjXtek9WEgKqR2FaavqMxKX2RtIyiI8KCqga0jAmsaJnRVWg6I9iQkZYRYcJhphloJGMgMlvGcZv+NgmoCumZh3Lu2WgZXS7ZiBKbPwRsuoq1pQBgcqMvI14wFM64e9PUGOEwn3MvXslwv0FEaKOEMjq+7KrpNtwWmijnN1WtjDL/GsJpm2qTp6HcWZoimZYJPCCUUQsLqqYw7jbdgDydJBMQYaU0hUnEyBOn2TdCZj8G8gNQ6wdomXaQlnGVExZyLrIWMWl57kbHfv4gMLUZWLOT/55Yb3cezgEUzrj7tEz490IIVDH03KNAAJooo4o2FzG5FhtRnnte17t6CId47x5vmxKpJH/jCtx0z13/XBm1bs/d7fC4k/Lc1eJGv8UukNcsjWEc2SqqRdi0wUA4DNAyYrRUSE9bBlKOAAi/TzOnQvZmTSUe+9YKa8lMbfKN++RGX0a8YCiecU94CgdSIfOQpD3PQISA524Ft4U2yvnSMkDENNxM3MzV5+NRaEd77gHOPRBQlRrm3ca9FcMbh2zLLK3Tjs7Jatz1AiDtQw4oG9xTpB7MsnK35LYC35WUcQ/Th8+YCqkXs/nnMeHgzx/i16nNwPaXADMXeJISBbTtxQuoJsl3CgDjQt1wds11z2c4RGihLLNl7ITDyG2jJUoo53XFx/DZfF9ajDHWuDciA6qR2jJRtIwSEwvL+Ahsl1+NKCYRUkiUAs3SOMaQpYsRvwb2rc6NQSqkl41FBJSqKIlOgJYpqcSHMKqnMsr7dDtWHdBCOzElHfr5g/w6tQm48GXAVT/lrXfeassQ0auJ6DEi2ktE74xY5o1E9DARPURE/5TvMLX9yGdoXCemTe5h/mPV9n4No7AgAht36lhL/vq0TE4XvOcVh3DusKzqjFsohnOPlB+IegCF9QAN3W78dRsGTyfJEM3SGCYyeO6+oqL2YYxkQxRI2xZKVVmh6n/vGfcwqkftx5J311Mhk2yGB+W5Twcr0otKyyR67kRUAvBhAD8A4ACAu4noJiHEw9oyOwG8C8C1QojTRBTdKTgrvKdwNC2zuXMAGF933sj05gkCsdeNtnGw0tuGSoXMa1DKuLvtnq/sK1TtOPdgKmSKgKrSeE+ZCmnaicnmAdosTWAWOZTu6x+2LeJYerykVJEVqqqIiVBWNE3YbMCj6urpj62GMOGwxGtd9Uye3BjcFoopP5DGc78awF4hxJNCiCaATwB4fdcyPwvgw0KI0wAghDiW7zB9JB1kAWBT5xCw+uJ+DaHYIKAtaZmSoeH04DbRQjm/C96J4LMB2WvTvEI1cmydNrB8AhibDf062nOPCKh6nnuyASIyrVA9OwHVUE2bTnRrwigEA9JVzpbRvoulZTzP3e53WAmHLZ1g3f4uOjc0LbQASGPcNwPYr/19QH6m4xIAlxDRHUT0DSJ6ddiGiOjtRHQPEd1z/LhdBV0aPfctnQPAmqFxD4PPuXeshcPIbaOFUn48ZBItY7CpROGw00+zgV7zrNCvg3nuuucuvcsozj2BlgEsZiFhhUQp0CpLWsZSPMxv76ftXP1uk1RIvQDNqbCeuyYcVo6jZXTP3QIBWiZtMLu5EBo7OZ8997Cf1X2YygB2ArgOwPUA/oaIZnpWEuIjQojdQojda9euNR1rYDBR3txoZwHTYg5YvdNq++c7CEALJV9+wNQACCFpmT547iG0DFnSMpFjO/4ov0ZoBlHEe99zt6dlAhx0CgjLyshWaQJlcu356jDP3WsqbpIKqfHcpQocLc+dQFq2TAznbtlHVQ8Kp65QjZCROJ/z3A8A2Kr9vQXAoZBlPiuEaAkhngLwGNjY546karONLTnJWDM07mFQAdUyOKBq7NxJA9wS5hkMkYhKM4RKpzPjqdV6ofCM+yWhX0cLh6mAatfsoq+0jH0qJIDoht5J+w2b/dgEVPWYTqmKstaJiSiBllFNsi0fUFbCYY3FcM/dst7gbCONcb8bwE4i2kFEVQBvAnBT1zKfAfByACCiNWCa5sk8B6qQRMvMtiXdf54018gbjqxQrZKlcJg0bu1B0TKmBjFJOOzEHtbnD+v+gwThMCATLWM8C7FMwWuW5YOmMW+8LhAhm2yVCqkHVCXnrgdUY2kZ5bnbcu6Q+zEQDgtr1o2Iat0CING4CyHaAG4AcAuARwDcKIR4iIjeT0Svk4vdAuAkET0M4KsAflMIcbJfgwZiaBlX3mwjPazQEJAVqqLsFTEZC4fJDIdWnhWqMbSMKXXkT8cjFjj+KLA2nG8HugKq+heRAVUzWsbkt1gJhwFoS+Muonq+Ju6XX7OmQgYK0ErMuevCYb7nHkbLjAb3awg9bpCelgk37kUVDktVxCSEuBnAzV2fvVd7LwD8mvzXVyQdY9+4p28q8D0FYq+7gjZKNsJhnueeY4VqqRzYto4Ab5sCiX1HTz0NbLk6cn39N4XSMmGeu1NO5dEOipZpldi4u/UF2JBnobMfi1TIQAFaqcqSv5qeeznOuOfluQOJVK6HSFqmmBWqhZMfSBIBGhOLcOGkmiZ/LyKQLWMTUO0ozz1H4bDYClW7BtmhQ3M7QGMuMg2S9xfctz/GqIDqErdiTGGFbWYhVgHVsuLc7Tx3PdPEg/rd1rRMpUc4rCxavL2w35i1iCmQ5x5vMzxE0TLncUD1nIKfLRP+/Zi7hGUaK2bu0gAQyJaxCqgOlpbJNaCqjF1McVsgoKp/4RUxdT2AYuSDe7aNwQiHtcpsoETTMqAqxeRCUyGNKlQ12q8UTIUkgIuYoh4WsXpDydDjBkk2w0OEuqcoaIVq8Yx7QnBkzF3GsmNe0fa9gkC2TJaAqight8lqAi1jFoQMydFWqMu+onHGXX8fFlDt9tw7rfTGfUB57h7nXrf13ENmPzapkF0B1bJexKSyZaKMe0bP3Vg4TKl7hgTaBYa0zECQJMA0KpawMjTukXCI0BYlbtZh6BUD8LxrzpbJaVCJtIyB544Yg5jCuCcKh3WPsdMIVzUMAXPuJr/FXjiM31ga97CgtEVANZCN5Xnu/CcRoYyYB2NtiqkwJQlgCGPhMC8wHk7LnLfCYecSkgSYxt3FoXGPgfLcuULVQjhM49xz4yE9WqbXc7fRc4+8EZVxr0UH2yPlB2ST5x5axqBptOksxFY4rCM5d7LMc3fDZj8WqZD6tpTnruu5M+ce8WAslYH1u4AjDxjtTyFQoZpGOKwRnfVUVOGwwhl3D1EBVXcZyzQ07tEgj3MvGXrFALybvC/CYZ0Qzh1m6Zqx6YMq7zstLYMuS1+q9dIybQPjPqCAKjkOFsWIdSpkaEDVhnMnbWOlKhxoFaoEVEQrXs5g4xXA4QesZBTCeqjGHntPujmClhka98HAiUkpGxvSMrFwNM7dLqAqK1RRhpPX1aO84hDPnUw9d2SjZSIrVAGmEDJ47nHXbRgE7DxGArCMkcyl++EBVZM2exot45QDnLuXLRO3vY1XAPUzwJl9qfepYJwKqSisYbbM2QXFFN+MiSWslIZpkFEg4hz1qqctY7iBjpYtk1tANcdUyLgbMZVxD+67Z5w9AdWmUUDVTFvGzmN0iNBEObzfa6r9hgRUOw2ASkaNM4h6aRm/QjWBlgHYuAPA4fsNRs9QR9kh0miZmGMfQ8ucz/ID5xwiuUshMCaWsTKkZSJBIDSh9NzNAnwAPO+6LUr5pRDE0jKmnnvMjZiGc494DyDGc08ZUIV55o/NA5QIaIhK74Mo7X617XgwSPn0xoFu+QEtoAqV5x6zzXW7+PXEY0b7BfTZh/Y7YgOq0R21hrTMABGpidJcRAkuVkpD4x4Fh4C2KKNEAiV0LDx3PVsmr4CqTKuMpGVy8nbrczztLkUXZgezZcI895A895QyuOYFWfaaJuy5986EUu3Xy5bRdu62/cB3SnRny5QRTIUsox2fWlkZ5QKxFfOWgbpwmFf4GLeCly0TwrkLu6yls41CGndQxBSrzgGz+pBzjwTJClUAKIlOBm2Zcr5T1TDDCYuS/Thvtz6f2J0rMlsGYM81LM/dJBXSsPGIVUCVCE1UQo9nuv2G0DJu27yXaQ8toxUx6RWqcRiZtjLu3mHWPPd4WkZx7mHZMsM894GBgPA7XmZDLDvhin9DMFpScaSEtkURk5Ytk+cVX6qE0jI2JfuR3m79jJlx7xljdaCpkLZcr0PIxrnL18D5ddupH2IKgfu0VEEJLkpwve8Ss2UAYHSGz5sh9HaLHiuTJs89SvJ36LkPBpG0jORU60NaJhKO5rmXRds6oNpEOd+pqlMOp2Vg0VQ6LqCaaNwp9D2AiIBqei7avB+sHddLAJpZOPewPHe37Wc1pYSjp9rKB0MZHfkduIgp6YExMuPHSgwQzJZJQ8tIzr0SElDFkHMfGIgAN+yOV7TM0LhHgrNlpOcu2FM2qlLVKlTzpWUqofIDMOapEwKqMcFUoFt+oOvL0ICqGS0ziE5MJLNlKCPnHqRlOsa0DN+n8g85uymj7Y2xkpQKCbDnnoGWSS35qwTgQvJ7WRWyeNa9mMYdEU9hmavaGHLukSDIKTvYcwcMa0Q6fciWAaRXHFWhmlNAtZGGc9c899AxdhcxGQRUYRo/sDvETMtk8dz5NavnTgh2YgKAKknjDqTn3DPRMtoxjDv4rZXIB81QOGyAiJzeSo6xbVgi/b0ExyEp+uVPkY2YGS2gOihaxpinjqVl4j33QCembvK+XOvlsTspDJTEoDoxQeW5W6dC+obRg41x1wvQ5OymIjTPHe0Uxt3Oc9eFw/xsmZiDH1OvkCVr6WyikMY9MltGiht1yCwf93sJBGjZMmxMjTJmOoOlZXIVDmsuJ3ZMCjTr6BljSEaPqXCYqZSCLeeOijUtE9qJycq4652YFC3T8badWKEKMC3TXAgNtsfBP84aLeNGLh4bGC9qJ6ZCGvfIwyy9qpYz9NwjQVq2jBUto2fL5HjBR6RC5iYcJoQMfsZro8SmQo7OAMunurZpUqFqKj9gB4cIDZHBcw/rxGTDufPW5KD4AajTMhWkoWVku0zDoKpfoer/jtjj2W5EjmXYiWmAiBRgUrQMDY17FAhanrsMbhlVqbp+tkyuzoxTDm/WkZdwWEptFIr5C5ObgKXj/gxDvaauUB2QcJjk3Mkyzz0yFTIHWsYLqEKgmoaWGVXG3YyaMRYOi9HlFzzgwqGQxj1SgEnewO7Qc4+EE5otY7CBTh/03IGYbJmchMNS6pHrxrTn901t5D0sHuW/PRnctKmQpt64sDrGipaxLWLKK6AarJHMI0cAACAASURBVFCVtIy85irE9Exi8w8VADc27vzKqZDys7gVOtGeO2xjH2cZhTTukQJM7TobLtNKuu8hEOmcu4VxV9oyyLETE8DT9ig9d4PNRHq7aT33AC3T7blv5Nf5w/xqqHEe4KBTwM0oHEaZaRkNnZaV565XqAJgbx2SkgGSH4yKljEMquoSCh4tE3ehx2jnxMpIn8MopnFHhEFqN3g6OugBFQgEQkt0ee4m5rPTgktl6IGqXBAVUDWmMqJomZSee8R7AL5xXzjEr8q4p2w9RzBN67QTDgNxERMJ1zgQCUTRMh2LVMhoWkZlzaSmZf7hx4CHb0q9b7/hiOa5xx36mHqFYRHTABEpwNRpoIlqIadQgwLTMirPnY2pqecu5E2e63GOzJYxD6iG0zLpmk3o6Y89v88z7kf41dhzNw+o2tMyqi+tufceLRxmWsSk3aeSXhl3uRbF89wTaZkZ//1tf5h63/pxTiUc1omuVxgKhw0QkSll7TpzjcU7D4NDQDiMvSezVEjfuOd6mCNpGfMK1dAbsWMeUO3ZzNhqHue89Nzb5rSMqXCYjcvoKOEwwEpfJlo4zFBbRr9Px1YDAKZcriL3PfckWkYrOgtppBGFAC2TpkK1HZ31NBQOGyBiaRmqFLJUeFDgPHf2wByPljFApwVX3uQDoWWiahoiEHkjpvTcY3+T47D3vmDJuWMwwmGcLaM8d/Ogqq6F7sGqQlX7veNrAABTLqc0lkXKY1cdA179h8DWa4BTT6Tet9BpGe+zmBWSdPkL6Lmbna1zBNEVquy5F7GabFDQhcO8VEhTWoYq3rZyQ1S2DGxombCAquLckwKqMbQMAExu8D13NRtI3WbPvPOVFS1DQCOD565iMMFsGfM894BwWGUMDVQxIY17RV57qeIV17wDaK8AX36flG2OrzIGfIdl5rFPwHHKANZb0TLqIVFEm1JI4x7pzbWbaKBaxIfswEB6EZNrIRzWacM19OBSIYKWMW9wEbG0TUA17Dqa3AAce4Tfd1Lyxtr2jHP2bfLcQWiKHDz3wGDs8ty9qlAizNMUJjts3BX3nqT142H1xfx66glg01WJiwsBbMZxbLz1N+Un/5SQLdMMPY9etW4B2YDi0jJhX0jPvYgnYlAgAC0RlB8w9dw5W6YftEyYnjuME91DG3enToUk7X3IAiPTfmOHtpnnLoeXfllLrlcVMQGw89z7IRwGYM6ZwqT03KfbJ/nDiQ3pNjZ7Eb+eTEfNCAjcUP6MNhY3BS3Tex51eqdoKKZxjxEOa6ISfnMPAUAKh3lFTBbCYZ1+ZsuEdGKCqZ57RPpgSs89IBwW9vuqE0Brmd8bFzENRjhMSf4CsMuW6YdwGIA5msKE9NynOlLGYXJ9uo3N7gBAqY27K4ArnCe9vzfiVLJwWBgtI1+LSMsU0gzGZcsw11jAMzFAZBMO65PnHkbLLJ3MTzgstfyA5rmHLVAd4649SlcG6KtwmGWau+a5m9My4cJhFnnu3cYdUx7nPt0+yfeqnuoYh8ooc+0rp5KXBR/nNTSH1tQFAIDtzpF44bB2I4KWCYk/FASpjDsRvZqIHiOivUT0zpjl3kBEgoh25zfEsP1EeJudJtMyxTsPA4POuTuWFapetkyeD9FuWuaprwF/cjEuaD6ej7ebWn4g/L2H6jggXN7eAPLcbY6wQ8SdmADLPHeVCpkxz72LljnjTGNd6xDwzz+JLfXHcBKrzDyE2pTXkCcRwsUs5lHfyKZoOx2NPvZCsGMRSsukH965hkTjTkQlAB8G8BoAuwBcT0S7QpabBPBLAL6Z9yC7EdlXU3ruRZxCDQp6DnTZ5RvftELVp2XyHFg5SMvs/xYgXLxw8SsWtEwI8pAfAPxc6+ay7xUbtdkzjB9YCoc1Mnju4RWq5vIDjhM0jmcgexs/9nlcvPRtnKBVZgOrTXp9kpNQbs6hTC6aay+HKNWwnY5EH/sUD+nztYjpagB7hRBPCiGaAD4B4PUhy30AwB8BqOc4vlBE8rDtBhpiGFCNAwHooIRFjKHW4imymefe7lNAtRqkZY4/CgDYvXQbUwIpIYBwdzd1tkwCLVMZ49fmojktA4v4wVmpUA1LhbQLqOqU3xwFUxhPICUlo1Cb9IPZCajWTwAA3LG1aE9fII17xMIxD373PA+obgawX/v7gPzMAxFdBWCrEOJzOY4tEpHpcSpbpoAnYlBQx2aeJlBrsRdk2mZPGfdcYxulCtMdihg99gjgVDDbOY4dnafTbyeSllGZLVk9d9nso7VsHFCFcVqnnVEhomx57qGpkJacu/b3MkYD39t57mmNO3Pz7vg6dKa2YAPFBFQ96eZoWqaIJiWNcQ/7Xd5RIiIHwJ8C+PXEDRG9nYjuIaJ7jh8/nn6UPduJCqg2JS1TxFMxGKhjs4AJVJXnbmJyGgtol8fltvIcmDQcbou59xN7gC3fBwBY5aYLogFxtEydb96EVCr92gn9fR4ts2TMuZv3g7XTNMlcoSpfe7RlSuadmPSfe4Kkp15leuYE+mjcG8pzXwNRm8IEVqKdmJhitNBjURCkMe4HAGzV/t4C4JD29ySA5wC4lYieBnANgJvCgqpCiI8IIXYLIXavXbvWetDR8gN11GGmf/G9inmaRK3FMqomVAHqZ9Cs8PQ6905MABujU0/y6/ZrAfh6JGkQKxyWQMkAXZ572GOiqmgZzbgbqUKmWhSA4Xnp2o8XULXRlnFDqAhr+QH/R9zrXIE/2PFR4NW/DwA4LkxpmanUxr3WkJ772BqIygQmaSX6eJ6ntEyas3U3gJ1EtAPAQQBvAvCT6kshxByANepvIroVwG8IIe7Jd6g+QtPj3A7gtljqtIAnYlBQBnnBmUS1dQCAYYXqyhk010jjnufAFG/daQEnH+f3214IAJgSBsY9SjisXU8V+EysUFW0THPJuIjJLq3T/CgrPXcAmTx3b89C2Oe5B7YrcGzkQuCCSbSpgj10gdnADAKqI42TaAsHGJuFqE2y5x5Jy0TPwM7rbBkhRBvADQBuAfAIgBuFEA8R0fuJ6HX9HmAYQvtqyhttSMvEQ1EN8zSJWpM999QXcKcFtJY8zz3X4+zRMm2/T+nqi9CBg0mRvn9mrHeWynNPqFAN0DKSq02plmjaDxbCshNT3hWqQsZBbDoxab/XdcFPjNUX4Xcv/xIexkVmA6tNyVhHskZ9rXESpzAFIgeiNolRakavF0evyfEX0aakOltCiJsB3Nz12Xsjlr0u+7DiEdpXU/JmdTFs1hEHdWwWaBLV5nxyWbYO2Q2nUVa0TI4D0z135Z2NTGOpNI1pN91UHEgQDkvjuScFVFW2TGuJrzmnksjje9sz7gebg/xApmwZNRBpFC0aZHf/XkV1uU7FnHaqyVTK5gIwGs/X1xoncVJMYYbgcfzl9lL4wjEprUWmZQpcodr1oea5F/FEDArKYC3SJAguJrGcniqQfSyV554ryjKTorXiF6rUprBYmsa0AS2DqL6jaT33iPcedFqmE174ErdxszR3S+EwIq9QLVOeu/rAM+7ZKlSFNhMx1bYH4Bv3FLx7rXkaJ8Ukzx5qyrgvhi/sBVR7Z2A9x6JAKKhxD0kpk3nMQ1omHjotAwAztGTsuXu0TJ7pMnqX+/oce1tOCculaUwaGPfIvqMpPfdgtkxMKmRzObJkPXrb5sJhtp2YAIJbqmbrxKR2bm3cg/epfm5Mq3UB+FK/KYx7ubOCJYzysVDGvRXhucektHqSvwWsjCymcUdIEFB57mLYiSkWynN3+IJfhYX0VEFd0TK8bq6HWU2zV04zLSNv5KXSNKaFCS0TJRxmkS0T9gNLFfbWm4s8TjnlTwPTfrBuZEVWwn7kwIVTs9SWUfID6gNZRJYxW0bAPzd8LAwHZuC5lzp1rKDKg/BomYj1YmkZyPEWD8U07nG0zLBCNRbqyCjjPkNL6T0oxblX+sC5613u63McPAOwXDajZQQivN3U2TIJFaoAe++tZWDpODCRPqU39LqNgbCsUPUcblvPXb56+/YCx6Y9VLtpGT88YSqiBsC7JtLoy5Q7K1gRNZ59jfC1Xon03ONomTAVtWKgeMa908Ya90QILSMDqigPtWVi4BUxSVpmGovpDY7i3Mt9yJbRaRnNc18sT2Ma86mtYmTf0XbdwnOP+H2VcebcF48D4+mNu2k/WF7HcAVoQcvyKI/TED3yA5a0TKATE4IzEVOKCoDmuacx7nWsoCZpGb6WIgOqcZXGXraM2VDPBRTPuN/xp/i7M29Fye3ySBTnLoadmOKgjs2iwxf8KlpM70F1e+55Dqzbc5cdepZLMyjD5c9SILLvaLthkS0TsVB1nGmZpWNGxp3ITF6ZFX/tKlQBoDmyFlg4Yrx+ngHVYEaM0Dh3s8whAIa0zApWUOUHlFyvEhVQbUcXo/nHonhGpXjGfXIjAK2Ti4Ly3EWpkCdiUPBoGeJ87WkY0DL1M0BlDG3ZQzXXw1yuskdcP8PTbultLZVlG7blkzErBxEdUDXMc49aqDoONBaBpRPAxLrU4wJMK1TthcMAoDm2Dlg8arx+tOduVv3dHWPQ5ekJFgVCtZQB1U4LJdHGilCee4Jxj5EfGKZCDhLSuM+0TwQ/V3nuwx6qsfCOTamEjlPBKDXMsmVGZry7MveH6OhMT0C1Likgr7ApAdF67ik9d/191IVUHQfmDwKiA4ynN+7m/WDthcMAoDW6Dlg4bLx+j1iWF1A149zRRb1wPIRCv0uF6jivmPSgl52yVlCFQwSqsSMTbdyVcFhYtow33MKheMZ9ahMAYKbTZdxbKwCAhigXsmvKoKCODYHgOjWMoGmW5z4607/WYyMzkpbxPfeGo1QY03HHscJhKTz3ROEwgI3Mqaf4/fiaiIXCtg0jd1XYNsiWqzRG1zGdJe+N9PvlVycHzl2/tPSZSPd3qUAEbL0a+MZfcjOXKMjfu4IRgADHcbAoRlCJLGJS2jLfe8Jh5xYmuaHuqm7PfYn/PiWmCvmUHRS8aTEBbmkENTTj24/pkJ67LyzVB8994TArQ0rPve1Ib6qVrk1AqHCYEHzD5yEcBrBxV1N5A1qGYCYGZt0gW742R+XYDHn3HirClnNHMMbQTcsYc+4AcP0nuOXeg5+KXkZ57jL+RiAsYhSV9nL48nG0TIFzIYtn3Edm0EANq7o99/kDEOURnMLkkJaJgee5E6FTqmGEWuk998Y8UJvsX9XeyAxwZp98z1x7W02VWxE3ZhdCqzqXTnBGhJz1xYEi/9AwvcV/b0zLDEA4TE45GqOy+bShce8ZYZYKVX272kzEqogJAMZmgZltwEJMLMHz3DkVkghYFKOodMxpGYUimpTiGXcinCqtxqpOF+82d0DedFTIKdSg4JV/w/fcUztQ7QZQGemdtueF0VXAsnxo16Rxd6S33U7nuffosTzxn8Ch+/j99NawVQKgNLTMzlf57w08d1PhML1c3wRqlfqIzOQx5N17aRm7IqZu4TB9VmXcclDHxAZgMeaB1WRHoI4qCLxP9txjaBmnHKoR1LdrfQAwO1vnCE47qzETYtzF1GbgYDGfsoOCXv7tlkcwgpaBcWfe2u1XlGlU0/fuoWVS8sYCIHWPzh8GPv5j7OkBwEwa466/j/iBW6/RxmmiSU6GbfYsaRmdcweMM2b6JRwm4FNdphRVAJPrgSMPRn8vZ3nLoubRMgtiFJubJ8N5u04z0msfZssMGOGe+0GIKZ4uF/FEDArq5nKI4JYMA6pdJfy5H2fdUMqAaqck95fSuAf03Pd8AYAAzjzDf6fx3CPeB1AqAzMX8PuUipCAeVVmFuEwAGhVpplHNvXcve3IN5a0TLdQWmAmkuXimdjANQZRvXVDaJlb3SuwdmkP8MAne5dvNyL74PYciwKhkMb9dGk1Zjsn/Cun0wIWDrPnjj4E+s4jaJlocEsjGCETWoaDkv2jZTTjPraad6mMezudcQ8Ihz32Bf+L2lRw+xFIFA5T+Pk7gF97JNWY/O0ZLR7dVSoBahUBAibWA4vHjPcL5JQto8EN0DJqXxbu+8R61piPSonUUiHVOD7aeQ2OTlwK3Pmh3uVXTkVKCHvCYQW0KYU07mdKq1FFy69anD8EQAw99xTQA1rKc0+dtSBzxXuEpfLChueyp/mS3wDWXgIAEFRGG056z10Jh3XawJO3+YUvKbx2IGWFKsCFMSkCtIFtG+q521eo8joCgpuLNCMCiRHoFQ7LKVtGm4l4Egk21MxkQqC4pTh3RcsALhwcndglbUUXFo9Fxk6sqaNzAIU07nVHan8r3Yw5bhfnKuM+ZN0j4c+KiT13NNORMkIw514Z7d9Udds1wHuOA694jz9eIjRQTZ8Kqca1copT3C68jr9IwbcD/b12Bi0cJgS452szXaaRt1/5mplz76FlgjEf/szGc+d06MhYgnQEloVPywDAUmW1vC5aweWX4jSC+pT2OwAU0rg3lXFX6XHyaezK6tUiivwMCo7nOXFAtUYpA6paE+FBZhAQAXVRM6RlyK9o3fFSflVB1aT9aXdE3r/PVDgsUuEyAQGvuDKWOo3U22+fhMP0mYj3ADLaokRKz33Fy5bhnS1VZvn7pa406hjP3b/WbQZ6dlFI497wPHc53WwwPePKgFwBH7IDQyBbxkuFTHGLqVREPVtmACAQ6qgaZMvICtUV1Yf1YuDlvw1c+ZNxa2n7097nfR2RWeGOG6VNn7Qb3SuujhsrQ/bKD0jjHhF0jBuHTmuILuEwwLKQaUIa96h0SHmt1JVwmByLZ9x1j7/T4mslol7Br2EqnlEpZCpkg2SQTU035cXrlrm/ZRFPxKCgjoxDBFGucSpkmhU1z93b1gAOs0MwMu6et6uCbWOzwMt+M/X+UgmHWYIAI1fVOqCqe8U2xr2birDMcwd689y7PWArP6EyynUQUYHi5hLaTg0CTqCuY9Hz3I/7yyovPkKX3z8WFuM8yyik597s5twbiwAIQjYvLuKJGBSCAdVR5txT0TLKcx8daAYBKeOeuohJBGmZ0Vmj/enG51ygZaxSIeVjSQhhScvwq3csLDl3p+tpps9EMh/b2kT0Q6u1gpYsflPHzyFizh0Ieu5L8gER4bkPaZkBo6mqFpWYVHMJqE5ACGW4CngmBgQ/FZIgSjVUqANXDzA9+nng+GO9K3rGvTbQ1mNEipZJKT+gCn8ULSNTKlPvT+/ElPMPJENaxjag6tMykJ67mXHvoSJy0nP3gt3aqzXFVxlLZdz1sSyUFS2jefyL0ouPzJbpm9hG31FMWqbbc28uANVxfwp1lsZVBPjZMoArKRZShrvdBD71NmDnDwBv+sfgihrn7nGyA/LcV4RBtowXUD3JBVfVMeP9+e/z/X0E02wZu2tZecUCkJ67Keeel3BYiJ67FtBXn1mhGjMjaS2jXRoJeNsEYoNfnQS++ntMx1THgK/9MS8QkS3jX+uW4zyLKKRxb3nGXePcaxOFnkINCl62DBFEmY8jdaThPPpdTh98+uvMs+rTcMW5V0a8h+ggjrMXUE2dLSO93eXTxpQM0OtZ5onu7JEkCNjRFwGvuDrGxrndDJW0jdqvvh37bBl/W94Dw/tOewDZoBIzI2kto+WMBOMnJPnzkWlgfgH4xoeD6yR01BoWMQ0Ize5smeYSUB3XdCCKdyIGBb1CVZSU5y4N98F7+bV+ple7QwU0yyP+tH2QnLuBJjkB7LkbUjK8btCzzBWE9PLKyN6JiWkZblRhUsjUmwppF1AlIk8yt9sDzkzLVGNmJK1ltJxaT+aTEADmD4Svo1r4daFvBXsDQCGNe4tqcEH+tKyxyJy7/H5o26Oh31zKc3cU5XLwPp62AsBTtwVX9LJlRjLMpc1BsKFlwJz7WHhJeez+ulL18oRpFpe1nrueLSOTDEyCqpGpkKadmKB57mpsXb/I+lKqxBRnNZXn7n/kUUTbXhhcdvVO4HlviTQaRaZlCmncHYe4y4rHuUvjPkCPsqjQsweEFAEjRXkcvBfY/mJgagtwrEs3RefcMTjqyzENqCrhsOVTVrSM4x0f41VTbNtQOMy6E5OWLVOVnawMgqpqhB4VoQLuGTox+RlWXdu25tzHo6+Juf1YqKztkW8WAsBPfRq44V5/2R/8XeB1IXozCA5vSMsMCERAnXTjzrRMN683RC9Ie/WNuzTcp58G1lzMueHdPUs1z921NDo24ApVk1RI6WUtn+TfYbo/77UPnntX9kgSMguHqQpVwCiomlsnJi07yM2blonKllk6CSwexbHRi7poGSm3XJvga1xx7Gt2xu6mb/LWA0AxjTuoy7gvymwZ+X0BT8Sg0C0cBkjj3lrhYOroKv63cjq4ovLupfzAoA4xAVhWnHsKQyCEgCMExw1sOHf96ZczCBadmDILh0njbuK599Aylpw7dFomyON7DyCjLWqIypY59jAA4PjYRV20DILHfu2zWcNdSTdHoMC2vZjZMo4D1DHin9zmErd/G6DmSVHhVewRedrs5Da4PyrAmupjs70BVS9bZhQCSwM7xkSEhqgCEF4nqDgIACNihSVhlSKk4f6APtEyzmA6MQWFw1RANb3n7jtJyr1uAyDzIibHr1Dt5q5VK0DrbkyVcZ7NdWd1STrx2MiFcMg3/l5AVeHyNwLrLk3xmwZXsJc3CmncCYRlGmGPXQjPcy9yZHtQ0LNBREULqNalcR+difDc9SKmwbnuRL4uN9orycZdAFU0+Q/5+2z22R8JC4tOTFa0DHnr29AyoXnuxtIDCEj++h5w0HO3ltStaoFiPdPl2MPA6CrMl1eDoBv3rrZ+z3tLqt1000lFQipahoheTUSPEdFeInpnyPe/RkQPE9EDRPQVIoqf62QFsVYzmstyuu4OaZmUCGTLyEYYTrse9NxHZ9m463l7LT+gioHSMjKgqo8hBkIIVIU/y7DbZ3+uId6mCS1jG1CV66s8dyAjLWNn3EFhtExwkCY0VQCViN917BFg7aUQoJA8d3N0P5SKhETjTkQlAB8G8BoAuwBcT0S7uhb7NoDdQojLAXwawB/lPdDAmACsYJSnmmq6WZ3QAqrFOxGDgjoyDpHnBVOY5y5coDHvr9iuA1QCSpVAGXm/4aiAKpAqY0YAqAnpuZfjvfwoEPXLbx9MhWowFVJmyxilQobkuVt57nq2DL/qQl7+IC2gsoC6ZyRLx4CpTQEFSrU/GwaoZxZTIKTx3K8GsFcI8aQQogngEwBery8ghPiqEEJdPd8AsCXfYQbhEKFONT6xzQX+MJAK2c+9Fxt6QBUqz70TwrkDQWpGNscGFBc8uGyZFdT8MSTAzcFzd6g/HOvAhMMCqZDKwzXn3D2+vxPdYzQOXKGqsmWCjlf2CtUIz70+B4xM9zwYTauDFYrMBqQx7psB7Nf+PiA/i8LPAPhC2BdE9HYiuoeI7jl+/HjYIqlABCyrPHfPcx8f5rmngF+hSnClsXY6K5rnvsrPD1/R0iFliz1AcsEDGy+hDmlYUtEy8I27recO6k+2DJml/lkLh3nrQz7Ayci49wiHSXkP43FoqZ/dRjJ7hWrIjEQIoD4PjEz11AiYpqEqdD+UioQ0xj3sV4UeJiJ6M4DdAP447HshxEeEELuFELvXro3XckgaEBcxLfsXbW1iKByWAgHKUxprp6177tN+s+CA576iee6De4ASAUtCeuCyKUsc2LirgKqZaJi/0/5cQ4OiZQJeseMYy/6KbtJdFgmaQhcO63a8Ag8gG1RCZiStFcBtseeO7ownsttXgdmANETaAQB6A8otAHq6zBLRKwG8G8DLhFCuU3/gkMyWaa3wkxoI0DJOIbP3BwOVgkZEcJwSzohxVJqngHqbGyA4JZ+WWdaNu5+G2M1n9hMEwiEh89XP7I9fGGyYKsIXObOBQ/5xyhM9GRsJYANlH1D1vOJqjDxuDLxDIIsEbdb3Aqp5C4dVQ2QVVIxoZLqn0K5bWz4tzvcK1bsB7CSiHURUBfAmADfpCxDRVQD+GsDrhBAR7VFyBAErQvKwSmw/kApZvBMxKJD2SgQcFGswunSQPffRaf4y1HOvBz33QY2XgMNYDeGUuYI2AQI6LWObLdOngGp3rnUCchEOA9gwNxaM9gtos7Pmkt0siKLz3NUgXdtcSC9QrAnK1eXMrjbVc40SmYm2KfRU6xYIicZdCNEGcAOAWwA8AuBGIcRDRPR+InqdXOyPAUwA+BQRfYeIborYXC4gEE4S90v1NFCGwmGp4BWRkPKK12Bk+RBz7rIHrfeqc+6tukfj2Gqe2MAhwIUDMb0tnXEP0DK22TL9oZ24QjU97LNlurziVduBk3uBL70XuOsvU+0X0PYtm+EYj8PbnujxgDMf3bBAsTLuIzM9wWjT6mCF875CVQhxM4Cbuz57r/b+lTmPKxYOAY9jO//x3X9hD21qE8RJ9tiGAdVo+NkyBIeAA2INrlt+FFhZzWmQAFAqM0Wj68u0654n3Mtn9nG8qiBn5gI4KYx7IFvG0nN35LHJGybCYT3piAYI5LkDwLrLgHv+N3BiDz+4r/n5WA+oh4qQRYKm8KgX0esB699ZIUztUlG0XkBVH4tlKqR8LaJNKSQ7TQQ8hc1AqQosHAY2XsH518MK1USE0TKV9hJw5hnfYwe4NPuBT/ot9wLZMoMVDgMAN63nDmicu30RUz+uIpOMjSxpvT20zPpd/HBuLQMLhzz9lSj0UBGWnLvO/Xd7wLllywQ8dz8poJeWMasOVjivaZlzEQRCm0r/p71zD5KjOA/479u70+nQnZ4nTi/0FmBeFkIQTPwoLAxGLkeEkCqwnVApKiR2XHHK5VThOKFIKlUpO4FUOSGxE5vEIXF4JS6TAmKwjR2b2AKBBBLhJWwskIQEQu/j9Ljt/NE9szOzM7Pz2Nm9PvpXdTe7M7PTX0/Pftv99fd9DXPP0DsWXgDY7ZPaKYI9dxB2qmF94PDuRs8d4Oqvau3wo1v1+xOjvrLsrM1dmOSiQAAAE8pJREFUlzQ+c4k2E42le8yEXCELKnekmmdIcvha+89ygTvd6BWbq4ycHT7hpYfTy46KWFS5e9cjJnGYNI4VomcKSE+k5x62uUcnQQtFw1pslrFSufsJmEbO1TsWrgGCUXA2NkVnaCQO06995Q46j7vHrKUw75yGh8qRPf4iwkWDa4rglVKfsUS/2P+L1POVUvTVj+FF0xahKrNMngnVaP7zvOVAYJQw90yQGkw7Vb/e8dNM16mJ6OjUE6OFbO6N5GAxicOiP0B5EWle/LvJWyYoC4V+SRpLStqnU6xU7ngJmBacr98uWgvYvSRWpwgmbhLRE6o+0WRKQ/N0j/7kcTj6BgzNB4oH1xSS15QzPk2XzeHXU8/XZpnjxXvtpsxK8rnn8LUuk7DKk90vqm9A292XvQ+mL9RtmVZ2PWCK8HrGBXru/vVCZhmJHCt8WR2TEfToGjuo0yT0DZh0yQ0EKWQC8jxsLNTtlmaF9HxWL7geFqzWvUzK2SnfKQQjBAV4E5MWd8U6mD4/fPLQfK1Mj7zeeE+HzTJmW/fS9x5Pd+nTPfexwtGpXplVJQ7LbpYpPqHq3bRQr/g3vqUXyH7gc9pzJrXswGUCEeC5xQiIHk0c1qhXCe0+e1m4Lib1gOeC2ZbEYd7nLewyWtlz92e+e/vhtIv8/aW+EO8QgsPimrG7/9flP4aP3dN88tB8HZm693n9fvoCgFh7ZlX4Q3vPLNDCX9v3cy8anUrw3rSXPB4bZToqsaacwbla8Z0yp3mVrYSyayKhxHz55Wg2yzSW2QuXVYjh0+HNlxoXGTuk60hzoF1NikWoTvbEYROOpCGWzT6pncI3ywQmDcemzNTuj1G8nvyup/TW9NyLBtcUwStmvM9T7kdSz6/XFX3144V93KG6L3Ie00CZVLNe5ya2rGlz9OjnZHIQechDpEzPPXC9aIBhKOd8UYZP1ykpNt0Bu7bonrsZ4emee1iWQmYZi60B1ppl4prJJQ5rTePWNO5R4iNvlDk7w8o9as+sFE9ReT3xDD333vqxUmYZkArNMtloyn+epxzvGnGFeUsPjr7VbIbzyzbXCfXci5tlVFAWiR4rod299U8f+CysvCzcc1cq/MNY0CzjfcqZZTpE0hCrMbPdYYEswhsqh/KnJD31Q/P0dueTer1Jk3Omk4nD/LaUHm0aaBVGr4yfe4kJ1apS/koO00DUjJGH1LwtvnLfl1p4yMcdSpplVNO+tpllPHZt1nNDJnVGNNCuVtDobnO+KgtF1sQNsWweQnWK6IQqpAxXvZ776Jta0Qe+rJ0zywSG71MGwwuIxFD3XCHLTKhKNc9Qo0fdWsuUyZOUGiCUQbmHUjofN2awUt4ycZ5sKaajrAwtaLwe3QcHdsCSS/wyw+kHSpplXM+9M4gQ+yvsVmJqjd8RDpgeEh/5voFG1Or0xhdJ+7lXJGCE0PC9f6ihbBJQoL1lSkyoVpk4DLL1VtsRkJdulnkzpexABHIps0zjAYtOEOe5F4nUavCR2+CyWxr7Vqwz11VNicOKpR9wE6odJWlFm3Z8ISY7EjTLZMnvsfhivV2wxt/VyZWYQsP3/qHWNnfl+bkX77lXuRITZLMOlJk/SpU9aHNPKTuU7hcKmmXM9WhOHNa2+3vhDXDxp6DWBzOXwJwVpsywHii8ElMJ81i3sXNClfghVplkS+8Ugn7GLc0yANfdBePHdbi3f34n/dwDw/cMyr2ulJlQLRPEJJVUMHi/e1oUUCZPkm+WiXNF8dI5j+6Dw3sA1Zhb8eULjH5LmGUa9W3O0ZLp2ctKbz+svk7b4EOmw/DdK5NbxkYfPDuVe8IQy7lCtkZi3qQ+8yJ+wjAP3SvqrC+k33M/mmxOgIC3TImee6DYtpLLLFNi/si37ccd7OnTHiWj++C+39J5/D/5WKggReDX+/hRPZkdeQYyyRFQtNE6t8UsE+RX/ib0Ni5xWJmybOwvWmyWiem5B4474gkOi2sFv2GdnFANtWWGnjuK0j33Wq06bxnI5v5X5lluaW47ZRgO7YLXnoC9zzZcXQOFN63CVEgO/3L+/6hZpl26PUo00K7wSkzBgC7LsFK5I/FDrFBODEcs8d4y+a7RlfQDvlmmlbdMXXvLlMktU6GfO2T7LS2TajbVWwZgxkJ4+VFtbgPYfGdT2b5Z5tBOP2FcUUF0EFNYNlrJWJJooJ0k6Iws1wE7rQFWKndJiEiw1zrWOaKJwyB/Zr6QN0XFhBTilEFtA06Rt5dxaoyXjlCtKnEY5DTLFCmnVa/47KvhhJkoXXC+7sFHyvabd9/LMGdVASkCsge9ZQLPn1dWFUQvG1ysO9d1LHavtlO5S/zQ1uK5j47h55apSbptNoVO9txDw/f+Iaif1AtPJDCgzJqafcX9svUPX+GPJ1836NbZAv+cgoJod+GEcs75Ne0qOmupXujmSHjZY4Vp33od3nrZ90ApJIO5XjTAsC2Jw1KImmUKxjBZnTjMygnVpARMNude7hRBb4WiS511MnFYyMTQP6TfHDucaHYZxqzGU9SUgK5bFSOTXFGZJV3wtEdZwsGp0+HyP9eTpAd26BTA4yf9/EJ++x7aqX9I56wsJENomb1I6ty2RKimEJ0XKjqh6hKHdRiXOKw4EtTuBe2enUwc5qG9ZUza35RJ1REx+b0DQVe5qcYTMuzW2YKykZEtV3268AY4/xMwOAKoUFBT3Ruaeel0Cyr34HyJihhN25I4LIWon7tQzL7vzDIdJmmI5RKHtca7M7p3WuwaVXk4xFELDt/709P+KqWYJyY4J+K7nYfqzTKtKRsZmTm98OCI3h7ZE9otUF65B80ykZFIYy6lKrOMajLLFLqOxWnErTTLJA2xXOKw1ngPvG1mGd/PHVKUO5yK6bkPxWc8zEJV+dwlx/0uGxkpxEdxN+Er94bdXSnFpbIJfvxNPXdR8IcybsK+bWuotsCfNzAUz+fufb4dUnUWO3vuZhv91XeJw1oT7wo5cf3cQ8P3aXP1zoSl9upKMU/2M9Y7vfwye1X03M226sRhXmGZ2tWbmwjc07qCq9X39EIt6/+y+KSu2SrV7FKYx0RVhHrELlPULOMSh3WYJH/hxpfGvoboFMEvV9HeU9SeWSUhD5PZy0Fq8OaLiXKNyH6OTJlbrsyJkDis5KMskK1hY8wyCsUSdsPS98L5Hy8mAOGRSlLisKq67i5xmKXKPSm6TfnHOyqOVeROHBZD1xKH9fZr970k5a5gRN5itL+kcpdqbKxFojKL3uek5HpN9E3V6QgCZhmpn2QRewrb2hsy6G1a4rAq52+CeqDlBHMCbkK1wySZE1zisNY0IgQb9yi/t0wnx0aR4bu3bmYMdaUYkQMc7S/uBgn6+alyQjWbt0y5yEiRhMRhcQyOhHruM47tppfx0sq9Ud9AnSV6rCqzTDjQLtU1NAWb04jbqdwTzTLmeGfFsQoJbEt5y3RjQhX00mr7tkN9vPnk+jhzOcDR0maZqlwhNVUnDvPKyqzLBkfg8G7/7fCxHfpFWeXuR+Sqpu9mnntRhObEYRQaJvjWMQuViqXK3RvSRXvueuuCmJKJSxyWO/2AUh0zfTW15fDpMH5MB9+cGIOH/wQO7tTH3niBXqlz+JRF5cq0PHGY97nMzTrvPNj9NJzQ0b3Dx17V+9vUc1eB/763VsVmmebEYWXzudunUyxV7nobfXjLJFt6pxDvLZPvGl1LHAYwfIbevvYEbL0X/vfL8JPbAai9+KA+NPuSkmXanThMfzCHyWPlOh2J+spjAJw69nMOMuivmVuUoCtk1JOtE2aZ4EMq4hKHWUFSAqZOBtfYSnzisHzX6FriMIBFa3Xv/Ue3wsav6n1b74HxE/S89BCb6yt5e2obJlQnSOKw4mXlYMkles3Z7d+FsUOce/AHbKytLidAQIaQt0wkcVhl3jKE74FLHJaAiHxYRF4Qke0iclPM8X4Rudsc3ygiS9staJDgLHwQm4dQnUJMi9dqkngfW6GHvG0WLIEmr4paD1z6R/DG87BnK7zrozo3yj3X07N7M4+MX1C6/auaUM13v8NmjNxl1XIos74BWPYB2PJNeOCzDNSPcnfvhkLlhmQItJ0ni3cPajlMVIWImGWSotozXMZ83j6d0lK5i0gPcDtwJXAWcJ2InBU57QZgv1JqJfDXwBfbLWhYJr2NDrNsTvLTKcKTTJ4nSr5rhPJ9V03c8P2sq+Bj98INj8A1/wyrPw4vPMCJlVfwT+NXtKXIar1lWp9bNiAvt3fIR/4KBmbA1nt5fMaHeb5Wzt4OYdOLL0vULFMvXUwszfncpZy3jIU6JUv6gYuA7UqpnwGIyF3ABuD/AudsAG4xr+8D/lZERFWUOCI4Cx/E5pntThEN/wZy2wBUdMxbIbFeFSJw+uWN91f9HXzwjxntmcPb275bupdVvVkmw4RqycjI3H7dMxfDjT+EY4e46+EDyCvJC2hnlsFslQoEA0USh1U2oUr4GRdPkLzXsdgDL4tyXwi8Gnj/GvBLSecopU6KyEFgDpC+4GVBvEbbcPtj9ARa8PDYSX3cyqboDNFhcU3gXzfu4KFt8SH9cby6f5R3L5pZhXhNeHJ+5q7NDPT1pJ47rnRwU9nWr4lQq2A2yntUP/G1jfT1pBdw7GQ99Jm81ATu37KLjT/Lr6RfPzjGzGl9xQoO4P3I/s6dTzIeWSXN29787W186b+fL11WlB1vjXLh0saEcE3gudcP86HbfpjrOvtH9WpVNpplsij3uFpFfwKznIOI3AjcCLB48eIMRcdz6Zmn8sxrBzkZM6abMTCFZcPFF2qY7JwypZc/vOIMrjhbJ4P6/XWreHFPi3VJI6waGWT9ucUTc+Xh3Ytmcs0Fixg9fjLT+ecsmMGlZ5YLYvrt9y2nr6f9X+ZLVg5z1eoFHB/PZotYu2QWa5fOKlTW735gBU/t2F/os6tGBnnP8jmFPhtk7dJZXL1mIWMndEzCJSvmcPYCnbZ55amDXHfRaRx8+0TpcuJYNTLIR89rpH2+7qLFDExJ7xwkcdqsU5hW8LPdRFoNEUXkPcAtSqkrzPvPAyil/iJwznfMOT8RkV7gdWBumllm7dq1atOmTW2ogsPhcLxzEJEnlVJrW52XZfD5BLBKRJaJyBTgWuD+yDn3A9eb19cA36/K3u5wOByO1rQ0yxgb+qeB7wA9wB1KqWdF5M+ATUqp+4GvA3eKyHbgLfQPgMPhcDi6RKbFOpRSDwIPRvbdHHg9Bvx6e0VzOBwOR1GsjFB1OBwORzpOuTscDsckxCl3h8PhmIQ45e5wOByTEKfcHQ6HYxLSMoipsoJF3gB+UfDjw1SU2qALuLpMTFxdJiauLrBEKdUyr3XXlHsZRGRTlggtG3B1mZi4ukxMXF2y48wyDofDMQlxyt3hcDgmIbYq93/otgBtxNVlYuLqMjFxdcmIlTZ3h8PhcKRja8/d4XA4HClYp9xbLdY90RGRV0Rkq4hsEZFNZt9sEXlERF4y22IrNFSMiNwhIntFZFtgX6zsovmyaadnRGRN9yRvJqEut4jITtM2W0RkfeDY501dXhCR8gu1tgkROU1EHhWR50TkWRH5jNlvXbuk1MXGdpkqIo+LyNOmLn9q9i8TkY2mXe42adQRkX7zfrs5vrS0EEopa/7QKYdfBpYDU4CngbO6LVfOOrwCDEf2fQm4yby+Cfhit+VMkP39wBpgWyvZgfXAQ+hVui4GNnZb/gx1uQX4XMy5Z5lnrR9YZp7Bnm7Xwcg2H1hjXg8BLxp5rWuXlLrY2C4CDJrXfcBGc7/vAa41+78CfNK8/hTwFfP6WuDusjLY1nP3F+tWSh0HvMW6bWcD8A3z+hvAVV2UJRGl1P+g8/UHSZJ9A/AvSvNTYKaIdGZtvgwk1CWJDcBdSqljSqmfA9vRz2LXUUrtVko9ZV4fBp5Dr2lsXbuk1CWJidwuSil1xLztM38K+CBwn9kfbRevve4D1knJhVttU+5xi3WnNf5ERAEPi8iTZk1ZgBGl1G7QDzhQbhHQzpIku61t9WljrrgjYB6zoi5mKH8+updodbtE6gIWtouI9IjIFmAv8Ah6ZHFAKeUtCByU16+LOX4QKLWQrW3KPdNC3BOcX1ZKrQGuBH5PRN7fbYEqwsa2+ntgBbAa2A3cavZP+LqIyCDwH8AfKKUOpZ0as2+i18XKdlFKjSulVgOL0COKd8WdZrZtr4ttyv014LTA+0XAri7JUgil1C6z3Qt8C93oe7yhsdnu7Z6EuUmS3bq2UkrtMV/IOvCPNIb4E7ouItKHVob/ppT6T7PbynaJq4ut7eKhlDoA/ABtc58pIt4KeEF5/bqY4zPIbjaMxTblnmWx7gmLiEwTkSHvNXA5sI3wAuPXA9/ujoSFSJL9fuA3jXfGxcBBz0wwUYnYnn8V3Tag63Kt8WhYBqwCHu+0fHEYu+zXgeeUUrcFDlnXLkl1sbRd5orITPN6ALgMPYfwKHCNOS3aLl57XQN8X5nZ1cJ0e1a5wCz0evQs+svAF7otT07Zl6Nn958GnvXkR9vWvge8ZLazuy1rgvz/jh4Wn0D3NG5Ikh09zLzdtNNWYG235c9QlzuNrM+YL9v8wPlfMHV5Abiy2/IH5Hovevj+DLDF/K23sV1S6mJju5wHbDYybwNuNvuXo3+AtgP3Av1m/1Tzfrs5vrysDC5C1eFwOCYhtpllHA6Hw5EBp9wdDodjEuKUu8PhcExCnHJ3OByOSYhT7g6HwzEJccrd4XA4JiFOuTscDsckxCl3h8PhmIT8P5g3x6eZ9nHAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "probabilities = posterior_decoding(casino_params, np.array(casino_observes))\n",
    "# calculation of accuracy\n",
    "estimate = [probabilities[i].argmax() for i in range(len(casino_observes))]\n",
    "same_indices = [i for i in range(len(casino_hmm)) if hidden_states[i]==estimate[i]]\n",
    "print(\"Accuracy: %f\"%(len(same_indices)/len(casino_hmm)))\n",
    "\n",
    "# calculation of plot\n",
    "plt.plot(np.arange(len(casino_hmm)), -hidden_states+1)\n",
    "plt.plot(np.arange(len(casino_hmm)), probabilities[:,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checking Accuracy of different versions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current stage: 0 of 1000\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:7: RuntimeWarning: divide by zero encountered in log\n",
      "  import sys\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current stage: 100 of 1000\n",
      "Current stage: 200 of 1000\n",
      "Current stage: 300 of 1000\n",
      "Current stage: 400 of 1000\n",
      "Current stage: 500 of 1000\n",
      "Current stage: 600 of 1000\n",
      "Current stage: 700 of 1000\n",
      "Current stage: 800 of 1000\n",
      "Current stage: 900 of 1000\n",
      "Viterbi:\tAccuracy median 0.800000\tAccuracy variance 0.004030\n",
      "Posterior:\tAccuracy median 0.820000\tAccuracy variance 0.002717\n"
     ]
    }
   ],
   "source": [
    "N = 1000\n",
    "accuracy_viterbi = np.zeros(N, dtype=float)\n",
    "accuracy_posterior = np.zeros(N, dtype=float)\n",
    "for k in range(N):\n",
    "    if k%100==0:\n",
    "        print(\"Current stage: %d of %d\"%(k,N))\n",
    "    hmm_test = generate_random_HMM(casino_params)\n",
    "    hidden_states = np.array([x[0] for x in hmm_test], dtype=int)\n",
    "    observes = np.array([x[1] for x in hmm_test], dtype=int)\n",
    "    # viterbi\n",
    "    estimate = viterbi_solver(casino_params, observes)\n",
    "    same_indices = [i for i in range(len(hmm_test)) if hidden_states[i]==estimate[i]]\n",
    "    accuracy_viterbi[k] = len(same_indices)/len(hmm_test)\n",
    "    # posterior\n",
    "    probabilities = posterior_decoding(casino_params, observes)\n",
    "    estimate = [probabilities[i].argmax() for i in range(len(observes))]\n",
    "    same_indices = [i for i in range(len(hmm_test)) if hidden_states[i]==estimate[i]]\n",
    "    accuracy_posterior[k] = len(same_indices)/len(hmm_test)\n",
    "    \n",
    "print(\"Viterbi:\\tAccuracy median %f\\tAccuracy variance %f\"%(np.median(accuracy_viterbi), np.var(accuracy_viterbi)))\n",
    "print(\"Posterior:\\tAccuracy median %f\\tAccuracy variance %f\"%(np.median(accuracy_posterior), np.var(accuracy_posterior)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}