{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bundle adjustment \n",
    "\n",
    "General formulation: Given a set of images, reconstruct a 3D scene:\n",
    "* Estimate feature in given images\n",
    "* Find correspondences by matching the features\n",
    "\n",
    "Given point correspondences estimate 3D location of the observed points + position ans orientation of the camera + camera internal parameters.\n",
    "\n",
    "And do that such that the **reprojection error** (my version -- projection error, since there is only projection from 3D into 2D is involved) is minimized\n",
    "\n",
    "$$ E = \\sum_j\\sum_i || y_{ij} - f(X_i, K_j, p_j) ||$$\n",
    "where $y_{ij}$ - 2D pixel coordinate of feature i in camera j, $X_i$ - 3D coordinate of a point that should ideally be projected in $y_i$, $K_j$ - camera's j internal parameters, $p_j$ - external parameters of camera j (position and orientation).\n",
    "\n",
    "To be robust against wrong correspondences put robust function \\rho. Check Huber, McClure etc.:\n",
    "$$E = \\sum_j\\sum_i \\rho(|| y_{ij} - f(X_i, K_j, p_j)||)$$\n",
    "\n",
    "In general, bundle adjustment is a maximum likelihood estimation problem. So, further is a simple example of solving the maximum likelihood problem.\n",
    "\n",
    "### Fitting a line using maximum likelihood estimation\n",
    "Given an array of 2D points fit the line.\n",
    "Let's denote the points that we observe with $y_i$ and the function that generates points $f(x_i)$ and $e_i$ additive gaussian noise.\n",
    "Then we have the following formulation $ yi = f(x_i) + e_i $.\n",
    "\n",
    "Since we are fitting a line than $f(x_i|a,b) = ax_i + b$. We also want to maximize the likelihood of the observed data, given we select the parameters for $a,b$.\n",
    "$$p(y_1..y_n|a,b) = \\prod \\frac{1}{\\sqrt{2\\pi\\sigma^2}} \\mathrm{exp}(-\\frac{(y_i-f(x))^2}{2\\sigma^2})$$\n",
    "$$-log(p(y|a,b)) = \\frac{1}{\\sqrt{2\\pi\\sigma^2}} \\sum_i-(-\\frac{(y_i-f(x_i))^2}{2\\sigma^2})$$\n",
    "$$ min_{a,b} L(y| a,b) \\propto \\sum_i (y_i - ax_i -b)^2$$\n",
    "\n",
    "**Regular solution **\n",
    "\n",
    "Now we need to minimize L with respect to parameters (a,b). For that we take the derivative with respect to a and b and set it to zero.\n",
    "$$\\frac{\\partial L}{\\partial a} = \\sum(y_i - ax_i -b)(-x_i)$$\n",
    "$$\\frac{\\partial L}{\\partial b} = \\sum(y_i - ax_i -b)(-1)$$\n",
    "\n",
    "After some manipulations we can get the system of linear equation.\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "    \\sum x_i^2 & \\sum x_i\\\\\n",
    "    \\sum x_i & N\n",
    "\\end{pmatrix}  \n",
    "\\begin{pmatrix}\n",
    "    a\\\\\n",
    "    b\n",
    "\\end{pmatrix} \n",
    "=\n",
    "\\begin{pmatrix}\n",
    "    \\sum y_ix_i \\\\\n",
    "    \\sum y_i\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "\n",
    "** Elegant solution **\n",
    "\n",
    "Consider $y = (y_1, .., y_n)$ be a vector and $x = [[x_1, 1], [x_2, 1], ..., [x_n, 1]]$ $nx2$ a matrix and $\\theta = (a,b)^T$ line parameters. Then given the fact that we fit a line, we minimize the following likelihood $ L = || y - x\\theta ||_{\\sum_i} $. Taken into account gaussian noise assumption\n",
    "$$L = (y-x\\theta)^T \\cdot (y-x\\theta)$$\n",
    "$$L = y^Ty - 2y^Tx\\theta + \\theta^Tx^Tx\\theta $$\n",
    "Then considering the fact that $(a^Ta)' = 2a$. the derivative over $\\theta$ will have the following form\n",
    "$$\\frac{\\partial L}{\\partial\\theta} = -2y^Tx \\cdot 2x^Tx\\theta$$ and $\\theta = (x^Tx)^{-1}y^Tx$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ML estimate of (a,b) for original data [1.74090899 1.03956472]\n",
      "ML estimate of (a,b) for contaminated data [1.37730457 1.57029592]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def mlEstimation(x,y):\n",
    "    x_sq = np.square(x)\n",
    "    sum_x2 = np.sum(x_sq)\n",
    "\n",
    "    sum_x = np.sum(x)\n",
    "    sum_y = np.sum(y)\n",
    "\n",
    "    sum_xy = np.dot(x,y)\n",
    "\n",
    "    m = np.ndarray(shape=(2,2), dtype=float, order='F')\n",
    "    m[0,0] = sum_x2\n",
    "    m[0,1] = sum_x\n",
    "    m[1,0] = sum_x\n",
    "    m[1,1] = x.shape[0]\n",
    "\n",
    "    b = np.array([sum_xy, sum_y])\n",
    "\n",
    "    sol = np.linalg.solve(m, b)\n",
    "    \n",
    "    return sol\n",
    "\n",
    "\n",
    "x = np.array([4.75, 5.30, 5.20, 2.75, 4.59, 1.20, 4.94, 0.22, 1.94, 0.32, 0.68, 5.76])\n",
    "y = np.array([9.39, 9.95, 10.52, 5.27, 8.96, 3.15, 9.71, 1.69, 4.19, 1.74, 2.23, 11.22])\n",
    "\n",
    "sol = mlEstimation(x,y)\n",
    "print \"ML estimate of (a,b) for original data\", sol\n",
    "\n",
    "\n",
    "x = np.append(x, 5.70)\n",
    "y = np.append(y, 2.10)\n",
    "\n",
    "sol_cont = mlEstimation(x,y)\n",
    "print \"ML estimate of (a,b) for contaminated data\", sol_cont\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABAkAAAElCAYAAACGUbtrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XlcVOXiBvBnNpYEV3BFEFldk1IUNXM38V67ptXVzLQyc8sUdzOXLLUQ1yy1zCWrW2a3bpJKivlTQRNFTAQBATU1xAVBEWY5vz/mMlc4B2UZZjnzfD+f++HOGebM+z7viff4nnPeVyEIggAiIiIiIiIicnhKaxeAiIiIiIiIiGwDBwmIiIiIiIiICAAHCYiIiIiIiIjovzhIQEREREREREQAOEhARERERERERP/FQQIiIiIiIiIiAsBBAiIiIiIiIiL6Lw4SENmY0aNHY+TIkRX+/S1btsDLy8usZVAoFPj1118r/Pu//vorFAqFWctAREREtqsmzj/MpXv37li4cGGFf1+n00GhUODgwYM1ViYie8JBAqIakJqaihdeeAEeHh5wcXFBu3btsHbtWgiC8MjPrl69Gh9//HGFv+vFF1/EqVOnqlNci/vss8/QokULaxeDiIjIorKysjB69Gg0bdoULi4uCAwMxFtvvYXLly+b7Tu8vLywZcsWs+2vPOY+/7Cnf6inp6dDoVAgKyvL2kUhqhEcJCAys7Nnz6Jz584wGAzYs2cPUlNTMWPGDCxatAgTJkwo93NarRaCIKBOnTqoU6dOhb/P1dUVnp6e5ig6ERER1ZDU1FR07NgRN27cwL/+9S+cP38eW7duhU6nw8qVK61dvErj+QeRfHGQgMjMJk+ejODgYHz33Xfo2LEjfHx8MGrUKGzbtg2ffvopDh8+DAA4ePAgFAoF9uzZgzZt2sDV1RU3btwQPW6QlZWFnj17wsXFBR06dMC3335bavS67O1+JZ9/5513UL9+fTRt2hRRUVGm94uKijBq1Cg0b94ctWrVwpNPPokDBw5Uqo5nzpxBx44d4eLigu7du4tG0o8ePYpevXqhbt268PT0xPDhw5Gbm2uq99ixY5GdnQ2FQmG6amCOchEREdmqiRMnws/PDz/99BOeeuopeHt7IywsDOvXr8f8+fNNv/fRRx+hefPmcHZ2RpcuXXD8+HHTeyV9/s6dO+Hr64u6devi1VdfRVFREQCgZ8+e+PPPPzFmzBgoFAr07NkTAPDTTz+hS5cucHd3R9OmTTFhwgTcvXvXtN+FCxeie/fuWL16NZo0aYK6devigw8+QFFREd58803Url0b/v7+iImJEZWlRHXPP/z9/QEAvXr1gkKhwOjRowEAer0e8+fPh5eXF9zd3dGzZ08kJSWZPicIAubNm4d69erB09MTH3300SPbIi8vD0OHDoWrqysCAwOxb9++Uu//9ddfGDZsGBo3bgx3d3f06NEDiYmJpvcDAgIAAL6+vlAoFKZHG5YuXYpWrVrhscceQ0BAANasWfPIshDZIg4SEJlRbm4uYmNj8fbbb4ue0Q8PD0dQUBC+/fbbUtsXLVqETZs24cyZM6hdu7ZonyUDBseOHUNkZGSpE4ny/PTTT9BqtYiPj8fChQsRERFh6lB1Oh0CAwPxn//8B6dPn8bgwYPx7LPPIicnp0J11Ov1eO655+Dt7Y2EhARMmTJFVKaCggKMHz8eJ06cwC+//IJLly6Z7qLo2rUrVqxYAS8vL1y9ehVXr15F165dq10uIiIiW5Wbm4sDBw5g2rRpknP41K1bFwDw1VdfYeHChVi2bBkSExPRvn17hIeH486dO6bfvXHjBrZu3YqffvoJP/zwA3788Uds3LgRALBr1y40adIEq1atwtWrV7Fr1y4AwP379zFv3jycPn0a33zzDWJjY7Fo0aJSZUhKSkJiYiJiY2OxatUqzJs3D4MHD0abNm2QkJCAAQMGYNSoUSguLi63ntU5/4iPjwcAfP/997h69SpWr14NwHieFB0dja+//hqnTp1Ct27d0K9fP1Mm27Ztw5o1a7Bx40YcPHgQcXFxOH369EPb4+2338bZs2dx4MABbN26FQsWLCj1fmFhIXr06IGYmBgkJCSgdevWGDx4MO7fvw8AiIuLAwAcP34cV69exfTp0wEAzs7O2LRpE86ePYv3338fc+fORXR09EPLQmSTBCIym/j4eAGAcOrUKcn3Bw8eLAwcOFAQBEGIjY0VAAgHDx4s9TuvvPKK8NJLLwmCIAh//PGHAEBIS0szvb9hwwYBgJCZmSkIgiB88cUXQrNmzUp9vnXr1qX2GRgYKKxdu7bccgcFBQlbt241vQYgxMTESP5udHS04OrqKty8edO0bdasWcLD/pzExcUJarVa0Ol0giAIwqZNmwQfH59yf7+8chEREdmjR50flOjcubMwY8YM02utVit4eXkJ69atEwTB2OcrFArh2rVrpt954403hKFDh5peN2vWTPjiiy8e+j1ff/214Ovra3q9YMECoV69esL9+/dN24KCgoRBgwaZXl+9elUAICQlJZnKYs7zD61WKwAQYmNjTe8XFhYKrq6uwpkzZ0p9LiAgQNi+fbsgCIIQGhoqzJo1y/TezZs3BVdXV2HBggWS35mXlyeo1Wrhl19+MW375ZdfRN/9IJ1OJ9SqVUv47bffBEEQhLS0tFLnYuUZN26cMGbMmIf+DpEt4p0ERGYkVGBiwrJCQkLKfS8tLc10i1+Jjh07PnKfbdu2LfW6cePGpa7IR0ZGon379qhfvz7c3NyQlpaGS5cuVai8qamp8Pf3R7169UzbQkNDS/3O5cuX8fLLL6Nly5Zwd3dHnz59oNPpcO3atYfuuzrlIiIisnepqano0qWL6bVarUbHjh2Rmppq2ubp6YlGjRqZXpft46UkJydjyJAh8Pb2hru7O8aMGSPqXwMCAuDs7Gx63ahRI7Rp06bUawC4fv16ud9j7vOPjIwMFBYWokuXLnBzczP9LyMjAxcuXABgzOzB85B69eqVOm8q68KFC9DpdKU+U/Y8RqvVYu7cuWjVqhXq1q2LOnXq4N69e488J9m9eze6d++ORo0awc3NDZs3b+Z5DNkltbULQCQnJZ3SuXPn0KFDB9H7KSkpGDBgQKltjz32WLn7EwShSksLajSaUq8VCgUMBgMA4Msvv8TixYuxdu1adOjQAbVq1cKQIUOg1WortO+KlGn06NEoLi7Gxo0b4eXlhczMTISHhz/0O6pbLiIiIlvl5+cHhUKB1NRUyfODynhYH1+ewYMHo3379tixYwcaNmyIQ4cO4Y033njkfh/cVtL3P+y7zH3+UVBQAMA4n1HJIxkl6tevLypbRZRc0HnYZ5YvX46tW7dizZo1CAoKgouLC0JDQx9a1gsXLuC5557DrFmzsGrVKtSpUwfLly9Henp6hctGZCt4JwGRGXl4eODpp5/G6tWrRXcVREdH4/z583jhhRcqvL/AwEDk5eUhIyPDtC0hIaFaZYyPj0fv3r3xyiuv4PHHH0fjxo1x8eLFCn8+KCgIaWlpuH37tmnb77//LvqOadOmoW/fvggODjZNWlhCo9FAr9ebtVxERES2ysPDA7169cKqVask7zrMy8sDYOxjS57NB4zP8Z84cQLBwcEV/q6yfWxubi4yMjLw7rvv4qmnnkJQUNAj7+yrCY/q51UqFZRKZamyt2rVCk5OTrh69Sr8/f1L/a9kkCAwMLDU5I63b99+6D/M/fz8oFarS31G6jzm+eefx9ChQ9G2bVs4Ozvj1q1bpvdLBkMeLOvJkyfh6uqKxYsXo2PHjggICEBmZmZlYyKyCRwkIDKztWvX4ty5c3jhhReQkJCA7OxsbN++HaNGjcKbb76J7t27V3hfbdq0Qbdu3TB27FgkJSVh//79ppmCq3KHAWDsHI8ePYr/+7//w9mzZ/HKK6888grEgwYMGIAmTZrg9ddfR3JyMnbu3ImtW7eKvmP79u1IS0vDnj178MEHH5R638fHB3/99RdOnDiB3NxcaLXaapeLiIjIlq1btw6pqano27cv9u3bh6ysLBw7dgyTJ0/G4sWLAQBTpkzB+vXr8dVXXyElJQUTJkxAYWFhqVWPHsXHxweHDh3CtWvXkJeXh3r16qFevXrYtGkTLly4gH/961/YsGFDTVWzXI/q5xUKBZo3b44DBw4gJycHBQUFqF27NiZNmoTx48fj+++/R2ZmJuLi4jB37lycPXsWADB+/Hh8/PHH2LlzJ5KTk/H6669DpVKVW47atWtjxIgRmDp1Ko4dO4b4+HjRBMx+fn7Ys2cPTp48iZMnT+KVV16Bi4uL6f3GjRvDyckJ+/btw/Xr13Hv3j34+fnhzp072LJlC9LT07FkyRLR4AORveAgAZGZtWvXznQVoF+/fggMDMTy5cuxYMECrF+/vtL72759O/R6PUJDQzFt2jTMnj0bAEo9N1gZb775Jvr06YPw8HD069cPTz31FB5//PEKf16lUmHXrl3IzMxESEgIoqKiTEv/lPjss8+Qnp6Odu3aYf78+ViyZEmp93v06IF//vOf6Nu3Lzw9PXHkyJFql4uIiMiWtWrVCidOnICXlxdeeeUVBAcHY+TIkVAoFJg2bRoAYPjw4ViwYAFmzpyJxx9/HElJSYiOjpZc/ag8CxcuxLFjx9C8eXM8++yzUKlU2LFjB/bt24c2bdpgw4YNpkEJS6pIP//hhx9ix44daNKkCSZNmgTAuCTkhAkTMH36dAQFBeGFF17ApUuX0KBBAwDGRxwnTpyI119/HT169EDHjh0fef6watUqBAUF4emnn8ZLL70kGiR455134Ovri+7du2Po0KF44403TN8HGM/BPvroIyxevBiNGjXChx9+iJCQELz//vuYOXMmnnjiCWRlZWHcuHHmiI7I4hRCVWZaIyKr+fLLLzF+/Hjk5eVBqeQ4HxERERERmQ8nLiSycTExMdDpdAgODkZycjLmzp2LESNGcICAiIiIiIjMjoMERDbu/v37mDFjBrKzs+Hp6YkhQ4Zg2bJl1i4WERERERHJEB83ICIiIiIiIiIAnLiQiIiIiIiIiP6LgwREREREREREBICDBERERERERET0X1abuNDZ2Rmenp7W+noiIiKru379OoqKiqxdDNnjOQcRETm6ypxzWG2QwNPTE5cvXzbLvgRBQH5+Ptzd3aFQKMyyTzlgLmLMRBpzEWMm0piLWHUy8fLyqqFS0YN4zlHzmIsYM5HGXMSYiTTmImapcw5ZPG6g0+kQGxsLnU5n7aLYFOYixkykMRcxZiKNuYgxE8fC9pbGXMSYiTTmIsZMpDEXMUtlIotBAiIiIiIiIiKqPg4SEBEREREREREAGQ0SqNVWm17BpjEXMWYijbmIMRNpzEWMmTgWtrc05iLGTKQxFzFmIo25iFkiE4UgCEKNf4sELy8vs00iREREZI/YF1oGcyYiIkdXmb5QFncSGAwG5OTkwGAwWLsoNoW5iDETacxFjJlIYy5izMSxsL2lMRcxZiKNuYgxE2lyy0UQgMOHgS1bjD+rcqneUpnIYpBAr9cjLi4Oer3e2kWxKcxFjJlIYy5izEQacxFjJo6F7S2NuYgxE2nMRYyZSJNTLtnZQKtWQJ8+wOTJxp+tWhm3V4alMpHFIAEREVFlVWZE/07RHexN32uxshEREZE8CAIwYACQkQEUFwMFBcafGRnAM89U7Y6CmsaZIIiIyOFkZxs77MxMwMnJ2Fn7+gJ79wI+Pv/7vSJdET458Qne/7/3kV+Uj7TJaWhep7n1Ck5ERER25cgRICsL0OlKb9fpgAsXjO93726VopVLFncSKBQKuLu7Q6FQWLsoNoW5iDETacxFjJlIk0MuFRnR1xv02HZ6G4LWBWHq3qlwUbtg/aD1aOLeRLQ/OWRCFcf2lsZcxJiJNOYixkykySWX9HRAo5F+z8nJ+H5FWSoTrm5AREQO5fBhoG9foKhI/J7GScCSf0Xjy6tzcCbnDOq51MOc7nMwKXQSXDWuZi8L+0LLYM5ERGQthw8b5yAoLha/5+QE7N9vmTsJHHJ1g+zsbNnMfGkuzEWMmUhjLmLMRJoccil3RN8rDvqXn8as039D+s10zO42GxemXMCMbjMeOkAgh0yo4tje0piLGDORxlzEmIk0ueTSrZvxkUZ1mQf91WqgZUvj+xVlqUxkMUig1+uRmJgoi5kvzYm5iDETacxFjJlIk0Mu/v5lRvM9k4F//gN4vSsMzY5icNM3kP5WOpb2XYq6LnUfuT85ZEIVx/aWxlzEmIk05iLGTKTJJReFwjjnkZ+f8c4BNzfjT39/4/bKPDlgqUw4cSERETmUkhH99OuXoH9qAfD4VkBpgCJ5GFpkLsG/jwVVqsMmIiIiehgfH+DcOeMkhenpxgGCbt0qN0BgSRwkICIih3Lr/k08vWQpzietBVRFUF3sBcX+ZfB3Da30iD4RERFRRSgUxrkHbG0lAymyGCRQKBTw9PS0+5kvzY25iDETacxFjJlIs+dc7mnvYXX8aiw/shx5RXno0KwDRjZehvre/REwRlHlEX17zoQqj+0tjbmIMRNpzEWMmUhjLmKWyoSrGxARkaxp9VpsPrUZi35bhKsFV+Fb1xfv934fL7Z9EUqFdafmYV9oGcyZiIgcncOtbqDX65GSkmL3k1qYG3MRYybSmIsYM5FmT7kIgoCdyTvR9pO2eHP3m9ALeqwduBYpk1IwvN1wsw0Q2FMmVH1sb2nMRYyZSGMuYsxEGnMRs1QmshgkMBgMSE1NtfvlMcyNuYgxE2nMRYyZSLOXXA5kHkDnzzrj+e+ex5X8K1jUcxEy3srApNBJcFI5mfW77CUTMg+2tzTmIsZMpDEXMWYijbmIWSoTWcxJQEREBACnrp7C7P2zsS9jHzRKDd4KfQvzesxDw1oNrV00IiIiIrvAQQIiIrJ7GTczMD92Pr7+42sooMDI9iOxuOdi+NbztXbRiIiIiOyKLAYJlEolvL29oVTK4ukJs2EuYsxEGnMRYybSbC2Xvwr+wnuH3sOGhA3QGXQIDwjH0j5L0b5Re4uVwdYyoZrF9pbGXMSYiTTmIsZMpFUlF0EAjhwB0tMBf39UeeUiW2WpY4WrGxARUY2qiQ77TtEdrDi6AiviVuCu9i66eHXB8r7L0cOnh3kKbSHsCy2DORMRyV92NjBgAJCZCTg5AcXFgK8vsHcv4ONj7dJZn0OubnDq1CnOfFkGcxFjJtKYixgzkVbZXLKzgVatgD59gMmTjT9btTJur4oiXRFWx6+G3xo/LD60GM3rNMeuF3bh6KtHrTZAwGPFsbC9pTEXMWYijbmIMRNplclFEIwDBBkZxsGBggLjz4wM4JlnjO/LgaWOFVkMEhgMBly8eJEzX5bBXMSYiTTmIsZMpFUmF3N22HqDHttPb0fwx8F4e+/bcFY547O/f4Yz489gSKshUFjxXkIeK46F7S2NuYgxE2nMRYyZSKtMLkeOAFlZgE5XertOB1y4YHxfDix1rMhikICIiGyPOTpsQRCw+/xuhGwIwah/j0Le/Tx82PdDpE1Ow2tPvAa1UhZT6xAREVE1pKcDGo30e05Oxvep4nh2RURENaKkwy4qEr9X0mF3717+5+MuxWH2/tk4lH0ILmoXzOo2C7O6zUI913o1V2giIiKyO/7+xrsVpRQXG9+nipPFIIFSqURQUBBnBC2DuYgxE2nMRYyZSKtMLlXtsJOvJ2Pu/rn4MfVHKBVKjH1iLBY8vQDNajerRslrDo8Vx8L2lsZcxJiJNOYixkykVSaXbt2MkxRmZJS+g1GtBlq2NL4vB5Y6Vri6ARER1QhBME5SKNVh+/sDycmlVzm4lHcJCw8uxJbTW2AQDBjaaiiW9F6CYI9gyxfeQtgXWgZzJiKSP6nVDVq2NK5u4O1t7dJZn8OtbqDT6XD06FHoyj746uCYixgzkcZcxJiJtMrkolAYO2Y/P2Nn7eZm/Onvb9xeMkBws/AmZuybgYC1AdicuBk9fHog/rV47Hxhp10MEPBYcSxsb2nMRYyZSGMuYsxEWmVz8fEBzp0D9u8H1q41/kxOltcAgaWOFVk8biAIAq5fvw4r3RRhs5iLGDORxlzEmIm0yuZS0mEfOWKcg8Df33jLn0IB3NPew5pja7Ds8DLkFeXh8UaPY1nfZRjgN8CqqxVUFo8Vx8L2lsZcxJiJNOYixkykVSUXhcI439HD5jyyZ5Y6VmQxSEBERLarbIetM+iw+eRmLPptEa7kX4FvXV+sH7Qe/2z7TygVsrjBjYiIiMhu8WyMiIgsQhAEfJ/8Pdqsb4NxP4+DVq/FmmfWIGVSCka0G8EBAgIAvPXWW2jRogUUCgX++OMP0/a0tDR07doVgYGBCA0NRXJyshVLSUREJF8VOiOz9Q5bpVKhQ4cOUKlUVvl+W8VcxJiJNOYixkykVTWX2MxYdPm8C4Z9NwxX8q9g4dMLkfFWBiZ3ngwnlVMNldYyeKyY17Bhw3D48GH4+PiU2j5u3Di88cYbOH/+PGbOnInXXnvNKuVje0tjLmLMRBpzEWMm0piLmKUyqdDqBocOHULLli3RvXt3/Pzzz2jbti0AoHfv3hg1ahRGjx6NnTt3YsWKFYiLi6vQF3OmYSIi+Tt19RTm7J+DvRl7oVFqML7jeMzrMQ8NazUEYFwBQWq+AkfBvrB8LVq0MJ1z5OTkIDAwELm5uVCr1RAEAU2aNEF8fDxatGjxyH0xZyIi4jlHxfvCCs1J0KNHD9G2nJwcnDx5Evv27QMADB06FJMmTUJWVlaFOmxz0ul0OHToEHr06AG1mtMslGAuYsxEGnMRYybSKppLxs0MzI+dj6//+BoKKPBSu5fwXq/34FvP1/Q7UksV+foaVz4ocxHZpvFYqXmXLl1C06ZNTfkqFAp4e3vj4sWLFT7n0Gq1pv+vVCqhUqmg1+thMBhE23U6XalJoVQqFZRKJXQ6HbRaLY4ePYquXbvC2dkZSqWy1L4BmMpZdvbp8rZrNBoYDAbo9XrTNoVCAbVaXe728spelTpJba9sne7fv2/KRa1Wy6JO1W0nAPjtt99MmcihTuZoJ0EQcPjwYXTt2rXU1VB7rlN126lkxvqSfkQOdXrU9orUqSSXbt26wcXFpVp1ys4G/vY3NTIzFdBoBGi1QIsWwO7dOrRsaT/H3v3793HkyBHT35XKtlNFVflsxpY67OLiYuTn56O4uBhKpdJuGrmm/2N8MBdBEGRRp+q2kyAIpTKRQ53M0U4GgwH5+fnQarWyqVN120mr1SI/Px8GgwGCIMiiTo/aXpE6leSi1Wol63T93nV8cOQDbEjYAJ1Bh2f8nsF7Pd9DSJOQUnUSBKB/fzUuXAB0OgWKi42fz8gQMGAAcPasAJXKPo69kkxK/q7UVIft6MquePGwGyGjoqIQFRVlep2Xl4fo6GjTa29vb4SEhCApKQkXL140bQ8KCkJwcDCOHz+O69evm7Z36NABPj4+OHToEPLz8wEA+/btQ1hYGBo2bIh9+/aVOvZ69eoFV1fXUt8JAOHh4SgsLERsbKxpm1qtxqBBg5Cbm1vqbkx3d3f07t0bly5dQmJiomm7p6cnunbtirS0NKSmppq1TgCqXad9+/bJrk5A1dqpU6dOKCgoMF1Qk0OdzNFO7dq1Q35+Po4cOYKCggJZ1Mlc7aTVamVXJ3O0U0JCArp161blOgkCMGlSb1y75ga9HiguNvYnGRkG9Op1HydOFKJRI/s49hISEkr9XalsO1VUhR43KPHgrX8JCQkYNWoUzp49a3q/U6dOWLFiheSdB1Id9pdffml6XRLIqVOnJAM5evSo5IF74MAByUB2795d7UbOycmRbOTs7GzJRk5JSZFsZNbJdurUqVMn0b7tvU7maKd27drhzJkzcHNzk+yw7bFO5mqnfv36QafTyapO5mgnDw8PdOvWzVSne/p7+DHnR/yU+xMK9YVoW68t/tngn2jr1layTsnJ9fHuu12h04mfqVOr9di16w7+/vd6dnHsHTlyBLm5uVVqp4kTJ/I2+HKUfdwgICAAN27cqPLjBpmZmabX1b0wERMTg379+sHFxcVhB1LLbi8sLDTlotFoZFEnc1yYiI6ONmUihzqZ68LEnj170L9//1J3X9lzncxxYSImJgYDBw6ERqORRZ0etb2iFyZiYmLQv39/uLq6VrlOR44oMGCAyjQ48CAnJwExMQJ69LCPY6+wsBD79u0z/V2pTDv5+vpW+JyjyoME7LBt/z9GdtjssNlhs8N+sE6P2l6VDvte0T18euJTLD2yFLmFuQhuEIwP+nyAvwf8vVQZy9Zp2zYFpkxR4e5dcYddq5aANWsEvPqqfRx7luqwHc2D5xwA0LNnT4wePRqj/zsPUmRkJOLj4yu0L3POSaDVahEdHY3w8HBTP0LMRQozkcZcxJiJNHPlsmULMHky8MB1MBM3N2DtWmD06Crv3qKqk0ll+sIqDxIAttNhGwwG5ObmwsPDA0oll9AqwVzEmIk05iLm6JmUN7lPSS716tfDN2e/wbsH30XW7Sw0c2+GhT0XYnSH0VArH/0k2+HDQJ8+MD1m8CAnJ2D/fqB79xqoWA2ozrHCCfXEJk6ciB9//BHXrl2Dh4cH3NzckJ6ejtTUVIwePRo3btxA7dq1sXXrVrRp06ZC++Q5R81jLmLMRBpzEWMm0ucdgmCeXHjOYWT2QQJb77CJiMh8HjahoLe3gF/Sf8Gc/XOQ9FcS6rrUxZzuczA5dDJcNa4V/g5BAFq1AjIygAcv/qvVxpOD5GTHmHGYfaFlMGciIttV0xMZ85zDqDJ9YYWGHz7++GNcvnwZOp0O165dQ3p6OgDj85hxcXE4f/48Tpw4UeEBAnPTarXYvXs3J4Mqg7mIMRNpzEXMUTMRBGNHnZFh7KQLCow/MzKAp1+KR88tPTHoq0E4f+M8ZnadiQtvXcDMbjMrNUAAGDvjvXsBPz/jCYGbm/Gnv79xuz111o56rDgqtrc05iLGTKQxFzFHzuRh5x0DBgj4+efq58JzjsqTzVpNVZm10REwFzFmIo25iDliJkeOAFlZpUfa4XEOuj7zkN3qB1y6qES/+v2wYcQG+DbwLW83FeLjA5w7J48mn/CdAAAgAElEQVQ1ix3xWHFkbG9pzEWMmUhjLmKOmonkeQeMrzMzgTNnamPAgOp/j72ccwiCgL0Ze/Hd2e+wafAmKBXia/qWOFZkM0hARETVl54OaDRAURGA2peBnguBDl8ASgNU55/DwqcWoV2TDHjV9jLL9ykUxucA7eVZQCIiIjKfUucdZWg0wNWrtcz2XbZ8zlGkK8KOMzsQFReFs9fPQqVQYVLoJIQ0CbFKeThIQEREJv7+QJHyJtBvGRC6FtDcBzJ7Ar8ug+p6Z3Qfq0NeXoa1i0lEREQy4O8vPaEgAGi1QJMmdy1bIAu7ce8GPj3xKdb9vg7XCq7B3ckdEWEReKvzW/Cu4221clVqdQNzMuckQoIgID8/H+7u7lDY2j0jVsRcxJiJNOYi5oiZ3NPew+r4NXhnz3IYnG4D1x4Hfl0GpA+AWq2Avz9w9qyAggLHyuVRqnOscEI9y+A5R81jLmLMRBpzEXPkTB4+qaCA+Ph81K4tv1wybmZgVfwqbE7cjHvae2heuzmmdJ6C1594HXVc6pT7OUudc8jmTgJX18pNmuUomIsYM5HGXMTkkkl5yxmW0Bl02HxqMxb9tghX8q+geX1f6PatQ+7B4XB2UqLYCWjZ8n+T+8glF3NiJo6F7S2NuYgxE2nMRUwumTzqnKOskkkFy65u0LIlsGcP8Nhj8silRNylOETGReKHcz9AgIAnmjyBiLAIPN/6eWhUmgrtwxLHiiwW4tTpdIiOjnbYCT/Kw1zEmIk05iIml0yys40j9H36AJMnG3+2amXcLggCvk/+Hm3Xt8W4n8dBq9dizTNrkP52Cv7c8xIO7Fdi7Vrj+sHJyYC3t3xyMSdm4ljY3tKYixgzkcZcxOSSycPOOR6mZFLB/ftR6ryjaVN55KI36LHr3C50/bwrum7uil3ndiE8IBwHRh3AibEnMKLdiAoPEFjqWJHNnQRERFTag8sK6XT/e+YvIwN4alQsmoycjeNXjsPNyQ0Ln16IaWHT4O7sbvq8rU7uQ0RERLblYecczzxj/Ef/o+4okNt5x93iu9iSuAUr41ci41YGnFROeD3kdUwNm4rWnq2tXbyH4iABEZFMSS4r1DgRur6zccl/L65e1WBy6GS80+MdNKzV0FrFJCIiIjv3sKUML1wwvi+nAYCHuVZwDeuOr8MnJz7BzcKbqO9aH+889Q4mhk5EY7fG1i5ehXCQgIhIpkotK1TvAtBrPtD+K0BQQJ38Et7vuxgzB7a0djGJiIjIzj1sKUMnJ+P7ch8kSL6ejKi4KGxP2o5ifTH86vlhcc/FGN1hNGo5mW8pR0uQzeoGOp0OarVadjNfVgdzEWMm0piLmBwyOXwY6P33HGjD3gM6bgBUWiDtGWD/Ujjd7ID9+yvfYcshF3OrTiZc3cAyeM5R85iLGDORxlzE5JDJ4cPGOQikljN0coJszzkEQcDBrIOIjItEdFo0AKBb826ICIvA4KDBUClVZv8+S5xzyGLiQgAoLCy0dhFsEnMRYybSmIuYPWeSX5SPGO1C6Cb6AZ3XAVdDgC2xwI5foM7tgJYtjTMOV4U951JTmIljYXtLYy5izEQacxGz90y6dQN8fY1LFz5IrYYszzm0ei2+OvMVntz4JHpv64096XswrPUwxL0Wh8OvHsaQVkPMPkBQwhKZyGKQQKfTITY21u5nvjQ35iLGTKQxFzF7zaRIV4Q1x9bAb40fFh9aBN8GzdD0/76HZls83HJ7wsnJuCRRyXKGlWWvudQkZuJY2N7SmIsYM5HGXMTkkEnJUoZ+fsY7B9zcIMtzjrz7eYg8GomWa1ripV0vIfVGKiZ1moS0yWn47vnv0MWrS41+v6Uy4ZwEREQyYBAM+OrMV5gfOx9Zt7PQzL0ZNv19E0Z3GA3V2+pKrVlMREREVFklSxnK8ZzjYt5FrI5fjU0nNyG/OB+N3Rrjg94fYFzHcajvWt/axTM7DhIQEdkhQTB2wmlpAnLr/YIvr81B0l9JqOtSF8v7Lsfk0Mlw1biafl9uywoRERGR7ZHbUoYJVxKwIm4Fvj37LfSCHm0822B61+kY3nY4nNXO1i5ejZHNIIG67AMwBIC5SGEm0piLmK1mkp3937WIi+Jh6DMbhua/QaFzwbgOM7E0fDbqudar0e+31VysiZk4Fra3NOYixkykMRcxW86k5MKENe4OsEYuBsGAX9J+QWRcJA5mHQQA9G3ZF9PDpqO/X3+rT6JoiUxksboBEZGjEASgZWgKsv3nQgj+ATAogcQxUP3fQgQ08kJysjxu63MU7AstgzkTEVVNyYWJzEzjHAPFxcYJCvfuNT5eICf3dffxZdKXWBG3Aim5KVAr1RjedjgiwiLweOPHrV28anO41Q0MBgNycnJgMBisXRSbwlzEmIk05iJmi5lcvnMZf//sdWSFtzEOEJwbAqz/A/jpM+hveeHCBeNIf02yxVysjZk4Fra3NOYixkykMRcxW81EEP5752KGcXCgoMD4MyMDeOYZ4/s1yVK55N7LxXu/vQefVT4Y+5+xuJp/FTO7zkTmlExsG7LNpgYILJWJLAYJ9Ho94uLioNfrrV0Um8JcxJiJNOYiZkuZ3Cq8hVkxsxCwNgC7r3wO5Z9PAZ/FAf/aBeS2Mv2ek5PxVsCaZEu52Apm4ljY3tKYixgzkcZcxGw1kyNHgKwsoOxE+jodLHJhoqZzSbuRhgm7J8B7pTfePfguXNWuWDlgJS5NvYTl/ZbDq7ZXjXxvdVjqWLHdh1+IiBzcPe09rD22FsuOLMPt+7fRvlF7vNx4GeZ+8AwMxeJnCoqLjc8KEhEREVVXejqg0QBFReL3Si5M2NsEhYIg4MilI1gRtwI/pvwIAQI6Nu2I6WHTMbT1UKiV/OcxwEECIiKbozPo8MWpL7Dwt4W4kn8FLeq2wNqBazGi3QgooMRnvsZb/R4c2VergZYtjZMJEREREVWXv7/xAoQUe7swoTPo8MO5H7AibgWO/XkMADA4aDAiwiLwlPdTVp+M0NbIYpBAoVDA3d2djVsGcxFjJtKYi5g1MhEEAT+k/IC5++ci9UYqPB/zxOpnVmPck+NKLbOzd694EqGWLY3ba7q4PFbEmIljYXtLYy5izEQacxGz1Uy6dTNOUmitCxPmyKWguACbT23GqvhVyLydCWeVM8Y9OQ5Tu0xFkEeQGUtrGZY6Vri6ARGRDTiYdRCzf52NY38eg5uTGyLCIhARFgF3Z3fJ37fmckRkPuwLLYM5ExFVjdTqBiUXJry9rV268l3Jv4K1x9bi04RPcfv+bXg85oGJnSZiQqcJaFirobWLZxWV6QtlcSeBwWDApUuX0Lx5cyiVspiL0SyYixgzkcZcxCyVSeK1RMzZPwd70vdAo9RgcuhkvNPjnUd2YAqF8TlASz8LyGNFjJk4Fra3NOYixkykMRcxW87Exwc4d846FyaqksuZv84gKj4KO5J2QGvQIqB+AJb1WYZRj4+Cq8a1hktc8yx1rNjWUVhFer0eiYmJNjcjqLUxFzFmIo25iNV0Jpm3MjFy10iEbAjBnvQ9GNFuBFImpWDNwDU2PcLNY0WMmTgWtrc05iLGTKQxFzFbz6TkwsTo0caflrpzsaK5CIKAmIwYDPhyANp/2h5bEregi1cX/PjPH5EyKQXjOo6TxQABYLljRRZ3EhAR2YucuzlYcmgJPj3xKbQGLZ7xfwZL+yxFh8YdrF00IiIiIrtRrC/G12e+RlR8FJL+SoJSocQLbV5ARFgEQpuFWrt4do2DBEREFpBflI8VcSuwIm4FCooLENosFMv7LkfPFj2tXTQiIiIiu3Gr8BY2JmzEmuNrcCX/CmppamFK5ymY0nkKfOv5Wrt4siCLQQKFQgFPT0+bmxHU2piLGDORxlzEzJVJka4IGxI2YMmhJbh+7zoCGwRiaZ+lGBI8xC7z5rEixkwcC9tbGnMRYybSmIsYM5FWNpes21lYFb8Kn538DHe1d9HUvSmW912ON558A3Vd6lq5tJZhqWOFqxsQEdUAg2DA12e+xvzY+ci8nYmm7k2x8OmFGBMyBmqlLMZnyQzYF1oGcyYisl/H/zyOFXErsDN5JwyCAe0btcf0sOl4se2LcFI5Wbt4dqMyfaFsJi5MSUmx2ck+rIW5iDETacxFrKqZCIKAX9J+wRMbnsDIH0bi1v1bWNZnGdImp2Hsk2PtfoCAx4oYM3EsbG9pzEWMmUhjLmLMRMwgGPDvc/9Gx/Ud0fmzzvj27Lfo17IfYl6OQeK4RLz8+MsOOUBgqWNFFoMEBoMBqampMBgM1i6KTWEuYsxEGnMRq0omxy4fQ+9tvRH+VThSb6RiRtcZyHgrA7O6z8JjmsdqsLSWw2NFjJk4Fra3NOYixkykMRcxZvI/hdpCfHriU7T6uBWGfDsEp3NPY1T7UUh6Mwl7Ru5B35Z9HfqxDEsdK/Z9SYuIyAak5KZg3oF52HVuF5QKJV4LeQ0Ley6EV20vaxeNiIiIyObl3M3Bx8c/xvoT65F7Lxd1XepiZteZCM4Lxsi/jYRGo7F2ER0KBwmIiKrozzt/YuHBhdicuBkGwYAhwUPwfu/30cqzlbWLRkRERGTzUnJTEBUXhW2nt6FIXwTfur54t8e7GBMyBs4KZ0RHR1u7iA5JFoMESqUS3t7eUCpl8fSE2TAXMWYijbmIPSyTW4W3sPzIcqw+thr3dffRw6cHlvddji5eXaxQUsvisSLGTBwL21sacxFjJtKYi5ijZSIIAg5lH0JkXCR+Pv8zAKBzs86Y3nU6hgQPgUqpAmB8/t6RcqkISx0rZlndYO/evZg7dy4MBgO0Wi1mzJiBV1555aGf4UzDRGRvCrWFWHt8LZYeXorb92+jfaP2WNpnKQb6D3To5+Oo6tgXWgZzJiKyPp1Bh53JO7EibgVOXDkBBRR4NvhZRIRFoFvzbjyXqmEWXd1AEASMGDECX3zxBU6dOoWff/4Z48aNQ35+fnV3XWF6vR6nTp3ijKBlMBcxZiKNuYg9mInOoMNnJz9DwNoAzPp1Fuq61MX2IdtxatwphAeEO1SnxmNFjJk4Fra3NOYixkykMRcxuWeSX5SPlXEr4b/GH8O/H46zOWcxvuN4pE5KxQ8v/oDu3t0lz6XknktVWCoTs92ncPv2bQDAnTt30KBBAzg7O5tr149kMBhw8eJFzghaBnMRYybSmIuYwWBAdnY2vk/+Hu0+aYex/xmLIn0RVj+zGikTUzCy/UgoFY53+xuPFTFm4ljY3tKYixgzkcZcxOSayeU7lzEzZiaar2yOafumoVBXiMU9F+Pi1ItYP2g9AhoEPPTzcs2lOiyVSbXnJFAoFPj222/x3HPPoVatWrh16xZ27doFJ6dHr1up1WpN/1+pVEKlUkGv15eqdMl2nU6HB5+MUKlUUCqV0Ol0pv1otVrT9gf3DQBqtbGqOp2uQts1Gg0MBkOpURqFQgG1Wl3u9vLKXpU6SW2vbJ0ezEUudapuO5Uw17FnC3UyRzuV/E7ZsthznarbTgcuHMCstFk4f/o8amlqYV73eZjaeSpqO9cGhP/lZk91etT2irRTyXfodDpoNBpZ1Kns9qrUCfjf35XK1okqpyqPOBIRkeUkXkvEirgV+OaPb6Az6NDKoxUi+0diZPuRcFG7WLt4VAHVHiTQ6XRYunQpfvzxR3Tr1g2///47/vGPf+DMmTOoX7++6feioqIQFRVlep2Xl1dqtkpvb2+EhIQgKSkJFy9eNG0PCgpCcHAwjh8/juvXr5u2d+jQAT4+Pjh06JDp0YaYmBiEhYWhYcOG2LdvX6mTxV69esHV1VU0Q2Z4eDgKCwsRGxv7v1DUagwaNAi5ubmIi4szbXd3d0fv3r1x6dIlJCYmmrZ7enqia9euSEtLQ2pqqlnrBKDadYqJiZFdnYCqtVOnTp1MmcilTuZop3bt2gEAjh49ioKCAlnUqart9PH3H2P7le04mX8SaoUab4a8iYhOETh7/CwO7z9sl3WqiXY6efIkunXrJqs6VbWdTp48CeB/f1cqWyequJJHHGNjY9G+fXtkZWUhODgYzz33HNzd3a1dPCIihyUIAvak78GKuBXYn7kfANCrRS9EhEVgYMBAh7z70p5Ve+LCEydOYNSoUUhOTjZt69SpEz788EP06tWr3M95eXkhMzPT9Lo6V590Oh0yMjLg5+cHJycnu7r6VJNX1IqKiky5qFQqWdSpuu2kUCiQmpqKli1bQqVSyaJO5mgnAMjIyICvr2+p2VLtuU6VbadL+Zfw7sF38dWZryBAwAutXsBon9Ho92S/UleL7alONdFOer0eGRkZ8Pf3h7OzsyzqVHZ7ZetUVFSE9PR009/aytTJ19eXE+pVgiAI8PDwwA8//IAePXogKSkJAwcORGZm5kPvYDTnxIV6vR5paWkICAgw9SPEXKQwE2nMRcyeMynSFWHHmR2IiovC2etnoVKo8GLbFxERFoEnmjxRrX3bcy41pTqZVKYvrPYgwV9//YWAgAD8/vvvCAoKQnp6OkJDQ3HmzBk0a9bMLIUkIqopOXdzsOTQEnx64lNoDVoM8BuApX2WIqRJiLWLRg6AfWHl7d+/Hy+++GKpRxz79u370M+Y88KEvQ5msU6sE+vEOpmzTtfyrmHjqY345MQnuHb3Gtyd3DH2ibGY8OQEeNfxtss6ybGdqnphotqPGzRq1AgbNmzAsGHDoFQqIQgC1q9f/9ABAnPT6XQ4fvw4QkNDTY1GzEUKM5HmiLnkF+UjKi4KkXGRKCguQGizUCzrswy9fI13QDliJhXBXMSYieXY0iOOgH09FsM6Wb5OoaGhiI2Nxb1792RTJ3M94njt2jUUFhY6/COOZevUv39/aLVam65Tk9ZNsD1tOzYlbEKRoQgeGg+Mbjoai/+xGI3rNkZ0dDT+wB9mbScPDw8+4vjfOh07dgy5ublVrlNFVftOgqoy59UTrVaL6OhohIeHQ6PRmGWfcsBcxJiJNEfKpVhfjA0nNuC9Q+/h+r3rCGwQiA96f4DnWj0nmtzSUTKpDOYiVp1MeCdB5djCI47FxcWIiYlBv3794OLi4hBXnypSp8LCQlMuGo1GFnWqbjsJgoDo6GhTJnKok7kmS96zZw/69+9famDVnutU3XbSarWIiYnBwIEDodFobLJO8ZfjsfLYSvw79d8QIKBDow6Y2mUqhgUPg0alqbHJkmNiYtC/f3+4urpavZ3MUaey2ytbp8LCQuzbt8/0d8Vm7yQgIrIHBsGAb/74BvNj5+PCrQto4tYEG/62Aa+GvAq1kn8KiexB8+bNcfnyZaSmppoecczIyEBgYOAjPys1gFMyh0RZ5d0RolarTSdhGo3GNIdLeYNDldmuVCpLzQnzqO3llb0qdapoGSuyvWSA4GFlt7c6PagydXpwdamy+7LXOgHVb6eSXNRqtWR57LFOj9pe0bKXzJ9lK3XSG/T4Ke0nRB6NRNxl49Xxgf4DMaPrDPRs0bPUxZXy6vSw7RWtU0nZbKWdHrbdUu1U9u9KZev0KDwzJiJZEwQBezP2Ys7+OUi8log6znWwtM9SvNX5LTymeczaxSOiSrCFRxyJiOTubvFdbEncgpXxK5FxKwNOKie8FvIapoVNQ2vP1tYuHlmALAYJVCoVOnTowFkvy2AuYsxEmlxzOXb5GGbvn42DWQfhrHLGjK4zMLv7bNR3rf/Iz8o1k+piLmLMxLKGDx+O4cOHW+372d7SmIsYM5HGXMRsJZNrBdew7vg6fHLiE9wsvIn6rvXxzlPvYGLoRDR2a2zx8thKLrbEUpnIYk4CIqIHpeSmYN6Bedh1bheUCiVGPz4aC3suRPM6za1dNKJS2BdaBnMmIirf2ZyziIqLwpdnvkSxvhj+9f0xrcs0vNLhFd51KSOV6QvFD0bYIZ1OhwMHDlRp5kY5Yy5izESaXHL5886fGPvTWLRd3xa7zu3CP4L/gTPjz+DzZz+v9ACBXDIxN+YixkwcC9tbGnMRYybSmIuYNTIRBAEHMg8gfEc42n7SFpsTN6NT007Y9cIupExMwfhO460+QMBjRcxSmcjicQNBEJCfnw8r3RRhs5iLGDORZu+53Cq8heVHlmP1sdW4r7uPHj49sKzPMoQ1D6vyPu09k5rCXMSYiWNhe0tjLmLMRBpzEbNkJlq9Ft+e/RaRcZFIvJYIpUKJYa2HISIsAl28utT491cGjxUxS2Uii0ECInJMhdpCrD2+FksPL8Xt+7fRrmE7LOu7DAP9B0rOuEtERETkiPLu52FjwkasOb4Gl+9cxmOaxzA5dDLe7vI2WtZrae3ikY3hIAER2R2dQYetiVux4OAC/Jn/J3zq+GDNM2swot0IqJSc3IaIiIgIALJvZ2P1sdX47ORnyC/OR2O3xni/9/t4s+ObFZrImRyTLCYuNBgMyM3NhYeHh+T6k46KuYgxE2n2kosgCPh3yr8x98BcpOSmwOMxD7zz1Dt4s+ObcFY7m/W77CUTS2MuYtXJhBPqWQbPOWoecxFjJtKYi1hNZJJwJQGRcZH47ux30At6tG3YFhFhERjedrjZz5lqCo8VMUudc8hikICI5O9Q9iHM+nUW4i/Ho5amFqaFTcP0rtNR27m2tYtGVGXsCy2DORORIzAIBkSnRSPyaCR+y/4NANC3ZV9MD5uO/n79+Simg3O41Q20Wi12794NrVZr7aLYFOYixkyk2XIuSX8lYdBXg/D0lqdx4soJTOw0ERlvZWBxr8U1OkBgy5lYE3MRYyaOhe0tjbmIMRNpzEWsupnc193HpoRNaLO+Df7+9d9x5NIRvNz+ZSSOS0TMyzEY4D/ALgcIeKyIWSoT2cxJwKUxpDEXMWYizdZyybyViXcPvosdSTsgQMDwtsPxXq/34Fffz2JlsLVMbAVzEWMmjoXtLY25iDETacxFrCqZ5N7LxSe/f4J1v69Dzt0c1HGug5ldZ2Jy58nwqu1VA6W0PB4rYpbIRDaDBEQkD9fvXseSQ0vwyYlPoDVoMcBvAJb2WYqQJiHWLhoRERGR1Z2/cR4r41Zi6+mtKNQVwqeOD1YOWInXQl6Du7O7tYtHMsBBAiKyCflF+YiKi0JkXCQKigvQqWknLOu7DL19e1u7aERERERWJQgCjlw6gsijkfgp9ScIENCxaUdEhEVgWOthUCv5zzoyH1lMXCgIAvLz8+Hu7m6Xz9vUFOYixkykWTOXYn0xNiZsxOLfFuP6vesIqB+AD/p8gKGthlq1jXisSGMuYtXJhBPqWQbPOWoecxFjJtKYi9jDMtEZdPjh3A+IjIvE8T+PAwAGBw1GRFgEnvJ+StYZ8lgRs9Q5h2yGnFxdXa1dBJvEXMSYiTRL52IQDPjmj28wP3Y+Lty6gCZuTbDhbxswpsMYaFQai5alPDxWpDEXMWbiWNje0piLGDORxlzEymZSUFyAzac2Y2X8SmTdzoKL2gXjnhyHqV2mIsgjyEqltDweK2KWyEQWqxvodDpER0dzYosymIsYM5FmyVwEQcDe9L14cuOTeGnXS7hx7wY+6P0B0t9KxxtPvmEzAwQ8VqQxFzFm4ljY3tKYixgzkcZcxB7M5Er+Fcz5dQ6ar2yOKXumoKC4AAueXoCLb1/Ep3/71KEGCHisiFkqE9ncSUBEtu/4n8cx+9fZiM2KhbPKGdPDpmN299lo8FgDaxeNiIiIyGqyCrPw2n9ewzdnv4HWoEVgg0As77scL7d/Ga4aXk0ny+IgARHVuNTcVMw7MA/fn/seSoUSYzqMwaKei9C8TnNrF42IiIjIKgRBQMyFGEQeiURMZgwAoIdPD0wPm45BgYOgVMjipm+yQxwkIKIa8+edP7Hot0XYfGoz9IIezwY9iw/6fIDWnq2tXTQiIiIiqyjWF+PrM19jRdwKnMk5A5VChe51u+PDf3yIMJ8waxePSD6rG+h0OqjVas58+QDmIsZMpJk7l9v3b2P54eVYfWw1CnWFeMr7KSzruwxdm3c1Q2ktg8eKNOYiVp1MuLqBZfCco+YxFzFmIs2Rc7lVeAsbEjZg7fG1uJJ/BW5Obng95HW81fkteLl5OWQmD+PIx0p5LHXOIZs7CQoLC+Hu7m7tYtgc5iLGTKSZI5dCbSHWHV+HpYeX4tb9W2jXsB2W9lmK8IBwu/zjzmNFGnMRYyaOhe0tjbmIMRNpjpZL5q1MrIpfhc9PfY672rto5t4My/suxxtPvoG6LnVLLWtHpTnasVIRlshEFg+66HQ6xMbGcubLMpiLGDORVt1cdAYdPj/5OQLXBWLmrzNR27k2tv1jG06NO4VBgYPscoCAx4o05iLGTBwL21sacxFjJtIcKZdjl4/hhe9egP9af6w5vgb+9f2x7R/bcGHKBczsNhN1XeoCcKxMKoO5iFkqE9ncSUBElicIAv6d8m/MOzAP53LPoYFrA6wasApvdnwTzmpnaxePiIiIyKIMggH/Sf0PIuMicfjiYQDAM/7PYHrYdPT27W2XF07I8XCQgIiq5FD2Icz+dTbiLsehlqYW5veYj+ldp6O2c21rF42IiIjIou5p72Hb6W2IiotC2s00aJQajOkwBtPCpqFtw7bWLh5RpchmkECtlk1VzIq5iDETaRXNJemvJMzZPwfRadFQK9WY0HEC5j89H43dGtdwCS2Px4o05iLGTBwL21sacxFjJtLklEvO3Rx8fPxjfPz7x7hReAP1XOphbve5mBQ6CU3cm1R4P3LKxJyYi5glMpHF6gZEVPOybmdhfux87EjaAQEC/tn2n3iv13vwr+9v7aIR2S32hZbBnInI3FJzUxEVF4Wtp7eiSF8E37q+mNplKsaEjIGbk5u1i0ckUpm+UBYTFxoMBuTk5MBgMFi7KDaFuYgxE2kPy+X63et4e8/bCFoXhC+TvkQ/v35IeCMBXw/9WtYDBDxWpDEXMWbiWNje0piLGDORZs+5CIKA37J+w+CvByP442BsPLkRIU1C8N3z3yFtchomd55cpRjL6EcAACAASURBVAECe86kJjEXMUtlIotBAr1ej7i4OOj1emsXxaYwFzFmIk0ql4LiAiz+bTH81vhh9bHVaN+oPfaP2o+9I/fiiSZPWLG0lsFjRRpzEWMmjoXtLY25iDETafaYi86gwzd/fINOmzqh59ae+Pn8zxgSPARHXj2CuNfiMKz1MKiUqirv3x4zsQTmImapTPiQBxGVUqwvxsaEjXjv0HvIuZuDgPoB+Lz35xjWehhn5CUiIiKHcafoDj4/+TlWHVuFi3kX4ap2xYSOEzA1bKqs76Yk4iABEQEwLtnzzdlvsPDQQly4dQFN3Jrg00Gf4tWQV6FRaaxdPCIiIiKLuHznMtYcW4MNCRtwp+gOGtVqhPd6vYfxHcejwWMNrF08ohoni0EChUIBd3d3XuUsg7mIMRMxQRAQcyEGM9NnIv10Omo718b7vd/HlM5TUMuplrWLZzU8VqQxFzFm4ljY3tKYixgzkWbLuSReS8SKuBX45o9voDPo0MqjFaL6R+Gl9i/BRe1SY99ry5lYE3MRs1QmXN2AyIEd//M4Zv86G7FZsXBWOWNS6CTM6T6Ho+REFsK+0DKYMxGVRxAE7M3Yi8ijkdifuR8A0KtFL0SERWBgwEAoFbKYwo3IMVc3yM7O5syXZTAXMWZidP7GeTz/3fPo/Fln/Jb9G0Y/PhoHhhzAh30/5ADBf/FYkcZcxJiJY2F7S2MuYsxEmq3kUqQrwhenvkC7T9ph4I6BOJh1ECPajUDCGwk48MoBDAocZLEBAlvJxNYwFzFLZWKWI7+oqAiTJk1CQEAA2rRpg5EjR5pjtxWm1+uRmJjImS/LYC5ijp7JlfwrGPefcWj9cWvsTN6JZ4OeRdKbSdg4aCOup1932FykOPqxUh7mIsZMHAvbWxpzEWMm0qydy417N/D+offRYnULvPrTq7iYdxERYRG4MOUCdjy3wyorOFk7E1vFXMQslYlZ5iSYPXs2lEolzp8/D4VCgatXr5pjt0RkJrfv38byw8ux+thqFOoK0d27O5b1WYZu3t0AAFqt1solJCIiIqo5GTczsDJ+Jb5I/AL3tPfgVdsLH/X7CGOfGIs6LnWsXTwim1LtQYK7d+/iiy++wOXLl00TKDRp0qTaBSOi6ivUFmLd8XVYengpbt2/hbYN22Jpn6UYFDCIk8AQERGR7MVdikNkXCR+OPcDBAh4oskTiAiLwPOtn+fqTUTlqPYgQUZGBho0aIAlS5bg119/haurKxYuXIg+ffo88rMPXr1UKpVQqVTQ6/WlnrEo2a7T6fDgHIsqlQpKpRI6nQ46nQ4NGjSATqczbS97ZVStNlZVp9NVaLtGo4HBYCh1K4dCoYBarS53e3llr0qdpLZXpU4lucipTg+qbJ0UCgU8PDxK7cfe6yTVTgYYsOOPHVgQuwCX8y/Du7Y3IvtFYmS7kXDSOInKLggCPD09Rbcu2VKdLN1OJf/9lOQjhzo9antF6lSSi16vh0ajkUWdym6vbJ30en2pv7WVrRPZF4VCAU9PTw60lsFcxJiJNEvkojfo8WPqj4g8Gom4y3EAgEEBgxARFoGeLXraXJvwWJHGXMQslUm1VzdISEhAx44dsXXrVowaNQqnT59G3759kZycDE9PT9PvRUVFISoqyvQ6Ly8PX375pem1t7c3QkJCcOrUKVy8eNG0PSgoCMHBwTh69CiuX79u2t6hQwf4+PjgwIEDyM/PN20PCwtDw4YNsXv37lIni7169YKrqyuio6NLlT88PByFhYWIjY01bVOr1Rg0aBBycnIQFxdn2u7u7o7evXsjOzsbiYmJpu2enp7o2rUrUlJSkJqayjqxTlarkyAIOJZ3DP/K/RcyCzJR16kunvN4DgM9BkKj1NhlneTYTqwT61RSp4kTJ3LWfQvg6gZEjuFu8V1sSdyClfErkXErA04qJ7zc/mVMC5uG1p6trV08IquqTF9Y7UGC3NxcNGrUCMXFxVCpVACA0NBQfPjhh+jZs+dDC5mZmWl6Xd07CTIyMuDn5wcnJye7ufpU01fUioqKTLmoVCpZ1MkcdxKkpqaiZcuWpuPV3utU0k6Hsg5hbuxcxP8Zj8c0j2Fal2mY2nkq3J3cH1knwHhXkK+vL5TK/81nau06WfPY0+v1yMjIQFBQkGn/9l6nR22vSJ1KcvH394ezs7Ms6lR2e2XrVFRUhPT0dNPf2srUydfXl/94raSioiJERERg7969cHJyQkhISKmLDlLMOUig1+uRlpaGgIAAUz9CzEUKM5FWE7lcK7iGdcfX4ZMTn+Bm4U3Ud62PiZ0mYmKniWjk1sgs31GTeKxIYy5i1cmkMn1htR838PDwQJ8+fbB3716Eh4cjOzsbmZmZCAoKeuRnNRrxc0AlJ1iigqqli6pWqyEIAtLT0xEYGGj6B47Uviu7XalUlvoH06O2l1f2qtSpomV82HaVSmXKpeR37L1O1W0nrVZr+g+r7L7stU5n/jqDOfvnYHfabqiVakzoOAHzn56Pxm6NJb9PquxarRapqanw8/OTLI+j/vdU8t+PQqGQTZ0etb0idSrJ5WFlt7c6PaiydVIqlaK/tUDl60QVY+3Jkg0Gg+nvJU9a/4e5/H97dx4QVbm/AfyZDUHBXMAtZJNEzQVNS0ArNXDLVrt509SWa6UtFt70mmmaEd0ETU2trK6ZP29l3W4qKYimKbiiouLCJmCu4AYKzHZ+f8xlAs9BtmHOzJnn849yYIbv+7wH5uWdc95XjJlIs2Uuxy4eQ3xqPL498i30Jj2CWwXj/UHvY0KvCWjm1sxGFTc+nivSmIuYvTKxye4GK1aswPPPP4/p06dDo9Hg888/5+KFRI3s9NXTmL1tNr5N/xYCBIzpPgbvD3ofwa2C5S6NiKhRcLFkIhIEAVtztyIuNQ6/Zv0KAAjvGI5pYdPwSMgj0Kj5xyRRQ9lkkiAoKAi//fabLZ6KiGpQeLMQH+z4AMv2L4PepEdkUCQ+HPIh7ulwj9ylERE1KkdYLLnieQwGg1PdFtPYt/pUzkUpbWpoP1Ww1bnnCG2yRT9VfM2ttdTUptLyUvxw/Acs3LMQhy8chlqlxuiuo/HGvW/gvjvvszy3yQyNWuN0517F5wVB4GLJldpU8T2MRiMXS67UJuDPc6axFku2ySSB3NRqNfz8/CQvG3VlzEXMmTMp0ZdgYepCfJzyMYr1xbin/T346KGPMCSo5sFxTZw5l8bCTKQxFzFmYj8GgwE5OTno1q0bYmNj67RYcuXFLisWo0xPT5dcjHLv3r2Si1Hu2LHDuhhlUlKSdTHKxMTEBi+wWVhYKLnAZkFBgeQCm5mZmZILbDakTQAa3KakpCTFtQmoXz/dd999aN68OZKSkhTTJlv0U8+ePeHn54eUlBSUlJTU2KY+4X2w5sQaLPh9AYoMRWiiboKR3iMR91Qc2ru3x7Zt25BwOEHWNtmqn8xmM4qLix2inxzp3Dt48KDi2lTffjp48CAAWH+v1LVNtdXghQvriysNE9WO3qTHFwe+wPs73seFGxcQ3CoYHwz+AKO7jYZaxT9KiJwZXwvrxhEWS3bmd5/YJrbJmdqUdy0PS/ctxVeHvkKxvhjtmrXD5L6TManPJLTyaOWUbaqpP9gmtqkx21SXxZIVMUlgMpmQnp6Onj17clGLSpiLmDNlYhbM+P7Y93hn6zvIuZKDdp7tMOeBOXih9wvQaWy78Jkz5WIvzEQacxFrSCacJKi7qKgoTJ061bpYct++fZGenn7btQk45mh8zEWMmUirKZcDZw8gLjUO3x/7HibBhO5tuiM6LBp/7f5XNNE2kaHixsdzRRpzEbPXmEMRb0OazWbk5+dXmW0h5iLFGTIRBAGJ2Yno+3lf/PXHv1rWIBj8AbJey8LLfV+2+QQB4By52BszkcZcxJiJfa1YsQL//Oc/0aNHDzz66KN2XyyZ/S2NuYgxE2lSuZgFMzac2oBBqwah7xd9sfboWgwKHIRNYzch/eV0TAydqNgJAoDnSnWYi5i9MlHEmgRESrHvj32YkTwDW3O3oommCaLDovGPAf9A66at5S6NiMghcLFkIuUoM5Zh9eHViN8djxOFJ6BVa/Fsz2cRHRaNXu16yV0ekcviJAGRAzhVdArvbH0H6zLWQa1SY2LoRMx9cC787vCTuzQiIiIim7puvI75v8/HirQVuHjjIu5ocgfeDn8br933Gnyb+8pdHpHLU8QkgVqtRkhICFeWvgVzEXO0TM4Wn8Xc3+biy4NfwiSY8EjII4gZHIO729xt1zocLRdHwEykMRcxZuJa2N/SmIsYMxE7VXQK8Snx+Nfxf6HcVA7/O/yxcOhCvND7BXg18ZK7PNnwXJHGXMTslYkiFi4kcjZXy67in7v+iUW7F6HUWIoBfgMQOyQWEX4RtX4OQQB27QKysoDgYCAiAqi0LTMROQG+FtoHcyaSjyAI2FWwCwtSFuCXk79AgIC+Hfri7+F/xxNdn4BWrYj3LIkcnsstXGg0GpGSklKvPSCVjLmIyZ1JmbEMC1IWoNPiTvhw54fo1KoT1v91PXZM3FGnCYK8PKBrV2DIEOC11yz/du1qOV4fcufiiJiJNOYixkxcC/tbGnMRc/VMjGYjfjj2A/p/2R8Dvx6I/578L0aFjMLWZ7diUbdFeCKEEwQVXP1cqQ5zEbNXJor4yRQEAZcuXYJMF0U4LOYiJlcmJrMJ3xz+BnN+m4OC6wXwu8MP8VHxGNdzHDTqum1fIgjA0KFAdjZgNAJ6veV4djYwbBiQkVH3Kwp4rogxE2nMRYyZuBb2tzTmIuaqmZToS/DVwa+waPci5F7NhbvWHS/d8xLe7P8mQrxDYDAYkHAkweVyuR1XPVdqwlzE7JWJIiYJiByVIAj45eQvmLl1JjIuZaC1R2vER8XjlX6vwF3rXq/n3LULOH3aMkFQmdEI5ORYPj9gQMNrJyIiIqqts8VnsWTPEqw4sAJXy67Cu6k33nvgPUzuNxk+zXzkLo+I6oCTBESN5Pe83zEjeQZSClLQVNcUswbOwrTwabjD/Y4GPW9WFqDTAeXl4s+5uVk+z0kCIiIisof0C+mIT43H/x35PxjMBnRu3RmxQ2Ixvtd4eOg85C6PiOpBEZMEGo0GoaGh0Gjqdtm20jEXMXtkcuTCEfwj+R/YmLkRWrUWr/R9Be/e/y7ae7W3yfMHB/95i8Gt9HrL5+uK54oYM5HGXMSYiWthf0tjLmJKzkQQBCTlJCEuNQ6J2YkAgPv978e0sGkY2Xkk1Krqlz1Tci71xUykMRcxe2XC3Q2IbCTvah5m/zYbqw+vhgABT9/9NOYPno/gVvX4q/02BMGySGHFmgQVtFrLBEF91iQgInnwtdA+mDNRw1TsqHQiU49cz7VYXxiHIxePQKPSYHS30YgOi0a/O/vJXSYR3YZL7m6wdetWrnx5C+Yi1hiZFN4sxJub3kTnpZ3xzeFv8FDQQ9j/t/349+h/23yCALBMAGzeDHTqZLm9wNPT8m9wsOV4fSYIeK6IMRNpzEWMmbgW9rc05iKmlEzy8oDOPa/gwVmxmHQsEDEZE3H0j1w8320qsl7Pwr9H/7tOEwRKycWWmIk05iJmr0wUcbuBIAgoLi7mype3YC5itsykRF+ChakL8XHKxyjWF+Oe9vcg9qFYPBT0kA0qvT1/f+D4ccusflaWZYIgIqL+VxDwXBFjJtKYixgzcS3sb2nMRUwJmeRczsU9/1iEq498CbjdAK7fCSR9BPWhSUjp2AL+o+v+nErIxdaYiTTmImavTBQxSUBkTwaTAV+kfYF52+fhwo0LCG4VjJWDV2J0t9G3vQfP1lQqywKFXKSQiIiIbGnPmT2IS43Djxk/whxiBs73AlKigWNPAyY3mMAdlYiUjJMERLVkFsz4/tj3mLV1FrKvZKOdZzssH7kcL/R+ATqNTu7yiIiIiOrNLJix/uR6LEhdgJ35OwEA3d2HImt1NMqOPQSg6uWK3FGJSLkUsXCh2WxGYWEhvL29oVYrYpkFm2AuYvXNJCk7CTOSZyDtXBqaN2mOt8PfxtT+U9HMrVkjVms/PFfEmIk05iLWkEy4oJ59cMzR+JiLmLNkctNwE98c/gbxqfHIvJwJnVqHcT3H4a2wt3D1VHcMGSK9q5KbG5CcXPdJAmfJxZ6YiTTmImavMYciJgmIGsv+s/sxY8sMJOcmw03jhlf7vYqZA2eiddPWcpdGRArA10L7YM5EYhdvXMSnez/Fp/s+RVFpEVq6t8QrfV/Bq/e+at22mTsqESmHy+1uYDAYsHHjRhgMBrlLcSjMRay2mZwqOoW//PAX9PuiH7bmbsWEXhNw6tVTiBsap8gJAp4rYsxEGnMRYyauhf0tjbmIOWomJwpPYNL6SfBb6Id5O+aheZPmWDxsMQreLMAHQz6wThAAjbOjkqPmIidmIo25iNkrE8WsScCtMaQxF7HbZXKu+Bzmbp+LlWkrYRJMGNV5FGKGxKB7m+52rFAePFfEmIk05iLGTFwL+1sacxFzlEwEQcCOvB1YkLoAG05tAAD09+2PaWHT8FiXx6BRa6p9rK13VAIcJxdHwkykMRcxe2SimEkCooa4WnYVH+/6GAt3L0SpsRQRHSMQ+1AsBvhxNR4iIiJyTgaTAesy1iEuNQ4Hzh2ACio83uVxTAufhvCO4bV+Hu6oRORaOElALq3MWIZP936KmJ0xuFx6GXf73I0Ph3yIhzs/DBVvsiMiIiIndL38OlamrcQnez5B/rV8eGg98ErfV/Bm/zdxV+u75C6PiBycIhYuFAQBxcXF8PLy4h92lTAXsYpMmjZritXpqzHntzkouF4Avzv8MO/BeRjXc9xtL7lTKp4rYsxEGnMRa0gmXFDPPjjmaHzMRUyOTAquFWDxnsX4PO1zXC+/jjbN2uC1e1/Dy31fhndTb7vUUBOeK2LMRBpzEbPXmEMxVxJ4eHjIXYJDYi5VCYKALQVbMHv7bBy7dAytPFohLioOk/tNhrvWXe7yZMVzRYyZSGMuYszEtbC/pTEXMXtlcvDcQcSlxuG7Y9/BaDaiq3dXxEfFY2zPsQ45vuG5IsZMpDEXMXtkoojdDYxGIxISEriwxS2YS1U783di4NcD8eS6J5F7NRfvDHwHOa/n4K2wtxzyBdSeeK6IMRNpzEWMmbgW9rc05iLW2JkIgoBfM3/FkG+GoM/nfbDmyBoM9BuIjc9sxNHJR/FCnxcccnzDc0WMmUhjLmL2ykQxVxIQVefoxaP4R/I/sOHUBmhUGgxrPQyfjf0Mfi395C6NiIiIqE7KjeVYc2QN4lLjkHEpAxqVBmN7jEV0WDR6t+8td3lEpACcJCDFyruah9m/zcbqw6shQMBf7v4L5gycg8zdmWjv2b7mJyAiIiJyEEU3i7Bi/wos2bsEF25cgJebF6aFTcPr972Ojnd0lLs8IlIQThKQ4hTeLETM7zH4dN+n0Jv0eCjoIcQOicU9He6BwWBAJjLlLpGIiIioVrIvZ2Ph7oX4+tDXuGm4iY7NOyIuKg4v9nkRzZs0l7s8IlIgxexuYDQaodVqufJlJa6Wyw39DSzcvRAfp3yM6+XX0ad9H4xrG4uWVyIRHAxERACAa2VSW652rtQGM5HGXMQakgl3N7APjjkaH3MRa2gmKQUpiEuNw3+O/wcCBPRp3wfRYdF4qttT0Gl0jVCxffBcEWMm0piLmL3GHIq5kqC0tBReXl5yl+FwXCEXg8mAlWkrMW/HPJwvOY/gVsH4oP/nWPLyU5iRq4abG6DXA4GBwKZNQKtWys+kPlzhXKkrZiKNuYgxE9fC/pbGXMTqmonJbMLPJ35GXGocUs+kAgBG3jUS08Kn4QH/BxTzhxLPFTFmIo25iNkjE8XsbrBt2zaufHkLpediFsz47uh36LasGyYnTIYgCFg2YhmOvZKBpS8/jZxsNfR6oKTEMkmQnQ0MGwZs3arcTOpL6edKfTATacxFjJm4Fva3NOYiVpdMbuhvYOnepei8tDNG/zAaaefS8GLvF5ExOQMbntmABwMeVMwEAc8VMWYijbmI2SsTxVxJQK5lS84WTN8yHWnn0tC8SXPMHzQfU/tPRTO3Zti5Ezh9Grj1Z8doBHJzgePHW2HkSFnKJiIiIrI6V3wOS/cuxfL9y3Gl7Apae7TGu/e/iyn9pqCtZ1u5yyMiF2WzSYK5c+fivffew5EjR9C9e3dbPS1RFQfOHsCM5BnYkrMFbho3vNn/TcwcOBPeTb2tX5OVBeh0QHm5+PE6HXDuXDM7VkxERERU1bGLxxCfGo9vj3wLvUlvuVVy8AeYEDoBTXVN5S6PiFycTSYJ0tLSsHv3bvj5ybfvvFbLiyKkKCWXzKJMzNo2C98f+x4qqDCh1wTMfXAu/Fv4i742ONhye4EUgwHw9S1r5Gqdk1LOFVtiJtKYixgzcS3sb2nMRaxyJoIgYGvuVsSlxuHXrF8BABEdIzAtfBpGdR4FjVojV5l2x3NFjJlIYy5i9sikwbsblJeX48EHH8T//d//YdCgQdiwYUOtriTgis5UG+eKz2He9nlYeXAljGYjRnUehZghMejepvpzTBCArl0taxBUvuVAq7VMIGRkAAq5rY+InBxfC+2DOZOcDCYDvj/2PRakLsCh84egVqnxRNcnEB0Wjf6+/eUuj4hchF13N5g9ezbGjRuHwMDAOj/WYDBY/69Wq6HRaGAymWA2m0XHjUYjKs9naDQaqNVqGI1GmEwmFBUVoXXr1tDpdFCr1VWeG/hzxuXWRR6qO67T6WA2m2EymazHVCoVtFpttcerq70+bZI6Xtc26fV6ay5qtdqp2lRUUoQFuxdgyb4luGm4iXDfcMQ+FIv+HSwvphWPq65NmzdrMXSogNxcyy0GBoNld4NffxVw4cJFtGzZEmq12iH6yRHOPZVKhcuXL6NFixZVFkZy5jY1tJ/MZjOKiorQtm1b6/M4e5tqOl6bNlXk4u3tDTc3N0W06dbjdW2TXq9HYWGh9XdtXdtEzsVsNqOwsBDe3t7W1xFiLlKu3LyCT3Z+gi+PfYkz18+gma4ZXrv3NUztPxVBLYPkLk82PFfEmIk05iJmr0waNEmQmpqKffv2ITY2tsavjY+PR3x8vPXja9euISEhwfqxn58fevfujfT0dOTn51uPh4SEoEuXLti7dy8uXbpkPR4aGgp/f3/s2LEDxcXF1uNhYWFo06YNEhMTqwwWBw0aBA8PjyrfEwBGjBiB0tJSbNu2zXpMq9Vi5MiRKCwsRGpqqvW4l5cXBg8ejIKCAhw6dMh63MfHB+Hh4cjMzMTJkyfZpga0KWxgGFZlrML7v72PYlMxOrp3xFTfqZgzZg7KysqqtLWmNm3alI+1a/Nx7lwztG9/A/ffr8Gdd/ZDQsKeKnmxn4AePXrgyJEj8PT0RElJiSLaZKt+ioyMtK4kq5Q22aKfvL29ERERoag21bef9u3bh8LCwnq3iZyLyWRCamoqRowYwUFrJa6eiyAAu3ZZ1kXyvDMPu0yf4MuDK1GsL0Z7z/aIGRyDl/q+hFYereQuVXaufq5IYSbSmIuYvTJp0O0GsbGxWLx4Mdzc3AAAZ86cQdu2bbFy5UoMHz78to/19fVFbm6u9eOGvPuk1+uRlJSEyMhIuLu7O827T439jlppaak1F51O59BtMplNWHN0DebumIuC6wXo2Lwj5tw/B2O7j4VGrbFZPwmCgISEBGsmjtBPjnDumc1mbNq0CVFRUVXuc3LmNjW0nwwGA5KSkjB8+HDodDpFtKmm47VpU0UuUVFR8PDwUESbbj1e1zaVlpYiMTHR+nulLm0KDAzkZfB2YMvbDQwGAxISEjBixAjr6wi5di55ecDQoUB26X4gLA7GkB8AtQmd77gbQ5sPRsxfY+Dp4Sl3mQ7Dlc+V6jATacxFrCGZ2O12gxkzZmDGjBnWjwMCAmq9JgEAyYZpNBpoNOKFW6pboEGr1VoHYRW3GlT33HU9rlarJWdoqjteXe31aVNta6zN8YoJgtvVLlebtFot1p9aj5nJM3Hs0jG08miFuKg4TO43Ge5a92rbVJvapY5Xvk3h1ueSu58qs3c/VeSi1Wol63HGNtV0vLa1q1QqqFQqRbXpdsdr26aKGpTUpgr1bdOtv1fq2iYicj4msxkRzyXgj/AFgP92y8HsSKj3RAOqSAyOXY8m2ibyFklEVEeKWC5SpVLBy8uryr3U5Pi57MzfiRlbZmBXwS54aD0wc8BMvB3xNu5wv6PRvqejZyIX5iLGTKQxFzFm4lrY39JcLZcyYxlWH16N+Vvj8McDJwGTFjg0Hkh9C7jQC2YAp90E5OX5ukwmteVq50ptMBNpzEXMXpk0eHeD+uJKw67r6MWjmJk8E+tPrYdGpcHf+vwNsx+YjfZe7eUujYjIrvhaWH9z587Fe++9hyNHjtR4BSNzJlspvFmIZfuWYenepbh08xI8VHfAuOclGHa+Blz3rfK1np7AkiXAxIny1EpEVFldXgsVsQKE2WxGXl5elXtFyfFyyb+Wj4k/T0TP5T2x/tR6/OXuv+D4lONY/vByu00QOFomjoK5iDETacxFjJnYX1paGnbv3g0/Pz+7f2/2tzSl53Kq6BRe2fAKOi7siDm/zUFTXVMsGroIP99fANWWj0QTBACg1wvw9Dyv2EzqS+nnSn0wE2nMRcxemShiksBkMuHQoUNVFqAix8ml6GYRojdHo/OSzlh1eBUGBw7Gvr/tw3ejv8Ndre+yay2OkomjYS5izEQacxFjJvZVXl6OKVOmYNmyZbJcgsr+lqbEXARBwM78nXjs34+hy9IuWHFgBXq06YHvRn+HrNez8Eb/NxD5gBcCA4FblyDRavG/43sUlYktKPFcaShmIo25iNkrE0WsSUCO6Yb+BhbtXoR/pvwT18uvo0/7PogdEovITpFyl0ZERE5q9uzZGDduHAIDA+v0uMq7SjRkF46K5zEYDE61C0dj7yxSoONATwAAIABJREFUORdnbxPUwE/Hf0Jcahz2nd0HFVQYeddIvB3xNsJ9w2E2myGYBBhMBqhUKmzerMXQoQJycwGdDjAYLBMEGzcacfSo7c49Z90BRmpHJalanLlNtthRCbBMTAmCoIg21XS8tjsqVTxWp9Mpok23Hq9Pm4A/z5m6tqm2OElANmcwGbAybSXm7ZiH8yXn0allJ3z+8Od46u6noFYp4uIVIiKSQWpqKvbt24fY2Njbfl18fDzi4+OtH1+7dg0JCQnWj/38/NC7d2+kp6cjPz/fejwkJARdunTB3r17cenSJevx0NBQ+Pv7Y8eOHSguLgYAJCUlISwsDG3atEFiYmKVweKgQYPg4eFR5XsCwIgRI1BaWopt27ZZj2m1WowcORKFhYVITU21Hvfy8sLgwYNRUFCAQ4cOWY/7+PggPDwcmZmZOHnypE3bBKDBbUpKSnLaNpWaSrHl8hYklyTj9LXTcFO5YVjrYRjlMwr+nv4Y6D8QFy9elGzTpk35WLs2H+fONUP79jdw//0a+Pn1w9GjlkwcrZ8A+c69Hj16AABSUlJQUlKiiDbZqp+MRiPKysoU1SZb9FNaWhoiIiIU1ab69lNaWhqAP3+v1LVNtaWIhQuNRiP27t2Le++9t9otp1yRvXMxC2b8cOwHzNo2C1mXs9C2WVvMfmA2/tbnb9BpHGPLL54r0piLGDORxlzEGpIJF9Srm9jYWCxevBhubm4AgDNnzqBt27ZYuXIlhg8fXu3jfH19kZuba/24oVcSHDhwAPfccw+aNGniVO8+NeY7amVlZdZcKrbTdZY25V3Jw6f7P8XKgytxtewqfJr6YHK/yZgUOgk+zXysX1/XNgHAnj170KdPH+v3krufHOHcEwQB+/fvR58+fapsK+vMbWpoPxmNRhw4cAD33XcftFqtItpU0/HatKkil759+8Ld3V0Rbbr1eF3bVFZWhv3791t/19alTYGBgbUecyhikoDktyVnC2ZsmYED5w7Ay80Lb0e8jan9p8LTzVPu0oiIHBZfCxsmICAAGzZs4O4GVC/pF9IRlxqHtUfWwmA2oHPrzogOi8azPZ+Fh85D7vKIiGzK5XY3MJlMOHHiBBe1uIU9cjlw9gAiV0cicnUkjlw8gqn3TUXOGzmYdf8sh5wg4LkijbmIMRNpzEWMmbgW9rc0Z8lFEAQkZidi6LdD0WtFL3xz+BuEdwzHL2N+wfEpxzHpnkk2myBwlkzsjbmIMRNpzEXMXpkoYpLAbDbj5MmT3B7jFo2ZS2ZRJp5e9zT6ftEXyTnJGN9rPE69egoLhy2Ed1Nvm38/W+G5Io25iDETacxFjJnI5/Tp0zVeRWBr7G9pjp6L3qTHqkOrEPpZKIZ+OxTJOcl4+u6nsffFvfht4m8YFTLK5usmOXomcmEuYsxEGnMRs1cmvKGU6uR8yXnM2z4PX6R9AaPZiIc7P4yYwTHo0baH3KURERERVXGl9Ao+O/AZluxdgrPFZ+Hp5omp903FG/3fQECLALnLIyJySJwkoCoEAdi1C8jKAoKDgYgIQKUCrpVdw8cpH2Ph7oW4abiJMN8wfPTQRxjoP1DukomIiIiqyL2Si0W7F+HLg1/ihuEG7vS6Ex899BEm3TMJLdxbyF0eEZFDU8QkgVqthp+fH9RqRdw9YTN1zSUvDxg6FMjNBdzcAL0e8O9UhqcXLMPyozEoKi1CN59uiBkcg0dCHrGu3OtMeK5IYy5izEQacxFjJq6F/S3NUXLZc2YP4lLj8OPxH2EWzOjVtheiw6LxdPen4aZxs2stjpKJo2EuYsxEGnMRs1cm3N2AAFiuIOjaFcjOBoxGACoT0PNbYPBs4I58dGzeEXMfnIvxvcZDo9bU+HxERFQzvhbaB3NWNpPZhPWn1iMuNQ4783cCAIYFD8O0sGkYHDjYKd/UICKyNZfc3eDgwYNc+fIWdcll1y7g9GnAaBSAzuuBl0OBxycCuhJotizA131O4bnezzn9BAHPFWnMRYyZSGMuYszEtbC/pcmRy03DTSzftxxdP+2Kx797HHv/2IvnQp/DkVeO4Nexv2JI0BBZJwh4rkhjLmLMRBpzEbNXJoqYJDCbzcjPz+fKl7eoSy5ZWYDafxfw/EDgmUeAVtnAjpnA4mx4HIpGQa67HSpufDxXpDEXMWYijbmIMRPXwv6WZs9cLt64iDnb5sBvoR8mJ0xG4c1CzBwwE6ffOI2vHv0K3dvYd8eL6vBckcZcxJiJNOYiZq9MFLEmATXM0YtHsfLGTJQ+sx4wa4D9LwHbZwPFHQAAejfLIoZEREREcjlReALxqfH45vA3KDeVI7BFIOY8MAfP9X4Onm6ecpdHRKQYnCRwYfnX8jHntzlYdWgVBAjwyn8KNzfMh+liZ+vXaLVAUJBllwMiIiIiexIEAdvztmNBygJszNwIAOjv2x/RYdF4vMvjTn8bJBGRI1LEJIFarUZISAhXvrxFdbkU3SxCzO8x+HTfpyg3lWNw4GDEDolFG2M/DN0F5F79c3eDoCBg82bLNohKwHNFGnMRYybSmIsYM3Et7G9pts7FYDJgXcY6xKXG4cC5A1BBhce7PI5p4dMQ3jHcJt+jsfFckcZcxJiJNOYiZq9MuLuBC7mhv4FFuxfhnyn/xPXy6+jdrjdiH4pFZFCkdWEfQbAsYpiVZbnFICJCORMERESOhq+F9sGcncf18utYmbYSi3YvQsH1AnhoPfBc6HN4M+xNBLfivY9ERPXlcrsbGI1GpKSkwGg0yl2KQ6nIpbS8FCv2r0DwkmDM2jYLPk19sPbJtdg/aT+iOkVVWflXpQIGDAAmTrT8q7QJAp4r0piLGDORxlzEmIlrYX9Lq08uggDs3An861/Aj1sKMC3x7+i4sCOiE6OhN+nx/qD3UfBmAT4d+alTThDwXJHGXMSYiTTmImavTBRxu4EgCLh06RJkuijCYZnNZvzn1H8w4cAEZF3OQptmbbB0+FL87Z6/wU3jJnd5suC5Io25iDETacxFjJm4Fva3tLrmkpcHDB0KZN88CITFwRjyHaAxIviOroiPisfYnmPhrnXunZV4rkhjLmLMRBpzEbNXJoqYJCCx5JxkTN8yHQfOHYCXmxfmPTgPb4a9ydV/iYiISFZms4CIiZtwtv8CCIFbLQdzBkO9Jxoa1TA8/4ZacVcyEhE5E04SKEzauTTM2DIDSTlJcNO4YZTPKKx4ZgU6tOggd2lERETkwsqMZViTvgbzt8bjjwczLNsupz8DpEYD5/rADCDXzbI20oABcldLROS6FDFJoNFoEBoaCo3GdbfBybqchVlbZ+G7Y99BBRXG9xqPOffPgaZYg3bN28ldnsPguSKNuYgxE2nMRYyZuBb2t7Tb5VJ0swjL9y/H0r1LceHGBbirvKDbFw3DzteBa35VvtbNzbJ4shImCXiuSGMuYsxEGnMRs1cm3N3AyZ0vOY952+fhi7QvYDQb8XDnhxEzOAY92vaQuzQiIqoBXwvtgznLI+tyFhamLsTXh75GqbEUHZt3xNT+U3F3+Yt4ZGhz6PXix7i5AcnJypgkICJyJC65u8HWrVtdauXLa2XX8O7Wd9FpcScs378c/Tr0w46JO7D+r+utEwSumEtNmIk05iLGTKQxFzFm4lrY39Iq55JSkIInv38SnZd0xrL9y9DVpyvWPLEG2a9n462wtxD1QHMEBgLaW65n1WqBoCDL9stKwHNFGnMRYybSmIuYvTJRxO0GgiCguLjYJVa+LDeWY9m+Zfjg9w9QVFqEbj7dEDM4Bo+EPFJlK0PAtXKpLWYijbmIMRNpzEWMmbgW9rc0o8mIxIJEvPOvd7D7j90AgJF3jcS08Gl4wP8B0XbLmzdbdjfIzbVcPaDXWyYINm9WzvbLPFekMRcxZiKNuYjZKxNFTBK4ApPZhG/Tv8Xs32Yj/1o+fJv74qvIrzC+13ho1LxPh4iIiOzvhv4Gvj70NRbtXoTsK9loommCF3u/iLfC3kJXn67VPs7fHzh+3LJIYVYWEBxsuYJAKRMERETOjJMEDk4QBGzM3Ih/JP8DRy8eRUv3lvg48mNM6TcFHjoPucsjIiIiF3Su+ByW7l2K5fuX40rZFbT2aI2n2j6F+Kfj4dvSt1bPoVJZ1h7g+gNERI5FEQsXms1mFBYWwtvbG2q1IpZZAACkFKRg+pbp2Jm/Ex5aD0ztPxVvR7yNFu4tavV4pebSEMxEGnMRYybSmItYQzLhgnr2wTGH7Ry9eBTxqfFYc2QN9CY9glsF463+b+HZns/i5rWbLpuLFFc/V6rDXMSYiTTmImavMYciJgmcmSCIL7XLuHQMM7fOxC8nf4FGpcELvV/AnAfnoINXB7nLJSIiG+JroX0w54YRBAHJucmIS43DpqxNAIABfgMQHRaNUZ1H8bZHIiIn4HK7GxgMBmzcuBEGg0HuUuokLw/o2hUYMgR47TVg8GMFaDHxefRc0RO/nPwFo7uNxrHJx/DZqM/qNUHgrLk0JmYijbmIMRNpzEWMmbgWV+pvg8mAb9O/Re/PeiNydSQSsxMxutto7H5hN35/7nc81uUx6wSBK+VSW8xEGnMRYybSmIuYvTJRzJoEzrY1hiBYVvXNzgaMuiLowz8E7l0Kg7YcTc8PxraZsbjXt1+Dv4+z5WIPzEQacxFjJtKYixgzcS1K7+9rZdfw+YHP8cmeT/BH8R9opmuG1+99HW/0fwNBLYOqfZzSc6kPZiKNuYgxE2nMRcwemShmksDZ7NoF5P5xA8awT4CIjwD368C53sCWWBgKIqF/RgXUbt0fIiIiogbLu5qHT/Z8gpVpK1GsL0Z7z/b4cMiHeOmel9DSo6Xc5RERkZ1wkkAGBpMBXxz8CoaX3wM8zwOXg4ANnwHH/gIIajTxtKxRwNV+iYiIqLHtP7sfcalx+OHYDzAJJnRv0x3TwqZhTPcxaKJtInd5RERkZw1euLCsrAxjxoxBRkYGmjZtinbt2mHFihUICAi47eNsuYiQIAgoLi6Gl5cXVA68wa4gCFiXsQ7vbH0HmZczgRttgN9mA2l/A0xu1q9zcwOSkxs+SeAsudgTM5HGXMSYiTTmItaQTLignn244pijJmbBjI2nNiIuNQ7b87YDACKDIhEdFo2oTlF1bptScrElZiKNuYgxE2nMRcxeYw6bXEkwadIkDB8+HCqVCkuXLsWkSZOQmJhoi6euNQ8PD7t+v7pKzknGjOQZ2H92P7zcvDD3wXn4dsqbyD3pCaPpz6/TaoGgIMsuB7bg6LnIgZlIYy5izEQacxFjJq7Fmfu71FCK1emrEZ8aj5NFJ6FT6zC+13i81f8t9GrXq0HP7cy5NBZmIo25iDETacxFzB6ZNHh3A3d3d4wYMcI6k9G/f3/k5OQ0uLC6MBqNSEhIcMiFLdLOpSFqdRQeWv0QDp8/jDfuewPZr2dj9gPvImmjJzp1slw54Olp+Tc4GNi8GbDFZJkj5yIXZiKNuYgxE2nMRYyZuBZn7e9LNy5h7m9z4b/IHy9teAnnS85jesR05L6Ri1WPrWrwBIGz5tKYmIk05iLGTKQxFzF7ZWLzNQkWL16MUaNG1eprK2/doFarodFoYDKZYDabRceNRiMq3xmh0WigVqthNBqtz2MwGKzHb90WQqu1NPXWQKs7rtPpYDabYTL9+Ta/SqWCVqut9njl2rMuZ+G9He/h+4zvoYIKY7uPxZz75yCgRQDUasvczJ13GpGeLiAlRYXsbOCuu9QYOFANk8kIg0Hc1rq2qXIutmgT0LB+kjpu736qYKtzzxHaZIt+qviaW2tx5jY1tJ8qPi8IAgRBUESbajpemzZVfA+j0QidTqeINt16vD5tAv48Z+raJqLGdKroFOJT47Hq8CqUGcsQ0CIA7wx8B8/3fh5eTbzkLo+IiByQTScJYmJikJmZiRUrVog+Fx8fj/j4eOvH165dQ0JCgvVjPz8/9O7dG+np6cjPz7ceDwkJQZcuXbB3715cunTJejw0NBT+/v7YsWMHiouLAQBJSUkICwtDmzZtkJiYWGWwOGjQIHh4eFT5ngAwYsQIlJaWYtu2bdZjWq0WI0eORGFhIVJTU63Hvby8MHjwYBQUFODQoUPW4z4+PggPD0dmZiZ2H92N7y98j8TCRJhgwsi7RuLZDs/C/Zo7MlIykIEMyTZ5ewP+/qFQqaq2CUCD25SUlNSgNp08edKm/WSLNtW3n/r162fNRCltskU/9ejRAwCQkpKCkpISRbTJVv1kNBpRVlamqDbZop/S0tIQERGhqDbVt5/S0tIA/Pl7pa5tIrI1QRCwM38n4lLj8MvJXyBAQL8O/TAtfBqe6PoEtGquW01ERNVr8MKFFRYsWIB///vf2LJlC1q0aFHj1/v6+iI3N9f6cUPefdLr9UhKSkJkZCTc3d1t9u6TyWTG77+bkZ0NdOpkWSdAp5N+9+mm6SY+2vkRFu1ZhJuGm7jvzvvw4eAPMShokGzvqJWWllpz0el0LvvOZ+XjgiAgISHBmokS2mSrKwk2bdqEqKgoa13O3iZbXEmQlJSE4cOHQ6fTKaJNNR2v7ZUESUlJiIqKgoeHhyLadOvxuraptLQUiYmJ1t8rdWlTYGAgFy60A1suXGgwGJCQkIARI0ZYX0cchdFsxE/Hf8KClAXYd3YfVFBhVMgoTAubhgF+AwCosGuXZQel4GDLuMZWa4E5ci5yYSbSmIsYM5HGXMQakkldXgttMkkQHx+PNWvWYMuWLWjZsnb76Np6pWGj0QitVluP1Xgh+YKZlwcMHQrk5lrWCtDrgcBAy3oB/v5/Pr7cWI7l+5fjg98/QOHNQnT17oqYITF4NORR2VfhbEguSsVMpDEXMWYijbmINSQT7m5gH44y5mgsxeXF+OrgV1i0ZxFOXz0Nd607JvSagDf7v4kQ7xAAtR/X1Jcj5iI3ZiKNuYgxE2nMRcxeY44GX2925swZREdHIygoCIMGDQIANGnSBHv27GnoU9dJaWkpvLzqdm9ddS+YmzYBw4YB2dmA0Wg5Dlg+HjYMyMgAzIIJa46swexts5F3LQ++zX3x5SNfYnyv8Q51GV99clE6ZiKNuYgxE2nMRYyZuBZH6e8/rv+BJXuXYMX+FbhWfg0+TX0w98G5eKXvK/Bp5mP9OkGwjHduN66xxfjbUXJxJMxEGnMRYybSmIuYPTJp8O4Gvr6+EAQB2dnZOHToEA4dOmT3CQKj0Yht27bV6R7Pyi+Yej1QUmL5NzsbeOAB4PRpywtp1e8DZOcI+Pi/GxH6WSgm/DwB18uv4+PIj3Hq1VN4vvfzDjVBUJ9clI6ZSGMuYsxEGnMRYyauxRH6O/1COib8PAGBnwTio10foZ1nO3z+8OfIm5qH2Q/MrjJBAFiumKxuXJOTY/l8QzlCLo6GmUhjLmLMRJpT5yIIwM6dwL/+ZfnXNnf42y0Tx/mL1s5u94L5xx+WKwtEfFNhGjod0w//Dg+tB2ZEzMD0AdPRwr3mNRiIiIiI6ksQBCRmJyIuNQ5JOZZFMh/wfwDRYdEY2Xkk1Krq3/fJygJ0OqC8XPw5NzfL5wcMaKzKiYhcTGPf32UHLjtJUNMLZpW1snwygCEzgS7/hdmswSMdJmH5mDno4NXBbvUSERGR69Gb9Fh7ZC3iUuNw5OIRaFQajOk+BtFh0ejboW+tniM4+M9bDETPr7d8noiIbMBe93c1MsVMElRekb02bveCaTQCHToAf5QUwDRwDtBrFaA2Q5UxGgG58/HznhBn6FsAdc/FFTATacxFjJlIYy5izMQ+ysrKMGbMGGRkZKBp06Zo164dVqxYgYCAALvWYY/+vlJ6BZ8d+AyL9yzGuZJz8HTzxJv938Qb970B/xZ1eycqIsLyJlbFmLWCVgsEBVk+bwv8ORBjJtKYixgzkeZ0udTm/q4GXrplj0xstgViXcm9orMgAF27Sr9gBnYrwoPvxmJl+hIImnJo8gdBlRyLYI97sXkz4OcnW9lERKQgcr8WOpuysjJs3boVw4cPh0qlwtKlS/HLL78gMTHxto9zppxzr+Ri0e5F+PLgl7hhuIE7ve7E6/e9jkn3TGrQ7Y1SV78GBYHjGiIiW/rXv4DXXrMseHcrT09gyRJg4kR7VwXAzrsbOAKz2YzCwkJ4e3tDra7dWowqleWFsfILZrn5JloM/QQX+n2EL45dQ+idoRjXLhat/KJw13Mqm+4nbA/1yUXpmIk05iLGTKQxFzFmYj/u7u4YMWKE9eP+/ftj0aJFtXqsodJ9hGq1GhqNBiaTCWazWXTcaDSi8nsoGo0GarUaRqMRJpMJRUVFaN26NXQ6HdRqdZXnBv58l+fWhaWqO67T6ZBakIq4lDj85+R/YBbM6NmmJ/4e8XeM7joaGmisbVCpVNBqtdXWXt3xO+80Ij1dQEqKCtnZwF13qTFwoBomkxEGg7itdW2TXq+35qJWq6HT6WA2m2EymaxfW1F7dcfr2qbb9ZPUcVv0U13apFarcfHiRbRs2dL6u8HZ22SLflKpVLh8+TJatGhRZQs3Z25TQ/vJbDajqKgIbdu2tT6Ps7eppuO1aVNFLt7e3nBzc3OONgUFQaXXQ+pPRkGvhykgAEIDfpfr9XoUFhZaf9fWtU21pYhJApPJhNTUVIwYMaJOAzR/f+D4cWD77wZ8ffgrbCyZi0v6cwhqFoQVDy/H092fvu1CQI6uvrkoGTORxlzEmIk05iLGTOSzePFijBo1SnQ8Pj4e8fHx1o+vXbuGhIQE68d+fn7o3bs30tPTkZ+fbz0eEhKCLl26YO/evbh06ZL1eGhoKPz9/bFjxw4UFxdbj4eFhaFNmzZITEysMlgcNGgQPDw8qnxPABgxYgRKS0uxbds2AIBJMCGtJA3bDduxM38nAKCPVx882uZRRLSPwJCeQ5CXl4dDhw5Zn8PHxwfh4eHIzMzEyZMn69Umb2/A3z8UKpXt2wRYBtAjR45EYWEhUlNTrce9vLwwePBgFBQU2LxNjdlP9W1Tv379RDt+OXubbNFPPXr0wJEjR+Dp6YmSSu+2OnObbNVPkZGR1tXrldImW/STt7c3IiIinKNNISHw8PFBs/Pnoa40ASJoNChp0wZbr14FEhLq3U/79u1DYWFhvdtUW4q43cBgMCAhIQEjRoyATqer9eMEQcC6jHV4Z+s7yLyciTbN2uDd+9/FpHsmwU0jtb2Bc6lvLkrGTKQxFzFmIo25iDUkE2e6DN7RxMTEYP369UhOTkbTpk1v+7W+vr7Izc21ftyQd5/0ej2SkpIQGRkJd3f3er37dL30Olanr8Yn+z5B1uUsuGnc8Ez3Z/B6v9fRvU13AM73LmFpaak1F51O57LvfFY+LggCEhISrJkooU226Cez2YxNmzYhKiqqyr3VztymhvaTwWBAUlIShg8fDp1Op4g21XS8Nm2qyCUqKgoeHh7O06acHGhHjrSsTaDTQWUwQAgKgnHDBuv9XfXtp9LSUiQmJlp/r9SlTYGBga51u0F9bM3diulbpmP/2f3wdPPE3Afn4q2wt+Dp5il3aURERHQbCxYswE8//YQtW7bUOEFQQWoCR6PRQKPRiI5XtyiUVqu1DsIqbjWo7rmljl8ouYBP932KZfuWoai0CC3dW+Kdge/g1XtfRdtm7bBrF7AmybK4csVigmq1WvIKlepqr0+balN7bY9XTBDcrnZna1NldWlTxeC8ciYVnLVNQMP7qSIXrVYrWY8ztqmm47WtXaVSQaVSKapNtzte2zZV1OA0bQoOBk6csCxSmJUFBAdDFREBncR96/Vt062/V+rappooYpJApVLBy8uryn1N1Uk7l4YZW2YgKScJOrUOr9/7OmbdPws+zXzsUKl91SUXV8FMpDEXMWYijbmIMRP7io+Px9q1a7Flyxa0aFH/hfzqqz79ffzSccSnxmN1+mqUm8oR1DII7z34Hp4LfQ7N3JohLw/o2tept9Tmz4EEZiKNuYgxE2lOnYtKZdnFoIE7GYif1j6ZKOJ2g9rIupyFWVtn4btj30EFFcb1HId5g+YhoEWA3WogIiKqjLcb1M2ZM2fQsWNHBAUFwcvLCwDQpEkT0X3ft5IjZ0EQsD1vOxakLMDGzI0AgP539sfIVtPQ4dpj6HyXxnq1QHW7LQUHO82W2kRE5OBccneDgoICdOzYUXTJyvmS83h/+/v4PO1zGM1GjLxrJGKGxKBn254yVWs/t8vFVTETacxFjJlIYy5izMR+fH19IdN7G1Y19bfBZMC6jHWIS43DgXMHoIIKT3R9AmMDozHz2XC8f8vVAvPnN/qW2nbBnwMxZiKNuYgxE2nMRcxemSgibZPJhEOHDlVZ+OJ6+XW8u/VdBC8OxrL9y9C3Q19sn7gdG57Z4BITBIB0Lq6OmUhjLmLMRBpzEWMmrqW6/r5efh3xqfEIXhKMZ356BscLj2NKvyk49doprHvqR8x8NhzZ2ZbJgZISy7/Z2ZbttKu7ZdTNzXI7qzPgz4EYM5HGXMSYiTTmImavTBRxJUFl5cZyLNu3DB/8/gGKSovQ1bsrYobE4NGQR53zfhYiIiJyWAXXCrB4z2J8nvY5rpdfR9tmbTF/0Hy83PdltG7aGgCwc2f1VwsUFlZ/O4Feb7nlgIiIyJ4UM0lgEkxYfWQ15u6Yi/xr+fBt7osvI7/E+F7joVUrpplERETkAA6eP4jF+xbju2PfwWg2optPN0SHReOZHs/AXete5WuzsixXC5SXi5+nSRPAy8syWXDrmgRBQX/uckBERGQvivjreXP2Zvw9++/IOZyDlu4t8XHkx5jSbwo8dB5ylyYrlUoFHx8fXkFRCTORxlzEmIk05iLGTFzLjvwdeD//fRw4dAAAMDhwMKaFTcPQ4KFQq6Tv4gwOtlwVIMVgAJYsAWbNqrq7QVCQZXcDZzmt+HMgxkykMRcxZiKNuYjZKxNF7G4wPWk6luxdgjfuewPTB0xHC3f7b4lERERUV9zdwD5smfNPeHcmAAAIDklEQVQnuz9BdGI0xnQfg+iwaPRu37vGxwhCzTsYAFW21EZEhPNMEBARkeOry2uhIhYufDv8bWwasQnzB83nBEElJpMJJ06c4GIflTATacxFjJlIYy5izMS1TOw1EUkPJ2HVo6tqNUEAWP7Y37wZ6NTJcqWAp6fl3+DgP68WqNhSe+JEy7/ONkHAnwMxZiKNuYgxE2nMRcxemShikqC5W3Ncyb8Cs9ksdykOxWw24+TJk8ylEmYijbmIMRNpzEWMmbiWptqmuF5wvc797e8PHD8OJCdbbi9ITrZcQeDn10iF2hl/DsSYiTTmIsZMpDEXMXtloog1CYiIiIgcXcXVAgMGyF0JERFR9RRxJQERERERERERNZwiJgnUajX8/PygViuiOTbDXMSYiTTmIsZMpDEXMWbiWtjf0piLGDORxlzEmIk05iJmr0wUsbsBERGRM+JroX0wZyIicnUut7uByWTCwYMHufLlLZiLGDORxlzEmIk05iLGTFwL+1sacxFjJtKYixgzkcZcxOyViSImCcxmM/Lz87ny5S2YixgzkcZcxJiJNOYixkxcC/tbGnMRYybSmIsYM5HGXMTslYkiJgmIiIiIiIiIqOE4SUBEREREREREAGRcuLBJkybw8fGx2fOVlJTA09PTZs+nFMxFjJlIYy5izEQacxGrbyaXLl1CeXl5I1RElXHMYR/MRYyZSGMuYsxEGnMRs8eYQ7ZJAlvjysXSmIsYM5HGXMSYiTTmIsZMXAv7WxpzEWMm0piLGDORxlzE7JEJbzcgIiIiIiIiIgCcJCAiIiIiIiKi/9G8995778ldhK2EhYXJXYJDYi5izEQacxFjJtKYixgzcS3sb2nMRYyZSGMuYsxEGnMRa+xMFLMmARERERERERE1DG83ICIiIiIiIiIAnCQgIiIiIiIiov9RxCRBZmYmwsPD0blzZ9x7773IyMiQuyRZvf766wgICIBKpcLRo0flLsdhlJWV4bHHHkPnzp0RGhqKYcOG4fTp03KXJbuoqCj07NkToaGhGDhwIA4dOiR3SQ5j7ty5/DmqJCAgAF26dEFoaChCQ0Px3XffyV2S7MrLy/Hqq6/irrvuwt13341x48bJXRI1Mo45xDjuEOOYQxrHHNXjmKMqjjnE7Dnm0DbaM9vRSy+9hEmTJmHixIlYt24dXnjhBaSmpspdlmxGjx6Nt99+GwMGDJC7FIczadIkDB8+HCqVCkuXLsWkSZOQmJgod1my+v7779GiRQsAwM8//4znn38eaWlpMlclv7S0NOzevRt+fn5yl+JQ1q1bh+7du8tdhsOYMWMG1Go1Tp06BZVKhXPnzsldEjUyjjnEOO6QxjGHGMcc0jjmkMYxR1X2HHM4/ZUEFy9eRFpamnUm5cknn0Rubq5Lz9bef//98PX1lbsMh+Pu7o4RI0ZApVIBAPr374+cnByZq5JfxYs1AFy7dg1qtdP/Wmiw8vJyTJkyBcuWLbOeL0S3unHjBr7++mvExMRYz5P27dvLXBU1Jo45pHHcIcYxhzSOOcQ45qDasPeYw+l/MgsKCtChQwdotZaLIlQqFfz8/JCfny9zZeToFi9ejFGjRsldhkMYP348OnbsiFmzZmHVqlVylyO72bNnY9y4cQgMDJS7FIczduxY9OjRAy+++CIuXbokdzmyys7ORuvWrTF//nz07dsXAwcORHJystxlUSPimIPqi2OOP3HMURXHHNXjmONP9h5zOP0kAQDRrBt3daSaxMTEIDMzEx988IHcpTiEb775BgUFBZg/fz7+/ve/y12OrFJTU7Fv3z5MnjxZ7lIczo4dO3D48GGkpaWhdevWmDBhgtwlycpgMCAnJwfdunXD/v37sXTpUowZM8blBzJKxzEH1RXHHFVxzPEnjjmqxzFHVfYeczj9JEHHjh1x5swZGI1GAJYX64KCAt7TQ9VasGABfvrpJ/z6669o2rSp3OU4lAkTJmDbtm0oKiqSuxTZbN++HSdOnEBgYCACAgJw5swZDB06FL/++qvcpcmu4veqTqfD1KlT8fvvv8tckbz8/f2hVqsxduxYAECvXr0QGBiIY8eOyVwZNRaOOaiuOOaoHsccHHPcDsccVdl7zOH0kwRt2rRB79698e233wIAfvzxRwQEBCAgIEDewsghxcfHY+3atUhKSqpyX5yrun79Os6ePWv9+D//+Q9at26NVq1ayViVvGbMmIGzZ8/i9OnTOH36NHx9fbF582YMHz5c7tJkdePGDVy9etX68dq1a9G7d28ZK5Kft7c3hgwZgs2bNwMA8vLykJubi5CQEJkro8bCMQfVBcccVXHMIcYxhzSOOcTsPeZQCQq4Tu7kyZOYOHEiioqK0Lx5c6xatQp333233GXJZsqUKfjvf/+L8+fPw9vbG56ensjKypK7LNmdOXMGHTt2RFBQELy8vAAATZo0wZ49e2SuTD4FBQV48sknUVpaCrVaDR8fHyxYsAChoaFyl+YwAgICsGHDBpdfXTcnJwdPPvkkTCYTBEFAUFAQPvnkE5f/4ygnJwfPP/88ioqKoNFoMGfOHDz++ONyl0WNiGMOMY47xDjmEOOYo2Ycc1hwzCHNnmMORUwSEBEREREREVHDOf3tBkRERERERERkG5wkICIiIiIiIiIAnCQgIiIiIiIiov/hJAERERERERERAeAkARERERERERH9DycJiIiIiIiIiAgAJwmIiIiIiIiI6H84SUBEREREREREADhJQERERERERET/8/8Zp0D6vT3zXwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x119ebe5d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## plotting\n",
    "fig=plt.figure(figsize=(16, 4), dpi= 80, facecolor='w', edgecolor='k')\n",
    "plt.subplot(121)\n",
    "plt.plot(x[:-1],y[:-1], 'bo');\n",
    "line = np.arange(0, 6, 0.1)\n",
    "\n",
    "f = sol[0] * line + sol[1]\n",
    "plt.plot(line, f, 'g-')\n",
    "\n",
    "plt.grid(linestyle='--', linewidth=1);\n",
    "plt.title(\"Original data\");\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.plot(x[:-1], y[:-1], 'bo');\n",
    "plt.plot(x[-1], y[-1], 'ro');\n",
    "\n",
    "f = sol_cont[0] * line + sol_cont[1]\n",
    "plt.plot(line, f, 'g-')\n",
    "\n",
    "plt.grid(linestyle='--', linewidth=1);\n",
    "plt.title(\"Contaminated data\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To sum up, outliers influence the estimation of the parameters. And robust function should help against them, however the solution will be not be that simple any more."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MAP\n",
    "How to add prior and make it a Bayesian estimation?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}