asr_lab3_4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ASR Labs 3 and 4 &ndash; Viterbi decoding\n",
    "\n",
    "In these labs we'll use the experience &ndash; and code &ndash; you have developed in previous labs to develop your own Viterbi decoder.  You'll want to refer back to [Lecture 5](http://www.inf.ed.ac.uk/teaching/courses/asr/2020-21/asr05-hmm-algorithms.pdf).  Remember that the Viterbi algorithm is used to find the joint probability of an observation sequence $X$ and the **best** path single path $Q$, allowing this best path to be efficiently discovered. \n",
    "\n",
    "Your decoder will find the best path through the WFST representations HMMs that you developed in Labs 1 and 2.  In the case that your WFST is just a linear chain, the best path will provide an *alignment* of the observation sequence to the HMM states.  If your WFST is a set of word or phone loops (or any other structure) then the best path will allow you to recover the most likely transcription for the observed acoustic features, subject to the constraints of your grammar and vocabulary.\n",
    "\n",
    "Observation probabilities $b_j(t)$ will be supplied for you &ndash; you do not need to use the observations $(x_1, \\dotsc x_T)$ directly.  These are supplied via the `observation_model` module.  You can use this as follows:\n",
    "\n",
    "```python\n",
    "import observation_model\n",
    "\n",
    "my_om = observation_model.ObservationModel()\n",
    "\n",
    "my_om.load_audio('filename.wav')  \n",
    "\n",
    "# or use dummy audio for debugging\n",
    "my_om.load_dummy_audio()  # will generate dummy observations as seen in Lab 2, \n",
    "                          # useful for testing\n",
    "    \n",
    "my_om.observation_length()  # returns the sequence length, T  \n",
    "\n",
    "my_om.log_observation_probability(hmm_label, t)  # returns log b_j(t) given HMM label in string form\n",
    "                                                 # raises IndexError if t > T\n",
    "                                                 # raises KeyError if hmm_label is not known\n",
    "```\n",
    "\n",
    "\n",
    "It's easiest to write your decoder as a Python class, and we will supply a template below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "partner = \"s2123998\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "import openfst_python as fst\n",
    "import numpy as np\n",
    "\n",
    "def parse_lexicon(lex_file):\n",
    "    \"\"\"\n",
    "    Parse the lexicon file and return it in dictionary form.\n",
    "    \n",
    "    Args:\n",
    "        lex_file (str): filename of lexicon file with structure '<word> <phone1> <phone2>...'\n",
    "                        eg. peppers p eh p er z\n",
    "\n",
    "    Returns:\n",
    "        lex (dict): dictionary mapping words to list of phones\n",
    "    \"\"\"\n",
    "    \n",
    "    lex = {}  # create a dictionary for the lexicon entries (this could be a problem with larger lexica)\n",
    "    with open(lex_file, 'r') as f:\n",
    "        for line in f:\n",
    "            line = line.split()  # split at each space\n",
    "            lex[line[0]] = line[1:]  # first field the word, the rest is the phones\n",
    "    return lex\n",
    "\n",
    "def generate_symbol_tables(lexicon, n=3):\n",
    "    '''\n",
    "    Return word, phone and state symbol tables based on the supplied lexicon\n",
    "        \n",
    "    Args:\n",
    "        lexicon (dict): lexicon to use, created from the parse_lexicon() function\n",
    "        n (int): number of states for each phone HMM\n",
    "        \n",
    "    Returns:\n",
    "        word_table (fst.SymbolTable): table of words\n",
    "        phone_table (fst.SymbolTable): table of phones\n",
    "        state_table (fst.SymbolTable): table of HMM phone-state IDs\n",
    "    '''\n",
    "    \n",
    "    state_table = fst.SymbolTable()\n",
    "    phone_table = fst.SymbolTable()\n",
    "    word_table = fst.SymbolTable()\n",
    "    \n",
    "    # add empty <eps> symbol to all tables\n",
    "    state_table.add_symbol('<eps>')\n",
    "    phone_table.add_symbol('<eps>')\n",
    "    word_table.add_symbol('<eps>')\n",
    "    \n",
    "    for word, phones  in lexicon.items():\n",
    "        \n",
    "        word_table.add_symbol(word)\n",
    "        \n",
    "        for p in phones: # for each phone\n",
    "            \n",
    "            phone_table.add_symbol(p)\n",
    "            for i in range(1,n+1): # for each state 1 to n\n",
    "                state_table.add_symbol('{}_{}'.format(p, i))\n",
    "            \n",
    "    return word_table, phone_table, state_table\n",
    "\n",
    "\n",
    "# call these two functions\n",
    "lex = parse_lexicon('lexicon.txt')\n",
    "word_table, phone_table, state_table = generate_symbol_tables(lex)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def generate_phone_wfst(f, start_state, phone, n):\n",
    "    \"\"\"\n",
    "    Generate a WFST representating an n-state left-to-right phone HMM\n",
    "    \n",
    "    Args:\n",
    "        f (fst.Fst()): an FST object, assumed to exist already\n",
    "        start_state (int): the index of the first state, assmed to exist already\n",
    "        phone (str): the phone label \n",
    "        n (int): number of states for each phone HMM\n",
    "        \n",
    "    Returns:\n",
    "        the final state of the FST\n",
    "    \"\"\"\n",
    "    \n",
    "    current_state = start_state\n",
    "    \n",
    "    for i in range(1, n+1):\n",
    "        \n",
    "        in_label = state_table.find('{}_{}'.format(phone, i))\n",
    "        \n",
    "        sl_weight = fst.Weight('log', -math.log(0.1))  # weight for self-loop\n",
    "        # self-loop back to current state\n",
    "        f.add_arc(current_state, fst.Arc(in_label, 0, sl_weight, current_state))\n",
    "        \n",
    "        # transition to next state\n",
    "        \n",
    "        # we want to output the phone label on the final state\n",
    "        # note: if outputting words instead this code should be modified\n",
    "        if i == n:\n",
    "            out_label = phone_table.find(phone)\n",
    "        else:\n",
    "            out_label = 0   # output empty <eps> label\n",
    "            \n",
    "        next_state = f.add_state()\n",
    "        next_weight = fst.Weight('log', -math.log(0.9)) # weight to next state\n",
    "        f.add_arc(current_state, fst.Arc(in_label, out_label, next_weight, next_state))    \n",
    "       \n",
    "        current_state = next_state\n",
    "        \n",
    "    return current_state\n",
    "\n",
    "def generate_word_wfst(word):\n",
    "    \"\"\" Generate a WFST for any word in the lexicon, composed of 3-state phone WFSTs.\n",
    "        This will currently output phone labels.  \n",
    "    \n",
    "    Args:\n",
    "        word (str): the word to generate\n",
    "        \n",
    "    Returns:\n",
    "        the constructed WFST\n",
    "    \n",
    "    \"\"\"\n",
    "    f = fst.Fst('log')\n",
    "    \n",
    "    # create the start state\n",
    "    start_state = f.add_state()\n",
    "    f.set_start(start_state)\n",
    "    \n",
    "    current_state = start_state\n",
    "    \n",
    "    # iterate over all the phones in the word\n",
    "    for phone in lex[word]:   # will raise an exception if word is not in the lexicon\n",
    "        \n",
    "        current_state = generate_phone_wfst(f, current_state, phone, 3)\n",
    "    \n",
    "        # note: new current_state is now set to the final state of the previous phone WFST\n",
    "        \n",
    "    f.set_final(current_state)\n",
    "    f.set_input_symbols(state_table)\n",
    "    f.set_output_symbols(phone_table)\n",
    "    \n",
    "    return f\n",
    "  \n",
    "# example\n",
    "f = generate_word_wfst('peppers')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABqQAAAA3CAIAAADlkLy4AAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOydZ0AT2ff3JwldKVIVkCJiWSxrwYYorJ2EDqKiiL0gYkcsu+gqYK/AYi9YAWmChSKoWFBWEUSKgjTpHUISkszzYv6/PNmETCaNFObzSpNz75x8uX3uPRcDgiCAgoKCgoKCgoKCgoKCgoKCgoKCIv1gxe0ACgoKCgoKCgoKCgoKCgoKCgoKinBAF/tQUFBQUFBQUFBQUFBQUFBQUFBkBHSxDwUFBQUFBQUFBQWlf1FdXW1hYTF8+PDi4mLh5kwmk+Pj43vN/+vXr8bGxiYmJt+/f6fT6Z6enoMHD379+jW7cXd397p16zQ1NUkkknDd6zOYf6xwc+akMMsTpUvhESNGLFmyRBQ5k8nksLAw5KURYKsdkqYVciSqmrP7I73CorVbKkC02CfepgeQxb8l2u4IHdG1OIBU9ZFDhw4NCQnh9C06jBAcGIVFWq/7SVMJI2/fjyoAKS/AIqrv/XYSK7oKDvSnJpQFkdZr2Wg20YosOtLS0tLT07du3RocHMzyVWdnZ3h4uIODA5J8PDw8zp8/397ezvgkMzNz1qxZveZvZmZWXFy8adOmzMzM+Ph4RUXFrKysgIAAdmMPDw91dfW8vDwlJSWBf6t4YP6xLF+JSGGWJ0qXwm/evOn1cwG1AgAgMzOzp6cHeWkE2GqHpGmFHImq5uz+SK+wnGo3T6oC/aZ2iwtEi31CaXoADn9Lrk0PIIt/S7TdETqiG08AUtVHXr9+HeZbdBghODAKc6rXaFOJHBh5+37OAEh5ARZRfe+3k1jRddxAf2pCWRDpbEE2mk20IvNHU1OTlZXVsGHDKisrDQwMnJ2dhw8fznItoZeXl56e3rp16+rr66urqzdt2gQAQHFx8bZt21xcXHA43J07d5jzUVJS2rp1q56e3u7du9++fTtkyJBLly4BAHD16lV5efkFCxZs3ry5oKAAAIDm5uZBgwYx5894qJKSEolE6urqcnZ2TktLW7x4sZmZGQ6Hc3BwYDZOT0/Pzc29detWREREnwqHGCQKM//YvlGY+YkAAEigwnv27MFgMBgMZtasWTBmkFwsWgFMsmMwGIZWAAAw5GLXCgCA5uZmPz8/5KWxtbWVWVvJL40AZ2EFqeYArLB8VHOpExamuLLXbpji2k9qtyQC/g9FRUVfX19dXd1du3aB/6WhocHDwwP6d1VV1caNG4uKivz8/ObNmxcREdHe3g6CYGNj44wZM0xNTSsqKpizevPmzeDBgyMiIkAQ7OjoCAsLmzZt2qZNm75+/QqC4L1796BsiUSinZ0d80NbW1sPHjzY1NQEgqCPj8+TJ09AEJw9e3ZLSwuzfVpampGRkY6OzsGDB0EJA0ZSCOhXQJKCIMiiqoCSgrCqskvKsJdeSZnLDNeCCgAAc1YMVdklBZEVVDFKylxUUlJSdu7caWFh8ddff7Fb8lSXZUkiAUGuMPjfeo02lUhALi/jx6LNJgR8k8io7xIil1SURoQdN8ihCeUkKdi/m1CEfXd/azaZf5e+vr6Tk5OZmRmdTmcxQysyfxw6dCgrK+v69esBAQEpKSl+fn6cLEtLS0NCQhj/tbCwOH36NOMPwZxPVFTU3r176+rqzMzMNm/eXFxczJJVeHi4qakpnU5/8OABp/xBEPT09MRisZGRkUuWLMnOzgZBkEAgFBUVMRv7+/tHRERUVFQYGBhAFUHSQKgw48cyPhG1wsxPlECF4+PjaTRaUFAQ5A8zzKN0CBatQCa5zM3NGVo1NTUdOHCARS5IKxAEmeXiqTQy7CW/NIKwwoL8VnMYYQWp5qD0CAuvKkvthimu/aR2SyD/f7GPpWQzG/HU9PD3t4SpIaDU/i1hJIWQqOEFKA3tDrykXMcToIz2kSzjLS8vr+rq6unTp7NbosMI/kCuMPhfBdCmEgnI5e37OQMo2QUYvknkWt/RSSw7Iuq4+3kTylPf3X+aTYRrJWhF5g9PT09oW4Ojo2NKSsqJEyc4WZ46dYpEIjH+W1RU5OvrO2XKlPPnz7e1tTHnExUVdfnyZRAEvby8bt++PXz48MOHD4Mg2NnZef369YULF27atCk/P//du3elpaWc8gdBkEKhJCUl4fH4LVu2pKSkgCA4ffr0hoYGZuNt27Y9f/4cBMHFixdXVFQIVxyhgFBhxo9lfCJqhZmfKJkKJyYmxsbGsn/OXtlZtAKZZAcAgKFVUVHR+/fvIblYtAJBkFkunkojw17ySyMEJ2FBfqs5jLCCVHNQqoSFUZWldsMU1/5TuyWN/xzjNTMz09XVtbKyamxshN8P+OjRo7KysmnTpl24cAE6FFBcXGxlZbVq1Spo+yUjq+nTp9vZ2f39998AAHR1dd24cWPRokVfvnxJTEzMzs62tLSEMoyNjd22bRvzI65fv56YmHjv3j0AALS1tdva2gAAaGpq0tTUZLYnk8mmpqZDhw61srJqbW1Ftp2x74CXlOVXs6gqoKTs+QNMqrJLyrCXXkmZywwETEEtKSlhzsre3h5SlUVSAAAQFlQxSspSVEaPHq2vr6+iooIkbT+RSEB4UphZAbSpRAJyebnW8X7YbCLvuAFxyyUVpREQTcfdz5tQgJe+u/80myy/y9DQEGFCtCIjwdTU9Pnz5y0tLWvXroUxS0lJ8fb2Li8vJ5PJ0CcjRow4f/58WloaDodbsWIFSz55eXlNTU2FhYVubm55eXlv376trKz09vZua2t78OABFCayrKwMOh7Ra/4AAMjLyxsYGIwZM2bu3LnR0dHFxcVUKlVbW5vZePLkydeuXauvr29oaNDX1xeZTvyDUGHGj2V8ImqFmZ8ogQqXlZV9+fLFyckpOTmZqzGLVgCT7CtXrmRoZWRkNGXKFEiuOXPmMGsFPRGSi6fSyGw/duxYCS+NAKywfFdzGGH5ruaAVAkLX1xZajdMce0ntVsSYSz7RUVFbd26tbGxccqUKd3d3cwrguzvGSA6OjpCQ0MdHBxAEDxw4AD0t0xMTGTJqru7e9GiRRUVFW5ubmfPnoUWekGmEwTPnz9vamoqKipiWRH//Pmzv78/CIJxcXEbNmwoKiqytLRksYfeOtbV1dna2lKpVGEuhAoMjKQg51/NUFUQSWHyh1RlkZTZ/urVq1IqKchUZpjptaCuXLmSJStI1alTpzJLCiIuqGKUlLmopKSkBAcHgyA4Z84cdkue6rIsSSQgyBXuVQG0qYQHubwgtzre35pN+CaRa33vS7mkpTSKruMG+3ETykff3R+aTZamj9PGKLQi80dNTc2ECRNMTU2/f/+up6c3cuRI9mNoBw4cgGZDI0eOrKqq2rJlC3w+UVFRdnZ2GhoaV65cefnypY2NzezZsykUCksShrzM+WdlZR0/fhwEwaCgIBUVFQKB0NjYSKPRHB0ddXV1MzIymI1BEKTRaAQCQVdXNyEhQfjqCAMkCjP/2L5RmPmJIAhKoMLz5s2D3LC1tWX5auTIkQAAhIeHgyDIVa6QkBCGViAIcpWL5eezy8WsFYu95JdGkLOwglRzkLOw/FVzUNqEhSmuaO2WCv6z2Mdcspnhqenh9W/JtekBpfZvCSOpRA0vQOlpd2AkZSkz/aqPZC4qxsbGAABAu6bPnDnDYokOI/gDocIirdcy3FQilLfvRxWgxDebME0iyFTfxS6XFJVGhB03yKEJ5WPCAPaDJhRh393fmk0kayUgWpElCcYxNBQRgSqMHFQrEYEKKyJQYcUCBvzfrUnR0dGtra3wO7ERIsSspBpUUqGDStqXoBKJGlRhkdIf5EWbRKGDSioKUFX7DFQfGBobG3V0dBj/9fHxuXjxIvLkJBJJWVkZAICGhgboLB4KC6jCPCGIXP1NK+SghVBEoMJKKf+32Mf8BwAAAP1bCg4qqdBh0QHtI7mCSiRS0G5PpKDycgXtZYQO2suIArSgIgdt9ySfX79+JScnZ2dnt7W1ycnJ9fT04HA4AAB0dXWtra3nzJkzaNAgcfso9bx58+bx48c/fvzA4XA0Gg2Hw9Hp9EmTJhEIhNGjR4vbO0mkpKTk0aNHX758AUGQSqVisVgAAMzMzBYsWDBr1ixxeyetVFdXJyUlZWdnd3R0yMnJUalUHA5HpVKHDx8+c+bMOXPmKCoqittHqYRKpT59+jQ5ObmpqQmq3QAA6OjozJ8/387ODmpRUYTF/9/ZJ3R+/vyZlJT04cMHMpkMgiAGgwEAAARBXV3dWbNmzZ07V0NDQ0SPllWIROLz589fvHhRW1srJycHXcSGwWDk5OSmTp06d+7cUaNGidtH6YPR4rS0tNDpdBwOB4KgkpLS7NmzXV1dVVVVxe2gOKHRaM+fP09KSoLuT8RgMBgMRkFBwdra2s3NTV1dXdwOSjdlZWXJyck5OTldXV1QwYM+HzhwIDSMQB6vHaVX8vLy4uLiCgoKsFgs1GACAKCurj579mx7e/t+Xrs58f79+8TExJKSEkYvAwCAjo6OjY0NgUBAh7a8AoLgu3fvnj59WlJSgsFgoNkCAAA0Gs3CwmL27NnW1tbQxAyFJz5//hwTE1NcXAw1niAIYrHYCRMmuLi4mJmZids7MZOXlxcdHV1YWMgszvjx411cXMzNzcXtnbRSVFR09OhRZWVlR0dHKysrliHQr1+/Xr169fjxY3V19b1796LdN3+8ffs2JCRk+vTprq6uzGWVSqXm5OQ8fPiwtrY2KCgIiviBAgBAZWXl/v37tbS0li5dOnHiRDk5OcZXhYWFiYmJr1+/3rVrl7W1tRidlDp+/Phx+PBhJSUlZ2dnKysrluFicXHxy5cvHz9+bGlpuWvXLnRcxBPPnj27ePGig4ODk5MT88utX79+PX78OCEhwcfHZ9GiRWL0UMYQyWJfXl5eUFCQjo6Oo6PjjBkzoNeMDGpqal69ehUfH6+npxcQEMD8Z0bhRFdXV3BwcFFREYFAmDdvHsu1MiQSKTs7OykpqaysbMuWLeg7HOQ8f/783LlzDg4Ozs7Ourq6jM87OjrS09Nv3749ffr0HTt2QEvV/Y0XL16cPHmSQCC4uLjo6ekxPu/s7MzIyLh169akSZN2796NTlP5ICcn5/jx44MHD3ZwcJg+fTrLzbOtra2vX79OTk7u6urauXPnuHHjxOWn9PL169fAwEBzc3N3d/fx48czl9KmpqaMjIz79++PHTs2ICBAXl5ejH5KFNAsy9LS0sXFZfTo0cztXk1NTVpaWnR09B9//OHr69s/m0Q+uHv37r1796ytrfF4/KhRo5jfV9Pp9IKCgrS0tLS0NDwev27dOrQtRUh5ebm/v7+ZmdmyZcuguyYhqFTqx48fb9++3dPTc+zYsf65waqqqsrf33/o0KGenp5jx45lfE6j0XJycm7fvk0kEo8fP66lpSVGJ6WRsLCwDx8+hISEMA+HeuX79+979uxZuXKlo6Nj3/gmM0REROTm5p48eZJlUMRMTU3Nli1bfH19bWxs+tA1CQXqtS9cuGBkZMTJhkQi+fv7Dx8+3NfXty99k15u37799OnTEydOcL3FNSUl5dy5c2fOnEFfoiDk7Nmz1dXVQUFBnMbeVCp1//79gwcP3r59ex/7JrMIPQrg2bNn161bB4U9hufr16/29vZQnGkUGPLz8+3s7F69esXVsrOzc8+ePfv376fRaH3gmLRz7ty5nTt3QjtPOXHnzh0vL6+enp4+80pCCA8P37p1K8vdfCw8fPhw2bJl7JG/UWCg0+mBgYFbt25tbW3lalxXV7d27dpz5871gWOyxIMHD5YvX861G3r69KmDg0NDQ0PfeCXhnD9/3sfHp7OzE97szp07Hh4eXV1dfeOV9NLR0eHp6RkWFsa1O6bT6Xfv3nVzc6urq+sb36SaDx8+ODg4VFVVwdhAo6afP3/2mVcSwufPnwkEQnl5OYxNYWGhnZ3d9+/f+8wrGeDkyZPQ1SUIoVKpW7dujY2NFZ1LskdsbOyuXbuQWFIoFDc3t/z8fFG7JOGUlJQ4OzsTiUQkxgcPHmS+gxuFEzdv3gwICICONSChoaHBzs6usrJSpF7JBrGxsQEBAUgsDx48GBUVJWp/+glCXuw7ceLEqVOnkNuTyeTly5e/e/dOuG7IEsXFxQ4ODkjWBRg8ePAAYX/Zn0lMTNyzZw8Sy+Tk5J07d4raH4ni2bNn27ZtQ2KZmprq6+sran9kid27d9+8eZOnJBcuXAgJCRGRP7JHamrqxo0bEQ7USkpKHBwcuru7Re2VhBMZGfnXX38hNH7//v2yZcvQV0owkMlkZ2fnnJwc5El+/Phhb2/f0dEhOq9kgJqaGnt7eyRrzTU1NXg8HuE0WDZoaGjA4/FIilB9fT1CSxQQBLOysvgY59BoNA8Pj9LSUlG4JHu0trY6Ojoi71aam5udnJyQr8jIJK6urshfEdHpdBcXFyR7cfozxcXFy5cv57VcVVVVoaWRK52dnfb29lQqFYkxnU53cHBAOymhgAsMDBTWJsGcnJwXL14cOXIEeRIcDrdo0aItW7Y4OjoqKCgIyxOZgU6nr1q16vr16zydRrGwsPj3338bGxvREH6cgPa0X758GcmxKXNz82fPnunp6Q0ePLgPfBM7PT09O3bsuHr1KpIIqcOGDXv58qWampqBgUEf+CbtJCQktLa28nqSYsqUKfHx8QoKCiYmJqLxS3bo7u7euXPnjRs3EMb31dTUHDRoUFxc3MyZM0Xtm8TS0NBw8uTJ0NBQhIdzDQwMWlpaCgoKxo8fL2rfpJTAwEBXV1eeCtWgQYNGjRoVEhJCIBBE55i0s3PnzsDAQCR98cCBA3V1dR89etR/IlX5+/sjDBU3YMCAoUOH3rlzBz0IiYStW7eGhobyOknBYDBTpkw5dOiQk5OTiByTJa5evbpo0aLhw4cjtFdWVq6tre3s7Oy3ATpzc3ObmpqQ9xcYDMbMzCwyMhKt9TDs2LHj+PHjvAZ0VlNTa25urqqqQm+PgSEqKmrixIljxoxBYozBYDQ1NdPT06dMmSJqx2QeYQaIOXHixNGjR3lNpaysvGPHjvPnzwvRE5khKirK2dmZj9AqO3fuvHr1KnS7DQo7T548cXd3R37dz969ey9cuCBSlySHlJQUR0dH5IHM9uzZExoaKlKXZAMQBC9duuTv789H2sDAwHPnzgndJdnj1q1ba9eu5SkM38KFC7Ozs0kkkui8knAuXLiwd+9ensLwrV69+uHDh6DILviSahobG0tLS+fPn89rwgkTJsjLy3/79k0UXskA7e3tJBJpxIgRCO3nzp374cMHKpUqUq8khO7u7ubmZoTzKAAArK2t8/PzKRSKSL2SAXJzc8eOHTtw4EA+0hoZGWGx2JaWFqF7JXu8ePHijz/+4CmJp6fnw4cPReSP5HP//v2lS5fylGTy5Mm5ubki8kcGaGhokJeX529jx6pVq6KiooTukiyRlJSEx+OR28+dOzczM1N0/vQfhLbYV11draWlxd8Fu3/88cfHjx+F5YksERsbu2zZMj4SysvLz5w5E1WVE4mJiQ4ODsjtdXR0Ojs7+8mcIT4+nqcX0YMGDerp6SGTyaJzSTZ4/fq1jY0Nf9dBDBgwYOzYsV+/fhW6VzJGWloaHxujHBwcnj17Jgp/pIL8/PzJkyfzlASLxU6bNi0nJ0dELkk19+/fX7VqFX9pN2zYcO/ePeH6IzMkJSXxWrttbW1fv34tIn8kiqdPny5cuJCnJPPnz3/x4oWI/JEZUlJSBNlsu3DhwpSUFCH6I5MQiURVVVXkL+ABAABBUEdHp729XXReSTilpaU8nd+CXs4ZGhrW19eLzCnpJiMjg4+3dBDKysogCNJoNOG6JEvQ6XSYi3dY+Pr1a1tbm7y8PCqp4Ahtse/169e2trZ8Jzc0NPz165ewnJENoMP/SkpK/CWfN29eWlqacF2SGbq7u5EvTL98+TI7O9vU1LS0tFSkXkkI7e3tyO/Ifv369fv374cPH15UVCRSr2SAzMzMefPm8Z18zpw56MSMKzgcjqcJA4S1tfX79+9F4Y/k09bWxnwROXL6s2jwvH//nu9T4RYWFujOPk68e/du9uzZCI1LSkpSU1PNzMz+/fdfkXolIfAkzo8fP1JSUkxNTfuJOILw7du33377jY+EFRUVeXl5o0aNKikpEbpXMkZNTQ2voWBMTU3Hjx+fl5fn5+d3//79uro6EfkmsfB6ezsejx89enRGRoavr++lS5fQQTs73759Y77eHTmhoaF37tzR1tauqKgQuleyQVdXF0/7oz09PceNG/fy5csJEybMmjULnf4IgpywMvr58yd/y+EJCQlqampDhw6tqqriesV1v6KhoYG/vcTl5eUZGRlWVlaVlZVC90o24OnAWnh4eHJysoqKSkxMjKam5pYtW1asWCE638QOT+JcvXo1JiZmwIABd+/e1dLS2rBhw+rVq0Xnm1RTVlZmbm7OR8Kmpqbu7u7hw4c/fvxY6F7JEhQKhde3IxkZGWvXrh01apSqqurXr19Hjx7N6wBa2mlubtbW1uYpyd27dwMDA0eOHKmvr19SUsJfqZZhenp6+HtLFxoaOn/+fD5Wq/sJzc3NyKOaVFRUODk5qamp4XC4yMjIqVOnhoeHi9Q98dLQ0ID8LV1NTY2rq+vAgQNxOFxUVNSkSZMuXbrEU9fffyASifyd4f327ZuHh4empubAgQMVFBTs7e3RONqcqK2t1dPT4ymJjo4OdHrp27dvt2/fzsvLE41rEgqdTud1rGJubv7kyRMAAPLz8x8/fnzr1q2RI0eKxjtppaOjQ11dnY+EsbGxmZmZWlpaz58/HzJkyKFDh9DAiCzU19fz9F75999/v3nzJgAA1dXVQ4cO5e+NiwzQ2tr68uXLL1++/Pjxg0wmgyCorKw8fPjw8ePHz5o1C2FwSaEt9pHJZEVFRT4SRkZGJiUlqaurR0dHr127lr8QdTIJhULh79KSlpYWX19fNTU1OTm5Xbt2eXt7Iw/j0h+g0+k8DWqtrKwePHgAHRbQ09MbN26cyFyTCHiKw2VtbX3r1q2Ojg4AANrb29GSBgPfNbqoqIhAIBgYGOjp6ZWXlxsbGwvdN9mgra1NTU2NpyRDhw7t7OxMSkrCYDCpqak6Ojpfv37tV5PetrY2XmNRGxoa1tbWlpSU4HC4mJiY33//PTU1VUTu9StCQkICAwMVFBTc3NzWr19va2vL36l/WYVKpSJfCR0zZoy6ujp0ZKS+vp7vw1nSQnd3t7KyMkJjCwsLNTW16upqAABqa2utra37VaPHEzgcjkqlysnxPF3S1dWVk5MrKysDAGDv3r2vXr1KSEgQgYOywMCBAzs7O3lKMnr0aGixD4vFOjs797c74rBYLK9h2W1tbf/55x8oTKe5ubmLi4toXJNiFBQU+AtJZGZmlpaWVldXV1dXx1Ng2f6DqqoqNFVEyKJFi+7evdvT04PBYObMmcPrywAZ4Pnz59evX1dSUrK1tXVycjI3N4fW2UgkUnFxcU5OzpYtW6hU6oYNG2bNmgWfldC2MAwaNKi5uZmPhBoaGkQisaam5v379zt27Hj16pWwXJJ2NDQ0+AvrO3DgwAEDBlRXV5eXl58+fZrXqz9lHl77yMmTJzNe9VhaWqIXUDIzYcIExlXRkyZNQm9NgoHvRhKa8Ofn57948WLKlCmOjo7Cdk1GUFVVbW1t5SnJ0KFDoRVYEATpdPqZM2f626RXQ0OD12JpbGwMbXWh0WhycnKnTp0SjWvSCt/3lqipqTU2Nv769SsmJsbT07OtrU24jkk7ysrKRCIRobGenh5j8cvQ0PDPP/8UmV8SgYaGBvLWb9CgQYzYSUZGRkFBQSLzS+rR19evqqqCMeBUT/X09Bh7IIyNjaFdKii9YmhoCC09swOCYK8K29jYQFvbhg4deuLECdH6J220t7ezT3N+//13aDONpqbm+fPn+9tQBwmGhoYwR+LodDqn5SrGvjNlZeX169ej5xTZ0dLS4jTUJJPJ7FfkTZs2DTp0YmhoeOzYMZH7J0nk5OS4u7t/+vTpypUr169f9/LyGjNmjKKiIolEIpPJSkpK48aNW7Vq1c2bN8PDwzMyMpYuXZqfnw+TodB29o0cOfLbt2/scWqam5sLCgp+/frV2tqqp6c3ZMiQcePGMZ9wYYROk5OTs7OzQ6+oZ6Cqqgo/1m9ra1NSUmLfUMkc5tbU1DQ6OlpULsocJBKpp6eHZZ+LhYXFwIEDW1tbDQ0Nw8LCxOWb2CGRSBQKhWXn1OjRo1VUVJqamgwMDP755x9x+SYVjBw5sqCggNM+dhqN1t7erq6uzn40Q05ODhqW0el0eXn5I0eOiNxX6URJSQnmrWx7e7uioiJLg6mgoMBYEZg6deqCBQtE66LkYWxsDDO6bW1tVVFRYdmRamhoCC1A43A4AoGAvv9gAXpRx3gLwgKRSKTT6b0eDGTsHtLS0rpx4wavx6tlHiMjo8rKSvajZ1DjqaGhwTJ9hXqrAQMG7NmzB3lccCnFyMiooqKCPRgxp54FEkdFRcXPz2/AgAF956i0MWnSpA8fPpiYmDA++fXrV0JCQmZmJgiCWCwWGqtjMBgcDjd37lx7e3toSUVHRwcaig8aNCgoKIhTg4ACAICWlhbzrRGVlZXx8fHQ1To4HG7gwIGQwnJycvPnz8fj8ZqampaWlpqamu3t7Z6enpqamuLzXZww9py2tbUlJyc/ffqUQqGAIKimpkYkEqFLBadNm+bk5GRiYmJsbAyNdsaMGcN1K1D/ZOLEia9evWI+gVtaWhoXF/f+/Xuo+CkrK3d0dGAwGEVFxUWLFi1atAhqSEeOHKmkpEQikczNzQMCAsT2AyQYDAbDGJ/T6fS3b9/Gx8eXlZXJy8vLy8tjMBgKhUKlUocOHero6GhlZWVsbDxgwAAsFmtnZ8dfXGkp5cSJE+Xl5REREZqamjQaLTMzMz4+vqqqSk5OTkFBAQTBnp6enp4eU7+bymIAACAASURBVFNTR0fH6dOn//nnn/X19fv27Zs8efLGjRt7zVNoi33Tp0+/e/fuunXroP9SKJTr168/ePDg1atXVCpVUVFRXV29qamJRqMNHDhwwYIFK1eutLe3BwCAcWh3+PDh169fF5Y/soGKikpnZydjVkCj0Z4/f56YmNjY2IjFYjU1Nbu6uqC/+rRp01xdXaHhiKqqKjTk1dbWDg0NRY9Fs6OqqtrU1AQpU1hY+OjRo3///VdBQUFBQUFJSam1tRUEwSFDhjg6OtrY2Kiqqg4YMACHw9nb2w8dOlTcvoscTU3N2tpaKF5kcXHxo0ePPn78yCKOnp6eo6Ojra2tkpLSwIEDsVjswoUL0eOl8MyePfvhw4fMw4jW1ta4uLjU1FQqlSovL6+qqtre3t7T06OoqDh//nwHBwdoGCEvLw/N0zQ0NAICAsaOHSuunyD5yMnJkUgk6H0SjUZLTU1NSEhoaGjAYrGDBg0iEolQgzllyhRXV9dhw4YBAABNdIcOHXr16lUxey8OMBhMT08PCIJQr0GhUJKTk5OSktrb2yHRoIvIaTSatbW1i4uLvr4+DoeDlk5MTEzOnTsn7l8gcVhbW7948YL5kFRubu6jR4++fv0KtaJycnLQzgsTExMnJ6cZM2ZAZtASqoKCgpOTEx6PF4/3Esy4ceM+fvwILfY1NzfHxsampaVBr0CgBRcqlaqkpLRw4UJ7e/uBAwcaGRl9+vTJ1NSUMTSVYcaOHZuTkwOFGWlpaYF6FhqNxtKzLFiwwMHBQVVVdcSIETk5OcbGxj4+PuL2XaKxtbXdv3+/u7s7AAA1NTVHjx4lk8mLFy9es2YNyyn77u7ulJSUbdu2GRoa+vv7a2hoKCgoYLFYKyurJUuWiMl9qUFRUbGjo6OtrQ16nbl48eKNGzeynJ4mEonPnj3z9fUdNmzY9u3bVVRU1NTU9u/fLyaXxcz06dOzsrImT558+vTp/Px8Nze3ixcvsmxWoNFoWVlZp06dIhKJBw8eVFNT09PTCw0NFZfPEs6kSZOOHz8O/fvHjx9///23urq6u7u7n58fSwSJ9vb2J0+erFmzZtKkSX5+fkZGRqqqqmpqateuXUOj7nJCW1u7vr4+JycnLCzMxsZmw4YNZmZmLDZlZWXx8fEnTpxYv369pqYmhULpPxvPyWSyr6/v1KlTd+/eDQBAfHz85cuX582bt3XrVua3TRA/fvyIjY09duyYj4/PggULrly5cv78eT8/v1OnTvUSdAIUHh4eHp2dnSAIxsTEmJmZDRgwYNWqVfHx8Y2NjZABjUarrKy8du2avb09FoudOXNmTk7O5cuXAQAYPHjwv//+K0RnZIOHDx/evHkT+ndsbCwejw8PD6+trWUxo9For1+/3rBhg4+Pz69fv+h0upGRkZKS0qZNm/rcZekgNjb2xo0bP3788Pb23rZt28ePH6G7j5mBDkE7ODikpaVNmTLF1NS0vb1dLN72McnJyREREWVlZWvWrNm6dWt2dja7OJWVlefOnbO3t3/27NmMGTOMjY1bWlrE4q10YW9vD4lJJBJDQkLc3Nzi4+OJRCKLWWdnZ0xMjLOz8+nTp0kkUmFhIRT9x83NTRxeSxORkZH3798HQTAhIQGPx4eGhv769YvFhkajvXnzZvPmzRs3bqyqqrK1tVVWVg4ODhaHvxLBsWPHXr16RaPRIiMj8Xj89evXGb02g56enrS0NG9v7127djU1NU2aNElDQ+PBgwdicVjCaWpq8vLygv6dn5+/dOnSgICAL1++sFt+//796NGjLi4u7969A0Fw6tSpAACMHTsWCsOMwkJHR8eKFSs6OjoOHz7s4eHx+PHj7u5uFpv29vaHDx86OjpevHgxODhYWVk5JSVFLN72MSQSydPTs6urKygoyN3dPSEhodeeJTo62snJ6ezZs+fOnVNWVo6PjxeLt9KFh4dHR0fHgwcPXF1dCwoKuNpnZ2fb29unpqaOHj3axMSkqampD5yUdm7durV161YPD4/i4mKuxllZWQQCQU9P7+zZs33gm2RSWlrq6elpZ2eXnp7O1bisrGzFihUTJkyws7PrA9+kl127dn379u3cuXPe3t7l5eVc7VNSUuzs7F6/fq2mpsbo91F65fbt27a2tsePH2fvuFkgk8lnzpwxNzdfv3593/gmCaxduzYtLQ0EwdbWVg8Pj7Nnz3IdCnZ3dx8/fnzlypXQ+ltiYuLWrVvZzYS52BcXFxcWFrZnzx4sFrthwwb2KRYznz59mjNnjpKS0vbt2xUVFY8cOSJET2QGMpns4OBAIpHWrFlz6tQprn/1goICJyen9PT0IUOGjB8/HtqvgcJOV1eXs7Ozh4fHz58/uVoePHhw7Nixx44d6xvfxA6JRHJ2dnZzcystLYW3JBKJhw8fHj9+fGBgYN/4Ju2cPXs2OTm5rKzMzs7u2bNnXO0TExPt7e3fvHmjpaVlYWHBPnlDYYFIJC5evHjDhg3Hjh0jkUjwxkVFRS4uLgsXLhwxYgSNRusbDyWQ6urq1atXL1269NKlS1QqFd44JyeHQCBMmTLF2tq6b9yTRtasWVNVVXX58mVvb++amhp449bW1u3btx85cmTWrFmDBw9GMtftt3h7e9vZ2b148QLejE6nx8TEzJw509LSsk/8kgjWrFmzaNEiJIubcXFx1tbW48eP7wOvZIAnT544ODgcPHiQ/cUnJygUio+Pj76+PvpGBCE7d+7cv38/coVJJJKlpeXVq1dF6pUk8/DhQ2dn566uLuRJli1btnr16v482uHKly9fFi5cGBoaijxJR0eHp6envr5+R0eH6ByTdiorK/F4fE5ODvIkN27cWLBgQVlZmcickiCCg4MvXboEgmBpaSkej8/Ly0Oe9uPHjwQCoaqqCgTBs2fPhoWFsRgIc7GPRqONGDFCRUUlOjoaiT2dTj9w4AAGgxk2bBjy9p0T5ubmHh4eAmbCDolEiouLo9Foy5Yt09PTe/XqVa8PZTc4f/48AABz587Nz883MjIyNjYuKSkBQZBIJK5du3bQoEFcF7YhTp06NWfOnJcvXyJ0mEwmb9q0SVtb+8ePHzz8Tjaqqqp+++03MzOzoqIiQfLpFRKJ9OjRI2FJCvKu6pkzZ/bt24e8wzt58uTSpUsFXzxl95wdJOKAbPowxGF/Cq/ihIaG7tq1i+ucn8HZs2fd3d352IpiaGgIs52q72t0H8jb1dVla2trb2/Pdf7PoLKyEnovwnXtlQUYeUVXu8XeYBKJRBsbGyQvuiF6enr++OOPP//8k4cfCYIgrLxcqzlLWrGL1tTUNGPGDOSb67u6uiZMmHDt2jWE9jCIqJqD4u5lvn37tmDBguPHjyN3+Pbt28OGDbtw4QLyJL0i0tot3r4pOzvbycmpoaEBocM5OTmzZ8/mteVkB2HHzUctFqI4nz59cnBwqKurQ/ijvn37NmvWLCGuLIu3LoOi1Pbs2bP8VUw3N7eoqCg+ErIj9jE5yFlhwcfkwcHBfCzb0Wi0HTt2JCYm8pqwV5BUcwbiHUmCIJienr5x40Ze5840Gu3hw4f79+/nKVWvyOQQHQTB3bt3x8bG8uo2jUZzc3NjqTu8IsN1vKOjg0AgVFdX8+p2XV0dHo+H4kfxTR9MvQXU5+3bt9u2bQNBsKWlBY/H19fX8/obq6qqCAQCtL9v7dq1LGuFQruNFwCAkJCQnz9/PnnyxNXVFYk9BoP5+++/9+7dW1FRkZaWJuDT37x50+vnnZ2d4eHhDg4OSDLx8PA4f/58e3s745PMzMxZs2bFx8crKipmZWWxxN1kPJTdAIPBgCCYkpJiZmZWXFy8adOmzMxM6BHq6up5eXnMt5RwAgTBvLy8o0ePWltbI/EfAAAFBYWwsDBra2v468O4kpaWlp6evnXr1uDgYJavBJQUAIDMzMyWlhZhSQrwqOqTJ09qa2uPHj3KfhMCJ3bu3Llq1SroFL0gsHvODhJxADZ9GOKwP4UncZ49e/b9+/cTJ04gjzrh5+fn6+vLx6XP8DE6e63RPJU9gMca3Qfydnd3Kygo3L59GwqJiAToWpjffvuNJRQLV2DkFV3tFnuDuXnz5kOHDtna2iLxHwAAOTm5pKSkxsbGFy9eIEwCASMv12rOkla8olGp1DVr1ly6dGnChAlcjSFUVFRev36dkpKSl5eHMAknRNRxA+LuZd6/f//HH3/w1GUsX758zZo1LS0tyJP0CqfaLZTGU4x9U01NzeHDhyMjI5HfWzJx4sSYmBhfX18B7zVG2HHzUYuFJU5dXd3evXtv3bqFPJD5qFGjYmNjt2/fzusN5pwQb10GRKZtenp6aWnpli1bkPjPwoMHD+Li4r58+cJHWhbEPiYHOCssYGuZlJTU0tKyevVqJP4zg8Vijx8/fvv27aKiIl7TsoOkmjMQ70jy58+fFy5cOHfuHK836mKxWHd3dwUFhYcPH/KUkB2ZHKLfvn1bS0uLj5tCsVhsZGTkyZMnOV0qjQQZruPQyJyPS4p1dXVDQkI43TuBkD6Yeguoz+nTp6ETchs3bjx+/LiOjg6vv9HAwCAwMBCKwBsUFMQS6FBoi32fP38+ePBgaGgor1f8BAUFLV++3MvLq7OzE95SSUlp69atenp68APo6urqTZs2AQBQXFy8bds2FxcXHA53586dpqYmKyurYcOGVVZWsmT19u3bIUOGXLp06erVq/Ly8gsWLNi8eXNBQQEAAM3NzYMGDUpLS1u8eLGZmRkOh+t1YMRikJubGx0dDbU70AU90LnR9PT03NzcW7duRUREIBEnNDR0zpw5UBwfnrh79+7x48cFkdTLy0tPT2/dunXQDVmQqiySAgDAUBWDwSCUFACA5ubmf//9VyiSAgDAk6o1NTVXr149evQoV0sW5s2bp6+v/+jRI66WMKoye84oqCz2SMQB/qtPZmYmQxxAgCLX0NDwzz//MMLTIsfa2trCwuLevXtcLZmrIQAAjY2NY8aMCQwMhEnCteyx1GhG2QMAQCg1WljyAgDg6+sbFhamrq6OxJjBiBEj7ty5g2ThAKG8zLVblhrMy5cvW1tb89oHKSkpnTlz5vTp01xXBBDKy9JAQQrDpBWvaEeOHNmwYYOFhQUSYwYqKirh4eEBAQE9PT3wlgg7boBDTeevHAJi7WWKiopevHixZ88erpYs7Nu3r7Oz8+XLl1wtEfbdQm88xdU3AQDg6+sbHh7O672xWlpaR44c8ff352opeMctxp7Fz8/vn3/+4bVn0dTUPHbs2M6dO7laMhcYAwMDZ2fn4cOHgyDIyZ7vgsd3XQZEoy2VSj1z5syJEye4WvYKFou9ePHi4cOHuVpyVZjXMXlfKixIa0mhUP755x++Y/DjcLiLFy/CjyEhuCrMteNmTiXe+h4YGHju3DkFBQUkxuwcPHjwzp07XV1d8GZ79uzBYDAYDAZ+TCUzQ/T29vbY2FgknUWvKCoqQks2XC05CSurdfzJkyfjx4+fOHEiV8teGTNmzMyZM2NjY7lachK21x6cxViM+jx+/HjGjBnq6uoxMTE2Nja//fYbAlV6YdKkSb/99tuzZ890dHRGjRr1n5VNXjcKcmL+/Pl8R/BpaWnR0tI6dOgQvFlUVNTevXvr6urMzMzY4902NDSw7Ci2sLA4ffo0Y5PzoUOHsrKyrl+/HhAQwJLVgQMHWM4yhIeHm5qa0ul0KNzGkiVLsrOzQRAkEAjM22sZD2UxoNPplZWV27Ztgw5ge3p6Qqv+/v7+ERERFRUVBgYGXG976OjocHV1hbeB4eXLlwJKCoJgaWlpSEgI478skoJMqpqbmyORFARBSFVhSQqCIE+qbt68+evXr/A2nKBSqQQCget5VXhVmT1nt29sbEQiDvhffQoLC5nFAfktcn5+fp8+fUKkBRs0Gs3JyYlrlDTmapiSkuLl5VVdXT19+nR2S5YaDVP2WGr05s2b2Y8mIanRopY3IyPj6NGj8DYw7Nu3LysrC94Gubwgt9otdQ1mZ2eng4MD3xEh3r9/z/VsC3J52as5p7TiFa2urm7FihUIJWInOjqaa2ibvu+4QXH3Mp6enpWVlfA2nCASic7OzlyLMU99txAbT3H1TQkJCWfOnIG3gWHbtm1cT6kL3nGDPNZiYYnz7Nkz5pacV/bs2fP27Vt4G5bmy8/Pj92Ga10GeSl4vNZlUDTaRkREQDc+CcKpU6eSkpLgbZAoDPIyJu9jhfluLc+ePZuQkABvw5W///6ba+wOJApz7bihVOIdSWZnZ/v7+yNQBY4XL14cPnwY3iY+Pp5GowUFBfV6qlT2hugHDhzgOsbmyvbt2z9//gxvAy+sjNVxGo1GIBAQnqHmBIVCIRAIXENpwQjLXrWZjcXbBrq7u5NIJIS/ER4ikWhvb0+j0dra2jw9PRmfC2dnX25u7vPnzzm9maHT6Z6enoMHD379+nWvBhoaGv7+/mfOnKFQKPAPMjMz09XVtbKyamxs5OrVo0ePysrKpk2bduHChfb29uLiYisrq1WrVkHr3MxZ2dvb29nZ/f33311dXTdu3Fi0aNGXL18SExOzs7MtLS0BANDW1oY2fTQ1NWlqarI/i8UAg8EYGhqeOXPmw4cPAABcv349MTHx3r17ZDLZ1NR06NChVlZWXI9OREREwBwcGDFixJIlS2CSW1tb5+XldXd3wz8FXtLY2Nht27Yx/ssiKQAADFVLSkqQSAoAAKSqsCQFAAC5qrW1tRQKhdOqeXV1tYWFxfDhw4uLi3s1wOFwy5cvR7L7HUZVZs/Z7ZOTk5GIA/xXHy0tLWZxAL6KXFNTU2tr6++//97rt1+/fjU2NjYxMfn+/XuvBlgsduXKlZGRkTCPAJgKDFQNR48era+vr6KiAp8KgC17LDV6+vTpUNkDAEAoNVoo8gIAEBYW5ufnB2Nw/vz5K1eucPp2z549oaGh8I/gSV742i11Dea1a9d8fHw4nWq5cOECBoOZN28ep+RTpkz58eNHR0cHzCOQy8tezTmlFa9op0+fhnmP3d3dvW7dOk1NTRKJ1KuBi4vLs2fP6HQ6/FP6uOMGxNrLfPv2TV9f39DQsNdvU1NTMf+jtraW3UBZWfmPP/54/vw5V6GQ991CbDzF0jcBAHDjxg3obTw7XV1ds2bNGjRo0IMHDzgl379/PxRbBx4BO25ATD3LlStXtm7d2utXZDIZurSUse2Fnb179/Las3Aq3izwUfD4rsuAaLR9+vTp4sWLe/2KMQjnOsfZtGkTr0MjTgojH5P3scL8tZYAAGRmZtrb2/f6FfM0B37K4+fnBx8Whl2ZXhXm2nFDqcQ7krxy5cqOHTt6/YqhEtfpjI2NDbQmBfMgBweH5OTk0aNHjxgxAt4lQPqH6DQa7du3bzNmzOj1W4awXKdCu3btunz5MrxW8MLKWB3PysqytbXt9Sgre6XmNAmSl5dftGgR12A7MMKyV21mYzHqU1tbq6WlpaiomJaWRiAQ5OTk2G1YhIKZzigrK8+aNevt27dqampYLJYxoxHOYl9cXJyJicnMmTN7/ZZTEBNmvLy82tvbuQZKyMvLa2pqKiwsNDIy4urViBEjzp8/n5aWhsPhVqxYYWpq+vz585aWlrVr17JkNWXKlLy8vLdv386ZM6etre3BgwdhYWEWFhZlZWWmpqYAAMydOzc6Orq4uJhKpfYaL6ZXgy9fvmhoaAAAIC8vb2BgMGbMmMmTJ1+7dq2+vr6hoYHr2fW3b9/Onj2b07ecwqMwg8fjU1NT4W1gJE1JSYGuHieTydAnLJICAMBQdeXKlUgkBQAAUlVYkgIAgFzV+Ph4Nzc3Tt/ChEtg4ODg8OzZM5hHQMCoyuw5uz2NRkMiDtCbPgxxAL6KXGJioouLC6dvkYQ8wOPx6enpMI8AmAoMVA2RA1P2WGq0m5sbVPYqKyu9vb2FVaMFlLepqUlDQwPmDNr37985jR4g1NXV5eTk4FejkMvLtXZLXYOZkZExZ84cTt8yh3ThBNfajVxe9mrOKa14Rfvx4wfMAV6uMUcwGMzMmTOzs7Phn9LHHTcg1l4mJiYGZlKKxWKhO7WdnJw4jSYXL14cHx/PVSjkfbcQG0+x9E3V1dWGhoaKioq9fpucnOzk5FRcXAwTR0JbW5tEInFas2YgYMfd629nQeji1NfXa2lpKSsr9/rtjRs3CARCQkICzHumQYMGAQAAf6aPv46bj4LHd10GRKBteXm5qakppxdIMKE/WVBWVsbhcIL33TyNyftYYf5ay2/fvo0ePZrTt8zTHPgpj6qqKoVCga/gSBRG2HGLcSRJp9NbW1s5heZkqIRkOmNpaclYHeuVsrKyL1++ODk5JScnw5hBSPsQ/c2bN1ZWVjDfQv/gOhXS19eH7t+DeRaMsLJXx+Pi4qCTreywVGr4SZCTk1NCQgLMgwBYYdmrNrOxGPVhzFzi4+MdHR17tWERCn464+zsDA0gbWxsXr169X+fCrJdkMG0adM47TkHQdDHx+fJkycgCM6ePbulpYWTmaWl5Y4dO2CeEhUVZWdnp6GhceXKFfZvR44cCQBAeHh4VVXVli1b2A1qamomTJhgamr6/ft3lqxevnxpY2Mze/ZsCoXCnOTevXvQP2g0mqOjo66ubkZGRlZWFuOKPcZDmQ1AEExISHBzc3N0dGxqagoKClJRUSEQCI2NjdB2Vl1dXeZd60QicfPmzSz7fhsbGzds2ACjBvuJCXbq6+s3btwIYwAj6YEDB6ASMnLkSBAEuaoaEhKCRFLwf6qyKMZQlVdJwf9tEmZR9eTJkw8ePGA5Vbps2TKu50yJRKKdnR2MwZIlS+D32cKoyuw5Q1Jme4TisNgwi8PylF7FCQ4OjomJYTmPvGLFiq6uLpjf1draevDgwV6PjDHw9PSE363NXA2NjY0BAPD09AQAgP2IFt81mmvZAznUaGHJm5ube+DAgYqKCuZH37lzJyYmhpMsNBotMjLy3r17ly9fhlEvMjIS/powhPIy126pazBBEPTx8Xnz5g3zJ21tbatXr+Yky+fPn2fPnv38+XMY6UAQbGlpWbNmDYwBQnlZGihIYU5pp06d2geiVVZW7tq1i+WIR1FR0b59+zj92LS0NCMjIx0dnYMHD8JoUlBQcODAARgDhB03yKGX4aMcgn3Vyxw6dCg+Pp7l6e7u7lwP4ba2tq5duxbGwM3NDT4HhH230BtPUfdN6enpQUFBLFfKXrp06enTp5ykyM7O3rFjR3d396lTp2AUi4iIgG8BBO+4QR5rMa/iJCcnnzhxAiqQDG7evBkfH8/pR61bt+7Lly8gCE6bNg1m+H3jxg34K02ZC4yent7IkSPZz0lxrcugBBe8Dx8+/PXXX1VVVcyPvnz58rNnzzhpwhiEI5njXL16FbLhBFeFeR2T96XCSFrLp0+fHj9+nOUq7fPnz798+ZKTJszTHK5TnosXL0Kec4Krwlw7bkaqvhlJNjY2+vn5scQd+vfff4OCgpAoBnKbznz+/PnIkSMwijG2Dtna2rJ/K71D9KSkJPaG9M8//ywoKOAkBbOwXKdCgYGB+fn5nL4FOQsr7XX8yJEjcXFx7CMiTjowq4pkEgSTFQQnYXvtwZmN+0af+Pj4M2fOsJQcf3//srIy5EIhmc5AWeXm5jIiRwlnsU9XV5dxVJ4dTqegWfD29nZycoJ5SlRUFPxMGDlCzEpwOjo6dHR0tLW1R44c6efn9/37dxAEs7OzT548CZMKyWIfCIJLliyB+VZWJQVBcPr06crKysbGxvPnz09OTqZSqSA3NSBYwiWws3v3bvh4TLxK0ffSWVpaqqioGBsbL1y4MCUlhUajgQjEYQ95wM6BAwegAtzHSFTxS01NVVRU1NXVHTduXEhICDS6PXToEEywyOjo6I6ODq6LfR8+fIBvFkSERMkLgqC+vr66uvqwYcNWr14NDao+f/4MMwJmicUGA5ImQhopKCjQ0NDQ0tIaNWrU/v37oTltUlLSzZs3OSVBGHOESqUuX74c5tEy3Mv8/vvvKioqRkZGeDw+IyMDYUMKguD9+/fhY3h5e3tDGwA5IauqRkREYDAYPT29CRMmXLx4sa2tDQRBf3//nz9/wqRav369tbU1/HumjIwMaLDOCcnvuENCQrBYrJ6enqWl5eXLlzs6OkAQ/PPPP2HG1f7+/n/99VdjY6OFhQXMm863b9+ePXtWJE5zRqIKXkJCgpyc3ODBg8ePH3/q1CloPnbw4EH24GIMOIX+7NU4PT09IiJCNL5zRKIUPn36NFR6J0+eHBERAfUpu3btghlO87TYl5SUdOvWLeH6LF6qqqq0tLSgueHu3buhNjAuLo6xIsYOi0rw05m2tjb4PSVCRKKKYnBwMKMhvXLlCtSQwve5zMJynQrdvHkzOTlZ6G73ikQJO3HiRGhEtGjRorS0NK4jImZVkUyClixZwndUbkngzz//xGKx+vr606ZNu3btWmdnJwiCy5Yto9FodDp96dKlnBIyC4VkOgNp3tXVxdgD0cvZYF7p6elpbGwcMmQIJwMkp6ABANDX13/37l1UVBSnp0D7ILBYLAAAa9asYXy1YMEC5v8icZiRlaqqKvKEIoJEIikoKFRXVzc2NhYVFUVGRkIDDk7Honmis7MTiaSqqqodHR18qyppkgIA0N3d3d3dXV5eXl5e/vr1a11dXUtLy9bWVhAE4e+qZwmXwM6QIUMePHjA6TwarwVVLNJRKBQikQiJ8/btW21t7UmTJnENFHL9+vVly5aFhYVBDvfKkCFDYmNjoY1LyBGk7AGSV/zy8vJwOFx9fX19fX1+fv6ZM2eGDh06YMCA1atXc0qyd+9exgHzyZMnc4qcqK+v/+XLF041mhMyJi8AADgcrq2tra2trbS0NCEhQUtLa8yYMZyCrQAAwIjFtn79+nXr1sHk3N3dzau8gMAK9wGNjY0KCgr19fVNTU3BwcFXr16FyCp9KQAAF3FJREFUNinAHDhliTnC6U+Pw+Gampr6Zy/T09NDJBIrKioqKipev36tpaU1ZcoU+MN6EG/evDl58iSMgY6ODvQ34vRcWR0OQW9E6urq6urqtm/fHhwcbGxsjMPhYILA1NXVTZs2bfLkyTt27AgLC+Nkpq+vHxkZKawRplh0Kysro9PpkDifP38ODAw0MTFRVFTkFLAPAID169fPnTs3Ly8P5hw0AAD6+vqfPn3iqemTsW4lNzcXAIDa2tra2trdu3efOnXKyMhIQUEBJnA2AyRznCFDhty7dw86MY0QGVP427dvjNKbm5t7+PBhExMTLBbLOHopIEOGDImPj+cUcaJXJLzjbm5uVlBQqKmpaWxsPHHixM2bNwcPHjxs2DBvb2+EOcBPZ9TU1MrLy5HUehkriiwN6V9//QXdYoHD4ZAk5zoVGjJkSFJSUmdnJ9esZExY5hHRmzdvtLW1J06cyDiPDA+SSZCqquqNGzcGDhzINTfJrNrV1dV0Ov3Xr1+/fv3Kyck5ePCgubk5lUqlUqnt7e1aWlpIMkEynVFSUuru7lZRUSESidAnQljs4wp0wtnIyAjmFDQEyPmUu7y8PPPFCEguSUCYlaRBo9HIZHJbWxvX+DICwqKDqqoq37JIuKRUKpVGo7W2tpLJZDqdDtOgM8IlGBsbcxoZw68V8lpQxS5dT08PlUpta2ujUCg0Gg1GnF6DFrEALw4nBCl7gARoCAP0GopCoeBwOJgaXVJSAgDA/fv3Ozs7Oa30Af+L1MCrDzIsLwAAPT09PT09LS0tXC8jYg7pIlwEVLiPodPpPT09ZDK5s7MTpkxCMUfGjx+PJPYfJ/pPLwM1pK2trdBLWmjNqFfodLqysrK8vDzfz+onwyFoENzd3S0vLw8zYYiNje3q6lq1atWjR49gcoNvPKWu46bRaFQqlUwmYzCYnp4eTmYmJiYvXrwIDQ11dXWFyY2Pjlu2uxUQBIlEIp1ORzIORzLHQRVmhlF6cTgcmUxGMnvnCh+jI6nruKFeG+HqCZLpDELFZL4okkgkDAZDJpMVFBS4JuE6FUJeFGVYWGhk3tbWhsVi4UdEEEgmQX0mbB8AbXskkUh0Oh06+MxTHwE/neklK8H2JP4f8Md4WU5Bc4LrMV5ZpaOjQ1dXV11d3czMbOnSpZ8+fQIRHONlPkMOg6yeSuPK9OnToRZ5xowZd+/ehQ74wKvBEi6BE1yP8Uo+lpaWkDhWVlaMyIbw4rDEI+CEuI7xShSpqalKSkqMwxfl5eUgt2O8EBJ7jFfSMDAwgA4LODs7Z2Vl0el0+GO8LCFdYJDVBrOgoGDQoEEaGhrm5uY+Pj6FhYUgt2O8nAImssD1GK8M8/vvvysoKBgYGNjY2DBC1XAtQllZWTAR6CC4HuOVVSIiIgAA0NTUtLCwOHLkSG1tLcjtGG91dfWkSZNUVVUZAfJ6hesxXsknJCQEAACWABHwx3gDAgIsLCxSU1PhcxbLMV6JIiEhAYfDQYEOAgICoDEe/DFeTqE/e0Usx3glitOnTzNKb1BQUH19PcjtGC/zNIfrlEcmj/Fqa2uzRCyBP8bLUAnJdKatrW39+vUicV2ygS4tYWlI4ftchrBIpkJ9eYxXopg4cSI0tbSxsXn06BEUFx5mRMReqWX+GC8GgxkyZMi4ceMYUSMZx3iXLVvGKSGzUEimM5DmRCLR29sb+kQ4O/tMTEwKCws5fYvFYuPi4rhmUlRUNHXqVKH4I13gcLi5c+dCQWcYi98mJibwF8nDCM6gublZTU1NOF5KGwsXLvTy8lq2bBmzAtA7cE5bKv7++2/oMnh4qqqqOJ2ukhYWLFiwfv36JUuWML9QlZOTI5FInA5BBAQEwBykYlBeXs73JiCZQVNTc/Xq1T4+Pr/99hvjw2HDhpWUlDB/wg7MmUqIkpIS6Kqyfo6tra2np+e8efMYG1GNjY1LS0s52dvb29vb23PNtr29XUVFRWheShIqKipOTk4+Pj6TJk1ifDhs2LC3b99ySoLFYhMTE7nm/OPHDyR37Mok8+bN27Ztm7u7O3OxoVKp8Klgzpsz6Ozs5HS/qmxjbGy8c+dOHx8f5obO1NS0pKSEU3QIfX39jx8/cs1ZBhrPkSNH7tu3b8OGDcw1DhJnxIgRvSYJCgoKCgrimrMMiCMgOjo6mzdv9vHxgaZVEIaGhuXl5ebm5r0mYR6Ec53jlJeXGxoaCsVVKcXc3Nzf33/Tpk3MFRlSmJMyzApznfKUl5ePGjVKKK5KCIqKinZ2dhs3bpw2bRpjq46hoWFqaiqnJMwqcZ3O9NsyCTWkGzduHDp0KONDQ0PDiooK5urPDLOwXKdC5eXlEyZMEIqr0sW8efN8fX3d3d0HDBjA+JBGo3GyZ6/U8JMgGo3G3+kxCWHMmDH79+/fsGEDc70zMDCoqqoyMjKC2aHPLBSS6Qyk+Y8fP8zMzP7vI2GtVg4fPlyQHOrq6rBYLNfX3f0KJycnAdewb968CX93Z38jNDQ0JSVFkBxIJJKnp6ew/JEorly58vjxY0FyoFAoSC6N6Z80NDTAX42NBC8vL5gbFfs5Li4u0MZ4vrl///79+/eF5Y9U4OLiImAOJ0+efP36tVCckQ0CAwM/f/4sSA6NjY19FjddKigvL9+2bZuAmSxZsgT+unkppba2dvPmzQJmsnz5cphLePotZWVlu3btEkpWK1asQBVmJz8/H/62d+QsXboU/ooe2YBKpXK9kxQhx44de/v2rVCykgEyMjKEtbvZ2dlZwLGoLOHn51daWip4PkLp6SSQO3fuPHr0CATB9evX19TUCJ5hWVnZzp07QRC8ceMG41wOl0PUCHF0dPz+/Xt2djbfOdy7d2/gwIE2NjZC8Uc2mDp16ps3bwTJISkpiXG3NAoAAI6OjtHR0YLk8Pjx4zlz5gjLH4nC3t4+JiZGkByePXuGVmFOaGtrNzY2co0rB0NnZyeFQhFR1DkZwMrKKjMzU5AcEhISFi5cKCx/pAITE5OioiJBcnj58uW0adOE5Y8M4Ozs/ODBA0FyiIqKIhAIwvJHBjAyMvr58yfMe2+utLS0yMvLy+S+XT09vbq6OoRhvHqlvb2dRqNJQnx3ScPExOTHjx+C50MmkykUCqowOxYWFvn5+YLn09XVhcViebqdQ0rB4XDq6upNTU2CZ/X+/fspU6YIno9sMHPmzNevXwueT11dnZ6eHtcQdf0HR0dHJGdEuBIfHy+T4yIbGxtor669vX1SUpLgGcbFxUFb/168eDFr1izoQ+EUx4kTJ86ZM2fPnj38Je/s7AwODvbz84O5LKwfsnHjxlOnTvGd/P379+bm5sybaVEMDAxAEITigPIBjUaDrmESrlcSgq6uroqKytevX/lLTqfTL1265OXlJVyvZIk1a9ZcuHCB7+QnTpxYv369EP2RMdavXw+9mOUv+efPnw0MDNTV1YXrlYSzY8cOKH4Nf8THx9va2iK8w66fMG7cuNLS0traWv6Sd3d3JyUl2dnZCdcraWf58uWXLl3iO3lQUBDMlbXSjre3d2hoKN/Jjx075uPjI0R/ZIk//vgjNjZWwEwuXbrENUZHv2X69OlPnz4VMJPQ0NAVK1YIxR/JZ9WqVWfPnhUwk6ysLAsLC3RNigEOhzMzM/vw4YOA+Zw+fVoSLn6VHGbNmpWSkkKhUATJhEqlJiQkzJ07V1heSQ76+vp1dXUUCmX+/PlxcXEwp56RQCaTU1NTZ86c2dXVRaFQ/v+MRvAdgxDZ2dkYDOb27dt8pN24caOOjk5bW5uwnJEZTp48+fDhQz4SkslkPB7f2toqdJeknYqKCnd3d/7OR589e/bu3btCd0ly+PXrF99nIcPDw69fvy5sj2QNNze3iooKPhKWlJTI6vlxIXL+/PnIyEg+ElIoFHt7e67Xd8gke/fuTU9P5yNhR0fHokWLoOt9UJjJy8tbs2YNf2n379///Plz4fojA9DpdEdHR/4OueTn569atUroLkkULi4uVVVVfCQsKipasWKF0P2RGSgUCh6Ph+7e4Y/29nYHBwepjisvUrq7u/F4vCDHHpubmwWPRyFdLF26VMB7Ap2cnNBJNwstLS3Ozs6C5FBZWbl69Wph+SMzxMXFnTt3TpAcwsPDHzx4ICx/JI2YmJiLFy+CIHjnzp0rV64IktXJkycTExNBEAwODmYeSQptsQ8Ewf379yspKb17946nVP/88w8Oh0tKShKiJzIDjUZbunRpbm4urwk3b96clpYmCpdkgOjo6L/++ovXVK9evVq3bp0I3JEsYmNjAwICeE315s0btIdDwq9fvwgEQmdnJ0+p2tvb8Xh8XV2diLySGeh0+sqVK7Ozs3lN6Ovr22/7IAqF4uzszOsl2nQ63cvLiw+p+wkXLly4cOECr6n4a377CT9//nR0dOR1cbmpqcnOzk7mQ51WV1fb29vzGpQQ6lmge1FROJGcnLx7927+0kJd0sePH4XrkowRExPDd+Q+aJaUl5cnXJcknJKSEjc3t56eHv6Sh4SEyNjNxcLi0qVLfEfuo1AoTk5OZWVlQvVIRhCkkhYWFrq6usr2+xJHR8e2tjYajebq6gpzBTw8ubm5S5cuBUGwubmZ5RJkYS720Wg0FxcXNTU1hLMmOp0eEhKCxWKPHz8uRDdkjLa2Nnt7+w8fPiC0p1KpO3bsuHTpkki9knaOHj169OhR5G1Hamqqm5tbP9nDcuLEicDAQOTiZGRkODs7w1xaj8LMhw8fXFxcoDvXkVBXV+fg4MDHin//pKOjw9HREXnYaRqN5u/vD71V67fU1dXZ2dkVFBQgtCeRSGvXro2OjhapV9LOrl27wsLCkNvHxMR4eXmhgb1hePXq1eLFi5EfWaiqqiIQCIWFhSL1SkJ49+6dq6trc3MzQvuamhp7e/v+tkrCHyEhIdeuXeMj4ZEjR9BVFST8+eef9+7d4yPh/v37o6KihO6P5JOcnOzn58dHwsTERCh+P0qv+Pr68nFZKJ1O9/HxSU1NFYVLMkBrayt/L5aam5vxeLzMH7t5+fIldBlUY2MjHo9H3o8zqK+vx+Px0HbdTZs2ffr0iflbYS72gSBIoVA2bdqEw+G2b98O/7cpLCy0t7dXUFCIiIgQrg+yR2dn54oVK/755x8qlQpvWVpa6ubmJuCdqv2EmzdvrlixguuxIDKZHBwc7OPjI8g5Dqnj3r17y5Ytq66uhjejUCgnT57csGFDP1kGFRZFRUV2dnavXr3iapmeno7H44VylVX/gUgkent7X7x4ket77/Ly8sWLF6NXloMg2Nzc7ObmdufOHa6r/AUFBQ4ODi9fvuwbx6Sas2fPbty4kes4taura9++ffv27UNX+riSl5dnZ2f3/v17rpZPnjwhEAgCHnaTLr59+2ZnZ5eVlcXVMiUlBY/Ho5tQkOPv7x8SEoLcnkql7tq16+TJk6JzSZag0+l+fn48baqiUCi+vr6hoaGi80rCuXHjxrp163gafl++fHn9+vVc55L9GQqFsnr16hs3biBPQiQSV61adefOHdF5JQOUlpbi8fj8/HzkSYqKivB4fFFRkei8khz++usvqAgVFhbi8Xie9vfl5eUxporXr18/c+YMi4GQF/sgbt26ZWBgoKGhsXXr1tTUVObb0FtaWh4+fLh06VJ5eflJkya9efNGFA7IHnQ6PTIyEo/H37lzp9fX2p8/f96xY8eaNWvKy8v73j0p5evXr0uXLj1w4ECvTUlDQ8Ply5cJBAJ0AL6/UVhYuHz58oCAgF73+zQ1NV27do1AIKALJfzR3t5+8OBBT0/PtLQ09jUpCoWSkpKyZMmSI0eO8HoyCwUEQTqdfu/ePTwef+vWrV5P8OXl5e3evXvVqlXoQiqDnp6eiIgIBweH6Oho9pPmdDr9w4cPmzdv9vHx4S90Wv/k48ePrq6uR48e7XVh5devXxcuXHBwcHjx4kVfeya1tLS07N2718vLKzMzk33WSiKRnjx54u7ufuzYMebBZz+hra1t//79y5cvT09PZ+9ZyGTys2fPPDw8goKC0M34vALdz4aky8jLy3NyckLfu/NKeHi4l5cXkrjG//77r4ODAxre9NWrVwQCAcnLj1+/fq1du/bUqVN94JW0Q6fTg4ODN27cWFtby9U4KysLj8cjP03Sn2ltbfXw8AgLC+O6e6anp+fSpUtubm4yv6ePAZ1OX7FiBbQs1tjY6OrqevnyZa7r8hQK5eLFi0uWLIGWhlJTU9evX89uhgH5vbsQHiKRGB4efv/+/ZycHBAEtbW11dTUGhoaOjo6FBQUbGxsvL29lyxZgsFgRPF0WYVMJsfHxz958qSrq0tJSUlPT49EIjU0NPT09IwbN87d3f23334Tt4/SR05OTnR0dHFxsaKioo6Ojry8fH19fXd3t6amJoFAsLOz6893TX7+/DkqKqqwsFBBQUFHR0dBQaG+vp5EIqmrq+PxeAKBICcnJ24fpZi6urqoqKhXr14BAKCmpqalpdXY2NjR0QEAgK2trZubm7a2trh9lGIoFEpCQkJycnJnZyfUYJLJ5IaGBgqFMmbMGHd39zFjxojbR4mjq6srNjb2+fPnJBJJRUVFT0+vs7Ozqampp6fH0tLS3d3d7P+1dzchqawBGIAd0xTF8qdxjCSMsznQpqhoF5m5KCorSAxa1CJoH0RRuJE2ta1l0KasFopSu7IIKgoJIVq0EArKnJlSQcufMb0L4XKo+3PuxZjM99n6Mb58qN/4zvj54wffGUvPycmJy+W6vb2VSqUkSQqFQpZlk8kkRVEWi8VsNuNc6L8KhUI7Ozunp6cEQSiVSpVKxbJsIpEQCoVGo3FkZESlUvGdkTfhcLiwshAEUV1drVarC5NDEERhZdFoNHxnLEl3d3cOh0Mul9tstvb29nd/ZprNZo+Pj51Op0wmm5+f12q1fOUsXcFg0OFwqNVqm83W1tb27oOR47ijoyOn06lWq+fm5vAyFggEsVhsaWnp/v5+dHS0q6tLIpG8G3B5ebm9vf3w8LCwsPDz509eQpai6+vrxcVFg8FgtVqbmprePZpKpQ4ODpxOZ0NDw8zMjEKh4CVkycnn8y6Xa319va+vb2BgoLa29t0Amqa9Xq/X6x0bG7NarWV1apRMJqempnp6ego77m1tbW1ublosloGBgY+ryePjo8fj2dvbm5iYGB4eFggEa2tr5+fnKysrlZWV7wZ/Vtn3p1AoFAgEwuFwJBKhKKqurq61tbWqqupTn7QcpNNpmqZlMhnqgCKiaZrjOIqixGIx31m+HIZhMpkMJueTJBIJhmEoipLL5Xxn+YYymQxN04Wqhe8sJSOVStE0XVVVVc69SXHl8/lwOJzL5XQ6XTlfRiqueDz+9PREUZRMJuM7y5dTWFl0Oh0mp1iCwaDL5bq4uBCLxRKJpKamhmEYjuOy2WxHR8fg4KBer+c7Y2m7ublxu91+v7+yslIqlWo0Gpqms9ns29tbZ2fn4ODgx46gzDEM4/F4fD4fQRAikUir1cZisZeXF47jmpqahoaGGhsb+c5Ykq6urtxu99XVlUgkksvlSqWy8GYnCMJkMlksFnwH/x84jvP5fF6vl2XZiooKkiQJgmBZNpvNajSa/v7+7u7uj41VOcjn8w6HIx6P2+12hUJR+JnX7u7u8/OzSCQiSTKfz7Ms+/b2RpKkxWIxGo1isTgajdrtdoPBMD09/ZeH/fSyDwAAAAAA4DtJJpPPz88kSX68owqK4vX1NRKJYIZ/XyaTYRhGpVLhynFxvby8RKNR3PFQdCzL5nI5iqL4DvJVnJycLC8vm83mycnJX0tPmqaFQuGvtyykUqnV1dWzs7PZ2dnW1ta/OyDKPgAAAAAAAAAAAD653e6NjQ2tVmsymVpaWurr6wsbR+RyudvbW7/fv7+/H4vFxsfHe3t7//lQKPsAAAAAAAAAAAD4Fw6HDw8PA4FAKBRKp9MCgUAikej1+ubmZqPR+JsbE6HsAwAAAAAAAAAA+CaE/z4EAAAAAAAAAAAASsEfbkNERGKwNHQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from subprocess import check_call\n",
    "from IPython.display import Image\n",
    "f.draw('tmp.dot', portrait=True)\n",
    "check_call(['dot','-Tpng','-Gdpi=200','tmp.dot','-o','tmp.png'])\n",
    "Image(filename='tmp.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observation Length:  30\n",
      "Backtraced State Sequence\n",
      "([0, 1, 2, 3, 3, 3, 3, 4, 5, 5, 6, 7, 8, 8, 8, 9, 10, 11, 11, 11, 11, 11, 12, 13, 14, 14, 14, 14, 14, 15], 'p eh p er z')\n",
      "Backtraced Phone Sequence\n",
      "['p_1', 'p_2', 'p_3', 'eh_1', 'eh_1', 'eh_1', 'eh_1', 'eh_2', 'eh_3', 'eh_3', 'p_1', 'p_2', 'p_3', 'p_3', 'p_3', 'er_1', 'er_2', 'er_3', 'er_3', 'er_3', 'er_3', 'er_3', 'z_1', 'z_2', 'z_3', 'z_3', 'z_3', 'z_3', 'z_3']\n"
     ]
    }
   ],
   "source": [
    "import observation_model\n",
    "import math\n",
    "\n",
    "\n",
    "def get_hmm_label(fst, arc):\n",
    "    arc_input = arc.ilabel\n",
    "    hmm_state = fst.input_symbols().find(arc_input)\n",
    "    return hmm_state\n",
    "\n",
    "class MyViterbiDecoder:\n",
    "        \n",
    "    NLL_ZERO = 1e10  # define a constant representing -log(0).  This is really infinite, but approximate\n",
    "                     # it here with a very large number\n",
    "    \n",
    "    def __init__(self, f, audio_file_name):\n",
    "        \"\"\"Set up the decoder class with an audio file and WFST f\n",
    "        \"\"\"\n",
    "        self.om = observation_model.ObservationModel()\n",
    "        self.f = f\n",
    "        \n",
    "        if audio_file_name:\n",
    "            self.om.load_audio(audio_file_name)\n",
    "        else:\n",
    "            self.om.load_dummy_audio()\n",
    "        print(\"Observation Length: \", self.om.observation_length())\n",
    "        self.initialise_decoding()\n",
    "        \n",
    "    \n",
    "    def initialise_decoding(self):\n",
    "        \"\"\"set up the values for V_j(0) (as negative log-likelihoods)\n",
    "        \n",
    "        \"\"\"\n",
    "        \n",
    "        self.V = []\n",
    "        self.B = []\n",
    "        \n",
    "        for t in range(self.om.observation_length()+1):\n",
    "            self.V.append([self.NLL_ZERO]*self.f.num_states())\n",
    "            self.B.append([None]*self.f.num_states())\n",
    "        \n",
    "        # The above code means that self.V[t][j] for t = 0, ... T gives the Viterbi cost\n",
    "        # of state j, time t (in negative log-likelihood form)\n",
    "        # Initialising the costs to NLL_ZERO effectively means zero probability    \n",
    "        \n",
    "        # give the WFST start state a probability of 1.0   (NLL = 0.0)\n",
    "        self.V[0][f.start()] = 0.0\n",
    "        \n",
    "        # some WFSTs might have arcs with epsilon on the input (you might have already created \n",
    "        # examples of these in earlier labs) these correspond to non-emitting states, \n",
    "        # which means that we need to process them without stepping forward in time.  \n",
    "        # Don't worry too much about this!  \n",
    "        self.traverse_epsilon_arcs(0)\n",
    "        \n",
    "    def traverse_epsilon_arcs(self, t):\n",
    "        \"\"\"Traverse arcs with <eps> on the input at time t\n",
    "        \n",
    "        These correspond to transitions that don't emit an observation\n",
    "        \n",
    "        We've implemented this function for you as it's slightly trickier than\n",
    "        the normal case.  You might like to look at it to see what's going on, but\n",
    "        don't worry if you can't fully follow it.\n",
    "        \n",
    "        \"\"\"\n",
    "        \n",
    "        states_to_traverse = list(range(self.f.num_states())) # traverse all states\n",
    "        while states_to_traverse:\n",
    "            \n",
    "            # Set i to the ID of the current state, the first \n",
    "            # item in the list (and remove it from the list)\n",
    "            i = states_to_traverse.pop(0)   \n",
    "        \n",
    "            # don't bother traversing states which have zero probability\n",
    "            if self.V[t][i] == self.NLL_ZERO:\n",
    "                    continue\n",
    "        \n",
    "            for arc in self.f.arcs(i):\n",
    "                \n",
    "                if arc.ilabel == 0:     # if <eps> transition\n",
    "                  \n",
    "                    j = arc.nextstate   # ID of next state  \n",
    "                \n",
    "                    if self.V[t][j] > self.V[t][i] + float(arc.weight):\n",
    "                        \n",
    "                        # this means we've found a lower-cost path to\n",
    "                        # state j at time t.  We might need to add it\n",
    "                        # back to the processing queue.\n",
    "                        self.V[t][j] = self.V[t][i] + float(arc.weight)\n",
    "                  \n",
    "                        if j not in states_to_traverse:\n",
    "                            states_to_traverse.append(j)\n",
    "\n",
    "    \n",
    "    def forward_step(self, t):\n",
    "        \n",
    "        # TODO - exercise\n",
    "        \n",
    "        for i in self.f.states():\n",
    "            if self.V[t-1][i] == self.NLL_ZERO:\n",
    "                continue\n",
    "            else:\n",
    "                for arc in self.f.arcs(i):\n",
    "                    if arc.ilabel == 0:     # if not <eps> transition    \n",
    "                        continue\n",
    "                    else:\n",
    "                        transition_prob = float(arc.weight)\n",
    "                        emission_prob = -float(self.om.log_observation_probability(get_hmm_label(self.f, arc), t))\n",
    "                        j = arc.nextstate\n",
    "                        combined = transition_prob + emission_prob + self.V[t-1][i]\n",
    "                        \n",
    "                        if combined < self.V[t][j]:\n",
    "                            self.V[t][j] = combined\n",
    "                            self.B[t][j] = i\n",
    "                        \n",
    "    \n",
    "    def finalise_decoding(self):\n",
    "        \n",
    "        # TODO - exercise\n",
    "        for i in self.f.states():\n",
    "            end_prob = min(float(self.f.final(i)), self.NLL_ZERO)\n",
    "            self.V[-1][i] += end_prob # -1 represents the last time step.\n",
    "            \n",
    "#         self.V[-1][:] = [x if x != math.inf else self.NLL_ZERO for x in self.V[-1]]\n",
    "    \n",
    "    def decode(self):\n",
    "        \n",
    "        self.initialise_decoding()\n",
    "        t = 1\n",
    "        while t <= self.om.observation_length():\n",
    "            self.forward_step(t)\n",
    "            self.traverse_epsilon_arcs(t)\n",
    "            t += 1\n",
    "        \n",
    "        self.finalise_decoding()\n",
    "    \n",
    "    def backtrace(self):\n",
    "        \n",
    "        # TODO - exercise \n",
    "        \n",
    "        # complete code to trace back through the\n",
    "        # best state sequence\n",
    "        \n",
    "        # You'll need to create a structure B_j(t) to store the \n",
    "        # back-pointers (see lectures), and amend the functions above to fill it.\n",
    "        best_state_sequence = []\n",
    "        \n",
    "        state = np.argmin(self.V[-1])\n",
    "        best_state_sequence.append(state)\n",
    "        b_idx = self.B[-1][state]\n",
    "        \n",
    "        for s in range(self.om.observation_length()-1, 0, -1):\n",
    "            state = self.B[s][b_idx]\n",
    "            b_idx = state\n",
    "            best_state_sequence.append(state)\n",
    "        \n",
    "        best_state_sequence.reverse()\n",
    "        \n",
    "        best_output = []\n",
    "        for start in range(0, len(best_state_sequence)-1):\n",
    "            arc_list = list(self.f.arcs(best_state_sequence[start]))\n",
    "            for stop in range(start + 1, len(best_state_sequence)):\n",
    "                arc_found = False\n",
    "                for arc in arc_list:\n",
    "                    if int(arc.nextstate) == best_state_sequence[stop]:\n",
    "                        arc_found = True\n",
    "                        if arc.olabel != 0:\n",
    "                            best_output.append(arc.olabel)\n",
    "                        break\n",
    "                if arc_found:\n",
    "                    # TO avoid repetition getting counted multiple times\n",
    "                    break\n",
    "        emissions = \"\"\n",
    "        for label in best_output:\n",
    "            emissions = emissions + \" \" + self.f.output_symbols().find(label)\n",
    "        emissions = emissions.strip()\n",
    "        return best_state_sequence, emissions\n",
    "    \n",
    "    def get_backtraced_phones(self, state_sequence):\n",
    "        \n",
    "        phone_sequence = []\n",
    "        for start in range(0, len(state_sequence)-1):\n",
    "            arc_list = list(self.f.arcs(state_sequence[start]))\n",
    "            for stop in range(start+1, len(state_sequence)):\n",
    "                arc_found = False\n",
    "                for arc in arc_list:\n",
    "                    if int(arc.nextstate) == state_sequence[stop]:\n",
    "                        arc_found = True\n",
    "                        phone = self.f.input_symbols().find(arc.ilabel)\n",
    "                        phone_sequence.append(phone)\n",
    "                        break\n",
    "                if arc_found:\n",
    "                    break\n",
    "        return phone_sequence\n",
    "        \n",
    "\n",
    "# to call the decoder (in a dummy example)\n",
    "# f will be a WFST that you have created in a previous lab\n",
    "decoder = MyViterbiDecoder(f, '')   # empty string '' just means use dummy probabilities for testing\n",
    "decoder.decode()\n",
    "# print(decoder.V)\n",
    "bt_states = decoder.backtrace()\n",
    "bt_phones = decoder.get_backtraced_phones(bt_states[0])\n",
    "\n",
    "print(\"Backtraced State Sequence\")\n",
    "print(bt_states)\n",
    "\n",
    "print(\"Backtraced Phone Sequence\")\n",
    "print(bt_phones)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises &ndash; Lab 3\n",
    "\n",
    "The `__init__()`, `initialise_decoding()` and `decode()` functions have been completed in the template above.\n",
    "You should aim to complete the `forward_step()` and `finalise_decoding()` functions for Lab 3.  Don't worry about implementing the back-trace in Lab 3.\n",
    "\n",
    "You should draw on your solutions to Lab 2 &ndash; the main difference now is that, rather than simply sampling a single path and computing its likelihood, you'll need to compute and store, at every time step $t$, and for every state in the WFST, the likelihood of the best path reaching that state after $t$ time steps.  For how to do this using the Viterbi algorithm, see Lecture 5, slides 11 onwards.\n",
    "\n",
    "Test your algorithms on an WFST that recognises the word \"*pepper*\" and one that recognises any word in the vocabulary.\n",
    "\n",
    "\n",
    "## Exercises &ndash; Lab 4\n",
    "\n",
    "Now complete the `backtrace()` function to allow the best path to be recovered.  As noted in the code, you'll need to create a structure to store $B_j(t)$, which stores the identity of the best preceding state reaching state $j$ at time $t$ .  You can follow a similar method to storing $V_j(t)$, given in the code already.  You'll need to add it to the `initialise_decoding()`, `forward_step()` and `finalise_decoding()` functions.\n",
    "\n",
    "Once you are happy that your function works, you should amend your code so that you can also recover the sequence of *output* symbols on your WFST's best path as well.  This should allow you to produce your first word-recognition result!\n",
    "\n",
    "### Working on real speech\n",
    "\n",
    "You can now find the first batch of recordings made by ASR students in `/group/teaching/asr/labs/recordings`.  You can test your decoder on real speech data by passing the full path to the WAV file when you create your `MyViterbiDecoder` object.  The transcriptions are also available in the same folder.\n",
    "\n",
    "When working with real speech, you may want to modify your code to print only the word output labels! \n",
    "\n",
    "If you are interested, the observation model (kindly supplied by Andrea) is a monophone time-delay neural network trained using the lattice-free MMI criterion, on the WSJ corpus of read speech.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}