hw5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\qquad$ $\\qquad$$\\qquad$  **TDA 231 Machine Learning: Homework 5** <br />\n",
    "$\\qquad$ $\\qquad$$\\qquad$ **Goal: Clustering**<br />\n",
    "$\\qquad$ $\\qquad$$\\qquad$                   **Grader: Aristide** <br />\n",
    "$\\qquad$ $\\qquad$$\\qquad$                     **Due Date: 21/5** <br />\n",
    "$\\qquad$ $\\qquad$$\\qquad$                   **Submitted by: Britta Thörnblom, 900718-1341, brittat@student.chalmers.se** <br />\n",
    "$\\qquad$ $\\qquad$$\\qquad$                   **Submitted by: Adam Tonderski, 930524-1037, tadam@student.chalmers.se** <br />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "General guidelines:\n",
    "* Since there is no theoretical part for this assigment, submit this ipynb only (with completed code/results).\n",
    "* All discussion regarding practical problems, along with solutions and plots should be specified here itself. We will not generate the solutions/plots again by running your code.\n",
    "* Your name, personal number and email address should be specified above.\n",
    "* All datasets can be downloaded from the course website.\n",
    "* All plots/tables and other relevant information should be included."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Practical problems\n",
    "\n",
    "The follwing code might be useful for this excercise.\n",
    "\n",
    "```python\n",
    "import scipy.io\n",
    "mat = scipy.io.loadmat('hw5_p1a.mat')\n",
    "print (mat.keys())\n",
    "X = mat['X']\n",
    "```\n",
    "\n",
    "## [K-Means Implementation, 20 points]\n",
    "\n",
    "a. Implement the basic (linear) $k$-means algorithm as described in the lecture, using the euclidean distance. Use (uniformly) random points from the data as initialization for the centroids. Terminate the iterative procedure when the the cluster assignments do not change.\n",
    "\n",
    "b. Run your implementation on the matrix $X$ in **hw5_p1a.mat** with $k=2$. Each row of the matrix is an observation, and each column is a feature. Store the cluster assignment both after 2 iterations, and at convergence. Produce a scatter plot of the data with colors indicating the cluster assignments at convergence and highlight points that have changed assignment after the second iteration.\n",
    "\n",
    "c. Implement the kernel $k$-means algorithm as described in the lecture, using the Gaussian RBF-kernel.\n",
    "\n",
    "d. Run the linear $k$-means **and** your kernel $k$-means on **hw5_p1b.mat** with $k=2$. For the Gaussian RBF-kernel, use $\\sigma=0.2$. Produce scatter plots of the data, with color indicating the cluster assignment at convergence, one plot for each of the algorithms.\n",
    "\n",
    "## [Discussion]\n",
    "\n",
    "In the first dataset we can see that the classes are linearly separable so the basic implementation works fine. A bunch of points (marked with white stars) changed assignment after the second iteration.\n",
    "\n",
    "In the second dataset the classes are arranged in circles, which makes it impossible to divide it into 2 clusters with a straight line. The basic algorithm might work more decently if we increase the number of clusters and let several clusters represent the same class (in this case the outer circle could be represented by many clusters arranged in a circle).\n",
    "\n",
    "However, a much better approach is to use the kernel version of $k$-means and replace the euclidean distance with something that is more suitable to our problem. As you can see in the last plot, a Gaussian RBF-kernel fits this role perfectly (since the data is circular in nature)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['__header__', '__version__', '__globals__', 'X'])\n",
      "4\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10e986e48>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAHVCAYAAAAkfVjMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+QHOV95/HPI/RjJZbyXBlXAcE6QOf4EuNcUtKm7gpK5fhsRXICPta54ATj2EoVMiius0G4zvySI0Wqu4AJiY9yxFmQhNs9nPKhMw7GYBs7JFXnaEWOxCY2IZDygS2XHWMwq52VEHruj0edHe329PSPp7uf7n6/qrZGuzua6emd6e/z4/t8H2OtFQAACMOyug8AAAAsIDADABAQAjMAAAEhMAMAEBACMwAAASEwAwAQEAIzAAABITADABAQAjMAAAFZXseTnnnmmfa8886r46kBAKjc448//k/W2teluW8tgfm8887ToUOH6nhqAAAqZ4z5dtr7MpQNAEBACMwAAASEwAwAQEAIzAAABITADABAQAjMAAAEhMAMAEBACMwAAASEwAwAQEAIzAAABITADABAQAjMAAAEhMAMAEBACMwAAASklm0fAQAt1O9LBw5IzzwjrVsnTU5KY2N1H1XjEJgBAMXNzEibNknHj0tHjkinny5t3y498og0MVH30TUKgRkAUEy/74Lyiy8u/Gx21t1u2iQdPkzPOQPmmAEAxRw44HrKcY4fl+6/v9rjaTgCMwCgmGeeccPXcY4ckZ59ttrjaTgCMwCgmHXr3JxynNNPly64oNrjaTgCMwCgmMlJafmQlKXly93vkRqBGQBQzNiYy77u9aTxcckYd9vruZ+T+JUJWdkAgOImJqTvftclgj37rBu+Zh1zLgRmAIAfq1dLv/7rdR9F4zGUDQBAQAjMAAAEhMAMAEBACMwAAASEwAwAQEAIzAAABKTwciljzJikxyStOvl4n7HW7iz6uAAAT9gnuVF8rGM+Kumt1tpZY8wKSX9pjHnIWvs1D48NACiCfZIbp3BgttZaSSc33tSKk1+26OMCAApin+RG8jLHbIw5zRjzhKTvS/qitfavfDwuAKAA9kluJC8lOa21r0r6WWNMT9IBY8yF1tpvDN7HGHOVpKskae3atT6eFgCQJIR9ksuY3275nLnXWtnW2heNMV+RtFnSNxb97i5Jd0nShg0bGOoGgLJF+yRHw9eDqtgnuYz57Q7MmRs3RVzgAYx5naRXTgbl1ZIekfRfrbV/Nuz/bNiwwR46dKjQ8wIARpifl84++9Q55kivV+4cc78vnXOO3+cu4zErYox53Fq7Ic19fcwxny3pK8aYv5U0IzfHPDQoA8Ho96XpaWn3bnc7P1/3EQF+1blPchnz2x2ZM/eRlf23kn7Ow7EA1enAcBggqb59ksuY3w5hzrwC7MeM7mEJCbqmjn2Sy5jfrnvOvCKU5ET3dGQ4DKjV5KS0fEjfb/ly9/sQHjNABGZ0T0eGw4BalTG/XeeceYUYykb3dGQ4DMjF5xrhMua365ozr1Dh5VJ5sFwKtapzCQkQsrikyOXLSYr0IMtyKXrM6J5oOGzYBYigjCbwXf2KpMhgEJjRTR0YDkOLlbHcL01SZNWZ3R1FYEZ31bGEBCiqrJ4tSZHBICsbAJqkrOV+UVJkHJIiK0VgBoAmKatn25E1wk1AYAaAJimrZ9uRNcJNwBwzADTJ5KRL9IpTtGc7mBT51FPSP/2TdOaZ0tNPS29+M8G5IgRmAAhZ3LKoMpf7rV4tveENLvizyUstCMwAEJooGD/2mHTvvdKyZUsDZFnL/VjPXDsCMwCEJFqj/MorS5O8FgfIMpb7sZ65diR/AUAoBnurwzKvpXJ3QWM9c+0IzAAQiqTe6qAyAyTrmWtHYAaAUCT1VgeVGSBZz1w75pgBwKcim0skbUk6qMwAySYvtWPbRwDwpei2iUlbkkqu4EfRbRjTNhyi+6XN+va921XLZNn2kcAMAD70+9I552Tf53txQFu7VnrHO6SXXz71fitXSjfdJF1/ff6AV9Z+y+zjPBKBGQCqNj0tbdsWPww9Pi7t27d0mdGwgBa3VEpKDvCjJDUcXvMa6Y47pOeey97bzdsg6ZgsgZnkLwDwIesyo8GlUbOzkrXuNmmpVJFlUkkZ3y+9JF19tbRzp2tcnH22azQUfdwyl3W1GMlfAOBDUuJWXBZ12qVRg4oskxqV8T0/726zVvnyue6ZeWpJ9JgBwI+sy4zSLo0aVGSZVNL65DiLe7v9vhuu373b3UaB3Ne655kZNyS+bVu+nnuLEJgBwIes2yZmDZRSsWVSSQ2HOIO93aSg6WPdc9Kw/qZNC42AjiAwA4Av0baJ+/ZJu3a528OH4zOTkwLaGWe4hCyf+yLHNRxWrRp+/6i3OypoSsX3cWae+hTMMQOAT6tXp9vkYVQhjwsv9L971OB+y88+K517rvThD8dnVEe93fvvT7epRZHdrqjPfQoCMwDUZXGgXBzQytjFaXHD4U1vSq7ylTZopm2QxMmaONdyBGYAqFORgObDqMZBFUFzctLtMx2n6vrcAWSGU2AEADBcUplQnwVEQqgeVuIxZCkwQo8ZADBcVZtajOq5l20wyS2SdU23JwRmAECyqoJmncP6aTLDKzo2AjMAYLS658LLFlBmOIEZAKoSQGIRhggoM5zkLwCoQgjJTU2Sdd/ooo2dkpPc2PYRAELC1ojZpG3E+G7skJUNAB0RUGJR8NJmR5eRRV13ZvhJBGYAKFtAiUXBS9uIKauxE0CSG5tYAEDZfG2N2AVpGzEtbuwQmAGgbD62Rlxs2P7ITZe2EdPixg7JXwBQBZ+JRW3O8E6bHV1VqVBPSP4CgND4Siyqs3RkFeuw05YArapUaA0IzABQFR+JRVVneEfB+LHHpHvvlYyR5uZcENy+vZxeetpGTCBZ1L4RmAGgSapMehocMl9cEavsXnraRkwAWdS+EZgBoEmylI4cNvScZkg6bsg8DuuwvSMwA0AI0s7fTk66IeQ4gxnecQli27dLd9whfehDS3++eEg6ach8UMOXJoWIwAwAdRsWROPmb9MkPSUliL3//dLgapzo529/u/S97y00BpKGzAc1fGlSiAjMAFCnPFnWo5Keknq7w5bIvvSSdOut0s03u++ThswH5V2HjaEoMAIAPuQt+JEmyzpOlPR0003udjB4p+3tLrZ378JxJxVFkVzQ7vUavzQpRPSYAaCoNEPRw+aQy8iyTtvbjRMlcg0bMj9xQrrySmnjxlYsTQoRgRmAH1UUnwhRmqHor399eODOkmWdVlKCWJKjR09tCLR0nXDoCMwAisuSvNQ2o4ai77tP+vCHhwfuf/zH5Drax4654fEsjZ3Fvd20Pee4hkAL1wmHjjlmAMUM9hhnZ11y0eys+37TpvZsrjDMqKHoBx9MDtyf/7wLor2eND7uKmuNj0tnnOF+/8EPSjt3Stu2udrQMzPpjivq7e7bJ/3Kr0irVo3+PyRyBYHADKCYvMlLbTFqlyNp9BzyYBDdtUv6xCekZctcA6dIYyfq7d57r/t30v1I5AoGgRlAMS3eFzeVUVs6/tIvpduecDDLeuVK6dVX4/9PnsbO2Jj0wAOuNx7HGDek3vZph4YoHJiNMa83xnzFGPN3xpgnjTH/yceBAWiIFu+Lm0o0n7t4KDrqgb773dn3Yi6jsfPcc8P/TsuWuSH1tNq6F3QgfCR/HZd0nbX2r40xZ0h63BjzRWvt33l4bAChS1siss1GZS9n3Z6wjExtX8G+y4l+FSkcmK21hyUdPvnvl40x35T0E5IIzEAXtHhf3EySspezLjsqo7GTNdjHLX+ztr69oDvE2GHl2fI8mDHnSXpM0oXW2h8v+t1Vkq6SpLVr167/9re/7e150VBdXffaVtHfk/WufsT1TKPGTp6e6fy8y+qO2y2q1zs1qA577muvlX73d+OD+/i4S15jaVUsY8zj1toNqe7rKzAbY8Yl/bmkPdbaxMyEDRs22EOHDnl5XpSozMDp+6IDtJHvxk6az12/L51zTnwAHxsbPp9sjMsov+mm/MfXYlkCs5cCI8aYFZL+l6SpUUEZDVHmPFKeov1AF/ku7pFmSH3Udo+rVrkKYYt1IdGvIoUDszHGSNov6ZvW2tuLHxJqV3bgTLPuleEwoByjgn1SktjRo8MLlXQl0a8CPnrMF0m6UtLXjTFPnPzZDdbaDLn3CErZgbPuda/MbcfjvDRP0b9Z3P8flST2kY9It9/e7US/kvnIyv5LSUNWraORyg6cZSwFSYulHvE4L81T9G827P9/7nPJ666vv17asYNEvxJ5zcpOi+SvwE1Pu7q8ZWVeZskO9SkpqaXM5w0d56V5iv7NRv3/z3xGuvRSF7SPHXNBe8UKGmoFZEn+oiQnlhpVYrDoPNKoSkllBYEiNZ3bXOmo67Wum6jo3yzp/x896oLysmUuKK9c6dYvf+5zBOWKsO0jlqqiYEQd+7zmHaJv+zBv3XP+yK7o3+xb3xq+FWS/f+r3x465r0suqXb0pMM5DwRmxKsicFa9z2ueue0uLO2qc84f+RT5m83MSB//ePbnrHLFRNsbwyMwx4zmyduSzjO3XfZ8ewjqmvNvqhB6cnn/Zklzy6NUVUCkpTkPzDGjvWZm3Id227bsm8fnmdvuwjBvXXP+TVTk/edT3r/ZqOIhK1dKa9bE/66q0RNyHhjKRoP4GFbOOkQf6jCv715bHXP+aYXQQ42OI6RpjTx/s6SGpjHS5s0usMepqoBIFxrDo1hrK/9av369BTKbmrJ2fNxalyN66tf4uPu9b/2+tb1e/HP2eu73VTt40D33+Li1xrjbXs/9vG1Ceq11vP98S3oNkrWrVp36fR3nvA3nOYakQzZljGQoG81RR0s6tGHewV7b7Ky7XM3Ouu83bWrXMq7QXmsbenJJSyGlpTWwV62SPvEJNxpQVdJV2cs1G4DAjOaIhpXjlDmsHA0Z7tvnkl/27av2QjWoS/Nvob3Wut5/PsU1NJMal8uXu3nnKhugoTWGa8AcM5qjjM3j06p6adcwbei1pRXaa63z/efT4rnpv/kbV+krTl3vqZBzHipAYEYxVSbmVFH4JHShJqOVIbTXWsb7r67EtsGG5vS09IUvxJ/nlSulc88t/3jihNIYrgHrmJFfmk3Xy+B78/gm6dKa41Bfq6/3X12fn8WSzrO0MITcgcIeZcqyjpnAjHxaWgSgEUK5oFehra81tM/PzIz09rdLL70U/3s+04VlCcwMZSOfsvdsxnBdmn9r62sN7fMzMSHdcYd09dXx2e51faZDWcNeMQIz8gktMadrujT/1sbXGuLn57nnli6XitRxTB2ul01gRj6hJeaEoqMtfGQU4ucnpGMKrcpaxZhjRj6hJubEqSpYtnU+FP6F+PkJ6ZhauHkMc8woX1OWLlU1HNbxFn5rVNWIC/HzE9IxhTjUXyECM/ILPTGnymAZWjJPiEIf5q96TjPEz08oxxTSsHoNCMwoJuTEnCqDZcdb+COFnshT14hHiJ+fEI6pLVXWcqJWNtqrymDZhjrKZQltM4o4odXl7rqO18umx4z2qnI4LLQWfkjDxk0Y5mfEIzyhDKvXgMCM9qoyWIaUOBPasHETgl7H5zSDFcKweg0IzGivqoNlCC38ELPDswa9Onr7oY14LOb7nIQ0ooIlWMeM9uvSphchrv/Msj62zrXgoa5D931cob7OlmMdMzCoS8NhIQ4bpx25qLu3H8KIx2K+z0nd5xipEJiRjCGvZgl1rjRN0AshSSy0RpzvcxLCOcZIBGYMF1oSEUYLea50VNALsbdfN9/nhHPcCKxjRrwmrD3FUk1e/8la8KV8nxPOcSOQ/IV4ISYRIb0mJryFtIlCKHyfE85xbUj+QnEMeTVbaHOlaYS0FjwUvs8J57gRCMyIF2oSEcoTQqJfiJnRdfN9TjjHwWMoG/GaPOQVQoBpGta2AqViKLslao0vTR3yIpM8O9a2hofGZafRYw5UFF9e+1pp/35p61bphRdqiC9NSiLq96VzzmlmL79OJPqFhdGLVqLH3HCDHZg9e6SLL5auu851/CrvwDQpiYjiCfkkJfrNzkpPPVXt8XQZoxcQ65iDdOCA+/xZK11zjXTaae7WWvdztoYdgkzyfJLWtkrSbbe5XhzKx77QEIG5FP2+Gx3cvdvdZq3F8cwzbtT4nnukuTn3s7k56e67pfPPJ74MRfGEfCYn3VDpMHNz2YvKFP0QdBWNS4jA7N3MjJvmvOUWaeNG6eabXXJzlg7HunXSyy9Ljz4qrVjhPo8rVrjvZ2eJL0MlBZi6y1GGLEr0W7Nm+H2y9NaiD8G2bdLOne4264egq2hcQgRmrwanh669dmFuOGsVyyi+bN3qvv/Yx9zt1q3El0RNLkdZt4kJ92YdJm1vjVKuxdC4hAjMXvmaG47iy+HD0qWXSh//uPTOd7rviS8jRMUT9u2Tdu1yt4cPk82axr/+164hE8da9zUKc6TF0LiEyMr2Kpob3rNHuvxyNzI4Nyfdd590443DN/2JMzEhXXih+/dFF4W/UikoTcokD0nSzlSSayFef33ym5A50uKozNV5BGaPBueG3/Medx1auTL/3DDxBZUaG3NzMLfcEv/7V18dveSMUq5+8OHvNIayPWJuGK125Ihb05yUbc0cKVAYPWaPoumhp592c8MPPyw9+aR0xRVMD9WG0obZrFvn5jTjeryrV7s1zcuWDS932tRSrm3Ce77xKMlZgiZVsWw1Shtml7R5iTHxCWBx5U75ENSD93ywspTkJDCjnaibnV/cxf3ECReU+/2l96eedhh4zwctS2BmjhntxLKd/OKWnF133fA1yGRbh4H3fGswx1wjpoJKxLKdYhZnBU9Pk20dOt7zrUGPuSY+SnciAaUN/Qol25oa3MPxnm8N5phrMDgVdOedrpTwvn0uwZWpIE+Skpg4yfnUnVhU9/OHjvd80Ej+Ctz09PBh6/l5tg32hgu5f3VlW+dNbOrafBHv+WBlCczMMdfAZ+lOJKC0oX91VaRKk9i0+LjigtTidddtw3u+FQjMNfBduhMJKG3YDlkTmwZ3uYpEiWubNrV7WJf3fON5Sf4yxtxtjPm+MeYbPh6v7SjdCWSUNbGJpUNoMF9Z2X8kabOnx2o9tnUEMsqaFc7SITSYl6Fsa+1jxpjzfDxWV7CtI5BB1hrc7HKFBvOWlX0yMP+ZtfbCUfftelY2AlVHBm/XsoaLSpsVztIhBKaW5VKjArMx5ipJV0nS2rVr13/729/28rxNw3U4UFUtMxl8A0jS7bf7eU7eWEuxdAgBCTIwD+pqjzm6Trz2tdL+/S7Z64UXuE7Urqri/4OBIm6INe9zEoCGY5crBIJ1zAEaXL2xZ4908cVuX4Dt29u/eiN4edbISkt7qVu2SA89FN9rjVu+M8zgc47qCXd5WVAaLB1CA3kJzMaY/ynpLZLONMY8L2mntXa/j8duiwMHll4jr7nGfTW52lcVI6ilP0eeDN7FvdSxMfeHXL3aHfDiYhZJwX/Yc6YpkJG3UQEgWF6WS1lrf81ae7a1doW19lyC8lJRta977nFVviR3e/fd0vnnN3P1RhUbcVSy2UfWNbKDvdTZ2YV9iq11f1Rr3c9ffNHdb34+OfjHPee55y59jsWPKbEsyCc2yEAg2F2qIoPVvlascNfMFSuaW+1rMDZde+3C0Pxg3Ch6nUvzHF5kXSObpfcb9VrXrXO1V9OIjiVNgQx2FPIjagFu2ybt3Olu2e4NNSEwV6Rt1b6ioXlr3XD8aae5W2vdz2+9tXhPd9RzeCveFK2R7fWk8XHJGHfb68Wvkc3S+416rZOT7sCTrFq18JzPPTf8OWZnF3rCoWzH2GRxIyBxoxNARQjMFWlbta+kofnzzpP27i3e0610+D8q/r9vn7Rrl7s9fDg+qzmpl7pY1GsdG5OuvDL5vpdcsvCco54jCvJJjYoHHnCtF4Zmk1G+E4EhK7tCbar29frXu6VecRtxvPii9KMfuTyoSJ5Et8o3+0ibwTs5mX4LsMFe68aN0tRU/FKp8XHpsssW3gyTk+6EDXP77dL117v7x+0otHatC/Rd2lkpL+bpERprbeVf69evt2iugwet7fWsPessa3/0I2uPHbN2xw53++UvWzs2Zu3ZZ1t7993WHjni/s+RI9bu3+/+z+7d6Z6n33fP8+ijS5+j13O/r010EsbHrTXG2tWr3e2aNe52fNz9/uDBhf8TvSDX3z31K+4F7doVf1/JPf7UVPyxzc1le56um5py5zPreQYykHTIpoyRDGUjk8HpuJtvll7zGunBB6XbbltYNnvDDX4S3YIe/l889P2pT0k//KH03//78KHwrHPZkrtPnKSeHEOz2XR9np5s9OAwlI1M4tZj/4f/sLBi6NOfdiOst99+aqLb3r3u+wcfzHadC3L4f/HC6h07Fg5o1FB4lo3s827EwNBsNlk3yGiTNGvlUTkCMzKJErL27JEuv9ytAJqbk+67T7rxRveZjq5zTz/teroPPyw9+aR0xRX5rnNBFW/ycSHzMZed1JNjZ6XssjSY2oKqccEiMCOTtAlZQfZ0i6r6Qpa3JzcqoG/Z4oYs2fDiVEG1ACtA1bhgEZiRSXTNTzNM3brrXB0Xsjw9uaSAfscd7jFCHrpkp6xqMOURLAIzMvE9TN0odV3I8rRw4gL6li3uNuShS+Y8q8OUR7AIzMislcPUi8X12pp2IVsc0Kenwx66ZM6zWnlzGFA6AnOH+BwhbN0w9aBhvbbPfa7Zy2pCH7pkzrNaXc5GDxyBuSMYIUwpqdd2ySUuOC+uqNWUC1noPf7QGw5t1MVs9AYgMHdAG0YIK8sHGtVr+3//r7kXstCHLkNvOLRVq4e/monA3AFFRghDSJCttLefptfW1AtZ6EOXoTccgIoQmFtgVPDMO0IYwvB3Um//LW+RvvMdV9HSm7b32rIOXVbZMgu94QBUxNhRe8SWYMOGDfbQoUOVP28bxQXP6DoWBc/pabfv+7BNjfbtW9oB7PfdfsqDATHS61U3/J107JKrPPbVr3psKMzPu42j637hIUjz5ipD1Bho2lQBkMAY87i1dkOa+9JjbrBRvckdO6Q3vlF6xzuyJxOHkiCb1NuXXDlQr/Pk9NqcpDfXxo3SJz8pvfvd5ZyPpk4VAJ6wu1SDJQXPuTm3ydG2bdL557uCT73ewnV0bCx5U6NQEmSjkeUk3jdMWrxzVNxOUW2X9Oaan5c+8AE3sjAzU+1xAR1Aj7nBRvUmpYVOzgc/6HYQjGYuog1nhwllqjUpHyhSSkMhS68thAw530a9uY4edV9NSesHGoQec4Ol6U1GXn5Z+vGP3bVUcrcvveSuq3Hbr4ayRW00srxmzfD71JqTNTPjJuNvucUN8d58czt6kmnfXG3Y35n9iBEYAnODJQXPtIZdV6OA2Ou5BDFj3G3S8HdZJiak558fHpxrW0kzOA977bXSxRdL113nvh/W4mmKtG+uphf+iBpW27ZJO3e62zY0rNBoDGU32OI8pWGZy0mSrqshFQX6F//CZV8HlZN14MDSYdxrrnFf8/PNLiE5+Oaanx/eyBg2XNGE4f02VN5BKxGYG24weD71lHTbbS7xK61Rw8AhJciG1FCQ5ILOBRdIe/ZIl1/uuvRzc9J990k33jh6cjx00Qn/9Kelq6+OD85xwxUhLIBPI5SlB8AiBOaA5O1kDAbPX/7l+F7lsB510woqhdRQ0Lp1bvL+0Uel97zHnfCVK933s7PNL0YiuRP+vvdJb3pTuuGKJvVCQ1l6ACzCHHMgfE11DVvp8+ijYcwXt0o0D7t1q/v+Yx9zt1u3Nq/FM0raJWRpeqGhSEpwa0OVNzQWPeYA+O5kxPUqgxsGboNoHvbpp6VLL5Ueflh68knpiiva2eJJM1zRpF4otbkRKAJzAKqa6gpqGLgtJiakCy90/77oonJaPE1IpIqEsgA+Daq8IVAE5gA0qZOBGGW2eKJEqte+Vtq/3w2Th5hIFWlaL5ShJASo0YG5SR2JJE3qZKBCg3Mce/YsrJPevr26RKqsH7JRvVBrXRGPkD60XR9KasuFtEUau7tUXRvflIENjbLpzHVkenr4i6tinXRcb/2FF9J9yOJ2iPr61/M/HsrRpgtp4LLsLiVrbeVf69evt0XMzVnb60XVnk/96vWs7fcLPXwtDh50xz4+bq0x7rbXcz/Hgk6dp127rD37bGvvvtvaI0fcz44csXb/fmvPOsvat77V2qmpct7wgx+yO++09vhxd5v3Q+b78VBcGy+kAZN0yKaMkY1cLtWkFRlpsaHRaIMju7Oz7goyO9uOCpixBtdJr1jhejQrVrjvJVefu6za3FFVM2tdJbPTTnO31rqfZ/2Q+X48FNfGC2lLNDIwtzVZKprquukmd9vK4dkCkq4j/b505ZUt24MgaZ30zTeXW5s7qmp2zz0LpeTm5qS773b7iGb9kPl+PBTX1gtpCzQyMFMXoJuSriNHj0qf+Yz0/ve7Kcy//Mtqj60UUSLV4cNunfTHPy6dOCG99a3l9zqTeusvvOCeL0tDIOnx2lIlrWmSLqRr1ri/MTtu1aKRgTmULQm7pu7d8dLsRHjsmOuIbdzYkuA8MSFddpkbDviFX5DOO6+aXuew3voHPuAaCw89lG0IvUtV0poi6UI6Nyf9yZ+w41ZNGhmYQ9qSsCtC2B0vyzaX1kq/+IstaehHcxy/+Zuud5ml15m3NbW4t37bbQsZ1RdfLO3YkW0IPa73/853uu/50NYj7kJ6+unuNkrgaH0iR5gau1xKil+Rwefbv37fBeUQlnMNru4Ytc3lypWuc9maJarRurr773fB8YYbpL17pb/4C+ld71r6hyiy3CnS77sgfNttroEQd0xZlm3xoQ3P4N/k8GHXU477cI2Pu6zU1nygqpVluVSjC4x0vS5AVULaHW+wUNOBA9L//t/Dj+2VV5qZvzJ0nXaW2ty+ipOsXi2ddZY7EB/bW/KhDc/g32T3bhLCAtDowIxqhJa8GV1HJiddZ3BYYG5iIuDIrYzT1uaOlicN/vyaa9xX1l5uF7a3bLu0VXkoQxiERs4xo1oXESTQAAAgAElEQVS+suB9J4+NjblOozHxv29aTlHqddpp1tX5Wp7U77uMuuPHiyVu1Z052GVRgsgtt7isyKS172TWBoHAjJF8fFbLSh67+GLpz//cja6uXDk6ETDk+OC13oOP5UnRH+13fkf6whek739fevvbsyduZQkM8GuwtXfttaPXvpNZG4a0JcJ8fhUtyYml5uZcdcZdu8qp0likFGYVlf/m5qy95x5rf+VX3Nc99yx93NDLee7a5Y4r7jwZY+3u3RkerN93L+7RR609dszaHTvc7Ze/nO6kJ5XQHBuLP8FZH4eyj+VLuhj0++73caILyu7d5ZV97RhlKMlJYG6BqgJO3s/q1JQ7priAMz4+/NqQxahz0ISywLnOU1KL7OBB97PNm91J2bLF2v/xP9K9MdJc0NO0BvMGBvgxqt56ptYeisgSmBnKbrgq60fnLRladvJYmnNQd1ngNEPomacMRg0RDxYn2bXLJW69613plkqNmqN+7LF0w9OU4qwXFdcaicDccHUHnDTKLqGa5hzUmVmedn490/Re2rnDvK2ppAv6yy9L996bbt6SwFAvKq41EoG54UJbyhSn7ETPNOegrvrqWUc0Uu8yVvZuTUkXdGulH/4w3XMTGOpFxbVGIjA3XBM29Cg70TPNOahyFcjgsPWOHdlHNFJ1csseIk66oF95ZfrnrjowhJx2X5ciUxqoBQVGGm5ycnjhpZA6JIMVu/JWYxxWIyHNOYjiw+LiHcuX+40PiwuErFjhlgHHKTSiUUXRj2HFTIyRpqbSP3faoihFxZUgPaU6S4dRca1Z0maJ+fwiK9uv0JcB+TDqNaY9B2WuAknK/PaekZ51OZTP9XRFl2KVgWVZCJwyZGU3ehMLLGjz3gBpN9Go+xxMT7vErlGba0QKbwAyM+PqZt97ryuBtnmzq5v9kz95ag/Rx2YWeZ+7KtPTw//gWUuQAiXozCYWWNDmkaq0m2jUfQ6SktAkN9r7yiseh9DTDBH72swiz3NXKZpz97HRBlAzAjOC14TMc2l0/f/f+A23TMprDBvVGvG5mUXW564SG22gRcjKRvCakHkuJWd+r1jhEpKzLicurCsFPliWhRYhMCN4Tdnwps76/0NXCXWlwEdT1+uyvAsxvCR/GWM2S/p9SadJ+pS19r8k3Z/krzCl3bK1DouXIa1a5X5+ww3S9deHc5xS9UlocXs4R3PYE2+ed+Pn99/v5pdvuEHau1f6i79wa1kLZZ8FqO4MwCzKSMpDsLIkfxUOzMaY0yT9vaS3S3pe0oykX7PW/t2w/0NgDk/ixT2Qa0S/L912m4sr1kpHj7oeaWjHWaVUGetfDyyDGqf+4e6806Xz79vnktQKp+sjRFVnZf+8pH+w1j578snvk/ROSUMDM8IymLgbiRKYiiTu+tTvS5/+tAvKg6N9oR1nUVlHLdJlrAeWQY1yk/LQeD7mmH9C0nMD3z9/8menMMZcZYw5ZIw59IMf/MDD08KX0DfCiDaBuPrq4VNwZR9nFVOBaTe7GJQ6Yz3vZhYoR1eS8pBLZclf1tq7rLUbrLUbXve611X1tEgh5OVIg735pGA4Oys99VQ5x5AnYGaVd/vOpmSsY5GuJOUhFx+B+TuSXj/w/bknf4aGCPnintSbX+y22/wGS6m6/a7zjlo0JWMdi7C8Cwl8BOYZSW8wxpxvjFkp6d2SHvDwuKhIyBf3UdW0Bs3N+Q2WUnXD/HlHLepcooUCmrq8C5UonPxlrT1ujPktSQ/LLZe621r7ZOEjQ2Wq2nkpj6RqWnGOH3dVGFeu9LPsq6ph/lFVw5JGLZJ27gp5CVznhVbWFMFgEwv8syqWgGYNFPMnl+HGLQeKY4xb47x8uZ9lX0kbU4yPuxUuPpJnk15n3tUzTVgCB3RFpeuY8yAwd1PeQBFXXCTrcHXe4FZGwBzGZyBNuyMXgGqwuxSCU2St9OKh2nPPlT70Iemll9I//+AuVFlUOcyfNCSdVdoduVqJ8Xs0HIE5g7Z+3qt4XUUDxeKNjN70plODZTSfOkyR+WCfATOJz79DyEvgRipyIuLKXG7fzvg9msVaW/nX+vXrbdMcPGhtr2ft+Li1xrjbXs/9vMmqel27drnHdwuOTv0yxtrdu7M/5tyctVNT1t5yi7Vr1sQ/dvQ1Pu7uG6K5OXd+xsbcl4+/w9SUe4ymnYt/fkOuW2ftV79q7QUXpD8Rc3PuvpK1d95p7fHj7lZyP+/3yz9+YAhJh2zKGElgTmHw8774q8mf9ypfV5mBIumxQ/87HTxo7Wte4/+4+/0GvmeLBtapqeH36fcDbo2gC7IEZrZ9TCH0kpV5Vfm6ylwrPWqtczRPnGdYuMxSnNG8e9Jced6/QyPXN0f1o611NaNPO83dWut+PupEpClzyTaLaAACcwpVztdVed2o8nWVGSiSKpeNjUmf/GS+6cWyS3GmqWrmY2583z5p1y53e/hwwFOtg4E1Shjo99PXjx5V5tJa9we95RZp40bp5pv911YFfEjbtfb51bSh7Krm66qex65jHjKaF969O3nkMYsyhm2rGOZPmndvxHywb/fc4170e95j7auvWnvihLu94gr387vucidj1674N0/0Rnj0UWuPHbN2xw53++Uvu/mC1avzD5MDBYk5Zr+qmK+rYx477+uKguuw62MdfDdqqmi0NHlu3Lu5OWuvusrdDvv96tWjk8IOHnQndvNm90bYssXae++1dtWq4Y/N/DMqkCUwM5SdQhXzdXXMY+d5XVXstJSH72HbKob5k+bdJek1r5GuvVa69daWT4dGb6ovfEF64gnpH/9x4cXOz0tPP+2Gsvt9d0Iuvli67rr4nUQmJqTLLpOuvNK9EX71V6Xf+i3p6FG3obe10iuvuPuyzSICxTrmlMpey1rXutMsr6tIkZAqLF7rXESR2tVpWevizN697t9Hj7qqZsZIv/Ebrub37/7uQlGTVi7HHXxT7dkj/fzPS1/5iqsic+SIK3r+27/tgvXq1Qv/75pr3Nf8/NJF8INvhOlp6XvfO/WNuWLFwv3YZhEBIjBn4PPCv1gVgWCYtK+rC9WkotoW3/qWdOJE/H187Lg1WH4zCshjY9INN7gO3gUXnJqtHVIDyKsoE3vwBb3tbafeZ+tWd0L27pUuv1xas8b1du+7T7rxRtdiGSZKKNuzZ+H/Wuv+AD/3c+6xH3yQbRYRFIayAxHy1ouRRleTSmFwmP53fsddv41x13Kf0xdxezzPz7uv22+XPvvZdi7PixW3xOmVV6TPf156/eulX/ol6TvfkX784+HZ1kmt1sWZ2sePS6++Kv3BH7g/wne+E/D6MXQVgTkQTVh3mrQsKalX34Slo3HBst93t5JbYeNrudGokYcHH2x3A+gUcUucJPdG+dGPpDe8QbrkEjekvXWr+93HPuZut24d3WqNWrzR//3oR90fdetW6eBB99itmhtAGzCUHZCqajLnNTk5fNRw2PUxbseksuZKi5RYTgqWy5ZJP/mT/obpR408SOVPawRT9z16Uw0G3b17F4aYP/7xhVbr009Ll14qPfyw9OST0hVXjG61Jv3fr37VtXyBwBCYA1PmPHZRSTstPfCAG2IdvNBbW12yWNEGQJXD9En5BKtWudfxpS/F/1/f89s+G0u5gn3aoDsxIV14ofv3RRdla7UW+b9ADdiPGZlFF+CoV792rRsRXBysr73WZRXHBaDxcTc07KMR4mPv4elpN7dc9rFKyXs8S26Z1O//vtvacvCcnjjhVgFt3Jg/rrzwgpu6jaZzB/V67m/60EPZe9KF95Je/KYicKJlsuzHTGBGIUlBcWzMZRzHvcWMcctMb7qp+DH4CKpJwTJtcM8iCmTDgvNgkHzsMenee905m5vLEfQGnvMtb4kPypIbrTHGDd1nCa4+GkZA22UJzCR/oZCkuVlr3dBsHJ9LwHwMQ1edfDcxIf3e7w0/P8ePu6B82WWuLsbcnHst1roGSFxtjSRRctuwoBzdZ25uIfkt7fO0dZMXoC4EZhSSFBSPHh3+/3wuAcubLb5Y1Zs+fO1rw89R1KDwFfTSbJgxzKjnafsyOqBqJH+hkKREpvFx6SMfcWtz4+YeffVC82SLD1NV8l2/74anh1mzxjUokoLe7Ky0f7/796gp2VFbYyYZFVzrLI4DtBE95gYKaV3wli3JFbKuv778XmgT1oAvduCAm8sdxloXbJNGAyS3/DdNvfJRj7NypWsMxBkVXJtQHAdoEpK/GqZw9msJx3Ls2Klzl6tXu7nTqo+pSYm9u3e7TUCGffy2bZP+8A9HZ3APSkq0SnqcNWtcj/qnfip/AldI70sgRFmSvxjKbpCQNpGIO5aIMW7PgSpqNwRTKCOjUVMAGze6fy9eOx53/0hSvfKkNeiPPCKddVby70ed09CL4wBNQmBukJA2kRhVKevzny//WKqsKuZblnnxwaC3f78bvo4zai54VPDMGlzjGkVxf/OmNp6AuhCYGySk7Ne6jyWk0YM8RvVgFx/7YFLawYP5E61GJbelTX5L2yhqcuMJqAuBOXCDvY3vfS+c7Ne6M3FDGj3IK8/wr88M9LzSNoqa3ngC6kJgDtji3ka0DW2cqrNfywwQaYY+6+6xR4oO02ZdnpW1p+37eKX0jaI2NJ6AOhCYAxXX24gCkTHuYlzWuuA0igaIYdIOffrssQ8LVqOCWF3DtHkTrXwdb9pGUSiNJ6BpCMyBSuptrFkjvfe9bvlLndmvvjNxswx9+uqxDwtWd9yxdBOJwSAWyjBt2tWOPo83baOo7ukOoLGstZV/rV+/3iLZrl3WGmOtu/Se+mWMtbt3132E/k1NWTs+Hv+ax8fd7wcdPGhtr+d+Z4y77fXcz9OYm3P3H3aO437e61nb72c/Vp/yvG6fx9vvDz9v0fnJcj+gCyQdsiljJD3mQHWxt5F16LNoj33UBhxxornRuoZp8/Z8fR5v2mmMsqY7gLYjMAcqhOzbquVpjBSpbZ2nfnQUxOpqOOVNqEpzvFkSw9I2iig8AmRHYA5UF3sbVTdGkoLVMFEQq6vhlLfnO+p41651eypnSQxL2yiqamMQoC0IzAHrWm+j6sZIUrAyJn44Owq6dTWc8vbUk473gQekSy4pN5GN6l9AemxikQEXl2pUuRnFsM0X4rKy4zZlqHrjjKTNKNJsNhF3vPff7zbNGFa3e9++Yj1eNrgAsm1iQWBOiYtLew0LrqHuVuX7vZi005UxbrvOm2469edpG6n9vhsiz9uQANqC3aU8C2XNKsoxbA401LlR31McWYfHsxQqofoXkB2BOQUuLt3RlOkKn42GLIlsWRupVP8CsltW9wE0AReXbpiZccOu27a5od1t29x87sxM3UdWrigxrNdzc8rGuNteb2kiW5pG6qCoNx6nrevxgaLoMafQxWIfXeNruqIpPe7F0g6PZ22kdnE9PlAUgTkFLi7+hBq4fExXNH3v4TTD41kbqV1cjw8URWBOgYuLH6EFrsFGwv/9v8MLjaSZruhKgmCeRmrX1uMDRRGYU+LiUkxogWtxIyFp1WCa6YquJAjmbaSGmuEOhIjAnAEXl/yKBi6fQ+BxjYQkaaYrupQgSCMVKBeBGZUoErh8D4EnNRIWW7Mm3XRF1xIEFzdS+31pejq83AGgiQjMqETewFXGEHiWXaVefTX+mBb33rucIFi04RRqQiBQF0pyohJ5azxPT/uv45z0mHEGjy+pHKbUvbKtRUtuVl3qlkYA6kJJTgQnb9JQGXO3Sb3bONEc+GWXje69d23utUjuQNUJgaGtCgCGITCjMnmShsqYu41rJCQNHEUNgLRBqEsJgkUaTlVmsoe2KgBIQmBGpbJmtpc1d7u4kfDcc9If/7F09OjS+0YNgKoyr5s03Fqk4VRlJntSI6Dfl+67T3rf+/w9H1AEgRlBK7O4y2AjYX5e+tM/jQ/MUQPg/vvLz7xu2nBrkYZTlZnsSY2Ao0elq6+W3vSmMM9xpZrUKmwxkr+QWR2f3Sr2Rh6ViDQ/L511lvTSS0v/r4+9hZu6d3HeBK68CYF5pEn4C/kcVyL6Q772tdL+/dLWrdILL4TbKmyYLMlfstZW/rV+/XqLZjp40Npez9rxcWuNcbe9nvt5G8zNWTs1Ze3u3e6231/43cGD1p5xhrVuRnrha3zcz+ufmnKPtfjxo+eYmir+HGVJOm9Jqno/9fvucePObVPOcanm5hZO0J13Wnv8uLuV3M/T/kExlKRDNmWMZCgbqXUhgWbYHHj02l9+eenvli+X3vzm4s/d5OpheaviVVVFLJoS2bjR9dTjhH6OT+F72OrAgaUf4GuucV/z8+2pKdsQhfZjNsb8R2PMk8aYE8aYdF10jBRVUdq9290Ou5BULetevG1SxWvv6t7FUVC/6SZ3W1bjbmJC+uQnhz9+Y85xtHH4Lbe4lsbNNxffOPyZZ9yLv+ceaW7O/WxuTrr7bun88xvUYmmHoj3mb0ialLTPw7FAYSf/JPXoZmdd8Ao1V6RoB6OM3uziY3rHO1zvO87y5dKWLaPLXpK7k+zd75Y+/OH4xm4jKrQNDlvt2SNdfLF03XXuIlFk2GrdOjcc9Oij0nve497UK1e672dnG9JiaZG0Y95JX5K+KmlD2vszxxxvcJpn8VcI0zxJc6CStWNjYc43+5jH9D3/O+yY/uiPsv188DW0ff7fl0afp6QJ/H4//yR5NAn/6KPWHjtm7Y4d7vbLXw7j4tMCYo65mULfOnBUxaz5efcV0nxz3nnxrL3ZLD2tpGP60Idc7/uhhxbmXLdscbdJr8Ha9s//+9Lo3bGiIec9e6TLL3e7rMzNuYXYN96YraTdoGgS/umnpUsvlR5+WHrySemKK9h0vgYjA7Mx5kuSzor51Y3W2s+mfSJjzFWSrpKktWvXpj7ALgk9+WdwTXG/H7/mVwqjERHJ09gZtvznjjtc4Cy6nnrUMT300KnHND2dbn475EZdaBq7hWuZQ84TE9KFF7p/X3RRw1os7TIyMFtr3+bjiay1d0m6S3LrmH08Zts0YevAqLfx3vdKn/lM/H1CaEREsjZ2svZm81y3sh5TmvtbG3ajLg3mx1OIhq22bnXff+xj0t697vsHHyw+Sd7YFku7MJQdkKZsHbh6tdvQ4QtfCLsRIWVv7GTtzVZxTEn3X7HCDVNPTITfqEsSctJjUBhy7oRClb+MMZdJ+oSk10l6UdIT1tpfHPX/qPw1XNXb4OVVZdWmIrIe5+7d0s6d8ZtaGCPt2uWW9VR5TEn3l9z2l6edJp04Eb/OOqS/R5ymVjyrVRWl8OBVZds+WmsPSDpQ5DFwqqYkppRZw9qnrMdZxXTC2Jj0wAPS5s3umI4dc4+9YkX8MQ2+hldeWTpkHR3r+LgLZCH/PeKEnvQYJIacW42h7AA15TPXlEZEluOsYjphZsaNQi5b5oLyypWuh/65zw0fFYlew44d0qc+5f5fnN/7Pfd4dfw98s4Rh570CFSNwIxCmtKISHucZY8ExCWXHTvmvi65JHnYdvVqt4nGK6/E//7IEen554sPtedRZI64CUmPQJUIzMCAft/l1Xzwg9IPfyideab0xjf663kWHbYNMYgVraHelKRHoCoEZuCkpMQ7X8PBRYdtQwxiRRsbTclXAKpCYAZU3c5ZRXu8IQYxH3PETclXAKpAYAZUXWawjx5vaEHM1/B6U/IVgLIRmAFVlxnsq8cbUhALcXgdaDICM6Bqk6pC6/EWFeLwOtBkhSp/5UXlL4SmKZXMQkYxKmC4yip/AW1Br6+4kIbXgSYjMCdgt5tuadsQM4BmIjAPwW433USvr340iNF1zDHHYLcboB5N2V0NyCrLHPOysg+midKsaQXg12CRl9lZt7HH7Kz7ftMml6AHdAGBOQa73bRPvy9NT7v9lqenuciHiAYx4DDHHKPujQKYY/OLfIF8qn4f0iAGHAJzjDorGRFE/KqqBnbc8za5cVXH+9B3g7jpfwN0F8lfQ9SRhELSmX/T09K2bfEX+/Fxad++dFnYWS7yTU9gqut96LPIS9P/Bmgfkr88iNa07tsn7drlbg8fLvdDzRybfz6GR2dmXKDatk3audPdnn22+/libUhgqut9GBV56fVco8kYd7tmjXT55e5505y/NvwN0G0E5gTRmtabbnK3ZfdWmWPzLxoejZNmeDTrRb7KoFZWQlud78PBBvFVV0knTrgAfdddyQ2iQTRw0XQE5oAUDSJYanLSDWHGSZMvkPUiX1VQy9KLz6ru9+Hq1dJll0mf/rQ0N+fOW5ZeLw1cNB2BOSBFgwiWGjY82uulq4Gd9SJfRVAre6g2hPdhkV5v3Q0LoCgCc0CKBhHEK5IvkPUiX0VQK3uoNoT3YZFebwgNC6AIlksFpu6NFNq6xCRvDeysS+eq2KWqiqHaut+HRZZOsVMYmo7lUvhnLDGJl+e8lLk3sa8lYCHzsXSK/aERkizLpQjMkMQa6lFCushnCVpNHgGhoYg2yRKYGcqGpHTzlk3vhRUR0naQaYdqm15Fru7hdKAuBGZIYolJ04wKWnWVIvUtpAZRFk0eqUD9CMyQVP/GHcguKWgxAlKfpo9UoH4sl4Iklpi0DSMg9aAcKHwgMENSGGtX4Q9FNupBOVD4wFA2/hnJNu1R59alXcZIBXwgMOMUTU22wakoslEPcjXgA+uYgRYLaf11F/jcUxrtwjpmAJIYAala1pEKllUhDoEZADxKm6vBsioMQ2AGAkaPqplGjVS0pQAMykFgBgJFj6q9KACDJARmBKvLvUV6VO3GsiokITAjSF3vLdKjajeWVSEJlb8QHMoa0qNqO0rgIgmBGcGhrCElNduOErhIwlA2gkNvkZKaVakzj4ESuBiGwIzgMP9GSc0qhJDHQAEYxKEkJ4JDWcMFlNQsR78vnXMO7zFUh5KcaDR6iwvoUZWDrHeEjMCMIDH/hjKRx4CQEZgRrKb2FrtcGKUpyGNAyJhjBjyKSyiKhuC7UBilKchjQNWyzDGzjhnwhMIozcE6YoSMoWzAExKKmoU8BoSKwAx4QkJR8zQ1jwHtxlA24AllNAH4QGAGPGFjAgA+EJgBT0goAuADc8yARyQUASiqUGA2xtwq6RJJxyQ9I+n91tqYlYHosq4V3CChCEARRYeyvyjpQmvtz0j6e0kfLX5IaJOZGbdZwLZt0s6d7vbss93PAQBLFQrM1tpHrLXRys2vSTq3+CGhLSi4AQDZ+Uz+2irpoWG/NMZcZYw5ZIw59IMf/MDj0yJUaQpuAIP6fWl6Wtq9293SeEMXjZxjNsZ8SdJZMb+60Vr72ZP3uVHScUlTwx7HWnuXpLskVys719GiUSi4gSzi6oxv306dcXTPyMBsrX1b0u+NMe+T9MuS/r2tY0cMBIsdfJDW4LRHJHrfbNrEphLolkJD2caYzZI+IulSa+2cn0NCW1BwA2kx7QEsKDrH/N8knSHpi8aYJ4wxf+jhmNASFNxAWkx7AAsKrWO21v4rXweCdqLgBtJg2gNYYOqYFt6wYYM9dOhQ5c8LIEzz8259+4sx5Yl6PeaY0XzGmMettRvS3Jda2QBqx7QHsIBa2QCCwLQH4BCYA9W1+tKARJ1xQCIwB4lCCwDQXQTmwFBoAQC6jeSvwFBoAQC6jcAcGAotAEC3EZgDExVaiEOhBQBoPwJzYMqsL82WegAQPpK/AhMVWliclb18ebFCC2R6A0AzEJgD5LvQApneANAcBOZA+Sy0kCbTm6IOABAG5pg7gExvAGgOAnMHkOkNAM1BYO6AMjO9AQB+EZg7gC31AKA5SP7qCLbUA4BmIDB3SFlb6rFFJQD4Q2BGIRQuAQC/CMzIjcIlAOAfyV/IjS0qAcA/AjNyo3AJAPhHYEZuFC4BAP8IzMiNwiUA4B+BGblRuAQA/CMrG4VQuAQA/CIwo7CyCpcAQBcxlA0AQEAIzAAABITADABAQAjMAAAEhMAMAEBACMwAAASEwAwAQEAIzAAABITADABAQAjMAAAEhMAMAEBACMwAAASEwAwAQEAIzAAABITADABAQAjMAAAEhMAMAEBACMwAAASEwAwAQECW130AQF36fenAAemZZ6R166TJSWlsrO6jAtB1BGZ00syMtGmTdPy4dOSIdPrp0vbt0iOPSBMTdR8dgC4jMKNz+n0XlF98ceFns7PudtMm6fBhes4A6sMcMzrnwAHXU45z/Lh0//3VHg8ADCIwo3OeecYNX8c5ckR69tlqjwcABhGY0Tnr1rk55Tinny5dcEG1xwMAgwjM6JzJSWn5kOyK5cvd7wGgLgRmdM7YmMu+7vWk8XHJGHfb67mfk/gFoE5kZaOTJiak737XJYI9+6wbvmYdM4AQFArMxpjdkt4p6YSk70t6n7X2uz4ODCjb6tXSr/963UcBAKcqOpR9q7X2Z6y1PyvpzyTd4uGYAADorEKB2Vr744FvT5dkix0OAADdVniO2RizR9J7Jb0k6RcS7neVpKskae3atUWfFgCAVjLWJndyjTFfknRWzK9utNZ+duB+H5U0Zq3dOepJN2zYYA8dOpT1WAEAaCRjzOPW2g1p7juyx2ytfVvK552S9HlJIwMzAACIV2iO2RjzhoFv3ynpW8UOBwCAbis6x/xfjDFvlFsu9W1JHyh+SAAAdFehwGytfZevAwEAAJTkBAAgKARmAAACQmAGACAgBGYAAAJCYAYAICAEZgAAAjKyJGcpT2rMD+TWPTfRmZL+qe6DCBznKB3OUzqcp9E4R+nUeZ7+pbX2dWnuWEtgbjJjzKG09U67inOUDucpHc7TaJyjdJpynhjKBgAgIARmAAACQmDO7q66D6ABOEfpcJ7S4TyNxjlKpxHniTlmAAACQo8ZAICAEJgBAAgIgTkjY8ytxphvGWP+1hhzwBjTq/uYQmSM+Y/GmCeNMSeMMcEvT6iSMWazMeYpY8w/GGP+c93HEypjzN3GmO8bY75R94ujv34AAAJ/SURBVLGEyhjzemPMV4wxf3fy8/af6j6mEBljxowxB40xf3PyPP123ceUhMCc3RclXWit/RlJfy/pozUfT6i+IWlS0mN1H0hIjDGnSbpT0hZJPy3p14wxP13vUQXrjyRtrvsgAndc0nXW2p+W9G8lbef9FOuopLdaa/+NpJ+VtNkY829rPqahCMwZWWsfsdYeP/nt1ySdW+fxhMpa+01r7VN1H0eAfl7SP1hrn7XWHpN0n6R31nxMQbLWPibphbqPI2TW2sPW2r8++e+XJX1T0k/Ue1Thsc7syW9XnPwKNvOZwFzMVkkP1X0QaJSfkPTcwPfPiwspPDDGnCfp5yT9Vb1HEiZjzGnGmCckfV/SF621wZ6n5XUfQIiMMV+SdFbMr2601n725H1ulBtGmqry2EKS5jwBKJ8xZlzS/5L0IWvtj+s+nhBZa1+V9LMn84IOGGMutNYGmb9AYI5hrX1b0u+NMe+T9MuS/r3t8ELwUecJsb4j6fUD35978mdALsaYFXJBecpae3/dxxM6a+2LxpivyOUvBBmYGcrOyBizWdJHJF1qrZ2r+3jQODOS3mCMOd8Ys1LSuyU9UPMxoaGMMUbSfknftNbeXvfxhMoY87poBY0xZrWkt0v6Vr1HNRyBObv/JukMSV80xjxhjPnDug8oRMaYy4wxz0v6d5IeNMY8XPcxheBk4uBvSXpYLlHnT621T9Z7VGEyxvxPSf9H0huNMc8bY36z7mMK0EWSrpT01pPXoyeMMe+o+6ACdLakrxhj/laucfxFa+2f1XxMQ1GSEwCAgNBjBgAgIARmAAACQmAGACAgBGYAAAJCYAYAICAEZgAAAkJgBgAgIP8f+2AZnwXovLYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.io\n",
    "import matplotlib.pyplot as plt\n",
    "mat = scipy.io.loadmat('hw5_p1a.mat')\n",
    "print (mat.keys())\n",
    "X = mat['X']\n",
    "k = 2\n",
    "num_points = int(X.size/2)\n",
    "cluster_assignment = np.zeros(num_points, np.int8)\n",
    "# Guess the mus\n",
    "changed_points = 1\n",
    "num_dims = 2\n",
    "\n",
    "mus_choice = np.random.choice(num_points, k)\n",
    "\n",
    "mus = []\n",
    "for i in range(k):\n",
    "    this_mu = X[mus_choice[i]]\n",
    "    x_val = this_mu[0]\n",
    "    y_val = this_mu[1]\n",
    "    mus.append(this_mu)\n",
    "\n",
    "def get_distance(point, mu):\n",
    "    # We can skip the sqrt\n",
    "    return np.dot(point-mu, point-mu)\n",
    "\n",
    "\n",
    "print(get_distance(np.asarray([0, 1]), np.asarray([0,3])))\n",
    "\n",
    "steps = 1\n",
    "while(changed_points != 0):\n",
    "    changed_points = 0\n",
    "    \n",
    "    # ASSIGN POINTS TO CLUSTERS\n",
    "    for point_counter, point in enumerate(X):\n",
    "        distance_vector = np.zeros(k)\n",
    "        for counter, mu in enumerate(mus):\n",
    "            distance_vector[counter] = get_distance(point, mu)\n",
    "        new_assignment = np.argmin(distance_vector)\n",
    "        old_assignment = cluster_assignment[point_counter]\n",
    "        if new_assignment != old_assignment:\n",
    "            changed_points = changed_points + 1\n",
    "            cluster_assignment[point_counter] = new_assignment\n",
    "    \n",
    "    # UPDATE MUS\n",
    "    for counter, mu in enumerate(mus):\n",
    "        my_points = X[cluster_assignment==counter]\n",
    "        mus[counter] = np.sum(my_points, axis=0)/num_points\n",
    "    if steps == 1:\n",
    "        cluster_assignment_legacy = np.copy(cluster_assignment)\n",
    "    steps = steps + 1\n",
    "        \n",
    "red_X = X[cluster_assignment==0]\n",
    "blue_X = X[cluster_assignment==1]\n",
    "\n",
    "plt.figure(figsize = (8, 8))\n",
    "plt.plot(red_X[:,0], red_X[:,1], 'ro', ms = 7)\n",
    "plt.plot(blue_X[:,0], blue_X[:,1], 'bo', ms = 7)\n",
    "\n",
    "not_changed = (cluster_assignment == cluster_assignment_legacy)\n",
    "star_X = X[not_changed==0]\n",
    "\n",
    "plt.plot(star_X[:,0], star_X[:,1], 'w*', ms = 7)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['__header__', '__version__', '__globals__', 'X'])\n",
      "4\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10d571898>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAHVCAYAAADo5GInAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+0HXV97//XmxySEE6vB4RqAkTAUtRbv7fqgbbfehF/RfwBSEotuqrQyiW16F0WqeKFEldybale1PZ++7VhAYrWCN58oYbKj6sixbVuKeewlog/QCDWkhgFQVie5JyEhPf3j88eM9mZmT2zZ/bs2bOfj7XO2mfvPfuc2WfPmfd8Pp/35/0xdxcAABh9Bw17BwAAQDUI6gAAtARBHQCAliCoAwDQEgR1AABagqAOAEBLENQBAGgJgjoAAC1BUAcAoCUmhr0D/TjiiCP82GOPHfZuAABQi3vvvfdn7n5kr+1GMqgfe+yxmp2dHfZuAABQCzP7UZ7t6H4HAKAlCOoAALQEQR0AgJYgqAMA0BIEdQAAWoKgDgBASxDUAQBoCYI6AAAtQVAHAKAlKgnqZnatmT1mZt9Jed7M7G/N7GEz+7aZvTz23Llm9lDn69wq9gcAgHFUVUv9s5JOy3j+jZJO6HxdIOnTkmRmh0taK+m3JJ0saa2ZHVbRPgEAMFYqCerufpekJzM2OVPS5zy4W9KUmS2X9AZJX3X3J93955K+quyLAwAAkKKuMfWjJD0au7+181ja4wcwswvMbNbMZh9//PGB7SgAAKNqZBLl3P0qd5929+kjj+y5+hwAAGOnrqVXt0k6Jnb/6M5j2ySd2vX4nTXtEzA08/PSTTdJjzwivfCF0urV0tKlw94rAKOurqC+WdJ7zex6haS4p919u5ndLukvY8lxqyR9uKZ9AoZiZkZatUras0fasUM69FDpwgul//2/pZNOGvbeARhllQR1M/uiQov7CDPbqpDRfrAkufvfS7pF0pskPSxpp6Q/6jz3pJmtlzTT+VHr3D0r4Q4YafPzIaA/9dS+x+bmwu2qVdL27bTYq0avCMZJJUHd3d/e43mXdGHKc9dKuraK/QCa7qabQgs9yZ490o03Su94R7371Gb0imDc1NX9DkChtbhjR/JzO3ZIW7bUuz9tRq8IxtHIZL8DbfDCF4bWYpJDD5WOP77e/WmzPL0iQNsQ1IEarV4tTaT0j01MhOdRjaxekbm5EPQXFurdJ2DQCOpAjZYuDeO5U1PS5KRkFm6npsLjdAdXJ6tXRJL+6Z+k5cvDuDvQFhZy2EbL9PS0z87ODns3gL5FGdlbtoQu93HLyO4nI73oaxYWQtCOj6knmZpifB3NZ2b3uvt0z+0I6gDqlJSRPjGRnZHez2vir5ufl3btSt5mclLasIFZB2i2vEGd7ncAAzc/L23cKF1+uXTqqaH1PDcnuYfbp54KwTdpjDuexZ73NZGTTpJ+/GPp9NPTt2HWAdqEoA5goGZmpBUrpDVrpPXrpZ07k7dLy0gvm8V+yCHSWWeFFnkSZh2gTQjqAAamu5WdJa3FXMXcfmYdYFwQ1AEMTFYru1tai7mKuf3MOsC4oKIcgIHJamV3S2sxr14dSrsWeU2SaHx9nGcdoP0I6kAWVgMpJWplZ3W9T07uy2RP+tNGrey07PciH8chhzQjy53DCoPClDYgTb/zqPBLWXPFly2TLr5YOvHEYvPUR72VXeSwIvgjwjx1oIz5+ZCynRSNqFZSSFuujaoIsEUOq7b83VCNvEGd7ncgSVaG186d0gc+IF15JYE9h0GPZdfRmi27hGu0jzfdFL5PEl96lxXm0C+COpAkK8Nr927p6qtDNRWaTbkMaiy7jvXS8wbYtIuL+D5m5RbEp+flmZvfhNwANA9BHUjSK8Nr9+7wRbNpaOpqzebptHn720PVuu6Li82bpTPO6F1/Xtp/el4Vc/MxnpinDiTJqlYSx8LcQ1PXeul5Om1OOSW5jO1pp+Wfpx+fnlfF3HyMJ4I6kCRerWTx4vTtaDYNTV2t2V5LuO7eHQJ5kl5d7lJyERwq4KFfBHUgTZThdf756YGdZtPQ1NWazdtpk2T37uzgfPLJ0v/8n2GoIJ4DQAU89IspbUAvWZOtmd42NHV+LFGy286dIVBX5dBDpYMPTk/sa8vcfJTHlDagKlWWNENl6vxYok6biy8OY+hVBfZo+CAtsa8pFfAwOmipA1ni85SOOSY8tnUrzaYGqbM1m9U7YBaq5O3cKS1Zkr3Oe7fJSWnDBgI40tFSB8qipNdIqLM1m9U7sHmz9Oij4eLivvukTZvy/1zyLVEVgjqQhJJeSJGnQt7GjdJtt/XOfI+Qb4mqENSBJJT0QoZevQNZy8UmYZoaqkJQB5JQ0gslpHXTS2Hs3X34+ZasANdOBHUgSVaZWPpKkUNaN7378Kep1VEzH8NB9juQhLnpaClWFR5NebPfqSgHJKGkFwZgfj4k0a1fH26LTHurSl018zEcdL9j+Jo6uDfohcAxVprS5U26SLsR1DFcTTnTpaGkFyrQpBmSpIu0G93vGJ74ma57zcpVq4bTNwkMQJO6vFkBrt0I6hieJp3pgAHK2+Vdx5g76SLtRvc7hofBPYyJPF3edY5EkS7SXgR15Fd1QhuDexgTWRXmJiakN74xHO51jrmTLtJOdL8jn5mZMLl1zRpp7dpwu3x5eLxfDO5hTPTq8r71VkaiUA1a6uhtUKm7rFOOMZLV5X3bbaMxEtXU2afYh6CO3ga5uAmDexgjaV3egxqJqjIIN332KQKCOnobdEIbg3sYc73G3PsZiaoyCDdpnj2yMaaO3qJmRBIS2oDSqp5mVnUJCGafjg6COnojoQ0YuGgkasMGad26cLt9e39d21UHYWafjo5Kut/N7DRJfyNpkaSr3f2Kruc/KenVnbvLJP2qu091ntsr6f7Oc//u7mdUsU+oEAltQC2qGomqOggz+3R0lA7qZrZI0t9Jer2krZJmzGyzu38v2sbd/yy2/fskvSz2I+bd/TfL7gcGrOqENtJogYGpOggPYswfg1F6PXUz+x1JH3H3N3Tuf1iS3P2vUrb/P5LWuvtXO/fn3H2yyO9kPfURl5TBE7X6SaMFSltYCGUkqlwzPenfdtEi6QMfCM9zbT5YeddTryKony3pNHc/v3P/nZJ+y93fm7DtCyTdLelod9/beWyPpG9J2iPpCnf/x5Tfc4GkCyRp5cqVr/jRj35Uar8xJPPzoYhN0tlm2TJp27Zw1gFQyiCunaMOti1bQvLdJz7BtXld8gb1uqe0nSNpUxTQO17g7tvM7HhJd5jZ/e7+SPcL3f0qSVdJoaVez+6iclkZPDt3SkcdJd15J2cFoKRBlICIxvyTrs2Z4tYMVQT1bZKOid0/uvNYknMk7Tcy4+7bOrdbzOxOhfH2A4I6WiIrg0cKgZ2zAlCJQZWAGGQ9KpRTxZS2GUknmNlxZrZYIXBv7t7IzF4k6TBJ/xJ77DAzW9L5/ghJvyvpe92vRYtkzXmPMPEVaDSmuDVX6aDu7nskvVfS7ZK+L+lL7v5dM1tnZvHpaedIut73H8R/saRZM7tP0jcUxtQJ6m2WNec9wlkBaIykNd6pR9VclYypu/stkm7peuzyrvsfSXjd/5H00ir2ASMimvN+6qmhqz0JZwWgEdJKzd58M/WomoqKcqjfSSdJW7eGbPcknBWAocsqNXv66SGwV1XWFtUhqGM4DjssZLlzVgAaqVcy3L//exhbf9e7pFe/Otz+8IdMXBk2VmnD8LDsKtBYvZLh7rordMVHXfP33BPG3JmnPlwEdRRXZYlXll0FGimr1OyyZdLnP79/Wgzz1JuB7ncUMzMTqk6sWSOtXRtuly8PjwNojayJKu7SQSnRgxmpw0VQR35VL9IMoLGy1nh/5zuZp95UBHXkd9NN0jPPJD/H5TnQOmlrvJ9yCvPUm4ox9TYa1LKmd93F5TkwZpLSXliKtbkI6m2TVi2ibErq/HzIjEmzbBmX58CYiLrm01aBI0lueAjqbRIf845UlZJ6003pmTFSGF/n8hwYG8xIbSaCepsMcumkXqurvfOd/DcDY6bMjNRBjRKOO4J6mwxy6aSsSauTkyFzBgByGNQoIch+b5dBLp2UNWmVzBgAOTEzdrAI6m0yyMCbNWmVzBgAOeUZJUT/6H5vk0GnpJIZA6CkQY4SgqDePoMOvNRqB1BCVnoOhWvKM3cf9j4UNj097bOzs8PeDQBAQQsLYbmI+MzbyNQUi8GkMbN73X2613aMqQMAakN6zmDR/Q4AqBXpOYNDUAcA1I70nMGg+x0AgJYgqAMA0BIEdQAAWoKgDgBASxDUAQBoCbLfAQBDx1Ks1SCoAwCGiqVYq0NQBwAMTXwp1khUF37Vqv3LxtKa742gDgAYmjxLsb7jHbTm8yKoAwCGJs9SrEVa8+OO7Pe2mp+XNm6U1q8PtwsLw94jADhAtBRrkmgp1jyteQQE9TaamZFWrJDWrJHWrg23y5eHxwGgQVavliZS+ownJsLzeVrzCAjqbRPvp5qbk9zD7VNPhcdpsQNokDxLseZpzSMgqLcN/VQARky0FOuGDdK6deF2+/Z9CXB5WvMISJRrG/qpAIygrKVYo9Z8d/b7xMS+1jwCgnrbRP1UUWpoHP1UAEZU1Jq/6abQNjn+eOapJzF3H/Y+FDY9Pe2zs7PD3o1mWlgISXHxuR+RqSnmfgDACDKze919utd2jKm3TZ6sEwBAK9H93kb0UwHAWCKot1VW1gkAoJXofgcAoCVoqQMARhqrt+1TSUvdzE4zswfN7GEzuyTh+fPM7HEz+1bn6/zYc+ea2UOdr3Or2B8AwHigKvb+SrfUzWyRpL+T9HpJWyXNmNlmd/9e16Y3uPt7u157uKS1kqYluaR7O6/9edn9AgC0G6u3HaiKlvrJkh529y3uvlvS9ZLOzPnaN0j6qrs/2QnkX5V0WgX7BABoOapiH6iKoH6UpEdj97d2Huv2e2b2bTPbZGbHFHytzOwCM5s1s9nHH3+8gt0GAIwyqmIfqK7s95slHevu/5dCa/y6oj/A3a9y92l3nz7yyCMr38HGY310ANhPr9Xbtm8fv1NmFdnv2yQdE7t/dOexX3L3J2J3r5b0sdhrT+167Z0V7FO7zMwcuJLBhReGCnHRMkYAMGZWrw6nwiQ7dkjXXSft3Dlep8wqWuozkk4ws+PMbLGkcyRtjm9gZstjd8+Q9P3O97dLWmVmh5nZYZJWdR5DhPXRASBRWlVss3Cq3LFj/E6ZpVvq7r7HzN6rEIwXSbrW3b9rZuskzbr7Zkn/1czOkLRH0pOSzuu89kkzW69wYSBJ69z9ybL7NDLyTK7MkwlC5TgAY6q7Kvb27aGFnjTWPg6nzEqKz7j7LZJu6Xrs8tj3H5b04ZTXXivp2ir2Y6Tk7VInEwQAMsWrYq9fH7rck4zDKZMyscNQpEs9KxPk4IPDZWnb+5MAIKdeyXPHH1/v/tSNoD4MRSZXrl4tTaR0qOzeLX3uc+NdPgkAYrJOmRMT4fk2I6gPQ5Eu9XgmSNLl5zhlgABAD2nJc1NT4fG2V5gjqA9D0f6hKBPk3HOlxYuTXzeu5ZMAoEt0ytywQVq3Ltxu397+6WwSq7QNR9bkyrT+oUMOkZ7/fOmZZ5JfNw4ZIACQUzx5bpzQUh+GfvuHxj0DBACQiZb6sHRPrjz++APnqXfPY3/Tm8Y7AwQAkImgPkxZ/UNJ89gnJqRPfUp6//sPfHwcMkAAAJkI6k2UtUjw+98fWva33prewgcAjCWCehP1msd+663jmQECAMhEUK9CnhruRV5DaVgAyK2fU3BbEdTL6mdZ1F6vibLcoy73OLLcAeCXWJl6f+buw96Hwqanp312dnbYuxEuD1es2H/sOzI1FaoddF8u5nmNFEq/Fvm5ADBm+jkFjyozu9fdp3ttxzz1MorUcC/ymnGvcwgAOfRzCm47ut/L6GfsO+9r8sxjB4Ax1m/6UZvH4AnqZfQz9l3kNeNa5xAAcshzOu0O4CtXSqef3t4xeMbUy1hYKD723c9rAAAH6HU6vfnm/QP4smXSzp1SUthr+umXMfU69DP2zXg5AFQi63S6eXMI6E89FVry7iGwp7Vj2zIGT/d7Wf2MfTNeDgCVSDud3nhjehJdkraUACGoV6GfsW/GywGgEkmn06wkuiRtKQFC9zsAoHWyVqpO0paFLgnqAIDWWb06faXqaOy9jSlNdL8DAFonSqJLWsF682bp0UfbmdJEUAcAtNI45iQT1JugzeWNAGCIxi0nmaA+bCwxBACoCEF9mObnQ0CPl0OK6h2uWtXs8kYAgMYh+32YWGIIAFAhgvow9bvEEAAACeh+H6Z+VnkDBoicTWC0sUrbMLFiGxokKWdzYoKcTaAJWKVtFLBiGxoinrMZrWg1Nxfur1oVrj8BNB/d78M2jtUR0Dh5cjbHaa4vMKoI6k0wbtURMHTdY+cPPEDOJtAGBPW6kIGEhkgaO3/22XA4zs8fuD05m8DoIKjXgapxaIisekdmya/pd0lKrmOB+hHUB42qcWiQrLHzQw4Jt2bh2nPx4hDQb765+CHKdSwwHGS/DxpV4zAI8/PSxo3S+vXhNmd6ela9o/l56e1vD5nvixdLu3dLBx0knX56CNJFdi0tk/6UU6TPfpZsemBQCOqD9sgjycVlpPA4GUgoamZGWrFCWrNGWrs23C5fnivyRvWOkhx6qPTFL0o7d4aALvU3rS3rOnZhQXrPe3LvLoCCCOqDdswx2c8ffXQ9+4F2KDmh/I1vDElxSZ59Nn1cvUinUlZvgBR2kfnvwGAQ1IFRUmI4Z2YmZLGnFZHcsye00pMUmdaW1RvQ/fsYfQKqRaJclZLSfR99NDR/0s6kd98tnXderbuJEdbnIkBJ+Zrdoi73JHmmtUWH/4MPpvcG5NxdAH2qJKib2WmS/kbSIklXu/sVXc9fJOl8SXskPS7pj939R53n9kq6v7Ppv7v7GVXsU+3S0n0vukhatiz9RPz5z0uf+hQZ8Minz0WAshr4efSa1tZ9+B9ySLiWnZiQnnkm+TW9LhSYEgcUV3pBFzNbJOkHkl4vaaukGUlvd/fvxbZ5taR/dfedZvYeSae6+x90nptz98kiv7NxC7rMz4fEpaRm0HOeE85qaf2ak5PShg1UlEM+fSwCND8vvetd0qZN+X7F0qUhGOdd1CXr8F+2TNq7V9q1K/fuSmJxGaBbnQu6nCzpYXff4u67JV0v6cz4Bu7+DXePotrdktqVHZbVDNq7V/rt305/LX2QKKLgIkBRovyXv5zvx09OSp/+dLjOXLcu3G7ZIj30UPrsuazD/6CDpEsvLbZmEYvLAP2rovv9KEmPxu5vlfRbGdu/W9KtsftLzWxWoWv+Cnf/x6QXmdkFki6QpJUrV5ba4cr1Guc8/PBwJmPddFQh5yJA8/PSa18r/eIX+X/0xIR0zjn7flSUXJdVRKbX4W9WbM0iFpfBKGnaMFGtiXJm9oeSpiW9KvbwC9x9m5kdL+kOM7vf3R/pfq27XyXpKil0v9eyw3n1Gud885ulr30t+bX91uDEeMuxCNANN+QL6Gb7d29HJ6S8xRDzDPMXWbOoz1xAoHZNrJxYRff7NknxydhHdx7bj5m9TtKlks5w91+OsLn7ts7tFkl3SnpZBftUr9WrwxkxSdT0Yd101OwrX8l+fnpa+sxn9nWzb9++/4ko7+y5Xod/0WvWXgVy6NhCEzR1mKiKoD4j6QQzO87MFks6R9Lm+AZm9jJJGxQC+mOxxw8zsyWd74+Q9LuSvqdRk2ecM+oyjQ9Wdp9FgRoddFCYcZnWHf7AA+nFEOMt5oLD/D1VfZEADEJTK4CX7n539z1m9l5JtytMabvW3b9rZuskzbr7ZkkflzQp6X9ZKFkVTV17saQNZvaswgXGFfGs+ZGSZ5yTddPRS4UDdG9+c3bG+7e+FboPk7oMZ2akK69Mf213iznnMH8u0UVCWvY7HVtogqYOE5We0jYMtU1py3OCbVqWBEZD0nFz//2VzuNaWJCOPDK9td1t2TJp2zZpyZL0KWqRrOloVYn+RGUvEoBB2LgxLLuQ9P81iJnKeae0EdTT5Jkoy2Ra9CPtuNmzJ/kMUSKCzsxIr399CPC7dkkHH5xeDEYKgf2SS6SPfSz9YmDZMunOOznEMd76KBlRCkG9jKxqGtGn5d57G5oV6Pbkk2GRn7RiRElKXvbHW7z33de7CM3SpeECIOnUYCZdfrn0kY/0tStAq9TZrssb1Kn9niRvBgSTaVHEzIx06qnFArpUeoAunsqxcaN0223ZXfLuoQs+KXv30EOlX//1vncFaJUqc0mqQlBPkicDwr2ZWRJopmj+S9GALlU6j2v16pAUl2XXLmnx4uTnyD4H9te0/GeWXk2SZ6Isk2lRRJkVVSqMpFFm+bJl2du5H7i2uhlrDwFNR1BPkmeiLJNpUURW70+apUsHUqDopJOkrVvTW+NSSKbrHlN3l97/fmqvA01GUE+Sp5qGe1hWdenS8EWVOGTJ6tlJ85a3hITL3/iNMBjevaLK/Hzy4zkcdpj0R39UbHek4RbVANAbY+ppsjIg4imPCwshq2jJEumDH5T+/M8J6DhQnsHsuMlJ6ayzkueuX3hh6Ad///tDk3rHjtDs/i//Rbr9dumVr8z1K045RfqHfyjWgUC6CNBsTGkrKs90N4I6kiTNf9mxI3nu2NTUvovJpGPNLH3O2V135QrsWfNs0wyiqAaA3pjSNiisC4l+JfX+rFwpnX568kTXW29NP9bSLsbdpTe8QXriiZ4Xl2nlWBctkp59NnmFN9JFgGYjqBfV1IK/GA1J81/Shnluu614cp0UuuRzXlymjTJlVaylIwpoLoJ6UXkWjwaKSJvomnWsZXnmmdAFn7PHKOnXN7GoBoDeGFMvqu6Cvxhf/Qx6R5Yty9UFD2A05B1TZ0pbEfPz0vXXSy9/eVgZY8mS8DhT2TAIWVMrr7oq+7UHHcTcM2AM0f2e18yM9OpXHzjGuXgxU9kwOFn94P/6r9I11yS/rkX5HaxuDORH93se8/NhYeq0pKXubnfOQqhD3Qs6DwGrGwMBS69W6bOfzS6/tXRpaDG94x2chbDPoC/uWp7fQUkIYB/mqVfpK1/Jfn5hIXR1Ritxxc9CUStq1SrOQuMk6eLuwgurvbhLm2je59yzpnUwURICKI6gXoWJiTDWyVkIUvmLuyLRtaK5Z3VcgxQ17JIQTbvIAfIgqOfx5jdLmzalP794cfiP//jHKUyDchd3/UTXkgs6N7WDaRAlIfIG6iZe5GCwqrqIG/rFoLuP3NcrXvEKr9X8vPuv/Ip7KMJ54Nc3vxm2+8IX3Ccnk7eZnAzPo/3WrXM3Sz4OzNzXr09+3c6d7lNTya+bmgrHYbTdF74Qfs8XvrDv8T419bCdn8/358jrnnvC6yYnw8cwORnu33PP/tvl/RjQHnmPjbp+ThJJs54jPjJPPY+lS6Wvfz1k50Rz0ycmQoGPb35z3+IZrLEOKXuZ1awmZp4W/sxMyB5bs0ZauzbcLl8eHu/TsLu50+RZATmveG/E3FwI0XNz4f6qVfuvWpvnY0B7FDk26vg5ZRHU84rGLq+9Nqxffd11oWLXK1+5b13rj388rLFexVkIo6vfi7te0fXBBwdy1uj3GqQO0b/dhg3SunXhdvv24l3gRQJ1Uy9yMBhVXcQ15WKQMfUiksYukwbfFi0KBWnMKJo9jvrNSu81iPyznw0kETNrqfcmdDD1kzLQPa75wAP5AzXLO4yXqi7imnIxSFAvIyvD6BOfYArbOOsnK71XdH3uc4udNXJm7GRdg2zeHK4VRikDPOk6+9lnw37Pzx+4fXegbvpFDqpV1UVcUy4GKT5TxhhU9ELNsooXPfRQ/uOtjyJI0TVAnqXem5oBnlWwxix5GfqkQjbUkBofVdVwGnQtKIrP1KEp/S1oj6wW/ktfmq8J2ecctXg3d1JwbMI0t16yxjUPOSTcHnRQ7xERlp4dH1XVcKq4FlTfCOplNKW/Be2SNogcP2ssLISvpUv3PR6dNSoogjRqdZSiXoZrrklffn5+XvqLv5BOPDFfoC45/R8jpKqLuCZcDBLUy2DwDXWan5duuy3c7t0bHoumUMfl6UHqMd4+Sp1Q8a7ytIAuhevsE0/cvzdi1PIFMDhVXcQN+2KQoF5GU/pb0H5R5OoesNu1K3zF+8R79SC5h771jHJpo9IJlTTSkCZ+nU3FOLQViXJV6M4w4pIfVcrK/orEE+V6Zey4S08/nfxc58JgVBaAy8pVjUxO7p/kxupvGEV5E+UoPlOFqL/lssvCLWcEVClrgDuyY0cYUN64MdxPK8V20UX7uu67xSpkVFnNbZCyhgnMpNe+9sCCNU0pEgIMAt3vQBNkjXFnRa6Iu3THHdI99+zrR45n7Bx9dNju859Pb9Z2DZY3Iemnl17DBH/8xweOb45SvgBQFEG9akNfogcjp9cAb1bk6tY97+wd7yiWSdY1WD7spJ9e+slVHZV8AaAfjKlXiYoVKCrPAK+UPsCdJhpjP+us3uPx3b+v4EXosK9ji/7bjUq+ABBH8Zm6NXVRajRb3gnh3bMsliwJ36e9NupHzjMeH88kK3iMNiGLvOgwAZNW0GYE9aqMWrUONEPeAd6kyLV7t/S+92X3I/fKJHvNa8LAcx/N6yZdxxYdJhiFfAGgHwT1qpB9g34UGeDtjlwLC9Kf/Vnyz40GlG+8sXgmWU6Duo6tqzu/6fkCqN+wh5KqQFCvCtk36EeZqoR5+pFXr5b+9E/7+/k9FLmOzXuyTOvOv+ii8HxdJ9o2nNxRTBOGkqpAolxVyL5Bv8omWGYVP5qZCZO1f/GL/V8zORmmwJU4W+VdpDDv28tbY2fQuafku46fUShIlDdRjqBeJc4G6NcgqhIO+EyV5zo2qkibZxfyVIercPcTjcLJHdUbhVW0qSg3DFH2zYYN0rp1B5ayAtIMoirhgEun5ak6V2QX8tTYqXD3E1Ftbjy1KSWqkjF1MztN0t9IWiTpane/ouv5JZI+J+kVkp6Q9Afu/m+d5z4s6d2S9kr6r+5+exX7lMsgBs7IvkFT1HCm6pXnUfbkAAAgAElEQVRFXmQXitTYGdSJtk0nd+TXppSo0i11M1sk6e8kvVHSSyS93cxe0rXZuyX93N1/TdInJf1157UvkXSOpP8o6TRJ/2/n5w3ezEzoZ1uzRlq7NtwuXx4eB9ogOlMlqehMFb8uTho1KLILq1eH0ao8BnWireFPhgbKOvZGbRXt0mPqZvY7kj7i7m/o3P+wJLn7X8W2ub2zzb+Y2YSkn0g6UtIl8W3j22X9ztJj6gycYRwMOHkzTwpJ0V3o/plpp6dB/ZuS7zq+mp4SVWdFuaMkPRq7v1XSb6Vt4+57zOxpSc/tPH5312uPSvolZnaBpAskaeXKleX2mEIxGAcDLJ2Wt/BM0V3o7s53l668MiwsV0flN6rNja+2FCQamXnq7n6VpKuk0FIv9cMYOEObdeeKbNki3XprpWeqItfFRU+W3WkpF19c74m2LSd3FNeGlKgqgvo2ScfE7h/deSxpm62d7vfnKCTM5Xlt9dqUFQHEZfUhVni26nVd3J1/WuZkGb02ulb5+McHXxCmDSd3lDOqBYiqGFOfkPQDSa9VCMgzkt7h7t+NbXOhpJe6+5+Y2TmSVrv728zsP0raKOlkSSskfV3SCe6+N+t3lh5TZ+AMbVRjrkivOeVLloTS9FWNSzZ9vBPt0sTjrbZ56u6+R9J7Jd0u6fuSvuTu3zWzdWZ2RmezayQ918welnSR9iXIfVfSlyR9T9Jtki7sFdArkWeCLTBqapxk3StTfdeuMB4+NxeuMVatCtfS/YiP38/NVfdzgSSjfrxVMqbu7rdIuqXrsctj3y9I+v2U135U0ker2I9CGDhDG8T7CL/97fSmc8W5IkkJZUuWpJ/wyuSfkteKOhU93prWTT8yiXIDwcAZmqKfM0N3H+GSJenbDiBXpPu6+L77pE2bkrctc01BXivqVOR4a+IiMOMd1KXsk2nTLsHQTv2cGZLmlGX1Cw6ogkb8unjjRum226rPPx1WXiv//uMp7/GWd1pn7dx95L5e8YpXeCXuucd9asp9ctLdLNxOTYXHs54DqrJzZziuwtDd/l9TU+7z88mv+8IXwjGZ9DrJfcmS2o/b+fn+3sqwfm4W/v3HV97jLetfcHIyPF8lSbOeIz6Ob0u912WWu/T008nPkR2PqvQ7YJzVR2gmnX669J/+U/5ckQqapYMq3FJ3QZgnn5ROPVXauXPfY/z7t1v34X/zzeFfKOt4a+qw0PgG9ayT6cJCen1KMnNQpX7PDL36CM86K/8xWuHA4KDyT+vKa52ZOTCgx/Hv3z5p09c2b5YefTT9eGtquZPxDepZJ9Osscm5uXBmYYANVej3zLB6dQi8SYqMnw9gYHBQ+ae9fm7ZzoboT5EW0CUS89om6/A/44zsw7+qf8Gqje966lnLMS1dmp1JfPPNrOiGavS7PFRVtRZasoB4FYsuZv0pIhScbJcyh39Ty52Mb0s96zJr6dLQ/b5rV/Lzu3aFLwbYUFaZAeMq+qQLdv83MSO8qs6GrD9FZNSW4US2suPiTSx3Mr5BvdfJVNpXPmgQFTWASJkzQ3ef9Px8mFuWN+oW6P5v4pxcqbriNFl/CklatoyCk21Txbh408qdlK79Pgyla7/HRU2PpJPp/Lz0rnelV9Qwk9atky67rJp9AfJIay73U7A65zoIvcrKRwvBDaMFv3596HJPOpUV+RddWJB+9VelX/ziwOeWLg1/iqmp8vuL5hilZUDqXE99tGVdZh1ySMgiHkRFDbTXIPuo05rLmzeHzJ6ifdA5u/+zWsM7d0rPe5500EFhVGpyst4WfFVZyNEs4yQTE805uaM6dU+XrAMt9V5G6VIOwzfI5Z2ymsuLF4fb3bsPfG5yUtqwIV/qeEr3f1ZrOE307+E+2HH4qv5Fs1aey/MnxOjqcfg3Ai31qrTxUg6DMei6kVnN5aRgHsmT8dNjYLDXeHOSPXvC2uef+ES5cfheHR9V/Ys2tZgIBq9p4+JlENR7mZ+XHnpIet/7pCeekI44QjrxxGZeymG4emVsXX99aFH322TNk56dpIJhotWrQyu2iLk56S//cv880zzXOPEgLklXXint3Zt9UVBFFnJTi4kARRDUs2R1pRLQ0a1XU+897wnHT79N1n6ay1Il87DcQ9JZEVmlHtKy0rv/5bq7+7MuCsq2tppaTAQoYnyLz/QS70qdmwtnl7m5cP/UU5MH8DDesgoauYcma/exFE2bzCOrUE2SxYsrq4Rx003FxtOlcBGQVuohqTs76V8uzSDq4jS1mAhQBEE9Ta9036OOoqIc9lc06ErFolN31MmyeLF0/vmhOVtBCnrRnv/nPEf6b/8t/RonqTs7T0W3yKDGuKNu/A0bwlS4DRsq+xMCtSCop+l1Ftu5s1grC+2X1tTLMjcnPfhg/t8RjzqXXx4qoiRZtiwMRlfUvMzqhIgsXRq+1q2TfvIT6c//vFgF3CIXDoMa446P5zc1CxrIwph6mjzjl1SUQ7fujK3t26Vrrknvh5ak//E/pLe8JX9zMD54/Ja31JL3kTXevGyZdPHFyfmjRbLSi6QMDGKMu6kV84AiCOppss5iEea5IIu79Nhj2QFd2tfr08+Ut5qKT/eaNpYW9IrsXp5/ucnJweSqDno2IlAXgnqa6CyWtbgy81zQrbu5l5UCHlem16emSbb9Xj/k3b34hcMzz4Q/3+LFIYh/4APh+0F1iVdVPx4YNoJ6lpNOkrZulY4+OjmwM88FcUnNvbw5F1m9PvGB3mOOCY89+uhQlkkb9PXDSSeFirennRaCeFRT58orpdtvl175ysH8XgrPoC0I6r0cdph0553J/Y6bN4dL+CatQ4nh6ZW+PTGR/nxar0/WxO26i6zXYH4+lLCPX0Pv3h2+TjlFuuuuwQR2Cs+gLQjqeXT3Ox59tPTDH4aTbbTuegtPsCgoq7lnJp1+eljKLKn1ntTrk9Tyjysz6NvEhdGVfV3kLr3hDaGwY9W7SuEZtAVBPdLrJBf1O0Ytp+4TLVk1yGruLVkiLVoUJm93F0NPy/zKO3G76KBvg9O8H3lkOBNOWOIBbUFQl/Kf5Hq1nKJtrr9eOu+8ge82GiarubewIG3aJN18c7h/+unSS16SvY5A3onbRQZ9G57m/cIX7j+W3m337sGNb9c0kQAYKIrPZJWD7S4uk6fltGtXqPFNtbnxk1R8ptuuXeFr0ybpb/9WOuGE9KiRp+KLVGzQN0+a9xD1Kso36PHtqEPussvCLQEdo4agXuQkl7fltLBAtblxFa/49nu/lz2lrVft97xlZ4sM+jY8zXvp0pDlnrZ4zMEHM74NZCGoFznJ5W05SY1o9WBIoubei1/cu/BM1nHSq+Xfz2ojWcfwwQeH7vchX4y+8pXSP/9zqFS3eHF47NBDWVgFyIMx9SJzWfKUvIo0oNWDIZqZCZOre+l1nJx0Urjw/Iu/kB54QPq1X5Omp6Wf/rS/Qd+sY3j3bulzn5M2bhx60tx//s/Sz37G+DZQlHnR9RQbYHp62mdnZ6v5YQsL0vLlyclvU1MHJg7Fk+qy0nQnJ0MXLGWoxs/8vLRiRb7leXsdJ9/8ZqjEsmdPCLqHHhpa1GWCbnQMR2XbkiQd+wCGxszudffpXtvR/V50EeW8q2QxuXV8FVlDNOs4+eY3pVe9KlRiidLBd+wovg57t+gYPvfcff3b3Rg+AkYS3e9S8bksQ1glCyMkT0Jlr5VJ5udDCz2tJ+2ZZ8pN2D7kEOn5zw8/JwnDR8BIIqhH0opa9ypKw+RWdOtVhOb006Wzzso+Tnq19qsIutRGBVqHMfUsSUVpstaabGjpTdSsaJ5GkvXrw/BOmsWLpc98plzORhX7CaAWecfUaamnKVp5q8GlN1GzKmqOvvCFoYs+LRmz35yN7gvPm28OPQcMHwGtQFBPU2SB5YaX3sQQlB2WyZp6ZhYqtBQ9ptJ6njZvDku5MnwEjDyCepoiRWmKXABgfJRZfDze2o+mni1aFL5uuaX4+qNZF55nnMGFJ9ASTGlLk1V5qzuJqOGlNzGiTjoptKLdwxj63r3h9uyzi68t0PCa7wCqQVBPk1V3u3s8s8gFAJDX/HxoRcfnqactNtQLF56ApPBvtXFjyEXduHHoVZErR1BPU6QoTZELACCvKlvXXHgCmpkJxR7XrJHWrg23y5e3a1FNxtSz5E12qiLbGejWq3V9zTXh+zyJbVmJdw248IyS8h94QHriCemII7KXmgeKGpd85lJB3cwOl3SDpGMl/Zukt7n7z7u2+U1Jn5b0HyTtlfRRd7+h89xnJb1K0tOdzc9z92+V2afK5U12oggNqpZVHMZduuMO6Z578k2dbPCFZ5SUv2tXOPFGli1jViiqMy75zGVb6pdI+rq7X2Fml3Tuf6hrm52S3uXuD5nZCkn3mtnt7h5dL/25u28quR/NUCbbGeiWZ1XAIk2NBl54JrWeIjt3hq82taIwPOOSVlI2qJ8p6dTO99dJulNdQd3dfxD7/sdm9pikIyXlWMIKGGNJreu0CpB5mxpFLzwHXCUxz9o3bWpFYXiOOSZMHtm168Dn2pRWUjaoP8/dt3e+/4mk52VtbGYnS1os6ZHYwx81s8slfV3SJe6e8CeXzOwCSRdI0sqVK0vuNjAi4q3ra64JXe5JBtHUyKqS+Bu/UUmwz7P2TZtaURiOmRnpz/4sOaBLjUgrqUzPoG5mX5P0/ISnLo3fcXc3s9RC8ma2XNLnJZ3r7s92Hv6wwsXAYklXKbTy1yW93t2v6myj6enpegrWU8sdTRBvXd9zTz0LsGRlFb3mNeEsWKIkcvSv9e1vhzVusqYVtakVhfplDfFI0nOeM/S0kkr1DOru/rq058zsp2a23N23d4L2Yynb/QdJX5F0qbvfHfvZUSt/l5l9RtLFhfZ+kKjljqapM4M9q1+8+6KiYApx/F8rrbR9XJtaUahf1qG8dKn0qU+165Redp76Zknndr4/V9KXuzcws8WSbpL0ue6EuM6FgMzMJL1V0ndK7k814pd2c3NhHLPfoh9AVYrUTigrT794txxz57v/tbIsWzaYt4bxknUo79olbd1a7/4MWtkx9SskfcnM3i3pR5LeJklmNi3pT9z9/M5jp0h6rpmd13ldNHXtC2Z2pCST9C1Jf1Jyf6oxLnMfMHrqymDPmk6XJsfgd9a/1pIl0pveFIqBHHmk9Ou/zogXyss6lNs4tFMqqLv7E5Jem/D4rKTzO9//g6R/SHn9a8r8/oEpMveBcXfUreqpk0nHcJ7pdN1ynCGz/rV275Ze/nLpssuK/VogS8PrLlWOinJJ8l7aMe6OUZd1DCcVq1m0SHr2WekXvzjwZ+U4Q45bqwnD1+C6SwNhnjbvtcGmp6d9dnZ2cL9gYSH0ASalS05NhWQg91BEOGubth0taJf5+d7HsPuBXf33359+huxxMZvnX4t/GwxC1CHVkLpLhZnZve4+3Ws7WupJ8lzabdzIuDtGW97cke7jOGlc/41vlG69VbrttsxhqHFrNaE5xqXgJ0E9Ta+EpHGpOYj2KnMMx8+QMzPh/yPnMFQDq9UCrUFQz5J1acfgIEZdFcdwn0tfjUurCagb66n3izXUMeqqOIarXPMdQGkE9X7VWQgEGIQqjmGGoYBGofu9DAYHMerKHsMMQwGNwpQ2AP1jjhpQC6a0NRUV6NAmzFEDGoWgXicq0KGNGIYCGoOgXpc+p/4AI4E5akAjkP1eF6b+AMBQzM+HIqDr14fbNq+eTUu9CnnGyZn6AwC1G7dRT4J6WXmPGKb+AECtxnHUk+73MuJHzNxcWNFqbi7cX7Vq/z4eKtABQK3GcdSToF5GkSOGCnQAUKtxHPWk+72MokcMU38AoDbjOOpJUC+jnyOGqT8AUIvVq0OKU5K2jnrS/V5GlePk4zTnAgBqMI6jnrTUyyhbIjOaCnfXXdLnPy8ddNB4zLkAgJqM26gnC7pUIQrORY6YaCrcM8+kj8uzIAYAQCzoUq+i4+RJkyeTRBn0jMEDAHJgTH0YsqbCxbV1zgUAYCAI6sOQNRUurq1zLgAAA0FQH4ZoKlwvbZ1zAQAYCIL6MGRNhZPaP+cCADAQBPVhSJo8eeih0rJl0po10oYNIeud6WwAgALIfh+WcZs8CQAYOIL6MFEyFgBQIbrfAQBoCYI6AAAtQVAHAKAlCOoAALQEQR0AgJYg+x0AMHaixTUfeSQU+WzLjGKCOgBgrEQrX+/ZE5bhOPRQ6cILQ02wUa/5RVBvkrZeOgJAQyStfD03F25XrQrFPEf5tEtQb4o2XzoCQENkrXy9Z490442jXROMoN4Ebb90BICGyFr5eseOULV7lJH93gR5Lh0BAKVlrXx96KFhGY5RRku9CYpcOjLuDmBM9Tr95Tk9rl4dRjaTTEyE50dZqaBuZodLukHSsZL+TdLb3P3nCdvtlXR/5+6/u/sZncePk3S9pOdKulfSO919d5l9GknRpWPU5R4Xv3Rk3B3AmOp1+st7eoxWvu7edmIiPD7qbSRz9/5fbPYxSU+6+xVmdomkw9z9Qwnbzbn7ZMLjX5J0o7tfb2Z/L+k+d/90r987PT3ts7Ozfe934ywsSMuX7z+mHpmaCmPq7tKKFdnbjPrRCAAJ5uezT3/R6tVJzy9bJn3gA9KLXrR/yz1q1Y/Kytdmdq+7T/faruyY+pmSrut8f52kt+Z9oZmZpNdI2tTP61slunScmpImJyWzcDs1te/SkXF3AGPqhhtC2yfJnj3SZZelnx537pTWr5fWrAltp5mZ8Hi08vVll4XbJgf0IsqOqT/P3bd3vv+JpOelbLfUzGYl7ZF0hbv/o0KX+1PuHn0UWyUdVXJ/RtdJJ0k//nH6pWPbUzYBIMHMjPSe96QH9R07pAceSD89RsZlQlHPoG5mX5P0/ISnLo3fcXc3s7S+/Be4+zYzO17SHWZ2v6Sni+yomV0g6QJJWrlyZZGXjo7o0jFJ3nF3AGiJaLZvWkCXwunvRS+S7rkn+fTYrQ1z0bP07H5399e5+28kfH1Z0k/NbLkkdW4fS/kZ2zq3WyTdKellkp6QNGVm0YXF0ZK2ZezHVe4+7e7TRx55ZIG32BKrV4dMjiRtSNkEgC5Zo46RiQnpox9NPz12a3vHZtkx9c2Szu18f66kL3dvYGaHmdmSzvdHSPpdSd/zkKH3DUlnZ70eHXnG3QGgRbJGHaX9T4vdp8c0be/YLBvUr5D0ejN7SNLrOvdlZtNmdnVnmxdLmjWz+xSC+BXu/r3Ocx+SdJGZPawwxn5Nyf1pt2jcfcMGad26cLt9O9PZALRSVqGYJUukT3963+kvOj1+8IPS4sXpP7PtHZulprQNS+umtA0CRWoAjIi001We2b7dxWfSpr5FrxnVsh55p7RRUa6NKFIDYET0Ol0VKRSTNQa/ZIn0yU+2/xRIUG8bFocBMCLynK56zfaNyxqD371b2rq1+vfQNAT1tmn7uoIAWiPv6So+23d+PjyeNLLIzF+Cevv0KlLDODuAhihaU6tXV33bF2vJg6VX2yYrXdRduvlmae3aA2smAkDNiiyDGu+qn5sLp7O5uXA/KlDDzF+CevtkFamRpF27kv8bAKBmRWpq5V3+Ytxn/tL93jZJ6wouWZK9GgLj7ACGoMgyqEW66rMqbrcdQb2NutNF77tP2rQpedu210wE0Gh5s9tJgsuHoN5W8UvVjRul227jvwFAI+VpWZMElw9j6uOAxWAAjDiS4PKhpT4OigxcAUBDFSlEM64I6uOC/wYALTDOSXB5ENTHCf8NANBqjKkDANAStNTHCcuxAkCrEdTHBcuxAmgA2haDZe4+7H0obHp62mdnZ4e9G6Njfl5asWL/9Q0jU1MsxwqgFklti2gSDm2LbGZ2r7tP99qOMfVxkFU0eWFBeuc7Q4EaasADGJA8C7KgPIL6OMgqmrywEErIsmobgAHKuyALyiGoj4Os9Q0jXDIDGKCia6ejPwT1cdBrOdY4LpkBDECRtdPRP4L6OOgumpyFS2YAA8ASFPUgqI+LqEzshg3S2WeHNdaTcMkMYABYkKUeTGkbRwsLISkuaYrbsmXStm3hPw0AKhbNU2cJimLyTmmj+Mw4iq/atmtX+C+LO+64AyeOUjECQAVYgmKwaKmPsyeflI45Rtq588Dn4kVpqBgBAENF8Rn0dttt0kEph0CUBU/FCAAYGQT1cZZn4igVIwBgZBDUx1meiaNUjACAkUFQH2d5Jo5SMQIARgZBfZzdf39y1/qv/Mq+iaNUjACAkUFQH1dRAtzc3IHPHXSQ9NKXhu+pGAEAI4N56uMqKwFu796QABdNJo2q0VExAkAfKHNRH4L6uCqaAEfFCAB9SCpzceGFlLkYFIL6uIoS4JK630mAA6DyLex4mYtIdMpZtWpffStUhzH1cUUCHIAMMzPSihXSmjXS2rXhdvny8HhelLmoH0F9XJEAByBFVYUkKXNRP7rfxxkJcAAS5Glh50mxYZSvfgT1cUcCHIAuVbWwV68OSXFJGOUbDLrfAQD7qaqQJKN89WPpVQDAfhYWQlJcPGs9El+VOa8oi55Rvv7lXXqV7ncAwH6iFnb3/PKJif5a2Izy1adUUDezwyXdIOlYSf8m6W3u/vOubV4t6ZOxh14k6Rx3/0cz+6ykV0l6uvPcee7+rTL7hJpRKgpoJfJoR1Op7ncz+5ikJ939CjO7RNJh7v6hjO0Pl/SwpKPdfWcnqP+Tu28q8nvpfm+IpFJR0aU8paIAoDJ5u9/LJsqdKem6zvfXSXprj+3PlnSru+8s+XsxbFVNZAUAVKZsUH+eu2/vfP8TSc/rsf05kr7Y9dhHzezbZvZJM1tScn9QF0pFAUDj9BxTN7OvSXp+wlOXxu+4u5tZal++mS2X9FJJt8ce/rDCxcBiSVdJ+pCkdSmvv0DSBZK0cuXKXruNQatiIivj8QBQqZ5B3d1fl/acmf3UzJa7+/ZO0H4s40e9TdJN7v5M7GdHrfxdZvYZSRdn7MdVCoFf09PTozcPr23Klopi6SYAqFzZ7vfNks7tfH+upC9nbPt2dXW9dy4EZGamMB7/nZL7g7qUWRCG8XgAGIiyQf0KSa83s4ckva5zX2Y2bWZXRxuZ2bGSjpH0z12v/4KZ3S/pfklHSPrvJfcHdSlTKorxeKBy8/PSxo3S+vXhlmvj8VRqnrq7PyHptQmPz0o6P3b/3yQdlbDda8r8fgxZvxNZWboJqFTe0SzSWNqPinIop59SUSzdBFQmPpoVif61Vq3aV9KVNJbxwIIuqF+e8Xj6EoFc8oxmkcYyPgjqqF+v8fj775dWrJDWrJHWrg23y5eHpgaA/eQZzSKNZXzQ/Y7hSBuPdw8BvVdfYhwDhRhjeUazSGMZHwR1DE/SePzGjb2bFPHXlBko5GIALbB6dTjkk0SjWTfeSBrLuKD7Hc1SpElRZqBwZoYufrRCntmlZcpKYLTQUkezFMmMzzNQmJSZnzddGBgRvWaXVr0+OpqLoI5mydOXGOl3oLDfiwGgwXrNLmV99PFAUEezFGlS9DvfnawhjKl+ykpgtBDU0Tx5mxRFWvVxFL8B0FLmPnoLnk1PT/vs7OywdwNNkJT9HrXq07LfFxZCUlx8TD0yNcWYOoDGMbN73X2613a01DHa+hkoJGsIQEsR1DG6uueZX3xx/oBM1hCAFiKoYzRVsToFWUMAWoagjtHDPHOMIAoYog4EdYwe5pljxLDsKepCUMfoYZ45RggdS6gTtd8xeqJ55kmYZ46GYdlT1ImgjtHD6hQYIXQsoU4EdYyePMtSAQ2R1bHkHr6AqlBRDqMrSidmnjkaLKuAoSQ95znST37CoYtsVJRD+zHPvH1aOO9r6VLpooukyy9Pfn7vXiZsoDoEdQDNMKbzvhhXR5UYUwcwfPF5X3NzYaB5bi7cX7Uq9GGPsBe+MKR9JGHCBqpESx2I9Or6bWHXcC3y/N1aXlCo31WCgaII6oDUu+t3TLuGS8v7d2v5vC8WBkRdCOpAr5JfW7ZQEqwfRUqpRfO+oufjWtI/zcKAqANBHejV9XvZZcPrGh7lLv8iXepj0j/NhA0MGkEd6NX1+8ADw+kaHvUu/wceSG55Swf+3Ua0f3qUr7nQTgR1oFfX74teJN1zT71dw6O+CsjMjHTllenPJ/3dRqx/etSvudBOVJQDskp+TU1JP/yhdNxx6c8PIsBu3CitWZN8ITE5KW3Y0Nx+3Pl5acWK9BJq0uD+bjXJeosj/tbQUHkryjFPHehVSz5+W1et+SZkg8/Ph4uL9evDbd654llj6ZK0bFmxv1u/+zFArLyGpqL7HZB6d/3W3TU87GzwMn3LWRckknTxxfn7pxvax92Eay4gCUEdiPRKTa4zdXmY2eBlx/OzLkgmJ6UTT6xnPwZo2NdcQBq634Em6h4SkKQlS/atDjJIZfuWq1rvvsF93FW9RaBqBHWgqaIu/w9+MARzM2nXLuljHwuJfTMzg/m9ZfuWq1rvvsF93FW9RaBqdL8DTfeJT+yfHDboLugq+paryEFoeB/3iM3Aw5hgShvQZMOY2tZril9dY9lN2Q+gAZjSBrTBMLqgm9K33JT9AEYI3e9Akw2rC7opfctN2Q9gRND9DjQZXdAARPc70A50QQMogO53oOnoggaQU6mgbma/L+kjkl4s6WR3T+wTN7PTJP2NpEWSrnb3KzqPHyfpeknPlXSvpHe6++4y+wS0EgtxA8ihbPf7dyStlnRX2gZmtkjS30l6o6SXSHq7mb2k8/RfS/qku/+apJ9LenfJ/QEAYGyVCuru/n13f7DHZidLetjdt3Ra4ddLOtPMTNJrJG3qbHedpLeW2R8AAMZZHYlyR0l6NHZ/a+ex50p6yt33dD2eyMwuMLNZM5t9/PHHB7azAACMqp5j6mb2NUnPT3jqUnf/cvW7lMzdr5J0lRSmtNX1ewEAGBU9g7q7v67k79gm6ZjY/aM7jz0hacrMJv0REdgAAAcZSURBVDqt9ehxAADQhzq632cknWBmx5nZYknnSNrsoerNNySd3dnuXEm1tfwBAGibUkHdzM4ys62SfkfSV8zs9s7jK8zsFknqtMLfK+l2Sd+X9CV3/27nR3xI0kVm9rDCGPs1ZfYHAIBxRplYAAAajjKxAACMGYI6AAAtQVAHAKAlCOoAALQEQR0AgJYgqAMA0BIjOaXNzB6X9KMKf+QRkn5W4c8bJt5LM/Femon30lxtej9VvJcXuPuRvTYayaBeNTObzTP/bxTwXpqJ99JMvJfmatP7qfO90P0OAEBLENQBAGgJgnpw1bB3oEK8l2bivTQT76W52vR+ansvjKkDANAStNQBAGgJgjoAAC0xNkHdzH7fzL5rZs+aWerUAjM7zcweNLOHzeyS2OPHmdm/dh6/wcwW17Pnift4uJl91cwe6twelrDNq83sW7GvBTN7a+e5z5rZD2PP/Wb97+KX+9nzvXS22xvb382xx0ftc/lNM/uXzrH4bTP7g9hzQ/9c0o7/2PNLOn/nhzt/92Njz3248/iDZvaGOvc7SY73cpGZfa/zOXzdzF4Qey7xeBuWHO/lPDN7PLbP58eeO7dzTD5kZufWu+cHyvFePhl7Hz8ws6dizzXtc7nWzB4zs++kPG9m9red9/ptM3t57LnBfC7uPhZfkl4s6URJd0qaTtlmkaRHJB0vabGk+yS9pPPclySd0/n+7yW9Z4jv5WOSLul8f4mkv+6x/eGSnpS0rHP/s5LOHvZnUuS9SJpLeXykPhdJvy7phM73KyRtlzTVhM8l6/iPbfOnkv6+8/05km7ofP+SzvZLJB3X+TmLGv5eXh37n3hP9F6yjrcGv5fzJP0/Ca89XNKWzu1hne8Pa/J76dr+fZKubeLn0tmfUyS9XNJ3Up5/k6RbJZmk35b0r4P+XMampe7u33f3B3tsdrKkh919i7vvlnS9pDPNzCS9RtKmznbXSXrr4Pa2pzM7+5B3X86WdKu77xzoXvWn6Hv5pVH8XNz9B+7+UOf7H0t6TFLPKlE1STz+u7aJv8dNkl7b+RzOlHS9u+9y9x9Kerjz84al53tx92/E/ifulnR0zfuYV57PJc0bJH3V3Z90959L+qqk0wa0n3kUfS9vl/TFWvasD+5+l0KDKc2Zkj7nwd2SpsxsuQb4uYxNUM/pKEmPxu5v7Tz2XElPufuerseH5Xnuvr3z/U8kPa/H9ufowH+Mj3a6gz5pZksq38P88r6XpWY2a2Z3R8MIGvHPxcxOVmitPBJ7eJifS9rxn7hN5+/+tMLnkOe1dSq6P+9WaFFFko63Ycn7Xn6vc+xsMrNjCr62Lrn3pzMccpykO2IPN+lzySPt/Q7sc5mo4oc0hZl9TdLzE5661N2/XPf+lJH1XuJ33N3NLHVeYueq8KWSbo89/GGFoLNYYf7khyStK7vPGftQxXt5gbtvM7PjJd1hZvcrBJRaVfy5fF7Sue7+bOfhWj8XBGb2h5KmJb0q9vABx5u7P5L8ExrhZklfdPddZrZGoTflNUPep7LOkbTJ3ffGHhu1z6V2rQrq7v66kj9im6RjYveP7jz2hEK3yUSndRI9PjBZ78XMfmpmy919eyc4PJbxo94m6SZ3fyb2s6PW5C4z+4ykiyvZ6RRVvBd339a53WJmd0p6maT/TyP4uZjZf5D0FYWLzbtjP7vWzyVB2vGftM1WM5uQ9ByF/488r61Trv0xs9cpXJC9yt13RY+nHG/DCh4934u7PxG7e7VCfkf02lO7Xntn5XuYX5Hj5BxJF8YfaNjnkkfa+x3Y50L3+/5mJJ1gIaN6scJBtdlDZsM3FMamJelcScNs+W/u7EOefTlgTKoTcKIx6bdKSszcrEnP92Jmh0Vd0WZ2hKTflfS9UfxcOsfVTQrjbJu6nhv255J4/HdtE3+PZ0u6o/M5bJZ0joXs+OMknSDpnpr2O0nP92JmL5O0QdIZ7v5Y7PHE4622PT9QnveyPHb3DEnf73x/u6RVnfd0mKRV2r/Xrm55jjGZ2YsUEsj+JfZY0z6XPDZLelcnC/63JT3duXgf3OdSVRZg078knaUwbrFL0k8l3d55fIWkW2LbvUnSDxSu/i6NPX68wknqYUn/S9KSIb6X50r6uqSHJH1N0uGdx6clXR3b7liFK8KDul5/h6T7FYLGP0iabPJ7kfR/d/b3vs7tu0f1c5H0h5KekfSt2NdvNuVzSTr+FYYAzuh8v7Tzd36483c/PvbaSzuve1DSG4f1ORR4L1/rnAuiz2Fzr+Otwe/lryR9t7PP35D0othr/7jzeT0s6Y+a/l469z8i6Yqu1zXxc/miwgyWZxTiy7sl/YmkP+k8b5L+rvNe71ds5tWgPhfKxAIA0BJ0vwMA0BIEdQAAWoKgDgBASxDUAQBoCYI6AAAtQVAHAKAlCOoAALTE/w+4PKVPLNGShAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.io\n",
    "import matplotlib.pyplot as plt\n",
    "mat = scipy.io.loadmat('hw5_p1b.mat')\n",
    "print (mat.keys())\n",
    "X = mat['X']\n",
    "k = 2\n",
    "num_points = int(X.size/2)\n",
    "cluster_assignment = np.zeros(num_points, np.int8)\n",
    "# Guess the mus\n",
    "changed_points = 1\n",
    "num_dims = 2\n",
    "\n",
    "mus_choice = np.random.choice(num_points, k)\n",
    "\n",
    "mus = []\n",
    "for i in range(k):\n",
    "    this_mu = X[mus_choice[i]]\n",
    "    x_val = this_mu[0]\n",
    "    y_val = this_mu[1]\n",
    "    mus.append(this_mu)\n",
    "\n",
    "def get_distance(point, mu):\n",
    "    # We can skip the sqrt\n",
    "    return np.dot(point-mu, point-mu)\n",
    "\n",
    "\n",
    "print(get_distance(np.asarray([0, 1]), np.asarray([0,3])))\n",
    "\n",
    "steps = 1\n",
    "while(changed_points != 0):\n",
    "    changed_points = 0\n",
    "    \n",
    "    # ASSIGN POINTS TO CLUSTERS\n",
    "    for point_counter, point in enumerate(X):\n",
    "        distance_vector = np.zeros(k)\n",
    "        for counter, mu in enumerate(mus):\n",
    "            distance_vector[counter] = get_distance(point, mu)\n",
    "        new_assignment = np.argmin(distance_vector)\n",
    "        old_assignment = cluster_assignment[point_counter]\n",
    "        if new_assignment != old_assignment:\n",
    "            changed_points = changed_points + 1\n",
    "            cluster_assignment[point_counter] = new_assignment\n",
    "    \n",
    "    # UPDATE MUS\n",
    "    for counter, mu in enumerate(mus):\n",
    "        my_points = X[cluster_assignment==counter]\n",
    "        mus[counter] = np.sum(my_points, axis=0)/num_points\n",
    "    if steps == 1:\n",
    "        cluster_assignment_legacy = np.copy(cluster_assignment)\n",
    "    steps = steps + 1\n",
    "        \n",
    "red_X = X[cluster_assignment==0]\n",
    "blue_X = X[cluster_assignment==1]\n",
    "\n",
    "plt.figure(figsize = (8, 8))\n",
    "plt.plot(red_X[:,0], red_X[:,1], 'ro', ms = 7)\n",
    "plt.plot(blue_X[:,0], blue_X[:,1], 'bo', ms = 7)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10d735358>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAHVCAYAAADo5GInAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+UHHWd7//XmwxJGGbXAWE1ASLgIurVe1cd2N3vehF/RfwBSJZ10bMKu3LJuug9LrKKF5Z4kusuqxd19379uuEAij8ieHNhDSs/roosnnMXmeEcEX+AQPxBYhQE4ThJT0LC+/vHp8tUOlXVVd3V1dXVz8c5c3qmu3qmerq63vX5fN6f98fcXQAAYPQdMOwdAAAA5SCoAwDQEAR1AAAagqAOAEBDENQBAGgIgjoAAA1BUAcAoCEI6gAANARBHQCAhpgY9g704rDDDvOjjz562LsBAEAl7r777l+6++HdthvJoH700Udrbm5u2LsBAEAlzOwnebaj+x0AgIYgqAMA0BAEdQAAGoKgDgBAQxDUAQBoCII6AAANQVAHAKAhCOoAADQEQR0AgIYoJaib2dVm9oiZfTflcTOzfzKzB83sO2b20thjZ5vZA+2vs8vYHwAAxlFZLfXPSDol4/HXSzqu/XWepE9JkpkdKmmNpN+XdKKkNWZ2SEn7BADAWCklqLv7HZIez9jkdEmf9eBOSdNmtkzS6yR91d0fd/dfSfqqsi8OAABAiqrG1I+Q9HDs5y3t+9Lu34+ZnWdmc2Y29+ijjw5sRwEAGFUjkyjn7le4+4y7zxx+eNfV5wAAGDtVLb26VdJRsZ+PbN+3VdLJHfffXtE+AcPTakk33CA99JD03OdKq1ZJS5cOe68AjLiqgvomSe82s2sVkuKedPdtZnarpL+LJcetlPTBivYJGI7ZWWnlSmn3bmn7dungg6Xzz5f+z/+RTjhh2HsHYISVEtTN7IsKLe7DzGyLQkb7gZLk7v8s6SZJb5D0oKQdkv68/djjZrZO0mz7V61196yEO2C0tVohoD/xxN775ufD7cqV0rZttNjLRq8IxkgpQd3d39rlcZd0fspjV0u6uoz9AGrvhhtCCz3J7t3S9ddLb3tbtfvUZPSKYMxU1f0OQAqtxe3bkx/bvl3avLna/WkyekUwhkYm+x1ohOc+N7QWkxx8sHTssdXuT5Pl6RUBGoagDlRp1SppIqWDbGIiPI5yZPWKzM+HoL+wUO0+AQNGUAeqtHRpGM+dnpampiSzcDs9He6nO7g8Wb0ikvSv/yotWxbG3YGGsJDDNlpmZmZ8bm5u2LsB9C7KyN68OXS5j1tGdi8Z6UWfs7AQgnZ8TD3J9DTj66g9M7vb3We6bkdQB1CppIz0iYnsjPRenhN/Xqsl7dyZvM3UlLR+PbMOUGt5gzrd7wAGr9WSNmyQLr1UOvnk0Hqen5fcw+0TT4TgmzTGHc9iz/ucyAknSD/7mXTqqenbMOsADUJQBzBYs7PS8uXS6tXSunXSjh3J26VlpPebxX7QQdIZZ4QWeRJmHaBBCOoABqezlZ0lrcVcxtx+Zh1gTBDUAQxOViu7U1qLuYy5/cw6wJigohyAwclqZXdKazGvWhVKuxZ5TpJofH2cZx2g8QjqQBYWA+lP1MrO6nqfmtqbyZ70v41a2WnZ70Xej4MOqkeWO8cVBoQpbUCaXqdRYa+sueKTk9KFF0rHH19snvqot7KLHFcEf7QxTx3oR6sVMraTghHFSoppysVRGQG2yHHVlP8bSpE3qNP9DiTJSvDasUN63/ukyy8nsOcx6LHsKlqz/S7hGu3jDTeE75PEl95lhTn0iKAOJMlK8Nq1S7ryylBMhVZTPoMay65ivfS8ATbt4iK+j1m5BfHpeXnm5tchNwC1Q1AHknRL8Nq1K3zRahqeqlqzeXpt3vrWULWu8+Ji0ybptNO615+X9p2eV8bcfIwl5qkDSbKKlcSxLvfwVLVeep5em5NOSi5je8op+efpx6fnlTE3H2OJoA4kiRcrWbw4fTtaTcNTVWu22xKuu3aFQJ6kW5e7lFwEhwp46BFBHUgTJXide256YKfVNDxVtWbz9tok2bUrOzifeKL0P/9nGCqI5wBQAQ89Ykob0E3WXGumtw1Ple9LlOy2Y0cI1GU5+GDpwAPTE/uaMjcffWNKG1CWMiuaoTxVvi9Rr82FF4Yx9LICezR8kJbYV5cKeBgZtNSBLPFpSkcdFe7bsoVWU51U2ZrN6h0wC1XyduyQlizJXue909SUtH49ARypaKkD/aKi12iosjWb1TuwaZP08MPh4uKee6SNG/P/XhIuURKCOpCEil5Ik6dC3oYN0i23dM98j5BwiZIQ1IEkVPRClm69A1nLxSZhmhpKQlAHklDRC/1I66aXwti7+/ATLlkBrpEI6kCSrDKxdJUij7RuevfhT1OromY+hoLsdyAJc9PRVCwrPJLyZr9TUQ5IQkUvDEKrFZLo1q0Lt0WmvZWlqpr5GAq63zF8dR3bG/Q64BgvdenyJl+k0QjqGK66nOjSUNELZajTFEnyRRqN7ncMT/xE17lk5cqVw+maBAahTl3erADXaAR1DE+dTnTAIOXt8q5izJ18kUaj+x3Dw9gexkWeLu8qh6LIF2ksgjryKzuhjbE9jIusCnMTE9LrXx+O9yrH3MkXaSS635HP7GyY27p6tbRmTbhdtizc3yvG9jAuunV533wzQ1EoBS11dDeozF3WKcc4yeryvuWW0RiKquv0U/wGQR3dDXJxE8b2ME7SurwHNRRVZhCu+/RTSCKoI49BJ7Qxtodx123MvZehqDKDcJ3m2SMTY+roLmpFJCGhDehf2dPMyq4BwfTTkUFQR3cktAGDFw1FrV8vrV0bbrdt661ru+wgzPTTkVFK97uZnSLpHyUtknSlu1/W8fjHJb2y/eOkpN9x9+n2Y3sk3dt+7KfufloZ+4QSkdAGVKOsoaiygzDTT0dG30HdzBZJ+qSk10raImnWzDa5+/ejbdz9r2Pbv0fSS2K/ouXuv9fvfmDAyk5oI4sWGJyyg/AgxvwxEH2vp25mfyjpQ+7+uvbPH5Qkd//7lO3/r6Q17v7V9s/z7j5V5G+ynvqIS0rgiVr9ZNEC/VtYCHUkylwzPelzu2iR9L73hce5OB+ovOuplxHUz5R0iruf2/757ZJ+393fnbDtcyTdKelId9/Tvm+3pG9L2i3pMnf/l5S/c56k8yRpxYoVL/vJT37S135jSFqtUMQm6WQzOSlt3RpOOgD6M4iL56iHbfPmkHz3sY9xcV6RvEG96iltZ0naGAX0tue4+1YzO1bSbWZ2r7s/1PlEd79C0hVSaKlXs7soXVYCz44d0hFHSLffzkkB6NcgakBEY/5JF+dMcauFMoL6VklHxX4+sn1fkrMk7TMw4+5b27ebzex2hfH2/YI6GiIrgUcKgZ2TAlCOQdWAGGRBKvSljClts5KOM7NjzGyxQuDe1LmRmT1f0iGS/j123yFmtqT9/WGS/kjS9zufiwbJmvMeYd4rUG9McautvoO6u++W9G5Jt0r6gaQvufv3zGytmcWnp50l6VrfdxD/BZLmzOweSd9QGFMnqDdZ1pz3CCcFoD6S1ninIFVtlTKm7u43Sbqp475LO37+UMLz/q+kF5exDxgR0Zz3k08OXe1JOCkA9ZBWavbGGylIVVNUlEP1TjhB2rIlZLsn4aQADF9WqdlTTw2BvayytigNQR3DccghIcudkwJQT92S4X760zC2/o53SK98Zbj90Y+YuTJkrNKG4WHZVaC+uiXD3XFH6IqPuubvuiuMuTNPfagI6iiuzBKvLLsK1FNWqdnJSelzn9s3L4Z56rVA9zuKmZ0NRSdWr5bWrAm3y5aF+wE0R9ZMFXfpgJTwwZTUoSKoI7+y12gGUF9Za7y//e3MU68pgjryu+EG6amnkh/j6hxonrQ13k86iXnqNcWYehMNalnTO+7g6hwYN0l5LyzFWlsE9aZJKxbRb0ZqqxUSY9JMTnJ1DoyLqGs+bRU4kuSGhqDeJPEx70hZGak33JCeGCOF8XWuzoHxwZTUWiKoN8kgV07qtrra29/OhxkYN/1MSR3UMOGYI6g3ySBXTsqaszo1FRJnACCPQQ0Tguz3RhnkyklZc1ZJjAGQF1NjB4qg3iSDDLxZc1ZJjAGQV55hQvSM7vcmGXRGKokxAPo1yGFCENQbZ9CBl1rtAPqRlZ9D4Zq+mbsPex8Km5mZ8bm5uWHvBgCgqIWFsF5EfOptZHqaxWBSmNnd7j7TbTvG1AEA1SE/Z6DofgcAVIv8nIEhqAMAqkd+zkDQ/Q4AQEMQ1AEAaAiCOgAADUFQBwCgIQjqAAA0BNnvAIDhYynWUhDUAQDDxVKspSGoAwCGJ74UaySqC79y5b5lY2nNd0VQBwAMT56lWN/2NlrzORHUAQDDk2cp1iKt+TFH9ntTtVrShg3SunXhdmFh2HsEAPuLlmJNEi3Fmqc1D0kE9WaanZWWL5dWr5bWrAm3y5aF+wGgTlatkiZSOo0nJsLjeVrzkERQb554N9X8vOQebp94ItxPix1AneRZijVPax6SCOrNQzcVgFETLcW6fr20dm243bZtbwJcntY8JJEo1zx0UwEYRVlLsUat+c7s94mJva15SCKoN0/UTRVlhsbRTQVgVEWt+RtuCI2TY49lnnoCc/dh70NhMzMzPjc3N+zdqKeFhZAUF5/6EZmeZuoHAIwgM7vb3We6bceYetPkSToBADQS3e9NRDcVAIwlgnpTZSWdAAAaie53AAAagpY6AGC0sXrbb5TSUjezU8zsfjN70MwuSnj8HDN71My+3f46N/bY2Wb2QPvr7DL2BwAwJiiLvY++W+pmtkjSJyW9VtIWSbNmtsndv9+x6XXu/u6O5x4qaY2kGUku6e72c3/V734BABqO1dv2U0ZL/URJD7r7ZnffJelaSafnfO7rJH3V3R9vB/KvSjqlhH0CADQdZbH3U0ZQP0LSw7Gft7Tv6/THZvYdM9toZkcVfK7M7DwzmzOzuUcffbSE3QYAjDTKYu+nquz3GyUd7e7/UaE1fk3RX+DuV7j7jLvPHH744aXvYO2xPjoA7Kvb6m3bto3dObOM7Petko6K/Xxk+77fcPfHYj9eKekjseee3PHc20vYp2aZnd1/IYPzzw8V4qJVjABg3KxaFc6FSbZvl665RtqxY6zOmWW01GclHWdmx5jZYklnSdoU38DMlsV+PE3SD9rf3ypppZkdYmaHSFrZvg8R1kcHgGRpZbHNwrly+/axO2f23VJ3991m9m6FYLxI0tXu/j0zWytpzt03SfqvZnaapN2SHpd0Tvu5j5vZOoULA0la6+6P97tPIyPP3Mo8iSBUjgMwrjrLYm/bFlroSWPtY3DOLKX4jLvfJOmmjvsujX3/QUkfTHnu1ZKuLmM/RkreLnUSQQAgW7ws9rp1ocs9yRicMykTOwxFutSzEkEOPDBclTa8OwkAcuuWPHfssdXuT8UI6sNQZG7lqlXSREqHyq5d0mc/O9bVkwBgH1nnzImJ8HiDEdSHoUiXejwRJOnqc4wSQACgq7TkuenpcH/DK8wR1IehaPdQlAhy9tnS4sXJzxvT6kkAsJ/onLl+vbR2bbjdtq3x09kkVmkbjqy5lWndQwcdJD372dJTTyU/bwwSQAAgt3jy3BihpT4MvXYPjXkCCAAgGy31YemcW3nssfvPU++cx/6GN4x1AggAIBtBfZiyuoeS5rFPTEif+IT03vfuf/8YJIAAALIR1Osoa43g9743tOxvvjm9hQ8AGEsE9TrqNo/95pvHMgEEAJCNoF6GPDXcizyH0rAAkF8v5+CGIqj3q5dlUbs9J8pyj7rc48hyB4C9WJp6H+buw96HwmZmZnxubm7YuxGuDpcv33fsOzI9HYoddF4t5nmOFEq/Fvm9ADBuejkHjygzu9vdZ7ptxzz1fhSp4V7kOWNe5hAAcunlHNxwdL/3o5ex77zPyTOPHQDGWa/5Rw0egyeo96OXse8izxnTMocAkEue82lnAF+xQjr11MaOwTOm3o+FheJj3708BwCwv27n0xtv3DeAT05KO3ZISXGv5udfxtSr0MvYN+PlAFCOrPPppk0hoD/xRGjJu4fAntaQbcgYPN3v/epl7JvxcgAoR9r59Prr05PokjSkBghBvQy9jH0zXg4A5Ug6n2Yl0SVpSA0Qut8BAM2TtVR1koasdElQBwA0z6pV6UtVR2PvDcxpovsdANA8URJd0hLWmzZJDz/cyJwmgjoAoJnGMCmZoF4HDa5uBABDNWZJyQT1YWOFIQBASQjqw9RqhYAer4YUlTtcubLW1Y0AAPVD9vswscIQAKBEBPVh6nWFIQAAEtD9Pky9rPIGDBA5m8BoY5W2YWLFNtRIUs7mxAQ5m0AdsErbKGDFNtREPGczWtBqfj78vHJluP4EUH90vw/bGBZHQP3kydkco6m+wMgiqNfBmBVHwPB1jp3fdx85m0ATENSrQgYSaiJp7Pzpp8Ph2Grtvz05m8DoIKhXgapxqImsekdmyc/pdUVKrmOB6hHUB42qcaiRrLHzgw4Kt2bh2nPx4hDQb7yx+CHKdSwwHGS/DxpV4zAArZa0YYO0bl24zZudnlXvqNWS3vrWkPm+eLG0a5d0wAHSqaeGIF1k39Iy6U86SfrMZ8imBwaFoD5oDz2UXFxGCveTgYSCZmel5cul1aulNWvC7bJl+QJvVO8oycEHS1/8orRjRwjoUm/T2rKuYxcWpHe9K//+AiiGoD5oRx2V/fiRR1azH2iEfueTv/71ISkuydNPp4+rF+lUyuoNkMI+Mv8dGAyCOjBC+hnNmZ0NWexpRSR37w6t9CRFprVl9QZ0/j1Gn4BykShXpqR034cfDs2ftDPpnXdK55xT6W5idPW6BlBSvmanqMs9SZ5pbdHhf//96b0BefcXQG9KCepmdoqkf5S0SNKV7n5Zx+MXSDpX0m5Jj0r6C3f/SfuxPZLubW/6U3c/rYx9qlxauu8FF0iTk+ln4s99TvrEJ8iARy69rgGU1cLPo9u0ts7D/6CDwrXsxIT01FPJz+l2ocCUOKC4vhd0MbNFkn4o6bWStkialfRWd/9+bJtXSvqWu+8ws3dJOtnd/7T92Ly7TxX5m7Vb0KXVCplLSc2gZzwjnNXS+jWnpqT166koh1x6WQOo1ZLe8Q5p48Z8f2Pp0hCM8y7qknX4T05Ke/ZIO3fm31+JxWWATlUu6HKipAfdfbO775J0raTT4xu4+zfcPYpqd0pqVnZYVjNozx7pD/4g/bn0QaKAomsARZnyX/5yvt8/NSV96lPhOnPt2nC7ebP0wAPp0+eyDv8DDpAuvrjYmkUsLgP0rozu9yMkPRz7eYuk38/Y/p2Sbo79vNTM5hS65i9z939JepKZnSfpPElasWJFXztcum4DnYceGs5krJuOEuRdA6jVkl79aunXv87/uycmpLPO2vu7ouS6rCIy3Q5/s2JrFrG4DEZKzcaJKk2UM7M/kzQj6RWxu5/j7lvN7FhJt5nZve7+UOdz3f0KSVdIofu9kh3Oq9tA5xvfKH3ta8nP7bUGJ8ZanjWArrsuX0A327d7Ozof5S2GmGecv8iaRb0mAwKVq2HpxDK637dKik/GPrJ93z7M7DWSLpZ0mrv/ZoTN3be2bzdLul3SS0rYp2qtWhXOiEmipg/rpqNiX/lK9uMzM9KnP723m33btn3PQ3mnz3U7/Ites3YrkEPHFmqhpuNEZQT1WUnHmdkxZrZY0lmSNsU3MLOXSFqvENAfid1/iJktaX9/mKQ/kvR9jZo8A51Rn2l8sLLzLApU6IADwozLtO7w++5LL4YYbzEXHefvpuyLBGAgaloCvO/ud3ffbWbvlnSrwpS2q939e2a2VtKcu2+S9FFJU5L+l4WSVdHUtRdIWm9mTytcYFwWz5ofKXkGOlk3HV2UOTz3xjdmZ7x/+9uh9zCpx3B2Vrr88vTndraY847z5xFdJKRlv9OxhVqo6ThR31PahqGyKW15zrA1S5LAaEg6bO69t9xpXAsL0uGHp7e2O01OSlu3SkuWpE9Ri2RNRytL9D/q9yIBGIgNG8LCC0kfsAFMVc47pY2gnibPRFkm06IHaYfN7t3J54d+AujsrPTa14YAv3OndOCB6cVgpBDYL7pI+shH0i8GJiel22/nEMeY66VoRB8I6v3IqqYRvVnu3behWYEOjz8e1vhJq0WUpN+L/niL9557uhehWbo0XAAknRrMpEsvlT70od72BWiUCht2eYM6td+T5E2AYDItCpidlU4+uVhAl/ofnouncmzYIN1yS3aXvHvogk9K3j34YOl5z+t9X4BGKTOZpCQE9SR5EiDca5kkgXqKZr8UDehSudO4Vq0KSXFZdu6UFi9Ofozsc6BDzRKgWXo1SZ6JskymRQH9LKhSZiCNMssnJ7O3c99/bXUz1h4C6o6gniTPRFkm06KArM6fNEuXDqY+0QknSFu2pLfGpZBM1zmm7i69973UXgfqjKCeJE81DfewrOrSpeGLKnHIkNWxk+ZNbwr5li96URgL71xQpdVKvj+PQw6R/vzPi+2PNNSaGgByYEw9TVYCRDzjcWEhZBUtWSK9//3S3/wNAR37yTOWHTc1JZ1xRvLc9fPPD93g731vaFFv3x5a3f/lv0i33iq9/OX5/sZJJ0mf/3yxHgTSRYB6Y0pbUXmmuxHUkSBp9sv27clTx6an915LJh1qZulTzu64I19gz5pmm2YANTUA5MCUtkFhXUj0KKnzZ8UK6dRTk6e53nxz+qGWdi3uLr3uddJjj3W/tkwrx7pokfT008krvJEuAtQbQb2omtb7xWhImv2SNspzyy3Fk+uk0CWf99oybZQpq2QtHVFAfRHUi8qzeDRQQNo016xDLctTT4Uu+LwdRkl/v4Y1NQDkwJh6URXX+8X46mXMOzI5ma8LHsBoyDumzpS2Ilot6dprpZe+NKyMsWRJuJ+pbBiArJmVV1yR/dwDDmDqGTCO6H7Pa3ZWeuUr9x/kXLyYqWwYmKxu8G99S7rqquTnNSm9g9WNgfzofs+j1QoLU6dlLXV2u3MWQgUqXs55KFjdGAhYerVMn/lMdvmtpUtDk+ltb+MshN8Y9LVd09M7KAkB7MU89TJ95SvZjy8shL7OaCmu+FkoakatXMlZaIwkXdudf36513Zp88x7nXpWtw4mSkIAxRHUyzAxEQY7OQtB/V/bFQmuZU09q+IipKhhl4So20UOkAdBPY83vlHauDH98cWLwyf+ox+lMA36urbrJbj2u5xzXTuYBlESIm+gruNFDgasrKu4YV8NuvvIfb3sZS/zSrVa7r/1W+6hCuf+X9/8ZtjuC19wn5pK3mZqKjyOxlu71t0s+TAwc1+3Lvl5O3a4T08nP296OhyG0XZf+EL4O1/4wt77e1XXw7bVyvf/yOuuu8LzpqbC+zA1FX6+6659t8v7PqBB8h4cVf2eBJLmPEd8ZJ56HkuXSl//esjOieamT0yECh/f/Obe1TNYYx3KXmY1q4WZp4U/OxuSx1avltasCbfLloX7ezXsbu40eVZAziveGzE/H0L0/Hz4eeXKfZetzfM+oEGKHBxV/J4+EdTzigYvr746LGB9zTWhZNfLX753YeuPfjSssV7GWQgjq9dru27B9f77B3PO6PUipArRx279emnt2nC7bVvxLvAigbquFzkYkLKu4mpyNciYehFJg5dJg2+LFoWCNGYUzR5DvWaldxtD/uUvB5OHmbXWex06mHrJGegc1rzvvvyBmuUdxkxZV3E1uRokqPcjK8PoYx9jCtsY6yUrvVtwfeYzi50z8ubrZF2EbNoULhZGKQM86Tr76afDfrda+2/fGajrfpGDkpV1FVeTq0GKz/RjHEp6oVJZtYseeCD/4dZLDaToIiDPWu91zQDPKlhjlrwOfVIhG2pIjZGyqjgNuBoUxWeqUJPuFjRHVgv/xS/O14LsdYpavJs7KTjWYZpbN1nDmgcdFG4POKD7kAhLz46Rsqo4lV0NqkcE9X7UpLsFzZI2hhw/ZywshK+lS/feH50zyqiBNGp1lKJehquuSl9/vtWS/vZvpeOPzxeo+53/jxFS1lVcDa4GCer9YPANFWq1pFtuCbd79oT7ohnUcXk6kLqNt49SJ1S8qzwtoEvhOvv44/ftjRi1fAEMUFlXcUO+GiSo96Mm3S1ovihwdQ7X7dwZvuJd4t06kNxD13pWtbRR6YRKGmpIE7/OpmIcmopEuTJ0ZhhxyY8SZSV/ReKJct3yddylJ59Mfiy6MBiVFeCyclUjU1P7Jrmx+htGUd5EOYrPlCHqbrnkknDLGQElyhrfjmzfHsaTN2wIP6dVYrvggr1d953i9THKrOY2SFnDBGbSq1+9f8GamtQIAQaC7negBrLGuLMCV8Rduu026a679nYjx/N1jjwybPe5z6W3ajvHymuQ89NVt2GCv/iL/Yc3RylfACiKoF62Ya/Qg5HTbXw3K3B16px29ra3FUsk6xwrr3sGeC+5qqOSLwD0gjH1MlGxAgXlGd+V0se300Rj7Gec0X08vvPvFb0GHfZ1bNGP3ajkCwBxFJ+pWl0XpUat5Z0P3jnJYsmS8H3ac6Nu5Dzj8fFEsqKHaB2yyIsOEzBpBU1GUC/LqFXrQC3kHd9NCly7dknveU92N3K3RLJXvSqMO/fSuq7TdWzRYYJRyBcAekFQLwvZN+hBkfHdzsC1sCD99V8n/95oPPn664snkuU1qOvYqrrz654vgCEY9lhSCQjqZSH7Bj3opyhhnm7kVaukv/qr3n5/N0WuY/OeK9O68y+4IDxe1Xm2Aed2FFWHsaQSkChXFrJv0KN+8yuzah/Nzoa52r/+9b7PmZoKU+D6OVflXaQw7+vLW2Rn0Lmn5LuOoRGoSJQ3UY6gXibOBujRIIoSDvo8lec6NipJm2cf8lSHK3P/k4zAuR2DMALLaFNRbhii7Jv166W1a/cvZQWkGERRwkFXTstTda7IPuQpslPm/ieh2tyYalBOVClj6mZ2iqR/lLRI0pXuflnH40skfVbSyyQ9JulP3f3H7cc+KOmdkvZI+q/ufmsZ+5TLIAbOyL5BTVRxnuqWRV67lxq9AAAgAElEQVRkH4oU2RnUebZB53YU0aCcqL5b6ma2SNInJb1e0gslvdXMXtix2Tsl/crdf1fSxyX9Q/u5L5R0lqT/IOkUSf9f+/cN3uxs6GdbvVpasybcLlsW7gcaIDpPJSnrPBW/Lk4aNiiyD6tWhdGqPAZ1nq3if4Yayjr4RmwZ7b7H1M3sDyV9yN1f1/75g5Lk7n8f2+bW9jb/bmYTkn4u6XBJF8W3jW+X9Tf7HlNn4AxjYNC5m3lSSIruQ+fvTDs9DepjSr7rGKt5TlSVFeWOkPRw7Octkn4/bRt3321mT0p6Zvv+Ozuee0TSHzGz8ySdJ0krVqzob48pFIMxMMjKaXkLzxTdh87ufHfp8svDynJVVH6j2twYa0hFopGZp+7uV0i6Qgot9b5+GQNnaLDOVJHNm6Wbby73PFXkurjoubIzLeXCC6s9zzbk3I5eNCAnqoygvlXSUbGfj2zfl7TNlnb3+zMUEubyPLd8DUqKAOKyehDLPFd1uy7uzD/t51wZPTe6WPnoRwdfEKYB53b0a0QrEJUxpj4h6YeSXq0QkGclvc3dvxfb5nxJL3b3vzSzsyStcve3mNl/kLRB0omSlkv6uqTj3H1P1t/se0ydgTM0UJWpIt3mlC9ZEmrTlzUsWfPhTjRNDQ+4yuapu/tuSe+WdKukH0j6krt/z8zWmtlp7c2ukvRMM3tQ0gXamyD3PUlfkvR9SbdIOr9bQC9Fngm2wIipco51t0z1nTvDePj8fLjIWLkyXEv3Ij5+Pz9f3u8FEo34AVfKmLq73yTppo77Lo19vyDpT1Ke+2FJHy5jPwph4AwNEO8h/M530lvOZaeKJCWULVmSfr7rJ/+UvFZUqugBV7Nu+pFJlBsIBs5QE72cFzp7CJcsSd92EKkindfF99wjbdyYvG0/FxXktaJSRQ64Gi4CM95BXco+m9bsCgzN1Mt5IWlKWVav4KDqZ8SvizdskG65pfz802HltfLxH1N5D7i88zqr5u4j9/Wyl73MS3HXXe7T0+5TU+5m4XZ6Otyf9RhQkh07wmEVBu72/Zqedm+1kp/3hS+EQzLpeZL7kiXVH7atVm+vZVi/Nwsf/zGW94DL+hBOTYXHSyRpznPEx/FtqXe7ynKXnnwy+TGy41GSXseLs3oIzaRTT5X+03/KnypSRqt0UIVbqi4I8/jj0sknSzt27L2Pj3/DdX4AbrwxfIiyDriajguNb1DPOpsuLKTXpyQzByXq9bzQrYfwjDPyH6JlDgsOKv+0qrzW2dn9A3ocH/8GSpu+tmmT9PDD6QdcTeudjG9QzzqbZg1Ozs+HMwsDbChBr+eFVatC4E1SZPx8EMOCg8o/7fZ7++1tiP4XaQFdIjGvcbI+AKedlv0BKOtDWLLxXU89azmmpUuzU4lvvJEV3VCKXheHKqvUQlPWDy9j0cWs/0WEgpMN088HoKb1Tsa3pZ51lbV0aeh+37kz+fGdO8MXA2zoUz/jxWV0SRft/q9jRnhZvQ1Z/4vIiK3CiW76HRevYb2T8Q3q3c6m0t7qQYOoqAG09XNe6OySbrXC1LK8QbdI938Np+RKKq84Tdb/QpImJyk42ThljIvXrN5J37Xfh6Hv2u9xUdMj6WzaaknveEd6RQ0zae1a6ZJLytkXIIe01nIv5arzLoPQra58tBLcMFrw69aFLvekU1mRj+jCgvQ7vyP9+tf7P7Z0afhfTE/3v7+okRFaB6TK9dRHW9ZV1kEHhTTiQVTUQGMNsos6rbW8aVPI6ynaBZ23+z+rNbxjh/SsZ0kHHBBGpaamqm3Bl5WEHE0yTjIxUZtzO8pU9XzJCtBS72aEruQwfINc3Cmrtbx4cbjdtWv/x6ampPXr82WOp3X/Z7WG00QfD/fBjsOX9RHNWnkuz/8QI6zbB6AGaKmXpYFXchiMQVeNzGotJwXzSJ58n27Dgt3Gm5Ps3h3WPv/Yx/obh+/W81HWR7SmtURQhZqNi/eDoN5NqyU98ID0nvdIjz0mHXaYdPzxtbySw3B1S9i69trQou61xZonOztJGaNEq1aFVmwR8/PS3/3dvnmmeS5y4kFcki6/XNqzJ/uioIwk5JrWEgEKIahnyepLJaCjQ7eW3rveFQ6fXlusvbSWpXKmYbmHpLMisko9pGWld37kOrv7sy4K+m1s1bSWCFDI+Baf6Sbelzo/H84u8/Ph55NPTh7Aw1jLqmfkHlqsnYdSNGsyj6xCNUkWLy6vDsYNNxQbT5fCRUBaqYe0OfCdH7k0gyiMU9NaIkAhBPU03dJ9jziCinLYR9GgKxULTp1BJ8vixdK554bWbBkZ6EW7/p/xDOm//bf0i5yk7uw8Fd0igxrjjrrx168PU+HWry/vfwhUgaCepttZbMeOYs0sNF5aSy/L/Lx0//35/0Y86Fx6aSiIkmRyMoxFl9W6zOqFiCxdGr7WrpV+/nPpb/6mWAncIhcOgxrjjo/n1zQJGsjEmHqaPAOYVJRDh86ErW3bpKuuSu+GlqT/8T+kN70pf2swPnb8pjdVk/aRNd48OSldeGFy/miRrPQiOQODGOOua8U8oAiCepqss1iEeS7I4C498kh2QJf2dvr0MuWtqtLT3aaNpQW9IvuX5yM3NTWYi5ZBT0cEqkJQTxOdxbIWV2aeCzp0tvayMsDj+un0qWqKba8XEHn3L37h8NRT4f+3eHEI4u97X/h+UBctZdWPB4aNoJ7lhBOkLVukI49MDuzMc0FMUmsvb8pFVqdPfJz3qKPCfQ8/PJxV0gZ9AXHCCaHk7SmnhCAeFdW5/HLp1lull798MH+XwjNoCoJ6N4ccIt1+e3K/46ZN4RK+TutQYmi6ZW9PTKQ/ntbpkzVvu+oa61VotUIN+/g19K5d4eukk6Q77hhMYKfwDJqCoJ5HZ7/jkUdKP/pRONtG66438QyLQrJae2bSqaeGlcySWu9JnT5JLf+4fsZ867guupR9YeQuve51obBj2ftK4Rk0BUE90u0sF/U7Rk2nzjMtWTVjL6u1t2SJtGhRmLvdWQs9LfEr77ztomO+dc7yfuih4Uw4YYkHNAVBXcp/luvWdIq2ufZa6ZxzBr7bqJes1t7CgrRxo3TjjeHnU0+VXvjC7GUE8s7bLjLmW/cs7+c+d9+x9E67dg1ufLuqmQTAIFF8JqscbGdxmTxNp507Q5Fvqs2NnaTiM5127gxfGzdK//RP0nHHpQeNPAVfpGJjvnmyvIepW1W+QY9vRx1yl1wSbgnoGDUE9SJnubxNp4UFqs2NqXjFtz/+4+wpbd1qv+ctO1tkzLfuWd5Ll4Ys97TFYw48kPFtIAtBvchZLm/TSapHswdDEbX2XvCC7oVnsg6Tbi3/XhYbyTqEDzwwdL8P+1r05S+X/u3fQqW6xYvDfQcfzMIqQB6MqReZy5Kn5FWkDs0eDM3sbJhb3U23w+SEE8J159/+rXTffdLv/q40MyP94he9jflmHcK7dkmf/ay0YcPwk+b+83+WfvlLxreBosyLrqdYAzMzMz43N1fOL1tYkJYtS05+m57eP3MonlSXlaY7NRX6YClDNXZaLWn58nyr83Y7TL75zVCIZffuEHQPPji0qPsJutEhHFVtS5J06AMYHjO7291num1H93vRRZTzLpPF5NaxVWQJ0azD5JvflF7xilCIJcoG3769+DrsnaJD+Oyz93Zvd2L0CBhNdL9LxeeyDGOZLIyMPPmU3RYmabVCCz2tI+2pp/qbr33QQdKznx1+TxJGj4DRRFCPpBW17laUhsmt6NCtCM2pp0pnnJF9mHRr7ZcRdCmNCjQPY+pZkorSZK01Wdfam6hU0TSNJOvWhdGdNIsXS5/+dH8pG2XsJ4Bq5B1Tp6WepmjprTrX3kSlyig5+tznhi76tFzMXlM2Oq87b7wx9BwwegQ0A0E9TZEFluteexOV63dUJmvqmVko0FL0kErreNq0KSzlyugRMPoI6mmKFKUpcgGAsdHP2uPx1n409WzRovB1003Flx/Nuu487TSuO4GmYEpbmqzSW51ZRHWvvYmRdMIJoRXtHsbQ9+wJt2eeWXxpgbrXfAdQDoJ6mqzC250DmkUuAICcWq3Qio7PU09ba6gbrjuBtlYrlE1cty7cDrsucskI6mmKFKUpcgEA5FRm65rrTkChi2v5cmn1amnNmnC7bFmjVtVkTD1L3mynMtKdgQ7dWtdXXRW+z5PYlpV4V4frzigr/777pMcekw47LHuteaCwMUlo7iuom9mhkq6TdLSkH0t6i7v/qmOb35P0KUm/LWmPpA+7+3Xtxz4j6RWSnmxvfo67f7uffSpd3mwnitCgZFnFYdyl226T7ror38zJOl93Rln5O3eG825kcpJZoSjRmCQ099tSv0jS1939MjO7qP3zBzq22SHpHe7+gJktl3S3md3q7tHl0t+4+8Y+96Me+kl3BjrkWRSwSEOjjtedSY2nyI4d4atBjSgM05gklvQb1E+XdHL7+2sk3a6OoO7uP4x9/zMze0TS4ZJyrGEFjK+k1nVaAci8DY2i152DLpKYZ/GbBjWiMExHHRWmj+zcuf9jDUos6TeoP8vdt7W//7mkZ2VtbGYnSlos6aHY3R82s0slfV3SRe6e8B+XzOw8SedJ0ooVK/rcbWA0xFvXV10VutyTDKKhkVUk8UUvKifY51n8pkGNKAzL7Kz013+dHNCleiSWlKRrUDezr0l6dsJDF8d/cHc3s9RC8ma2TNLnJJ3t7k+37/6gwsXAYklXKLTy1yY9392vaG+jmZmZagrWU8sdNRBvXd91VzULsGTlFL3qVeEc2E9F5Oij9Z3vhEVusmYVNagRhWHIGuORpGc8Y/iJJSXqGtTd/TVpj5nZL8xsmbtvawftR1K2+21JX5F0sbvfGfvdUSt/p5l9WtKFhfZ+kKjljpqpMoM9q1u886KiaAJx/KOVVts+rkGNKAxD1sG8dKn0iU806pze7zz1TZLObn9/tqQvd25gZosl3SDps50Jce0LAZmZSXqzpO/2uT/liF/Zzc+Hgcxeq34AJSlSOqFfebrFO+WZO9/50coyOTmY14Yxk3Uw79wpbdlS7f4MWL9j6pdJ+pKZvVPSTyS9RZLMbEbSX7r7ue37TpL0TDM7p/28aOraF8zscEkm6duS/rLP/SnHmEx9wOipKoM9azpdmjxj31kfrSVLpDe8IdQCOfxw6XnPY8QLJcg6mBs4ttNXUHf3xyS9OuH+OUnntr//vKTPpzz/Vf38/YEpMvWBcXdUrOyZk0mHcJ7pdJ3ynB+zPlq7dkkvfal0ySXF/i6Qqe6Vl0pGRbkkea/sGHfHiMs6hJOK1SxaJD39tPTrX+//u/KcH8es0YQ6qHPlpQEwT5v4WmMzMzM+Nzc3uD+wsBD6AJOyJaenQzaQe6ghnLVNww4WNEur1f0Qdt+/q//ee9PPj92uZfN8tPjYYCCiLqm6VF4qyMzudveZbtvRUk+S58puwwbG3THS8qaOdB7GSeP6r3+9dPPN0i23ZI9CjVmjCXUyJhU/CeppumUkjUnJQTRXP4dw/Pw4Oxs+HnlHoepYrhZoCoJ6lqwrOwYHMeLKOIR7XfhqTBpNQOVYT71XrKGOEVfGIVzmmu8A+kdQ71WVlUCAASjjEGYUCqgXut/7weAgRly/hzCjUEC9MKUNQM+YogZUgyltdUUFOjQIU9SAeiGoV4kKdGggRqGA+iCoV6XXuT/ACGCKGlAPZL9Xhbk/ADAcrVaoArpuXbht8PLZtNTLkGecnLk/AFC9MRv2JKj3K+8Bw9wfAKjWGA570v3ej/gBMz8flrSanw8/r1y5bxcPFegAoFpjOOxJUO9HkQOGCnQAUK0xHPak+70fRQ8Y5v4AQHXGcNiToN6PXg4Y5v4AQDVWrQo5TkkaOuxJ93s/yhwnH6MpFwBQiTEc9qSl3o9+a2RGU+HuuEP63OekAw4YiykXAFCZMRv2ZEGXMkTBucgBE02Fe+qp9HF5VsQAAIgFXapVdJw8ae5kkiiDnjF4AEAOjKkPQ9ZUuLiGTrkAAAwGQX0YsqbCxTV0ygUAYDAI6sMQTYXrpqFTLgAAg0FQH4asqXBS46dcAAAGg6A+DElzJw8+WJqclFavltavD1nvTGcDABRA9vuwjNncSQDA4BHUh4mSsQCAEtH9DgBAQxDUAQBoCII6AAANQVAHAKAhCOoAADQE2e8AgPETra750EOhymdDphQT1AEA4yVa+nr37rAOx8EHS+efH4qCjXjRL4J6nTT0yhEAaiNp6ev5+XC7cmWo5jnC512Cel00+MoRAGoja+nr3bul668f6aJgBPU6aPiVIwDURtbS19u3h7LdI4zs9zrIc+UIAOhf1tLXBx8c1uEYYbTU66DIlSPj7gDGVbfzX57z46pVYWgzycREeHyE9RXUzexQSddJOlrSjyW9xd1/lbDdHkn3tn/8qbuf1r7/GEnXSnqmpLslvd3dd/WzTyMpunKMutzj4leOjLsDGFfdzn95z4/R0ted205MhPtHvJFk7t77k80+Iulxd7/MzC6SdIi7fyBhu3l3n0q4/0uSrnf3a83snyXd4+6f6vZ3Z2ZmfG5uruf9rp2FBWnZsn3H1CPT02FM3V1avjx7mxE/GAEgUauVff6Llq9OenxyUnrf+6TnP3/flnvUqh+Rpa/N7G53n+m2Xb9j6qdLuqb9/TWS3pz3iWZmkl4laWMvz2+U6MpxelqampLMwu309N4rR8bdAYyr664LjZ8ku3dLl1ySfn7csUNat05avTo0nmZnw/3R0teXXBJuaxzQi+h3TP1Z7r6t/f3PJT0rZbulZjYnabeky9z9XxS63J9w9+id2CLpiD73Z3SdcIL0s5+lXzk2PGMTABLNzkrveld6UN++XbrvvvTzY2RMZhR1Depm9jVJz0546OL4D+7uZpbWl/8cd99qZsdKus3M7pX0ZJEdNbPzJJ0nSStWrCjy1NERXTkmyTvuDgBNEU33TQvoUjj/Pf/50l13JZ8fOzVgLnqWrt3v7v4ad39RwteXJf3CzJZJUvv2kZTfsbV9u1nS7ZJeIukxSdNmFl1YHClpa8Z+XOHuM+4+c/jhhxd4iQ2xalVI5EjSgIxNANhP1rBjZGJC+vCH08+PnRres9nvmPomSWe3vz9b0pc7NzCzQ8xsSfv7wyT9kaTve8jQ+4akM7Oej7Y84+4A0CRZw47SvufFzvNjmob3bPYb1C+T9Foze0DSa9o/y8xmzOzK9jYvkDRnZvcoBPHL3P377cc+IOkCM3tQYYz9qj73p9micff166W1a8Pttm1MZwPQTFmFYpYskT71qb3nv+j8+P73S4sXp//Ohvds9jWlbVgaN6VtEChSA2BUpJ2v8kz37Sw+kzb1LXrOiNb1yDuljYpyTUSRGgCjotv5qkihmKwx+CVLpI9/vPHnQIJ607A4DIBRked81W26b1zWGPyuXdKWLeW/hpohqDdNw5cVBNAgec9X8em+rVa4P2lokam/BPXG6VakhnF2AHVRtKhWt676hi/WkgdLrzZNVraou3TjjdKaNfuXTASAqhVZBjXeVT8/H85n8/Ph56hADVN/CeqNk1WkRpJ27kz+MABA1YoU1cq7/sWYT/2l+71pkpYVXLIkezEExtkBDEORZVCLdNVnldxuOIJ6E3Vmi95zj7RxY/K2DS+ZCKDm8ma3kwSXC0G9qeJXqhs2SLfcwocBQD3laVmTBJcLY+rjgMVgAIw6kuByoaU+DoqMWwFAXRUpRDOmCOrjgg8DgCYY4yS4PAjq44QPAwA0GmPqAAA0BC31ccJyrADQaAT1ccFyrADqgMbFQJm7D3sfCpuZmfG5ublh78boaLWk5cv3Xd4wMj3NcqwAqpHUuIhm4dC4yGRmd7v7TLftGFMfB1k1kxcWpLe/PRSooQY8gEHJsyAL+kZQHwdZNZMXFkIJWVZtAzBIeRdkQV8I6uMga3nDCFfMAAap6Nrp6AlBfRx0W441jitmAINQZO109IygPg46ayZn4YoZwCCwBkUlCOrjIioTu369dOaZYY31JFwxAxgEFmSpBFPaxtHCQkiKS5riNjkpbd0aPmgAULZonjprUBSSd0obxWfGUXzVtp07w4cs7phj9p83SsEIAGVgDYqBoqU+zh5/XDrqKGnHjv0fixeloWAEAAwVxWfQ3S23SAekHAJRFjwFIwBgZBDUx1meeaMUjACAkUFQH2d55o1SMAIARgZBfZzlmTdKwQgAGBkE9XF2773JXeu/9Vt7541SMAIARgZBfVxFCXDz8/s/dsAB0otfHL6nYAQAjAzmqY+rrAS4PXtCAlw0lzSqRkfBCAC9oM5FZQjq46poAhwFIwD0IqnOxfnnU+diQAjq4ypKgEvqficBDoDUfws7XuciEp1zVq7cW+AKpWFMfVyRAAcgy+ystHy5tHq1tGZNuF22LNyfF3UuKkdQH1ckwAFIU1YlSepcVI7u93FGAhyAJHla2HlybBjmqxxBfdyRAAegU1kt7FWrQlJcEob5BoLudwDAvsqqJMkwX+VYehUAsK+FhZAUF89aj8SXZc4ryqJnmK9neZdepfsdALCvqIXdOb98YqK3FjbDfJXpK6ib2aGSrpN0tKQfS3qLu/+qY5tXSvp47K7nSzrL3f/FzD4j6RWSnmw/do67f7uffULFqBQFNBOJtCOpr+53M/uIpMfd/TIzu0jSIe7+gYztD5X0oKQj3X1HO6j/q7tvLPJ36X6viaRKUdGVPJWiAKA0ebvf+02UO13SNe3vr5H05i7bnynpZnff0effxbCVNY8VAFCafoP6s9x9W/v7n0t6Vpftz5L0xY77Pmxm3zGzj5vZkj73B1WhUhQA1E7XMXUz+5qkZyc8dHH8B3d3M0vtyzezZZJeLOnW2N0fVLgYWCzpCkkfkLQ25fnnSTpPklasWNFttzFoZcxjZTweAErVNai7+2vSHjOzX5jZMnff1g7aj2T8qrdIusHdn4r97qiVv9PMPi3pwoz9uEIh8GtmZmb05uE1Tb+Voli5CQBK12/3+yZJZ7e/P1vSlzO2fas6ut7bFwIyM1MYj/9un/uDqvSzIAzj8QAwEP0G9cskvdbMHpD0mvbPMrMZM7sy2sjMjpZ0lKR/63j+F8zsXkn3SjpM0n/vc39QlX4qRTEeD5Sv1ZI2bJDWrQu3XByPpb7mqbv7Y5JenXD/nKRzYz//WNIRCdu9qp+/jyHrdR4rKzcB5co7nEUeS+NRUQ796aVSFCs3AeWJD2dFos/WypV7S7qSxzIWWNAF1cszHk9XIpBPnuEs8ljGBkEd1es2Hn/vvdLy5dLq1dKaNeF22bLQ0gCwrzzDWeSxjA263zEcaePx7iGgd+tKjGOcEOMsz3AWeSxjg6CO4Ukaj9+woXuLIv6cfsYJuRhAE6xaFY75JNFw1vXXk8cyJuh+R70UaVH0M044O0sXP5ohz/TSfupKYKTQUke9FMmMzzNOmJSZnzdbGBgV3aaXlr0+OmqLoI56ydOVGOl1nLDXiwGgzrpNL2V99LFAUEe9FGlR9DrfnaQhjKte6kpgpBDUUT95WxRFWvVxFL8B0FDmPnoLns3MzPjc3NywdwN1kJT9HrXq07LfFxZCUlx8TD0yPc2YOoDaMbO73X2m23a01DHaehknJGkIQEMR1DG6OueZX3hh/oBM0hCABiKoYzSVsTgFSUMAGoagjtHDPHOMIioYogIEdYwe5plj1LDsKSpCUMfoYZ45Rgk9S6gQtd8xeqJ55kmYZ466YdlTVIigjtHD4hQYJfQsoUIEdYyePKtSAXWR1bPkHr6AklBRDqMryiZmnjnqLKuCoSQ94xnSz3/OsYtMVJRD8zHPvHmaOO1r6VLpggukSy9NfnzPHmZsoDQEdQD1MK7TvhhXR4kYUwcwfPFpX/PzYZx5fj78vHJl6MIeZc99bsj7SMKMDZSIljoQ6db128Su4Srk+b81vaBQr8sEAwUR1AGpe9fvuHYN9yvv/63p075YGRAVIagD3Sp+bd5MRbBeFKmkFk37ih6Pa0r3NCsDogIEdaBb1+8llwyva3iUu/yLdKmPS/c0MzYwYAR1oFvX7333DadreNS7/O+7L7nlLe3/fxvV7ulRvuhCIxHUgW5dv89/vnTXXdV2DY/6IiCzs9Lll6c/nvR/G7Xu6VG/6EIjUVEOyKr4NT0t/ehH0jHHpD8+iAC7YYO0enXyhcTUlLR+fX27cVstafny9Apq0uD+b1XJeo2j/tpQS3kryjFPHehWSz5+W1Wt+Tpkg7da4eJi3bpwm3eueNZYuiRNThb7v/W6H4PEymuoKbrfAal712/VXcPDzgbvp2s564JEki68MH/3dF27uOtw0QUkIKgDkW6ZyVVmLg8zG7zf8fysC5KpKen446vZj0Ea9kUXkILud6COOocEJGnJkr2LgwxSv13LZa13X+cu7rJeI1AygjpQV1GX//vfH4K5mbRzp/SRj4TEvtnZwfzdfruWy1rvvs5d3GW9RqBkdL8Ddfexj+2bHDboLugyupbLyEGoexf3qE3Bw1hgShtQZ8OY2tZtil9VY9l12Q+gBpjSBjTBMLqg69K1XJf9AEYI3e9AnQ2rC7ouXct12Q9gRND9DtQZXdAARPc70Ax0QQMogO53oO7oggaQU19B3cz+RNKHJL1A0onuntgnbmanSPpHSYskXenul7XvP0bStZKeKeluSW9391397BPQSKzDDSCHfrvfvytplaQ70jYws0WSPinp9ZJeKOmtZvbC9sP/IOnj7v67kn4l6Z197g8AAGOrr6Du7j9w9/u7bHaipAfdfXO7FX6tpNPNzCS9StLG9nbXSHpzP/sDAMA4qyJR7ghJD8d+3tK+75mSnnD33R33JzKz88xszszmHn300YHtLAAAo6rrmLqZfU3SsxMeutjdv1z+LiVz9yskXSGFKW1V/V0AAEZF16Du7q/p829slXRU7Ocj2/c9Jt2UWOEAAAcdSURBVGnazCbarfXofgAA0IMqut9nJR1nZseY2WJJZ0na5KHqzTckndne7mxJlbX8AQBomr6CupmdYWZbJP2hpK+Y2a3t+5eb2U2S1G6Fv1vSrZJ+IOlL7v699q/4gKQLzOxBhTH2q/rZHwAAxhllYgEAqDnKxAIAMGYI6gAANARBHQCAhiCoAwDQEAR1AAAagqAOAEBDjOSUNjN7VNJPSvyVh0n6ZYm/b5h4LfXEa6knXkt9Nen1lPFanuPuh3fbaCSDetnMbC7P/L9RwGupJ15LPfFa6qtJr6fK10L3OwAADUFQBwCgIQjqwRXD3oES8VrqiddST7yW+mrS66nstTCmDgBAQ9BSBwCgIQjqAAA0xNgEdTP7EzP7npk9bWapUwvM7BQzu9/MHjSzi2L3H2Nm32rff52ZLa5mzxP38VAz+6qZPdC+PSRhm1ea2bdjXwtm9ub2Y58xsx/FHvu96l/Fb/az62tpb7cntr+bYveP2vvye2b27+1j8Ttm9qexx4b+vqQd/7HHl7T/zw+2/+9Hxx77YPv++83sdVXud5Icr+UCM/t++334upk9J/ZY4vE2LDleyzlm9mhsn8+NPXZ2+5h8wMzOrnbP95fjtXw89jp+aGZPxB6r2/tytZk9YmbfTXnczOyf2q/1O2b20thjg3lf3H0sviS9QNLxkm6XNJOyzSJJD0k6VtJiSfdIemH7sS9JOqv9/T9LetcQX8tHJF3U/v4iSf/QZftDJT0uabL982cknTns96TIa5E0n3L/SL0vkp4n6bj298slbZM0XYf3Jev4j23zV5L+uf39WZKua3//wvb2SyQd0/49i2r+Wl4Z+0y8K3otWcdbjV/LOZL+34TnHippc/v2kPb3h9T5tXRs/x5JV9fxfWnvz0mSXirpuymPv0HSzZJM0h9I+tag35exaam7+w/c/f4um50o6UF33+zuuyRdK+l0MzNJr5K0sb3dNZLePLi97er09j7k3ZczJd3s7jsGule9KfpafmMU3xd3/6G7P9D+/meSHpHUtUpURRKP/45t4q9xo6RXt9+H0yVd6+473f1Hkh5s/75h6fpa3P0bsc/EnZKOrHgf88rzvqR5naSvuvvj7v4rSV+VdMqA9jOPoq/lrZK+WMme9cDd71BoMKU5XdJnPbhT0rSZLdMA35exCeo5HSHp4djPW9r3PVPSE+6+u+P+YXmWu29rf/9zSc/qsv1Z2v+D8eF2d9DHzWxJ6XuYX97XstTM5szszmgYQSP+vpjZiQqtlYdidw/zfUk7/hO3af/fn1R4H/I8t0pF9+edCi2qSNLxNix5X8sft4+djWZ2VMHnViX3/rSHQ46RdFvs7jq9L3mkvd6BvS8TZfySujCzr0l6dsJDF7v7l6ven35kvZb4D+7uZpY6L7F9VfhiSbfG7v6gQtBZrDB/8gOS1va7zxn7UMZreY67bzWzYyXdZmb3KgSUSpX8vnxO0tnu/nT77krfFwRm9meSZiS9Inb3fsebuz+U/Btq4UZJX3T3nWa2WqE35VVD3qd+nSVpo7vvid03au9L5RoV1N39NX3+iq2Sjor9fGT7vscUuk0m2q2T6P6ByXotZvYLM1vm7tvaweGRjF/1Fkk3uPtTsd8dtSZ3mtmnJV1Yyk6nKOO1uPvW9u1mM7td0ksk/W+N4PtiZr8t6SsKF5t3xn53pe9LgrTjP2mbLWY2IekZCp+PPM+tUq79MbPXKFyQvcLdd0b3pxxvwwoeXV+Luz8W+/FKhfyO6Lkndzz39tL3ML8ix8lZks6P31Gz9yWPtNc7sPeF7vd9zUo6zkJG9WKFg2qTh8yGbyiMTUvS2ZKG2fLf1N6HPPuy35hUO+BEY9JvlpSYuVmRrq/FzA6JuqLN7DBJfyTp+6P4vrSPqxsUxtk2djw27Pcl8fjv2Cb+Gs+UdFv7fdgk6SwL2fHHSDpO0l0V7XeSrq/FzF4iab2k09z9kdj9icdbZXu+vzyvZVnsx9Mk/aD9/a2SVrZf0yGSVmrfXruq5TnGZGbPV0gg+/fYfXV7X/LYJOkd7Sz4P5D0ZPvifXDvS1lZgHX/knSGwrjFTkm/kHRr+/7lkm6KbfcGST9UuPq7OHb/sQonqQcl/S9JS4b4Wp4p6euSHpD0NUmHtu+fkXRlbLujFa4ID+h4/m2S7lUIGp+XNFXn1yLp/2nv7z3t23eO6vsi6c8kPSXp27Gv36vL+5J0/CsMAZzW/n5p+//8YPv/fmzsuRe3n3e/pNcP630o8Fq+1j4XRO/Dpm7HW41fy99L+l57n78h6fmx5/5F+/16UNKf1/21tH/+kKTLOp5Xx/fliwozWJ5SiC/vlPSXkv6y/bhJ+mT7td6r2MyrQb0vlIkFAKAh6H4HAKAhCOoAADQEQR0AgIYgqAMA0BAEdQAAGoKgDgBAQxDUAQBoiP8f+0KlT++nnVEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.io\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "mat = scipy.io.loadmat('hw5_p1b.mat')\n",
    "X = mat['X']\n",
    "k = 2\n",
    "num_points = int(X.size/2)\n",
    "# Guess the mus\n",
    "changed_points = 1\n",
    "num_dims = 2\n",
    "sigma = 0.2\n",
    "kernel_factor = 1/(sigma*sigma*2)\n",
    "\n",
    "\n",
    "cluster_assignment = np.random.randint(0, k, num_points)\n",
    "\n",
    "def gaussian_RBF(point1, point2):\n",
    "    diff = point1-point2\n",
    "    # this diag thing feels like bullshit\n",
    "    stuff_to_put_up = -kernel_factor*np.diag(diff.dot(diff.T))\n",
    "    return np.exp(stuff_to_put_up)\n",
    "\n",
    "def compute_new_assignment(points, assignment_vector, num_classes):\n",
    "    assignment_matrix = np.zeros((num_points, num_classes))\n",
    "    assignment_matrix[range(num_points),assignment_vector] = 1\n",
    "    class_counts = np.sum(assignment_matrix, axis=0)\n",
    "    \n",
    "    # Precompute all pairwise kernels\n",
    "    kernel_results = np.zeros((num_points,num_points))\n",
    "    for i in range(num_points):\n",
    "        kernel_results[i,:] = gaussian_RBF(points[i], points)\n",
    "    \n",
    "    distances = np.zeros((num_points, num_classes))\n",
    "    for i in range(num_classes):\n",
    "        ass_vec_class = assignment_matrix[:,i,None]\n",
    "        first_term = 1\n",
    "        second_term = -2/class_counts[i] * kernel_results.dot(assignment_matrix[:,i])\n",
    "        third_term = 1/(class_counts[i]**2) * np.sum(ass_vec_class.dot(ass_vec_class.T) * kernel_results)\n",
    "        distances[:,i] = first_term + second_term + third_term\n",
    "        \n",
    "    return np.argmin(distances, axis=1)\n",
    "\n",
    "steps = 0\n",
    "while(changed_points != 0 and steps < 200):    \n",
    "    new_assignment = compute_new_assignment(X, cluster_assignment, k)\n",
    "    changed_points = np.sum(new_assignment != cluster_assignment)\n",
    "    cluster_assignment = new_assignment    \n",
    "    steps += 1\n",
    "    print(steps)\n",
    "        \n",
    "red_X = X[cluster_assignment==0]\n",
    "blue_X = X[cluster_assignment==1]\n",
    "\n",
    "plt.figure(figsize = (8, 8))\n",
    "plt.plot(red_X[:,0], red_X[:,1], 'ro', ms = 7)\n",
    "plt.plot(blue_X[:,0], blue_X[:,1], 'bo', ms = 7)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}