TP Forward-Backward.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "PGM_Lab_1.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "QXRWTxuzbttZ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.linalg import toeplitz"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y32y45tscuB0",
        "colab_type": "text"
      },
      "source": [
        "## Question 0:\n",
        "\n",
        "$s = \\log(\\sum_{i = 1}^{I}e^{v_i - v^* +v^*}) = v^* + \\log(\\sum_{i = 1}^{I}e^{v_i - v^*}) $\n",
        "\n",
        "\n",
        "With $v^*  = max{v_i} ,i \\in I$ . **This** was seen in the course. In the end the arguments in the exponential are positive so the log can be applied."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "J_r8zqkhTGG7",
        "colab_type": "code",
        "outputId": "dc979969-e805-41e4-e3ca-6d037411b0d3",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 33
        }
      },
      "source": [
        "def log_sum_exp(v):\n",
        "  v_max = np.max(v)\n",
        "  return v_max + np.log(sum(np.exp(v-v_max))) \n",
        "print('For v = [-1234,-1235] we get: ' + str(log_sum_exp([-1234,-1235])))"
      ],
      "execution_count": 41,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "For v = [-1234,-1235] we get: -1233.6867383124818\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "f5XpC6PXVXHL",
        "colab_type": "text"
      },
      "source": [
        "## Question 1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0L4rPbXTgyt8",
        "colab_type": "text"
      },
      "source": [
        "###1. \n",
        "Let $N$ be the different states.\n",
        "\n",
        "Let $X_k$ be the true robot position after k iterations\n",
        "\n",
        "Let $Y_k$ be the state detected by the sensor.\n",
        "\n",
        "\n",
        "Let $E$ be a variable such that: $P(E=0) = 1- \\epsilon$\n",
        "and $P(E = 1) = \\epsilon$\n",
        "\n",
        "Let $U$ be a random variable following a uniform distribution on \n",
        "${1,...N}$\n",
        "\n",
        "Let $W$ be a random variable such that: $P(W=1) = w$\n",
        "and $P(W = 0) = 1 - w$\n",
        "\n",
        "\n",
        "So we have according to the text:\n",
        "\n",
        "$X_{k+1} = R[(E X_k + (1-E)(X_k + 1)),N] + 1$\n",
        "\n",
        "$Y_k = W X_k + U(1-W)$\n",
        "\n",
        "$R(a,N)$ being the rest of the euclidiant division of $a$ by $N$ ie. we have to work with the modulo\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EpUGLZCPVzF9",
        "colab_type": "text"
      },
      "source": [
        "### 2.\n",
        "The conditional probability table is given by $A_{i,j} = \\mathbb{P}(X_{k+1} = i|X_{k} = j)$ so with 1. we conclude that it is of the following form :\n",
        "\n",
        "\\begin{pmatrix} \\epsilon & 0 & ...&1 - \\epsilon \\\\\n",
        "                1-\\epsilon & \\epsilon &...& 0\\\\\n",
        "                .&1-\\epsilon\\\\\n",
        "                .\\\\\n",
        "               .\\\\\n",
        "               0& .&.&\\epsilon\\end{pmatrix}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "h1T35FklWnq4",
        "colab_type": "text"
      },
      "source": [
        "### 3. \n",
        "We write M the number of total observations.\n",
        "* (a) $\\mathbb{P}(x_k|y_{1:k})$\n",
        "* (b) $\\mathbb{P}(x_k|y_{1:M})$\n",
        "* (c) $\\mathbb{P}(x_{k+1}|y_{1:k})$\n",
        "* (d) $\\mathbb{P}(y_{k+1}|y_{1:k})$\n",
        "* (e) $\\mathbb{P}(x_0|y_{1:k})$\n",
        "* (f) $argmax_{x_k}\\mathbb{P}(x_k|y_{1:k})$\n",
        "* (e) $argmax_{x_{1:k}}\\mathbb{P}(x_{1:k}|y_{1:k})$\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RvrbwCHYYcpi",
        "colab_type": "text"
      },
      "source": [
        "###4."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "K1aYtShJmuGb",
        "colab_type": "code",
        "outputId": "9af6903e-28fc-4d10-f0b8-af89c83af3b0",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 283
        }
      },
      "source": [
        "def simulation (w,eps,N, n_iter):\n",
        "  X = np.random.choice(np.arange(1,N+1), p = np.ones(N)/N)\n",
        "  List_X = [X]\n",
        "  W = np.random.choice(np.array([1,0]), p = np.array([w, 1-w]))\n",
        "  U = np.random.choice(np.arange(1,N+1), p = np.ones(N)/N)\n",
        "  y_0 = w*X + U*(1-W)\n",
        "  List_y = [y_0]\n",
        "  \n",
        "  for i in range (n_iter):\n",
        "    \n",
        "    E = np.random.choice(np.array([1,0]), p = np.array([eps, 1-eps]))\n",
        "    X = (E*X + (1-E)*(X+1)) % N \n",
        "    W = np.random.choice(np.array([1,0]), p =np.array([w, 1-w]))\n",
        "    U = np.random.choice(np.arange(1,N+1), p = np.ones(N)/N)\n",
        "    y = W*X + U*(1-W)\n",
        "    List_y.append(y)\n",
        "    List_X.append(X)\n",
        "    \n",
        "  List_X = [u + 1 for u in List_X]\n",
        "  List_y = [u + 1 for u in List_y]\n",
        "\n",
        "  plt.plot(np.arange(len(List_y)), List_y, label = 'evolution des observations')\n",
        "  plt.plot(np.arange(len(List_X)), List_X, label = 'evolution des vrais états')\n",
        "  plt.ylabel('états')\n",
        "  plt.xlabel('itérations')\n",
        "  plt.legend()\n",
        "  plt.grid()\n",
        "  plt.show()\n",
        "  \n",
        "simulation(0.8,0.3,50,100)"
      ],
      "execution_count": 42,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXmcHGWd/99PVV9z554kJCHhJhcJ\nuTgUBlBWEQV/IrorEm5c12tXWdT1xkUUBQRRuUQOV1lARAEPZDMcEpIQcpALQi4yuTMzmbuvquf3\nR3VVV18zfUxPd88879crr8n0dFc/NdP1fOv7/XwPIaVEoVAoFCMXrdQLUCgUCkVpUYZAoVAoRjjK\nECgUCsUIRxkChUKhGOEoQ6BQKBQjHGUIFAqFYoSjDIFCoVCMcJQhUCgUihGOMgQKhUIxwvGUegHZ\nMG7cODl9+vS8XtvT00NNTc3gLqjMUec8MlDnPPwp9HxXr159WEo5fqDnVYQhmD59Oq+//nper21u\nbqapqWlwF1TmqHMeGahzHv4Uer5CiF3ZPE+FhhQKhWKEowyBQqFQjHCUIVAoFIoRjjIECoVCMcJR\nhkChUChGOEXNGhJC7AS6AAOISikXCiHGAI8B04GdwKVSyvZirkOhUCgUmRkKj+AcKeU8KeXC2Pdf\nBV6QUh4PvBD7XqFQKBQlohR1BBcBTbH/PwQ0AzeWYB0KRdkRCfaw6smfoPUdAUB6A8z56A3U1o8u\n8coUwxlRzJnFQogdQDsggXuklPcKIY5IKUfFfi6Advv7pNdeB1wH0NjYuOB3v/tdXmvo7u6mtrY2\n31OoSNQ5Vy6T1t/FiW1/x5QCTVjX5tNH3UDD8e9Jee5wOedcGGnnXOj5nnPOOatd0ZjMSCmL9g84\nKvZ1ArAOOAs4kvSc9oGOs2DBApkvy5Yty/u1lYo65wpl7W+l/Ha9vOu/Pi3/sfWQ3Ll5tZTfrper\nnrk37dOHxTnnyEg750LPF3hdZrFXF1UjkFLuiX09CDwFLAYOCCEmAcS+HizmGhSKiuDgFnjm3+mY\nsITbo5egawLN47V+ZhilXZti2FM0QyCEqBFC1Nn/B84HNgB/BJbGnrYUeLpYa1AoKoJwDzy+FLzV\nbDrjNgx0PLpA1y1DIM1owW/x2vZWWtp7Cz6OYnhSTI+gEXhFCLEOWAk8K6X8C3AL8H4hxFbgfbHv\nFYqRy7NfgUNvwcfuoy9gNYr0aBqax8rlkEbhhuDzv13D/S/vKPg4iuFJ0bKGpJTbgVPSPN4KnFes\n91UoKoo1j8K6/4Gzb4RjzyWycT8AuibwiMHzCIIRg2BEhZgU6amINtQKxbDkwEZ49ssw4yzLEACG\naWUKeXSBpsUuz0EwBIYpiZrFyxBUVDaqxYRiQA51hVixvbXUyxhW9HV30P3Ip5D+evh/94OmAxAx\nTMAKDekeHzA4oaGoKR0jo1AkowyBYkAeenUnl/9qpZ3uqygQaZpsvPcqqrp2svvcu6Cu0fmZ4xFo\nAt2jPALF0KAMgWJAuoIRQlGTnrCKMQ8GK5+8nYWdf+f26CW0TTgt4WdRV2jI47E1gsJ+71Ja3kA0\n5m0oFMkoQ6AYkGDE2kCO9IZLvJLKZ9v6V5m34QcsF6dwt3FRyuYcNWyPQEO36wjMSEHvaRsX5REo\nMqEMgWJAglHrjvRIb2Eb0kinq6MN31NX0yHq6Lrg50g0Ikbi5myYlmHQNYGu26GhwjwCO9xUao2g\nJxRl5Y62kq5BkR5lCBQD0hcLCXX2KUOQL9I0efu+K5lk7ufwB37BqHGTgdTN2TYMXl2g6TqGFAVr\nBOXiEfx+zR4+ee9yuoLqc1RuKEOgGJBgNBYaUoYgb1Y+fisLuptZdcxnmXnaB/DoAoCImRgasg2D\nrlk/N9ALNgSGYXsEpdUIekJRTBkPNSrKB2UIFANiFyJ1KEOQF1vXvsz8TT9iXdVillz2PQC8mnXp\nGUmhIfuu3atbP4+iIwoMDdnGJjkMNdTYeki0xAZJkYoyBIoBCUWURpAvnUdaqX76atrEKI6++hE0\n3aoXsO/4kzdFe7OMewRa4R5BmWgEtiGKRJVoXW4oQ6AYkD7lEeSFNE3euW8pE8zDHPnQPYwaN9H5\nmdcODWXwCDy2IRA6YphoBHaxXHI4TFF6lCFQDIgd0+3oU+mjubDisR9was/LrD7+C5y06H0JP3Pu\n+M1kQ2CiawJrZlNMI5AFZg2ViUZgG6KIqmcoO5QhUAyI0ghy5+03mjl1y09YU30GS/7lWyk/tzWA\n5E0xakrHSACYaIPgEcRi8yXWCOxzLfU6FKkoQ6AYkKDSCHKio/UAdX+8hlYxhmOueRihpV5mdtZQ\nskdgGBKvyxAY6IgCPYJyCQ3ZBkB5BOWHMgSKAXHSR5UhGBBpmmx/4ArGyja6PnI/DWPGp32efdcf\nSQkNJXkEovDQUNQoF7G4PLKXFKkoQ6DoF8OUhKO2RqAMwUCs+J/vMb/3Vd448T844dSmjM+z00dT\nWkyYJh49flkag5A+ajgeQWnvxG0DoHoelR/KECj6JRSNb0LKEPTPllV/Z8HWO3mj5r0s+eTX+32u\nniE0FDWkkzEElkegycHRCJJrFoYaex3JXpCi9KjBNIp+sTOGxtb4aO0JEzFMR+gc6fT1dPH2ij8j\nTQNpGkx69dsc1MZx3LUPpdUF3NgeQbr00QRDMAgagVFuGkFUeQTlhjIEin6xheLG+gCtPWE6+yKM\nrfWXeFWlJxoJs/2nF3BKeL3zWFB6afnoHzhq1NgBXx8Xi1NbTCSEhkThhsAJyZSJRlDqEJUiFWUI\nFP0SNwR+Nu2zwkPKEMCqh27k9PB6XjvhK4yd2QTA6IlHc9zEaVm93r7rT/YIIoaZEhoqNH3U8QhK\nHJtXYnH5ogyBol/squKJDQFANZ4DePPF37Nk94OsHH0Bp/3LN/M6hhACXRMpd8dGSh2BjlZw+qjp\nHLuUqIKy8kUZAkW/2BrBhDrLEHSMwBTSrs4j7Fz3Ikgwo31Me+kGdunTmHPtvQUd1zIEaTQCV2jI\nEouHh0agCsrKF2UIFP0SSvIIRlzmUDRExy/+iTl9W5yHemQA7RMPU1VTV9ChvZpI2RSjSaEhOQga\nQbRMms7Z5xpWHkHZoQyBol/s6WQT62OhoZE2rvJv32BK3xZ+GvgM559zLgDjp57E0ZOPLvjQHl1L\n02sotaDMaxb2O4+6xGIppdPHaKiJewTKEJQbyhAo+sUODY2vswTijr7ChMuKYuNTsPJenvJfxPqJ\nl/DFJYsG9fAeTaT2GjKk05kUYqEhBkcjADsrqVSGoDxCVIpUVEK4ol/sMZW1fg91fg9HRkoH0tZt\n8PTn4aiF/Mj4Z8cQDiYePTU0lCwWS80zaBoBlHYTtg2SCg2VH8oQKPrFDg0FvDr1Vd6RoRFEgvD4\nUtA9GJc8yIEesziGQNPSiMWJBXtS6OiDpBFAaXUCJ0SlxOKyQxkCRb/YoaEqr86oau/IyBr6y1dh\n/5tw8S9p8zRiSornESRPKEv2CAYhNFQuHkHEVBpBuaIMgaJf7IIyv1ejYSR4BG8+AasfhDO/CCd+\ngENdIQAmFMUjSJM+mtJrqPDQkHvjLeUmbI+oDCuPoOxQhkDRL6GIgRDg92iMqvYO64Kyvn2bkX/6\nAkw9Dc61CsUOdgWBInkEmpa++6i7T5HQ0AsWi8skNKQ8grJFGQJFv/RFDPweDSEEDVXeYTuTINjb\nxf77P0lXRIdLfgW6F8DxCMbXBgb9PdOJxVFTOp1JYXiJxSprqHwpuiEQQuhCiDVCiGdi388QQqwQ\nQrwjhHhMCOEr9hoU+ROMmFR5dQAaqnx09kWQcvhdyOvvu54Zxk6+zueg4Sjn8UPdliEYVzf4H9N0\noSHDTJxQJoU+fDwCQ2UNlStD4RF8Edjs+v6HwO1SyuOAduDqIViDIk+CEYNAzBCMqvYSNkyn/9Bw\nYdXTP2dx+7M8pH+MZ/pm0x2K10oc6gpR6/dQ7Rv8khuPrqWKxYZEd4WGpOZBo7CNs+w8AmUIyo6i\nGgIhxBTgQ8D9se8FcC7wROwpDwEXF3MNiv4JR01HEE5HMGo6hqChygqXDCfBeNfm1cx64zts9M1h\n9Ie+DUBLe6/z80NdoaLoAxDzCFJCQ4ktJhgEjyBSLmKxqXoNlSvF9gjuAP4TnFuascARKZ2RSy3A\nUeleqBga/vvZTVz+q5UZfx6MaQQAo2KGYLjoBL3dHcjHlxIUfiZc8SjTxjcA0NLW5zynqIZAz5A1\nlKQRFFpHUA4egWFK7IiiCg2VH0VrMSGEuBA4KKVcLYRoyuP11wHXATQ2NtLc3JzXOrq7u/N+baWS\nyzmv2Rrk3S4z4/P3HugjHIXm5mZ2tFob0ovLV3FgjD5Iqx0ccv07S9NEf+123mO08KcZ36Th7Z10\nhqydatmq9XgOWkZv18FeptVpRfkMdR4J0h2RCccOhsPs37eX5uZWAMLdveik//tke87bd8SrwVes\nXMWBhqH/27lTRvcdOKiu5ywZqvMtZq+hM4GPCCEuAAJAPfBTYJQQwhPzCqYAe9K9WEp5L3AvwMKF\nC2VTU1Nei2hubibf11YquZzz3VteJdzZmfH5d295lXpNo6npNMbv7eBHq15h+gmzaJo9cfAWPAjk\n+nde9fufsij8EsunXcNFV3wZACkl//nKX6gaexRNTTMB6Gn+KzOPmUJT06xBX/Oju1ZhdgRpanqv\n85hY9lemT5vqvP/yt5/E02OkPbdsz/nV3s2wYzsA805dwLypowZl/bnQHYrC838FYNTosTQ15de3\naaRdz0N1vkULDUkpvyalnCKlnA58Evg/KeWngGXAJbGnLQWeLtYaFAPTHTLoDRsZY8fBiEnAa31M\nbI2gs8I1gh0bVzBn3U1s8M9j8dIfOo8LIZgyuprdMY0gGDHoCkaLFhrS02gEyRPK0DyFZw0Z7qyh\n0oRl3HOK1fD68qMUdQQ3Av8hhHgHSzN4oARrUMToiWXIuDNl3CRmDVkplJXceK6n6wj6k1fSLWqY\neOWj6J5Ep3jq6Cpa2i2NwKkhKJpGkJo1lNx0Ds2DLiSmkb8xcL9HqcZERszyEKwV6RmSNtRSymag\nOfb/7cDioXhfxcDYhqArGHU2ejfBaNwQ1Ph0dE1UrFgsTZPN913NfGMvm89/hNkTp6Y8Z8roat54\n9wgAB4ttCJLqCKSUKRPK0KzfvWFE0fT8YvvlUEfg9krUqMryQ1UWj3C6XYYgHX3heGhICMGoCu43\ntPLJ21nY+XdWTr+e2Wd+OO1zpoyuoqMvQmcw4qoqLpYh0JLCNjL2uNsjsMJxRjT/37lhlD5rKNEQ\nqNBQuaEMwQgmapiEYrHbrmD6jSbkCg2BpRNUYr+hbetfZd6GH7A+sIAll9+c8XlTRlcDsKe9z6kq\nnlBfHEPgTeo+am/S7tCQiHkE0QIMQaJHUCKNIOE8lUdQbihDMILpCcfjzpk8AndoCKCh2ltxYnFX\nRxv+p66iQ9Qx5apH+g2xTBldBcDutl4OdYXQBIytGRqx2N6w3RPK0KzorRHNfzKckRCfL5FG4AoH\n2V1IFeWDMgQjmB6XQNwVSt3cDVMSMSQBT5JHUEEagTRNtt53BRPNAxz+wC8YM6H/+sWpYyyPoKW9\nj0NdIcbU+BPF20HEqycOprFDOO4WE0K3DUH+An2kDArKbAPk92gJ3oGiPFCGYASTYAjSeAR26wlb\nIwAqTiNY+fiPOLX7RVYd81lmnvaBAZ8/utpLtU+PGYJg0YRisD2C1JBJOo/ALMQjMOKZSKUyBLZH\nUO3TVYuJMkQZghFM9wCGoC8SH1NpY3kElZE+unXty8zfdCvrqhaz5LLvZfUaq5agipb23qK2lwCr\nxUS6u/W0GoFRmEYQiLUJKZVGYJ9blVdXWUNliDIEI5ieUFwj6EwjFtseQVWCRuCjKxTNOg3xzZYO\n5nznr+zvCBa42tzoaD9MzdNX0y4amHbVwzmlXlpFZVZoqFgZQwBeTUvbByi5oAwK9AhME7/XNigl\n8ghiSQlVPl1lDZUhyhCMYAbyCOx5xX5XaGhMtRcpydor2LSvg65glLcOdBW42uyRpsn2+5cy3jxM\n+wX3MHr8pJxeP3V0FS1tvRzqDhUtYwisO3+rGVtie2b3hDLN1ggK9Aj8jkdQqoKymEfgUx5BOaIM\nwQjG1gh0TdDdr0YQv5seG7tDbu3JzhDYz9vf0TfAMwePFb+7mfk9r7D6+C9w0uL35/z6KaOr6QpF\niRiyuB6Bnhi3dzwCt0ag2x5BAYbAiBuCUrV3sI1ctdejKovLEGUIRjC9YWvzn1DnT1tHEIqmMQQ1\nVvVxa3d2hqAt9rx9QxQaevuNZk596zbWVJ/Bkn/5Vl7HsFNIoXhVxRDPDrLDNfGCMpdHYIeGCvAI\nDFM6f0OjRJuwHQ6q8umq11AZogzBCKY7phE01gfSi8Vha9OwhUZwewShrN7D9ggOdA6BIehto+FP\n13JYG8sx1zyM0PL7eNtFZVBcQ2B7BHY6pR0ySRCLY7OTC+015GgEpfIIYueoxOLyRBmCEUxPKIom\nrM2uv/TRKl/cIxgT8wjacgwNFd0jkBL+8FnGmK08d+IPaBgzPu9DDZVHYIvCRopH4DYEtkdQiFhc\nBhqBEReLpSzt7GRFKsoQjGC6Q1FqfB7qA960oaFgmtDQ6GovQsDhLENDrbE2DcXOGpq6+w/w9p+5\nxbyMww2zCzrWqGovtX5rAy5qaEi34/a2R5CqEdihocESi0tXRxAPDVnfK6+gnFCGYATTE4pS4/dQ\nF/D0mzXkriz26Bqjqry0ZRkasj2H/cUMDb37Gsdsfxh58kd4IHK+s+nli11L4Pdo1PmL16DXa3sE\nZmaNwA4NyYI9gtKmj9rvW+1VhqAcGZI21IrypCccpcavUx/w0B2OYpoSzRWWSFdZDJZOkI1YLKWk\ntTvstK4OJjWwy4kj78KvP2R9TUMwMBFxwZ2w5lV8BRoCsFpN9IYNhChOewmIawH2JmnH0RM8Ar3w\nOoKIYbpCQ6UqKIuHhkANsC83lCEYwXSHDGr9HuoCVm1ATzhKXcDr/NwxBL7EzXtMjS+r9NHuUJSw\nYXLSxDq27O9if0eQ6eNqcl9oNAyPXwF9R+Cs/yRqSu5u3sai6aM549hxIDTW9R3NTG8tQMEeAcAN\n/3Ri1jpIvnjt0FDs7tjeHN0ageYZHI3Ao4tYt9PSbMDhaKIhUP2GygtlCEYw7tAQWEVlaQ2BJ9EQ\njKv18db+gQvEbK9h1uQGtuzvYl++huDv34Y9q+HSh2HmRXT3hrn9789z7aQZnHGONds32NxMyC6A\nGwRDcEJjXcHHGAj7zj85NKSnEYtlgRqBrgmngK0U2AYoHhpSHkE5oTSCEUzcEFibf7JOEIyYaCKp\nCRqWR5DN3bLtNcycXA/A/s48iso2/wle+zksvh5mXuSsC3BmKdiEDdsQ5Bl+GmLsO397U7Q9A69r\nQplup4+aBXoEmrAG4ZS6oMznSfheUR4oQzCC6QlHqfHp1DoeQeJdpx3TT46Tj63x094bGfBito3F\nLNsQdKQXmA90Bjn71mUse+tgwuN7tm+m+7Hr6Ro7F86/yXncboZnewA2odjjg6ERDAW2KNyfR6A5\nHkH+hsDyCLSUbqdDiW3s7HYlSiwuLyrjilEUhZ6QkRIacpM8lMZmbK1VS9A+wFwCO3V06phq6gKe\njG0mXnzrELtae/n3x9ay54j1nFCwl97fXIYpJX8+6Qfgiadx2iGrsJHeI6gUQ6AnFZSlazqneeys\noUJaTJgxj6B0GkHUNNE1EW91oUJDZUVlXDGKotAdilLr91AfMwTJHUj7wmZC51Ebe2LXQNXFdmho\nbI2PSQ2BjCmky7e30lDlJWpIPvc/bxCOmqy9/3Mcb7zDVyLXc9g7MeH5tiGwW2DYDKZGMBR4k1pM\nxLOG3E3nCg8NRWNisUcvnUYQMeLhKVBZQ+VGZVwxikEnYpiEo2b/GkHUSOg8auNUFw+QQtraHabG\npxPw6jTWB9IWlUkpeW17K+85fhw/umQua949wt0/u5Ulh5/ktcZ/5nm5iGA4ccN3NIJIZXsEHqfp\nXOasIT0WGqLQrKESawQRw8Sna845J3tzitJSGVeMYtCxO4+6Q0PuttQQG1yfRngdFwsNHR5AMG7r\nCTm9iSY1BNK2mXi3rZd9HUFOP2YsF8yZxL+fqnNN++285TmJBVf/lCqv7mgCNnbFc4pYHK1MsTju\nEaSpLB6EFhNujaBkWUOG5ZX4dNsjUIagnFCGYIRib/q1fp0qr46uiTRisZlSTAZuj2Dg0JD93In1\nAQ51h1JEwuXbWgE47ZixEAny+dbv4/X6GHPFb/D6/AS8uuMBOOuKeQjhJENgh4oqxyNIFIvTTSiz\n6wgYlKyh0moEHl1zzrlU61CkpzKuGMWgY08nq/F7EEKkbTORqRJ4VLUPTQw8k6C1O+y0rZ7YUIWU\ncKgr0Xi8tr2V8XV+jh1fA3/5KtqBNwlceh/jpxwHMIBHkPh43COojI91PH3UWreRZjCNp8AWE1JK\nDFcdQanuxMNRiVcTKjRUplTGFaMYdHrC8dAQQK0/1RD0RYy0YrGuCUZXD1xd3NoTcjKMJjUEgMQu\npFJKlm9v5bRjxiI2PAmrH4QzvwQn/JPznIBXSzUEGeoI7O8rxyMYeDBNoR6BfUyvLvDoJawjME28\nHs0VGlIeQTlRGVeMYtBxNIJYgU9dmg6k/fUGGlvrc9JD0yGlpK0nzJhYhlFjvWUI3HMJdrb2cqAz\nxPsndMCfvgjTTodzv5lwnIBXd+oD3OuCdKGhmCHQK+Nj7WTQJBsCt1jssQypzNMQxGsTNDyl1ghc\nHoHSCMqLyrhiFINOXCy2Nvq6gIfONJXF6bKGYODq4s6gNepxXD8ewWvbWwkQ4p82fdWqE/jYA85o\nRpu0oaEBPIJMay434mKxmfDVHRrSHY8gv8E0buOil1AjiBgmXl1zzk2FhsqLyrhiFIOOPZ3M7rtf\nn0YjCGUoKIOBO5DaRsIWi0dVe/F7tASPYPm2Vn5Y9Si+1i3w/+6FhqNSjlPl0+lLSh/ty1BH4GgE\neoVkDWUKDbk8Ao9dUJanIbCH3ugxsbh03UeTs4ZUaKicUIZghOJOH4VMoSEzbfooWEVi/WkE9rwC\nO31UCMFEVwqplJLR7zzJRfIFxHu/DMe9L+1x/J7UrKFQhhYT4YrzCFJnFgtBQitwPWYIMPOrLI64\nWlvrmihZRW/EMPFoWkrthKI8qIwrRjHoxNNHbUPgSakj6IsYVPkyh4Y6+iIZe8bYE8zsrCGwUkjt\nNhNbN7zOjdF7OTBmITR9LeM6q3y6ownYOJXFRob00UrRCJI2xYghnWpjG6egLF+PwDXsxqtrJR1V\n6dM1p6FeWHkEZUXRrhghREAIsVIIsU4IsVEI8d3Y4zOEECuEEO8IIR4TQvgGOpZi8OkJRRN6v9jp\no1LGO2EapszsEcTu9NszeAXJoSGAibE2E73dHfieupI+EcB76YMpuoCbqn6yhsJR01mv/b1XFwl3\n1OVMcosJI9aPx42m6xhSFJw1VGqNwC4o8yqxuCwp5q1TCDhXSnkKMA/4gBDiNOCHwO1SyuOAduDq\nIq5BkYGekNV51O4sWuv3YpjS2XTj08kyh4Ygcy2BnVGUbAgOdATZeN+1TDNa2HfeXYyZOK3fdVoF\nZenrCCBRMA5FzYrxBiDedC7qajrn0VONmIGetyEoF40gYsqEgjLVfbS8KNpVIy26Y996Y/8kcC7w\nROzxh4CLi7UGRWbs6WQ2yR1InXnFGeLtjiHIIBi39oSp9XsSDMmk+gAXsYxFHX9lxbRrmP3eiwZc\nZ7qsIbd47DYE4ahZMTUE4MoassXiWIplMlF0RN5ZQ4kaQalE2kjUxKvFPQLVfbS8KOpVI4TQhRBr\ngYPA88A24IiU0r69aQFSU0UURac3HKU6rSGwREn7LtzfTx0BZO5A2toddp5jcxy7uMnzIBv881i8\n9Jas1umPtZhwh4CCSZu/TShqVEyfIUjfa0jXUi9JA63g0JAey+HPJzT02Kp3+UXztrzeP74OK300\nORymKA+KOqpSSmkA84QQo4CngJOyfa0Q4jrgOoDGxkaam5vzWkN3d3fer61Usjnn3fuCmFHpPG/n\nIWujaX51JS2jdPZ0Wxvsjq1v0dydugl0h60LecXaTTQc2Zry820tfXiiOMfXo73Mf/3LBPUaWmZ/\nicOvvJLVuezbbXkcz/9fM77Y3eSBQ/G5Bi++8g/GVWl0d3eze08QI2JWzN/bNm7vbN9Bs76Hlj0h\nzKiRsv556PR0daQ8ns3f+d1Oy6C/tXkTbYejdHXn/vt59PUg7UGTk9md0+vcdHb10ib6eOmlFxHE\nzrl5T87HGWnX81Cd75DMLJZSHhFCLANOB0YJITwxr2AKkPbTIKW8F7gXYOHChbKpqSmv925ubibf\n11Yq2ZzznZv+QYNPp6npNABqd7Zx++rlHD9zLmefMJ43WzrglVc49ZQ5NM1sTHm9aUr05j8zZtI0\nmppOTPn5LWtfYsaEapqaFoKU8OQ1ENoPl/+RD8x4b9bnstO7g8ff3sSi085kdCwcdfuGV6C9A4D5\nCxdz7PhampubGT2ujlFGN01NZ2d9/FLj/ftzTJk6jaamk/jjgbVU97al/O1am3VqqwMsSXo8m7/z\nmy0d8OornDJnDnvYx55Q6vEH4pdvL6dLBgu6jrwr/4+jJo2hqWkevr//mclTp9LUdHLOxxlp1/NQ\nnW8xs4bGxzwBhBBVwPuBzcAy4JLY05YCTxdrDYrM9IQMp70E4Mwk6LY1gqgtFqf/iGhOv6H0oaG2\nnnjDOVY/CBuegHO+DjkYAev9rVCPWycIRkzs6ZkJoaFIZWkEQEImT9SUCfOKbQx0RN6hIev3o8c0\nAiOPkEwoaqbUcuS8Dpf+4dU1FRoqM4rpEUwCHhJC6FgG53+llM8IITYBvxNCfB9YAzxQxDUoMmBP\nJ7PJpBFkyhqCWFFZGrHY7jM0ttYH+9bDn78Kx54H7/lyzuus8ukJ6wHLKNT5rZYYCWKxYVZM51Eb\nr6YlTChLTh8FMNERstA6gvw1gnDUTKnizpWIYToZQx5dqKyhDGxZ9XfMv30bj7Suq0bDYO/03zJ5\neqrXPZgUzRBIKdcD89M8vh38buCyAAAgAElEQVRYXKz3VWRHTzjqVBVDataQnZmTrvuozdja9NXF\nnX1Roqak0R+Gx5dC9RirhUQaIXQg0nsEBg3VXssQuB6vSI9AFwkTytJlDZlCy9sjiLjSR/OtIwhF\nzZS+Tvmsw9Z4vLqmsobScHj/bsY9ezUmGvsCxwIQkRFGD0HLlCHRCBTlR08o0RDU+DwI4fIIov2n\nj4JVI7Bxb2fK41a4SHLe1u9D+y644lmoGZfXOm1DEEwyBBPq/eymL7GOwDBp8Hnzep9S4dHim2J/\ndQSFegR2w7d8CrlCUYNgxEBK6dSd5ErU5RF4SzgXoVwxolH2P3gZx8ke9n78WU6ZvQSwNIKJU48r\n+vtX1u2TYlAIRQ0ihqTWH7/T0DRBrS/egdRJH+0nHXNcrZ/DaVpRt/aE+bT+PFP2/hXO+yYcfXre\na61yDEF84whGTRqqrA0/USMwKi80pMeLvDKlj5pCR8gCNQKnoCwPjyBiYsrCpopFXEbO69FUaCiJ\nlQ9/jdmhtayf+w2OiRmBoUR5BCOQ3ljn0Wpf4p/fPaUslIVGMCu6ifOMO9h0800Jj/sNk2943qJr\n6rnUnfHFgtZqGwI7VGWakrDLECRrBBUXGnIVeRmmmT40VEBBmVsj0PPVCIx42+90YnY2RAzTqSHw\naIJICUdV/ugvW3jv8eM5/dixRTn+6mfvJ7Du12TtO0nJkvBGVo06n0Uf/UJR1jQQWf1VhRAfF0LU\nxf7/DSHE74UQpxZ3aYpikdxwzsbdgXSgymI69/HRt29ktncvEpHwT9d13qx7L95L8tMF3Njv77S+\niImWcUMQ3yDD0QoUi11TwyIZNAJD5B8aSigoK8AjAFJafWSLYUqkxDEiXl0jUqDmkC+GKfnFi9t4\nem3uNQzZ8PYbLzJn5X/SEDmccl1k/Cc0Vo86n1nX3o8o8HrJl2w9gm9KKR8XQrwHeB9wK/ALYOh9\nGEXBJI+ptHF7BH39eQRGFJ64Co8RZPS/LmP0+OJlNCRrBLaBqk8XGqpAQ2AJuLGZxaZMa3hNoaPl\nawiMePdRXbOMTi6xfimlY2zzFYztMJDHJRaXqvlde28YKWG/ay7GYNHRdojaP11DmxhD3edeZMrY\n1PqbciXbq8b+FH4IuFdK+SyguoZWKMnTyWzcraiDEQOPJtKHApZ9H959FS68HYpoBCA1fdT+mjY0\nFDUrqsUExMIktlgc69mfTCHpo+5eQ96Yt5GLVxA1JfbT8/UI3HOT7bWUSiOw0533dwyuIZCmyfb7\nL2e82cqRC++hoYKMAGTvEewRQtyDVRT2QyGEHyU0VyzJ08lszows58TDT7L+Fg9nhw3O9Ep45L7E\nF5sG7HgRTl0Kp3yi6GtNTh+1N6P6QGpoKBQ1Kk4jcM8IiJqZ0kfz9wiSNQLnfbK0l8kFe/mQPILT\nq5VOLLYLIPcNsiFY8dvvc1rvq7x24lc4beF5g3rsoSBbQ3Ap8AHgx7F2EZOAG4q3LEUxSZ5OBsDe\ntVy1/yYO66Npj47Br0F1lQ6hrtQDzP0EfPCHQ7LWQGxjDzpxauvrqOqYIYh9L6UlIldSG2ogNjUs\nHhpKV1AmB1EjsN8nW0JJTf3ywRabvU7WkMjbqBSK7RF09EXoCxuOx1kIW15/gQVv38GamjNZ8sn/\nKvh4pSBbQ3CPlPLT9jdSyn1CiB8BfyvOshTFJEUsDnbA40vRaicw4fqXmVBTnGyKfPDoGj49PpzG\n/lrr96CJ+CZjSDAlFacRWOmj8WFA6UJxptDxmpnHgvaHe0KZnZqaS3zevfnn22bC1insc/NoGt15\nZkEVSpurAHJ/Z5AZ42oKOt6Rw/sZ9cz1HNLGccy1D5dM7C2UbFc9y/1NrG3EgsFfjmIoSPAIpISn\n/w06WuCSB6GMjICN36s56aPutFafR3PuWO0b10oLDXlcLSYyeQSm0NHIVyxOrCNwP5YNyW2+81tD\nzBg5WUOiZFlDra66l0J1AtMw2PXA5YyR7fR85H4aRudXNFkO9HvVCCG+JoToAuYKITpj/7qw5guo\nZnEVSm/YriPQYcU9sPlPcN63YVp5JoFVeXVnE4o3w9Pxe3THMNj7SqV5BB5dOAPmI0b6ymKpefLP\nGnIJtfax8w0N5esR2OfnTcgaKpVG4PYI+vp55sCs/t13OaVvBWtO/grHzz+r0KWVlH6vGinlD6SU\ndcCtUsr62L86KeVYKWXmieOKsqY7FMWjCfwH1sDfvgEnfBDO+Hypl5WRKp/ueATu+ga/R3NCQ3aB\nkq8Cs4bsjdnIIBZLoaMXKBYneAS5GIJI4R5BJEks9pSw+2hbT5ijRlUBBQrGu5Zz6ta7eMlzBosv\nvXGQVlc6stIIpJRfE0KMBo4HAq7HXyrWwhTFoycUZaIviHj8SqibBBf/HPLsITMUBDx6StZQlVfH\n79WcjSpSoaEhPanXULoWE7KQ0FAajSA3jyD9WNCc1mAkpo96NeEY8KGmtTvMUaOr6ApGOJCvIeg5\nDE9cxSG9kQfG/AdnVagu4CYrQyCEuAb4ItYgmbXAacByrPnDigqjOxjhZu0X0LUPrvqL1R20jAn4\n9JSsoYBXx6fHNYJIhYaGEnsNmc5m6cYU+YeGCvUIkns55UPEyRqKVxaXyiNo7Qlx4sQ6JjYE8vMI\nTBN+fx30tvLjhh9j+OoHf5ElINusoS8Ci4DXpJTnCCFOAm4u3rIUg82aZ+/D3PIcAJ/obmOJXAsf\nuAWmLCzxygamypuaNRTwxDSCqK0R2KGhyjIE7jCJYWRIH9V09Dw9gnhYRjjHzkUsDiVVbueD45WU\nQ0FZT5gxNT4mNlRxIJ/q4ld+AttegAtvZ8vyGTRW2OctE9meRVBKGQQQQvillFuA4paUKgaPHS8x\nd9V/MrVzDRO6tzCZg6ye+AlY8plSrywrAl49pbLY79Ws0JDtEcT2yUrzCKwGbPHuo+nSR2WBHoEm\nrO6yXj0PjSAhfXRwNAJrHsHQG4KoYXKkN8LYGj8T6/25ewQ7XoZlN8PsS2DBlYQiJv5+2rRXEtl6\nBC2xsZN/AJ4XQrQDu4q3LMWg0XUAnriad5nEg7N+xfcutWYUTy3xsnKhymUIQhEDIawNP11oqOI8\nAtf4yEwTyhBa3h6BVa1s/U7y0wgK9wgiyRpBnl1QC6Wt18oYGlvrQwKHukMZazdS6D4IT14NY46B\nD98BQhCMGhXX0iQT2YrFH4399zuxIfQNwJ+LtirF4GAa8OTVyFAXnwl/hw+OKW8tIBNVXpdYHDUJ\neHSEEPi9Oh19VrdUe2OptAvTSh/tv8WE1Dx5GwJ3bUJeWUODERpK0gg8JfII7GKysTV+PJqGlHCo\nK8TkWBZRRkwDnrzGKry87PfgrwOsjKpK80Azka1Y/IhdWSylfNF+DPh0vy9UDAkrHv8xYs9q5/tG\nEYQjj0PXXtj5Mu3vv4O3/jSBK+sD/RylfPF7dfrC8VbIdodOv0dzxMxKFYs9mtVryG7VnK7pHCJ/\njcA9/lJ3WkzkpxHkHxpK1AjsUZWFTDzLB7u9xJgan1VDg5VCOpAh2PK/3+SkHS/CR+6CibOdxyux\n220msg0NqcriMuX1Z+5lycabaKWBCF5MKfFogPG29YQzv8SOKRcBy5nYUJmGoMobLxzrCxtOIzq/\nR3OJxdZzKy40FBNO3V1Ck5GaB13mK9SaTrM52yDkMi/Y/r1X+/T8m86ZSVlDrp5H6c63WNjFZONq\nfY5XNJBgvOHlp5m5+ee8XPM+3js/8b43FDXw9zO4qZLo1xAIIb4GfB2oEkLYw2kFEAbuLfLaFAPw\n7ttrOXnVN9nsm8lxNzTj9fn510dXs2bHAV77jwuc5+1bvxeASQ0DuMBlSpVPSwwNxS4+n8ddR2CH\nhirMEMQmlLnTPFMoIDTk1gjsFg/5aAT1Aa9T1Z0r7swl9zqsSuq8DpkXdnuJMTU+tJgn0p9gfHjv\nLia98Hm2ycn8svbfeK/Le7HmNAwfjyDbyuLbgH8DfhL7fj7w1BCsT5GBvp4ujN8tJSx8jFn6KF6f\nH4Apo6to7bPcbhu7p8rECg0NBTw6UVMSMUyCrrnEVvpohYvFsTbUTvhkkNNHjTShoXzqCOoCnrw9\ngkhS0zlbNI4McZuJtp4wmoBR1T5GVXvxe7SMHkE0EubAry+jSga5pe5rdBjehJ9boa3Ku/HIRMaz\nEEK8JxYCAqjHKiL7ZOz7LuDuIq9N0Q9v3n89Rxu72N10B41TjnUenzqmmohpZUTY7O8IUuXVqa+q\nzBHV7uE0wUi8dbClEVgbZFwjqCxX3U4fdc8NSEHz4CH/HP5ksThXjcCna1T59IKbzrkLytyPDxWt\nPWFGV/vQNYEQot+islUP3cis8Ho2zv8W/skzU/os2b+LSvu8ZaI/c2ZijaMEWCyl/DcgCCClbEdN\nKCsZq/5wN4vbn2XFlCuY2/SxhJ9NGW2Ff1ra4w219nUGmdgQGFJhbjBxD6cJRaysISChjqBiNYJY\n9op95+1Jl8qoedCExDRy34gN03Ti8PGCstzqCPweq69T3m2ok/QP++tQZw61docYWxvftibWB9jf\nkdp4bn3zkyzZ/SArR13Aoos/R8CTagTtz91wqSPIeBZSyleBn8e+jcS8AwkghBgPed6iKApi1+bV\nzFrzXTb65rDoih+l/HzK6Gog0RAc6AhWbFgIXHOLwyZ97qyhWB2BlLJyNYLYpmhn5KT3CKzzN4xo\nzsePuD2CPArKwlGraCrgzd8jcEJDWqJHMNSGoC1WVWwzsSGQMrv4QMs2pjZ/iV36NOZca8mgfm+q\nEXQMQYV93jIxkEawNvbfO7E0gQlCiP8GXkG1mBhyers7kI8vpU8EaLziN3i8qU6Z3Vlxd1uv89i+\njmDFZgyBlTUEVgtqK33U9gisr2HDJGJaffPSbqRljL1eW4hNLxZb8WkjGsn5+G6NIN86Ap9emEeQ\nOrw+d89kMGjtDjO21u98P7EhwIGOkKOnRSNh2h76NH4ZQvvEw1TVWPUCfo+ekjrrnosxHMjKnEkp\nfwP8J/ADYB9wsZTy8WIurJIJRgweeGXHoN7xSNNk433XMM1oYe+5dzFu8tFpn1fj91DnjXsEpik5\n0FnhhsBnfUz7wgbBaGL6KFh3rVETfLpWceEvOxRkb7LpqlxFzCOI5mEIErKGnMri3DQCvzexr1PO\na0jTdA5KEBrqCTPW5RFMqg8QNkyn0GzVr77MyZGNbFp4E0efOM95XsCbmjo73DyCrNXDWH+hLUVc\ny7Ch+a1D3PTMJqaPrea8kxsH5ZirnrqTxR1/Y/m0azn9rIv6fe64ao2WdssjaO0JEzVlZYeGPHGN\nIBhJTB8F66KMGLIiL0rHI4j05xFYl6kRzT00VLBGEMvScusxuZLcYsI2SLnUMxRKxDDp6LP6DNnY\nN0f7OoK0rHiK0/c9zIoxH2HJh69PeG3Aa829cFdpxw3BCPIIFLlhb8Kb9nYO8Mzs2L5hBXPXf58N\n/nksXnrLgM8fVyXYE/MInNTRCvYIAklZQ+7KYrAuyqhZeUNpIFUjSNeGWui2Ich9bnE0jUaQSx1B\n2DBjYrFekFhsZ+qAK310CD2C9thd/xi3WByrq3nrrU0c/fKX2aYfwynX/jLltfaNR7qW3JV485GO\n4XEWZYYdltk4CIagu7Md75NX0C1qmHjlo+iegZ24cVUaLUf6ME3piGGTKtgQOBpBJEkj8MQv0IhZ\nmRelLaDam2y6wTS2R2Dm4RGkazERyUUjiJj4PBoBr1aQWOzWbpz00SGsIzjcbfcZSswa8hLlmBc/\nj0ca+P/5EQJVqcPs7c+VWycYMVlDivxxPIJ9hRkCaZpsue9qJpv7OHD+Lxg3MbueoeOqBOGoyaHu\nkJMeV9GhodjG3x0yiBgynj7qeAQGEbMyQ0PxUEPmrCERMwRRIw+xOKGOIKYR5DSPwOqwac2Hzl8s\ndmsf8fTRoQsNxRvOxQ3B+Do/N3p/x3yxlbeW/IApx81O+9qAK1nBRoWGFANiewTvtvXSGcz94rVZ\n+eRtLOx6gZUzPsOsMy4Y+AUxxlWJ2Dp62dcRxKOJhGyJSsP2CI7E2gjboSFHI4jYoaHK+zinpI/2\nExrKyyMwTVfXz/yyhuw6grBhYubRPjpqyISQVynE4tYeq8DSXUegv/0c1+jPsbrxEhZccGXG19qf\nt2Ca+c2VePORjqKdhRBiqhBimRBikxBioxDii7HHxwghnhdCbI19HV2sNZQCKSUt7X1MH2vl82/Z\n15XXcd5Z9w/mbbiF9YGFLPn0f+f02vFV1p+1pb2P/Z1BJtT504uQFULcEFhGNV5ZbH0NVXBoyJMS\nGkrnEdh1BIV6BHloBK46AiCvWcNR00wolCtFZXFrd7wFNQDtO+EP/wqT57Pg2p9nfiHxZIWE+c0R\n5RFkSxT4spRyJlZ7in8TQswEvgq8IKU8Hngh9v2woaMvQncoyj/NmgjAxr0dOR+j80grgT9cRYeo\nY8pVD6PpuX3YxjoeQR/7K7yGAOJx2HbbI3BVFoOdPiqHhUeQdkiK4xHkmz6af68hdx2Be525EDGk\n03EU3F1Qh84jaOsJo2uChiovREPw+JVWeezHfw2e/r1lJzSU4BEojSArpJT7pJRvxP7fBWwGjgIu\nAh6KPe0h4OJiraEU7G6zwkLzp41iXK0v58whaZq8c/8VTDQP0vbBexgz4aic1+DXBeNqfbS097K/\nwmsIwLrTFyLuEfhTsoaMmEdQeXdnXscQZPYINFsszsMjiBqpGkHuLSZ053eeTwppxEjvEQylRtDa\nE2J0tQ9NE/C3b8LeN+Diu2H09AFfm14sHl6hoSHpQiaEmI7VsXQF0Cil3Bf70X4gbaK9EOI64DqA\nxsZGmpub83rv7u7uvF+bD6v2W3HcA9s3MTFgsHLrXpqb27N+fdebf+TD3S/x7OjLqOnzsT+PtXd3\nd1Ov66zftoc9R0yOqw4N6e+gGPg02L7nAADvvLWZ5iNb2d1lbUpvrHuTUMSgq6O94s5z4yHr8/L2\nth0ArFvzBp3bEw1ax569nApseHM9uw7FW4dk89nu7O6lTfTS3NzsVNBu27GD5uY9Wa2vpy/E4YP7\n2BmyfvcvvvIqE6pz2/z27g8SCZnOWvf3WH+39Rs2UtP2Vk7Hyvd6fmtnkAAmG//3v5m16R52T/kw\n2w7UwYGBj/VOu7Xpr1y9huC71pa5abvlna5Y/g/8RZypMFT7V9ENgRCiFngS+JKUslMk9vSWQoi0\ntwVSynuJzTxYuHChbGpqyuv9m5ubyfe1+bD1pe2wdjMXv+8s9nm38cAr2znjPWdlFbbYuuYljj78\nEOuql/DBz9+Zc0jIprm5mZOPrue1ba2EjDCLZh1P01nH5HWscqHm5ecR/gDQycL5p3D2CePZcbgH\n/tHM8SeejPHOOiY3jqepqbLmJfneOQyrVzBh8lGwYydLFi1i5uT6hOesNfZDC5x4wgmccOrZzuPZ\nfLb9q5YxeeIomprmA+D9+3NMmTqNpqaTslqf+X9/Yca0qZwybRS8uYZ5CxZxQmNdTuf4P+++The9\nNDWdBcTan7y8jONPPImmBVNyOla+1/Ndm19lYdV+Zm37JRy1gKlX/Iqpnuz6Zo7f2wErXuGEk2fT\nNNsK+a6Nvg1vb+V95zQVVX8bqv2rqH6NEMKLZQR+I6X8fezhA0KISbGfTwIOFnMNQ01Ley91AQ8N\n1V5mTq4nYkjeOdg94Os62g9T88draBOjOPrq3HWBZKaOrnYmMjVWeGgILMHYzhqq8qamj1Zq1pCu\nJYaG0mUNFRQaconF9vvlOpjG79XigmkeKaTRpElkpeg+2tXdzZeP3AxCi+kC2TdPtjWCUFL6qFcX\nFZ2E4aZoHoGwbv0fADZLKW9z/eiPwFLgltjXp4u1hlLQ0t7ndACdOcm6s9u0r5OZk+uJRsKsefZe\njO7DKa+r3b2ME83DbPvw45w0bmLB67DbUUNlF5PZ+L0a+45YxXEp6aOxyuJKjNfasfNQP91HhW41\nnTPz6D5qmInFXB5Ny1osjsbaKrg1gnymlEUMM2EWczxraOgMwVU99zGVrfDPv4NR03J6rd+Vpmxj\nDa6vPE0qE8UMDZ2JNdz+TSGE3cX061gG4H+FEFcDu4BLi7iGIWd3ey9Hj7WqE2eMq6HKq7NxbweX\nLJjCqge/wul7H0r7OkMKXj/5RpYsPG9Q1uE2BJVcTGZT5dWdcZXJTedCEavXUCV6BI5Y7BSUpZ6D\npufvEUQMmVCtrGsi6w3YThW1W0xAfh6BVVDmqiOIrSc8RGLx6mfv55P8jY3TlzLrxA/m/Pr0BWVG\nRd54ZKJohkBK+QrWfON0DM5uV2bYNQRnHjcOsC66kybVsWlvJ+uWPc7pex9i5egLmXnlz1Jeq+se\nltTkFnvtD9srAZhQX7nFZDZVrna/8crixDbUvgLDaaUgOTSkpwsNxdqKyLw8AjPJIxBZewT2pm+3\nmLDWmbtHEE0y0l6P3fyu+B7B7q3rOGnlf7HFN5MT/uXWvI4RcLU4sRlO84phiLKGRgrtvRF6w0bC\nJjxzUj2vrVnH0ftuZJs+g7nX3kOgurboa7E9grE1vmHhwrr7vtubklcXCGGFVSJmZeZ0e/XEDdab\nNjRkG4I8NuE08flsNQJ3GwV38V6uRExJtbvFhJ3GmkeVci4Ee7uJ/PZywsLLqMsfceZ650rASR9N\nrCPwD5NZBKBaTAwqdo+hqa6wzKyJVfyIO/DKKL5PPjQkRgCsjXNcrb/iawhsEgxBrLJYCIFP1+iL\nGBjSmkdQacQ9gsxtqHVbIzCHViNw58oHvFrCY7kQNcwEA2eHicJ5trXOlvX3Xc8x5k7ePft2Jk49\nLu/jeHQNjyaSKotVaEiRAbvHkNsj+OCBexmtbWX14p+w4PhThnQ9p04bldBbpZKx20pAPDQE1ibV\nHbI2yIr0CLTEQq3+NIJ8QkNW1lCiRpCtR2Bv1D6XRpBPK2qroCxuCIQQsRBV8QzB63/8BYvbn2H5\n5Ms5/ZyPF3y85AltKjSkyIg9HvIo2yPY8hyj196DsfBqFlxwzZCv555PL6i4iV2ZsN1zTST27Pd7\ndTqD1gZZiR5BNk3nNI/lEci8KotTNYJs0zbdU7jcqbq5r0GmtM7w6KJovYZ2bXmDmau/zSbfbBZd\n+ZNBOWbAq6dUFg+HkKtN5V05ZUxLex/1AY/Vz6R9F/zhMzDpFPQP/KAk6xkuRgDiHkHAqyecl9+j\n0RUzBJV4h+ZJFovTtZjIMzRkmhJTkncdgRMa8upp++1kS8Q0UwyBV9PyamA3EL3dHcj/XUpQ+Bl/\nZfq53vlgGYJkjaDyPm+ZGD5nUga0tPdaYaFoGJ64EqTMqqmVYmDsrKHkYeE+j0ZXrNV3Jd6heZLF\n4jRejR4LDZFjaMiQiSMi7ffLXiMYPI8guT7C69GK4hFsvP96phm72XPunYyfPH3QjutPGsxj1REM\nn+1z+JxJGdDS3sfUMVXw/Ldgz2q46GcwprJbO5QLdoZGVZIh8Ht0xyOoxDqCZLE4/cji/ArK7Dt/\nt0bgyckjiGsEmmYJ8/k1nZMJTefsdQy2RrDyqbtYdOTPrJx6FXPO+uigHjt5VGcoaqisIUUqdg3B\n+1gJK34BSz4DM/sfMq/IHqetRJI77vdodPZFnP9XGvbdeihqxfLThfOc8aQ5hobsO39PUmgo1zoC\n+/dqCab5VRYnz2L26hrh6OB5BDs2rWLO2u+x0XcKi6740aAd1yZ5VOdwE4uHz5mUmNaeMOOie/nI\nrv+GyafC+28q9ZKGFVUxAxDwpAsNVa5H4M6pTycUA+hafllDdsGWnlxQlrVYbKePxo1wPh5B1Eij\nEeiD5xH0dB1Be+IKekQ1jVdlN9c7VwKeZLF4eLWYqLwrp0zZc/gId3vvRMujqZViYAKORpDqEdit\nJyrxwkzO8U+HVqhHoOfnEYSjyR6Bnp9HkMbIeXQtbfaSaUpW7WxzWmYPhDRNNt93DVOMPex//88Y\nNzG3PkLZEvAmpY8OszqC4XMmJabvmRuZq+2g7f13wOijS72cYYc7a8iNe/OvRI9A04SjC2T0CGI3\nFTJHQ2A4oaHEhm+5VxbHDEEhHoGW7BFoaQfT/G3TAT7+y+W8tDW1MWM6Vv3+DhZ2Ps/K6dcz+8wP\n57y2bAl49dTQkMoaUrhZ/ez9nHb4KV5r/GcaF3+s1MsZlgQyicWui7ESDQHEN+p0nUehhBqBq8UE\nWOGRXJvOGbEU1mQj59XTh6jWtxwB4H9f3z3gsbe9+RqnvHkzb/pPZXGOc71zxV1QJqVUoSFFIrvf\neZMTV36DLZ6TWXD1T0u9nGFLIEP6qN8Ve65UV93eJDP1tvc4BWU5egSGnTWUqBEYWcbm43UEbo8g\nt9CQHf5JKSjTRFqPYGNstOvzGw848yfS0dXRhu/3V9Ih6ph81SNF0QXcuAvK3F1ZhwvD50xKQLCv\nh/BvP01U6Ixa+mjeTa0UA5Mxa2hYeASJM4WT0WOGIFePIGKmDrux2lDn2GJCj2cN5eoR2N5Huqyh\ndBrBpn2dzJpcT9gw+dO6vWmPKU2Tt++/iknmfg5/4BeMbcxtylk+uA1BMKIMgcLFuvs+w7HGDnad\n9ZOCmlopBiZTQZk/qe9QJWLn2GfUCOyCMjO3u/F4HUFi99FcQkNeXVgD34lthjl6BHb4J9nIedMU\nth3sCnKoK8THTp3CyZPqeXx1S9pjrnzixyzoWsaqYz7LzNM+kNN68sWtj7grrocLlXnllAGv/+ke\nlrT9keWTLueUcz9Z6uUMewIZ0kfdm3/lewTpDYGm6xhS5K4RGKlisUfLQSxOmsKVj0dgh3+SPQKP\nntrzaPO+LgBmTq7nkgVTWN/SwVv7uxKe07nvbeZv/CHrqhaz5LLv5bSWQgh4dEJR09IHlEegANj1\n1lpmvv5NNntnseiqwWlqpegfRyz2JX5k3Zt/pYp3A4WGAAz0vD2C1ME02WsE7s0uOXMmG+zNPrmy\nOF3W0Ma9HQCcPKmei8EQP0gAACAASURBVOdNxqMJnlgdF407j7Ry6lu30i4amHZV4XO9c8HvtOE2\nU7KphgOq+2iO9PV0YT52OSHhZ+wVg9fUStE/TvpoPx5BpV6Y9ibZ3yD0KDrCzK37qL3h60kagZGD\nRuBL+v3m2nQu6ngEaQrKkjyCTXs7mTK6ymraCJx30ni2vrGMdTWbEYCx8n5my1a2fegxTho/Kad1\nFErAacNtpBTaDQeUIciRN++7joXGu2w45wHmHjWj1MsZMdT6PdT4dBqTBu0k1BFUYBtqiGsDyeET\nNwYayPyEWk+SRhDJQSPwJ3lcOXsEpp01lBQa0lLFYlsotrmx5jmOMW6D5vhznhlzORcufn9OaxgM\n3N1XHY9gGNURKEOQA6v+cDeLjzzH8ilXcnqTqhcYSgJeneYbzmF0tTfhcfuOVRc4omalYW/U/XkE\nhtAReWoEhbShdhva5OraXNaQTix2h4Z6QlF2HO7holOOsh7Y+QozNtxB54wPsXfmVQD4axqoPdiT\n0/sPFu4JbcNRI1CGIEt2bn6dWWu+y0b/nKI0tVIMzPi61PRc+2Ks5JuzeEHZABqBzE8j8CbNC862\n11BqaMjyCKSUWc+6iGsEaQrKXFrFlv1dSGkJxXQfhCeuRow5hvpP3kO9v8553q6DzVm972DjntCm\nQkMjlJ6uI4jHr6BPBGhUukBZYbvnlWwI7LBJpvRRsAxBzh6BrRGkFJTlGxrSMKUVcuovjOXGNgTJ\nYTsrayi+jk37rEKymRNr4PefguARuOxJcBmBUmJ7BJZGMPw8guFzJkVCmiab7r+WqUYL+867i3GT\nVR+hcsK+K8uUelkJ6FmEhkx0RJ4eQUKLiRzrCNxx8HicPPt1pGt8B6kFZZv2djKq2svk9T+D7c3w\nwR/BxNlZv0+xcZ+7bQiSGyBWMsPnTIrEqqfuZFHH31hx9LXMfq+aL1Bu2HealVzb43HOoZ/QkNBy\n9ggiGVpMZG8IEjUCdwpl9mvIXFCWaAg6uGTMdkTzLTD3E3Dq5Vm/x1AQcKePVnC320yo0FA/bN+w\ngrnrv8+bgfksvrw0c4cV/WNvTpXspWcjFhfmEcR/OXqsoCybOH84aiaEdOwUylwMQTRDQZnXNbw+\napgc3r+bLwV+COOOhw/dBmU2b9vvSfUIVGhoBNDd2Y73ySvoErVMurL4Ta0U+WFfoN4KDg05LSb6\nMwRCR8j8NAJ3WMb+PWWjEySHhvyuOHm2ZCoo82hWiwkpJTsOdnKruJOA7IVLHwZ/bdbHHyocjSCh\noGz4eATKEKRBmiZb7ruayeY+Dp7/c8ZNnFrqJSkyMByyhuzNOXmzdGOiIwahstguLssmPJQ8oN0Z\nYJ9DCmmmFhP29xFDcuiZ73KGvonWs38AE07O+thDiT9dQVklf+iSGD5nMoisfPI2Fna9wKoZn2XW\nGReUejmKfrDTGyvZS9cH6DUEsTqCHEND0XRN53LyCJI1gthmmENRWdQpKEvVCADWNf+e01oeZOWo\nC2g866qsjzvU2GJxKBKvI6jUAsZ0DJ8zGSTeWfcP5m24hfWBRSz+tJo7XO7EPYLKDQ15swwNabka\ngjTFXLo9IzmLNhPpWkxAbh5BfA2poyobaePYl7/ELn0ac669N+tjloJAUq8hn0er2ALGdKjAt4vO\nI60E/nAVHaKOqVc/MqRNrRT5Yd+lDguPoJ/cfEssznVUZfo6AiCrxnPJdQTOXXEOHkHYPZhGSmh9\nB4ww43v3cKfvZwQI0/WJh6mqKY96gUwkpo8Or3nFoAyBgzRN3rn/CuaaB3nngqFvaqXIDyd9tIKv\ny4EmlEGeHkGaoTD2ew0UGjJMSdSUaec95NJmIqHp3F+/Dq/9HICPAGjw+qk/ZOGJ87I+XqnwxGZL\n272GhpNQDEU0BEKIXwEXAgellLNjj40BHgOmAzuBS6WU7cVaQy6seOwWTut+ideO/QKnLTm/1MtR\nZEm8srhy3XRvFi0mpNDRcvYIMmsEA4nFznSyNG2+c/EIbM+j+p1nLCMw71Nwwj8RMUx2hetZuOC8\nrI9VSoQQzpSyZBF9OFDMs/k1kDw+6KvAC1LK44EXYt+XnLffeJFTt/yYdVVLWPyp75R6OYoc8A8H\nsVgfWCwuxCNIpxEM5BHE++m4Q0P5ZQ0dLfZT99cvwVEL4MI7YOZFeOd8lOMqxAjYWPMYTIJRY1hl\nDEERDYGU8iWgLenhi4CHYv9/CLi4WO+fLR1th6j90zW0idEcffXQDrtQFM5wCA1llT4qdDTyGxOZ\nziNINy/YTbpWy/l4BDLcx8+9PwVNh4//GjyV26cr4NFcHsHw2ieG+vJplFLui/1/P9A4xO+fgDRN\ntt+/lPFmK0cuvJdR4yaWcjmKPBBC0FDlpcZbuaEh3QkNZT4HqXkK8AgS21BDFh5BJLVoKpBHi4mT\n1/+AWdoujIt+DqOmZf26csTv1WMFZUosHjSklFIIkfHTKIS4DrgOoLGxkebm5rzep7u7O+Nru9f9\nngt7/8EzY5ZS262zP8/3KDf6O+fhyFdO9eA3eiv2nPftDQOw+91dNDfvS/ucQDhKnRFJOMeB/s7b\ndoQRwEsvveg89tZ+S2d4beUq9tRl3sz2dlub/Ttvb6G56x0gblg2vf0Ozca7A55Xx5ZlXHTkjzxf\n82G8B2rgQOa1ZkspP9vRUB8t+4L0Rqzfw1CsY6jOd6gNwQEhxCQp5T4hxCTgYKYnSinvBe4FWLhw\noWxqasrrDZubm0n32i2vv8CxbY+ypuZMPvT5OxD9CHWVRqZzHs5U8jmvCm2BHds4/tgZNDUdn/Y5\nb7x+G55eEs5xoHNeEdyCd+eOhOdENh2Ata8z/9QFzD6qIeNrN+zpgFdeYf7c2TTNsjxlKSXi+ec4\naurRNDWd2O85vfv2WsYu+yVbfDNp+uIDeH2psyTyoZR/57Eb/0FtwIMZjNJQ5aWpaXHR33Oozneo\nd78/Aktj/18KPD3E7w/AkcP7GfXM9RzSxnLMNQ8NKyOgqDzs0JDez+fQFLmHhgxTpqSkZps1lK6x\nmhCCgEcfsNdQsLcb43dLCQsvo5c+OmhGoNQEvJrTfXS4hYaKdjZCiN8Cy4EThRAtQoirgVuA9wsh\ntgLvi30/pJiGwc4HljJGttH94ftpGDN+qJegUCTgzaLFhNR09BzF4ohhphzTNgwDTSmzBWFf0obn\nj22G/bH+vuuYYe7k3bNvp3HKsTmtuZwJeHVCEYNwdPiljxYtNCSl/OcMPyppztjK33yH0/pe47WT\nbuS0U88u5VIUCsCVPtpPZbHM1yNIGRqfWx1BcnZMwKP3mz666umfs7j9WZZPXsrp53w8p/WWO36P\nNmwLyoaXWRuAzSv+ysJtP2N17dks+URZlDAoFK6Csn4yn/LwCKKmTClSs1NUB64jSN9z3+/VMjad\n27V5NbPe+A6bfHNYdOWPc1prJRDw6gSjsRYTlZyvnIbhdTb90HZwD2P//Bn2axM4/poHlS6gKBs8\n+sB1BFLkbggMQ2YODWVpCJLHMfo9WlqPoLe7A/n4UoLCz/grHh2Wc71tb0hVFlcqpsmBh66kQXYR\nvPhX1I8aW+oVKRQO2UwoIw9DEO1HLDYGaDrntJhIKrC074qT2XjfdUwzWthz7s8YP3l6TuusFAIx\nb2g4hoYqtulcJBKhpaWFYDDY7/MaGhrYtG4l4j1fZoNeT7VvFJs3bx6iVZaGhoaGYX+OyZTjOQcC\nAaZMmYLX6+33edlMKJOaB11mX8gFVp+fZN1BdyqLs2wxkYVHsPKpO1nc8ReWT7uG088avnO9A16d\n3rBB2Bh+HkHFGoKWlhbq6uqYPn16v7NXe9sPUtW3hw45A8/YGdQG+r8ohwNdXV3U1ZV3W9/BptzO\nWUpJa2srLS0tzJgxo9/n6lm0mEDzDI5HkGX30XhlceKaAl6dnlC8+d2OTauYs/YmNgTmsXjpD3Na\nX6Xh92hxEV1pBOVBMBhk7Nix/Q/gNiIEgvuJCi/7xQRq/BVr9xQVhhCCsWPHDuixQrxN9GCnjxqG\ndIRoG1s8zr6OIDEE4vdoGKEe+nq6aD+0D+2JK+gR1UwcAXO97dkXMLzmFUMFewRA/0ZASjiyCyFN\ndsqJ1Nf4+3++QjHIZPt582SVNeTBQ66hoUHQCNwegZRc3/pDFnU+D7dCFVAvBZvPf4TZEyu7j1A2\nBLypsxmGC8PrbNwIAbWNHPGMp0/6aKgq35BQbW3tgM+544476O3tdb6/4IILOHLkSEnWki9NTU28\n/vrrRTt+ttx8880J359xxhklWomFRxu4jgDNgyYkppG9V2D0oxEMNKoyFDXwaCLRkLz2CxZ1Ps/m\nSR/ltWO+wGvHfIHN73uI2Wd+OOs1VTKBhE6sw2vrrGiPYED8dbSaJj5dUO2rbFfujjvu4LLLLqO6\nuhqA5557rsQrKj1SSqSUaDmmAt988818/etfd75/9dVXB3tpOWFrA/21mECzPr/RaCQlkycT/WkE\n2YSGEja73avg+W/CSRdy8icetG60RhgJ09q8lb2fJDO8zFoSUcOkLyppqPYWJSz06KOPsnjxYubN\nm8f111+PYRj88pe/5IYbbnCe8+tf/5rPfe5zANx2223Mnj2b2bNnc8cdd6Qcr7m5mQsvvND5/nOf\n+xy//vWvufPOO9m7dy/nnHMO55xzDgDTp0/n8OHDGY+7c+dOTj75ZK699lpmzZrF+eefT19fX8p7\n7tixg9NPP505c+bwjW98I+Fnt956K4sWLWLu3Ll8+9vfBqCnp4cPfehDnHLKKcyePZvHHnss5Zhr\n167ltNNOY+7cuXz0ox+lvT0+hO6RRx5h3rx5zJ49m5UrVwLw4osvMm/ePObNm8f8+fPp6urK+P47\nd+7kxBNP5PLLL2f27NncdNNNGX/fF198MQsWLGDWrFnce681HP2rX/0qfX19zJs3j0996lNA3AuS\nUnLDDTcwe/Zs5syZ45yb3fjrkksu4aSTTuJTn/oUUkrneDNnzmTu3Ll85StfSfldZIPtEfQ7ZU2z\nPFojGsn6uNGC6giMeFiotw2euBLqJ8NFPxuRRgCUR1D2fPdPG9m0tzPl8ahpFX9U+XS0HD+8MyfX\n8+0Pz8r4882bN/PYY4/xj3/8A6/Xy2c/+1l+85vf8LGPfYzTTz+dW2+9FYDHHnuM//qv/2L16tU8\n+OCDrFixAiklS5Ys4eyzz2b+/PkDruULX/gCt912G8uWLWPcuHEJP0t33IULFzJlyhS2bt3Kb3/7\nW+677z4uvfRSnnzySS677LKE13/xi1/kX//1X7n88su5++67ncf/9re/sXXrVlauXImUko985CO8\n9NJLHDp0iMmTJ/Pss88C0NHRkbLeyy+/nLvuuouzzz6bb33rW3z3u991DFRvby9r167lpZde4qqr\nrmLDhg38+Mc/5u677+bMM8+ku7ubQCCQ8f2n/f/2zj08qupa4L81k5AhoCFQDFBeQZEEyINoEmIa\nQOmFiN5iQN4I3MtD7C1gaSl4y+3XfmgLlooiFnyA8pZKAR9QBSPIQ5D3m5RnkGARhCZCEBgm+/5x\nzoyZZPKCJEPO7N/3zZc55+xzzlpnT2bNWnvttZs359ixY8yfP5+OHTty4cIFUlJS+N3vfuf1vAHm\nzZtH/fr1+f7770lMTKR3795MnTqVWbNmsXfv3mJyr1ixgr1797Jv3z6+/fZbEhMT6dSpEwB79uzh\n0KFDNGnShNTUVLZs2UJ0dDQrV64kKysLEbnlUJ07A6W0TBQp5BGUF5ePmcXuwWNXGbWGbrhz5QsK\nYNUzcPkcDP8EaoeX+/5Ww+Fj/WarYC1tinDTpRChwkagPGRmZrJr1y4SExOJj48nMzOTkydP0rBh\nQ1q1asW2bdu4ePEiWVlZpKamsnnzZjIyMqhTpw5169alV69ebNq06bbl8HVdd6gjMjKS+HhjYfAH\nHniA7OzsYudv2bKFAQOMslBPPfWUZ//atWtZu3YtHTp0ICEhgaysLI4dO0ZMTAzr1q1j4sSJbNq0\nibAw71LGeXl55Obm0rmzUcdp6NChbNy40XPcfa9OnTrx3XffkZubS2pqKuPHj2fmzJnk5uYSFBRU\n4v0BWrRoQceOHQE8z3v79u1ezxtg5syZxMXF0bFjR86cOeM5v7RnOWDAAOx2OxEREXTu3JkdO3YA\nkJSURNOmTbHZbMTHx5OdnU1YWBgOh4Phw4ezYsUKT9iuoiS1rM/0PnF0aFbKl6zN+M1WUBGPwNcY\nQUVCQ8E22PoqHP0Yur9gLDUZwBQeLHZYLDRkCY+gpF/u15wuvruSzz3hd1f6PZVSDB06lD/96U/F\njvXv35+//e1vREVFkZGRUf7skaAgCgplc5Qn9bA0QkJ+KP9rt9t9hobAd3aLUornnnuOp59+utix\n3bt3s2bNGiZPnkzXrl09v8bLQ9F7iQiTJk3iscceY82aNaSmpvLJJ5+UeP/s7Gzq1Knjta9///6s\nXLmSEydOeJ73hg0b+PTTT9m6dSuhoaF06dLltp5n0Wd58+ZNgoKC2L59O5mZmSxfvpxZs2bx2Wef\nVfjaQXYbTz7QtNQ2Yjf+VSvqEZScNVT2PIIOZMGnf4C2PSFpVLnva1WsHBqyljZFcATbqR1UNfHM\nrl27snz5cs6fN9bWuXTpEqdPnwYgIyOD999/n6VLl9K/f38A0tLSWLVqFVevXiU/P5+VK1eSlpbm\ndc0WLVpw+PBhrl+/Tm5uLpmZmZ5jd911lyd2Xhhf161IFkxqairvvvsuAIsXL/bs7969O/PmzePK\nlSsAnD17lvPnz/P1118TGhrK4MGDmTBhArt37/a6XlhYGOHh4R5vZ+HChR7vAPDE3Tdv3kxYWBhh\nYWGcOHGCmJgYJk6cSGJiIllZWSXe3xcZGRmsXr3a63nn5eURHh5OaGgoWVlZbNu2zdM+ODgYp7P4\nF2paWhrLli3D5XJx4cIFNm7cSFJSyYuPXLlyhby8PHr06MGMGTPYt29fiW1vG9MjUBXIGnLexhhB\n8PVLTP7+z8bykj97NWDHBQrjNVis5xFoANq2bcvzzz9Pt27dKCgoIDg4mNdee40WLVoQHh5OdHQ0\nhw8f9nyRJCQkMGzYMM/2iBEjio0PNGvWjL59+9K+fXsiIyO9jo8aNYr09HSaNGnC+vXrPft9XTcu\nLo6LFy+WS49XXnmFgQMHMm3aNHr2/KE8QLdu3Thy5AgpKSmAMaC6aNEijh8/zoQJE7DZbAQHBzN7\n9uxi15w/fz6jR4/m6tWrtGrVirfffttzzOFw0KFDB5xOJ/PmzQOMjKj169djs9lo164djz76KCEh\nIT7vb/eRMRMeHk6bNm04evSo5zmkp6czZ84coqOjadOmjSeU5H6WsbGxJCQkeBm/jIwMtm7dSlxc\nHCLCiy++SKNGjcjKyvL57C5fvkzPnj25du0aSileeumlcj3zW0FMQ3DTdbsegTmhrEj66E3nDQ7M\neIK2+cYA/ksUoMQGfVeCo+SVzAIJL4/AYjOLxZ39cCfz4IMPqqL550eOHCE6OrrMc++00gPVgdb5\nzqG8n9Oy2LFqFol7f8vXQ7+kSWQUUPYyhl3/soGoxnfz2sAEr/2Rz61mzMP3Mb7bD8tNbn1zHCln\n32F7+GO4HEZRxrvbd6dd6uPcSfhzqcozl66S9qLxI2zzxIdpGn5rY0IV4Xb1FZFdSqkHy2qnPQKN\npgbg9ghcFfQIfM1WDrKJV2ho//rlHiOQNG7JbctqVUK8xgisFRqyln+j0VgU92BxxbKGioeGwBgn\ncA8Wf5NzgmafP8spW0tiR75ROcJaFO8JZdb66rSWNhqNRXHPIyioFI/AhtOlcN64zqX5g6mlnNj7\nz8cRWnXlRayAzhrSaDR+RezumcU3y2j5A06X8lm2wvAICtj59niinYc5kjiF5vfHV5qsVqWW3eZJ\nnqpVWsnwGogeI9BoagCe0JCr/IbAVVBQ8hhB1j9I+X4RXzZ4guTH9RyB8iAiOILsFChluUrG1jJr\nmoBmyZIlfPXVV/4Wo0qwuWcWVyA0dLNA+axo2tR2gQlXZ3Dcfi9xI/5aaTIGAo5gm+XCQqA9gjuC\nunXreiZOlcTLL7/MqFGjPGUMevTowZIlS6hXr161y1LVzJkzh9DQUIYMGVLucxYsWIDT6aR5c6Mu\nfnZ2Nl988QUDBw4s9bzytvM37tBQaR5B3qULnJ39BPc4cwDIVBC63w7/9P43X+rMowDBMWABjtp1\nfF1KUwIhQXZcNSDlvqJoQ1BDsFoZapfL5XNyGMDo0aMrfL0hQ4Z4zSPIzs5myZIl5TIE5Wnnb2z2\n0j0CVVDAybeG0P7GEfbUT0eJHQFaR9wFod5rcXzzbT5Xo56k3X3tq1psy+EItpU5K7smYj0fpxrR\nZaiLl6HOysryKsuQnZ1NTEyMR+aJEyeSkJDAe++9x5tvvkliYiJxcXH07t3bs/DO73//e6ZPnw4Y\nhePcZZ7d5SMK43K5mDBhAp07dyY2NpbXX38dMMpDb9q0ifj4eGbMmEF2djZpaWkkJCSQkJDgKcxX\ntN2hQ4c8fRobG1tmobrqwhbkLjHh2yP4cunzdLj6Bbvuf5akcUtIHruQpLELCe/3V/jPV7xekf/1\nFu1S0qtTfMvgCLbr0NAdyz8mwbkDPg/Vdt0E+y2o2SgGHp1a4mFdhtp3GeqoqChu3LjBqVOniIyM\nZNmyZfTr189zvEGDBp76RBcvXmTkyJEATJ48mblz5zJmzBiv602dOpVTp04REhLis8zz3LlzCQsL\n4/PPP6dWrVqkpqbSrVs3pk6dyvTp0/noo48Ao/z1unXrcDgcHDt2jAEDBrBz585i7caMGcO4ceMY\nNGgQN27cwFWB2j5VSWmDxVk7M3ng6MvsqZNK8oDJxY5rKo+Q4IqXtK8JWM+0VRO6DLXvMtQAffv2\n9XgKRQ1B4fcHDx4kLS2NmJgYFi9ezKFDh4pdKzY2lkGDBrFo0SKCfCyOvnbtWhYsWEBqairJyclc\nvHjR5694p9PJyJEjiYmJoU+fPhw+fNjn80xJSeGPf/wj06ZN4/Tp09SuXdtnu+rGbo4RKNdNvj13\nhsMvpNJ6w0i+/sP9NP9wABdsDWg1Yj5SwdXaNBWjdrDNcpPJwCoeQSm/3L+voho0ugx1yWWo+/Xr\nR58+fejVqxciQuvWrT3HCpeQHjZsGKtWrSIuLo533nmHDRs2FLvX6tWr2bhxIx9++CEvvPACBw4c\n8DIISileffVVHnroIa9+LnqtGTNmEBERwb59+ygoKMDhcPh8FgMHDiQ5OZnVq1fTo0cPXn/9dR55\n5BGfbasTzxiB8xrn3h7MvTeOsjMkGUft2uRIEPd0/xVN6jf0s5TW55ku9+EqKH1Rn5qI9UxbNaHL\nUPsuQw1w7733YrfbmTJlipcHUJTLly/TuHFjnE6n173dFBQUcObMGR5++GGmTZtGXl5esYym7t27\nM3v2bE9Z6aNHj5Kfn1/seeXl5dG4cWNsNhsLFy70hHyKtjt58iStWrVi7Nix9OzZk/3795f5DKsD\nW5DhEdy981XaX9/Lgfj/w/XQr0n85XskPbuUltFl1hXTVAKd72/II1ER/haj0rGGR+AHdBnqkstQ\ng+EVTJgwgVOnTpV47ylTppCcnEzDhg1JTk4uZuhcLheDBw8mLy8PpRRjx44tli47YsQIz0CwiNCw\nYUNWrVpFbGwsdruduLg4hg0bxs9//nN69+7NggULSE9P93gmRdtdv36dhQsXEhwcTKNGjbwWufcn\nNjM0dJ/rBDvC0knKGOvTg9JobgVdhtqCaJ3vHCqrDHXO8YM0XZRKtq0594zfTGjdML+WZPYXgaZz\ndZWh9ktoSETSReSfInJcRCb5QwaNpibRuGUU2yIGYB+4hNC6eqEYTeVS7aEhEbEDrwH/AeQAO0Tk\nA6WU7zQOjUaDPSiIjs/M8bcYGoviD48gCTiulDqplLoBvAv0LOMcjUaj0VQR/jAEPwbOFNrOMfdV\nmJowvqEJXPTnU1NTuGOzhkRkFDAKICIioliGRN26dcnJySEsLKzUPH2Xy+Uz7dLKaJ39j1KKvLw8\n8vPzqyy758qVKwGXORRoOleXvv4wBGeBZoW2m5r7vFBKvQG8AUbWUNGRc6fTSU5ODmfPFjvVi2vX\nrpU4eciqaJ3vDBwOB3FxcQQHB5fd+BYItAwaCDydq0tffxiCHUBrEYnEMAD9gQqXfgwODiYyMrLM\ndhs2bChXPR8roXXWaDQVodoNgVLqpoj8AvgEsAPzlFLFi8xoNBqNplrwyxiBUmoNULML6ms0Go1F\n0LWGNBqNJsCpESUmROQCcPoWT/8R8G0lilMT0DoHBlpn63O7+rZQSpVZlrZGGILbQUR2lqfWhpXQ\nOgcGWmfrU1366tCQRqPRBDjaEGg0Gk2AEwiG4A1/C+AHtM6BgdbZ+lSLvpYfI9BoNBpN6QSCR6DR\naDSaUrC0IbD6Ajgi0kxE1ovIYRE5JCLjzP31RWSdiBwz/4b7W9bKRkTsIrJHRD4ytyNF5Euzr5eJ\nSC1/y1iZiEg9EVkuIlkickREUqzezyLyS/NzfVBEloqIw2r9LCLzROS8iBwstM9nv4rBTFP3/SKS\nUFlyWNYQFFoA51GgLTBARNr6V6pK5ybwK6VUW6Aj8D+mjpOATKVUayDT3LYa44AjhbanATOUUvcB\n/waG+0WqquMV4GOlVBQQh6G7ZftZRH4MjAUeVEq1xyhH0x/r9fM7QHqRfSX166NAa/M1CvC9YPgt\nYFlDQAAsgKOU+pdSarf5/jLGl8OPMfScbzabDzzhHwmrBhFpCjwGvGVuC/AIsNxsYimdRSQM6ATM\nBVBK3VBK5WLxfsYogVNbRIKAUOBfWKyflVIbgUtFdpfUrz2BBcpgG1BPRBpXhhxWNgSVtgBOTUBE\nWgIdgC+BCKXUv8xD54AIP4lVVbwM/AYoMLcbALlKqZvmttX6OhK4ALxthsPeEpE6WLiflVJngenA\nVxgGIA/YhbX79Re7RAAABB9JREFU2U1J/Vpl32lWNgQBg4jUBf4OPKuU+q7wMWWkhVkmNUxEHgfO\nK6V2+VuWaiQISABmK6U6APkUCQNZsJ/DMX4BRwJNgDoUD6FYnurqVysbgnItgFPTEZFgDCOwWCm1\nwtz9jdtlNP+e95d8VUAq8DMRycYI9z2CET+vZ4YQwHp9nQPkKKW+NLeXYxgGK/fzT4FTSqkLSikn\nsAKj763cz25K6tcq+06zsiHwLIBjZhb0Bz7ws0yVihkbnwscUUq9VOjQB8BQ8/1Q4P3qlq2qUEo9\np5RqqpRqidGnnymlBgHrgSfNZlbT+RxwRkTamLu6AoexcD9jhIQ6ikio+Tl362zZfi5ESf36ATDE\nzB7qCOQVCiHdHkopy76AHsBR4ATwW3/LUwX6/QTDbdwP7DVfPTBi5pnAMeBToL6/Za0i/bsAH5nv\nWwHbgePAe0CIv+WrZF3jgZ1mX68Cwq3ez8AfgCzgILAQCLFaPwNLMcZAnBie3/CS+hUQjEzIE8AB\njIyqSpFDzyzWaDSaAMfKoSGNRqPRlANtCDQajSbA0YZAo9FoAhxtCDQajSbA0YZAowFEZKCINPe3\nHBqNP9CGQGN5ROQL829LERno4/hw4B6l1Fe3cO3/9XUvjaYmodNHNQGDiHQBfq2UerwC5wSpH2rb\n+Dp+RSlVtzLk02j8hfYINJZHRK6Yb6cCaSKy16x1bxeRP4vIDrO++9Nm+y4isklEPsCYzYqIrBKR\nXWZ9/FHmvqkY1TH3isjiwvcyZ3/+2aylf0BE+hW69oZCawssNmfOIiJTxVhbYr+ITK/OZ6QJbILK\nbqLRWIZJFPIIzC/0PKVUooiEAFtEZK3ZNgFor5Q6ZW7/t1LqkojUBnaIyN+VUpNE5BdKqXgf9+qF\nMRs4DviRec5G81gHoB3wNbAFSBWRI0AGEKWUUiJSr9K112hKQHsEmkCmG0btlr0Y5bsbYCz6AbC9\nkBEAGCsi+4BtGIW/WlM6PwGWKqVcSqlvgM+BxELXzlFKFWCUBWmJUWb5GjBXRHoBV29bO42mnGhD\noAlkBBijlIo3X5FKKbdHkO9pZIwt/BRIUUrFAXsAx23c93qh9y7APQ6RhFFZ9HHg49u4vkZTIbQh\n0AQSl4G7Cm1/AjxjlvJGRO43F3wpShjwb6XUVRGJwlgW1I3TfX4RNgH9zHGIhhgrjG0vSTBzTYkw\npdQa4JcYISWNplrQYwSaQGI/4DJDPO9grGPQEthtDthewPfShx8Do804/j8xwkNu3gD2i8huZZTD\ndrMSSAH2YVSI/Y1S6pxpSHxxF/C+iDgwPJXxt6aiRlNxdPqoRqPRBDg6NKTRaDQBjjYEGo1GE+Bo\nQ6DRaDQBjjYEGo1GE+BoQ6DRaDQBjjYEGo1GE+BoQ6DRaDQBjjYEGo1GE+D8Pw9RfRSRVFwtAAAA\nAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7Iabso2U0DQ0",
        "colab_type": "code",
        "outputId": "ca73f692-1644-408f-eb77-dac4121539f0",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 134
        }
      },
      "source": [
        "N = 50\n",
        "eps = 0.3\n",
        "col_1 = np.zeros(N)\n",
        "col_1[0] = eps\n",
        "col_1[1] = 1-eps\n",
        "row_1 = np.zeros(N)\n",
        "row_1[0] = eps\n",
        "A = toeplitz(col_1, row_1)\n",
        "A[0][A.shape[0]-1] = 0.7\n",
        "print(A)"
      ],
      "execution_count": 43,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[[0.3 0.  0.  ... 0.  0.  0.7]\n",
            " [0.7 0.3 0.  ... 0.  0.  0. ]\n",
            " [0.  0.7 0.3 ... 0.  0.  0. ]\n",
            " ...\n",
            " [0.  0.  0.  ... 0.3 0.  0. ]\n",
            " [0.  0.  0.  ... 0.7 0.3 0. ]\n",
            " [0.  0.  0.  ... 0.  0.7 0.3]]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hkB1pmwi0UJb",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def forwardpass(A,w,observations):\n",
        "  ###Initialisation\n",
        "  B = np.ones(A.shape)*((1-w)/N) + np.eye(A.shape[0])*w\n",
        "  pi = np.ones(N)/N\n",
        "  log_alpha_kk1 = np.log(pi)\n",
        "  log_alpha_kk = log_alpha_kk1 + np.log(B[observations[0]-1])\n",
        "  \n",
        "  for i in range(1,len(observations)):\n",
        "    v_max = np.max(log_alpha_kk)\n",
        "    log_alpha_kk1 = v_max + np.log(A.dot(np.exp(log_alpha_kk-v_max)))\n",
        "    p_yk_xk = B[observations[i]-1]\n",
        "    log_alpha_kk = np.log(p_yk_xk) + log_alpha_kk1\n",
        "  result = np.exp(log_alpha_kk)\n",
        "  return result\n",
        "\n",
        " \n",
        "def backwardpass(A,w,observations):\n",
        "  ###Initialisation\n",
        "  B = np.ones(A.shape)*((1-w)/N) + np.eye(A.shape[0])*w\n",
        "  log_beta_kk1 = np.zeros(N)\n",
        "  log_beta_kk = np.log(B[observations[len(observations)-1]-1]+ log_beta_kk1)\n",
        "  for i in range(1,len(observations)):\n",
        "    v_max = np.max(log_beta_kk)\n",
        "    log_beta_kk1 = v_max + np.log(A.T.dot(np.exp(log_beta_kk-v_max)))\n",
        "    log_beta_kk = log_beta_kk1 + np.log(B[observations[len(observations)-i-1]-1])\n",
        "  return np.exp(log_beta_kk1), np.exp(log_beta_kk)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YBrLg8qcFmtD",
        "colab_type": "text"
      },
      "source": [
        "###5.\n",
        "\n",
        "\n",
        "$\\mathbb{P}(x_{k}|y_{1:k})$ est égal à $\\alpha_{k|k}(x_k)$ à une constante près mais comme c'est une densité de probabilité on peut simplement renormaliser $\\alpha_{k|k}.$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vUbeq1r5tAiH",
        "colab_type": "code",
        "outputId": "d8623328-5c93-4d06-d773-a87e8c74fa92",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 251
        }
      },
      "source": [
        "alpha_kk = forwardpass(A, 0.3,[1,2,3,4,5,6])\n",
        "PxkY1k = alpha_kk/np.sum(alpha_kk)\n",
        "print('Voici P(xk|y1:k):')\n",
        "print(PxkY1k)"
      ],
      "execution_count": 65,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Voici P(xk|y1:k):\n",
            "[4.72276797e-08 1.96827423e-07 7.86771010e-06 3.92005741e-04\n",
            " 1.95999466e-02 9.79997961e-01 4.48901837e-08 4.48901837e-08\n",
            " 4.48901837e-08 4.48901837e-08 4.48901837e-08 4.48901837e-08\n",
            " 4.48901837e-08 4.48901837e-08 4.48901837e-08 4.48901837e-08\n",
            " 4.48901837e-08 4.48901837e-08 4.48901837e-08 4.48901837e-08\n",
            " 4.48901837e-08 4.48901837e-08 4.48901837e-08 4.48901837e-08\n",
            " 4.48901837e-08 4.48901837e-08 4.48901837e-08 4.48901837e-08\n",
            " 4.48901837e-08 4.48901837e-08 4.48901837e-08 4.48901837e-08\n",
            " 4.48901837e-08 4.48901837e-08 4.48901837e-08 4.48901837e-08\n",
            " 4.48901837e-08 4.48901837e-08 4.48901837e-08 4.48901837e-08\n",
            " 4.48901837e-08 4.48901837e-08 4.48901837e-08 4.48901837e-08\n",
            " 4.48901837e-08 4.48901837e-08 4.48901837e-08 4.48901837e-08\n",
            " 4.48901837e-08 4.48901837e-08]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3I92Eam-HIlz",
        "colab_type": "text"
      },
      "source": [
        "$\\mathbb{P}(x_{k}|y_{1:M})$ est proportionnel à $\\gamma_{k|k}(x_k) = \\alpha_{k|k}(x_k)*\\beta_{k|k+1}(x_k)$ donc on peut procéder de la manière suivante de façon analogue à juste avant. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eS2h4QnVHOss",
        "colab_type": "code",
        "outputId": "004294d0-34bc-4635-ce78-0058e4237b11",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 251
        }
      },
      "source": [
        "alpha_kk = forwardpass(A, 0.3,[1,2,3,4,5,6])\n",
        "beta_kk1 = backwardpass(A,0.3,[7,8,9,10,11,12])[0]\n",
        "PxkY1M = alpha_kk*beta_kk1/np.sum(alpha_kk*beta_kk1)\n",
        "print('Voici P(xk|y1:M):')\n",
        "print(PxkY1M)"
      ],
      "execution_count": 57,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Voici P(xk|y1:M):\n",
            "[4.43996873e-08 1.85041401e-07 7.39659180e-06 3.68532446e-04\n",
            " 1.84263023e-02 9.21315705e-01 4.10777703e-02 1.84263023e-02\n",
            " 3.68532446e-04 7.39659180e-06 1.85041401e-07 4.43996873e-08\n",
            " 4.22021605e-08 4.22021605e-08 4.22021605e-08 4.22021605e-08\n",
            " 4.22021605e-08 4.22021605e-08 4.22021605e-08 4.22021605e-08\n",
            " 4.22021605e-08 4.22021605e-08 4.22021605e-08 4.22021605e-08\n",
            " 4.22021605e-08 4.22021605e-08 4.22021605e-08 4.22021605e-08\n",
            " 4.22021605e-08 4.22021605e-08 4.22021605e-08 4.22021605e-08\n",
            " 4.22021605e-08 4.22021605e-08 4.22021605e-08 4.22021605e-08\n",
            " 4.22021605e-08 4.22021605e-08 4.22021605e-08 4.22021605e-08\n",
            " 4.22021605e-08 4.22021605e-08 4.22021605e-08 4.22021605e-08\n",
            " 4.22021605e-08 4.22021605e-08 4.22021605e-08 4.22021605e-08\n",
            " 4.22021605e-08 4.22021605e-08]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HJ8Kgs69JAO8",
        "colab_type": "text"
      },
      "source": [
        " $argmax_{x_k}\\mathbb{P}(x_k|y_{1:k})$ se calcule en prenant l'argmax du résultat obtenu deux cellules plus haut.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HVTnFxxBJ9RK",
        "colab_type": "code",
        "outputId": "a4d9191e-69ec-4f83-b198-eb259b8c5c6a",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 33
        }
      },
      "source": [
        "print('argmax_xk P(xk|y1:k)=' + str(np.argmax(PxkY1k)+ 1))"
      ],
      "execution_count": 58,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "argmax_xk P(xk|y1:k)=6\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mbbkAvxFcbl8",
        "colab_type": "text"
      },
      "source": [
        "### 6."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-rzJO2bHQLsc",
        "colab_type": "text"
      },
      "source": [
        "###1. \n",
        "\n",
        "Here is the new scenario: \n",
        "\n",
        "\n",
        "Let $N$ be the different states.\n",
        "\n",
        "Let $X_k$ be the true robot position after k iterations\n",
        "\n",
        "Let $Y_k$ be the state detected by the sensor.\n",
        "\n",
        "\n",
        "Let $E$ be a variable such that: $P(E=0) = 1- \\epsilon$\n",
        "and $P(E = 1) = \\epsilon$\n",
        "\n",
        "Let $U$ be a random variable following a uniform distribution on \n",
        "${1,...N}$\n",
        "\n",
        "Let $W$ be a random variable such that: $P(W=1) = w$\n",
        "and $P(W = 0) = 1 - w$\n",
        "\n",
        "Let $V$ be a random variable following a uniform distribution on \n",
        "${1,...N}$\n",
        "\n",
        "Let $K$ be a random variable such that: $P(K=1) = \\kappa$\n",
        "and $P(K = 0) = 1 - \\kappa$\n",
        "\n",
        "\n",
        "So we have according to the text:\n",
        "\n",
        "$X_{k+1} = (1-K)(R[(E X_k + (1-E)(X_k + 1)),N] + 1) + KV$\n",
        "\n",
        "$Y_k = W X_k + U(1-W)$\n",
        "\n",
        "$R(a,N)$ being the rest of the euclidiant division of $a$ by $N$ ie. we have to work with the modulo\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sbQvrdxfPSCB",
        "colab_type": "code",
        "outputId": "49dcb1ac-b813-4036-c8e1-8e04bfc63784",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 283
        }
      },
      "source": [
        "def simulation_2 (w,eps,N, n_iter,kappa):\n",
        "  X = np.random.choice(np.arange(1,N+1), p = np.ones(N)/N)\n",
        "  List_X = [X]\n",
        "  W = np.random.choice(np.array([1,0]), p = np.array([w, 1-w]))\n",
        "  U = np.random.choice(np.arange(1,N+1), p = np.ones(N)/N)\n",
        "  y_0 = w*X + U*(1-W)\n",
        "  List_y = [y_0]\n",
        "  \n",
        "  for i in range (n_iter):\n",
        "    \n",
        "    E = np.random.choice(np.array([1,0]), p = np.array([eps, 1-eps]))\n",
        "    K = np.random.choice(np.array([1,0]), p = np.array([kappa, 1-kappa]))\n",
        "    V = np.random.choice(np.arange(0,N), p = np.ones(N)/N)\n",
        "    \n",
        "    X = (1-K)*((E*X + (1-E)*(X+1)) % N ) + K*V\n",
        "    \n",
        "    W = np.random.choice(np.array([1,0]), p =np.array([w, 1-w]))\n",
        "    U = np.random.choice(np.arange(1,N+1), p = np.ones(N)/N)\n",
        "    \n",
        "    y = W*X + U*(1-W)\n",
        "    \n",
        "    List_y.append(y)\n",
        "    List_X.append(X)\n",
        "    \n",
        "  List_X = [u + 1 for u in List_X]\n",
        "  List_y = [u + 1 for u in List_y]\n",
        "\n",
        "  plt.plot(np.arange(len(List_y)), List_y, label = 'evolution des observations')\n",
        "  plt.plot(np.arange(len(List_X)), List_X, label = 'evolution des vrais états')\n",
        "  plt.ylabel('états')\n",
        "  plt.xlabel('itérations')\n",
        "  plt.legend()\n",
        "  plt.grid()\n",
        "  plt.show()\n",
        "  \n",
        "simulation_2(0.8,0.3,50,100,0.1)"
      ],
      "execution_count": 59,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXecHHX9/5+fmdlyPZdCKqTQQuql\nJ8bAIQoIKlJEMFQplp+C+jWCiiJiAUVQkC9NkC5NigJfisAZAgnpIQkJpCek55Kre7s75fP7Y3b2\ndu927+buttzl5vl45HHZ3dmZz+zezXte7yqklHh4eHh49F6UfC/Aw8PDwyO/eIbAw8PDo5fjGQIP\nDw+PXo5nCDw8PDx6OZ4h8PDw8OjleIbAw8PDo5fjGQIPDw+PXo5nCDw8PDx6OZ4h8PDw8OjlaPle\ngBv69+8vR4wY0an3NjY2UlRUlNkFdXO8c+4deOd8+NPV8122bNkBKeWA9rbrEYZgxIgRLF26tFPv\nraqqorKyMrML6uZ459w78M758Ker5yuE2OZmO8815OHh4dHL8QyBh4eHRy/HMwQeHh4evRzPEHh4\neHj0cjxD4OHh4dHL8QyBh4eHRy/HMwQeHh4evRzPEHh4eHQbIobJM0t34I3QzS2eIfDw8Og2zP/k\nAD957kPW7KzL91J6FVmtLBZCbAXqARMwpJRThRB9gaeBEcBW4Hwp5aFsrsPDw6NnUB/Wk3565IZc\nKIKTpZQVUsqpscfXA29JKY8F3oo99vDw8CAUNQFojP30yA35cA2dBTwS+/8jwFfzsAYPDwBW//d5\n1iz4V76X4REjFDWSfnrkBpHNoIwQYgtwCJDAfVLK+4UQNVLKPrHXBXDIedzivVcDVwMMHDhwylNP\nPdWpNTQ0NFBcXNzZU+iReOfsjrrdH3P6+p/xsXoM1SfemqWVZY/D8Xt+cWOUFzfqXDrGz8lH+Vq9\nfjiec1t09XxPPvnkZQnemLRku/voZ6WUO4UQRwBvCiHWJ74opZRCiJSWSEp5P3A/wNSpU2VnO/D1\ntm6F4J2zG2oO7KGp6kr8wsCv0CM/r8Pxe34/tA42bmbYiKOpPHFUq9cPx3Nui1ydb1ZdQ1LKnbGf\n+4AXgOnAXiHEYIDYz33ZXIOHR0ss02Tb3y6inzzELjEQBc8f3V1ojNguoUbPNZRTsmYIhBBFQogS\n5//AqcAa4F/ApbHNLgVeytYaPDxS8cHjv2RieAkrxl7HvoJRKNIzBN2FpliQOOQFi3NKNhXBQGCB\nEGIVsBh4RUr5GnAL8AUhxAbg87HHHh45o//211jnG8P0836MFKpnCLoRjhJoiHiKIJdkLUYgpdwM\nTEzxfDVwSraO6+HRHn4ZoT4wFKEotiHwXEPdBkcJhDxDkFO8ymKPXodPRjDVIACW0FA9RdBt8OoI\n8oNnCDx6HX4ZxdJsQyAVTxF0J+LBYk8R5BTPEHj0OgIygowpAukpgm5Fk+4pgnzgGQKPXkeAKDKm\nCFBUVE8RdBsaI16MIB94hsCjV2EaBn5hgq8AAClUVKw8r8rDoSnquYbyQbYriz26Cx8+A2tfAGBf\nfYSdR5zEpK9em+dF5Z5wUwNFgNAC9hOK5imCboJlSUK93DX0yd56/vLWBiK6xeDoVs7f8yf2Hv0o\nA488NqvH9QxBb2Hp32HPh9B3JMV7N1G6dzP0QkMQaWq0DYGjCBQvRtBdCBsmUoKmiF7ZdK62SefK\nR5ZyKBTluD7wy7obKbbqaZJq1o/tuYZ6C2YUjpoJ317AysBUNNk7+71HwyGg2RAgvBhBd8FJHe1f\nHEA3JRGj93wvliX5n2dWsqumiYcvm8Y/hz7FUXIX2yf+mKFHte65lGk8Q9BbMKOg+gEIWxq+XmsI\nGgFQ/AmKwIsRdAtCsUDxgJJA0uPewL3zN/Gfdfv4+ZknMGXvs7D2efjcDdSUT8jJ8T3XUG/B1EG1\n2/o2SQ2fjOZ5QflBjykCxV8IgFA0fMJEWhZC8e6L8klIt91BjiFojBqUF/nzuaSssfb9VzHn32a3\nN5FQEdZ5vdzPcRuKYfsiOO50mP1DmD8/J+vxDEFvIUERNFkafnqpIYjYhkBLUAQAlmWheoYgrzip\nowOKA0mPDzd2bf2YYW9cSYQA1b7BAAwsUhjRL4CwDBh9Jnz5z5DD30fPEPQWTB0UWxGELA2f7H3B\nOAAjZgjUQLMiADCMKKrm/TnkEydAnKgIDjci4RCNj8+lRFpEL/43JxwzLt9LArwYQe/BjDa7hiyN\nQC9VBGbUNgS+oG0IUOyMDNPonTGT7oQTLD6cYwQrH/gOxxob2Dj7NoZ1EyMAniLIOnVhnZ8+v5qb\nzxpH33z6OxNcQ42miqpJDD2K5js8fbDpMCNNAPgCsawh1VEEh9/dZ3dgxfZD/HP5p9x81jjsybTp\naakIDodW1MtefZDi5fcjsFCkyQxzE4sGzWXmqRfle2lJeIogy6zbVccrH+5m1ac1+V2IqYPqRzct\nmqStDCKxwGlvwow6hqDIfiLmLrM8RZAV3lq3j8cXbXd1UU9MH7Uf92xDsGHFfMZ/8BMKrHrCWhkh\nX18WHfF1pl7x53wvrRWeIsgyhmWPZDbMlKOZc0fMNRTWTaLYFz89EoaS/C4r11gxQ+AvsA2BiCsC\nzxBkg9omPf6zJNh6GH0iLdNHe3J1cc2BPRS/9E0OinJKv/s2R/UflO8ltYmnCLKMbto56oaZx1x1\nKeOuobBuEYkZgmjMTdKbkLp9zoFgcrDYMnv23Wd3pS5sG4K6pvY/Xyc43L/Ydlf21MZzzTOxD1L3\nlb/Rp5sbAfAUQdYxY4pAt/KoCCwTkDFDYBKV9tfupFL2JuKGIKYInBiBafTM4HntoQNsvf9CyiK7\nAZCoNFT+ivEnnp3nldkkKoL2aIqaBH0KRX4NIXpu47kPHv0Zs8JL+GDsz5kxuTLfy3GFpwiyjG46\nrqE8KgIzdpGLuYYi2HdcRiScvzXlC8M+58NBEViWyeYHLmZMaBkHC0ZysHAkw8wdNH70Zr6XFqcu\nZgAcZdAWjVGDQr+GoggKfWqPdA2tnv8CM7bex9LSzzP9vB/nezmu8RRBljEs2wDo3cIQ2K6haOxr\n16O9zzWEESYsfQRjxToirgh6Xoyg6cN/8rnQ+yw6fh4zv3EDAPt+NRIlkufEhATqwraBrXOhCEIR\nk0K/nc5bFNB6XLB4z46NDHv7GrarRzLm6od6VKV6z1lpD8UJEuv5DBabsT9C1UeTbsZjBEa09ykC\nxQgTEc0ps8LJGuphimDte69w+qEnWVZcyYwLfhZ/vlEpwRetzePKkumIaygUNSny24a5KKDR0IPq\nCKKRMLWPzMUvoygXPEZhcVm+l9QhPEWQZbpFsDhJEZhEZMw11AsVgTDCRAg0P44rgu5tCHZv+5j6\nxy6i1DgEwChZyw4xmOOvfjjpzrNJK8Wvdx9D0OwachcsLogpgkK/2u2DxQsfuJaRO18BQEPneGpY\nNuPPTDmuIs8r6zieIcgy8fTRfAaLWxgCxzVkRSP5W1OeUMww0QRFoKhOjKD7uoYi4RANj81lsP4p\n68srAdim+jkw5DTOLC1P3lYrpU94Zx5W2ZqwbhIx7BsgN66hpqhJUaDZNdSdW0wsfuEuZu18mNWB\nSYQL7KygzUfOZPoZl+d3YZ3EMwRZxlEC3dE1ZOq90DVkRtBFa0XQnQvKVj7wHWYYG1gx+26mJVSk\nVlVVtdpW95dRFFqXw9WlJzFA7MYQNEZN+hTaRrrIr3KgoXtmcm1avYgJK29iTbCCMfP+c1j0qPJi\nBFlGj8cIuodrKJJQR2Dqvc81pJphdKXZEDiKQHbTGMHSf93DjOoXWTRoLpNctCWwguUUy8YcrKx9\nEi/+7mIERlwRFHZTRVBXU03g+cuoEyUM/uYTh4URAE8RZB0na6i7xAgSFYHUe6ZrSFoWq//wBeT0\nq5j4uQs69F7NimAoKWIEOTYEH+2q46cvrOaWc8ZzwuDS+PMbVy0g+OIVFMsGACpkiI8C4123JZAF\nfSgUEcJNjQSdWok8UZtQROYmfTQUNSmMBYuL/Vq3azonLYtN91/MOGsfm854itEDh+V7SRnDUwRZ\nJq4I8hojaHYN2QVlsUyZHuoaqj24jwnhpTRtXtzh97Y0BEqsI6vMcYzgwQVbWLWjhu88vix+kayt\n3kvhi5fjl1E+7n8aH/c/jSVDLmLgFU+5bg6oFPYFoOHQgayt3S3OeQ0sDbiqLA5FjHj6aGFA7XYF\nZR88eROTQu+x7LhrGT3j1HwvJ6N4iiDLGN2qoCy5xYQ0eqYiOLRvB30AaXX84u2zIoSU/vHHzcHi\n3F106sM6r67ezZTh5azcUcO8Z1fxv9+YxNa/XcwJVjVbz3q+0xWpWlE/ABpq99N/yPAMrrrjOK6h\nI8sL+fRQ225IKSUh3aTIqSPw264hKWW7XUtzwUcL/4+pG+5kecmJzLjwF/leTsbxFEGWaS4o6x5Z\nQ4muIXqoImg4sAsA0Ym7eJ+MYqrB+OO4IuiEUeksL3+4mybd5IYzT+CnXxzN62v38tTtP2Ri0wes\nGPMTjutCWwJ/sa0IQrXdQBE4hqBvYbuuobBuISUUJNQRWJJ41lE+ObBnO0e8/m12K4M49qpHelSh\nmFsOvzPqZmQrWHzlI0t5esl2dxu3cA35/PaFsKcqgvChWHqk1fG7eL+MYKmtg8VWDusInlm6g2OP\nKKbiyD5c8dmRXDtyFxc0PMrSks8z/Ws/6dK+C8psRRCpy78hcALEw8oLCEXNNv8GnMBwc/qo/TPf\nMwkMPcreh75BkQxhnPcwJWV987qebJF1QyCEUIUQK4QQL8cejxRCfCCE2CiEeFoIcVhPRjHiBWWZ\nVQQfbK5m5Q6XhUOJWUOGSTDgR5cqmD3TEJi1doM10Ym7eD9RpNasCJysD9kJo9IZNuytZ8X2Gs6f\neiRCCET9bq6tvYX6klGM+daDXb7bLCwbAIDeeDATy+0SdWGDoE+JzxdoK4W0KdZXyAkWOz/zHTBe\n8tAPGRtdzZrJNzFy7Iy8riWb5EIRXAskJjbfCtwhpTwGOARckYM15A0j3n00s4ogYlpE3crmRNdQ\nrMOjjobooYqAxn0A9qDvDhKQUaSaKlicG0Pw9JIdaIrg7MlDbaX27GUoRpg+lz5FYXGfLu+/pNw2\nBFY3MAS1IZ3SoI+yAvszbiuF1FEETrC4OKYI8plCuuKNx5m1+3E+6PdVpp313bytIxdk1RAIIYYB\nZwJ/iz0WwOeA52KbPAJ8NZtryDfN6aOZUwRSSnTTIurW3ZTkGrIo8KlEhB/RQxWBL9Q5QyAtiyBR\npK8g/pyi5S5rKGpY/Gv5dr4wui/9CxR485ew4wP4yl0w4LiMHKO4pA+GVJBNhzKyv65QF9YpLfBR\nWqDFHqf/vkJxReC0mLDfk8vMIWlZ6NEIejTCjg2rOPr9eWzQjqXiqntytoZ8ke2soT8DP6F5DlY/\noEZK6Xy7nwJDs7yGvBLPGsqgIjAtiZSgd0YR6CZBn4qOr1sagrqwzhf//C53XljBlOGp/bHByH4A\nhOzYRSIaDRMQEqE1GwJVzV3TuTWLXmeBORf/ZhNujj0549sw7pyMHUMoCnWiGCWc/w6kdWGdsoJm\nRdCWa8hxARXGg8WOIsiNa6ih7hC77/wCxxobADgSqKGYooueiLcsP5zJmiEQQnwJ2CelXCaEqOzE\n+68GrgYYOHBgynJ6NzQ0NHT6vZlgx077Yrt77/6MrSNi2MZl977U+2x5zkN2ruU44L0PlrBnv+0W\niaIRaazL62eTil0NFjtrmnjm7WXUj0w92nBEtBqAaFPzebr5niNN9ZwGHKht3jZUs48hwJ5dn2b9\ns6hbt4DJwuTjQWfhKyjG0ErYHfg8MsO/2yMoxqzbm/fv9tN9TZQFBOtXrwRg4bJVWLtSX3KW7bUN\n8brVKwhtU9lRb9/kLF6+Cpnwnmz8PUvLIrDwD8zQN/J6ydmYMdehOWQGJVt388nW3Rk9XkfI1fUr\nm4pgNvAVIcQZQBAoBf4C9BFCaDFVMAxI2SFLSnk/cD/A1KlTZWVlZacWUVVVRWffmwle2LMCdu2i\nrLwvlZXTM7LP2pAO/3mDkrJyKitbB7BanfPCj2ADzJ5TSXDDWsoKfBifBijyCybn8bNJxbrddbDg\nXQr7D6GycmzKberfOQQCCnwqn4mt3833fGDXNvgABgw+khmxbQ/s2gYrYeCA/vHnssXiA8thL5Sd\neSODjjwagK44hNKd8/oFZZQoYcbn+buVi99h1NA+fP7EE/j5grcYNupYKmekrm2oWbETVqxkzqwZ\njBpQzI6DIXjvHUYcczyVU4+Mb5eNv+dF//gtM/WFLBz1PU679LcZ3XdXydX1K2sxAinlT6WUw6SU\nI4ALgLellHOBd4DzYptdCryUrTV0B5oLyjIXI4iYtlzuTLA4HDUJagqG8KOY3a+pl/M57alNXePQ\n1FhPibCLk5QOuoYiYXs0p+JPjBHE7oVykDXkxCE0f6CdLbtGxFdKgVGX1WO4obbJjhF0JFhcFHCy\nhmzXUCgLrqFoJELtwf3UHtzPmvf+zZT1f2Jl4SxmXPTrjB+rp5CPyuLrgKeEEL8BVgAP5mENOUM3\nMz+hzKlNiLjdp5NmqfoJGyYFfhVD+FGt7mcInOyq3WkMwcG9n8aDSh2NEegRuxlboiFQNTt7ORfp\nozJmeH2+7BqCqL+MAeGtWT1Ge1iWpD4WIwhoCn5VabPNRFOLYLFjEDKdNSTrdnHwjhMZJO04Uxmw\nSxnIyCsfQ1HVjB6rJ5ETQyClrAKqYv/fDGTGR9IDyMY8AidI7D5YrAMCFNVOH9VUTMXXPQ1B7JzS\nKYL6/Tvi/1dkx+4W9ZgiUP3NwT81h4oAw/683fYN6iymv4xiqz6rx2iPhqiBJaE06EMIQWmBr21F\n0CJYHNAUVEVkto7A1In+4xJKrDpeOuLbDCwvASEYMft8yvoOyNxxeiBer6Esk40JZU7aqPv00Sio\nfhCCsG4rAlMJEDS6zyQrB8dg7qsPY5gWmprsvQwdtENKNRSjdPDirUdaGwJNy10dgdPGwqnsztpx\nCsopFSEMPZp1o5MOJ0PISR0tLdDabDMRihrxiz+AEIJCv5rZyuL//IrA7iX8WP8eV35lHhOP7Hrd\nxuGCZwiyjOPzjmYwRuDEBtzHCHSIpUmGDYuAT8FSfGjdUBE4xs2SsL8hwuCygqTX9do9AFSrR6B2\n0DVkxgyBL5DoGoplJuVCEcRiBD5f6myoTCEK7Kll9TXVlA8YnNVjpcO5+3fiA6VBX9vpo9HmwfUO\nxZ0YYF9fe5DGOruYTlUE/Yv8dtO6rQtg4V9ZMehrvLp9Nn8cVNLOnnoXniHIMtmYR+CojA4Fi1Uf\npiWJGpbtGlIDaHS/qVyJQfVdNeFWhsCq34MhFRp9/SjQO1Y0ZcQMgZaQF66qjmsoB/nqZpSoVPFn\nuWmZVhRrRV2zL2+GwIkHlAZtQ1BW4KMmlP7GozFqxN1CDoV+tUN1BJ8s/y/DXvoag0Sa+pihU7lb\nvZxjBti1NB7NeIYgy+hm5mMEcUXQQddQxLD/qAr8KpYSwNcNFUFiUD1VnEBt3MdB0ceOcXRUEcQm\nsvkCzYZAUVVMKXKiCISpY6CRbWeNr8RuPNeYxw6kjhuo1FEEBT62Hwyl3T5xXrFDUUBzXVlcc2AP\npf+6gjpRyprx/w8QvLVuH4PKglw+ewQIFY7/Ih/esZzPHtO/vd31OjxDkGWa21BnPmuoQ8Hi2CwC\ngKCmYKl+fN1QESR+TrtrW/ewD4T3U6v2RSoaKh27i7eitmHxBZMnd5moOVEEwtIxRPb/5IIlTgfS\n6qwfKx0tXUNlBVo76aNmvAW1Q5HLKWWWabLtbxdxgjzE1rOeZ/rkkwB46YXV/GnlLi6e+AU0VWF/\nfYR99RHGDi3r7GkdtnhtqLNMNuoInIul6/TRmGuoSbf/qII+FakG8HdLQ9D8OaVSBMV6NSF/PyzR\ncUUgY4rA36JlgIHaqU6mHcYy0HNw71XYJ9aBtCF/hqA5WJwcI5Ay9d9BKGLEh9I4FAVUV+mjHzz6\nMyaGl7Bi7HUcFzMCALOO7kdDxGD1TjspYu0u++fYIaUp99Ob8QxBlmlOH82cIogkBIvT/WElEXMN\nhfVm1xCqH7/sfq4hJ5YS9CnsrmttCMrMg0SCA2xF0ElDEChoqQgU6GAqamcQZtRWH1mmJGYIjDx2\nIK1r0hECSmL1AGUFPgxLpi0QSxUsLvS37xpaPf8FZmy9j6Wln2f6eT9Oem3mKFsZvb/JNohrd9lF\ndmM8Q9AKzxBkGaOjgV0XJLpPXE0+i2UNOUU7AU0FLYgfA5nh9thdxTm3o/oWtlIEpmFQLmuxigZ2\nyjXkGIKWQ91NoXaqpXVHUSwdQ2Q3YwigpI/tA5eh/HUgrQsblAQ0lFg6qKMM0qWQhlIEi4sCWpvB\n4j07NjLs7WvYrh7JmKsfajXLoX9xgNGDSlgYNwS1HNW3MB7A9mjGMwRZJhvB4kRD4CpgnCJYjBZA\nERJd716qwPm8jupbyO6a5BjBoX07UYVElAxEKj40Onjx1sMYUsHXosWDiZobRSANjBy4hlRNo44i\nlHD+DIHTXsLBufimixOEUgWL/SqhNIogGglT+8hcfFJHueAxCotT+/1nHd2PpdsOEjFM1u6qY9xQ\nTw2kwjMEWSYb8wgS1YWrgHHcNdQcLEaLdSGNtD1UPNPsr2+79bVj5IaVF7K3PoKZYEBr9n8KgK/P\nEFA0tA4qAmGEiaTI2THJnSIwcxAsBmgQxaiR5lbUhxqjcddgLqhr0pPuvJtbUaf+nENRkwJfi/TR\nmCKwYr8DOzevpW73x2xYMZ8V91/N8cZ6Ppn5e446riLtOmaN6kdYt3j3kwNsqw4xdogXKE6FZwiy\njGMAMjmhrOOKINk1FPSpCJ9d3RoNp0/pyzRrdtYy7bf/YeO+9O0PHOV0VN9CTEtyoKHZcDRW21XF\nhf2G2IqggzECYTQRSTEZ1UJF5CBrKJeGoFEpwRe1g6OWJTnjzne59bX1OTk2NM8icIgPp0mhCKSU\nNEaNVorAmVLWpJssfuEuhj76Gb7y8U849qUvM6P6JRYd8XUmf/HyNtcxY1Q/FAF/W7AZ8OID6fDS\nR7OMc9GW0h4o45TQd4XEKmVXsQczCv4iwgmuoUgeFMHeWPB3b12EY45IXdnpnM+Rfe3Mnt21YQaW\n2kYrcmgXAKX9h3FIUTusCBQzQpTWDd9MoXS4gV1nENLISYwAIKyVEjBsg7tuTx27a8O8s34fN345\ndWvvTFPbpDOyf3Mspq0OpBHDQsqYyzIBJ2awbsUCJqy8ibWBiWwceAbDR4xACxYxfeYZ7a6jrMDH\nuKFlLNpsB869jKHUeIogyyTGBjJVS5B48Y90yjWkosQUgRFJ3dwtG7hpjWFYFqoiGNLHXt+ehFoC\ns95uL9F34DBQfB0OFitmmKiSRhHkIEagWjpWjhRB1F9GoWlnyTjB0q3VIXbV5Mbw1zUZyYogmD5Y\n7GQGFbUKFquU0sgRr11NnShh0BVPUnb0TCpOuYBxs7/sulvorKPt7KEjSgIcUZLdPk89Fc8QZBnD\nlPhUWwVkyhAkZw11wDXk1BH4FYRm/0Ho0dwpAseN1Zbx0k2Jpoh4a4nEdtRKw17qKLKzflQfqpBY\npvsLuGJG0EUqRZAj15A0sJTcGALDX0aJtBXBwk3V8dRMxyhkm9oWMYKSoBZ/viVJ84qlhJ3LYcu7\nDKtZxp989zBYHuDgGffTb+CwTq3lM0fbWVSeGkiP5xrKMrplURzQ0E0jYwHjxDtq164h1U8koaBM\njXXANKK5UwQRF60xdNPCryqUF9p97BNTSH1N+zmk9KUUEIp9kdH1CAHV3UxZzQxjpIsR5MA1pEoD\nXRS1v2EGsILllMoGdMNg8ZaDnFUxhNfW7OH9TdWcO6VzF1S3RA2LJt1MUgSaqlAc0FIGi5sNgQZV\nt8B/bwFgGoAKi477MTOnf6HT65k2opwiv8qU4eWd3sfhjmcIsogzZL7Ap1IfNjIWMO5csNjfHCzW\nVNRYCqWZwxiBK9eQKdFUgRCCwWVBdiUYgsJoNfU+W+YTaxZn6FHXw8U1K4KhtFYEllA7PNugM6hS\nx1RyEyMQheWoQvLBx9upjxjMPqY/tU06izZXI6W0O3JmiZZ9hhxKg6nbTDjVw8Oq34P/3grjzoMp\nlyGR7AoHmDG6a+NLCv0ab/7oJPoV56cld0/AMwRZxLlgO0GwjCkCszOKwEfYMFEVgU8VqD7b9WLk\n0jXkwhDopoUvNoNgUFkwKUZQalSzq2S8/SDWVtvQ3beG0KwIIX/rO3JTaLmJEUgjZzECpdDuQLpg\n9QbArrI91Bjl1dV72H4wxPB+2VMmLWcROJQW+FLGCJqiJkPZz5iFN8LAsfCVu8BfiID4NLquMqRP\nQfsb9WI8Q5BFnEBxgS/DhqCTrqGwbhHUFIQQaLGe/Kbedl5/JokP1DHSX3R1U8YNweCyApZstbM9\npGVRbh1ie4HdPkHEDYH79ftkBFNpHSy0FUH2XUOaNJA5UgT+Yls5rfpkC8cPHEP/4gCzYr7y9zdV\nxw3Bnh0b2fPJ0pT7GDJ6BkcMHdnhY9eF7c+yrKUiKEg9kyAUCnG3/05blZ3/KPjdKTyPzOEZgizS\n3DfHNgSu20a3Q6ddQ5HmPuxqLGvI0vOgCNqJETjB9UFlQfbWhbEsyeqqZ5goooi+I+wNYxdU03Cv\nCPwyiqWmdg3lQhFo0sDKkSEIlNqGQIRrmFVh///oAUUcURJg4aZqLpx+FDs3r6Pk0VOooDHlPlYv\nnswRP32nw8d23D8tWzmUBn18eqh13UrJf39JhbKJg6c9SN9+R3f4eB5dxzMEWaSVIshYjKATdQSq\nj7DebAh8AdsQ5FQRuEwf1eKKIIhuStatW82I+T9ikzqKiV/+HgBCcxSB+xYZfhnB0lorAinUnMxv\nVsmdIigste/++9AQT58UQjDr6H68t7GacFMjTU/MpQRYc8qjBIqTA6nmKz8hYDZ06th1LVpQO5QV\n+Fi3O1l5Lf33fcysfoFFg77BzKnndep4Hl3HMwRZxHEFZTxGYNiZNVHTat8QSJnUfTTosy+ymt92\nDVl6DusIXDTgixrNrqFBpUH+v5LhAAAgAElEQVQCRBHPXYZAEpz7OMHCYgAU1VEEHTAERLHUVK4h\nDZ/MvjLSMOKxjWxTVBYzBKKBmSP7xZ//zNH9eGnlLt7/6xV8ztzEyjn3UTHnrFbvX/lmH0oiezt1\n7NqmNMHiFjMJtq5bypilv+Aj/zimXvGXTh3LIzN4hiCLxIPFsbvwjBWUmRaFAZVoyGrfNeT00FF9\nhHUrbpT8QdsQyDwogrbmKBhWs2toSJ8Cfqk9xhi5iRWfuZtJo5qrYp0YgdUBRRCUUaTWOmgoc5Q1\n5MthjKC43DYEXyjaRNmml+LPn2JG+YH2Jp9r/D8WDrmUWadckPL9luJDle7dbrUH97Nx0b+Q0sLY\nWcuXlUOUb24EtblUaVrDHg7ou1n68nZAcsTyPxMSBRxx+ZNoPi+jJ594hiCLOK6hYNw1lBlFoJt2\nbUJNSG9fEZixC6WjCDTHNWQH5KSRhzoCl1lDo3b9m3HaWywafDEzT70oaTvFcQ25jBGYhoFfGOBL\nEyzuYJVyZ9AwIEcFZcGCIvbTl5Oi8+Gf8+PP9wd+oMHqwCSmXX5b2vdbih/NpSEIhxo4cPepTDHt\nfj5Tgcv8wIvJ250BnOEHYrHpqFTZcOpjjB0y3PV5eWQHzxBkESdY7FR1uh4t2Q5RwzYE4EJlJBiC\nJt2Ml/EHYjECmUvXkCtDYFcWs3ctha//GOPIzzDjsjtabSdU+w7ScukaioQbKYR4RXUiUmiouVAE\nGEg1d3e+RT9cghmpRm1RMxDSDUYPOBbNl16ddKSp34cPXM10czNLJ9/KgONnAlBe6GsVLJZIdtU0\n4Xz9JeUDGHtEphJEPbqCZwiyiN4iRqBnUBEUxQxBu72GzNhdXcw11K8o5hqKpY/GDUUOcBMj0E2L\nAb4IPP3/IFiGdv7DKf3qHY0RRJpCFAL4UriGlOwrAtMwUIXMWYwAoLCsP7YGaPG8i/daLuc9LH7h\nTqYfeoWFQy9j1le+3ea2Ahg6wMXBPXKOZwiyiJMlFHcNZSxGICmN9W7piGsokhAsVlSVqFQhh64h\np36grbiGYVhc1XA7NG2Fy16GkoEpt4vHCEx3d62RsJ0iqfjzEyPQ9Yg9pDKHhqBLKD587cy03rR6\nERNW/po1wQqmX/6nHC3MIxt4TeeySFwRxIPFmcsaCmgqPlW0Hyxu4Rpy1gIQxY/IpSJw4Rr6cugF\npoXehc//CoZ/Ju12TozArWtIdwxBKkWQA9eQHrWD8iKHrqGuIFU/vnZcQwfeuRsTlcHffAJV8+4p\nezKeIcgiRjxryP6YM1dHYOHXhJ1C2iHXUHMdAYAufIhcKoL2XEPbFnJF+O+sKJoDn/l+m/tS44bA\nXUAzGrbTQ9VAatdQR1tadxRTbzbIPQGp+tt1DalGiBqlrNNdQT26D54hyCJOlpAzYCOTbaj9qoJf\nUzocLHZcQwBRfDlVBBG9jcrihn3w7GXsFgN5ash10E5TNEWLBYtNd4bAiNgVrUqK9gVS0VDI3AS5\nlMePK4Ie4hpS/fiF2Wabb8WKpuzm6tHz8AxBFnEu0kF/5l1DvpghcBsjkEqsjiBJEfhRclBR65B2\nHoFpwHPfhHAtP9XmYfrb7xvfUUVgxLqsaqn62IjsKwI9lp0ltJ5x4XRcWHobdSaKpWPkqImeR3bx\nDEEWcSqJg5qS9Lir6KaFT1PwdcA1pMfyAgIJhsAQPhSzG7SYeOe3sPVd+NLtrJPD4wVlbaHGCpCk\nS0VjRG1FoKV0DWk5cA3Z30NPiREQM7RObCMV9gzmHqJwPNoka4ZACBEUQiwWQqwSQqwVQtwUe36k\nEOIDIcRGIcTTQhy+2tKJCcRbTGQoRuC0mPBrSptVukBcEUSlvYZgC0WQix47DikNwcf/Bwtuh8mX\nQsU3YpXF7f9aOopAulQEVswQ+IIp2i8LFS3LwWKnp5MT5O72xAyW0YYhUK2oZwgOE7KpCCLA56SU\nE4EK4HQhxEzgVuAOKeUxwCHgiiyuIa9kLWvItPBristgsX2hj0gtaS0ApuJHtXKnCByXUMRpQ31w\nC7zwLRg8Eb74ByA2mEZxYwjsLqJuYwTOAB5fiiE2UvGhZjtGEDNYPc011KYhkEbOBu14ZJesGQJp\n47Qv9MX+SeBzwHOx5x8BvpqtNeQbs2WLiYwFi+05yAFXwWL7AhS2HEXQ/JUbih/Vct9PpqvEs4ZM\nC/QwPHup/cL5j8ZbP0RNC5/mxjUUuwC5NAROu21/qmlmOckachRBDzEEPtvQ6m1UnnuK4PAhq5Ee\nIYQKLAOOAe4GNgE1UsYTlD8lzRAiIcTVwNUAAwcOpKqqqlNraGho6PR7u8rqT+2L1JqVywD4eOMm\nquSOLu3TkhLTkuzasZ2mRpNoiFbnl3jO/fevYBywfM16YCCbPllPVd1GAIoNKLaacvb5hML2xbC+\nMcyuhy5myO5VrB73c6pXbQW2AnYbjl07dlBV1Xbny3D9IU4H9u7eSVVVVbvfc/3O7QCsXLWGwMbk\n76CprgFNWLzz9tuIFmpECEFRURGqqtIVLLOAdac9g6n1YcWKFV3al0NpaWnG9tUSdeg01g16hqYD\ndew5mPoY6uf/gBRq1taQimyec3fE7fmapkljYyNSds7rkFVDIKU0gQohRB/gBWB0B957P3A/wNSp\nU2VlZWWn1lBVVUVn39tVdi/eDmtWc+LsWfDu2xx51AgqK4/r0j7Dugmvv8axx4xirzxA1LCorEwu\nvEo659UHYC0cfcJ4WLGPqRUTqBx9BAArlhQTaDqUs8/Heus1wOQrygKG7H4DPvsjxn/+J/HXTUsi\nX3uVY0aNpLLy2Db3VVu9F5bBgH7lzKysbPd7XrT9LTgIJ1Z+joKikqTXFm75P6iHOXM+26oL5pYt\nWygpKaFfv35dmvPb1FBDQd0WmkpHUlDcp9P7SaS+vp6SkpL2N+wEjbUHKGrcQaT8OAIFqcdaRnet\nQVeCFA06JitrSEU2z7k74uZ8pZRUV1dTX1/PyJEdnygHOcoaklLWAO8As4A+QsRzzoYBO3Oxhnzg\nuIJ8qoKqiIwEix33ihMsdltZ3BRzDQUSXEOmEnDdYTITRE2L0WI7P7PuhxFz4OSfJ73uuLm0DmQN\n4dK15XRZTTnoPtYRNFUn03A43GUjACTcqWVvaHwmEbF1StnW75dst97DI/sIIejXrx/hcOeLQ7OZ\nNTQgpgQQQhQAXwDWYRsEZxTRpcBLqffQ83GCwz5VoCkiI+mjTnDYfbDYvrhFYoYgMVhstxHITdaQ\naUkKrUbu8d1BLUVw3kOgJgtSPcHItYcWTx91acj0MBHpQ0nh4hGxdaQbe9lVIwDYA4IAIXpIxrZw\nDEH631mBpKcYtsOdrv6OZvO3cjDwjhDiQ2AJ8KaU8mXgOuBHQoiNQD/gwSyuIa84CkBT7Yt2JrKG\n9ASV0SFFYLZOH7XUQLuNxTJFVDe5zXcvR4r9/L/oNViFrdtQOobSjSLQnDRMy12QVxhNRNIFNuOK\nIHsD7OMX1DzeQRcXF7e7zZ///GdCoVDcYH31nPOoqalJua1AIjt5Pm7W0lkqKytZunRp1vbvlt/9\n7ndJjz/zmfS9s/JNNrOGPpRSTpJSTpBSjpNS/jr2/GYp5XQp5TFSyq9JKXOXv5hjnAu/pgg0VWSk\nxYRuOCrDbWWxfaEPWfZX3VIR+HOkCFh4F6epS7lLvZilcnRKA5Zo5NpDjd3FC5eKQDHDREmdsSNi\nhsBtlXKniLlYMqIusohjCByD9cJzz9CnT+qYhpCS3lCTKqXE6oRbt6UheP/99zO1pIzj6lsUQnxN\nCFES+/8NQojnhRCTs7u0no9hNl+0NVXJUIzAvgPuaB1BKIUikGoQv4ue811m63sE/3szr5rT+VfB\n2UDqfkPOvAY3lcVCUYhKFekyRqAYYaIikObFmGvIrZupU3TeNfT4448zffp0Kioq+Na3voVpmtx7\n773ccMMN8W0efvhhvve97wFw++23M27cOMaNG8ef//znVvurqqriS1/6Uvzx9773PR5++GHuvPNO\ndu3axcknn8zpZ3wZgBPGTeDAgQMp9yuQbNvxKSeccAJXXXUVY8eO5dRTT6WpqfX85y1btjBr1izG\njx+ftG6AP/7xj0ybNo0JEyZw4403AtDY2MiZZ57JxIkTGTduHE8//XSrfa5cuZKZM2cyYcIEzj77\nbA4dOhR/7bHHHqOiooJx48axePFiAP773/9SUVFBRUUFkyZNor6+Pu3xt27dyvHHH88ll1zCuHHj\nuPnmm5k3b17Kz/urX/0qU6ZMYezYsdx///0AXH/99TQ1NVFRUcHcuXOBZhUkpWTevHmMGzeO8ePH\nx8/NSXg477zzGD16NHPnzo0ryeuvv54xY8YwYcIEfvzjH7f6LLqK26yhX0gpnxVCfBb4PPBH4B5g\nRsZXdBhhWBZCgKoIfIrIiGsoGlMEflW4dA3ZF7dGw74AJdYRoPkJCB1pWa3SJjNG/V547nKMshH8\nZM/VjIxNrUplwJwJbm4UAYCB1gFFECGapoi9vRiBw03/XstHu+pcHa8llhFFsaJI3/IkYzBmSCk3\nfnls2vetW7eOp59+mvfeew+fz8d3v/tdnnjiCc4991xmzJjBX/5iD31/+umn+fnPf86yZcv4+9//\nzgcffICUkhkzZnDSSScxadKkdtd4zTXXcPvtt/POO+9QXBiEmg2O/Uq93zED8fcbzoYNG/jHP/7B\nAw88wPnnn88///lPLrooebTotddey3e+8x0uueQS7r777vjzb7zxBhs2bGDx4sVIKfnKV77C/Pnz\n2b9/P0OGDOGVV14BoLa2ttV6L7nkEu666y5OOukkfvnLX3LTTTfFDV8oFGLlypXMnz+fb37zm6xZ\ns4bbbruNu+++m9mzZ9PQ0EAwGEx7/KOOOooNGzbwyCOPMHPmTPbv38+sWbP44x//mPR5Azz00EP0\n7duXpqYmpk2bxrnnnsstt9zCX//6V1auXNlq3c8//zwrV65k1apVHDhwgGnTpnHiiScCsGLFCtau\nXcuQIUOYPXs2ixYtYsqUKbzwwgusX78eIURaV11XcPvX7zhizwTul1K+Aml0tkcc3ZT4YhdYTVUy\nUlAWD6h2SBEImmI3/omKgFh1biTS+g4uI5gGPHc5hOvYc/p9NFAYH7GZat2JMRU3GEIFy52iUa0w\nenuKIJsxgk6+76233mLZsmVMmzaNiooK3nrrLTZv3syAAQMYMWIEixYtorq6mvXr1zN79mwWLFjA\n2WefTVFREcXFxZxzzjm8++67HT5uSxdWy/2effZXefeDFSAEI0eOpKKiAoApU6awdevWVvt77733\nuPDCCwG4+OKL48+/8cYbvPHGG0yaNInJkyezfv16NmzYwPjx43nzzTe57rrrePfddykrK0vaX21t\nLTU1NZx00kkAXHrppcyf3zyb2TnWiSeeSF1dHTU1NcyePZsf/ehH3HnnndTU1KBpWtrjAwwfPpyZ\nM+3RmwMGDGDUqFGtPm+AO++8k4kTJzJz5kx27NgRf386FixYwIUXXoiqqgwcOJCTTjqJJUuWADB9\n+nSGDRuGoihUVFSwbds2ysrKCAaDXHHFFTz//PMUFrqZMdcx3CqCnUKI+7Azf24VQgToDc7BLmKY\nVjzw6VNFRkZVRlsEi121oVb9RAxbnQS05q/Nmd8bjYQJpskV7xJv/xq2vQdn30djn+OBfRS3MVkt\nUe24wURDuHQNaWYEQ0ltCBTVXYygrTv39mio3k1xZA/GgHFtzgpuiZSSSy+9lN///vetXjvvvPN4\n5plnGD16NGeffbbr+IOmaUk+71Rph+25sJqzYQWBQPPnqqpqSteQvWnr9Ukp+elPf8q3vvWtVq8t\nX76cV199lRtuuIFTTjmFX/7yl22uqa1jCSG4/vrrOfPMM3n11VeZPXs2r7/+etrjb926laKi5L+J\nCy64oNXnXVVVxX/+8x8WLlxIYWEhlZWVXUrjbPlZmqaJpmksXryYt956i+eee46//vWvvP32250+\nRircXszPB14HTovVBPQF5rX9Fg/DkqiKYwgypAgS3Ce+WCaS1ZaBMfX4LIKApiT9gYiYIoiGQ11e\nVyvWvwLv/QWmXA4TL4hf+EscRZDis4grApduKgMV4VIRaFZ6Q0DcNZTNwLkTI+hYsPiUU07hueee\nY9++fQAcPHiQbdu2AfClL32Jl156iX/84x9ccMEFAMyZM4cXX3yRUChEY2MjL7zwAnPmzEna5/Dh\nw/noo4+IRCLU1NTw1ltvxV8rKSmhvr4+YZ0y5X5fevFF5syYhNtLyOzZs3nqqacAeOKJJ+LPn3ba\naTz00EM0NNjdaHbu3Mm+ffvYtWsXhYWFXHTRRcybN4/ly5cn7a+srIzy8vK42nnsscfi6gCI+90X\nLFhAWVkZZWVlbNq0ifHjx3Pdddcxbdo01q9fn/b4qTj77LNbfd61tbWUl5dTWFjI+vXrWbRoUXx7\nn8+Hrre+uZgzZw5PP/00pmmyf/9+5s+fz/Tp09N+dg0NDdTW1nLGGWdwxx13sGrVqrTbdha3iuA+\nKWVcz0kpdwsh/gC8kfEVHUboZnMnTTtrKHOKwK/ZisB5LqikaYFgRuOD6xMzhgCUWH8fI5rZKWW7\nt3xE+bNXExhcgTj9FnuNMUPQliKIZw1p7i4utiJwZwh8VpiQv2/K1xQnayibwWIna0jpmCEYM2YM\nv/nNbzj11FOxLAufz8fdd9/N8OHDKS8v54QTTuCjjz6KX0gmT57MZZddFn985ZVXtooPHHnkkZx/\n/vmMGzeOkSNHJr1+9dVXc/rppzN48CCqnmwONLfc7+WXX86kcaNZu8/d5/+Xv/yFb3zjG9x6662c\nddZZ8edPPfVU1q1bx6xZswA7oPr444+zceNG5s2bh6Io+Hw+7rnnnlb7fOSRR/j2t79NKBRi1KhR\n/P3vf4+/FgwGmTRpErqu89BDDwF2RtQ777yDoiiMHTuWL37xiwQCgZTHT9VSJNXnffrpp3Pvvfdy\nwgkncPzxx8ddSc5nOWHCBCZPnpxk/M4++2wWLlzIxIkTEULwhz/8gUGDBrF+/fqUn119fT1nnXUW\n4XAYKSW33367q8+8Q0gp2/0HLG/xWAU+cvPeTPybMmWK7CzvvPNOp9/bVX7y7Co57TdvSimlPOuv\nC+RFf1vU5X2+vma3HH7dy3L1pzXygfmb5PDrXpa1TdGkbZLO+aXvS/nH4+T/PLNSzvzdf5K2W/LS\nPVLeWCq3f7Kyy+tyaGqslxt/PVEe+uVguWXD2vjz723YL4df97K89f/WyeHXvSyXbj3Y6r0LNx2Q\nw697Wb63Yb+rY+341XFyyZ/OkVK2/T0vf+1hKW8slQsf+GHq119/TMobS+WGlQtavfbRRx+5Wkt7\n1O/fIeXO5dKyrIzsT0op6+rqMravluh6VMqdy2X9gZ0pX4+EQ1LuXC4bDu7J2hpSkc1z7o505HxT\n/a4CS6WLa2ybt15CiJ8KIeqBCUKIuti/emAfh3FFcKYwLBlXBD41M5XFuplcRwBtD4N3XEPhFoPr\nARR/5hXBh/dfxdHmFn6of5f6YHM/QWduQkYVgdBQ2okR7NiwimPfv45PtOOYdNFvUm7jZA1ZZhZT\naaWFlN2/jsChpWuoJTJeIOeFCg8H2vwWpZS/l1KWAH+UUpbG/pVIKftJKX+aozX2WAyrOVisKZmp\nI2iZNZT4XEoSXEOBNK4hPZIZQ7Dk+b8wveZV/ll0Ae9YkwhFmy+srmIECQV4brDacQ01NdZj/OMS\ndKFResmTqfsMkSvXkET2oHYM8WBxuhYTTg+iHmLYPNrGlTmXUv5UCFEuhJguhDjR+ZftxfV07CEr\nsWCxlpkWE9F4sFjE1UbbiiAaVwRJNQSA6rPHNhrRrqePbvrwfcavupk1gQqGn2vfeTfpze0fOhQj\ncJk+agoVRaY3BB8++F2Gm9vYUfkXBh2VvpupM1A+q4qAnmYIYmtNYwik1cN6J3m0iatgsRDiSuBa\n7G6hK4GZwELsITMeaUgMFtsFZZnvPgpuXEM+GqMGRf7kr1sL2IrAbGMKlRtqDx0g+MJl1IkSBn/z\nCfZZdnZOOJUhCLRRUJbg9nKDKTREG4Zg5MEFrCitZErluW3uJ54+muXKYtlz7ABCCCwpSO8a8hTB\n4YRbc34tMA3YJqU8GZgEZL687TDDsGSzayhDMYKk7qMxQxBxoQhqm3TKCpLz19WYa8jUO68IpGWx\n+YGLOcI6wMEz7qffwGHxWEQommAInBhB3DXUullcc0GZS9eQ0FDamDXsJ4oRKG93P4pmr0lmURGI\nHuYagpgJSOsa8hTB4YTbbzEspQwDCCECUsr1wPHZW9bhgW5a8Zx4TVXQMxgjSAoWtxsj8FPXZFBa\nkKwIfDFFYLUxjrA9Pnj8RiaF3mfZ8T9k9PQvAFDotw1BKtdQSZsFZe7bUANYQkVpI0YQkFGkVtDu\nfnISLO6BLZvtzqKpf7d6syJ48skn2b59e76XkVHcGoJPY7MFXgTeFEK8BGzL3rIOD4zYbGGwXUOZ\nyRpqVgQBJ1jswjVUF9YpbaEItFjWkNXJrKENi19j6qa/srz4JGZc0DxkJugYgmgq11BbLSbct6EG\nWxGoaVxD0rIIoCO1NEVkCaixGEE2FUF3CBZ3pA01gERw7oWXp+5tIztXINeRtWSbe++9l0cffbRD\n73nwwQfZt28fRx11FGBXID/55JPtvs/tdvnCbbD4bClljZTyV8AvsGcInNX2uzwMq1kRZKqy2LmA\naoqIp1m2pwhMxUfUsCgNJhsCX8DOopFGJ2IE9XsY+uZ32SYHcvQVDyc1rXNcQ4mGIGLY/3eCxanc\nWR0NFluKljZYHI2GUYREuFAESo4UQb4NgRtaGoLnn3wwZRtqidNWu3u7hswULkiHb3/721xyySUd\n2t8VV1zBD37wg/jjXmUIhBCPOf+XUv5XSvkv4KGsreowQTcTYwRKhnoNSfyq3SrC7yprSMeI5QS0\njBHEFUFHXUOmAc99E81s5BrrR5SVJ1fs2u0vRErXkBOwTj2PIBYsdtlioi1FEAnH4h6+9hWBomVf\nEQhkpxvP5aMN9cknnwwIxkw9KWUb6rv/914Atm3fnpc21OvXr09qy7B161bGjx8PwIgRI7juuuuY\nPHkyzz77LA888ADTpk1j4sSJnHvuuXFD96tf/YrbbrsNsBvHOW2enfYRiZimybx58+LrvO+++wC7\nPfS7775LRUUFd9xxB1u3bmXOnDlMnjyZyZMnx2cQtNxu7dq18e90woQJ7TaqyzZuW0wkddsSQqjA\nlMwv5/DCsKykgrKMDKYxrXhswF3WUJRo7Gtu6RryBx1F0EFD8NZNsO09nh/2C/btHZVyk6BPTQoW\nR2LrDrSx5uaCMnd3zlLRUEl9xxcNNwIgfG4UQcwQtNfA7v+uhz2rXa2tJYFoyB7t6G/R3G/QePji\nLWnfl6821P379yeya238tZb7nTZ1Cl+YcgzFwyfnpQ316NGjiUajbNmyhZEjR/L000/z9a9/Pf56\nv3794v2JqqurueqqqwC44YYbePDBB/n+97+ftL9bbrmFLVu2EAgEUrrCHnzwQcrKyliyZAmRSITZ\ns2dz6qmncsstt3Dbbbfx8ssvA3b76zfffJNgMMiGDRu48MILWbp0aavtvv/973Pttdcyd+5cotFo\nm8olF3iVxVkksY5AU5SMZQ05cQe3weKotF01LRVBIGhfJDvkGlr3Mrx/J0y9gvkFJ7fap0OhX22V\nPhpQFRTFnt+cMkbgDK93rQh8aRVBtMm+6xOxzKi2cLqBZjVGQOdaUeerDbW93maD3HK/X/nSGXYb\naoW8tKEGOP/88+NKoaUhSPz/mjVrmDNnDuPHj+eJJ55g7dq1rfY1YcIE5s6dy+OPP46mtb4/fuON\nN3j00UepqKhgxowZVFdXp7yL13Wdq666ivHjx/O1r32Njz76qNU2ALNmzeJ3v/sdt956K9u2baOg\noP0blmzSpiKQUv4e+H2swdxqYJSU8iYhxFHAoFwssCejt2xDnSFF4KgMt66hiIwpgmDy1+2PuYZw\nawiqN8GL34Ehk+D031P78IpW+3Qo8KmtXEOJSiZl1lC8jqADiiCNIdAjtiFQ/O5jBLK9BnZt3Lm3\nR3S33VCsYPDoDr1P5qkNtU36OoLEiWv5akP99a9/na997Wucc845CCE49tjmosHEFtKXXXYZL774\nIhMnTuThhx+mqqqq1bFeeeUV5s+fz7///W9++9vfsnr16iSDIKXkrrvu4rTTTkt6X8t93XHHHQwc\nOJBVq1ZhWRbBYOobkW984xvMmDGDV155hTPOOIP77ruPz30uf2VZaW+9hBCfjbmAAEqxi8gc51k9\ncHfKN3rEMSyZkD4q4lkxXSGayjXUjiKIpFEEQlGISB+4ab+sN8Ezl4KiwvmPghagrslIqwhauoZa\nGYKULSYsNEW4vqC15RpyDIHqwjWkarEZS92wjiBfbaiBpMH0Lff771deY86MSa6/q0y3oQY4+uij\nUVWVm2++OUkBtKS+vp7Bgwej63rSsR0sy2LHjh2cfPLJ3HrrrdTW1sbXk7jOe+65J95W+pNPPqGx\nsTHp8wLbhTV48GAUReGxxx6Lu3xabrd582ZGjRrFNddcw1lnncWHH37Y7meYTdpSBBb2OMqrgelS\nyslCiBUAUspDQqSZ++cRx0gIFvtUBdOyZwcoHWxFnEjUsOJKwJ0iiBK2bEPQMkYAEEVDmO3ECKSE\nV/4H9q6Buc9CHzt1rrZJZ9SA1ANtWrmGEgxYII0iSFQ7bpCKDy3NzGUjNnVNzaQi6BKdqyPIVxvq\nIUOG8MqT98afb7nfS+dewKRxo9kSdvd9ZaMNNdiqYN68eWzZsiXtsW+++WZmzJjBgAEDmDFjRtIF\nGexA8EUXXURtbS1SSq655ppWmVJXXnklW7duZfLkyUgpGTBgAC+++CITJkxAVVUmTpzIZZddxne/\n+13OPfdcHn30UU4//aNDTlUAACAASURBVPS4Mmm5XSQS4bHHHsPn8zFo0CB+9rOfufocs0ZbrUmB\nitjPD7BbTy+PPR4ArHDT3jQT/3pqG+rpv31T/uTZVVJKKf/69gY5/LqXZVg3urTPqx9dIk+9/b9S\nSikbwrocft3L8t6qjUnbJJ3zbwbL5fd/Rw6/7mUZ0c1W+ztw45Fy0Z0Xt3p+8Yt3yw2/niQ3/rpC\nbrlpvJQ3lkr51m+Stpl40+vyhhdWp1znNx5YKM/53/fij7/92FL5+T9VSSml/Oytb8kfPrWi1Xtu\nfGmNHHfja6lPPAWL/vpNWXvjYCll6+959fwXpbyxVK59/9V293No/267TfWTv231WqbaUId3rpGN\nuz/JyL4cst2SuXHXx7Jp59qUr9Xv2yatTzPbVtsNXhvq9GStDbWU0pm8fCfwAnCEEOK3wALgd1mx\nTIcRiYrACRp3NWCsm7LDWUMhU6HAp8a3T9ofPhQz2TW0fvGbVCy/AVUa1AUGsUsZyD3Gl9Hn/CS+\njWVJ6lK0rXAo8KmtCsri61aVeFvq5HOzXFcVQ0wRpIkROG0ztIALRRBLH3U7/7hz9Iw6giSEsDOd\nUhFzdfWUttoebeMqfVRK+YQQYhlwCra+/aqUcl1WV3YYkDyhzP7ZVUOQmDXkGJe0QWgpwdIJmWqr\n9hLxNQo/itkcLK7e+yl9X72afUp/+n7/bUaW9+exhVu59aW1nBs2OSKWYdMYNbAkafdb4NfSuob8\nmkpET92G2m1VMTiuodQxAjPmGvK5MATxoGAWDYHoiS0m2gkWS88IHDa4rSNA2v2FUs9S80iJHSxu\nzhoCutxvKPGCKoTAr6W+uwbs9hJAyFTS3rkbwo9i2YrANAx2PzSXY2Q9O8/5F0PL+wNQXmSHg2pC\nOkeU2FkQdWH7opleESjJdQSJsY00weKOxghQ0weLzaijCFLPIEjejVNHkGVD0NMunC4UgcfhQfeu\nD+/h2He4zS0mnOe6QsuLZUBNHXgF7IZzQKOhtGovEV+j4osrgsUPz2NcZCWrJ/6Coyd8Jr5N30Lb\nEBxsbHYh1YZsI5Nuv4V+LW36qL3m1hdwPWGimysUDVVIrBTFODJWLe1PM4wmEc3JGrJSGxWZrgNn\nB7DvrXvahTO9IRA90dV1GNPV31HPEGQR3eqAG8cliVlDkD4nH4gbgnpDtKkIVKmz6u1nmPXpQyzu\ncwbTzrk2aZs+MUNwKMEQ1IVtQ9BW+mhLQxBop47AMC3XNQQAQrGPreut6yCsWIwgEEyd1ZSI4hSw\npVAEwWCQ6urqLv+h9URFIBHpL/XS6nTLDI/MIqWkuro6bc2CG1y7hjw6hmlJpCSp6Rx03RAktpiA\n9gyBfbFu0JWUqaMApuKjPLqb4vk/YJM6iglXP9Bqm74x19ChUHMLhtqmmCJoI1gcNSxMS6IqolX9\nQ01TateQ26piAGJpn4aeog4irghctKFWFHSppjQEw4YN49NPP2X//v3u15UCq2YvulpPoLrr0+Ac\nwuFwl/742yNSX43PDKHUpkg7rtuHYhlotbkNFWb7nLsbbs83GAwybNiwTh/HMwRZwrngJw6mAbpc\nVGYHi5svlj5VSW9cYoqgThdpK4BNJcBQuZc6CgnOfZxgYev2wH0K7QvBoVCCImhqWxEkziQoDmit\n6h/STShzO7gegFiPIENv3SPI6Z/kRhEAmCgpDYHP52PkyJHu15SG6I2zWTZ0LhVX39XlfTlUVVW5\n6iHUWRbd+11G736G4E0HWr228g8/pDCyh6N/sTLFO7NHts+5u5Gr8/VcQ1nCueD7EobXQyYUQbIf\nPV3gFWh2DenpXUOmYrcH2PSZPzJ01NiU2wR9KoV+NTlG0NR2jMCZSeAMsHfTYkI3LXwdKLYTcUOQ\nokWGESYqVdQUfWNSYaIi0sQIuoq0LPzCBLVn1WBK1Y8vTcGeYumYIvV379HzyJohEEIcKYR4Rwjx\nkRBirRDi2tjzfYUQbwohNsR+tj9LsAfSsoGaYxC6nD7a0jXUZrDYvlhHpZbehXPi91gy6fdMOvWi\nlK87lBf6kxVBLGuouI1eQwDhqNVq3eljBB0NFtvnZBqtFYEwwkRwf+E1RWrXUCYwnPUpPezCqfpR\nhcQ0Wn8uqhXF8JoLHDZk0zVkAP8jpVwuhCgBlgkh3gQuA96SUt4ihLgeuB64LovryAt6iwZq8ayh\nrqaPGhb+hICqX1PSzyyOKQKd9IbghBmnpXy+JX2L/MnB4iadkqCGmuYOvuW4SnvdanzNqVRM1LQo\n8bn/lRSaowhaxwiEESbagQuViYpoY/5xV9CjYXw0K5geQ2y6mx4No2rJLkNVGhiKZwgOF7KmCKSU\nu6WUy2P/rwfWAUOxJ5s9EtvsEeCr2VpDPnEu+GpC0zmAqJHZ9FE3weIoWloXjlvKi/wcTAgW1zXp\nbe6zeYB9CteQmtp4GVbHKoudOQJmiqZ5ihkh2hFFkCZGkAn0aGx9Ws+6cDqGKxpt7XpTrSimZwgO\nG3ISIxBCjAAmYfcsGiil3B17aQ8wMBdryDWOC0jLsCJolTXkIliso6WNEbilvNBHTSg5fbStfRYk\nKAIppbumc0bHKoudC5WVQhGoZpio0v50Mgc7RpAl11AshiF6WIxAxAyXkWKmtSYNLC9GcNiQ9awh\nIUQx8E/gB1LKusTeJFJKKYRIeYsshLgau/MpAwcOTNlD3A0NDQ2dfm9X2NNoX+g2fvIxVQ2b2FRj\nux2Wr1yFubNzH7slJbop2bljO1VVewCorw1TE5FJ5+icc59DH1IB6FLjk7UriexQU+/YBaFDEfbV\nGvHjbN/ThCJa92N32Fxrn+/iZStp3GYfd+f2rVRV7WLXzihRw+Kdd95J6lVT2xDikAi5/r5qd+1h\nCvDhqhWI4kFJ7ytsaiBiaa73dSwK4VB2fleaavfxRWDvgeqM7j/bv9t1Bw4C8MHC9ygoOyLptZFW\nhJBu5fxvK19/z/kiV+ebVUMghPBhG4EnpJTPx57eK4QYLKXcLYQYjD3trBVSyvuB+wGmTp0qKysr\nO7WGqqoqOvverrBhbz28O5/xY8dQOXEI/XfWwqIFnDBmHJVjOzfTJ2KY8PprHHfMKCorjwHgqR3L\naDrQQGXlSfHt4ue8wYBVoKNyypxZDCtvv8o2HR+aG3hz2yfMnnMiPlXh9yvmM7x/IZWVU1NuP2Rv\nPSyczzGjxzDzuAHwxhuccNwxVM4ZxWpzA3LTJ3z2xJOSU2EXv83QwX2prKxwtaaVxm7YCaOPP46d\ntVbS97x64a+wRIHr7/7T/2oUBv1MzcLvys7N62AFDB56FNMyuP9s/24vqVkD+2HC+PEMHXVC0mu7\nqwz8haXMzvHfVr7+nvNFrs43m1lDAngQWCelvD3hpX8Bl8b+fymH6cjL9MHizscInH26riy2EmIE\nGXANQXMtQa3LGEFT1Iyvr72uqYmjPd3QVoxA66AP2xTZcw2ZsQlwoqfFCBzXUIr0XA0D2dOyoDzS\nks0YwWzgYuBzQoiVsX9nALcAXxBCbAA+H3vcLVj40Dw2/Cb1HW5HcWIBWotgcVfqCJwLZ2IbBp+L\nXkOG0Cj2d038JTaeAztG0JZxSYwROBlCiU3noLUh+P/tnXl4XGd97z+/s8wiyZK8O14SO2CSkMUJ\ncZyF5aoQ9pTQtKWEQGmT1txe2lJa2oZuD32eckvLUi6U9hJISgJpCgUKaZsLhDQmkGbfHDsJISGL\nHduxjS1ZsjQzZ3nvH+ec0SxnRrNIGp2j9/M8fiSNNDPvqyO/3/PbHc9vq6Asah/tx0wWs/wirtl6\njMCfw6whLzxIjYRlDZl28PvzYoTAxkElLOahacycuYaUUj+icd/d183V+3bKQ9/7Chc+fw0QdOFs\ntRCpEU5tsLhcUNaNRRDdWU/7+oNUzAavGWYNZTK5rqaiQXXjOcfzmSx5zYPFrVgEXq0QqLYKyozw\njtV3S0B1Gb6tiuViuVbwZQ6FIKwjMKzW17MQEKuJEChXC0GK0JXFwN6ndvHSOz+Ep4JDaPL4sa5f\nMyoom55HIFWPd0KcRRBk4DQ4wEKLYDZ6s1Q2njtWripuLJZlIXBihKDBiM1221BHIyb9mIIyW5Xw\n27EIxMRoMOSmW8pZQ1bCLIJQaLVFkH4WvRAUJico3vRufDG4d917gsdmQwjCWEDk8y67hrqIEZRd\nLLVN52ZIH83PghBUNp4rzyLoa3ywGYaQtQymSl65ZqDWNVRbS1DZtrsVyhaBVy8EGVXCN1vf91y6\nhqL0VjNhMYLIgnGd6vTRpLbM0DRm0QvBzmt+k5d4z/Dsaz6NvSbotVM4Ptb16043nQtbTBjRPILO\nLQKnxtcefT5TQVmuhQ6cM1HZeG6mPkMRfZmgFXWxxiLIxsQIyrUGbdQRRANl4iyCDCVU2xbBHAlB\nGMw2EiYEph0KbU2dRimsK5A2fr+ahc2iFoJ7v/Fpto3ewl3rr2TLa9+BmQvK6GfFIqjNGrK6H0zj\nuNFrVncf9VUDgQktgr5890KQs036w8ZzM3UejcjbJpMtxgi8yIJqwyKIDirl1WcNZVUJZbVhEcyh\nEERZTYadNCEIDnq/JiurVAwthIQJm6Yxi1YInnrkTrbs/CsezZ7Dtl/7OAB23yAApcnxrl+/Lmso\ndBE1dOO0QCmcxGXXuIagQRC6LASd1w9UMhw2nouG0syUkpoLLYJoz+XBNGHPoUqLYLpbaxtCEI2Y\nrLEIlO+TEwcWiBD4oY/dspN1B21lIiGojhG4UcsJ7RpKDYtSCMaOHib/rV9nTJaw7qobyxlCmVAI\n3EL3FkHDOoIuLIKoT1FtHUHwvTiLIDgg+/tmZ5BH1HiuHddQodIiqGg6V7vmklcfCJ8JM/Rh18YI\niuHgemW3vm8l1pwFi1W4vrRYBE7kGkpYFpSmMYtOCHzP45kvvJtV/mGOvPULLFu1rvy9bH4JAO7U\nRNfvM910TsofRbrrNTQdLK7uPgpQjJnb6zlFPCUMzpIQRI3njk01H1wfMbNraHrN0660dlxD4fvX\nCkEhEAJZKBZBeJBaCXOlWKGQqhohiCyCpBXIaRqz6ITgnq/8BWdP3sUDp/w+p553cdX3sgNDAHiF\n7l1DTszBZhtGd3UEbnVKKgSD4CHeIigWC01bULfL0j67bBHYppCzm//5RHOLowO/Wfpo7US3Voiy\ncFSNEDhTxwEQu/XYiJoPIWjDQlkIRK6hWiGILAIjYa4uTWMWlRDsuvPf2fbTz/HAwAjnv/NP6r6f\n7w+EwC/OgkVQU1AGgdtjVrKGKmIEthW1t47p718qUJqFzqMRSytiBEN5u6phXBx9GZNCXB1BTPqo\n49WL3EzYdrwQFAuTABiZNoTAMDGYGyGILJayBZMQrOj3WxMjKFdKa4sgNSyamcWH9j3LCbe+n73m\nOk7Z/iUkZkh6vj9wDalZEAKvJlgMQUZMVy0mYg7LyO8eZ2k4xSLGLMwiiFjWn2G84PKziWJLr1nn\nGjIbp4/WxlRaIcoainoqRbjFUAjajBGYc2QRRFlNyQsWh64hr9Y1FFkEybJwNI1ZHBaB5zD+5feQ\nUwXUL9/AwGD8dMxMNkdJWVDqXgjiDjbblO4KymoOVGgeLHac4qy7hgCe+9lkS6+Zz1jN6wgqRLG2\nErsVrAYWQSkUAjPTeraUkrmzCKL1RQdrUshk42MEXlhgZiZM2DSNWRxC8P2PcPLUTr665kNsPO3c\npj96XPIYzmTXb1lOH6042CzD6NI1FGYNxaSPluKCxaUiLiZD+dkx/KLGc88faVEIbDPoNVSbPtok\na8iKsdQaYUcHUc3QeS/MGjLbOHiVYWGq7oYGNSQ8SO2EuYbKv98aoZ12DWkhSAupdw2tOHQX7P57\nrndfT+G0y2b8+QI5DOd41+9bbjpX0UTNMqXL9NGwjqCqoCx4/bjRj55TxG8yuL5dosZzMzWci8hn\njOpeQ026j7plkWsnWGzhK0Fqg8WREGQXikUQCkHCLALTsnCVUa5HiYiEwMxoIUgL6bYIfvY0pz7x\nGcaXb+Gj7rs5a93wjE8pGnkst3shiEuHzJhGQ9fQ6GQJpZqLRJxFEOdvj/DcUuAamqUYQdR4Dpo3\nnIvoy1h4vuJ40cUypNwBtWnWUBsWAYCLifJr71gDIbDaChZbmHMVLA7nHFgJazoHwZhTqRGCpBbI\naRqTXiFQCr79fpSY/Nvmj1LC5sx1QzM+rWj0Ybmz4xoSma4jgMgiqD+wRydLXPDXt/Gfj+6v+14l\ncUVXzYLFvlvEFYuc3fmIykqixnMwc1UxUH7fsSmnSrws08CQ6hhBXLptK7iYdRZB5Bqyc21UVIs5\nd0LglSgpMzZBYaHjiFVnEUSVxkmLeWgak7y/zFYRgbd9lt2n/zF3He7jxGV9TbtlRpTMPLY31fXb\nB731q3+9lhGfNbT36BQFx+fJA83rF8ptqI3WgsXKLeHP4hSp4YrfXyuuob5wOM3oZLUQQLDu+PTR\n9uYmuGKW77gj/FIg5Jk2hCCIEcyNEIjn4CbUC+tiITUWVxQ81sHi9JBeIQBYsZnRpWfy6AtjnLl+\nZmsAwDX7yPqzYBF4fl1xlG1K7J37ofHgDmv/WKHue5UE/fqlashMs2AxnjOr4wSjxnMwc3sJmJ5J\nMDblVGU6QX3X1Ci43q5F4MUcVH6Y1WK3ESPAsDCZm2Cx+E5wZ51AHOyGrqGkxTw0jUm3EADjJcXe\no1Oc1YJbCMC1+8mqxhaB7ysOzHBgQ9BErXb+rmUasS0mDo4Hr3fgWCtCUH3Jojvo+F5DJVQbc3tb\nIYoTtGIRVLqGsnatRWBWWQRRH6V2KoshdA3VWAQqjBFk27IITCzmptcQfoItArEwai2CqC5CC0Fq\nSL0QPDsW3Cm3ahH4Vh951fhA/s9H9/Oaj9/O6GR96+NK4g5ty2huEewbbe6SKrn1rzltEVS/7pEX\n97KytBfPHmj6mu0SxQkGW0hJjVxDcRZBMFlttiyCGiEIfdjtuIbEsLHmyiJIsmtI7BjXkLYI0kbq\nheCZY8F/7jNatAiU3U9fE4vg+SOTlFx/xrv3YNpW9d1txoqvIzhY4RpqljlU8lSdrz0b09LZ9132\nXXcFeTXF0Jv+tOk62yWqJWgtfbQyWFwdsK6drNZJiwkIYgRS2zW0I4vAwhCFH+di6xLDd/AS6hry\niLEIwhhBVHCmST6pF4Jnx3xOXtHfcgqlyi4hKw5OqX5OK1BuwXzk+AwWge/XpUI2sggOHgvea7Lk\nMV5s7J5wPL/e1x5XnPXQP3NG8WEe3fLnvOSsi5qus12i6uJ2YgSTJa8+WGxWz1qOq7toBS/GdYFb\npKQsDLP1bCkxgoPajZl21i3iO7hJFQLDxlA1v5NyXYQOFqeF9AvBMb9ltxCAZPoBmDwen8ETuYRG\nJ5sfGHEWQaNeQ4cmpkWnWfyh5Pqx2TfR9wAe+a+v8abxb3Dv8Fs477IPNF1jJyzti1xDrVsEMN0l\nNSJT4xqKa6jXCn6Ma0jcKQrSZmzECNbqzYEQGMrFI3k1BACe2Jh+zU2PG7Q3txI2X0HTmFQLwaHx\nIkcKqqX6gQgjG42rjJ9bHAnATBaB6/t1d7e2KeVJXJUcHC+wbjgofmqWORRlDVViGhLm5Hvse/bH\nbLzj9/gJJ3HW9i80XV+nbFzex2DOaqmgLF9RvxAnYNW9hjq1CMy6gTKGV6BEm4eUOXcWQeAamp1a\njvnGExujVmi9oFBRkx5SLQS7XggO87PWz1xRHGGEw2kKEw2EIHQNHZ1JCDwVEyyujxEopTg0XmTL\nhkCsDow1jk/EBYshOFQLU1Mc/8oVCIonz/pjcn2zGySOuOKCk7jtD0Zami3cl2kiBGa8RWC3aRF4\nYtXFCMQrUmrbIggONn/OhCChFoFhY8a4hkoJ3Y8mnlQLwc69Ywhw+trBlp9j5QIhKDaYWzwWWQQz\nZA25fr1ryDbrB9OMF10Kjs/pa4cQgX2jTVxDXr1rCIJD9aUP/hWb3Z/w9EUfp3/Zuphnzw62abBy\nSWu+4cqK5rjYRmwb6jarb32x6gbKmF4BR9rzX89ljMBQbmKFwDdszFqh9R2chLq6NPGkWggefWGU\nE/qF/mzrZqwdCkFpKn5u8ehUazECx6sPFgcFZdUWQRQoXjecZ8VAtmmMIC4lFeBS40dcbnyfu094\nD+e84d1N1zWfZC2DaHbNTJXFboeVxb5Yda4LwyvhdGoReLNfS2AqF99IpmvIFxurxiIQr5TYdFhN\nPKm+mle+chN35eJdPI3I9IcD7BtZBC1mDQWuodpgcX2MIKohWLUky9qhHPubpKWWXJ++TM0le3E3\nf+p/nscyZ7H1yk81XdN8IyL02SbH47KGGqSPmu3GCAyLbM0ELcsv4BrtWQRGGCOYs2Cx9M/6684H\ncRaB4ZdwEmrhaOJJtRBc9NIVlPa2t8VsXzjAPmZuccHxKDjBgXV0RtdQXPpofdZQVFW8ajDLmqEc\nzxxu3PnUqa0jKByDr76HTP8Qm7d/bUFmceQz8UKQrY0R+IqMacw4/rIWFRMsNv0SbrsV1UYkBHNh\nETiz2vNpPlFGvUVg+KXEuro08aTaNdQJuf5ogH39lLLIGoCZhcCJSR+1Y+YRRBbByoEcJwzlW88a\nCrurcvRZjHdcjz10QtP19IooTjBjjMCt783UCr5YdXestl/Ea9MiEDNyDTW/rp1gKRdlJPOeS5k2\nNrVC4OBqIUgVWghqyIWuIVWstwiiuMCawRxHj89QR+DHtJiI6TV0aLxIxjIYzFusGcoxXnCZaFBU\nVpU1dPc/wOM3w8UfgZNmt2hsNokyh7IzpY/69VlWrRA3R8BWpc6FoEOLYLzgsHPvaOz3ghjBwrPW\nWkEZmboeTEmulNbEo4Wghv6BwCJQpXoXTVRMtmlFPxNFN77RW4jr1Tedi7KGKttIHBwvsmpJFhHh\nhKGgZL9RwLicNfT83XDrX8Cpl8BFv9PeBueZqJZgpvTRUkyNRCv4Uu/DtlUJ3+wwRtBhsPiGu57j\nF//xv5ks1T/fUi4qoQenMjPYca63doPxmgXNnAmBiFwnIgdFZFfFY8tE5FYR+Un4MX6KfA8xLYsp\nlUHihCB0DW1aGQT+mjWei+0UGgpDZcD40HixnI65ZjAQgv0Nagkcz2epGoN//TUYPhHe/g/Qpk99\nvomqi2dyDbkNMqJmIpgjUH1QZVQR32qvD46EPnzf6yxY/MLoFI4XdLqtxcJBmcl0pQSuoRohUC5e\nQmMemnjm0iL4EvCmmseuBm5TSm0Gbgu/XnBMSQ5xYmIEoWto0/JACJrVEsTVEURFWJVxgoPjBVaF\nQrB2hupi13G4Ys9fwtRReMcNkGu9YrpXNLQILAPXV/ihKMbFVFohzjWUwUG1aRFMxwg6swiiNOA9\nR+pnWZh4szoXYl4xM2TERVW4NE1/dgceaXrPnAmBUuoO4EjNw5cC14efXw+8fa7evxumJI8ZM8A+\nqiHYuCIUgiYppK6n6lIhI9eHU/GfKnANBXevqwaDw6uRa+h9/lc5eeIBeOunYM2ZrW6np5QtgkY9\nksI4QaMaiZlQhl3nw86qUtsWgWF2V1l8KMz+ihMCW7kJFoJg3ZWFdpZyEhvz0MQz347L1UqpaDDv\nAWB1ox8Uke3AdoDVq1ezY8eOjt5wYmKi7eeuV1m8ydG65+18soQpcODp3QDced/DlPbE/wqPTxU4\nfPDFqtd45rngP9MP7vgRSzKC4ytGJx0mDu9jx47DAAxm4MHHf8oO84Wq1xt7+m62y79xz8DrmBpb\nB0321Mme54qxnwV3ys8+/RQ7nOfKj+95Nvhd3LbjDvptYf+LBUpTftvr9iYLmMot71n5Pj8nDmMT\nhbZea+y55zkTePKJxzkw2f6hvedwIAB3PfokGyv2CXAhLseOt7eeVpiP63x8NEia2HH7beWJb5v8\nEpMlryd/Ywvpb3s+mK/99iyCpZRSItKw+b5S6hrgGoCtW7eqkZGRjt5nx44dtPvcJ+4coN9wOb/m\ned87+ihLDx7gDf/jIv78zts4YeNmRi44KfY1jB/eyoZ1axgZmb5zf+Ge5+DxXZx/wYWsGszxwugU\nfO+/OP+sUxnZdiIAJz76Q6Q/y8jItunn/fQJluz4DE+ZL2HL+79MLt+8OKmTPc8V3x99lDv3Pc8Z\nLz+Vka0byo/vuetZeGI32y64iJVLsnz52fsoWgVGRl7d1uvf/fi/YE96DAwMMDIyQmFyAu6AoRVr\nuLCN38Fjd03BM3DyyRs541WtPw+CflHjt/6/4Iv+5YyMbK36vne7y5LhZW2tpxXm4zrfvf9OGIVt\nW89laHlw37Z/h0umb5BX9uBvbCH9bc8H87Xf+c4aelFETgAIPx6c5/dviZLZR8arN/HHJh2G8nZ5\niHuzxnPxweJqd8jBY9PFZBFrBqtrCQpTx5m68V0A5K+4cUYRWGg0ixHA9O+iFNOSoxUC19B0jKA4\nFbj0xO7MNaQ6iBGMTjrlXkl7jlQHiz3XxRQFZjJdKRKu23Gmq7ctEuzq0sQy30JwM/De8PP3At+e\n5/dvCdfsI+PXZ3+MTpUY7suQtUwGslbTYLEXO7M4zBoKD43KYrKItcO5KiHY+YX38VLvaX76qk+y\n7uTTOt9Uj8iHLTFq6wiyVvVkNddTdZlFLWFWB4tLxeC6iZ1v62UMK8oaal8IoglzKwYy7DlafQNR\nPkATmjVkWIEQuBWDmmwcVEKFTRPPXKaP3gTcBZwiIntF5CrgY8DrReQnwMXh1wsOz+oj59cHbEcn\nHYbDgSxL++2mjeeCwTT1BWUwPZ83OkCqLIKhHGNTDpMll/u+9fdsO/Lv3LX2vZx98eXdbapHRBZB\nNmZUJUwLgeN1VlmMYWOKwvcDMSgVjocPtykEXVgEUZuQc05cynjBLWeXAeVJd5JQIZA4IVBu21lZ\nmoXNnMUIlFKNTq7XzdV7zha+3U+eetfQ6KTDKWuCXkRL+zJNs4Ycv75AKqojcCosAhFY3j99dxUV\nlT1w748476GPyWALSQAAEvVJREFUsDu7hfN+/RPdbaiH9DXKGjJrhMBX9HdoEQBlIXBCi8Boc7C6\nER7UqnbsZQtElt25Jy3l1sdeZM/RSYb6gtTe8gGa0DvoSAg8p9YiSKawaeLRlcUx+Jl+8qreIjg2\n5TCcD/5jLO3LNOw35PkKpahvOldTR3BwvMjy/kyV5bBmMM8SJtlw6/sYlwFWX3njgmwm1yozxwiC\nA9ztsLI4utOO7uSdQiDgZqb1wfUAZhd1BJFl94oTg/rIyhRSzwn+RqIDNWkYVnDn74ZCoHyfjHig\nLYJUoYUgjkw/GfEoFafFwPF8xotuOVC8rL+xEEQdRuOazsF0HcGh8QIrl1TfuZ4wmOVv7c+znoMc\nfvPnWbFmA0km16SyGCjPJIib39ASRrVF4JYCi8Bs1yKwunANHSvSlzHL1mJldbHjBH9DRkLvoA07\nFIJSsI9S+FESauFo4tFCEINkg//Qk+PTTcSOhe0lhsIYwXCf3bDxXNRCos41FB6GjhsJQbFu2tf6\nH1/Lm837eOjUD3La+W/sdis957Q1SzhxWR/rl1b77KPf4+GJQExdT7U9phLqLQK3GNyNW9n2LAIj\nPNg6EYJDE0G/qKG8zZKcVRUwdhNuEZiRa8gN9lG+OUrofjTxaCGIwQwH2E9VDLCP+gyVLYK+TMPG\nc9G0rfp5BNW9hqKGc2We+2+s2/4STnsb573zz2dpN71l8+ol3PFHP8fygWrBe8nKAWxTeGxfMAmu\n5PnlGEpbGJFvPzjAvfCO1cq0Fyw2Q4sAvxOLoFCuDt+wtK/KNeSHQpDUGEFkEfihayjpMQ9NPFoI\nYjBygRBUzi2OMoSGyllDmfDxevdQef5ug15Djufj+8HQ+rIQjL8I//rrsHQjXPq5Bd9MrlsylsHm\nVUvYvS8Q22CiWwcWQTntM3ANeaXgELZzbcYIukgfrbTsNizLs6fCNRT51o2E3kHXWgRO5BqydIwg\nTSSzN+4cY+WDmQSFCotgLOwzNNwX/MdYFgrBkckSqwar/dFRemht+mhGFflf5rdwbvkGd2csrjaO\nccG+5fCdIXjuR1AYg/d8E3KDc7OxBcYZ6wb5/uMHUUp1nD46ne0THOB+GCOws+0V3kXpo51YBIfG\ni7zmZaEQLO3jB08eQimFiJRHX5p2MmMEZgOLILIUNOlAC0EMdj6IETiT0x1II4tguCJGAPGN56Ks\noNqCspfd/xHOtL/GxGjgtjjLgvx+Ew5I4HO99O9h9emzvJuFy+lrh/ja/Xs5cKzQcdO5cowgDBb7\nYXA2k2vPNWRZ1YLSKlMlj/GiW64F2bCsj4Ljh3GDXPkANRLqSrEyoRCEgjZtESRzP5p4tBDEkMnX\nzy0uC0FF1hAQGzCOzRp68Aayu26C1/whA6/9szlZd9I4fW1g+ex+4RiOpzpKH40qggnnCCgnsAgy\nuTYtguh12hSCqJhs5cC0awiCVhOrluRw3cg1lMw7aCts1aHcYJ9RPYGpLYJUoWMEMWTDcZXu1LHy\nY6NTDiKwJDcdLIb42cVRMLgcLN7/CPznh+DkERj58NwtPGGcdsIgIrBr31jsaM9WiNIYlQosAhVa\nBNkOLQLajBEcKleHTweLAfaGmUOq7BpK5h10ZBGoMEbg6hhBKtFCEEM+HFfpF6ddQ2OTJQZzdnnG\nQBQriGs8F1kEtikwNQpf+1XoXwG/eC0YZt3PL1b6sxabVvSze9+xcDBN+3+O5fz8MFiMGwlBZ8Fi\nfK/5D9ZwsNwvKjgY15eFILBMoiBrUoPF066hYB9eKLTaIkgXWghiyJcH2FfECKacslsIgqyXRo3n\nyjECAb71WzC2F375S4EYaKo4fe0Qj+wJ6jUyXbiGIosAt0BR2UibxWmdxghqO8jmM2bQfO5IZBEk\n25Vih+tWXiQEyXZ1aeLRQhBDvi+IEVCsDhZHqaMRjRrPRa6hTU9+EX58C7zho7BhW93PaYI4QXRX\n3ZFFEN1ph3fyhlug2MFgdcM08ZUgbQrBoYkipiFlVyEEVkFUVBYFWZPaJiSyCKgRAjOjhSBNaCGI\nwTBNjqscVIyrHJ2qF4JlDRrPuZ7PBcZjnLzzU3D6L8D575vzNSeVKGAMdBQjmA7yBgeuuAWKdHbo\nuhgdWARFVgxkMCoyxDYs6yvPJYhcKmZCXUN21KojFLQoC8pqc96DZmGjs4YaMCU5jAohODblcOKy\nar/zcIPGc2riAJ+1P0thcBP5t3029cVh3XD62qHy551kDZm1FoFXxOnAIgDwMNu2CCpnTkeck93H\nyPj/5f5PfpbBqb1AxZ11wrDtqPVG8Hfuh66upO5HE48WggYUJI/hTLcKGJ0slWsIIpb1Z/jp4Ymq\nx5xSkcF/304fBY6//Z/Ih32LNPEs68+wdijHvrFCZxZBuRBsWghKXQgBqr1g8aHxImuGKoRg8ghX\nPP2HuOYRjkwE3Ugfy5zJS1at62hNvUYMg5KywAu7j4YWjrYI0oUWggYUjD4sN7AIfF8xVhMshrAV\ndU0dwf3X/R4XOru4/7yPs3XTlnlbb5J5+doh9o0V6grwWqGclhneyZt+EbdTIRCjI4vgrPWhVeP7\n8M3tZAuHyP7md+hfd25H61hoOFiIV+Ma0hZBqtAxggaUjDx2OLd4vOjiK+pjBP12VeO5h757PRce\n+GfuWXEZWy/ZPu9rTipRnKB2ZkErlNM+wzt50y/iGJ0dUh5mWwVlrufzs+MV/aJ++El46lZ4019D\nSkQAwBELCV1Dquwa0hZBmtAWQQMcM0+fG6Q1jpWrijNwYBfc8bfglnj70UlOscfZ9cnPYIiw+fjD\nPGm/jLN/43O9XHriOGNdcEfdyTwCM0pjDA9wyy/hGl3ECNpwDR05XkIpgoZzT98Ot38UznwHbL2q\no/dfqLhY5WB8FCuwtRCkCi0EDXCtfrKlfUAwtB5gpTEBN10OxTEYPokVyufk7BQqjBc/kzuNlVdc\n03Yx02LngpOX8cbTV3P2icNtPzdq5iYqEoIiBauzpn2BRdC6EERpr+vNUfjGVbDyFPj5T6cuOcDF\nxoiEIIwRZLJaCNKEFoIGBAPsgxTA0UkHwefs+/8YxvfDld+F9eeSA17S22WmgiU5m8+/Z2tHz42y\nhiQ8wG1VxO+wwZsnZllQWuHQeBEbl/Pu/yC4RfiVr0CmvR5HScAVC4lmOXuRELTXwkOzsNFC0ADf\n7idHUDU6OuXwfvPbDL2wA97yCVifHv9v0onSGyOXjq1KeGZnd6s+ZllQWuHgeIEPW//MwMEH4Zf+\nCVZs7uh9Fzqu2BihAOAW8ZRMD/LRpAJ9NRugMgMMqCl2/+/XsMHzuMTaTeHUy8id9xu9Xpqmgihr\nSHwXz3XJqyn8DoXAExMjxiJ44JZryT5yA6JU1eOnug5brMdwz9uOdcZlHb1nEnDFxlBhwZ5XooSN\ntgfShc4aasDSLW/l8eyZiPLIGvDg0MVk3v6Z1Pl/k45VYRHce8OHWco45saLOnotPyZY/MT9t3Hm\nPX/IcOkAoryqf7ZpcN/yn8d640e73sdCxhOrHCPAK+GIvn9MG/qKNuDUba+Hba/v9TI0MxA1RVs3\n/hBnjj7BfcNvYOslnbX08KVaCI4e2s/wf2znsLGcJb/zI9YvWzkra04antiYUQsPv4RDMqetaRqj\nLQJNojEtC18JW/zHec48kTO2X9d259EIX6ZbTHiuy55rr2CZGuX4pdcxtEhFAMA3rArXkBOkk2pS\nhRYCTeJxMZlQeazLv0y+v/OWHr5YGKFFcO/1V3NW4QEeOuNP2Hz2q2drqYnEE5uTSk/z5F+dx+bR\nH+KItgjShpZ2TeJ5YMOvctBcx6Wbu2vp4YuJoTx23v51zn/+i9w3/Ea2/eIHZ2mVycU963Ke2hl8\nvscaYnLDq1nf2yVpZhktBJrEc+Fv/B07duzo+nV8MRkuvcjADz7Ac+ZJnLH92o7dTGni3LdcBW9J\nV7W0phr9V67RhPhisl7tx1Qe1ru+0pWbSaNJEloINJoQJcE86acu+hs2vPTMHq9Go5k/eiIEIvIm\nEfmxiDwlIlf3Yg0aTS3q3Cu55+V/xjlvfG+vl6LRzCvzHiMQERP4HPB6YC9wn4jcrJR6bL7XotFU\ncvbFl/d6CRpNT+iFRbANeEop9VOlVAn4F+DSHqxDo9FoNPRGCNYBeyq+3hs+ptFoNJoesGDTR0Vk\nO7AdYPXq1R2nB05MTMxKamGS0HteHOg9p5/52m8vhOAFYEPF1+vDx6pQSl0DXAOwdetWNTIy0tGb\n7dixg06fm1T0nhcHes/pZ7722wvX0H3AZhHZJCIZ4J3AzT1Yh0aj0WjogUWglHJF5LeB7wImcJ1S\navd8r0Oj0Wg0AT2JESilbgFu6cV7azQajaYaXVms0Wg0ixxRNeP3FiIicgh4rsOnrwAOz+JykoDe\n8+JA7zn9dLvfk5RSMw7TSIQQdIOI3K+U2trrdcwnes+LA73n9DNf+9WuIY1Go1nkaCHQaDSaRc5i\nEIJrer2AHqD3vDjQe04/87Lf1McINBqNRtOcxWARaDQajaYJqRaCtA/AEZENInK7iDwmIrtF5APh\n48tE5FYR+Un4cWmv1zrbiIgpIg+JyH+EX28SkXvCa/3VsH1JahCRYRH5uog8ISKPi8iFab/OIvLB\n8O96l4jcJCK5tF1nEblORA6KyK6Kx2KvqwR8Jtz7ThF5xWytI7VCUDEA583Ay4HLReTlvV3VrOMC\nf6CUejlwAfD+cI9XA7cppTYDt4Vfp40PAI9XfP03wN8ppV4KHAXSNm39/wDfUUqdCmwh2Htqr7OI\nrAN+F9iqlDqDoB3NO0nfdf4S8Kaaxxpd1zcDm8N/24F/nK1FpFYIWAQDcJRS+5VSD4afjxMcDusI\n9nl9+GPXA2/vzQrnBhFZD7wV+GL4tQCvBb4e/kiq9iwiQ8BrgGsBlFIlpdQoKb/OBC1w8iJiAX3A\nflJ2nZVSdwBHah5udF0vBW5QAXcDwyJywmysI81CsKgG4IjIRuAc4B5gtVJqf/itA8DqHi1rrvg0\n8EeAH369HBhVSrnh12m71puAQ8A/he6wL4pIPym+zkqpF4BPAM8TCMAY8ADpvs4Rja7rnJ1paRaC\nRYOIDADfAH5PKXWs8nsqSAtLTWqYiFwCHFRKPdDrtcwjFvAK4B+VUucAx6lxA6XwOi8luAPeBKwF\n+ql3oaSe+bquaRaClgbgJB0RsQlE4Eal1DfDh1+MTMbw48FerW8OeCXwNhF5lsDd91oC//lw6EKA\n9F3rvcBepdQ94ddfJxCGNF/ni4FnlFKHlFIO8E2Ca5/m6xzR6LrO2ZmWZiFI/QCc0Dd+LfC4UupT\nFd+6GXhv+Pl7gW/P99rmCqXUh5VS65VSGwmu6X8ppa4Abgd+KfyxtO35ALBHRE4JH3od8Bgpvs4E\nLqELRKQv/DuP9pza61xBo+t6M/CrYfbQBcBYhQupO5RSqf0HvAV4Enga+NNer2cO9vcqArNxJ/Bw\n+O8tBD7z24CfAN8HlvV6rXO0/xHgP8LPTwbuBZ4C/hXI9np9s7zXs4H7w2v9LWBp2q8z8JfAE8Au\n4MtANm3XGbiJIAbiEFh+VzW6roAQZEI+DTxKkFE1K+vQlcUajUazyEmza0ij0Wg0LaCFQKPRaBY5\nWgg0Go1mkaOFQKPRaBY5Wgg0GkBE3iUiJ/Z6HRpNL9BCoEk9IvLf4ceNIvKumO9fBaxSSj3fwWv/\nSdx7aTRJQqePahYNIjICfEgpdUkbz7HUdG+buO9PKKUGZmN9Gk2v0BaBJvWIyET46ceAV4vIw2Gv\ne1NEPi4i94X93d8X/vyIiPxQRG4mqGZFRL4lIg+E/fG3h499jKA75sMicmPle4XVnx8Pe+k/KiK/\nUvHaOypmC9wYVs4iIh+TYLbEThH5xHz+jjSLG2vmH9FoUsPVVFgE4YE+ppQ6T0SywJ0i8r3wZ18B\nnKGUeib8+kql1BERyQP3icg3lFJXi8hvK6XOjnmvywiqgbcAK8Ln3BF+7xzgdGAfcCfwShF5HPgF\n4FSllBKR4VnfvUbTAG0RaBYzbyDo3fIwQfvu5QRDPwDurRABgN8VkUeAuwkaf22mOa8CblJKeUqp\nF4EfAOdVvPZepZRP0BZkI0Gb5QJwrYhcBkx2vTuNpkW0EGgWMwL8jlLq7PDfJqVUZBEcL/9QEFu4\nGLhQKbUFeAjIdfG+xYrPPSCKQ2wj6Cx6CfCdLl5fo2kLLQSaxcQ4sKTi6+8CvxW28kZEXhYOfKll\nCDiqlJoUkVMJxoJGONHza/gh8CthHGIlwYSxexstLJwpMaSUugX4IIFLSaOZF3SMQLOY2Al4oYvn\nSwRzDDYCD4YB20PEjz78DvA/Qz/+jwncQxHXADtF5EEVtMOO+DfgQuARgg6xf6SUOhAKSRxLgG+L\nSI7AUvn9zrao0bSPTh/VaDSaRY52DWk0Gs0iRwuBRqPRLHK0EGg0Gs0iRwuBRqPRLHK0EGg0Gs0i\nRwuBRqPRLHK0EGg0Gs0iRwuBRqPRLHL+PzQRlhC8LV/yAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F1lQIWRPS_CD",
        "colab_type": "text"
      },
      "source": [
        "The transition matrix has changed so here we have to define a new one.\n",
        "\n",
        "In this Question the transition matrix is equal to:\n",
        "\n",
        "\n",
        "$B = \\begin{pmatrix}\n",
        "        (1-\\kappa) \\epsilon+ \\frac{\\kappa}{N}&\\frac{\\kappa}{N}& \\frac{\\kappa}{N}&.&(1-\\kappa)(1-\\epsilon) + \\frac{\\kappa}{N}\\\\\n",
        "        (1-\\kappa)(1-\\epsilon) + \\frac{\\kappa}{N} &(1-\\kappa) \\epsilon+ \\frac{\\kappa}{N}&\\frac{\\kappa}{N}&.&\\frac{\\kappa}{N}\\\\\n",
        "         \\frac{\\kappa}{N} &(1-\\kappa)(1-\\epsilon) + \\frac{\\kappa}{N}&(1-\\kappa) \\epsilon+ \\frac{\\kappa}{N}&.&\\frac{\\kappa}{N}\\\\\n",
        "         .&.&.&.&\\frac{\\kappa}{N}\\\\\n",
        "         \\frac{\\kappa}{N}&.&.&(1-\\kappa)(1-\\epsilon) + \\frac{\\kappa}{N}&(1-\\kappa) \\epsilon+ \\frac{\\kappa}{N}\\\\\n",
        "     \\end{pmatrix}\n",
        "     $"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AdcwhByVXR13",
        "colab_type": "code",
        "outputId": "46358e2b-ac41-4a12-8b82-1b47ab94154a",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 134
        }
      },
      "source": [
        "N = 50\n",
        "eps = 0.3\n",
        "kappa = 0.1\n",
        "col = np.ones(N)*(kappa)/N\n",
        "col[0] += (1-kappa)*eps\n",
        "col[1] += (1-kappa)*(1-eps)\n",
        "row = np.ones(N)*(kappa)/N\n",
        "row[0] += eps*(1-kappa)\n",
        "A_2 = toeplitz(col, row)\n",
        "A_2[0][N-1] = (1-kappa)*(1-eps) + kappa/N\n",
        "print(A_2)"
      ],
      "execution_count": 67,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[[0.272 0.002 0.002 ... 0.002 0.002 0.632]\n",
            " [0.632 0.272 0.002 ... 0.002 0.002 0.002]\n",
            " [0.002 0.632 0.272 ... 0.002 0.002 0.002]\n",
            " ...\n",
            " [0.002 0.002 0.002 ... 0.272 0.002 0.002]\n",
            " [0.002 0.002 0.002 ... 0.632 0.272 0.002]\n",
            " [0.002 0.002 0.002 ... 0.002 0.632 0.272]]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mhZ9_FBQS1cY",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def forwardpass_2(A,w,observations):\n",
        "  ###Initialisation\n",
        "  B = np.ones(A.shape)*((1-w)/N) + np.eye(A.shape[0])*w\n",
        "  pi = np.ones(N)/N\n",
        "  log_alpha_kk1 = np.log(pi)\n",
        "  log_alpha_kk = log_alpha_kk1 + np.log(B[observations[0]-1])\n",
        "  \n",
        "  for i in range(1,len(observations)):\n",
        "    v_max = np.max(log_alpha_kk)\n",
        "    log_alpha_kk1 = v_max + np.log(A.dot(np.exp(log_alpha_kk-v_max)))\n",
        "    p_yk_xk = B[observations[i]-1]\n",
        "    log_alpha_kk = np.log(p_yk_xk) + log_alpha_kk1\n",
        "  result = np.exp(log_alpha_kk)\n",
        "  return result\n",
        "\n",
        " \n",
        "def backwardpass_2(A,w,observations):\n",
        "  ###Initialisation\n",
        "  B = np.ones(A.shape)*((1-w)/N) + np.eye(A.shape[0])*w\n",
        "  log_beta_kk1 = np.zeros(N)\n",
        "  log_beta_kk = np.log(B[observations[len(observations)-1]-1]+ log_beta_kk1)\n",
        "  for i in range(1,len(observations)):\n",
        "    v_max = np.max(log_beta_kk)\n",
        "    log_beta_kk1 = v_max + np.log(A.T.dot(np.exp(log_beta_kk-v_max)))\n",
        "    log_beta_kk = log_beta_kk1 + np.log(B[observations[len(observations)-i-1]-1])\n",
        "  return np.exp(log_beta_kk1), np.exp(log_beta_kk)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mzYkCsovPTDL",
        "colab_type": "text"
      },
      "source": [
        "Now, at each step the robot can be kidnapped with a probability k and so he lands anywhere in the map with a uniform distribution.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "d7Tdo0A3Ymlt",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "\n",
        "\n",
        "$\\mathbb{P}(x_{k}|y_{1:k})$ est égal à $\\alpha_{k|k}(x_k)$ à une constante près mais comme c'est une densité de probabilité on peut simplement renormaliser $\\alpha_{k|k}.$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nU3itbI0YpCB",
        "colab_type": "code",
        "outputId": "cf981ea5-ab66-4fa4-f560-baea3d017ce9",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 251
        }
      },
      "source": [
        "alpha_kk = forwardpass_2(A_2, 0.3,[1,2,3,4,5,6])\n",
        "PxkY1k = alpha_kk/np.sum(alpha_kk)\n",
        "print('Voici P(xk|y1:k):')\n",
        "print(PxkY1k)"
      ],
      "execution_count": 69,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Voici P(xk|y1:k):\n",
            "[1.50709059e-04 1.50855273e-04 1.58388833e-04 5.37037245e-04\n",
            " 1.95367573e-02 9.72835153e-01 1.50706793e-04 1.50706793e-04\n",
            " 1.50706793e-04 1.50706793e-04 1.50706793e-04 1.50706793e-04\n",
            " 1.50706793e-04 1.50706793e-04 1.50706793e-04 1.50706793e-04\n",
            " 1.50706793e-04 1.50706793e-04 1.50706793e-04 1.50706793e-04\n",
            " 1.50706793e-04 1.50706793e-04 1.50706793e-04 1.50706793e-04\n",
            " 1.50706793e-04 1.50706793e-04 1.50706793e-04 1.50706793e-04\n",
            " 1.50706793e-04 1.50706793e-04 1.50706793e-04 1.50706793e-04\n",
            " 1.50706793e-04 1.50706793e-04 1.50706793e-04 1.50706793e-04\n",
            " 1.50706793e-04 1.50706793e-04 1.50706793e-04 1.50706793e-04\n",
            " 1.50706793e-04 1.50706793e-04 1.50706793e-04 1.50706793e-04\n",
            " 1.50706793e-04 1.50706793e-04 1.50706793e-04 1.50706793e-04\n",
            " 1.50706793e-04 1.50706793e-04]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DaBUGzHAY5ft",
        "colab_type": "text"
      },
      "source": [
        "$\\mathbb{P}(x_{k}|y_{1:M})$ est proportionnel à $\\gamma_{k|k}(x_k) = \\alpha_{k|k}(x_k)*\\beta_{k|k+1}(x_k)$ donc on peut procéder de la manière suivante de façon analogue à juste avant. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ne5_T-Z5Y5nj",
        "colab_type": "code",
        "outputId": "be781b6d-f107-4b8b-e169-9293262332db",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 251
        }
      },
      "source": [
        "alpha_kk = forwardpass_2(A_2, 0.3,[1,2,3,4,5,6])\n",
        "beta_kk1 = backwardpass_2(A_2,0.3,[7,8,9,10,11,12])[0]\n",
        "PxkY1M = alpha_kk*beta_kk1/np.sum(alpha_kk*beta_kk1)\n",
        "print('Voici P(xk|y1:M):')\n",
        "print(PxkY1M)"
      ],
      "execution_count": 70,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Voici P(xk|y1:M):\n",
            "[1.41776525e-04 1.41914073e-04 1.49001119e-04 5.05207019e-04\n",
            " 1.83788126e-02 9.15175162e-01 4.08039881e-02 1.83788126e-02\n",
            " 5.05207019e-04 1.49001119e-04 1.41914073e-04 1.41776525e-04\n",
            " 1.41774393e-04 1.41774393e-04 1.41774393e-04 1.41774393e-04\n",
            " 1.41774393e-04 1.41774393e-04 1.41774393e-04 1.41774393e-04\n",
            " 1.41774393e-04 1.41774393e-04 1.41774393e-04 1.41774393e-04\n",
            " 1.41774393e-04 1.41774393e-04 1.41774393e-04 1.41774393e-04\n",
            " 1.41774393e-04 1.41774393e-04 1.41774393e-04 1.41774393e-04\n",
            " 1.41774393e-04 1.41774393e-04 1.41774393e-04 1.41774393e-04\n",
            " 1.41774393e-04 1.41774393e-04 1.41774393e-04 1.41774393e-04\n",
            " 1.41774393e-04 1.41774393e-04 1.41774393e-04 1.41774393e-04\n",
            " 1.41774393e-04 1.41774393e-04 1.41774393e-04 1.41774393e-04\n",
            " 1.41774393e-04 1.41774393e-04]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "308wZLiVZELV",
        "colab_type": "text"
      },
      "source": [
        " $argmax_{x_k}\\mathbb{P}(x_k|y_{1:k})$ se calcule en prenant l'argmax du résultat obtenu deux cellules plus haut.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "714VYbaxZERv",
        "colab_type": "code",
        "outputId": "76a9a091-5ab3-490b-bd2e-4b1341b4e1de",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 33
        }
      },
      "source": [
        "print('argmax_xk P(xk|y1:k)=' + str(np.argmax(PxkY1k)+ 1))"
      ],
      "execution_count": 71,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "argmax_xk P(xk|y1:k)=6\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}