DishonestCasino.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import random\n",
    "import math\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Входные данные:\n",
    "\n",
    "pi - распределение;\n",
    "\n",
    "transit - матрица переходных вероятностей;\n",
    "\n",
    "emis - матрица эмиссий;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi=np.array([2.0/3,1.0/3])\n",
    "transit=np.array([[0.95,0.05],[0.1,0.9]])\n",
    "emis=np.array([[1.0/6,1.0/6,1.0/6,1.0/6,1.0/6,1.0/6],[0.1,0.1,0.1,0.1,0.1,0.5]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Генерация следующего состояния\n",
    "def next_state(weights):\n",
    "    choice = random.random()\n",
    "    for i, w in enumerate(weights):\n",
    "        choice -= w\n",
    "        if choice < 0:\n",
    "            return i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Генерация последовательности использования костей\n",
    "def hidden_seq(L):\n",
    "    out=[0 for i in range(L)]\n",
    "    out[0]=next_state(pi)\n",
    "    for i in range(1,L):\n",
    "        out[i]=next_state(transit[out[i-1]])\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Генерация последовательности выбросов\n",
    "def obs_seq(hidden, L):\n",
    "    out=[0 for i in range(L)]\n",
    "    for i in range(L):\n",
    "        out[i]=next_state(emis[hidden[i]])\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Печать скрытых состояний для алгоритма Витерби\n",
    "def print_LF(seq):\n",
    "    for i in range(len(seq)):\n",
    "        if (seq[i] == 0):\n",
    "            print('F', sep = ', ', end = '')\n",
    "        else:\n",
    "            print('L', sep = ', ', end = '')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Алгоритм Витерби в логарифмическом пространстве\n",
    "def viterbi_log(obs_seq,pi, transit, emis):\n",
    "    T = len(obs_seq)\n",
    "    N = transit.shape[0]\n",
    "    delta = np.zeros((T, N))\n",
    "    psi = np.zeros((T, N))\n",
    "    delta[0] = np.log(pi[0]) + np.log(emis[:, obs_seq[0]])\n",
    "    for t in range(1, T):\n",
    "        for j in range(N):\n",
    "            delta[t,j] = np.max(delta[t-1] + transit[:,j]) + emis[j, obs_seq[t]]\n",
    "            psi[t,j] = np.argmax(delta[t-1] + transit[:,j])\n",
    "\n",
    "    states = np.zeros(T, dtype=np.int32)\n",
    "    states[T-1] = np.argmax(delta[T-1])\n",
    "    for t in range(T-2, -1, -1):\n",
    "        states[t] = psi[t+1, states[t+1]]\n",
    "    return states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Алгоритм Витерби с применением масштабирования\n",
    "def viterbi_scaling(obs_seq,pi, transit, emis):\n",
    "    T = len(obs_seq)\n",
    "    N = transit.shape[0]\n",
    "    delta = np.zeros((T, N))\n",
    "    psi = np.zeros((T, N))\n",
    "    delta[0] = pi * emis[:,obs_seq[0]]\n",
    "    m = np.zeros(T)\n",
    "    m[0] = max(delta[0])\n",
    "    delta[0] /= m[0]\n",
    "    for t in range(1, T):\n",
    "        for j in range(N):\n",
    "            delta[t,j] = np.max(delta[t-1] * transit[:,j]) * emis[j, obs_seq[t]]\n",
    "            psi[t,j] = np.argmax(delta[t-1] * transit[:,j])\n",
    "        m[t] = max(delta[t])\n",
    "        delta[t] /= m[t]\n",
    "\n",
    "    states = np.zeros(T, dtype=np.int32)\n",
    "    states[T-1] = np.argmax(delta[T-1])\n",
    "    for t in range(T-2, -1, -1):\n",
    "        states[t] = psi[t+1, states[t+1]]\n",
    "    return states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Печать графика использования кости (желтый - честная, красный - нечестная)\n",
    "def print_chart(path):\n",
    "    for i in range(L):\n",
    "        if path[i] == 0:\n",
    "            plt.axvline(x = i, color ='yellow' )\n",
    "        else:\n",
    "            plt.axvline(x = i, color = 'red')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Печать графика предсказанного Витерби использования кости (желтый - честная, красный - нечестная)\n",
    "def print_viterbi(path, orig, s, L):\n",
    "    L = len(path)\n",
    "    if L < 1000:\n",
    "        print_LF(path)\n",
    "    print('\\nЧастота совпадения:', s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Печать сравнения логарифмического и масштабированного Витерби\n",
    "def print_viterbi_diffs(path_log, path_sc, orig, s_log, s_sc, L):\n",
    "    L = len(path_log)\n",
    "    if L < 1000:\n",
    "        print('Original \\n')\n",
    "        print_LF(orig)\n",
    "        print('\\n')\n",
    "    print('Viterbi_logarithm')\n",
    "    print_viterbi(path_log, orig, s_log, L)\n",
    "    print('\\n')\n",
    "    print('Viterbi_scaling')\n",
    "    print_viterbi(path_sc, orig, s_sc, L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Нахождение точности предсказания Витерби\n",
    "def find_accuracy(orig, viterbi, L):\n",
    "    s = 0\n",
    "    for i in range(L):\n",
    "        if orig[i] == viterbi[i]:\n",
    "            s+=1\n",
    "    return s / L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Обработка n последовательностей размера L\n",
    "def process_sequences(L, n):\n",
    "    s_log = 0\n",
    "    s_sc = 0\n",
    "    for i in range (n):\n",
    "        hidden=np.array(hidden_seq(L))\n",
    "        observed = obs_seq(hidden, L)\n",
    "        path_log=viterbi_log(observed,pi,transit, emis)\n",
    "        path_sc=viterbi_scaling(observed,pi,transit, emis)\n",
    "        s_log += find_accuracy(hidden, path_log, L)\n",
    "        s_sc += find_accuracy(hidden, path_sc, L)\n",
    "    s_log /= n\n",
    "    s_sc /= n\n",
    "    print_viterbi_diffs(path_log, path_sc, hidden, s_log, s_sc, L)\n",
    "    return s_log, s_sc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Подсчет результатов для заданных входных данных\n",
    "def calculate_results(Lens_obs, Count):\n",
    "    t = len(Lens_obs)\n",
    "    results_log = np.zeros((2, t))\n",
    "    results_sc = np.zeros((2, t))\n",
    "    for i in range (len(Lens_obs)):\n",
    "        print('Входные данные: ', Count[i], ' последовательностей; ', Lens_obs[i], ' элементов в каждой;', )\n",
    "        start_time = time.time()\n",
    "        s_log, s_sc = process_sequences(Lens_obs[i], Count[i])\n",
    "        results_log[0, i] = Lens_obs[i]\n",
    "        results_log[1, i] = s_log\n",
    "        results_sc[0, i] = Lens_obs[i]\n",
    "        results_sc[1, i] = s_sc\n",
    "        print('\\n')\n",
    "        print('Время:', (time.time()-start_time) / 60, 'минут \\n\\n\\n')\n",
    "    return results_log, results_sc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Входные данные:  40  последовательностей;  300  элементов в каждой;\n",
      "Original \n",
      "\n",
      "FFFFLLLLFFFFFFFFFFFFLLLLFFFFFFFFFFFLLLLLFFFFFFFFFLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLFFFFFFFFFFFFLLLLLFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFLLLLFFFFFFFFFFFFFFLLFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLLLLLLLLLFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLFFFFFF\n",
      "\n",
      "Viterbi_logarithm\n",
      "FFFFFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL\n",
      "Частота совпадения: 0.7545833333333332\n",
      "\n",
      "\n",
      "Viterbi_scaling\n",
      "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLLLLFFFFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLLL\n",
      "Частота совпадения: 0.7781666666666668\n",
      "\n",
      "\n",
      "Время: 0.012573607762654622 минут \n",
      "\n",
      "\n",
      "\n",
      "Входные данные:  40  последовательностей;  1000  элементов в каждой;\n",
      "Viterbi_logarithm\n",
      "\n",
      "Частота совпадения: 0.7695000000000001\n",
      "\n",
      "\n",
      "Viterbi_scaling\n",
      "\n",
      "Частота совпадения: 0.8035500000000001\n",
      "\n",
      "\n",
      "Время: 0.052584278583526614 минут \n",
      "\n",
      "\n",
      "\n",
      "Входные данные:  40  последовательностей;  2000  элементов в каждой;\n",
      "Viterbi_logarithm\n",
      "\n",
      "Частота совпадения: 0.7559875000000001\n",
      "\n",
      "\n",
      "Viterbi_scaling\n",
      "\n",
      "Частота совпадения: 0.8002625000000002\n",
      "\n",
      "\n",
      "Время: 0.09787732362747192 минут \n",
      "\n",
      "\n",
      "\n",
      "Входные данные:  40  последовательностей;  3000  элементов в каждой;\n",
      "Viterbi_logarithm\n",
      "\n",
      "Частота совпадения: 0.7474999999999999\n",
      "\n",
      "\n",
      "Viterbi_scaling\n",
      "\n",
      "Частота совпадения: 0.79455\n",
      "\n",
      "\n",
      "Время: 0.12051969369252523 минут \n",
      "\n",
      "\n",
      "\n",
      "Входные данные:  40  последовательностей;  10000  элементов в каждой;\n",
      "Viterbi_logarithm\n",
      "\n",
      "Частота совпадения: 0.759255\n",
      "\n",
      "\n",
      "Viterbi_scaling\n",
      "\n",
      "Частота совпадения: 0.7978050000000002\n",
      "\n",
      "\n",
      "Время: 0.4827980558077494 минут \n",
      "\n",
      "\n",
      "\n",
      "Входные данные:  40  последовательностей;  100000  элементов в каждой;\n",
      "Viterbi_logarithm\n",
      "\n",
      "Частота совпадения: 0.7558069999999999\n",
      "\n",
      "\n",
      "Viterbi_scaling\n",
      "\n",
      "Частота совпадения: 0.7946280000000001\n",
      "\n",
      "\n",
      "Время: 4.262119710445404 минут \n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#results_log, results_sc - массивы точек формата (размер посл-ти, точность)\n",
    "Lens_obs = [300, 1000, 2000, 3000, 10000, 100000]\n",
    "Count = [40, 40, 40, 40, 40, 40]\n",
    "results_log, results_sc = calculate_results(Lens_obs, Count)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2f060327b70>]"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHrJJREFUeJzt3X+UHWWd5/H3t393Oh3SIY3mByFhgSxROQQDEwTFH6DIeGDYddZkM0dRZ7PHFQ6yciSIgRCHw4J7zIwuCoyoDDMSkVU2ywTBgcCq/DCJ8itgoAlIOsGkAwmhu5P07e7v/lF1k+rb9/at7r63b3fV53VOn6771HPrPpXqfOqpeqrqmrsjIiLpUVXpBoiIyNhS8IuIpIyCX0QkZRT8IiIpo+AXEUkZBb+ISMoo+EVEUkbBLyKSMgp+EZGUqal0A3JNnz7d586dW+lmiIhMKJs3b97j7q1x6o674J87dy6bNm2qdDNERCYUM/tT3Lo61SMikjIKfhGRlFHwi4ikjIJfRCRlFPwiIimj4BcRSRkFv4hIyiQy+O+/H159tdKtEBEZnxIZ/MuWwTXXVLoVIiLjUyKD/8ABePhh0PfIi4gMlrjgd4dMBnbvhuefr3RrRETGn8QFf0/PkemHH65cO0RExisFv4hIyiQ2+Ovq4NFHg9M+IiJyRGKD/0Mfgs5O2Lixsu0RERlvEhv8n/gEmOl0j4hIrsQG/7vfDQsXwr/9W2XbIyIy3sQKfjM738y2mlmbma3IM3+OmW0wsz+Y2bNmdkFk3tXh+7aa2SdK2fh8ssFfXw/nngtPPAFdXeX+VBGRiaNo8JtZNXAL8ElgAbDUzBbkVPsGcI+7LwSWAN8L37sgfP0e4Hzge+HyyiY6uPuxjwWDu7/5TTk/UURkYonT4z8DaHP3be7eA6wFLsqp48CUcPooYGc4fRGw1t0PufurQFu4vLKJBv/ZZwe/dZ5fROSIOME/C9geed0elkWtAv7GzNqB9cBlw3hvSUWDf9IkeM974LnnyvmJIiITS5zgtzxluU/BWQr82N1nAxcAd5lZVcz3YmbLzWyTmW3q6OiI0aTCDh0KftfVBb/nz4etW0e1SBGRRIkT/O3AsZHXszlyKifri8A9AO7+BNAATI/5Xtz9dndf5O6LWltb47c+j2iPH4Lgf+01OHhwVIsVEUmMOMG/ETjRzOaZWR3BYO26nDqvAx8DMLOTCYK/I6y3xMzqzWwecCLwu1I1Pp/c4D/ppODBba+8Us5PFRGZOIoGv7v3ApcCDwIvEly9s8XMVpvZhWG1rwL/xcyeAe4GLvHAFoIjgReAXwJfdve+cqxIVr4eP+h0j4hIVk2cSu6+nmDQNlp2bWT6BeCsAu+9AbhhFG0cluh1/BD0+EHBLyKSldg7d7M9/uZmmDlTwS8ikhWrxz+R5AY/BL3+cgV/dzc8/XRwhDFpEjQ1Bb8nTYLGxuB5QSIi40kqgn/+fPjZz0r7Ob298KMfwapVsHPQdUpHRHcGub9HOi9ap7a2tOslIsmXmuB/6y3YswemTx/d8t3hF7+Ar389OIo480z4zneCAO7uDp4L1N09cDrf746O/POGq7a2vDsXHbWIJE/igj/3Bi4YeGXPaIL/scfgqqvgqafg5JPhvvvgwgtLF4zuwRfFx9lxFKvT3Z1/5zKSL6YptJOIe1RSrI6OWkTGVuKCP9vjj4ZJNPjPynvt0dCefRauvhrWr4dZs+COO+Czn4WaEv/rmR0Jy3LJZI7sGEa7c3nzTXj99cH1fdC92UOrqRn9Ka+h5jU2QlXiLmMQGblEBn9t7cBe+HHHBWXDHeB97TW49lr453+Go46Cm2+GSy8NgmSiqq0N1uWoo8qzfPfgLum4p7zi7Fxyy6LfqxxXY2N5dy7RI0yR8S6RwZ/7n7CmBk44AV56Kd4y9uyBG26A730v6Cl+7WvBKZ6WltK3N2nMgpAt586xt3d44ylD1XnzTdi+ffC8kRy1lHusRUctUiqJDP7szVtRcR7W1tUFa9YEPfuuLvjCF+C662D27PK0VUampgamTAl+yiF71DKc8ZShdjxvvjm4bKRHLeXcueQeKUtyJTL48x12z58P//qvQW8x99x8JgM/+AFcfz3s2gUXXxz0+E8+eWzaLONL9Kjl6KPL8xnZo5ZS7FyyRy2584Z71FJdPfpTXkO9X0ct40eqgj+TCc7bn3BCUNbfD/feC9dcA21t8MEPBpdqnnnmmDZZUmgsjloOHRr+lWCF5r311uA62SvohqOhobw7Fx21xJOa4I8+s+eEE4Jv5brqKti8Gd77Xrj/frjgAv3RSDKYBSHb0FDeo5YDB0qzc9m7F9rbB9YZ6VFLKW6MHGpeEo5aUhP82Us677svuOHqoYdgzhy4805Ytiz4gxGR+GpqgmdhNTeXZ/nZo5aR3MNSaOeSW2ekRy3l3LnkG6MstdQE//TpMG1acC5/2jT49rfhS18KNqKIjD/Ro5Zp08rzGX19wxtrGarO3r2wY8fgOv398dtz+unwu9+VZ12jEhf8hw4Vvqb62muDjXPFFeW7jl1EJo7q6vIftfT0xD8qGe0jZeJKXPAX6vEDXH752LZFRNLNLDh1U19fvqOWkUjAMMVAha7jFxGRQCKDX7fPi4gUpuAXEUkZBb+ISMoo+EVEUkbBLyKSMgp+EZGUUfCLiKRM4oL/0CFdxy8iMpTEBb96/CIiQ4sV/GZ2vpltNbM2M1uRZ/4aM3s6/HnJzPZF5t1kZs+HP58pZeNz9fcHj4pV8IuIFFb0WT1mVg3cApwHtAMbzWydu7+QrePuV0TqXwYsDKf/EjgNOBWoBx4zswfcfX9J1yKUyQS/FfwiIoXF6fGfAbS5+zZ37wHWAhcNUX8pcHc4vQB4zN173b0LeAY4fzQNHkr2e0wV/CIihcUJ/lnA9sjr9rBsEDM7DpgHPBIWPQN80swmmdl04CPAsSNv7tAU/CIixcV5LHO+LyMs9IVoS4B73b0PwN0fMrPTgceBDuAJoHfQB5gtB5YDzJkzJ0aT8lPwi4gUF6fH387AXvpsYGeBuks4cpoHAHe/wd1PdffzCHYiL+e+yd1vd/dF7r6otbU1XsvzUPCLiBQXJ/g3Aiea2TwzqyMI93W5lcxsPtBC0KvPllWb2dHh9CnAKcBDpWh4Pgp+EZHiip7qcfdeM7sUeBCoBn7o7lvMbDWwyd2zO4GlwFp3j54GqgV+bWYA+4G/cfdBp3pKJfvFybqBS0SksFhfveju64H1OWXX5rxeled9Bwmu7BkT6vGLiBSXqDt3FfwiIsUp+EVEUkbBLyKSMgp+EZGUUfCLiKSMgl9EJGUSGfy6jl9EpLBEBr96/CIihSUq+LN37ir4RUQKS1Twq8cvIlKcgl9EJGUU/CIiKZPI4K+J9eg5EZF0Slzw19WB5fvOMBERARIa/CIiUljigl83b4mIDC1xwa8ev4jI0BT8IiIpk6jgP3RIwS8iUkyigl89fhGR4hT8IiIpo+AXEUkZBb+ISMokLvh1Hb+IyNASF/zq8YuIDE3BLyKSMgp+EZGUiRX8Zna+mW01szYzW5Fn/hozezr8ecnM9kXm3WxmW8zsRTP7jln5np2pG7hERIor+uR6M6sGbgHOA9qBjWa2zt1fyNZx9ysi9S8DFobTHwDOAk4JZ/8GOAd4tETtH0A9fhGR4uL0+M8A2tx9m7v3AGuBi4aovxS4O5x2oAGoA+qBWmDXyJs7NAW/iEhxcYJ/FrA98ro9LBvEzI4D5gGPALj7E8AG4I3w50F3fzHP+5ab2SYz29TR0TG8NYhQ8IuIFBcn+POdk/cCdZcA97p7H4CZnQCcDMwm2Fl81Mw+NGhh7re7+yJ3X9Ta2hqv5XnoOn4RkeLiBH87cGzk9WxgZ4G6SzhymgfgYuBJd+90907gAWDxSBoah3r8IiLFxQn+jcCJZjbPzOoIwn1dbiUzmw+0AE9Eil8HzjGzGjOrJRjYHXSqpxT6+6G3V8EvIlJM0eB3917gUuBBgtC+x923mNlqM7swUnUpsNbdo6eB7gVeAZ4DngGecff/W7LWR2QywW8Fv4jI0Ipezgng7uuB9Tll1+a8XpXnfX3Afx1F+2Lr6Ql+K/hFRIaWmDt3FfwiIvEkJvjr62HlSjj99Eq3RERkfIt1qmcimDwZVq+udCtERMa/xPT4RUQkHgW/iEjKKPhFRFJGwS8ikjIKfhGRlFHwi4ikjIJfRCRlFPwiIimj4BcRSRkFv4hIyij4RURSRsEvIpIyCn4RkZRR8IuIpIyCX0QkZRT8IiIpo+AXEUkZBb+ISMoo+EVEUkbBLyKSMgp+EZGUUfCLiKSMgl9EJGViBb+ZnW9mW82szcxW5Jm/xsyeDn9eMrN9YflHIuVPm9lBM/urUq+EiIjEV1OsgplVA7cA5wHtwEYzW+fuL2TruPsVkfqXAQvD8g3AqWH5NKANeKiUKyAiIsMTp8d/BtDm7tvcvQdYC1w0RP2lwN15yj8NPODu3cNvpoiIlEqc4J8FbI+8bg/LBjGz44B5wCN5Zi8h/w5BRETGUJzgtzxlXqDuEuBed+8bsACzGcD7gAfzfoDZcjPbZGabOjo6YjRpaE+2P8muzl2jXo6ISBLFCf524NjI69nAzgJ1C/Xq/xPwC3fP5HuTu9/u7ovcfVFra2uMJg3tgn+5gG89/q1RL0dEJIniBP9G4EQzm2dmdQThvi63kpnNB1qAJ/Iso9B5/5Jzd/Ye3Mvurt1j8XEiIhNO0eB3917gUoLTNC8C97j7FjNbbWYXRqouBda6+4DTQGY2l+CI4bFSNXooB3oPALDv4L6x+DgRkQmn6OWcAO6+HlifU3ZtzutVBd77GgUGg8uhOxNcNKTgFxHJL3F37maDf+/BvRVuiYjI+JS44O/q6QLU4xcRKSRxwX+4x39APX4RkXwSF/xdma7DvzN9ea8eFRFJtcQFf7bHD/D2obcr2BIRkfEp0cGv0z0iIoMlLvizg7ugAV4RkXwSF/wDevy6pFNEZJDEBX92cBfU4xcRySdxwR/t8Sv4RUQGS3Twa3BXRGSwxAV/V08XU+qnUFtVqx6/iEgesR7SNpF0Z7ppqm2ivrpeg7siInkkLvi7Ml001TVRZVXq8YuI5JG44O/OdDOpdhINNQ0KfhGRPBIb/FPqp+hUj4hIHskb3M100VTbREtDi3r8IiJ5JC74sz3+qQ1TdTmniEgeiQv+rp5gcDfb48/5CmARkdRLXPB3Z7qZVBP0+DP9mcNfvi4iIoFkBn94qgd0966ISK7EBX/2Ov6WxhZAz+sREcmVqODv7e+lp69nYI9fl3SKiAyQqODPPqCtqbbpcPCrxy8iMlAig39S7SRaGnSqR0Qkn8QGvwZ3RUTyixX8Zna+mW01szYzW5Fn/hozezr8ecnM9kXmzTGzh8zsRTN7wczmlq75A2W/b7epTqd6REQKKfqsHjOrBm4BzgPagY1mts7dX8jWcfcrIvUvAxZGFvFPwA3u/iszmwz0l6rxuaI9/trqWppqmxT8IiI54vT4zwDa3H2bu/cAa4GLhqi/FLgbwMwWADXu/isAd+909+4h3jsq2e/bbaptAgge26CrekREBogT/LOA7ZHX7WHZIGZ2HDAPeCQsOgnYZ2Y/N7M/mNm3wiOIsoj2+AFaGvWgNhGRXHGC3/KUFXoAzhLgXnfvC1/XAB8ErgROB44HLhn0AWbLzWyTmW3q6OiI0aT8coNfPX4RkcHiPI+/HTg28no2sLNA3SXAl3Pe+wd33wZgZvcBi4E7om9y99uB2wEWLVo04qeqRQd3AVoaWti+f/tQbxm17kw3v339tzTUNNBU10RTbROT6yYfnq6tri3r54uIDFec4N8InGhm84AdBOH+n3Mrmdl8oAV4Iue9LWbW6u4dwEeBTaNudQH5evzP7X6uXB8HwJUPXcn3N32/4Py66jqaaptoqgt3CPmmayfn3WkUmm6qa2JS7SSqLFFX44rIGCka/O7ea2aXAg8C1cAP3X2Lma0GNrn7urDqUmCtR56D7O59ZnYl8LCZGbAZ+MeSr0Uo7+BuGa/jf3Xvq/zg9z9g2fuWccmpl9DV00VXpovOns7C05kuunq62NW56/B0trynr2dYnz+pdlLxnUWxnU6e6brqOoLNJSJJFOurF919PbA+p+zanNerCrz3V8ApI2zfsGR7/I21jUBwqmf/of30e39Zesff/H/fpMqquOncm5g1Je9497Bk+jKHdwZFdyB5dibZ6T3dewbV6/f4V9FWW3WsI4/h7liaapuorirb2L6IxJSo79ztznTTUNNwOOSnNkzFcd4++Pbhp3WWyktvvsSdz9zJ5X9xeUlCH6C2upap1VMP33xWKu7Owd6Dg3YQsaYjO5C9B/eyff/2ATudg70Hh9WWhpqGshylNNQ06ChFJKZEBX9XT9fh0zzAgEczlzr4Vz26ioaaBlacPehG5nHHzGisbaSxtpHpk6aXdNl9/X10Z7qHtzPp6aIzM/CoZcf+HQN2Mp09nfQdvjisuCqrGnBkUXBnEXMcRQP0kmSJCv7u3u7DA7tA2R7b8Nyu51j7/FpWnL2CY5qOKemyJ5rqqmqa65tprm8u6XLdnZ6+nlEfpbzT8w5vdL4xYCeTPSUYVzkG6CfXTaaxtlED9FIRiQr+7PftZpXrmfzXPXodzfXNXPmBK0u6XDnCzKivqae+pp5pjdNKuux+7+dA5sCoj1KyA/TRsZSRDNAXPaU1gqMUDdDLUBIV/NmvXcwqx6OZN+/czC/++Auu//D1JQ8kGRtVVhUEZKSTUCrRAfo4O5OhBuhzy73gfZOD1VTVFBxg1wC9JDr4y/Fo5pUbVjKtcRpfWfyVki1TkmMsBuhjX+2VZ8fy1oG3Dg/QZ8s1QJ8+iQr+rkzX4V4+UPLv3f3t67/lgbYHuOncm5hSP6UkyxSJo9wD9KW4jHjH/h2DykczQF/Ky4g1QD9QooK/O9PNrOYjl1ZOrptMlVWVLPhXbljJu5rexZdP/3LxyiITRHVVNVPqp5S8MxMdoB/NUcr+Q/sPD9Bny0c6QD/U3fBpGqBPVPDnDu5WWRVH1R9VksHdR159hA2vbeAfzv+HspwbFkmacg/Qd2e6R32Usqtz16DyTH9mWG0p5QD91IapzGieUdJ/q3wSFfzdmW4m1UwaUNba1Mrurt2jWq67s3LDSmZPmc3y9y8f1bJEZPSqrIrJdZOZXDe55Mse7gB9vsuIuzJddHR3DNoBFRugP2PWGTz1t0+VfJ1yJS/4awcG/8zmmex8p9DDROP5ZdsveXz749z6l7fSUNMwqmWJyPhW7gH6oXYgYzV2mJjgd3e6Ml2DTsPMap7F49sfH9Vyv7HhG8ybOo/PL/z8aJspIikVHaBvpbWibUlM8Pf09dDv/QV7/O4+okvF7vvjffz+jd/z44t+TF11XamaKyJSMRNvOLqA3EcyZ81snsmhvkMjGuDt6+9j5YaVzD96PstOWVaSdoqIVFpievwAF//7iznp6JMGlM1sngnAznd2DvvKgnu23MOWji2s/Y9rqalK1D+ViKRYYtJsWuM0fv6Znw8qjwb/e495b+zl9fb3ct2j1/G+Y97HX7/nr0vWThGRSktM8BeSDf4d+3cM6313PXMXL7/1Mvd95r4JeYOGiEghiU+0GZODmyGGc0lnT18P1z92PYtmLuLC+ReWq2kiIhWR+B5/Y20jLQ0twwr+O35/B396+0/c9qnb9NAoEUmcxPf4IbykszNe8B/IHODvfv13nD3nbD7+7z5e5paJiIy9xPf4AWZNmRW7x3/rplvZ+c5OfvIffqLevogkUnp6/DGCv7Onkxt/cyPnHn8u58w9ZwxaJiIy9tIR/JNn8sY7b9Dv/UPW++5T36Wju4NvfuSbY9QyEZGxl47gb55Jn/fR0dVRsM6+g/u4+fGb+dRJn2Lx7MVj2DoRkbGVmuAH2PFO4Wv51zyxhn0H97H6w6vHqlkiIhWRquAvdJ5/T/ce1jy5hk8v+DQLZywcy6aJiIy5WMFvZueb2VYzazOzFXnmrzGzp8Ofl8xsX2ReX2TeulI2Pq5iwf+t336Lzp5Orv/w9WPZLBGRiih6OaeZVQO3AOcB7cBGM1vn7i9k67j7FZH6lwHRbvMBdz+1dE0evndPfjeQP/j/3Plnvvu777LslGUsaF0w1k0TERlzcXr8ZwBt7r7N3XuAtcBFQ9RfCtxdisaVSm11Lcc0HZM3+G/89Y309PVw3TnXVaBlIiJjL07wzwK2R163h2WDmNlxwDzgkUhxg5ltMrMnzeyvRtzSUZrVPPgmru1vb+fWzbdyyamXcMK0EyrUMhGRsRXnzt18t68W+sbgJcC97t4XKZvj7jvN7HjgETN7zt1fGfABZsuB5QBz5syJ0aThy3cT1w2/viH4IvUPrSzLZ4qIjEdxevztwLGR17OBQrfBLiHnNI+77wx/bwMeZeD5/2yd2919kbsvam0tz3dR5gb/tr3buOMPd7D8/cs5bupxZflMEZHxKE7wbwRONLN5ZlZHEO6Drs4xs/lAC/BEpKzFzOrD6enAWcALue8dCzObZ7K7azeZvgwAqx9bTU1VDV//4Ncr0RwRkYopeqrH3XvN7FLgQaAa+KG7bzGz1cAmd8/uBJYCa909ehroZOA2M+sn2Mn8j+jVQGNpZvNMHOfPnX+mK9PFXc/exRWLrzh8qaeISFrEejqnu68H1ueUXZvzelWe9z0OvG8U7SuZ6LX8a55cQ2NNI1eddVWFWyUiMvZScecuHAn+X7b9kp9u+SlfWfwVWpvKM54gIjKepS74b/zNjRxVfxRfPfOrFW6RiEhlpCb4Wye1Um3VHOo7xJUfuJKWxpZKN0lEpCJS8Q1cANVV1cxonsGBzAEu/4vLK90cEZGKSU3wA6w6ZxXvmvwumuubK90UEZGKSVXwf/G0L1a6CSIiFZeac/wiIhJQ8IuIpIyCX0QkZRT8IiIpo+AXEUkZBb+ISMoo+EVEUkbBLyKSMjbw8fmVZ2YdwJ9G+PbpwJ4SNmci0Dqng9Y5HUazzse5e6xHDo+74B8NM9vk7osq3Y6xpHVOB61zOozVOutUj4hIyij4RURSJmnBf3ulG1ABWud00Dqnw5isc6LO8YuISHFJ6/GLiEgRiQh+MzvfzLaaWZuZrah0e4bLzI41sw1m9qKZbTGzy8PyaWb2KzN7OfzdEpabmX0nXN9nzey0yLI+F9Z/2cw+Fyl/v5k9F77nO2ZmY7+mg5lZtZn9wczuD1/PM7Onwvb/1MzqwvL68HVbOH9uZBlXh+VbzewTkfJx93dhZlPN7F4z+2O4vc9M+nY2syvCv+vnzexuM2tI2nY2sx+a2W4zez5SVvbtWugzinL3Cf0DVAOvAMcDdcAzwIJKt2uY6zADOC2cbgZeAhYANwMrwvIVwE3h9AXAA4ABi4GnwvJpwLbwd0s43RLO+x1wZvieB4BPVnq9w3b9d+AnwP3h63uAJeH0rcCXwun/BtwaTi8BfhpOLwi3eT0wL/xbqB6vfxfAncDfhtN1wNQkb2dgFvAq0BjZvpckbTsDHwJOA56PlJV9uxb6jKLtrfR/hBL8g58JPBh5fTVwdaXbNcp1+j/AecBWYEZYNgPYGk7fBiyN1N8azl8K3BYpvy0smwH8MVI+oF4F13M28DDwUeD+8I96D1CTu22BB4Ezw+masJ7lbu9svfH4dwFMCUPQcsoTu50Jgn97GGY14Xb+RBK3MzCXgcFf9u1a6DOK/SThVE/2DyurPSybkMJD24XAU8C73P0NgPD3MWG1Qus8VHl7nvJK+3vga0B/+PpoYJ+794avo+08vG7h/LfD+sP9t6ik44EO4Efh6a0fmFkTCd7O7r4D+J/A68AbBNttM8nezlljsV0LfcaQkhD8+c5hTshLlcxsMvC/ga+4+/6hquYp8xGUV4yZfQrY7e6bo8V5qnqReRNmnQl6sKcB33f3hUAXweF5IRN+ncNzzhcRnJ6ZCTQBn8xTNUnbuZiKr2MSgr8dODbyejaws0JtGTEzqyUI/X9x95+HxbvMbEY4fwawOywvtM5Dlc/OU15JZwEXmtlrwFqC0z1/D0w1s5qwTrSdh9ctnH8U8BbD/7eopHag3d2fCl/fS7AjSPJ2Phd41d073D0D/Bz4AMnezlljsV0LfcaQkhD8G4ETw6sE6ggGhNZVuE3DEo7Q3wG86O7fjsxaB2RH9j9HcO4/W/7Z8OqAxcDb4WHeg8DHzawl7Gl9nOD85xvAO2a2OPysz0aWVRHufrW7z3b3uQTb7BF3XwZsAD4dVstd5+y/xafD+h6WLwmvBpkHnEgwEDbu/i7c/c/AdjObHxZ9DHiBBG9nglM8i81sUtim7DondjtHjMV2LfQZQ6vkwE8JB1UuILgS5hXgmkq3ZwTtP5vg0O1Z4Onw5wKCc5sPAy+Hv6eF9Q24JVzf54BFkWV9AWgLfz4fKV8EPB++53+RM8BY4fX/MEeu6jme4D90G/AzoD4sbwhft4Xzj4+8/5pwvbYSuYplPP5dAKcCm8JtfR/B1RuJ3s7A9cAfw3bdRXBlTqK2M3A3wRhGhqCH/sWx2K6FPqPYj+7cFRFJmSSc6hERkWFQ8IuIpIyCX0QkZRT8IiIpo+AXEUkZBb+ISMoo+EVEUkbBLyKSMv8fvV4sEGTnP6QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Построение графиков зависимости точности от размера последовательности\n",
    "#Зеленый график - для Витерби в логарифмическом пространстве\n",
    "#Синий график - для Витерби с масштабированием\n",
    "plt.plot(results_log[0], results_log[1], color = 'green')\n",
    "plt.plot(results_sc[0], results_sc[1], color = 'blue')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Задание начальных условий для работы forward/backward/aposterior/baum-welch\n",
    "L = 300\n",
    "hidden=np.array(hidden_seq(L))\n",
    "observed = obs_seq(hidden, L)\n",
    "viterbi_path = viterbi_scaling(observed,pi,transit, emis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Алгоритм прямого распространения\n",
    "def Forward(obs_seq, transit, emis, pi):\n",
    "    L = len(obs_seq) \n",
    "    M = len(emis)\n",
    "    forw = np.zeros((M, L))\n",
    "    scaling = np.zeros(L)\n",
    "    forw[:,0] = pi*emis[:, obs_seq[0]-1]\n",
    "    scaling[0] = max(forw[:, 0])\n",
    "    forw[:, 0] /= scaling[0]\n",
    "    for i in range(1, L):\n",
    "        for j in range(M):\n",
    "            forw[j,i] = emis[j, obs_seq[i]]*sum(forw[:, i-1]*transit[:, j])\n",
    "        scaling[i] = max(forw[:,i])\n",
    "        forw[:,i] /= scaling[i]\n",
    "    proba = sum(forw[:,-1]) \n",
    "    return forw, proba, scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "forw, prob_forw, m_f = Forward(observed, transit, emis, pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Алгоритм обратного распространения\n",
    "def Backward(obs_seq, A, E, X):\n",
    "    L = len(obs_seq) \n",
    "    M = len(E)\n",
    "    back = np.zeros((M, L))\n",
    "    scaling = np.zeros(L)\n",
    "    back[:,-1] = np.array([1,1])\n",
    "    scaling[-1] = 1\n",
    "    for i in range(L-2, -1, -1):\n",
    "        for j in range(M):\n",
    "            back[j,i] = sum(back[:,i+1]*A[j]*E[:,obs_seq[i+1]])\n",
    "        scaling[i] = max(back[:,i])\n",
    "        back[:,i] /= scaling[i]\n",
    "    proba = sum(back[:,0]*X*E[:,obs_seq[0]-1])\n",
    "    return back, proba, scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "back, prob_back, m_b = Backward(observed, transit, emis, pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_posterior(F, B, P, m_f, m_b):\n",
    "    L = len(F)\n",
    "    Result = np.zeros(L)\n",
    "    for i in range(L):\n",
    "        Result[i] = np.log(F[i]) + sum(np.log(m_f[:i+1])) + np.log(B[i]) + sum(np.log(m_b[i:])) - np.log(P) - sum(np.log(m_f))\n",
    "    return np.exp(Result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Вычисление апостериорных вероятностей\n",
    "post_prob = find_posterior(forw[0],back[0], prob_forw, m_f, m_b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ниже приводится график, построенный по следующим принципам:\n",
    "    \n",
    "    1) Желтые промежутки - использование честной кости;\n",
    "    \n",
    "    2) Красные промежутки - использование нечестной кости;\n",
    "    \n",
    "        а) На первом рисунке - действительное использование;\n",
    "        \n",
    "        б) На втором рисунке - предсказанное алгоритмом Витерби;\n",
    "        \n",
    "    3) Черный график - постериорная вероятность использования честной кости;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2f05e8aa898>]"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztvXmcI1d57/19Smr13rN29+ybZ2zPeME2c41hwDZm8xCwzUtI7BsMN0CcsCUGcgOEXN68JDcBEgLvTQyEACEOL2AIBkxiIMELxAY7Hi/YHq/jGXumZ+uenp7eF0l13j+qjqRWl6QqSa2q0pzv59MfdatLpVPbr576nec8R5RSGAwGg6G5sMJugMFgMBjqjxF3g8FgaEKMuBsMBkMTYsTdYDAYmhAj7gaDwdCEGHE3GAyGJsSIu8FgMDQhRtwNBoOhCTHibjAYDE1IMqwvXrlypdq0aVMVn3wamPb/r/Z2OIvSn6n1q88EOnx8br3P5YI0M/fd7e4bFT5cav1+tsHPehaTdncbpz2+OGj760G990G57SumUcfL7/f4vi58nqfl1h1k2+t1vgfC7za24wpTYB588METSqneSstJWOUHdu7cqfbs2VPFJy8HHvH/rwsugLsp/Zlav/pu4AIfn/ucz+WCNDP33XrFFT5cav259fik1HoWkwvcBj7i8cV3E6z99eBy6rsPym1fMXfTmOPl93u81u/5WZ/nabl1+21TqXYFXUdg/G7jBW5DgiMiDyqldlZaztgyBoPB0IRUFHcR+aqIDIrI4yX+LyLyf0Rkn4g8KiIX1b+ZBoPBYAiCn8j9a8CVZf6/G9jm/twAfKH2ZhkMBoOhFiqKu1Lq58DJMotcDdysHO4DlorI6no10GAwGAzBqYfnvhY4VPD3gPuewWAwGEKiHuIuHu95puCIyA0iskdE9gwNDdXhqw0Gg8HgRT3EfQAni1uzDjjitaBS6ktKqZ1KqZ29vRXTNA0Gg8FQJfUYxHQb8D4R+RbwEmBUKXW0DuuNPFNKcc8vM+y/z+bkSUVLi7BunXDBBQm2b0+E3TyDwXAaU1HcReSbOMMBVorIAPB/Ay0ASqkvArcDrwf2AVPAby9WY6OCUoq/mpvjL2dmOPUe72UuuSTB5z/fzoUXGpE3GKLK2Jji61+f48EHs/T3W7ztbS2cfXZzXLMVxV0pdV2F/yvgvXVrUQy4YXqaL6fTvCGZ5L2fS3H+mxIsXy7MzcHBgzZ33pnhL/9ylssum+DOOzvZGV6VB4PBUIL9+212757kmWds+vqEkycVn/nMLF/7WjvXXZcKu3k1Y0aoBuR76TRfTqf5o1SK2zo6uHJXC2vWWLS1CT09wrnnJvj9329lz54uli8X3vKWKebscEo8nM787GcZ/vEf5xgctMNuiiGCZLOK3/iNKYaGFHfc0cnx4z0MDHRzySUJrr9+mr17s2E3sWaMuAfAVooPTE9zvmXx521tiHglCjmsXWvxhS+08/zzin88NtfAVhr+9m9nufzySd7xjmle8pIJDhwwAm+Yz1e/mubBB7PcdFMbV1zhPFn391t897sd9PQIf/AHja6MV3+MuAfgP0ezvKAUH21tpaWMsGuuvDLJJZck+JsBI+6NYt++LH/wBzNcfXWSO+/sZHhY8Ud/FP8LtVlQSnHvvRkOH67/Dfezn51l3box3vCGSaanyz8tf+ELs7z4xQmuvbZl3vsrV1p89KOt3HFHlmeeiXf0bsQ9AN8YTNMBvLGlpeKyACLCf//vLTwzbfPcoXifKHHhi1+cI5GAz3++nVe+Msl73tPKrbdmTPQeET796UO8/OWTbNgwzg9/mK7bevfsyfDBD86wapXF7bdnePe7S9/Q9w/YPPywzXXXtXg+fWvB/+53y7fv6FGbm2+eI52Opu1qxN0nSim+dyLNVS0tdPqI2jW7dzuPfD+6N7NYTTO4zM4qvvrVNG96k9MPAvD+96cQga98xTw9NZL9AzYf/OD0vJvqY49N8Md/fIBrrkmyfbvFjTfOMDNTH2H8yEdm6OsT7ryzk/e+N8U3v5nmZNr7hv7dOxzRfvObvYO09estLr44wa23lr5mH3ggw7Zt47z97dO84x3TFJZOz2adp5NsNlzRN+Luk4MHFUNpxaWJYGlSW7cm2NpuGXFvAI88kmVkRPGbv5m/aNeutXjxixP8/Odm/9fKnj0ZPv3p2YqR6knb5nXvmeSzn53jRS8a5/nnHZG95RZnVPo//EM7n/lMG/v329x6a+3R+6lTirvuyvK7v5uip0d45ztTzM3BNwe9133nf2U491yLTZtKy98b35hkz54sp055b+snPzlLa6tw440pvv71NHfc4TyZj45mePnLH+blL5/kD/9wpuZtqwUj7j556CHn4F0UUNwBLluS4IEm6H2POvfd5+zjSy6Zf4x27UrwwANZ5uai+fgcB555JstrXzvFhz88wzXXTFFukp+vpNPsO2Tzta+1MzEBX/ua89T0/e+f4NJLl7BypcWrX51k2TLhzjtrv+nedVcG24bXvtZ5Sr7gggTnnWfxL0Pe4v7os1kuuqj8daz//9hjC6/bgQGbH/wgwzvf2cL//t9tdHTkLZyvfe0Y9903zu7dST73uTnuuiu8oMKIu08eeihLAji/CnHf3pFgaEQxPGx838Xk/vuzrF0rrF07/7TetSvJzEz+Bl2JuTnF1VdPctVVk+aG4PIXfzFLNqt4//tT3H57hscec87l0VHFl4/O8a25uZzg/zCd5oKzLN7+9hRXXJHg5pvnePbZLHv3TnHNNSsBSCSEyy5L1EXcf/rTDF1d8JKX5K/NXbuSPDyRXXATOmHbHBlSnH9++etY///RRxeeM7fdliabhXe9K0VHh/D61yf53vfSZLOKf/7n41x4YRe33tpBezt873v161cIihF3nzz0UJbtnRbtAfx2zfYOZzc/+aQR98Xk/vuz8y5wza5dznv33utP3D/0oRluuy3DD3+Y4UMfCvfROgrYtuLHP87wa7/Wwsc+1opIXrQ+8YkZfueZaa6bnuaX2SzDts292SxvvNSxxq6/PsWBA4ovftGJ3q+4YmluvVdckeTAAZWzbarl7rszvOIVSVpa8tfmhRcmGM3CgSJxf8x2vuv888tL39q1wrJlkruJFfLgg1lWrhS2bXPW8aY3tXD8uOKWW9I8+OAE11/fT1ub8MpXJvnRj/I3r1OnVNknnnpjxN0nDz+c5cKu6oYlb+9wPvfUU0bcF4uREcX+/TYXX7xwNPCqVRb9/cKTT1YW92zWGY7+1re28Na3tvD1r8+F3jEWNo8+anP8uOJ1r0vS32+xa1eCW29NMzur+Kd/SrN7eZIe4Ka5OX6ayWADb7jUOQ4ve5lz7t9ySxrLgjPPzM9Mfam7zC9/WX30nskonn3W5oIL5l+bF17oSNsj2fnH/FH37xe9qPy1LCKcd57lGbk/+GCWF784kcu02bnTWdeXv+zcwHbvXu6+Jtm3z2bfvixPPJGlv3+Miy+e5NChxuiAEXcfzMwojh5VnNle3e7a0Ca0teJLXAzVobMyzjrL+xht3mz5ihAfecTm1Cnnwty9O8mpU857pzM/+YkjvtrT3r07yaOP2tx8c5rhYcUH1qb47VSK76TT/CKbxQJedKYjeFu2WLS3w+HDijPOaKe1NX98zjzT+f2556rfv4cOKdJp2Lp1/nE/77wECeBhD3HvWy709VW+ls8/P8Fjj2WxC0aYz8wo9u6153n2W7ZYpFLws59laW0Vtm5tB8gNjvrlL7N88pOzJJPwq19l+bu/m612cwNhxN0HR444B3dta3W7KyHCWRstE7kvIgcPOvt2wwbvY7Rpk8WBA5Uj8DvucITsiiuSuYtTv3e68sgjWTZtElavdvbtjh2OsH31q3N0dsKrliW5MpkkDfxLOs0my6I15US1iYRwzjnO8tu3d8xbb3u7sGaNsH9/9dfFvn2OeJ9xxvzj3tYmbO+0FkTuz9s2W9f7u463bbOYmIDh4fx58/jjWTIZ5ol7MimceaaFbcP27Z0kk862b91qYVmOHfiNb6S54YYUl12W5LbbGnM+GXH3wcCAc/KtSwX32zVnbUrw7LNG3BeLF17Q4u59jDZvtjh40K5osdx1V4YdOyxWrXJ+duywQs14iALPPWezbVtezM4+25GN++7Lcs45CSyRXKLBEaU40yqOop2/i8UdnKi3lsh93z7ns8WRO8AZbRYv2PPXfVQp1vT6u471WAkd3IFjUQFccMH879Mlvs89N7+NqZSwcaPwne84HbBvfnMLV12V5KmnbPbtW/xR00bcfaCHSq+rMnIHWL1SOHbMiPticfCgTXs7rFjhfeFu2mSRycy/UDU33TTL97/vdBA+9dT8NLmLLkrw1FPxtNMeeCDD0aO1n3P79tnzIuMtWyx00ti55zrvrxZhuetBLxR378hdr6u2yN057qtXLzzua1otDhd1YB6xbVav9Hcd63UW7kMd6BU/IW7f7vx97rmd897fujXByZNOG847L8Eb3+h0NP/wh8O+2lALRtxLkMkoXvGKCV7zmsncY3kt4t6/XBgfp2LNC0N1HDxos2GDVbKY26ZNzvvFvvtDD2V53/tmuPbaKR57LMuhQ2re4JZNmywOHVJkMvE6bg88kOFlL5vk6qvL56RXYmREMTKi5ol7KiVs2eL8rS0XEeE8V9TPKhL3yy5L0toKl1zSs2D9Z5xhcfiwqnqkqr7xWNbC4742JZxUihl3+6eUYhQCR+5Hj+bbdviwordXSBU9xe/YofdHsbg772/aJCxZImzaZPGjH3Xwrnet8reBNWDEvQT33JPlnnuy/PSnGb7ylTTd3dCdrN6W6V/h7OrBwXiJRFw4eFCV9NuBnGAXi/uf/MkMK1YInZ3CO94xTTbLAnHPZmFgIF7H7Z3vnCaRgAceyPKvw9XbSs895zy1FNseuuNaR+4A57nh/JlFY0EuuijB5GQPZ53lHbkrtfC4+GXfPtvTkgEncgfHKgI46lo0QSP3I0fybTt82Gbt2oU68IY3tPBnf9bKa16zbN77Ol2yMDvnyitb6O5e/DkejLgXUOjH3nprmrY2eMMbnINQPDAmKH3LnRPi+HFjzSwGL7xgl/TbATZudI5fYa2TbFZx990Z3vrWFl7zGme4OTj+vKbUTSHKnDxp89hjNv/rf7Wydq1w8/Hq6+poT7u4w1KLu47cAV6aSNACnGMtvFYSCe9jo9dbre9+9KgqeW2udaPrI66oa5Ff0+cvSGtrc3LdCyP3I0dsz+/r7BT+5E/a5mUDQf6mWGnQ1GJgxN3l6FGb3t5xPv/5WadI2PfSXHllMueRTU3VFrn1r9DiHo8IcHDQZnIyHm2dmVUcP65yAu5Fa6uwcqVw7Fh+m/bvt5medi68wukQCyN3LfRxqir58MNOWy++OMHOnQn2TlXfdi262obRvOtdKf7iL1pZsyYvlNe2tPBcdzerPcS9FFooCwXUL5mMYxmtXOkt1jpyP6zF3X1d0+u/fatXS1HkruZtcyUuuCBBaytcfnnjZ2Mz4u7yqU/NMjKi+NjHZnjqKZuBAcXllyd55Sudi/7gwRrFfbmzq+Mg7uPjivPPn+C3fmsq7Kb44pibqqY90lIsXy65zi2Axx93LtrzzrNy2Q+WBevX5y/edesEy4pX5P7ww84TyIUXJjjnnATPTtvMVem7Hzyo6OtzbKtCzj47wUc/On/CGkuE9QGEHfId4IXphn7Rx7KUuK9NFdky7uvqEst7sXp1PnKfm1MMDpZ+UvBi/XqLkZGeXFptIzHijiNmf//3c1x+eYJTp+DjH3eGnG/ZYrF1q8XVVyf5zncW+oVB0LZMHKZ9+9znZjl+XPGDH2TYsyf6aYAjY87Ft3x5+Yt2obhnEXHytvUIx3XrZN4w9lRKWLtWYifu69YJK1c6qZwZBc/Y1bX/xAm7pHjWg44OSKUcKyko+oZQqn1Lk9DG/Mi9FVjW43971qyxctkyWuSDWrTt7Yu3/8phxB1n5OjMDNx4o/OY+eMfO4K2ebOTffH973fy67/ub4KOUrS1Cj098Yjcv/CFOa64IsHSpfB3fxf9Ouha3JctCybujz1ms2WLRWen0N9vsXq1eJaB3bzZipUt89BD+XRO7YnvzVaXzjk8rEqml9YDEWHFCqkqcj9xQou7t4yJCGssKxe5H1GK1SJlp8csRkfuSqlcSrRXh2oUMeIOPP10fuj6jh0JJiac9ws71upBf78VeXEfGXFKLeze3cJ/+29J9u6NvqhpcV+6tLK4j4zk9//evVnOOSd/jP/qr9r48IdbF3xu1SorNllOSikOHLBzHZ5nn21hAXurjNwXW9yBmsW9XPvWiOS89mGl6A1oG/X3W8zNwehofoxEJfsvKsSjlYvMU0/ZJJNOz70ejODlM9ZKX59E3pZ5+mknwjv7bIuzzrJ4+umFZVOjxqlxf5H7smUy7/H/8GF7XqT+W7+V4vWvX/iE5nwu2vtAMzoKs7PODQmcjI9NbVbVtkyUxb2SLQOwTIRT7vk7phRLAlZ1XbLEWX5sTOUy3VatMpF7bHj66Sxbtli0tEhO3OsdtQP090vkI3dd/+assxxxHx9nXoZJFBnxKe7LlwunTjkpkLOzitFR54ZbCR3xR/0mB+RGQRcK0KqUMFiFuCulIi3uJ07Yuc+XYokIY+5xG1WKhcOoytPTkxf30VFnPVrwo44RdxxB04+xuihScepXPSi2BaLIU0/ZtLQ4Nze9T7RtFVVGxhSWBV1d5ZfTHa6nTimGhpzj0OsjLW7ZMiGTgcnJmpu66OjgQUfuAL0twlAVN6bJaZibgxUrFlcmarFl2tuho6O8uI+6v48pRU/AyL3HvRuMjSnGxpzO37Y2I+6xIJt16kHrYkg6cl8Mce/qEiYmoi3uTz+dZetWy610l8i9F2VOjSuWLhXPIeiFaHE/eVLlPHQ/kbt+IoiDNaMj9/7+/HZVK+7Dpyp72vVgxQqL4eHgT0bDw6Vz3DU9IowqZ93Vifv8yD0uUTsYcefIEcXcnFPgB6Cvz+Lmm9v5vd9L1f27uruFiQkW/fH+oYfm16AOwlNP5W9069cL7e3xiNwrWTJQLO7ONvm1ZYDIP3VBYeReKO4WJ5TCDiqeo40R9+XLnSej8fFgnztxorK4LxEhC0wCY+7fQSgW954AaZRhc9qLu+5x7y0oJnT99SnWrVucyF0pmFrEsUH33ZfhxS+e4JOfrG5CgEOH7NxIT8tyUgOjnuM9MqYqZspAPgIfGakuco+DuB87pkgm5/c/9KUcgTsVWNwre9r1oNqBTE5/QPnrdIn7esS2URA4cs93qDoCv2RJhQ9EiNNe3PUJtdgnMEB3t/M6Pr54IqHntvyzP5sNLMqzs4rJyfnZBytXVueHNpJT48Ejd+25+5mRJ262TH//fIuq1x2UFdSaaZwto8U92PnqJ3LXYn7I3fZqbZnRUWPLxA59QjVG3J3vWExxv+22DNu2WczMwL//e7DRpV43umo7uxrJSEBx/8AHZvjGN+Zobc3fcP18Lg6R+/Hjap7fDjWIe4NsGb3+oDfPkZHKx13bMANutlDQbBndSW889xiSH+W2+Aetq8v5Dj1Iqt4cOGDz1FM27353ilQqPwWZX/J5w/nTIhbiPqZYurTycloIBgcVDz1k09vrb7RivGwZe16mDDieO1QfuVcq61Ar1doyk5OKzs7yy2hxP+SKe1DP3bKE7m7juccSfUIt9gkMix+566yWiy9OVDV9mdeIPy3uUc3xVkr5tmWSRfX4/Q5S6+qCRCIetoxn5J6qPnLv6WFerZ3FIH9d+P+MbSumpysfw1ptGXCsGScVsgkjdxG5UkSeFpF9IvIRj/9vEJG7RORhEXlURF5f/6YuDsPDjTmBIf+It1jifvhwvrDR1q1Wrha3X7xtGYt0evGeNmrhhG2zdnycuXTlAUyae+7p5Oabndnp/WYBiTh1veMQuXsNOtK2TNCBTGOTjREzLdBBymrrpIRK4l4cuVcr7qOjTp57U4m7iCSAm4DdwA7gOhHZUbTYnwDfVkpdCFwLfL7eDdVMT2f5t39L1219TqdMYx5gdISyWLnuurDRmjWSE/cgEbdX/0MtJVkXm7uz2VwZ13IDWQrZtSvJW97ilBjwkymjWb5c+M//zPCf/xndKpmZjGJmJn+eaVotoYfgkfukj8i4HnS4BVeDibuzbFBbplpxP3JEYds0nS1zMbBPKbVfKTUHfAu4umgZRb6vYglwpH5NnM+f//lBrrpqinvvrc9F1ojh1ZpqHj+DMDBg09fnzO+4davF1FSw0gGlOlQL/xcl7s0458D7fjPFNdf4r9rZ1ib89Ked3HNPBWUoIJWCvXtt3vzmqarHECw2+umqWNwBei0rsLhPzaiGiHsq5dheQUYA62Ur3dR1f3mt4q4nxm6qyB1YCxwq+HvAfa+QPwXeKiIDwO3A++vSOg8+/OH1bNwovO1t03WZtLiR4q47VBfTltHlSPX0XkGsmeFhp4OqcHh1tWlqjeAX2SyXJhL87Ufay86f6sWrXpVk2zb/U59pC2doSPHEE9HbF5A/r7zEvQcYDxy5q1xUvZiICB0dwSJ3PUtYpZtPQoQuyJUgCJotA04JgkOHdF2ZKlYQEn6uCK+9V3wUrgO+ppRaB7we+GcRWbBuEblBRPaIyJ6hoaHgrQV6epL88R+3sX+/XdPgmuFhm4svnmDPnmwDxd15XUxbRk8koAufBdlHJ04svNFFNXKfUoqHsllelmzMDDff/nYHH/qQM2r5rruiac3o88qrxk6nCJOBI/fG2DLgfE+QaR3z4l55WW3NdOKIfVAKrZhmi9wHgPUFf69joe3yTuDbAEqpX+JMgLKyeEVKqS8ppXYqpXb29vZW12LgzDOdZu/fX72433JLmgcecLJLGpEGCU6nbWvr4kXuAwMqN7K2mtxsr6eYqIr7M7ZNBnhxojETD19zTQt//dftbNwo3H13NMW9XOTeJULQPvFGRe7g2CtBRm5rW8bPzUdbMdVYMjBf0JvNc38A2CYim0UkhdNhelvRMgeBVwGIyHYcca8uNPeBLupVy+w43/hGvlNW2yWNQNeXqTczMzbDw3lbRg/HDy7u808JfZOImrgfdz3U1VVesNWya1eSBx+MZiG1cp57J1QRuTfGcwdqsGUqL1uruDdt5K6UygDvA34CPImTFbNXRD4hIle5i30I+B0R+RXwTeB/qEVMjF6zRkilqo/cjxyxuffeLOef72z+Yo4YLaa7e3G+78gRp5aMtmWSSWfwRVBxL36KaWlxpgeMmrgPuqdXX4PFfe1a4dixaOb96/PKK1ipxpaZnPafhVQrji3jf3m/njvAJnf2pR0BZ2HSFNYtipO4+zIslVK343SUFr738YLfnwB21bdppbEsYfNmq2px14N7PvnJNh591Obtb69tftQgdHXJooj70aPOXKerV+dPvqC52SdPKs/BXFEcpXrcFar+Ki/Yaunvt5idhdEs+BgU21C05+5VUqEaW8aJ3Gtvlx+CR+7Oqx9x/8f2dj7V1sa6KgOBt761hePHFdPTijVrmkzco8iWLdWLu84H37DBYvfuxgk7aFum/kI5Nub4wIVRRlBxn5hQno/0S5Yszg2pFgbdmex9lIapK7qU7vE5m6U0xu/3S9nInWC2jFKq4ZG7vi79EMSWaRdhYw1PeL29Fp/8ZFvVnw+L2JYfqE3c8yM5G40Tudd/vePjjg9cKM7LlgmnTvm7oJ2p5/DsQAuaydAIBpWiP+BM9vWgv985Z47PRWt/QH78RKkO1Rkg61PgZwGlvM+HxSB4h2qwwWunI7EV982bLU6dIjevYRCOHLFpbw8nZ3WxPPexMUfcCzt/gkTu+sLyuliCXniN4LhS9DXYkoF85H5sLnq57hMTCpESN2j3JujX1tZRfuNSIYPbMiLQFr+AumHEVtx1x181xZycwT5Ww6M+cGq1HDpkc/BgfcWhVOTuX9xLP+YG9UMbwaBt0x/C8dNFuY6no7U/wAkaOjvxnG5QH1a/1oy+lzcqMu7oCNahOjXlbGsY13BciK24615rv7ZDIUeO2LmUwUbzgQ+ksCy4/vr6hsLacy/sTGv6yD2Um7OQSMCxCNoypfpMwLFlACZ8ins+cq9P2ypRTeTeqKeKuBJbcdcdh9WI++HDijVrwtn0s89O8O53t/KLX2TJZusnEGNjWTo7IZGYH7lPTcGcDyEq52FGLXJXSjmeewi2jGUJfX3C8QjaMuPjpcdsBLVlwojcZ2bwfU04tdyNuJcj9uIe1HNXSrnD9MM7MTZvdiYEPnKkfoI5Pp5dELUFmWRCi3epDtUoifspIE3jc9w1/f0S0Q5VVXJmqaC2TOM9d+d7pqf9Le9noo7TndiKe7W2zMiIkxUSVuQOsGmT890vvFC/6G9sLLNgaLS+Ad56a5qZmfL7qVx97KjZMroueVjivmqVxbEIeu71tGWmVOmb/WIQtOyvsWUqE1txr9aWOXrUWb5wsE+j2bhxMcS9dOT+nvfM8A//MFf28+Ui944OmJujLlU468GoKzxLQxL3vj5hMJK2jKqbLaOXa1z5Abd9PhtoIvfKxFbce9zanUHFfWzMFYal4Yt7LVUtixkfz+T2iaZwdqJ77ilfD6VcfWz9XlSi95zwhCTu3d3CRATLy0xMeOe4QxXZMg2O3IPOxjQ5qUyOewViK+6JhFM7JajnXq5yXqNob3c65eoduRfbMoUVHn/xi/KVDPOpkN4dqoXLhI22Fjwq2zaEri5hso6d4fWinOdefbZM4wqHAb4HyxlbpjKxFXdwou+gkXu5UXyNZONGixdeWNwO1a1bLW6+uZ0///NWBgYUA7OlbyaVOlSdZerW3JrICU9IkXtnJ8wqSEeoeNi+fVlOniztuUc9WyboOWZsmcqchuIefuQOTqdqPW0Zrw5VEeH661O87nVOCaFfjpX2EvzZMtEQMy1QXSGJu/a1A4y5WVSUUrz2tZN0dQnveEfKcxl9zw6eLVOPFlYmaOQ+M+M8ARtKE2txX7JEarBlFqNF/lmzRjh2bHE7VDVnneUUuHphprrIPaq2TFiBW07cIxK5Dw8rDhxQfOxjrZxGy9AdAAAgAElEQVRzjncxs4QI7QTIlgGSSafkcyMIGrnPzSlaWxexQU1ArMW9lsi9kRN0eOFUWqQuky3Pztqk02pBh6qmq8uZhPhEmfS9qSlobZ0/CEoTuQ7VkG2Z3HSJERF3be/pFNtSdIoEqi3T0cC6LVrc/V7Pc3OQSpnIvRynobg7EUnYd/2eHkGpYDO+l0KXHig1BZiIsHKlcCJdOnIvl30Q9JF5sZnAqVXtbUAsPlqIFmFCrarQdYp0FlYpOvF/Q5pQis62xonn+vVCb69w552VpzBUSpFOOwGLoTSxFnfHlgn2GSejoPGlYovRQlxNVctivIqGFeOIe7nIvfR8mVHsUO0kvKJR+qkvOpG7np+g/P7wOxuTUoqfZbOcf2bj6tUnEsLVVyf5t39LMztbvo1pd4ZMI+7libW468g9yJRnzii+RWyUT7SFovPua8Gr3G8xlcW9dGZE5DpUlQrNkoG8LRMVz/3gQZuODhZMbl5MG06d9ko8ats8Z9u8+VWNncjmTW9qYXycitH7nDsez9gy5Ym1uC9ZItg2gSac1pF72OjyCfUQ9/Fx52LoKpP4vXKlVTFyL5UZEbkOVcLLlIFoRu4bNlQuYd0qQvlxyg7fTaexgGsub+xEbbt2Od/3xBPlEw10ITwTuZcn1uJeTfQ7Ph5+GiTko+x6iPuMmwVTLjWsUuRe3nOPpi0TFtHz3FVFvx2cPoo5Hzekx7JZzrYsepc3Vh56epy+sMHBSuLuvBpxL0+sxV2LWaWiWIVEJXLPe+61r0t7lOU6iVeuFE5mVMlp1srbMnqZaESqE0qFHLnn2xEFnMi98v5IifiyZY4rxaoQ9q+IW075ePn9amwZf8Ra3PUUWzMz/j9Tboh2I1mMyL2tTHbDypWCAkZKintpWyaREFpbIxS5E14aJEQrzz2TUQwN+ZufwG/kHlatfHDmqB0crCTuxpbxQ6zFXUfu09Pxi9zr6bnPumUFKkXuACdKXNyTk+WHmjvToIUvZhC+LZNKQYJo2DK6nIY+n8rRCr4890HbDq2ccl+fVBT3Wffxw4h7eWIt7tVF7tHw3PWjfT1SIfPiXj5yh9LiXi4VEqI1G1PYtoyI0JWIhi2jB+WVy5TSpHx0qE4rxTiEMj8t4Noyfj338K/jKBNzcdeeu7/llVJla143kkRC6Oqqty1TepmVK51DXTpyL19CNUoTdoRtywB0JSQS4q7PH1/iDsxWaPOg+/++kG2ZcunN2pYJeyBi1Im1uLe3O69+O1RnZyGTCb+ujKanR+pky+iTvXLkPuRx0Ty5P8voKJx9dunTIXKRe8ht6Er4H8q/mASplZSisi1z3J3lKszIfW6ufKKBidz9EWtx15G733kXJyaiURFSU03hMy/8eO7akx33EPdb73SG/F19delBK6mU5C6qMMkqxSzhR+6dEbFlgkTufvLcc5F7iOIO5dMhTSqkP2Iu7s6r38g9KrXcNU7kXvt6tC1TTtz1U86UhyB9784ML3lJgrVrS58OqVT+cThMwp6FSRMdW8Z5rZctc9z9f3jZMlrcK9syRtzLE3Nx9++5Dw3ZXHWVIw1hTrFXSP1sGZuWFsGySm9XS4vQIvlJGDRZpXjoqSyvelX50YipVL6mR5iEPQuTJnq2TH06VMOefLyvz5GkcrnuxpbxR6zFPYjn/uCDWR5/3OZDH0rx2tc2dlh1KXp66pcK2dpa+VB2WAsj92GlUKryhOEtLdGwZcIu96vpjEzkHqxDdQ7KdlYeV4puoN3YMrEn1uIexHM/eNA5oW+8sTUyE+vWy3OfmVG0tfkQ94QsEHfdwdrbW36fRCVyj4otszQJR2ybYbt+E65UgxZ3Px2q2rUrV5ZrSCl6Q7JkIP8EUq5elLFl/BFzcXde/UTuBw/aJBKVI9RG0t1dP1umXKaMpsOSBbZMXtzLnwotLdHw3Gfc9jZwHglP3r26lRng3UEGWSwC4+OK9nZIJv3ZMlC+MuS4UpSY86Uh+Jk7wNgy/vAl7iJypYg8LSL7ROQjJZb5DRF5QkT2isg36ttMbxIJoaXFn+f+wgs269aJ50xDYdHWJrnRdrXg15ZpTyy0ZXIea1+lyD0atoxuQirkyP3C7gTvSqX413Q6UMnpejM25s+SgfzkJuVKEIRdTtmyhPb28mm3JnL3R0XzWUQSwE3Aa4AB4AERuU0p9UTBMtuAjwK7lFIjItK3WA0upq3Nf+S+YUO0HlRaW50oxLZV2c7QSszM2P5smbKReyXPPRq2TE7cQ22FwxbLYhoYBZaG1IYg5TS0LVPuHj0F9IR843RKXZT+v/Hc/eFH7S4G9iml9iul5oBvAVcXLfM7wE1KqREApdRgfZtZmrY28em5R0/cta1Ua0RcS4eqFvdKEz1EJRVSR51RGJy4xhXBoyH67mNjpefOLUY/7ZQ73SaVokwViobQ2WlsmXrgR+3WAocK/h5w3yvkTOBMEblXRO4TkSvr1cBKtLdXjtyzWcXAgIqcuGufvFZrZnZW+fPcE0LxfXBQKVYslYqerbFlFrLa7Xg8EqotowLbMuVy3aNQ2qGzs3yROmPL+MOP2nkd6eI9nwS2AZcD1wFfFpEFT6oicoOI7BGRPUNDQ0Hb6klbm1T03I8dmyOTgY0bwxeEQoIOwiqFf1vGI3K3bXqXVd4vUbFltDBF4bqOQuQexJbJee5llgnbcwct7qX/r4OMlsbOAhg7/Ij7ALC+4O91wBGPZX6glEorpQ4AT+OI/TyUUl9SSu1USu3s7e2tts3zaGurXPL30CEnNF63rlkjd5+2TIlUyD4f4h4ZW8Z9jYK4RyNy99+h2urTlgmznDJUrmM0NwfJJDX1U50O+FG7B4BtIrJZRFLAtcBtRct8H3glgIisxLFp9tezoaXwE7mPjjqZvVEZmaqpV+Tu33Nf2KE6qJTPyF0iEblHyZbpFqELJ989LMbGgkfupWwZpRRTQEfkI3dlLBkfVFQEpVQGeB/wE+BJ4NtKqb0i8gkRucpd7CfAsIg8AdwF/E+l1PBiNboQP577+HgWiE5NGU29InfftoxHKuSQUvQu8zGLT8qJmMJM+4NodagCrLEsjoa4T8bH69ehOgvYEHrk7qdD1Yh7ZXyNw1dK3Q7cXvTexwt+V8AH3Z+G0tYmnDoVT3GvX+SuAkXuSilEBKUUJ90O1UroiymTgTCtzihF7gCrRUKL3NNpxewsvucnqOS5R6a0Q8UO1fLlrQ0O0TKhq8CP5z4+7tgyUanjrtFVHOvjufvJlnEiM31x60itq73yd7S0OOsP25qJkucO4Ubu+rwpN3duIbk891ITtrivYYt7pYlhjC3jj9iLe3t7Zc89upF7sJmkShFkEBPkK0PqSK3DhzjoiynsdMgoZcuAUz1xMKTIPT9Ji7/lK5Uf0JZdHPLcTY57ZWIv7n5GqI6PZ2lpid6jXD5yb1CHqlt6QV/EWuQ72/2lQoJjBYSJvrdEo64ndIlT+nex+iLKrVcHBb7F3X2NeuTe2ekMTLRt73Yaz90fTSDu/iL3qEXtUL/IPcgIVSgQ91zkXvk7dKQUduQ+hyNSEhHPvQvH2qpDiaAFPJHNYo2NcWfGu46jtmX8Bi1x8tyBktbM7KyxZfzQBOLux3PPRs5vh/pE7pmMIpvFd8lfyEfs+tWPLaMj99DFXanIZMpAXggXo7b7/+vu7FtLdHTo86bcxOiFVMpzz4m7/yYuCroyZKlcd2PL+CP24u7Pc880beQeJHorGbn76FDVkVIUbJmoZMqAY8sAlCk/XhXTSvEdV9RHStw4qo3cS+W5a1smCnnuQMlcd9Oh6o/Yi3tbm5PBkc2WFp2o2jL1iNzznWqL3aEaDVtmVqnIdKZCPsqdrHPk/ng2mxP1J7JZz2Wq7VAtdQinIhK558XdeO610ATiXnkgUFTFvR6Ru/5sIFumpg7VwE2sK9pzjwpdi2TLDLrr25VI8LRtY3usP2jkXqnkb3Q8d+e1vLhH73qOGrEXdz2PajnfvZk99/wFvtgdqs5r2PVlombLaCGs92TZuhTzZckk08ALHuKus8SCZstUsmXCFnc9DaaxZWoj9uKuo99ygx6iGrnri7I2zz2ALVNTh2o0bJmodaguduR+WSIBwJMe1kz+xu5vnQmcEq+VIvco5LlDpQ7VBjYopsRe3Nevdy6uAwdKDySZmIimuFuWM01gLSNU86MU/UfukzVE7saWmc9iee5Dtk0bsN0Vd6/Kk0FHqIoIKUrnuU+5N85E6LaMH889etdz1Ii9uG/f7kY2T3p3OimlIpstA/6nCSxF/tG88vYtaxG6gcfdEZXVjVAN15aZJVq2zGJlywwqRZ8IS9z1j3qKezBbBhzfvWTkTviWDJhsmXoRe3Ffv17o7IQnn/SO3Gdsm2w2enVlNK2ttU2SHcRzT4rwqmSSH7uTOk/hPKa3+rhQIlNbJmLZMlrc6x65K0WvZdGFc4y8xd15DTLyOiVScsBVFGq5Qz7P3WTL1Ebsxd2yhLPPTvDEE97iPp6NZl0ZTa2Re35WGn/b97pkkheU4hnbZsq9mP2M9oxKbZnIdai6r4vhufeJYInQQ/0i93K2zCTh57hD5Qw4Y8v4I/biDrB9u1XSlom6uNcauWcyzoXqW9zdnMY7MplAEzPkR6hGIFsm1BbMp8X1settywzZNr3usVki4inu+doyASJ3ytRzVwqfg10XlXID5hybVeU6XQ2laRpxHxhQjGcWngxRF/daI3dddqTSBNeaTSIkgcNKMRVgpnsdKUXBlolStgw4PnU9bRmlVC5yB0fcT9XLcxcpKe5povFUVK7UxeiE835/f1NI16LSFHto82ZnMwZmF1oz4676dXU1tEm+aW2tLVsmqLiLCMtFGFaKSaWqiNyraWX9iJotA07xsHraMpPADNDrztFaKnIPmgoJTuReKs89Kv0ZliUkEt7n2uBJ5xrv74/WORBFmkLclyxxDvSYRwmCGTczpN3HKMwwcKpa1h65t7T4P5TLRRi27UC2TGSyZSIiQIV0umV/64XOca9ky8zOOhNFJxL+z+32Mm2dI9xZtgopNSH74Ennvb6+aF7PUaIpxF3P/j7qYcvozqOo9q7XHrk72+c3cgdY4UbusbRliJbnDk7GTD0j95NuQLKiUNw9lnOmVwy27v4yk4tE6akolfI+144PO/vZ2DKVaYo9pCP3UY8+1Vn3RI7aRB0aP/XoyxHUlgFHNE66qZB+85qNLVOaTuqbCpmr+eMjcg96XveLcKxEW9MReipKpaSsLWMi98o0l7h7RO7aXwwa4TQKJ3JvXIcqVBu5O6+hl/yNYIdqvSP34unutLgXz8o0MxM8cl9lWQwpRdajvVF6Kmpp8bZljp9UiMDKlUbcK9H84m4i9wXEvkM13CYsoN6ee67mj47cgQwwXbScE7kHW3e/CDZwopS4R+SpyPHcF74/eFKxYoUEOt9PV5pC3Lu7QcS7Q3XOFfdoe+61zcQEwSP3aWA4gLiLCMlkuOKulIpc+QGof7aMV+QOCwcyVWPLrHIzcI57ibtSEepQ9bZljg/bxpLxSVOIu2UJ3d3xtGVC8dzdC3ycYBUAnU6u8GyZLKCIXuTe5eahlxr5GZQFkXtJcVe+p9jT9LvrOubRqZomOvu21Lk2OKJMGqRPmkLcwbFm4mjLhOW5a4IMN3d8UN+L153cU1jEIvfdySSTwFfqtHNykXtFca8icnfX5Rm5E519W8qWOT6s6OtrGtlaVJpmLy1ZIox62DJ5cW90i/yxdKkwPl59RFyt567ZZvk/BVIpCTUVMpfWGl4TPHldMskrEgn+opac1gL82zLVdagCnhkzURnEBE45jVLZMsaW8UfTiHtPj3fkrgWhJSpmYhEbN1rYNgwMVCvu1Xnumtclk74/VyqDoVHoyD1q92kRYXcyyYBSTNfBmpnCuTC10ObEvWi5mZngkXsX0A4c97Bloj6ISSnF+FQ+gcJQnqYRd8eWWfj+rG2TSomvyodhsGmTcwief770ZCPlqNWWWRUocg/XlpmNqC0D0Ofux6F6iLuboqrP2R73dawOkbuIsErE05aJSm0Z8B7ENAfYdr4ksKE8TSXuXtkys7btq9Z5WNRL3N0Je3yx0r2AXxXkQxhbphy6yFep0Z9BKC4LoftMi2vCVJMKCU7NmsGiddlKkSE6+9YrWybX0dwRjRtQ1PH/TB5xch2qRWfnnFKRFvd16wTLqk3cEwl/Ndk1rSI82tXF1gBRO0THlomKABWSE/c6Ru6aNnfdxUlVs7P+p9grpIuFI2r1PTsq+zZ3rhXsiFxfhBF3XzSRuJcexBTlwv6plLB2rdQk7tUM6DgvYNQOpet9NIpc5B4R66AQbcvUTdwLtlEH5wsj9+C2DDiDrkaKnjB0kNwSkX3rZQHmpoU0towvohvSBmTJEmFWeVwAEbdlwOlUrVbc02nVsNF6HR3CxISJ3L3orWfkDp7i7hW5V5Pi2+ExojYdMcvLywLM1dzpjMYNKOr4Uj0RuVJEnhaRfSLykTLL/bqIKBHZWb8m+kNXhtw+Ps5IwQU2G3FbBhzfvdGRezWsWCEMD4cn7iNuB8OSiESXhXTiZKEM1cNzL7JlxJ3tqTjRspraMuC0daroJqSD5OiI+0ILcMpE7oGoqHoikgBuAnYDO4DrRGSHx3LdwO8D99e7kX7QhYQOKMWvsvnykHO2HdkBTJp16ywOH15YGMoPjRT3lSuFEyfCE/dBN5Trj6C4iwh9InWL3IurdbYBM3XqUPWK3HPiHpF96zVgznSoBsNPSHsxsE8ptV8pNQd8C7jaY7k/Az7NwqfHhnD11S187gwnr+BgQfTkeO7Rjtx7eiCbhZkqgr4wIvdqbkL1YNC92vsCdgQ3ij6PLJRq8KrW2SYy78JSSjE3V50t0ymyMHLX40ECr21x8MqWmTQdqoHwc5WsBQ4V/D3gvpdDRC4E1iul/rWObQtEe7tww2rnoXKeuCsV+chdz+864ZHKWYlMRgWahakWVq60mJ2Fqdqdh6oYTKdJ4VRJjCJ9ZSbCCEJxhyo4vnthf5IWvqoidxyLp7Dsb9SyZbxqyxhbJhh+VMFLGXN7XUQs4LPAhyquSOQGEdkjInuGhob8t9In7QmhV4RDhZ57DDpUtbiPe0w2UolGR+4AJ0IqHnZ8bo4+ie6AtN462jLF4l4cuQ8OulPx9QY/t7XlU2jNRM2WKZfnbjpU/eHnzBgA1hf8vQ44UvB3N3AucLeIPA9cAtzm1amqlPqSUmqnUmpnb29v9a0uwwbLmhe5Rz3PHQrFPdqeuxb34ZDEfTCdzuWTR5FekbqOUC2kOHI/dMg5x9evryJbpuB7NFEbIKY990IL0ETuwfCjeg8A20Rks4ikgGuB2/Q/lVKjSqmVSqlNSqlNwH3AVUqpPYvS4gpsEPHw3KMrCABdXc5r1MVdd1oPp8PxZQbTafoj6reDUyZglnxaYTUod/rDSpH7wIBzDNatqyFy97BlouO5g1JOmWeN6VANRsUzQymVAd4H/AR4Evi2UmqviHxCRK5a7AYGZb0bues7fpxsmWo994bbMh6DxRqBtmWiihbk4hmTgpDGEbRKkbsuNFeNuOt2ThW8F0VbBmCuII6YVM4Ue0Fr2J+u+BqhqpS6Hbi96L2Pl1j28tqbVT0bLIsJnAp6S4mbLRP8s6eLLaOUYnBujr6olvdkvt3RU6VIFk/UoWlj/k1jYMCmowOWLg3+HZ3u62SEbZncnL0Fp9qUUnS0BSu1cToTbdWrgvXugR9wrZk4Re7jVUTEjRT35cvDE/fxbJZZpSJty3hFxEEpnqhD0yYyL899YECxbp1VldB1lovcA69tccjN2WsXeO5ARxW1dE5XonulVMnSovKocfDcu7ud1+o99zo3qATJpLB0aTjZMsd1jnuEo7acuNfguRdP1KFpZf4I1YEBu6rO1MJ1e3ruEdm/OVvGI3I3+KPpxL3LPTnHY+S5d3XFw3MHJ9d9OATPfcgdndobEfHxwisLJSglbRmPDtVq/HYoiNxjYMsUe+4mcvdPtFWvCrS4T7h/z8VgEFMqJaRS0ffcwbFmwrBlJt2SEt1RFvc62DJ6hqSlRdtZ2KGazSqOHFFVi7tup2eee1VrrD95z93YMtXSvOKuFLZSpGPQoQqO7x71VEiAzk6YqqKdtTLlil6QCb0bTbv7Wkvkfl82iwAvLirJXBi5j04oslno7a1uX+gO1XmRu/salWyZlhZvW6azvcQHDAuIvuoFxE0ZZ0Kp/Akb8doy4Pju1dkyjRX39nZhOoQ095y4N/6rfVOPVMhfZrNst6wFlS9byRcOm3AfDXRHfFC8IncdIUclF8nLljGRezCir3oB6S6wZXQHVNRtGXB89+oi98Z67u3tMG2HZ8tEOXKv1XNXSnFfNstLPSZSaXMHSAFMTDnr7+xcsJgvdPA76RW5V7fKupMT94I2TpoO1UA0nbi34WzUhFK5EyM2tozHBN+VCCdyD8GWccW9Sj1rCLVmyzw7Pc2IUrzUI/2pDSejxVaKiWln/bojPigiQgfRtmV0tkzhYOgppehoj0b74kD0VS8gIkIXjrjnI/fob2a1nns63WhxJ1xbJiLi40WtHaovzDiu+pkeufy5qfaASdf3qVbcwcmYmdehGtVsmULPHWPLBCH6qlcFXSKMK5XLLoh6njvEx3Pv6BCmw+hQdTsaq6hw2zBqtWX0TFPLPG5ghZNka1umFnEvjtyjVlvGcxCTUnQaW8Y3TSvuE+QfNeMQucfHc5dQ6rlP2TYdVnUjMhtFiwhJqo/cy4l74STZtXru4EbuRbZMkugM7S8exDSqFOPAymXRaF8ciL7qVUF8bZngn2u85+7U+8jWkO5XDVO2TYdHR2PUKI6Ig3DKb+ReB1umQ2R++QGlImPJwMLaMve7++bicxo0HLsJiL7qVUG3iCPuuQ7V6N/ta8lzb9RMTOBE7lBbul81TGWzdES4roymWDSDMJLJ0EI+m6UQHbnPKMVkjR2qAGtEeK5w3gOi47dDYSqks606///ic6N/g48K0b9aqkDbMjpyj0Oee3s7ZKqIiMOI3AGmTeTuSYfH/KR+GclkWFZipikduc9SeyokwK5kkmdsmyFX4NNEJ1MGFg5iui+b5RzLoqeGG9rpRvRVrwq63Mg9P6Fu9Dcz5zEG/FwYqZBgIvdS1GrLeFky4KRCghO5T0w5ka0WwGrY5d4o73VTTOeUikxnKswfxKSU4v5slpfE4OYeJaJ/tVRBF07hsOPuRdbfH6UHTm/0RMez5RdbgNOhWvfmlMRE7uWp1ZYpJe6thZH7tKrJkgHYmUjQCtzjetlRtWXSSjEDnFSKrTG4uUeJptxbOnKPk7h7jcjzg4nco0UtkftIJrOgYJimMHKfnFa5qRmrpVWECxMJ9ujInWjZMoXZMiPu/ix14zN4E/2rpQq6Xc/9qG3TnUjQ2Rn9iE93+gaP3I3nHiVqitzT6dKRu/uq89w7O2s/5qtFGHaPYzpi2TL6PJvIqpy4l7rxGbxpSnHvEsEGnrdt+lNROmVLE7fIvZayttUwGZPIvZ1F8tzndajWlimjWSbCSV2MDGiPkHimUkJXF5xMK06ZyL0qon+1VIF+Yt1n26yK8JybhVTjudu2QqnTJ3LvjEvkXsW+UUpxqowtU5gKOVEHWwYcsdRR8UmlWBEx8VyxQhjOKGPLVElzirt7Ejxn26yKTeQePFsm4w56anT5ATCeeyk6RKraN+PZLFlKC5iO3KfB9dzrE7lP44x6HbbtaIp72tgy1RL9q6UKtLjPQGxsmVzkHiDqcxMdTpvIPRaeO9XZMuVKD0B+esFB22Ziirp47svddY4oxXAkI3drnribyD0YTSnuWwoiPBO515dctkwDxT2tFBmlYhO5T+LYLEHQpQdKZsuIsEKEAbe2TL1sGYBBpRgjL/ZRYfly4WQm77mbyD0Y0b9aquACy0KfBnER92o894xbrqDZUyFzk0bHIHJf53bmH6yzuOt1H7btuuS5Q17cdRmC6EXueVumC0hGrH1RpynFXUS4xBWC7hgIAlSXLRNO5O68NjJy1zZHHCL389zz7bFssCpw0z6mEVxnWQzYNpPT9cuWASfxAKIp7iMZxzIylkxwon+1VMkn2pxhHxfV4/m1AVST5x6G597aCoL/yH1cKQ7atdUIzol7DG7U51Yp7jPuPmorI2JrLYt9to1t11ZXRhOHyF0BB2zbiHsVNG39zFcnk9g9PUhHlKdUzlNN5J4OQdxFhDbLX6fh09ksZ09M0ApM9/RUXSs8Z8vEIHJfIsJGER4LeEPLiXuZZda5g/MA+vtr3xfaY9/n3oiWR2z/rliRf7Lwmp3KUJ6m3mNRmXjAD1VF7iF47gDtlr90v8/OOd3Ds+SHkH9zbo7+sTGOBBC/ONky4FgzixG5ryvY/ssuqz0u0/7+sxGN3Jcvd9pzzNgyVRGPq+U0IBe5B/hMGJ47QLvlz3M/XCDgL7jL35JOM6gU/2Paf5esHiK/pJEV0mrg3ESCp2ybTICnMC3uXrXcNWsLBG7z5tqPeUKEHvKdv1ETdx25g0mDrAYj7hGhqjz3sMQ94S9yP2Lb9LkX5QuueOnP/Ucmw0mf0fsT7nJnx8RiWyVCBhirQtzLRe5r3Mi9va1+T6VaNJPkR3ZHhUJxN2mQwTHiHhGqynN3PfdGzsQE/iP3o0rlspY+OTvLpRMT82b/Oe5T/PZms6xOpVgek1ISS1whOhXgM348962WxVmWxbc/Vb+bnBb3FSUmCQmTFSvy5/WWmFhyUSIez7mnAbVF7ovQoDJ0+PDcs27J5fMSCf4jk+H+Ag/64kSC/8pmOaYU23183+O2zbk9PTW1uZHoKHM0wLHUqZCtZZZpE+Gp7m64tJbWzafXssC2y35vWCxdCjeuTbFx2OJ3YzJeJUr4uh2KyJUi8rSI7BORj3j8/4Mi8oSIPCoid4jIxo9+QBEAABAXSURBVPo3tbmpznMPp0N1aTJfKrYUQ0ph45SV3VgUdb3UjeaP+7BlbKXYm81ybkwsGSiI3APaMq0hRM9/3dbG9S0tvL81evIuInx2azs3trZGqtZ8XKgo7iKSAG4CdgM7gOtEZEfRYg8DO5VS5wP/Any63g1tdixLSErQbBnntdHivqXdYn8FYT7qCttqyyop7sd8iN8B22YaOLceid0NYkkVkfuMbdMWgvVwfiLBzR0d/GEExd1QG37OpouBfUqp/UqpOeBbwNWFCyil7lJK6XTk+4B19W3m6UFKol84DOCMNosRlS/o5MVRV/xXi7DBFbsdrnhdlEjQgj/PXY+ePDNGkXs1tkxY4m5oXvy4tWuBQwV/DwAvKbP8O4Ef1dKo05VWS2KRCrml3RGh/bbNi0uMGi2M3H87lcrlaP/V7CybLIs+EY75sGW0tbEiJmmQAEvc16C2THsMRuAa4oOfK8ZLOTzPWhF5K7ATuKzE/28AbgDYsGGDzyaePqQC2jJz6XA89zPaHKF+LpstLe6ucK8SYVMyyUuTSWaU4m2pFC0irLIsX7bMqZjluEO8bBlD8+LnbBoA1hf8vQ44UryQiLwa+BhwlVLKU6OUUl9SSu1USu3s7e2tpr1NTasVrPzApJuy0tXV2IivMHIvhR5VWJi33SbCJlfAVon46lAddV/jJO5JEToJHrkbcTfUEz9n0wPANhHZLCIp4FrgtsIFRORC4O9xhH2w/s08PUhZEihyn5h2xKPR4t6VEPpE5uWsFzOmVM6e8KJfhIdtm0smJsr2M4wqRZL4lB7QLBEJHrmbjBBDHal4xSilMsD7gJ8ATwLfVkrtFZFPiMhV7mJ/hTPA7Tsi8oiI3FZidYYytErQyD0ccQdnUEm5yH1SKTrLiJWuzX1/NpvrNPXilFIsieAAm0osrUbcY3YDM0QbX8+6SqnbgduL3vt4we+vrnO7TktSVjDPfcLNTwpD3HtFGCgn7lBW3N+YTPIlt7DYPtvmnBLe/agr7nFjiQinlOKAbXPz3BwfrZCrPW3bTqdxjeWRDQaNCRUiRKsEy5aZmFKIQFtb4w9jtwhjZf4/qRTlMtPf0NLCcHc3QFl7Z7SCvRNVloowCnxxdpY/nZ3l9yoUSjORu6HemLMpQqSsYHnuznRr4ZQ27hFhvExbJ5Wio0K7llsWS6GiLRPHolE6cte5/P+YTucyiLww4m6oN+ZsihCpKiL3eky3Vg3dlK96WMmW0WxNJCpH7jEV91Gl5m1bORvL5Lkb6o0R9wjRGjByr9dcmtXQI8IMkC7R3kq2jOYMd+q4UsRV3JcWiPt5bkR+pMyxNZG7od6YsylCpKxghcOcyH3RmlOWHldwS1kzlbJlNFstixdsu+RNIq62zHL3KeyoUrzCzdGvaMvEcDsN0cWIe4RoleB57p2dIdkyrhCV6lT1a8usFyGLd50ZWynGIZaR+6sLBl29NJFAyJdk8MJE7oZ6Y86mCJEKOEJ1YipcWwa8I/e0UqTBly2zpMx6xnHqXMRR3C8sEOqzLItekZKRe1Yp0koZcTfUFXM2RYjWgHnuk9Ph2TLd7qtXp+qk++onctc3Ca/15OrKxFDcRYTfcGeO2ppIsEakZOSuj7kRd0M9MWdThIhTtky5yH3Sfc9P5F5O3PUIzzh67gBfb29nb1cXy0RYbVkcKRG5z7jbacTdUE/M2RQhgmbLTEwTvudeTtyDRO4e/xuNceQO0CLCDje9cXWZyF0PbzKpkIZ6YsQ9QsQxcvcS5XrZMgfdSLcvpuJeyGrL4qhSnvVmTORuWAzM2RQhWi1I42SJVCKtFHNpwvPc62XLuK9e4v7TTIblIpzTBKK33t1fS8fGGCqyZ2bcVyPuhnpizqYI0ZNwJ3nwsayOjsMcoQq12zKlbhJKKf49k+HVySSJJojcr0ul+EN3FvSnSol7E2ynIToYcY8QK1qci/ukj8qAE0qX+w1HEBIidFAicndf/Yh7iwhtLLxJPGHbHFGK18Zoko5y9IjwO664HygWd2PLGBYBczZFCC3uwx5R7OEiQdDi3unH+1gkekTKR+5B1lP03uNZZ4LYlzRRJ+NGV7yfLzqWU0bcDYuAOZsixIqkt7j/z5kZ1o2PzxP4sCN38BZlCGbL5NZTtM0n3L+boTNV0yrCGpEF4j6kt9WN7A2GemDEPUJ4Re4jSvEZd1KLPW40C+F77uD45bXaMuB0qhaLuxa85U0k7gCbLGuBLaMnCl9lxN1QR4y4Rwgvcf+XdDr3+68KxF2P3uzuJjR6wDO1rypbxiNyXy6Sm46vWdhsWQsi92NuXZmeJrKgDOFjxD1CLE0KwnxxP2jbWDhzlj6SzbLPFXhdJ3zTpvAO4VqPKBQccU+Bb2EuJe4rm0zYwYncDyk1b7DaMaVYlUrFbp5YQ7Qx4h4hEiIsE5kn7kdsm1UiXGRZfC+TYdvEBM9mszxr2yzrEVasCE8Qzk8kOKLUguwevxUhNV7iPmTbTSnuWyyLLNA7NsYxd78dU4pVbh0ag6FeGHGPGCuKxP2oUqyxLLYUZFI8kM3ybDbLtg1WqNGenoTisSJxHwk476lXx+wJpehtQnF/S0sLH0qlGCdvsx2zbeO3G+qOEfeIsdwjcl8jwvWpFK9zc74fcSP3bRvCPXznuR7xYwV9AQAv2DabAqT16chdFWz3UJPaMt0ifKi1FcjPHattGYOhnhhxjxgrRBguiISPuJH7uYkEP+7s5ELL4v5slkNKhS7ua0ToBt4/M8NX5vJVcQ4EFXecsgu69K1SyvHcmzTve5U7AOw5dwaqE0bcDYtAc149MWaFCCfdCHZOKYaUYk1BBPuiRIKfZ7Mo4MyQxV1EONuN3t817dQ2nFGKI0qxOWDkDvnMm3EcsW9GWwac/abnjh1SCoVJgzTUHyPuEaPQc9f5z2sKhPKCgnS587aFnzp3c3s7v+baRSdtmxfcp44g4r6+aOSmLqzVjLaMZqtl8Zxtmxx3w6JhxD1i9IowgTMCVU/uUBi5/18tLfxGSwvf6+jg3K3hi/vZiQTvdYVpr23nBDqILbPdXfYJ97N6dGozi/sZrrg/4fZXrHN9eIOhXjRHVaYm4kI3Mv+vbDY3UKkwcl9vWdzS0RFK20qhJ6R4wrWLIFjkvtmySAFPukI3dBqI+1bLYhb4m9lZ1ohwYVi1mw1NixH3iPHSZBIB7s1k0HH5uoiL3AYRunAi93aghflPG5VIinCmZfGkG7k/ks0iwLYm7VAFuDyZJAU8bNt8IJXCivgxNsSP5r16YspSd3KKe7NZvpNOc0kiwYqIi5y408k9ms3y75kMOywrsFjtSCRytszPs1nOtSyWR3y7a+GsRILPtbVhAdcbv92wCDTv1RNjdiWT3JnJ8Ihtc21MRi6+JpnkZ9ksj9g2N1bhH293SxmMKcUvMhkubZI67uV4d2srQ93dOSvOYKgnRtwjyA2pFCvdSSzeEhNx/3hrK5cnEmy3LH6rija/LJFAAb83Pc0kcOlpInjN/HRiCJfmD49iyEWJBM93d3PCHcAUB1Ii3NHZySzO7EpBeU0yyWuSSb6ZTtMvwqtOg8jdYFhMfCmHiFwpIk+LyD4R+YjH/1tF5Bb3//eLyKZ6N/R0IyUSG2HXWCK0V9kxKCJ8qb2dG1MpHuzqinw/g8EQdSpeQSKSAG4CdgM7gOtEZEfRYu8ERpRSW4HPAp+qd0MNzc8my+Kz7e2sNcJuMNSMn6voYmCfUmq/UmoO+BZwddEyVwP/5P7+L8CrxBSnNhgMhtDwI+5rgUMFfw+473kuo5TKAKPAino00GAwGAzB8SPuXhF48dxqfpZBRG4QkT0ismdoaMhP+wwGg8FQBX7EfQBYX/D3OuBIqWVEJAksAU4Wr0gp9SWl1E6l1M7e3t7qWmwwGAyGivgR9weAbSKyWURSwLXAbUXL3Aa83f3914E7VeHMCwaDwWBoKBWTiZVSGRF5H/ATIAF8VSm1V0Q+AexRSt0GfAX4ZxHZhxOxX7uYjTYYDAZDeXyNFFFK3Q7cXvTexwt+nwHeUt+mGQwGg6FaTEKxwWAwNCESljUuIkPAC1V+fCVwoo7NCROzLdHEbEs0MdsCG5VSFTNSQhP3WhCRPUqpnWG3ox6YbYkmZluiidkW/xhbxmAwGJoQI+4Gg8HQhMRV3L8UdgPqiNmWaGK2JZqYbfFJLD13g8FgMJQnrpG7wWAwGMoQO3GvNHFI1BGR50XkMRF5RET2uO8tF5H/EJFn3ddlYbfTCxH5qogMisjjBe95tl0c/o97nB4VkYvCa/lCSmzLn4rIYffYPCIiry/430fdbXlaRF4XTqsXIiLrReQuEXlSRPaKyB+478fuuJTZljgelzYR+S8R+ZW7Lf+P+/5md0KjZ90JjlLu+/Wf8EgpFZsfnPIHzwFbgBTwK2BH2O0KuA3PAyuL3vs08BH3948Anwq7nSXafilwEfB4pbYDrwd+hFMx9BLg/rDb72Nb/hT4Q49ld7jnWiuw2T0HE2Fvg9u21cBF7u/dwDNue2N3XMpsSxyPiwBd7u8twP3u/v42cK37/heBd7u/vwf4ovv7tcAttbYhbpG7n4lD4kjhZCf/BFwTYltKopT6OQurfZZq+9XAzcrhPmCpiKxuTEsrU2JbSnE18C2l1KxS6gCwD+dcDB2l1FGl1EPu7+PAkzjzK8TuuJTZllJE+bgopdSE+2eL+6OAK3AmNIKFx6WuEx7FTdz9TBwSdRTw7yLyoIjc4L7Xr5Q6Cs4JDvSF1rrglGp7XI/V+1y74qsF9lgstsV9lL8QJ0qM9XEp2haI4XERkYSIPAIMAv+B82RxSjkTGsH89tZ9wqO4ibuvSUEizi6l1EU4c9K+V0QuDbtBi0Qcj9UXgDOAC4CjwGfc9yO/LSLSBXwXuFEpNVZuUY/3or4tsTwuSqmsUuoCnDkwLga2ey3mvtZ9W+Im7n4mDok0Sqkj7usg8D2cg35cPxq7r4PhtTAwpdoeu2OllDruXpA28A/kH/EjvS0i0oIjhv+fUupW9+1YHhevbYnrcdEopU4Bd+N47kvFmdAI5rfX14RHQYibuPuZOCSyiEiniHTr34HXAo8zf7KTtwM/CKeFVVGq7bcBb3OzMy4BRrVNEFWKvOc34RwbcLblWjejYTOwDfivRrfPC9eX/QrwpFLqbwr+FbvjUmpbYnpcekVkqft7O/BqnD6Eu3AmNIKFx6W+Ex6F3atcRS/063F60Z8DPhZ2ewK2fQtO7/6vgL26/Tje2h3As+7r8rDbWqL938R5LE7jRBrvLNV2nMfMm9zj9BiwM+z2+9iWf3bb+qh7sa0uWP5j7rY8DewOu/0F7Xo5zuP7o8Aj7s/r43hcymxLHI/L+cDDbpsfBz7uvr8F5wa0D/gO0Oq+3+b+vc/9/5Za22BGqBoMBkMTEjdbxmAwGAw+MOJuMBgMTYgRd4PBYGhCjLgbDAZDE2LE3WAwGJoQI+4Gg8HQhBhxNxgMhibEiLvBYDA0If8/sKyn/Z5zk+QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(L):\n",
    "    if hidden[i] == 0:\n",
    "        plt.axvline(x = i, color ='yellow' )\n",
    "    else:\n",
    "        plt.axvline(x = i, color = 'red')\n",
    "plt.plot(range(L), post_prob,color = 'black')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2f05e6d14a8>]"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztvXmcHFd57/17qrfZNVpmRstotWRLsq1Ytq5sENjGBmw5YJlLyLVvMNwAceJgEoO5AUIub17CJQZC4L2JgRAMxHABs9ggEgMJXgAb7Fi2heRFsseSJY1GmhmNRrNPb/W8f5w63T091d1V1dVdp1rn+/no05qe6p5T26+e8zvPeQ4xMzQajUbTWBhBN0Cj0Wg0/qPFXaPRaBoQLe4ajUbTgGhx12g0mgZEi7tGo9E0IFrcNRqNpgHR4q7RaDQNiBZ3jUajaUC0uGs0Gk0DEg3qDy9ZsoTXrFnj4ZMHAcy42L7ZenXzmYBwu2sAcC6Aljrto5f2+UWztY8zM9Y+B9SOSng9RoX75waVj8U8QnQvesbpPjYDOM/TX3jqqadOMXNXpe0CE/c1a9Zgz549Hj55JYC9Lra/yHp185mAuBLum/ktABfVaR+vrP2fKIncx717rX0OqB2VuBLejlHh/rlB5WMxjxDdi55xuo8XAXjE018goiNOttO2jEaj0TQgFcWdiL5KRENE9GyJ3xMR/R8i6iOifUR0sf/N1Gg0Go0bnETuXwdwbZnf7wSwwfp3C4AvVt8sjUaj0VRDRXFn5l8COF1mk10A7mHB4wA6iWiZXw3UaDQajXv88NxXADhW8HO/9Z5Go9FoAsIPcSeb92xXACGiW4hoDxHtGR4e9uFPazQajcYOP8S9H8DKgp97AQzYbcjMX2bmbcy8raurYpqmRqPRaDziR577bgC3EdF3AFwKYIyZT/jwvcozPc149NEMDh0ycfo0IxYj9PYSLroogk2bIkE3T6PRnMVUFHci+jbE1IwlRNQP4P8BEAMAZv4SgAcAXAegD8A0gD+sVWNVgZnxmc+k8Ld/O4szZ+y3ueyyCL7whWZs3apFXqNRlfFxxje/mcJTT2XR02PgHe+IYePGxrhnK4o7M99U4fcM4L2+tSgE3HLLDL7ylTTe9KYo3vveOLZsiWDRIkIqBRw9auKhhzL4279N4oorJvHQQ63Yti2wicAajaYEhw6Z2LlzCi++aKK7m3D6NOOzn03i619vxk03xYNuXtXoGaouuf/+NL7ylTT+4i/i2L27BddeG8Py5QaamggdHYQLLojgz/4sgT172rBoEeFtb5tGKmU7vqypIb/4RQZf+1oKQ0Nm0E3RKEg2y/j935/G8DDjwQdbMTjYgf7+dlx2WQQ33zyD557LBt3EqtHi7gLTZLz//TPYssXAJz7RBCK7RCHBihUGvvjFZrzyCuNrX0vVsZWaf/iHJK68cgrvetcMLr10EocPa4HXzOWrX03jqaeyuOuuJlx1lehZ9/QY+MEPWtDRQfjzPw9/cTMt7i741a+yOHKE8ZGPJBCLlRZ2ybXXRnHZZRH8/d9rca8XfX1Z/Pmfz2LXrigeeqgVIyOMv/iL8N+ojQIz47HHMjh+3P8H7uc+l0Rv7zje9KYpzMyU7y1/8YtJXHJJBDfeGJvz/pIlBj7ykQQefDCLF18Md/Suxd0F3/pWGi0twJvfHKu8MQAiwn//7zG8+KKJl18O94USFr70pRQiEeALX2jG614XxZ/+aQL33ZfR0bsifPrTx/Ca10xh1aoJ/PjHad++d8+eDD7wgVksXWrggQcyuPXW0g/0Q4dMPPOMiZtuitn2vqXg/+AH5dt34oSJe+5JIZ1W03bV4u4QZsb996dx/fUxtLZWjtolO3eKLt9PfpKpVdM0Fskk46tfTeMtbxHjIADwvvfFQQTcfbfuPdWTQ4dMfOADM3Meqvv3T+Iv//Iwbrghik2bDNx++yxmZ/0Rxg9/eBbd3YSHHmrFe98bx7e/ncbp0/YPdCnab32rfZC2cqWB7dsjuO++0vfsk09msGHDBN75zhm8610zEHklgmxW9E6y2WBFX4u7Q44eZQwPMy6/3F2a1Pr1Eaxfb2hxrwN792YxOsr4b/8tf9OuWGHgkksi+OUv9fGvlj17Mvj0p5MVI9XTp01cc80UPve5FH7ndybwyitCZO+9V8xK/+d/bsZnP9uEQ4dM3Hdf9dH7mTOMhx/O4o//OI6ODsK73x1HKgV8+9v23/3QQxlccIGBNWtKy9+b3xzFnj1ZnDljv6933plEIkG4/fY4vvnNNB58UPTMx8YyeM1rnsFrXjOFD35wtup9qwYt7g55+mlx8i6+2H0O7BVXRPDkk9qWqTWPPy6O8WWXzT1HO3aI46+zlrzz4otZvPGN0/jQh2Zxww3TcyLVYu6+O42+PhNf/3ozJieBr39d9Jp++MNTuPzyBViyxMDrXx/FwoWEhx6q/qH78MMZmCbwxjeKXvJFF0Vw4YUGvv99e3Hfty9b8T6Wv9+/f/59299v4kc/yuDd747hf//vJrS05HsDX//6STz++AR27ozi859P4eGHgwsqtLg75Omns4hEgC1b3Iv7pk0RDA8zRka071tLnngiixUrCCtWzL2sd+yIYnY2/4CuRCrF2LVrCtdfP6UfCBaf/GQS2Szjfe+L44EHMti/X1zLY2OMr3wlhe98J5UT/B//OI2LLjLwznfGcdVVEdxzTwovvZTFc89N44YblgAAIhHCFVdEfBH3n/88g7Y24NJL8/fmjh1RPPNMdt5D6NQpEwMDXPE+lr/ft2/+NbN7dxrZLPCe98TR0kK47roo7r8/jWyW8Y1vDGLr1jbcd18LmptF6nRQaHF3yNNPZ7Fpk4HmZud+u2TTJnGYX3hBi3steeKJ7JwbXLJjh3jvscecifsdd8xi9+4MfvzjDO64I9iutQqYJuOnP83gd383ho9+NAGivGh9/OOz+KM/msFNN83gN7/JYmTExGOPZXNJBzffHMfhw4wvfUlE71dd1Zn73quuiuLwYc7ZNl555JEMXvva6JwMtq1bIxgbAw4fnivu8qG0ZUt56VuxgrBwIeW2L+Spp7JYsoSwYYP4jre8JYbBQca996bx1FOTuPnmHjQ1EV73uugcO/bMGS7b4/EbLe4OeeaZrOdSArLOzIEDWtxrxego49AhE9u3z58NvHSpgZ4ewgsvVBb3bFZMR3/722N4+9tj+OY3U4EPjAXNvn0mBgcZ11wTRU+PgR07IrjvvjSSSca//EsaO3dG0dEB3HVXCj//ubBI3vQmcR5e/Wpx7d97bxqGAZx7bn4178svF9v85jfeo/dMhvHSSyYuumjuvbl1q5C2vXvnnnMZif/O75S/l4kIF15o2EbuTz2VxSWXRHKZNtu2ie/6ylfEA2znzkXWaxR9fSb6+rJ4/vksenrGsX37FI4dq48OaHF3wOws48QJxrnnejtcq1YRmprgSFw03pBZGeedZ3+O1q41HEWIe/eaOHNG3Jg7d0Zx5ox472zmZz8T4is97Z07o9i3z8Q996QxMsJ4//vj+MM/jON730vj17/OwjDy4rlunYHmZuD4ccY55zQjkcifH3k/vfyy9+N77BgjnQbWr5973i+8MIJIRARlhezbl0V3N6G7u/K9vGVLBPv3Z2Ga+Yf77CzjuefMOZ79unUG4nHgF7/IIpEgrF/fDAC5yVG/+U0Wd96ZRDQK/Pa3WfzjPyY9768btLg7YGBAnNxiL9cpkQjhvPMMHbnXkKNHxbFdtcr+HK1ZY8zrotvx4INCyK66Kpq7OeV7Zyt792axZg1h2TJxbDdvFsL21a+m0NoKXH11FNdeG0U6DXz/+2msWWMgkRBRbSRCOP98sf2mTS1zvre5mbB8OeHQIe/3RV+fEO9zzpl73puaCJs2GfMi91deMec9CEqxYYOByUlgZCR/3Tz7bBaZzNzEimiUcO65BkwT2LSpFdGo2Pf16w0YhrADv/WtNG65JY4rrohi9+76XE9a3B3Q3y8uvt5e93675LzzInjpJS3uteLIESnu9udo7VoDR4+aFS2Whx/OYPNmA0uXin+bNxuBZjyowMsvm9iwIS9mGzcK2Xj88SzOPz8Cw6DcAOTAwPwe7oUXip+LxR0QUW81kXtfn/isnWCfc46Ruy4kJ04wli93dh/LuRIyuAOERQUAF1009+9J6/WCC/L7GI8TVq8mfO97YgD2rW+N4frrozhwwERfX+1nTWtxd4CcKt3b6/1wLVtGOHlSi3utOHrURHMzsHix/Y27Zo2BTGbujSq5664kfvhDMUB44MDcNLmLL47gwIFw2mlPPpnBiRPVX3N9feacyHjdOgMR6xBdcIF4f9kywqJF4tjPF3f7yF1+V3WRuzjvy5bNP+/Llxs4fnzu+R4YMHM9kErI7yw8hjLQK+4hyqSJCy5onfP++vURnD4t2nDhhZHcQPOPfzziqA3VoMW9BJkM47WvncQb3jCV65ZXI+49PYSJCVSseaHxxtGjJlatMkoWc1uzRrxf7Ls//XQWt902ixtvnMb+/VkcO8ZzJresWWPg2DFGJhOu8/bkkxm8+tVT2LWrfE56JUZHGaOjPEfc43HCunXiZ2m5yAFIYP64xxVXRJFIAJdd1jHv+885Rwiw15mq8sFjGPPP+4oVooyv/O7pacbYGFxH7idO5Nt2/Dijq4sQj8/9js2b5fEoFnfx/po1hAULCGvWGPjJT1rwnvcsdbiH3tHiXoJHH83i0Uez+PnPM7j77jTa24H2du+2TE+PONRDQ+ESibBw9CiX9NsB5AS7WNz/6q9msXgxobWV8K53zSCbxTxxz2aB/v5wnbd3v3sGkQjw5JNZ/Ou/ereVZE2kYttDCriM3IF8hH7uuXMzUS6+OIKpqQ6cd5595M48/7w4pa+vtIdebKvICNxt5D4wkG/b8eMmVqyYrwNvelMMf/M3CbzhDQvnvC/TJQuzc669Nob29tqv8aDFvYBCP/a++9JoasqndHkdTJV0d4sLYnBQWzO14MgRs6TfDgCrV4vzV1jrJJtlPPJIBm9/ewxveIOYbg4If15S6qGgMqdPm9i/38T/+l8JrFhBuOce73V1pKddPGApxV1G7gDwqldFEIsB558//16JROzPjfxer777iRNc8t6UIizFWYq808i9qUnkuhdG7gMDpu3fa20l/NVfNc3JBgLyD0Uvkx+rRYu7xYkTJrq6JvCFLyRzRcKuvTaa88imp6uL3Hp6pLiHIwIcGjIxNRWOts7OMgYHOSfgdiQShCVLCCdP5vfp0CETMzPixiucw1AYuUuhD1NVyWeeEW3dvj2CbdsieO45722XoittGMl73hPHJz+ZmCOUN94Yw8svtzuOjIF80FQooE7JZIRltGSJvVjLyF2OmUmRl+87YdkyKorcnQ/IAqIUQiIBXHll/Vdj0+Ju8alPJTE6yvjoR2dx4ICJ/n7GlVdG8brXiZv+6NFqxV0c6jCI+8QEY8uWSfzBH0wH3RRHSMGudNMuWkS5wS0AePZZcdNeeKGRy34wDGDlyvzN29tLMIxwRe4yt3vr1gjOP19kaXkto3D0KKO7m+ZVQt24MYKPfGTugjWGQVi50p2kyAHwwnRDp8hzWUrc5YMjb8uIV7vB11IsW5aP3FMpxtBQ6Z6CHStXGhgd7cil1dYTLe4QYvZP/5TClVdGcOYM8LGPiSnn69YZWL/ewK5dUXzve/P9QjdIWyYMy759/vNJDA4yfvSjDPbsUT8NcHRU3HwyW6MU88U9CyKRty1nOPb20pxp7PE4YcUKCp249/YSliwRqZyZDPDii97af+qUWVI8/aClBYjHUbI8bznkA6FU+zo7gaamuZF7IgEsXOh8f5YvN3JevRR5txatl5IlfqDFHWLm6OwscPvtopv5058KQVu7VmRf/PCHrfi933O2QEcpxBqr4Yjcv/jFFK66KoLOTuAf/1H9OuhS3CvdtMXivn+/iXXrDLS2Enp6DCxbRrZlYNeuNUJlyzz9dD6dU3riXtcEHRnhkumlfkBEWLyYPEXup05JcbeXMSLC8uVGLnIfGGAsW0Zll8csRkbuzJx7SNgNqKqIFncABw/mp65v3hzB5KR4v3BgzQ96egzlxX10VJRa2Lkzhv/yX6JV+bX1Qop7Z2dlcZfbAkLwCgf/PvOZJnzoQ4l5n1u61AhNlhMz4/BhMzfguXGjmCXp9TzWWtwBVC3u5dq3fHneMx8ZYXR1ubune3oMpFLA2FjhgGw4ZDMcrawxBw6YiEbFyL2cjGDnM1ZLdzcpb8scPCgivI0bDZx3noGDB+eXTVUNuaBCpch94UKa0/0/ftycE6n/wR/Ecd1183to4nNqHwPJ2BiQTIoHEiB6jGvWGJ5tGZXFvZItA4hzJ6+P8XHGggXu9kVuPz7OuUy3pUt15B4aDh7MYt06A7EY5cTd76gdEBkzqkfusv7NeecJcZ+YwJwMExVxY8ucOSNSIJNJMaFFjoVU+tzoaH3LtXpFzoIuFKClS70FFcystLifOmXmPl+KBQsI4+Piu8fGGB3z51GVpaMjL+5jY5z7zjCgxR1C0GQ3VhZFKk798oNiW0BFDhwwEYuJh5s8JtK2UpXRUYZhAG1t5beTA65nzoglEwE46qYvXEjIZICpqaqbWnNk8CAjdwDo6qLc/rphagpIpYDFi2srE9XYMs3NQEtLeXEfGxP/Hx/nnFg7RT4MxscZ4+Ni8LepSYt7KMhmRT1oWQxJRu61EPe2NsLkpNrifvBgFuvXG1alu0juPZU5c4bR2Um2U9ALkeJ++jTnPHQnkbvsEYTBmpGRu5xXAXgXdym4tY/cDYyMuO8ZjYyUznGXdHQQxsbEd3sT97mRe1iidkCLOwYGGKmUKPADAN3dBu65pxl/8idx3/9WezthchI1794//fTcGtRuOHAg/6BbuZLQ3ByOyN1JettccRf75NSWkX9HdfKRe6G4Gzh1il1fE/US90WLRM9oYsLd506dqizuCxYQslnRCxkfd2+pFIu724dDkJz14i5H3Lu68ift5pvjVRUJK0VbG4EZmK7h3KDHH8/gkksmceed3hYEOHbMzM30NAwxGKd6jvfoKFfMlAHyEfjoqLfIPQzifvIkIxqdO/7Q3S0ETg4sOkWu+VsPz138PfcPn0qW0YIF4nVgwAQzXItzfkBVDsi6+nignPXiXq/oBADa28XrxETtREKubfk3f5N0LcrJJGNqam72wZIl3vzQenLmjPvIXdoUTlbkCZst09Mz16KSgYtba6Z+towUd3fXq5PIXYr5sWM852enyO3HxrQtEzrqFZ0A+aqStRT33bsz2LDBwOws8O//7m52qd3N7HWwq564tWXe//5ZfOtbKSQS+Qeuk8+FIXIfHOQ5fjsQHnF3+/B0ct6lGMs67G6zZeQgvfbcQ0h+llvtT1pbm/gbcpKU3xw+bOLAARO33hpHPJ5fgswp+bzh/GURFnHv7Ky8nRSCoSHG00+b6OpyNlsxXLaMOSdTBshnBHkV90plHarFqy0zNcVobS2/jRRjuSi1W3E2DEJ7u/bcQ0m9LmCg9pG7zGrZvj3iafkyuxl/UtxVzfFmZse2jFzbUuJ0klpbGxCJhMOW8Tty7+jAnFo7tSB/Xzj/jGkyZmYqn8NqbRn5GZEK2YCROxFdS0QHiaiPiD5s8/tVRPQwET1DRPuI6Dr/m1ob6nUBA/kuXq3EXS4ptmKFKHgma3E7xd6WMZBO1663UQ2nTplYsWICqZTzYlCPPtqKe+4Rq9M7zQIiEnW9wxC52006kuLudiJTvcRMCrSbstoyKaGSuBdH7l7FfWyMPWXbBElFcSeiCIC7AOwEsBnATUS0uWizvwLwXWbeCuBGAF/wu6GSmZks/u3f0r59nxiUqU8HRkYotcp1l4WNli+nnLi7ibjtxh+qKclaax55JJur1FduIkshO3ZE8ba3iRIDTjJlJIsWEX71qwx+9St1q2RmMozZ2fkrhiUSomid28h9asp576YaWqyCq+7EXWzr1pbxKu4DAwzT9Pb5oHCiatsB9DHzIWZOAfgOgF1F2zAAOVSxAMCAf02cyyc+cRTXXz+Nxx7z5yarx/RqiZfupxv6+010d4v1HdevNzA97a50QKkB1cLfqYS8Bm67LY4bbnBetbOpifDzn7fi0UcrKEMB8bgovvXWt057nkNQa2Tvym45yK4uw7W4T09zXcQ9Hhe2l5sZwHLbSg91OWBerbjLAdmGitwBrABwrODnfuu9Qv4awNuJqB/AAwDe50vrbPjQh1Zi9WrCO94x48uixfUUdzmgWktbRpYjlct7ubFmRkbEAFXh9GqvaWr14Ne/zuLyyyP4h39oLrt+qh1XXx3Fhg3Olz6TFs7wMOP559U7FkD+urIT944O99fd1BTnoupaQkRoaXEXuctVwio9fCIRQlsbciUI3GbLyM9Iz77R8tztjl7xWbgJwNeZuRfAdQC+QUTzvpuIbiGiPUS0Z3h42H1rAXR0RPGXf9mEQ4fMqibXjIyY2L59Env2ZOso7uK1lraMXEhAFj5zc4xOnZr/oFM1cp+eZjz9dBavfnV9Vrj57ndbcMcdYtbyww+rac3I68quxk5rK7leNnF6uj62DOC+fXlxr7ytjLZbW0uv5VqOwmi/0SL3fgArC37uxXzb5d0AvgsAzPwbAE0AlhR/ETN/mZm3MfO2rq4uby0GcO65otmHDnkX93vvTePJJ0V2ST3SIAExaJtI1C5y7+/n3MxaL7nZdr0YVcX9xRdNZDLAJZfUZ+HhG26I4e/+rhmrVxMeeURNcS8XuYu6Ru6+r16ROyDsFTczt6Ut4+ThI8XZq19eKOiN5rk/CWADEa0lojjEgOnuom2OArgaAIhoE4S4ewvNHSCLelWzOs63vpUflJV2ST2Q9WX8ZnbWxMhI3paR0/Hdi/vcS0I+JFQTd1lb2816mH6wY0cUTz2lZiG1cp57ays8RO718dwBVGHLVN62WnFv2MidmTMAbgPwMwAvQGTFPEdEHyei663N7gDwR0T0WwDfBvA/uIaJ0cuXE+Jx75H7wICJxx7LYssWsfu1nDFaTHt7bf7ewICoJSNtmWhUTL5wK+7FvZhYTGRaqCbubmrD+MmKFYSTJ9XM+5fXlV2w4sWWmZpynoVULaJ9zrd36rkDyC3Isnmzt6y4wrpFYRJ3R4YlMz8AMVBa+N7HCv7/PIAd/jatNIZBWLvW8CzucnLPnXc2Yd8+E+98Z3Xro7qhrY1qIu4nToi1TgsjWbe52adPs+1kLhVnqcrqhz099Z2H19NjIJkUA3ROZsXWE+m525VU8GLLiMjdh4Y5wH3kLl6diPvXvtaMT32qCb293oT57W+PYXCQMTPDWL68wcRdRdat8y7uMh981SoDO3fWT9gBacv4L5Tj48IHLowy3Ir75CTbdukXLKjNA6kahoZMx7Vh/ESW0h0cNNHZWR+/3ynlI3d3tgwz1z1yl/elE9zYMs3NhNWrve9HV5eBO+9s8vz5oAht+YHqxD0/k7PeiMjd/++dmBA+cKE4F64fWQmx9BxsB9C8dOlrzdCQmGbvZiV7P5A9BRWXS5TXVakB1dlZcZ6dkEwCzPbXQy1wP6DqbvLa2UhoxX3tWgNnziC3rqEbBgZMNDcHk7NaK899fFyIe+Hgj5vIXd5YdjeL2xuvHgwOsqNyvX4jI3e54pFKTE4yiEo/oAHnE4XceNp+0Nrq3pYhAprCF1DXjdCKuxz481LMSUz2Meoe9QGiVsuxYyaOHvVXHEpF7s7FvXQ3160fWg+Ghsx5BbLqgfybakbuwiO3W25QnlenPbByD/ta0NLibkBVjgcEcQ+HhdCKuxy1dru6DCAid5kyWG/e//44DAO4+WZ/Q2HpuRd60I0fuQfxcCZEIu7KOtSLUmMmQGG5aWftduNp+4GXyL1evYqwElpxlwOHXsT9+HHG8uXB7PrGjRHcemsCv/511rH/6YTx8ey8GXgLFwpRTqUq/51yHqZqkTszW557/c+hYRC6uymXZ68SExOl52y4tWWCiNzdjAmIWu5a3MsRenF367kzszVNP7gLY+1asSDwwIB/gjkxkZ0XtblZZEKKdym/ViVxPzMBpNP1z3GX9PSQsrZMqewht7ZM/T138XdmZpxt72ShjrOd0Iq7V1tmdFRkhQQVuQP5SRVHjvgX/Y2PZ+bNwJMPwPvuS2N2tvxxKlcfWzVbZui0OG5BifvSpUbD2zLlHva1wG3ZX23LVCa04u7VlpH1v+s9bb2Q1atrIe6lI/c//dNZ/PM/p8p+vtzN3NICpFJARpFZmWOWQBXm9NeT7m5yvfBFPZiYYN9sGTeThPxA2j9usnl05F6e0Iq7LN3pVtzHx4MVBiAv7tVUtSxmYiIzr5xp4epEjz5avh5KufrY8r1pRUqqTFld96Ait1rVB6qWyUn7HHfAS7ZMfSN3t6sxiaJmOnIvR2jFPRIRtVPceu7lKufVi+ZmMSjnd+RebMsUVnj89a/LVzLMp0LaD6gCwLQii1RMTpcubVsP2trUm9QFlPfcvWfL1K9wWOHfrYS2ZSoTWnEHRPTtNnIvN4uvnqxebeDIkdoOqK5fb+Cee5rxiU8k0N/PudVk7Kg0oAqoFLnXV3iKaW0VMzjTaXUEvq8vi9OnS3vuqmfL5CN3Z9trW6YyZ6G4Bx+5A2JQ1U9bxm5AlYhw881xXHONKCH0m9+UVmdHtowikbu0ZepZqrkQ+XfdTLqpJcyMN75xCm1thHe9K267jfvIuL557m7bNzsresCa0oRa3BcsoCpsmVq0yDnLl5OvU9jtBlQl550nClyVs4EqDagCwLSPefnVMDlTX+EpJi/uahyPkRHG4cOMj340gfPPty9mFokQmpvdZMsA0ago+VwP3EbuqRQjkahhgxqAUIt7NZF7UFGfRFRahC+LLSeTJtJpLrk+ZFubWIT41KnSf2t6Gkgk7Jchy0fuVTfVF4K2ZWq9XKJbpL0nU2xL4aZmej1XYQLy59Lp/ZxKAfG4jtzLcRaKu4hIgn7qd3QQmP3p2svSA6VWmiEiLFlCOHWqtDqXyz7IdZlVidytqDJu70DUHClEqmTMyDpFMgurFK2tzh9Ik5NkCAeXAAAgAElEQVT1nQG6ciWhq4vw0EOVlzBkZqTTwZ3/sBBqcRe2jLvPiIyC+peKLUYKsZeqlsXYFQ0rRoh7uci9dKSW6zIrFLkHWTTKbeZJrZF226pV5Y+H09LNzIxf/CKLLVvqV68+EiHs2hXFv/1bGslk+TamrRUytbiXJ9TiLiN3N0ueiVl8NWyUQ6SFIvPuq8Gu3G8xlcW9dGZEPs9dDTET4h7cw1naMqp47kePmmhpwbzFzYtpahJZPpXYt8/Eyy+beOtb67uQzVveEsPEBCpG7ylrPp62ZcoTanFfsIBgmu66xzJyDxpZPsEPcZ+YEDdDubzvJUuMipF7qQFK9fLcgx0zUTFyX7WqcgnrRIJywliOH/wgDcMAbrihvgu17dgh/t7zz5fvIspCeDpyL0+oxd1L9DsxEXwaJJCPsv0Q99lZcTOUSw2rFLmX99wVy3OfDTbHWT3PnSv67YAQQycVQvfvz2LjRgNdXfWVh44OMRZWqbRDPnKvQ6NCTKjFXYpZpaJYhagSuec99+q/S3qU5QaJlywhnD7NJUuqlrdlrG2UidxL11CpB+ply5gV/XZA2BhObJnBQc6tOFVPiGQ55fLHVdsyzgi1uMsltmZnnX+m3BTtelKLyL2pqXzkzly6/G85WyYSISQSCkXuM8FOPVcpzz2TYQwPO1ufwGnkHlStfECsUTs0VEnctS3jhFCLu4zcZ2bCF7n76bknk0LcK0XuQOlc90or3be0EKYUidxltkxQxONAJKKGLTNhTfqR11M5Egk48tyHhszAyimLipvlrzPZ+9DiXp5Qi7u3yF0Nz1127f1IhcyLe/nIHSgt7uVSIQFrNSZFsmWCtmWICG1tatgyE1bvoVymlCQerzygOjPDmJhAIOvTAnC0ypW2ZZwRcnGXnruz7Zm5bM3rehKJCIHw15Ypvc2SJeJUl47cy5dQbWkhdfLcZ4OvCNjWRkqI+7grcUfFHHIZNXd3B2vLlEtvlrZM0BMRVSfU4t7cLF6dDqgmk0AmE3xdGUlHB/lky8iLvXLkPjw8/++98EIWY2PAxo2lLwf1Ivdg2yDK/gbbBiAfuTu5poXnXn4bGTUHGbmnUuUTDXTk7oxQi7uM3J2uuygjLRVsGcBb4TM7nHju0pOVtXUKue8+MeVv167Sk1bicYKDsbiak2VGMhV85O5mKn8tcRO5O8lzz0fuwYm7aEfpbqJOhXRGyMVdvDqN3FWp5S4RkXv13yNtmXLiLns5divd3H9/BpdeGsGKFaUvh3gcSCkwoCqD5aDFXR1bRrz6ZcvINMTgsmWkuFe2ZbS4lyfk4u7ccx8eNnH99eJOCHKJvUL8s2VMxGIEwyi9X7EYIRabX1I1m2U8/XQWV19dfjZiPA6osDbFJAe7CpNEPVvGnwFVGTEHF7kLSSqX665tGWeEWtzdeO5PPZXFs8+auOOOON74xvpOqy5FR4d/qZCJROVT2dIyP3IfGWEwV14wPBYjpBQYUJ3iYMv9SlpbVYnc3Q2oplIoO1g5OCjmgQS1EIa2Zfwj1OLuxnM/elRc0LffnlBmYV2/PPfZWUZTkxNxp3niLgdYu7rKHxNVIndVbJnOTmBgwMTISLBPvPHcOFLlbaVtlylTl2t4mOtedqAQ2QMpN4dA2zLOCLm4i1cnkfvRoyYikcoRaj1pb/fPlimXKSMR4j73vby4l78UYjEg5aL6Zq2YtdpQLu2zHtx6awKzs8Ctt7qYZFEDJqYZzc1ANOrMlgHKV4acmCi96Es9cLLcnrZlnOFI3InoWiI6SER9RPThEtv8PhE9T0TPEdG3/G2mPZGI8JGdeO5Hjpjo7SXblYaCoqnJWa2PSji1ZZqb59syTj3WeFwNW0ZaxkHf2Fu3RvCe98Txr/+adlVy2m/Gp5xZMkA+0i1XgkAsPB3csTUMsr1OC9GRuzMqms9EFAFwF4A3AOgH8CQR7Wbm5wu22QDgIwB2MPMoEXXXqsHFNDU5j9xXrVKroyKng5smlx0MrcTsrOnClpn7nlNbJhZTw5bJi3ugzQAArFtnYGZG5GR3dgbThokp5+U0pC1TblB1etr5w6JWtLSUH6zWnrsznKjddgB9zHyImVMAvgNgV9E2fwTgLmYeBQBmHvK3maVpaiKHnrt64i6tBSf1PspRzYCqFPdKCz3E42rYMrINKsxOXL5cHLMTJ4Lr0oxPObdRZG+n3PVW77VT7Wht1baMHzhRuxUAjhX83G+9V8i5AM4loseI6HEiutavBlaiubly5J7NMvr7WTlxlz55tdZMMsmOPffiB+HQEGPxYqro2WpbZj7LlonraWAgQFtmkl3bMuVy3aemgh+srrQcoLZlnOFE7ezOdPGRjwLYAOBKADcB+AoRzeuoEtEtRLSHiPYMDw+7bastTU1U0XM/eTKFTAZYvTp4QSjE7SSsUji3Zewid7OiJQOoY8skWZ0bW4XIfWLauS2T99xLbxO05w5IcS/9e9n+WH1XAQwdTsS9H8DKgp97AQzYbPMjZk4z82EAByHEfg7M/GVm3sbM27q6ury2eQ5NTZVL/h47JkLj3t5Gjdyd2jL2qZBOJqyoMkNVJc9dicjdxYCqvN4qi7sfLfOOXRBSSCoFRKOoapzqbMCJ2j0JYAMRrSWiOIAbAewu2uaHAF4HAES0BMKmOeRnQ0vhJHIfGxOJvarMTJX4Fbm7E/e57w0NscPInZSI3FWyZdrbRWXPgYEAPfdJ95F7KVuGmcuuyFUvKkfurMTDXXUqKgIzZwDcBuBnAF4A8F1mfo6IPk5E11ub/QzACBE9D+BhAP+TmUdq1ehCnHjuExNiCSFVaspI/Ircq7NlnE1aEQOq5Wc31gOVBlQBYPlyAydOBHdMJqb9G1BNJgHTROCRu5MBVS3ulXE0D5+ZHwDwQNF7Hyv4PwP4gPWvrjQ1Ec6cCae4+xe5s6vInZlBRGBmnD7NFTNlgPzNlAEQpNWpUuQOiElxQUXuaatCptP1CSp57lJQ1fDcy4u7kwSCsx21TGgPOPHcJyaELaNKHXeJjD798dydZMuIyEze3DJSc1KEKxYT35+upqE+oJLnDgQbucvLptzauYXk89xLL7UIBC/udvZhIdqWcUboxb25ubLnrm7k7m4lqVK4mcQE5CtDyujIiceai/q8NdE3VMqWAeSan8FE7kmXFlWl8gPSsgtDnrsqPTeVCb24O5mhOjGRRSymXlcuH7nXb0AVyN/EUuSdRGoy7SwdtOduvUbVKOyZK/1bq7GIct8rYwLn4i5eVY/cW1vFfAyzRHaW9tyd0QDi7ixyVy1qB/yL3N3MUAUKxd15pJYbjPPWRN9IAYjHxCLVKtDWJqwtP2oEFfN8NgtjfBwPlSjjKP+k06AlTJ47MH/tAUkyqW0ZJzSAuDvx3LPK+e2AP5F7JsPIZuHJlpGvTmwZGbkHLu7MSCh0Y0shqkVt9//PUuH70vYjHUmXFTIr5bnnxd1FI2tAcRBSjLZlnBF6cXfmuWcaNnJPWve90wFVwGvkLl5VsGXiMXXOpcxUKVd/3AszzPieJeqjJY6518i9VDAhbRkV8twBlMx11wOqzgi9uDc1Aem0qB9TClVtGT8i92RKDqrVekBVDVsmyYy4In47kI9yyw0AeuHZbDYn6s9ns7bbeB1QLRW5y4d90JF7Xty1514NDSDuThYgUFPc/YjcZ5Pyu+o0oOq+ib4iPXdVyEfu/or7kCXcOyIRHDRNmDbRu9vIvVLJX3U8d/FaXtzVu59VI/TiLtdRLee7N7Lnnrdlaj2gKl6DLvurmi1TyULwyrB1nK+IRjED4IjNcZ91HbmL10q2TNDiLoMQbctUR+jFXUa/5SY9qBq5y5uyKs+9ClvG3YCqGraMagOqtY7cr4hEAAAv2Fgz+cjd2XdGIgBR5chdhTx3oNKAah0bFFJCL+4rV4qb6/Dh0hNJJifVFHfDEMsEVpNGl7RuVKe1ZYD8TexpQNV1C/1FNVumVp77sGmiCcAmS9wHytgyTmeoEpGoEVQiz316mpFIIPClKJ157urdz6oRenHftMmKbF6wH3RiZmWzZQDnywSWYjYXuVfev4ULCe3twLPPigehpxmqAdsySQBxB4tB14taZcsMMaObCAusfP4xO3H3UERNLu1ohwoLdQA6W8YvQi/uK1cSWluBF16wj9xnZ01ks+rVlZEkEtUtki0jdye2TDRKuPrqKH7603SuvCuRM3FQprYMs1KRuxR33yN3ZnQZBtogVsuxFXfr1c3M63i89PWmQi13YH4Psxhtyzgj9OJuGISNGyN4/nl7cVe1royk2sg9ZaltzOEg4zXXRHHkCOPFF01MT4ub2clsT1Vqy6g3oCpea+G5dxPBIEIH/Ivcy9kyU1PB57gDlTPgtC3jjNCLOwBs2mSUtGVUF/dqI/eMld/vXNxF2PvggxlXCzPkZqgqkS0TaBPmEIsJH9tvW2bYNNFlPXQXENmKe762jJvIvVw9d3Y827WW5MZ3bFaHETarGj0M1WkYce/vFye9GNXFvdrIPWM90yotcC1Zs4YQjQLHjzOmp52vdC8jJRVsGZWyZYDK9cfdwsy5yB0Q4n7GN8+dSop7Oq1GRJwLJGzaOTYm3u/paQjpqikNcYTWrhW70d8/35qRtdyd1CwPgkSiumwZt+JORFi0iDAywpiaYveRu5dG+kgKag2oAuLa8tOWmYKIyrsMcV2XitzdpkICIiouleeuykClYRAiEXtxl+WVe3rUugZUpCHEfcECcaLHx226rrPiYmhuVvNiEFUtqykcJl5jMeenUoi76cqWUSZbRrEBVaDymp9ukTnulWyZJETpYzepi83NpduaSuUf4kFTamxgaEi852RR97OdhhB3ufr72Nj8i0FeICpEJHZUH7mL/XMauQPA4sUicg+lLQO1PHdAZMz4GbmfNkVAsrhQ3G22SzIj4fJY9PSUXlxEpYHKeFzYRMUMDorjrG2ZyjTEEZKR+5jNHZBMigtZtYU6JE7q0ZfDrS0DCHE/fVqkQjrNa1bKllEoWwaovHKQW+Rk61YHkXvCpRj39BBOnrRvazqthi0DiIdMOVtGR+6VaTBxt7kBku4HneqJiNzrN6AKeI3cxWvgJX8VHFD1O3Kfto6xPDVS3ItXZZr1ELkvXWpgeJhtq6iqlD8ei9nbMoODDCJgyRIt7pU4C8S9wSN3y3N3I+6hH1BVzJbx23OXkXuLjNwBZADMFG0nInd3393TQzBN4NSpUuKuxn1SKmVzaIixeDG5ut7PVhpC3NvbxUxLuwHVVEqIuyoRSTHVR+7ePPeZGWBkxLm4ExGiFKy4M7Ny5QcA/7Nl7CJ3YP5EJi+2zNKl4paX3nUhqRQrNKBqb8sMDpraknFIQ4i7YYiaKWG0ZYLx3MVpn5hwVwEwTsHaMlkADPUi97Y2wpkzXHLmp1vmRe6lxJ0ZTR4idwA4eXL+oKrIc3f3fbVCDKjaZ8voNEhnNIS4A8KaCaMtE5TnLnEz3TwWcOSesrJIVBtQ3bkziqkp4O67/Tk6uci9krjDS+QutreP3NW3ZQYHGd3dDSNbNaVhjlJlca93i5zR2UmYmLCPUpzgRdwXLcpvu2GD80sgblCgqZAyx161yP2aa6J47Wsj+OQnq8hpLcCxLeNxQBWAbcaMKpOYAFHWoVS2jLZlnNEw4t7RYS/usqusipdYzOrVBkwT6O/3KO4Zb5675JprnC9IGqNgJzHJyF21bBkiws6dUfT3M2bKrOXrlGmIG1PuZk7ci7abhfvIva1NrF42ODjfllF9EpOoK5NPoNCUp2HEXUTu899PJk3E4+So8mEQrFkjTsErr5RebKQc1doyMpJzQtwI1pZJSltGsQFVADmrYNhjD6yQaWa0IF+ts8N6HbeL3F0+6IgIS5eSrS2jSm0ZwH4SUyoFmGbwK0WFhYYSd7tsmWTSdFTrPCj8EndrwR5HyBzhq6928SEAcdK2TCmkVTCU9nYeC5lG3m8HAFmoMWnnuXs4Fl1dRm4av8Q0GZmMSgOq820ZN8tCagDnfXLFKeW5p1KstLj39hIMozpxF2tjOr/gEwnCvn1tWL/e3XFRxZZRWtx9yJiRkbukyTq3xUlVSQBNHhIF2trmz6iVUbIq4m43iSm/LKQWdyc0kLiXnsSkSlfTjnicsGIFeRf3jDtLRnLhhe6idkDYMmpE7uqdT2nLDPllyxQ8rGUuwLzI3eNs3dZWwujo3OtNRslO1wWoNXbZMqos4B0W1A1pXbJggVj0ojitUHVbBhCDql7FPZ3hus3WazEIkzpyt6Wry8fIHbAVd7vIPeFBjFta5s+oldlaqkTu8TjN89ylLaPCOq9hwJHqEdG1RHSQiPqI6MNltvs9ImIi2uZfE50hK0Nu2jSB0dH8DZZMqm3LAMJ3r8aWqZe4L44RRgIU91Gr1sKCNvVu7tZWkYUy7IfnXmTLEBHiyNdvl8x6jtzzFodERsnqiHs5WyaIFoWPiqpHRBEAdwHYCWAzgJuIaLPNdu0A/gzAE3430glykPDwYcZvf5tfci+VMpWdwCTp7TVw/Pj8wlBOqKe4L4kRTgUo7kNWKNezWL3zSUTo7iZ/bBnkK0JKmiDEvBAvtWUA+8g9L+5qHFvhuc99Tw+ousNJSLsdQB8zH2LmFIDvANhls93fAPg05vce68KuXTF8/vMir+Do0Xz0JDx3tSP3jg4gmwVmPQR9dY3coyJy9/IQ8oMh627vXqTm+ezuNmoyoAqIQdXCG4uZkYI3W6a1lWwid7Xmg9hly+Q9dy3uTnByl6wAcKzg537rvRxEtBXASmb+Vx/b5ormZsItt4gwZq64s/KRu1zfddLDBJhMhl2twlQNS2IGksjXPqk3Q+k04gAWKLpkoojcfbJliiL3BOYOqErd8xa5iwViCsv+qpYtY1dbRtsy7nCiCnbKmDvqRGQA+ByAOyp+EdEtRLSHiPYMDw87b6VDmpsJXV2EY8cKPXf1B1SluE9kK2xoQ709dwCBWTODqRS6Sd0JaV1d/tkyxeJeHLnnluJb6P7algOShdaMarZMuTx3PaDqDCdXRj+AlQU/9wIYKPi5HcAFAB4holcAXAZgt92gKjN/mZm3MfO2rq4u760uw6pVxpzIXfU8d6BQ3NX23KW4BzWoOpROo1tRYQeEuA/XyJYpjtyPWZlDKz1USJSRb6E1o9pylNJzL7QAdeTuDieq9ySADUS0lojiAG4EsFv+kpnHmHkJM69h5jUAHgdwPTPvqUmLK7BqFdl47uoKAiAmlQDqi/sSKe5m9daDF4bSafQY6j6oOzoISa6uLDIzO4rc+61z0OthLdF85D7fllHHcweYxViURA+ouqPilcHMGQC3AfgZgBcAfJeZnyOijxPR9bVuoFtWrhSRu3zih8mW8eq513NAFQjellEVKTrFKya5IQ1Rt75S5N5v/d+LuMt2ThcMnqhoywBzM2ampsQSe01NJT6kmYOjGarM/ACAB4re+1iJba+svlneWbXKwOSkWCy7szNstoz7z54ttgwzYyiVQrcqoaUNObuDOVfsyy3FC3VImjD3odFvmmgB0Nnu/m+0torXwshdNVsmt2ZvwUQmueavqmMuqqG26nlg5Upx4vv7Rbc1TJH7REZtW2ZRgOI+kc0iyay0LZOLiKv4juKFOiRNRHPy3PuZ0WsYnoRO2jL2kbvrr6sJuTV7U4Weu7Zk3KDuneKRzk6rPOp43pZRpatZinYr+vLuufvcoBJEidCJYGyZQZnjrnDUlhP3Ko5P8UIdkgTmzlDtN02s9HgsZA/D3nNX4/ja2TIyctc4o+HEvc2amj4xER7PXbbZk+eerZ/nDgBLDCOQyH3YUp8upcVdvFYl7vK7HAyo9nrsxeQjd/VtmWLPXUfuzlFb9TyQE8pJ8bPw3NW+IOJxQjzu0XP3WBXSK4somPoyU1baRLvS4l69LTNoZcF0lpnElGXGgGXLeEG20z7P3dNX+k7ec9e2jFcaWNwZpslIp9UfUAWE7656KiQAtKK6yNQr05boFUe0KtHcLF6rOT6PZ7MgAJcUrb5SGLmPMSML770YOaA6N3IXr6pYmNIeKrZlZNs1lVFf9Vwic8YnJ7ngglV/N9vbvdoy9RX3ZqKqUv28khP3AP62U/xIhfxNNotNhpFbN1WSQL5wmNUp9dyLsYvcZYSsSjKSnS2jI3d3qK96LsnljE+K+hkAlLdlANHj8Ba519dzbwYwE6Ato3LkXq3nzsx4PJvFq2zWTGwiyg2oypr6XoNY2cOYmwopXlWzZQqzZYTnHlCDQkjDiXtTE2AYMnIXF0ZobJmM+8/V23NvJgpE3KctcVe5V15ttsxLMzMYZcarbNKfmiAmOJnMOXFv8/igIyK0tKhty8h2zM9zV6N9YUB91XMJEaGtTYh7PnJXfze9eu7peos7qrMdvBIGz73aAdUjs8JVP9dmoDS31B4A6aZ4FXdAZMzMHVBVP1tG2zLuUF/1PNDWRpiY4NySe6pEI+UIi+feEmDkTsiLnIpUa8vIlaYW2oh24SLZ1UbuAOZF7qrVlrGfxKQHVN3QsOI+OZl/6ochcg+N504USD33adNEi8cZmfUiFiNEyXvkXk7cCxfJrtZzB2TkPteWiUbVmdpfPIlpbIwxMZFfcU1TGfVVzwPhtWXcf67u2TKwilvVOXqfNk202Aw0qkaL4T1yP+M0crfeqy5yp6LyA6yMJQPMry3zxBPi2GzfXqfp2A2A+qrngfZ2ssRdDqiq/7T3nOeeQd1WYgJE5A7U33efzmbRonBdGUlLxHvPZjSTQQziAVqMjNxnmTHlgy2zfDnh5ZcL1z1Qx28H5mfLPP54FkTA9u3qP+BVQf27xQPSlpGRexjy3JubgQy7j4iDiNyB+qdDhidyp6o894UlVpqSkXsS1adCAsCOHVG8+KKJ4WEh8Om0WmNTxZOYHn88i/PPN9DRoU4bVUd91fOAEHcuWFBX/d3MeYwuPxfEJCZAR+6laIlUZ8vYWTKASIUEROQ+CSAOIFZF5L5jh3hQPvaY8AJTKVZmMBWYmy3DzHjiiSwuvVT9h7tKqH+3eKCtTRQOGxwUN1lPj0L9zRIkrH53svxm8xADqr43pyQ6ci9Pi1GdLVNK3BNFkXs1lgwAbNsWQSIBPPqo8LJVtWXSacbsLHD6NGP9+oaUq5rRkEdLRu5hEvdcpKK6LaMj97JUM6A6msnMKxgmKYzcp5jR5rF9kkSCsHVrBHv2yMhdLVumMFtmdFQcz4UL1WlfGFD/bvGAGFAFTpww0d4eQWur+hGfHPR1HbkHMIkJ0JF7KaoaUE2nS0fu1qvMc2/1IWVx2TLCyIg4j+m0WtkyskTC5CTnxF2u1aBxRkOKe1sbwTSBV14xQxG1A+GL3Oud6z4Vksi9ucpUyJKee6Etg+oyZSQLFxJOn7aKkU0Czc3qiGc8Lmaanz7NOHNGR+5eUP9u8YCsDNnXZ2LpUoVGicrgxXM3mcF89kTurWGJ3D0cG2bGmTK2TGEq5KQPtgwgxFJGxadPMxYvVks8Fy8WPQtty3ijQcVdXAQvv2xi6dKwRO7us2VknbF6lx8AtOdeihbDW0nkiWwWWdhPYALykfsMIDx3nyL3mRkgmWSMjJjKi7u2Zdyh/t3iASnus7PhGEwFCiJ3F1FfEOKuPffyeB1QLVd6AMgvzDFkmpgEfPHcFy0S3zE6yhgZUTFyN3TkXgUNKe7r1uV3S0fu/pLLlqmjuKeZkWEOR+QeIUxB2CxukKUHSmbLEGExEfp9tmUAYGiIMT6eF3tVWLSI5njuOnJ3h/p3iwcuusiAvEfCIu5ePPeMJSCNngqZWzQ6BJF7b4JgAjjqs7gDQC8RjpumL3nuQF7cZRkC9SL3vC3T1lbf67wRaEhxJyJcdpkQgvZ29QUB8JYtc7bYMtLmCEPkfqGVdrs/664K3IyDZQR7DQP9pokp+JctA4jEA0BNcZeWkbZk3KP+3eKRj39cTPu4+GI/OrC1x0ueexDingBAcB65TzDjqGlW3rAMOXEPQeR+gUdxn7WOUVMZ0V5hGOgzTZjwZ0WqMETuzMDhw6YWdw80bP3M178+CtPsAFE4Fl30ErnLFcjqKe5EhCY4GzQ8mM1i4+QkEgBmOjo81wrP2TIhiNwXRAmribDf5QMtJ+5ltuklypX77fHhWEiPva8va/2s1vGVD5u+PhPnnqtW28JAQx8xVRYecIKnyD0Azx2w1lF1sN3nrJJ+SQCjVlu/nUqhZ3wcAy7EL0y2DABcGInUJHLvLdj/K3woKCQHKF96Sc3IXT58Tp7UtowXwnG3nAXkIncXnwnClgGsdVQdRO7HCwT8iLX9vek0hpjxP2acD8mOWJ9dUM8KaVVwQSSCA6aZe/g6QYq7XS13yYoC4V/rQ+ASiRA6OoCjR0U7VRP3wvZocXePFndFCEueO+A8ch8wTXRbInTEEi/5uf/IZHDaYfT+vLXdxpZwWGxLiZABMO5B3MtF7sutyL0Z/vVKpWhGo/mZ3apQKO46DdI9WtwVoZo893quxAQ4j9xPMOMyaxD0zmQSl09O4uUCQR90KH7PZbNYFo9jkUoFx8uwwBLeMy4+48RzX28YOM8w8F0fH3JS3Bcvtl8kJEgWL85f14VzVzTOCEc/9yygusjd//aUo8VB5J5lxiAzLoxE8B+ZDJ4o8KC3RyL4z2wWJ5mxycHfe9Y0cUFHR1VtricyV33MxbmUqZCJMts0EeFAe3s1TZtHV5cBwMxdfyrR2Qncfnscq1cb+OM/Dsd8FZVw9DgkomuJ6CAR9RHRh21+/wEiep6I9hHRg0S02v+mNjaePPeABlQ7iXI+eCmGmWECWEaE1UUDoa+yovlBB7aMyYznsllcEBJLBiiI3F3aMokSS+zVkr/7uybcfHMM73ufeupORPjc55px++0JpWrNh4WK4k5EEQB3AdgJYDOAm4hoc9FmzwDYxsxbAHwfwE+hjE0AAA+/SURBVKf9bmijYxiEKKmf5w4A6wwDhyoI8wlL2JYZRklxP+lA/A6bJmYAXNDqR2Z3fVjgIXKfNU00BZANtGVLBPfc04IPflA9cddUh5OraTuAPmY+xMwpAN8BsKtwA2Z+mJllOvLjAHr9bebZQZzCMaB6jmFglDmX3mjHCUv8lxFhlSV2my3xujgSQQzOPPc+63vODVHk7sWWCUrcNY2LE7d2BYBjBT/3A7i0zPbvBvCTahp1tpIwKBSpkOssETpkmrikxKzRwsj9D+PxXI72Z5JJrDEMdBPhpANbRlobi0OSBgkAC6xXt7ZMcwhm4GrCg5M7xk45bK9aIno7gG0Arijx+1sA3AIAq1atctjEs4e4S1smFZDnfo4l1C9ns6XF3RLupURYE43iVdEoZpnxjngcMSIsNQxHtsyZkOW4A+GyZTSNi5OrqR/AyoKfewEMFG9ERK8H8FEA1zOzrUYx85eZeRszb+vq6vLS3oYmYbgrPzBlvba11TfiK4zcS3GSGQuJ5uRtNxFhjfXZpUSOBlTHrNcwiXuUCK1wH7lrcdf4iZOr6UkAG4hoLRHFAdwIYHfhBkS0FcA/QQj7kP/NPDuIG+Qqcp+0xKPe4t5GhG6iOTnrxYwz5+wJO3qI8Ixp4rLJybLjDGPMiCI8pQckC4jcR+6K5Zlrwk3FO4aZMwBuA/AzAC8A+C4zP0dEHyei663NPgOgDcD3iGgvEe0u8XWaMiTIZeQekLgDlTNmppjLrhYUtX73RDabGzS14wwzFgSQIlgtnV7EPWQPMI3aOOrrMvMDAB4oeu9jBf9/vc/tOiuJG+48d1khMAhx7yJCfzlxR/ml4N4cjeLLVmGxPtPE+SW8+zFL3MPGAiKcYcZh08Q9qRQ+kkggXmY/ZkxTDBpXWR5Zo5HoUEEhEuQuW2aSGURAU1P9T2M7EcbL/H6KuWzN8TfFYhixZluWs3fGKtg7qtJJhDEAX0om8dfJJP6kQqE0Hblr/EZfTQoRN9zluU8yo60lmNLGHUSYKNPWKWa0VGjXIsNAJ1DRlim39JyqyMhd5vJ/LZ3OZRDZocVd4zf6alKIuIfIva05GOFrR/mqh5VsGcn6SKRy5B5ScR9jnrNv5Wwsneeu8Rst7gqRcBm5TwFoawlG+DqIMAsgXaK9lWwZyTnW0nGlCKu4dxaI+4VWRD5Q5tzqyF3jN/pqUoi44a5wmLRlgqDDEtxS1kylbBnJesPAEdMs+ZAIqy2zyOqFnWDGa60c/Yq2TAj3U6MuWtwVIkHu89xbg7JlLCEqNajq1JZZSYQs7OvMmMyYAEIZub++YNLVqyIREPIlGezQkbvGb/TVpBBxlzNUJ4HAPPdykXuaGWnAkS2zoMz3TEDUuQijuG8tEOrzDANdRCUj9ywz0sxa3DW+oq8mhUi4zHOfCtCWkUtG2A2qyrIITiJ3+ZCw+55cXZkQijsR4fetlaPWRyJYTlQycpfnXIu7xk/01aQQnrJlAhxQBewjbjlz1knkXk7c5QzPMHruAPDN5mY819aGhURYZhgYKBG5z1r7qcVd4yf6alIIt9kyk0Dwnns5cXcTudv8bizEkTsAxIiw2UpvXFYmcpfTm3QqpMZPtLgrRJjy3MuJsl+2zFEr0u0OqbgXsswwcILZtt6Mjtw1tUBfTQqRMIA0RJZIJdLMSAHBee5+2TLWq524/zyTwSIinN8AorfSOl6d4+MYLrJnZq1XLe4aP9FXk0J0RKxFHhxsm6vlHpDnXnZA1YUtU+ohwcz490wGr49GEWmAyP2meBwftFZBP1BK3BtgPzXqoMVdIRbHxM192kFlwFwt94BsmQgRWlAicrdenYh7jAhNmP+QeN40McCMN4ZokY5ydBDhjyxxP1ws7tqW0dQAfTUphBT3EZso9niRIEhxb22uT9vs6CAqH7m7+Z6i957NZgEAlzbQIONqS7xfKTqX01rcNTVAX00KsThqL+7/c3YWvRMTcwQ+F7kHZMsA9qIMuLNlct9TtM+nrJ8bYTBVkiDCcqJ54j4s99WK7DUaP9DirhB2kfsoMz5rLWqxx4pmgeA9d0D45dXaMoAYVC0Wdyl4ixpI3AFgjWHMs2XkQuFLtbhrfESLu0LYifv30+nc/39bIO5y9mZ7QNkygBBlu9Q+T7aMTeS+iCi3HF+jsNYw5kXuJ626Mh0NZEFpgkeLu0J0RgmEueJ+1DRhQKxZujebRZ8l8LJO+JrlwZ3CFTZRKCDEPQ44FuZS4r6kwYQdEJH7MeY5k9VOMmNpPB66dWI1aqPFXSEiRFhINEfcB0wTS4lwsWHg/kwGGyYn8VI2i5dMEwuJsLgzOEHYEolggHledo/TipASO3EfNs2GFPd1hoEsgK7xcZy0jttJZiy16tBoNH6hxV0xFheJ+wlmLDcMrCvIpHgym8VL2Sw2GEag0Z5chGJ/kbiPulz31G5g9hQzuhpQ3N8Wi+GOeBwTyNtsJ01T++0a39HirhiLbCL35US4OR7HNVbO914rct8QcOrchZZHvL9gLAAAjpgm1rhom4zcuWC/hxvUlmknwh2JBID82rHSltFo/ESLu2IsJsJIQSQ8YEXuF0Qi+GlrK7YaBp7IZnGMOXBxX06EdgDvm53F3al8VZzDbsUdouyCLH3LzMJzb9C876XWBLCXrRWoTmlx19SAxrx7QsxiIpy2ItgUM4aZsbwggv2dSAS/zGbBAM4NWPyICBut6P09M6K24SwzBpix1mXkDuQzbyYgxL4RbRlAHDe5duwwMxg6DVLjP1rcFaPQc5f5z8sLhPKignS5CxVInbunuRm/a9lFp00TR6xehxtxX1k0c1MW1mpEW0ay3jDwsmnqHHdNzdDirhhdRJiEmIEqF3cojNz/ayyG34/FcH9LCy5QQNw3RiJ4ryVMz5lmTqDd2DKbrG2ftz4rZ6c2srifY4n789Z4Ra/lw2s0ftEYVZkaiK2WYP9nNpubqFQYua80DNzbEuDMJRvkghTPW3YR4C5yX2sYiAN4wRK64bNA3NcbBpIA/j6ZxHIibG1rC7pJmgZDi7tivCoaBQF4LJOBjMt7FRe5VURog4jcmwHEMLe3UYkoEc41DLxgRe57s1kQEPiAcS25MhpFHMAzpon3x+MwFD/HmvDRuHdPSOm0Fqd4LJvF99JpXBaJYLHiIkfWcnL7sln8eyaDzYbhWqw2RyI5W+aX2SwuMAwsUny/q+G8SASfb2qCAeBm7bdrakDj3j0hZkc0iocyGew1TdwYkpmLb4hG8YtsFntNE7d78I83WaUMxpnx60wGlzdIHfdy3JpIYLi9PWfFaTR+osVdQW6Jx7HEWsTibSER948lErgyEsEmw8AfeGjzqyMRMIA/mZnBFIDLzxLBa+TeiSZYGj88CiEXRyJ4pb0dp6wJTGEgToQHW1uRhFhdyS1viEbxhmgU306n0UOEq8+CyF2jqSWOlIOIriWig0TUR0Qftvl9gojutX7/BBGt8buhZxtxotAIu8QgQrPHgUEiwpebm3F7PI6n2tqUH2fQaFSn4h1ERBEAdwHYCWAzgJuIaHPRZu8GMMrM6wF8DsCn/G6opvFZYxj4XHMzVmhh12iqxsldtB1AHzMfYuYUgO8A2FW0zS4A/2L9//sAriZdnFqj0WgCw4m4rwBwrODnfus9222YOQNgDMBiPxqo0Wg0Gvc4EXe7CLx4bTUn24CIbiGiPUS0Z3h42En7NBqNRuMBJ+LeD2Blwc+9AAZKbUNEUQALAJwu/iJm/jIzb2PmbV1dXd5arNFoNJqKOBH3JwFsIKK1RBQHcCOA3UXb7AbwTuv/vwfgIS5ceUGj0Wg0daViMjEzZ4joNgA/AxAB8FVmfo6IPg5gDzPvBnA3gG8QUR9ExH5jLRut0Wg0mvI4minCzA8AeKDovY8V/H8WwNv8bZpGo9FovKITijUajaYBoaCscSIaBnDE48eXADjlY3OCRO+Lmuh9URO9L8BqZq6YkRKYuFcDEe1h5m1Bt8MP9L6oid4XNdH74hxty2g0Gk0DosVdo9FoGpCwivuXg26Aj+h9URO9L2qi98UhofTcNRqNRlOesEbuGo1GoylD6MS90sIhqkNErxDRfiLaS0R7rPcWEdF/ENFL1uvCoNtpBxF9lYiGiOjZgvds206C/2Odp31EdHFwLZ9PiX35ayI6bp2bvUR0XcHvPmLty0EiuiaYVs+HiFYS0cNE9AIRPUdEf269H7rzUmZfwnhemojoP4not9a+/L/W+2utBY1eshY4ilvv+7/gETOH5h9E+YOXAawDEAfwWwCbg26Xy314BcCSovc+DeDD1v8/DOBTQbezRNsvB3AxgGcrtR3AdQB+AlEx9DIATwTdfgf78tcAPmiz7WbrWksAWGtdg5Gg98Fq2zIAF1v/bwfwotXe0J2XMvsSxvNCANqs/8cAPGEd7+8CuNF6/0sAbrX+/6cAvmT9/0YA91bbhrBF7k4WDgkjhYud/AuAGwJsS0mY+ZeYX+2zVNt3AbiHBY8D6CSiZfVpaWVK7EspdgH4DjMnmfkwgD6IazFwmPkEMz9t/X8CwAsQ6yuE7ryU2ZdSqHxemJknrR9j1j8GcBXEgkbA/PPi64JHYRN3JwuHqA4D+HcieoqIbrHe62HmE4C4wAF0B9Y695Rqe1jP1W2WXfHVAnssFPtideW3QkSJoT4vRfsChPC8EFGEiPYCGALwHxA9izMsFjQC5rbX9wWPwibujhYFUZwdzHwxxJq07yWiy4NuUI0I47n6IoBzAFwE4ASAz1rvK78vRNQG4AcAbmfm8XKb2ryn+r6E8rwwc5aZL4JYA2M7gE12m1mvvu9L2MTdycIhSsPMA9brEID7IU76oOwaW69DwbXQNaXaHrpzxcyD1g1pAvhn5Lv4Su8LEcUgxPD/MvN91tuhPC92+xLW8yJh5jMAHoHw3DtJLGgEzG2vowWP3BA2cXeycIiyEFErEbXL/wN4I4BnMXexk3cC+FEwLfREqbbvBvAOKzvjMgBj0iZQlSLv+S0Q5wYQ+3KjldGwFsAGAP9Z7/bZYfmydwN4gZn/vuBXoTsvpfYlpOeli4g6rf83A3g9xBjCwxALGgHzz4u/Cx4FParsYRT6OohR9JcBfDTo9rhs+zqI0f3fAnhOth/CW3sQwEvW66Kg21qi/d+G6BanISKNd5dqO0Q38y7rPO0HsC3o9jvYl29Ybd1n3WzLCrb/qLUvBwHsDLr9Be16DUT3fR+Avda/68J4XsrsSxjPyxYAz1htfhbAx6z310E8gPoAfA9Awnq/yfq5z/r9umrboGeoajQaTQMSNltGo9FoNA7Q4q7RaDQNiBZ3jUajaUC0uGs0Gk0DosVdo9FoGhAt7hqNRtOAaHHXaDSaBkSLu0aj0TQg/z9ohPHmxRoF0QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(L):\n",
    "    if viterbi_path[i] == 0:\n",
    "        plt.axvline(x = i, color ='yellow')\n",
    "    else:\n",
    "        plt.axvline(x = i, color = 'red')\n",
    "plt.plot(range(L), post_prob,color = 'black')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 226,
   "metadata": {},
   "outputs": [],
   "source": [
    "d1 = 0\n",
    "d2 = 1\n",
    "states = (d1, d2)\n",
    "hidden_n = 2\n",
    "obs_val = 6\n",
    "eps = 0.001"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вспомогательные функции для алгоритма Баума-Велша"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 227,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getAlpha(seq, transit, emis, PI):\n",
    "    T = len(seq)\n",
    "    alpha = [[0.0 for _ in range(T)] for _ in range(hidden_n)]\n",
    "    for i in range(hidden_n):\n",
    "        alpha[i][0] = pi[i] * emis[i][seq[0] - 1]\n",
    "    for t in range(1, T):\n",
    "        for j in range(hidden_n):\n",
    "            for i in range(hidden_n):\n",
    "                alpha[j][t] += emis[j][seq[t] - 1] * alpha[i][t - 1] * transit[i][j]\n",
    "    return alpha\n",
    "\n",
    "def getBeta(seq, transit, emis):\n",
    "    T = len(seq)\n",
    "    beta = [[0.0 for _ in range(T)] for _ in range(hidden_n)]\n",
    "    for i in range(hidden_n):\n",
    "        beta[i][T - 1] = 1.0\n",
    "    for t in range(T - 2, -1, -1):\n",
    "        for i in range(hidden_n):\n",
    "            for j in range(hidden_n):\n",
    "                beta[i][t] += beta[j][t + 1] * transit[i][j] * emis[j][seq[t + 1] - 1]\n",
    "    return beta\n",
    "\n",
    "def getGamma(alpha, beta, T):\n",
    "    gamma = [[0.0 for _ in range(T)] for _ in range(hidden_n)]\n",
    "    for t in range(T):\n",
    "        den = 0.0\n",
    "        for j in range(hidden_n):\n",
    "            den += alpha[j][t] * beta[j][t]\n",
    "        for i in range(hidden_n):\n",
    "            gamma[i][t] = alpha[i][t] * beta[i][t] / den\n",
    "    return gamma\n",
    "\n",
    "def getKsi(seq, transit, emis, alpha, beta):\n",
    "    T = len(seq)\n",
    "    ksi = [[[0.0 for _ in range(T - 1)] for _ in range(hidden_n)] for _ in range(hidden_n)]\n",
    "    for t in range(T - 1):\n",
    "        den = 0.0\n",
    "        for k in range(hidden_n):\n",
    "            den += alpha[k][t] * beta[k][t]\n",
    "\n",
    "        for i in range(hidden_n):\n",
    "            for j in range(hidden_n):\n",
    "                ksi[i][j][t] = alpha[i][t] * transit[i][j] * beta[j][t + 1] * emis[j][seq[t + 1] - 1] / den\n",
    "    return ksi\n",
    "\n",
    "def getPI(gamma):\n",
    "    pi = [0.0 for _ in range(hidden_n)]\n",
    "    for i in range(hidden_n):\n",
    "        pi[i] = gamma[i][0]\n",
    "    return pi\n",
    "\n",
    "def getA(ksi, gamma, T):\n",
    "    transit = [[0.0 for _ in range(hidden_n)] for _ in range(hidden_n)]\n",
    "    for i in range(hidden_n):\n",
    "        den = 0.0\n",
    "        for t in range(T - 1):\n",
    "            den += gamma[i][t]\n",
    "        for j in range(hidden_n):\n",
    "            num = 0.0\n",
    "            for t in range(T - 1):\n",
    "                num += ksi[i][j][t]\n",
    "            transit[i][j] = num / den\n",
    "    return transit\n",
    "\n",
    "def getB(seq, gamma):\n",
    "    T = len(seq)\n",
    "    emis = [[0.0 for _ in range(obs_val)] for _ in range(hidden_n)]\n",
    "    for i in range(hidden_n):\n",
    "        den = 0.0\n",
    "        for t in range(T):\n",
    "            den += gamma[i][t]\n",
    "        for v in range(obs_val):\n",
    "            num = 0.0\n",
    "            for t in range(T):\n",
    "                if seq[t] - 1 == v:\n",
    "                    num += gamma[i][t]\n",
    "            emis[i][v] = num / den\n",
    "    return emis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "metadata": {},
   "outputs": [],
   "source": [
    "def twoDimDifMax(m1, m2):\n",
    "    maximum = 0.\n",
    "    for i in range(len(m1)):\n",
    "        maximum = max(maximum, oneDimDifMax(m1[i], m2[i]))\n",
    "    return maximum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "metadata": {},
   "outputs": [],
   "source": [
    "def oneDimDifMax(m1, m2):\n",
    "    maximum = 0.\n",
    "    for i in range(len(m1)):\n",
    "        maximum = max(maximum, abs(m1[i] - m2[i]))\n",
    "    return maximum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 230,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Алгоритм Баума-Велша\n",
    "def baum_welch(seq, transit, emis, PI):\n",
    "    while True:\n",
    "        alpha = getAlpha(seq, transit, emis, PI)\n",
    "        beta = getBeta(seq, transit, emis)\n",
    "        gamma = getGamma(alpha, beta, len(seq))\n",
    "        ksi = getKsi(seq, transit, emis, alpha, beta)\n",
    "        PI_STAR = getPI(gamma)\n",
    "        transit_STAR = getA(ksi, gamma, len(seq))\n",
    "        emis_STAR = getB(seq, gamma)\n",
    "        if twoDimDifMax(transit, transit_STAR) < eps and oneDimDifMax(emis[1], emis_STAR[1]) < eps and oneDimDifMax(PI, PI_STAR) < eps:\n",
    "            return transit, emis, PI\n",
    "        transit = transit_STAR\n",
    "        emis[1] = emis_STAR[1]\n",
    "        PI = PI_STAR\n",
    "    return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 231,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Приближающие итерации по алгоритму Баума-Велша\n",
    "def run_on_seq(seq, transit, emis, PI):\n",
    "    emis_new = [[0.0 for _ in range(obs_val)] for _ in range(hidden_n)]\n",
    "    transit_new = [[0.0 for _ in range(hidden_n)] for _ in range(hidden_n)]\n",
    "    pi_new =np.array([0, 0])\n",
    "    for i in range(100):\n",
    "        transit_new, emis_new, pi_new = baum_welch_train(seq, transit, emis, PI)\n",
    "    return transit_new, emis_new, pi_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 232,
   "metadata": {},
   "outputs": [],
   "source": [
    "transit_new, emis_new, pi_new = run_on_seq(observed, transit, emis, pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 233,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.95576286, 0.04423714],\n",
       "       [0.10228545, 0.89771455]])"
      ]
     },
     "execution_count": 233,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "transit_new #Приближенная по Бауму-Велшу матрица переходов"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 234,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.1726212 , 0.19541371, 0.14307666, 0.1799557 , 0.14178638,\n",
       "        0.16714634],\n",
       "       [0.06147175, 0.13153735, 0.08716328, 0.1233187 , 0.09026933,\n",
       "        0.50623959]])"
      ]
     },
     "execution_count": 234,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "emis_new #Приближенная по Бауму-Велшу матрица эмиссий"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Viterbi = O(L * M^2)\n",
    "Forward algorithm = Forward algorithm = O(L * M^2)\n",
    "Baum Welch = O(L * M^2)\n"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}