{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-9673c77d0e8030f6",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "# K-means: Profiling and optimization demo\n",
    "\n",
    "Here is a more complex example, which shows a simple implementation, diagnosis, and reimplementation development workflow. In short, it walks you through the following:\n",
    "\n",
    "1. A brief mathematical summary of Lloyd's algorithm for the k-means clustering problem.\n",
    "2. An initial implementation, which was written by a student in [Georgia Tech's MS Analytics program](http://analytics.gatech.edu) as part of a homework assignment.\n",
    "3. An invocation of Python's `line_profiler` tool to identify potential execution time bottlenecks.\n",
    "4. A reimplementation of one those bottlenecks in Intrepydd.\n",
    "5. An optional exercise for you to try improving a different bottleneck in the same code.\n",
    "\n",
    "This example follows an earlier [toy problem](./003-profiling.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-541e8cffa94a67c8",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**Scaffolding.** This example includes some auxiliary scaffolding code to help generate a sample dataset and visualize it, so you can focus on implementation basics. Take a moment to skim and run these cells now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import pandas\n",
    "\n",
    "from auxiliary import gen_mixture_of_gaussians, visualize_clusters, count_matches"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Dataset.** The following cell generates a sample dataset, which is a two-class mixture of Gaussians."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 200000 2-dimensional points.\n",
      "There are 2 classes: [0 1]\n",
      "\n",
      "Coordinates of the first few points:\n",
      "[[-0.38123772 -0.51561695]\n",
      " [ 0.69589787  0.26799208]\n",
      " [-0.2817528  -1.09621906]\n",
      " [ 0.53102266 -0.16483183]\n",
      " [-0.45889088 -1.4770781 ]]\n",
      "...\n",
      "\n",
      "Class labels of those first few points: [0 1 0 1 0] ...\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAFgCAYAAABKY1XKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvXtwJNd93/s53TM9b8wsFgtwueRyQRK0RIuiI4uxBcsivVRKguJrVmiXYusm2dhVccU2K6xbN3sjJ7Z8i76JmGyuFVZWVy5XEmdvYodWxYzplAmnLkVzJQqUIjkxJWtJcZZc7IO7Cywe8370dPe5f5zueWGAAbB4DXA+VVuz09M9aAx6zq9/r+9PSCnRaDQajWYtjN0+AY1Go9HsfbSx0Gg0Gk1ftLHQaDQaTV+0sdBoNBpNX7Sx0Gg0Gk1ftLHQaDQaTV+0sdBoNBpNX7Sx0Gg0Gk1ftLHQaDQaTV9Cu30CG+WTn/yk/NM//dPdPg2NRqPZKGK3T+B2GDjPYmFhYbdPQaPRaA4cA2csNBqNRrPzaGOh0Wg0mr5oY6HRaDSavmhjodFoNJq+aGOh0Wg0mr5oY6HRaDSavmhjodFoNJq+aGOh0Wg0mr5oY6HRaDSavmhjodFoNJq+aGOh0Wg0mr4MnJCgZn+SzU4zM3OGXO4Smcw4k5OnmZiY2u3T0mg0Ptqz0Ow62ew009NPUSzeIBodpli8wfT0U2Sz07t9ahqNxkd7FppdZ2bmDIZhYVkJACwrgW2r7b28C+2FaDQ7j/YsNLtOLneJcDjesS0cjpPLza7YV3shGs3uoI2FZtfJZMZpNCod2xqNCpnMiRX7tnshQggsK4FhWMzMnNmhs9VoDibaWGh2ncnJ03iejW2XkVJi22U8z2Zy8vSKfbu9kFotT7F4jcuXv8a5cye1h6HRbBPaWGh2jGx2mnPnTvLcc+MdC/vExBRTU2dJpY5Sqy2TSh1laupszzxEuxdSq+UpFK7gODahUFSHpDSabURIKXf7HDbEhz/8Yfntb397t09Ds0GCXINhWITDcRqNCp5nr2oU1vM+xeI1HMdGCEinjxOJpLHtMqnUUU6demVbf5fuBDugk+6afgz0DG5tLDQ7wrlzJykWbzQrnoB1LexrLcyXL3+NUChKKnUHkUgaACkltdoyTz/97rb8Hr2MXq2WAyTR6KHbMoSafc9AGwsdhtLsCBupeApYrfIJ4NSpV7jnnh9naOiupqGA1RPjW0WvBHu9XqBeL+qku2Zfo42FZkfYSMVTQL/Kp40kxreKXkbP8xp4ntOxrZ8h1GgGDW0sNDvCZhb2ft7IRhLjW0Uvo2cYYQyjs791LUO4WqJfo9nL6A5uzbbSnnOwrCGEgFptmUzmRN8kcCYzviLP0b0IT0xM7WheYHLyNNPTT2HbNPMTkcgQoAxge86ilyFsz3l0htZ0fkOzt9HGQrNtdC+MwSL6qU99cV0LY6+FebvDTGsRGD7bLuO6SxhGhNHRB/nEJ34LCKqhZtc0hL2kTYrFEi+88Bmi0YyupNLsWbSx0KyLzegxbVTzqRu1z9l1LcLbTbvhS6WOdRiu4HzWc1653CWi0eHm81otT7U6j5Qemcy92tPQ7Fm0sdD0ZbOhk+6FEVYmfvsZoe0MM23EAN6u4QvoDq2Vy3NICaFQrJnE38z7ajTbjU5wa/qyWT2mfhVQuykKuNGfvZnS3150J/odpwZIksmx23pfjWa70cZC05fNLpT9KqB2UxSw1892XZsXXvhMzyql9Zb+9qt06q7gsqwEyeTYjvaKaDSbQRsLTV820yMB/Utbt+puHTZejtr9s+v1PKXSHLZdRgiTa9e+yfPP/xRf+tLDZLPT6yr9Xa+3MjExxalTr/D00+/y5JO/j2laO9orotFsBi33oenLRnSd1psHyGaneeGFz2DbJUKhGJaVwraLOE4V04wwPHwf9XphXcn0zehOtcuP1Gp5crlLSOkihAEIDCMMSIQwSaXuYGrqLKA8klu3LuC6dQzDYnT0B5vnd+7cSRYWsth2HtetY5oRLCvNyMjEOiVNdjeJr9l2BlruQxuLfc5WTZVbz4K23kU72K/RsKlW5/E811+oTYQQeJ6HYRik0/dgGOENLfwB/XSnstlpXnzxF6hWl/A8e8XrhmFhGCZSumQy9zbfa63f8Y/+6O9Sr+cAgRAGUnqAJBI5xOnTNzf8mWv2HQNtLPZENZQQwgS+DbwnpfzJ3T6f/cJWNoCtpyppvRVDwX6p1CHC4Si53KXma0KEME0JCMrleQ4ffqBvddB6qq56I5HS6bFd+BIeDQCKxWtUq8t9f0fPs5ESTNP0fxcT13XwvPqKn7DTo2H1KFrN7bJXchZPA2/u9knsN3Y6gZzLXcJ1Gywuvs38/HdZXHwb122sWLTb8wXRaBrDMDHNGIZhAh5CGAhh4LpqkfW8Blevvr5qPqI7p1Kv51lYeItS6caq+YuZmTNEo4cwjBChUAzTjNC68ZPNf4YRxnFsbDtPNjvdce61Wp7FxbdZXn6Xa9debx7reR5Sgucpz0K9d4v15ja2ShZEj6LVbAW7biyEEHcBfx34N7t9LvuNrUwgrwfLGqJQuILr2v5dtU2hcAXLSnXs1724q31ruG4Dz3PxPAcpPUwzQr2eJ5+/jBDGqgtde/K5VsuRy83ieQ2SyWOrLozBZ2OaEaT0fINl0R4pMM2Ib7ggHh/h5Zc/S62WY37+u9y6dcE3jrZ/jIFtlwFw3TqOUwUgmRzjyJEHO372eiqxzp9/ZssWeD2KVrMV7IUw1L8C/g8gtdoOQohfBH4R4Pjx4zt0WoPPerSV+hGEL1ZL6rYjBKg7ctH2T/rbW7TLeHhew19wJcHl6HkNhPAYGjpGoXAdEKRSd67atNbe6X316ut+iOtOolFVjmrb8PLLn+0Iw1jWEI1GhURijELhip83AdO0mh6N69qYZoR0+jieJ1lYuEAyeQe2XW72R7iu6/9WYTzPBTwMI4IQIKWL6zZWVDZ1h82CSiwpPTzPoVB4j8uXzxONDpPJ3ANsvgmw18+D7b1p0CGv/cmuehZCiJ8E5qWUf77WflLK35FSflhK+eEjR47s0NkNPrcr4R2EL5aWslSrS9h2hXo9x8JCtuddbr1eIJ2+B9MMI6WDaYZJp++hXi+ueG8pJYuL32d5+R2kdIlEhrGsGIZhYBgWoVDUTxB7DA0dby780HuhC8pRU6k7GBl5X8f+rttgYeFCx116uTxPrbaMYYRIpe5GCBMpG0gp/eR2xPc0VEipVLqBYViYZtRXmG0VhphmpGkohDAJhSyCEFY4HGdm5kxHKCnwrOr1PPPz32Np6SKe10BKF9uuIKWLlB7V6gLF4o01f+/1sNnS582gQ177l90OQ/0Y8FNCiFngeeCkEOI/7u4p7R9uV8L75Zc/S7F4k0LhPX8xBBDYdr4jjBHE1ovFmxSLN0gkxhgdfYjDhx/AMMIdi1JQhVQoXMUwrOZ2286TTKrjRkc/QCIxytNPv8tdd30E0wx3nNdaC12vhTFY6NvDMNFohmRyjFTqqF++ayGEiWGEiMcPI0TLQyoWr+N5NpGICrMpI9ZCVRS2jMfhww8wOvoQqdRRCoWrKxbOEyceo1bLsbT0Lq5b6/oNvI73L5WuU6/n+/7e7Z9vd55jJ+d+6JDX/mVXw1BSyl8FfhVACPEY8A+llH9rN89pv7FZbaVsdpqFhQuoQjW1EHpeA8MI4br15l1ucCepQkngOBWWly8C+D0LEApZZLPTTExMMTNzhnq9AAikbDR/npQepdIckUi6Y1HcqPJssH+xWMK2czhOHc9zSCTGOvYLh+PUass8/vizfsXYIWxbeTnV6iKx2GFsu+SHmwRHjjzI4uJFVA+GSTP6BB0VVaIt5lYstowUtEJJs7OvkkiMUqnMs57K9aWli4RCcSKRoabCbfA3ag/3nDjxGG+8cW5F9dvU1FmmpnZGkHGnQ16anWMv5Cw0e5DgDlFKz+8ZCAyGg2Ulmwv6zMwZXNemXJ6nu2dHHWtSrS41S3ZVUrjedves8hogaTQqHXe9q0mCr7XQTUxMcf36KV577fN4noNpRhHCoFpdxLKSzfBU+/kHC3ooFMVxanieQ6l0s5ngDnDdGkKE/QXeALzm72maUd+YWkgpm0ZtaOiejvNrLZzSLxM2kdJr5knWpvX59iqLfu21zxONjpBKHQI68xynTr1y28ahl3GanX21IzexFXkyzd5kt8NQTaSUr+oei71DLneJVOooQRxeoRZ1y0o3F/Rc7hLVag4wkNKlvZpICAPDCFOrtcJWmcx406gIITruxKX0mqEyoBn7tqwkjmNTqy1SqSwCa5eVvvnmHwKGX8kkiMUOA5Ji8XpHGObEice4evV1lpffYXHxbUwz7PdWqN9TVWjVCYdjOI6NECE/ce0QCkUwjJBvaJJkMieIx0dIJse4desC+fwspmnhOJ1GIFg4M5lxDCPUrMRau19LcOTI+4lGDzXDOb3CPZ7nYNu55lG1Wp5i8RqXL3/ttifydeciFhayfPWrv8nSUnZFiG2nR91qdgbtWWh6EtwhptPHKZXmcJyq31ltMjIy0by7t6whcrnL/lHdnoVs9ksEd9Sf+tQXuXz5q0jp+UYjWCQNEokjzY7rc+dO+nfpDoXCVdTib7K0dJEXX/x5QBCNZpoL1Ysv/gKJxCjl8hzl8px/127RaNSw7RJCGDhOjULhPUZHH2yGbJTnIP39VibiAWq1HPH4YeLxUWq1BdLpccLhOJXKPKXSPI5TJ5ebJRJJUyi8h/I4QnieTal0jUplnnT67mY3erBwvvjiz/sd5Cs/u3ZCoSi1Wp5y+SZLSxc5d+4kt25dIJU61rGfaUabxqlWy/v5FXX87c7J6G5GtO08YFCr5UkkxjpCbDsV8tLsLNpYaHoyOXmaF1/8Ber1gn+3bWAYBrHYoeY+2ew05fI8qy10Qohmv0SlMo9tl3jppV8mFIrgunYzaW6aURKJEYaHJ5rHBrHvpaVs82dLqTquK5VFpHQolzslNMrluea5SOngOG7bcyUnUqnMMTs7x+zsn/kek9GRO+nxWwBQKLyHECaOU2N5+SKe5/q/nyQUihGJpCmXb7Qd18pjeJ7N8vIlxsY+wOOP/1Zz4Xziid/l5Zc/y9LS2yiJELMpEdL+mUYimRULf72exzTDHbmYWCxDpbKAbZcpl28ipSpnTqXuuO05Gd25CNet+/0xLc8puCHY6VG3mp1BG4sDzto18UGeQpWwCmFgWcnmXarKAWTwPJtqdWnFe6tYfA0pXQqFEoYR9lVd1cIfiQzjumUcp+4vxqFmIjzwbNSiFGq+n3rsJdHROt/VnqswWffzzm2rvW/QZKc+D8c/Xr2/41RwnErPIwOEECwtXeSll36543MOPmtVJfbzVKvLTZkRgEjkMLXaUrOAwPNUR7xlJahUFgiHk83Ev2lafPSjv8rs7Kt+UjxKKnUHAIuLb+M4tWZRwkYX8+5chGlGcJw6oVCrO13nJvY3WkhwH9OvOWotUbyZmTMUizeQ0mF5+ZK/sBqEQhGOHHkQ2y6Ty836jXNX/bt+tytxLZrehUrohv0ejCChK5v7qt4FiMWGeeSRX+bChT9kYeECnuf5uQ9VleW6DdYK2ex1DCOMYYSJRIZ45JFfaiaII5EhcrkrvsfVwDDChEIRIpE0+fxsr3ciEklx9OiHeoZ7AnFFKR3y+Suo9GSniu5GDEb3tVIqzVOp3CSZvIN4fHRdSr+awRYS1MZin7IeBdi11FpzuUsIYVIoXG2TtACQDA/fj2UNsbBwARXGUSqxgC/X0cA0owwN3enPiAhyAYJwOAZAo1El6PYWIpAEb92tG0ao2SUd6CtFo4dWhJ4GjUBeJPAeQqE4qdRRisUbeF6DdFpVT5XLczhOrcsAdyN47LH/k0cf/dyKV1oJ6ZtNVV8pXQzDxDAshofv55d+6Y0NnXu38nCrGkrnJtaJNhY7iTYW62M9st3PPTdONDrcVZEkqdWWyWROcO3aN/2FzW4ri1WVN6nUMUzTavZitEtyBwJ6SuJbdoRVAh2lVvinXbwvKKNlxWut999vtOdjgnyF8rSECOG61bWPNkJ87GO/vqrB+PKXn8RxGigZklCzAktKl5/92T++Tbl6LeexQQbaWOyZ0lnN1nLr1gUKhWtN9dd6Pb+iOWotGYig6zcIXbRUWEN+ZZTNxz/+LCMjSiTPcWq+B2I0jcRKQ4H/Pt15Atn12LZdSvXg7UdDAeD5womtZLwyuM4auRmFquQy+MY3vtCzlHhiYorh4Qf8bnQVJlRFBWLTXdXnzz/Dl7/8JFeufI1yeYGlpd7SL5r9hzYW+5Bsdpp6Pe+rv4Zw3Qb5/BUqlfmOBORaMhATE1OMjDzYHEhkmlFCoRiqzyLZDGc9+OBP++GNMKYZ9fMRKq8QJILXZg3Pdu189T5HfX6tHpcee/jhrFot11OP6fz5Z/xmSa9tfxvXtYlE0mvKvvcim532mx09DMPC8xzK5Xlc19ZyHgcAbSz2IWpWw0hT+VQlmaFSWehojuqnHfXxjz9LKnUHmcy9HDnyINHoYV/KO8LMzBmy2WlmZ18lmbzDr4pxCYUiDA3d2cxNrEWQp+jJaobhgBmM7gquXvuAKu3N5d5laSmLlA6GYfHVr/4zv5y4G49KZX5N2fdeqAFPKu8hBH6eyqBazWk5jwOALp3dh+Ryl0gmRwmHo5TLc81Qk+s6zTvAdnnvtaQzAunvGzf+B/V6ESEEtl1maSnbbCoD4RsJpQ4rpapaUlPiVo4sDVgZotJsFte1CYViTS8yFEr0nNAXoLrlV5d970Uud8nvkWk0PR7VdFnTJbMHAG0s9iFBTXyggxQ0dJmmtaKTN5udbjaGSQkjIw/w+OPPds2KgC9/+clmgtTznKZ0dtCFbdvFrg7oYJ6FZuvoLgBoRzarnTwP6vXlPu/V+T7hcJz5+QucO3dy1cR1JjOO4zjNuelCGL6nEdJyHgcAHYbah7TnInp18gbJzfPnn+H5559gfv47OE4d13W4desCL774Cx0hiT/5k1/2BfZsHKfiJ1+lLxQIvRvbJEr2YpMGY7XDDrT9WTsGpyb02esUJYR8fpZaTcmfl0rz2HZ+zTkUk5OnCYctYrFRX0fLwTAMPvrRX9XVUAcAbSwGkelpOHkSxsfV43RnrLk9F9Fo1DBNi3T6OJGI8jSCu8hAmbWFi5SCer3QDFedP//MiqYwVaWjyjz7l7O2T87bIN2HHGhDsT7U32Z9iR2VoL7pj6NdIB4f8UNSBYrF98jnr/DCC59pGozguhoZmSAUSmBZCSKRNLOzrzb32aq54Zq9hw5DDQAdde3uEJP/ZY6JhUMwPAw3bsBTT8HZszA1tWL/SCSJZSWbhgLwG/TqbYaiFd6QsoHres2E5de//i+24De4jay0NhDbSjD4KRJJE4+P+jPPVce3EEqepT1sGXgQwQyQcDje9EKuXz/Vc57GZsULNXsL7VnsMv3uxFaMqbx5kekfWSJ7wlGxpUQCLAvOnOm5fyiUpFS4QfnSXyK/+x3si2/iFXOYZgTTjPo/ZaWGkpSSbHaaRqO8Ex+DZpcQwsC2SySTYzQaFUqlOYJGQfAIhaIrejJWm4b3jW98QU/J28doz2IX6TXApvtOrFsa2qo62BGDmfvmmLiVJnskz1c+do2FxFuIfxpDCINIZLg5ACflRaEkscMNjFiITAkmX5fM/I0xlmJLqxqDfP4yf/Inv7yPO6c1EMwcsZASPydVRYiw37wnSSTGVjRzrjYNr14vksncu2K7LqvdH2hjsYusMAQ9ShhXfDEjEcINm/lUhd/+2AXmU1XlF0gwpfQVUOuEw1FVDTU3R7IRpiYkT7/ykHqPchm+Di9+ZLX4tgpL9RawU6jZDI22KXoHqgFi36Aq264BJpnMPX7jXs1XrD1GNJrGtssdpbFKF+wtpHQwzQiJxBiGESISSdFoVPSUvH2KDkPtIrncJcLheMe2vpIcqRQls45tuizGfEMhALNVxiql12rGqtdphCFTaUlJE4/DnJr90D42VNErq9y5j+FCvOD6CiAtxdgOZNs/zZ7Gtkt+JdR7pNPjfnm06sJXye8clcoizz03zpe+9DCFwlVfL0zgujb5/GVqtWV+9Ef/Nz0lbx+jPYtdZD3ziicnTzM9/RS2DeFKg0ZunloS4mUoJ/ALjdQCrzp3w344oaZmQcfUF//EQpLf/tgFFpI1hASBIGbeyaFD97K09A699ZmUh6EWD6/5siegHPWIl6ASd5HdpbO9ZDp0onrPU60uEYsdIh4foVyeY3n5HcLhBIYRwnVtf5zqW3ieTTw+im0Xcd06hhEmmRzj0Uc/x513PqKn5O1TtLHYBFulujk5eZrp//wL2MtXCFcbNGJhvENDTH7it5r7tHdR5669TqZuUYtCXIapyTouHkEjhZQSwxCYZoRwOKbUY++4nxPnL/Ot989TtdzmQu6EJI3CdZLJO9ruJLtRq3zztZbwLNKESpxgTEKneGwvtMEYACT5/JW2mSECz3P8GRuOLxvjACa2XeTw4QfUUb5SMaxUBAgKOLRC7eCjjcUGWU9Ser1MXAReksw8BLmkaCafJ34ImGjbL/gCjo/D8DDnPpKlGLVJ1sPkzDqegODLLaVHLDbME0/8bmtuBQ9TX/5LhARDGBAO43p1wKNUuo5a8VX1y0rWiCPpIOY+Q7TNLhH+c+U5lMtzRKNpTFONxG1v/FstL7GV3xXN7qPnWWyQ9cyJWDcnT6o+iUTrvbKpOWYeLpEbz6y8E/P3z55wmH7oCoYn8HApxF1cU2IYYUZG3s/HP/5sx5fxuefGKRTewzDCCEFfzaZVWetS6edZtO+jGQCEH9IMGjAFphkGDH+4VYSRkfetOSFvte+KaVrE44cPorcx0N8A7VlskNXKBjdVHnjpkmqs88keyTP90DxGwyMavZfijSzT//ZJ+NoQE5EfhMceg3PnmJi1QN7NzL03yMVd7s48xOT/8mzPL9z5889QKFzD8xy/einECg+i17yhNXWI1mC1wwb6a3IQkf4NRfsfzmz24CSTY80hWast9r2+K67bYHn5IocO3a+9jQFDG4sNsp6k9LoZH+/wLGbum8NwwDJjUChgvTePHYaZv1Jh4is34Nw5OH4cXnuNibccJkIh+Jt/E/7lf+z59ufPP8P58890yFx3D9QxXDCloGH6K3wzt7ABQ9GrgKqnAdIMHuoPqW4ylAS9ZY2SyZzo60n3+q6USjf6lotr9iY66rxB1hoYtGFOnwbbVn0PUpKL1Qg3JIyNqdJWwyAsTXJxWxmUfB5efRUMA2Ix9fgHfwDPPNPz7b/+9X+x9jwECZ4BDUOu2L4hVvMktPDsvsE0Q4yOPsThww8QCkW4dq3/4KTVviup1NGO/XTj3mCgPYsN0lGddLvlgVNTStPpzBmYnSXTSFA8lsRKp+HKFQiFaJhuq0cil1OPoVDr0XHgC1+ARx5R73PpEtkPDTHzY6xfqmM9C/pgpbY0W4zjVLl16wKuayOlh2n2Tlp3Vwo+/PApZmdfbX5XQiELx+nMl+nGvcFAJ7j3EO3VI+HZazQ8Gy8EU989zsStNPz5n6ueiljbFDopoV6HEyfAssje02D6A1cwXMlips9Y040ko/tdJtqDOAC0qhgMI0Qmc4JIJN0s8Ah6ggzDIhyO90x+d1zjq+yzjxnob4k2FnuM5p3Z9QtkLueZ/P4IE4VRqFTg7bfJ3ieZ+aggl5Zk8oLJGcHEOwLuvx8SCX77YxdYTNSQAlxDrn15bsQA7ESVk85zDAACIUxfL0zNY08kRpHSI5M5sa5KwZb3obyNEyce872PfV8dNdBXtTYWe5np6WaIihMnyN4vmY69iuFBuAGNsMo5PPymxewDFrdSNcoRB+GB6alOa8+k9yW63j/7TpXE6gqqPY8QJoZhtvVigGlGkNLhyJEHqdcLRKPDCNH6wwUNe08//W7P9zxgnsZAX9E6Z7GXmZpqzqgAmDl3EmN2BGtuGaSL5ZqURxK89qMl0iWwTdWhLQ2/qTs4sLt7eiP3B+vpvL5dj0B3fQ8Eat52kG8I5p+4qGZQVf20sJDFtvO4bh3TjGBZaUZGJlZ9z/WIaWr2BtpY7BLrlQxp369YvMnQoWNw9J7m69WF7+MhsBrKyzBc8ELgmp3vE2oInPAOeJHrNUzr8Vg0e4peUQjPc7CsNEtLb5PLRbDtvO+BWDhOHce5yYkTf2/V99zSviXNtqKNxS6wpgzCRTqqmqY/PI+RyhCNDlMqzZPPX6bRqDZF3FzXRgiD+SMenuu7Ex6dRdHSz19sdmHeTCltv54+bSQGjvYeHcOwMM0QrtvAtvOEQjHULHaBlC6uW/X3FHzjG1/gzjsf6XkztKV9S5ptRfdZ7AKrTRqb+a+fVSNSb9yA4WFmxi5izC9gVZWIWyp1J1JKSqXrOI7dvNOT0sOVDtIX+FvxV/VUaGpH2QpjoENQexbPc2g0qk0p82g0g22X6CU5XKvl+PKXn+T8+VY/UCAweOvWBfL5WYrFOaSUlMtz5POzzM9/T8/w3mNoz2IXWNX1vnYBrBPNju5c0iFaN1WDXjpNNJpGSrWCem1Cbn0X1W5l2EFgkM71QNIpGVMq3Vh7b8/jtdc+z513PgLQ9KxTqWOYZphKZYFGQ3nL0egIyeSolgLZY2hjsQt0uN75PMzN0XCrZHIeNFpS4ZlKhGLUxqorw1Cr5QEn0HVr0c8QDNrCO2jnq6GfK2kYJp7nNOdxtye11ejWJOXyTdLpEzrZvUfRYahdoCmDsDSHvHKZslEhn3SZH/E49+Pvko1fV/u9M4aHhx0LIaWkWLy+umHYLzkAbSj2JUHn9+XLX+Pq1dc75qfU63kKhWvUajmKxWv+TZFCJ7v3Dtqz2AF6VT5NTZ1l5ouf4VbCox72iJYlyTIUE5LpD92At+JQq2PVJQuHqohr/wPPECQq/oS8/YI2DgeCQKNMCBOQ5HKzmKaF5zXwPBchDIQwcRybQuEKcJxoNK2T3XsI7VlsM0HlU7F4o6vyCU69mOFILkw6J0lVhEp2N8Bw4OXxd5j+0E0cQ3Ikb5GuRjA8iekK5dyLAAAgAElEQVRAqI+Kx8CgDcWBQ/VqNJDSxXFqeJ4LSKR0saw0QqgeoXL5pp7hvcfQxmKbWbXyaeYMjI+Ti9UJBx65X90UbsDSYTBsB2nCUqpBLm7jSSj787cHHm0oDiRCBNUWAdL3KsJI2SCdPo5pWjQaNVKpo/u1k3sg0WGobWbNpqPTXyTzB69STIFlqy9QLQL5IdVYt3gYwMVwVbOd0nuCamzFjxk8Bq06S7MlqA5wda1bVrK5TQgT160TiaQRIrS5yZOabUV7FttMJjNOo1Hp2NaMw05NMfnePXimwA5DNQLLGWUo2m++PNPvk/AXV6kXWc3A0rpLUL0aleZMb8OwdOhpD7OrxkIIcbcQ4s+EEG8KIb4nhHh6N89nO+g5AKaYY/KPF2F8nIn5Iaa+PkSqYlAaUgrkgSfRjmv4ooCAE97532Nb2KoKLtn2T7PnUUlugee5GIYFGEjp0miUWFx8i1JpjuvXv7Xbp6npYrc9Cwf436WU7wd+FPgVIcSDu3xOW8rExBRTU2dJpY5Sqy2Tsi2mXpJMXLDV/G3bZuJajFN/dpxkGYQHhgfmaovffvMqbneBX9kwrNnTCEIhi1hsmFAoArgYhtF8zTRjuG6Dr371Nzs6vjW7z67mLKSUN4Ab/v+LQog3gWPAhd08r61mYmKqlaQ7eRIW7GaXdvPx6FEyGUGhcQknRKe2kgQRSDvtN2MBrd9zvbM3+okQ6nzIniWTuZdC4Qqu2yCTuYdIJM3Nm3/RfN0wBBDCdR2+8Y0v8Oijn9u9k9V0sGcS3EKIE8BfAb7Z47VfBH4R4Pjx4zt6XrdNMJPiwgWo18keyTHzUZPcYUGmGmPynTEm5BBcuMCJVI3ZH6K10LUtihJVUuvtlxBUL9bbcKiNwUAihIkQglAoiuPUKZXmiETSbT0YomPfer3AuXMnD8JQpIFgTww/EkIkgfPAP5VSvrDWvgM1/Gh6WgkD2jbMzZG912P6E6q6KexAIxbGCwsevmAxO1Lmyt2ymZdoJ1JVa6NtKRnyfb9Q9lOsXe97aPY8mcx95POX/LneEQxDfQEcpw5Ihocf2E9DkQb6qtztnAVCiDDwh8Dv9TMUA8eZM8pQ3LgBjQYzP6oMhdVQYSWr0sB1bF77oRLFpFSGoAf1GDjCv9IG+nJbJxud4qcZWHK5d4nFRvxRrRLPk7iug5Qu0eih3v1Jml1ht6uhBPBvgTellL+1m+eyaaanVR5ifFw9TrdJKl+4oBRjXZfs/XD1Llg+pPon6hG1SzWqKp1yadZc/JwohOoqAa5ZB7vvMGvWhSSVuoNHH/0clhXH8xpYVpxIJE06fQ+1Wp7FxbeZn/8uxeI15uc3ns4M5NCfe25cy57fBrsahhJCfBT4GvBdWprH/1hK+dJqx+ypMFQQZrIsiMehUlGexNmzahzqoUNQLJIdd5n+FBRS/gQ7P8ySLCmdJylY912y6a6cgncg2anZ4JptxzSj/NqvVTu2nTt3koWFLNXqPCAQwvBLbQ0+/ekX1h2K2mMzvgf6itwTOYuNsKeMxcmTKsSUaFP2K5fJPmgx81OHyX3nPJlFj1waikM9FnmdqN082ljsK8bGPoiUYNsFMplxTpx4jNde+zye52EYIaT0AI9EYpTh4Yl1d3efO3dyxSQ+2y7vVof4QF+Re6YaaiC5dEn1SrSRvafB9ImLGMX7iYo4C4dKlNL+i71KRLXBuD1WS4brz3RgMIwwt25dAARDQ8cpFm/wxhvnCIWiuK6D59mYZoRkcgzLGtqQZLme8b117HqCe6AZH1ehpzZm7r2BYfrCgWN3YMfoWMwENLuNrQaazdBtCPo91+waphnps4fAMEIIEUIIg0plvpnMBhgauovR0Yc4fPgBIpGNS5avKbej2RDaWNwOp0+rHEW5HOgqk4vbhEeOqtfTadyIiWiT7mjaDQGigV7YtgrR9k+zJwgMwFrEYodwnBquW8d1bRxH5S7C4ThSSvL5Webm3mBh4fuUy3Mb1o3qKbejtac2hTYWt8PUlEpmHz0Ky8uqC/vOB2nEW51zphVHhM2W1lPbglaPo6t2NsJWGwOtKbXNyI6JeL2oVpe79nep1fJUKvPNedwqHFWnUlng4YdPbSgxvUJuR8uebxqd4F6DXhPu+l1k2f/0DNNvfB6j7hA2o5QORah4eZJmBje3SDmuPm+zAW5gU3Teoj9b/fnoPMe2YxgRlAGw132MECFM0wI8TDOBlEqR1jQjWFaakZH1J7f3IAN9hWnPYhVWm3C3Zo329DQTv3aOqW+OkGpEqYk6I1fLfOyOv8nwsYeoJ0xCDRhegtEF9F3tetkJQ7HWds2m8Lw6rYr4tRHCwDSjCCGQ0sMwIjhOsW3WhU21Or+pPgvN1qCNxSqsOeFuNc6cUT0X0SgqKSGoRiUXvv9fyOUuIV2XVBEidbW76dAMrXRLkms0+wHPW88MYIGKcHgYRpjDhx/AMARSgmEoPSnDMJEyMEAK3Wy3s+jS2VXYTMldtv49Xv5knoWhOoanNJ0KKUDAULmG8CCfBpFXBiNdhKWMmmEBvrJsEJIaaId1i9Fhun2Oculct4EQLuXyfHO753kIYfh9FrJZXdXebNc5217nI7YL7VmswkZL7rLZaaZ/vMBSwkZ4yvkuJ9XAIs+AnH0TTyj9p6VhuDUCOX/YUciGWA3SObShWI2tDNmt9vnqz32XEQghME11D5tMjmGaYaR0MM0wyeQYR46ocTeb8vw1t4U2Fquw0ZK7mf/6WYyGRAqJhJ7qsV4g9YGadueFlNaTG4ZKrL8+lIbtMxj6c99VhDAIh6MIEaJazWEYEUzTIpU6xpEjHyCVOoZpWs3vXy53iXA43vEeutlue9HGYhXWXXI3PQ0PP0zu6ncIV2xMl1XVY5tImr0XnqGeSwP911gvW2kwtCe3J1BhJmU0XLfG6OiDa37/dLPdzqNzFmvQMeGuF4GQ4M2bZD4kKCYliRLYw6sfAqjFyTcSCN9QaDQHHNdt+AlxSbW6CLBqmezk5Gmmp5/CtukQCDxx4jE9MGmb0MvUeuklRR5UPzkOk/89hGcKDA+MPjIewvOVZvGT2vrOdvPoxrp9g2rgk0Sjh1hcvMjzz/8UX/rSwz2rnHp5/g8/fIo33ji3sXJ3zbrRTXnrYTUp8nIZjh2DbBZsm+x9kplHbG4dgUrcr2zqYY4TRZX8RviS4wdh+t1OoT/HgcY0LaT0ECIESIQwSaXuWFfX9R5TmO3FQF+d2rNYD4EHkUio8qVEQj2v15XhGBsDKZl4R3Dq90L8w38JH/hOqyRWZbxVL0V6GcKuqnwyXPWS6bDe3iVNPwbr3kfTgcB1bTxPTcqTUuI4VfL5q7zwwmf6egg66b296JxFD1bIfNS/x0T8rs6d4nGIRJSHYVlw991qtoXrkv1wmmv35EkVoBprzbFIFZQBsWpgR/x1zVBhKNMDVydbNQealqVvaUoJhAhh26W+fRSZzPgKz0InvbcObSx8AgNx69YF6vU80egIyeQoCwtZvvz4AhF7kSOlOJPvjDFxK608igcfVMqzZ86QtS8w85MWubEUtfIiVt0gUfBIVNQI1VxahZ6SRbg1qpLawlX5C8/UhkKj6Y1ESodQKNbso1jNWKyW9NYKs1uDNhZ0doOqfgqPanUeKV1qtUWkKbBNl6JVZ/qhK/CtUSaWLWUopqbI3o9//CGi4TiF6hy2PzM7UodSkmZ5bCXZSmhLs1OyXLMF6M9x3+F5DRoNj0LhGrVaS6W2l9Dn1NRZf9ssmcyJpqHQFVK3j05w05kYm5//LkKEkNLF8xx/MIuJdG1GCzFst0pKJjn1id9XEuWsTKwtXv0LHFxCLhxehPlRQPqhJlPP0N42tKE4AAjGxh7i/e//ad5449yas7Wz2Wm+8pXPcuvWBQzDIpk8immG9QzuTaIT3HQmxkwz4ldjGEjpNnVpzHAMHniA8Ps+SG78UNNQMD1N7s3XCX//HXj7bcjnScRHEYDjew5BqWyipKqfNFuMDuEdIEwWFy/y2mufp9GwV5X7CKIFi4sXARMpXYrFq/4NoJYF2QzaWNDZDZpMjgEenucihInnuYDnb+9KmPkltZmCQSNiQKMBV64QDSWI1QwsG2oxOLwEsSoYHsRLu/Ir7m90n8UBwsFxqjiOjW3nOl5pr3wKtKOkdDAME8NQWjvl8pyukNok2ljQqQNlWUMkEqMYhkE4nMQwDBKJUSxraKU+1Gc/CzdvMvlKHc+1sQ0XaQjsW9cJD4/y5CsZnv5iiL//2/DEi5AqAobyNDTbgDYYBwRVi27bJWq1fHNr+41cEC0IIgUQSInUdYXUJtEJbvBjl63E2PDwBD/5k7/DxMRUWxKtlTCbmJhSXsWFC2CaTFyx4KU6Mx9pkMtA5pZg8n+OMpFV84QxDCYuekxcVE9feAK++zA6dLIdaDnzA0U+fxkpjzdzEcGNXFBGm0yOkc9fwfNANfmFdIXUJtEJ7s1y8iR885vgeeqf3TU60jD8Fm5UZ/ck5A6pHotyEjwB1Th6YdsO9Gd6gDAwDJPjxz/aUeXUXuHoeQ2KxRt4ns3IyIN8/OPP7lY11EBfmdqz2CyXLsHRo3D1qspVdKNuZVRZ7adUt3a0CgsjSmk2k1fS5HYIHQzUaDaBEAZChAF3RTlsd7Tgrrt+RJfM3ibaWGyS7IeGmBm7SC7hkZmXTM7QDDO1MzOpDIXl2xPpq80WUuCEGPB7DY1mOxGsnohSnd1CgBC9m/X6qkZrNoQ2Fpsgm51m+sPzGDfrREsexZTyHnhppcGYPwKOpXorTL9jWxq+oRisCODgoPMW+wJVxRQGwHXt5mhVEH7/kzImyeTRdVc39Wrk0wZlfWhjsQlmZs5giBCWowr8rYbERnkR7cYiez/YUZWfMDylLusZLe/CcNW0PI1G044gFjtMJnNPc0ugHlupLLK0dNEftWqRSIxhGCFSqaM9DQHQ3BaJDFEqzRGNHtJzuzeBXqrWQzC74tIlGB8n99cuEF0sQyikEtnVKuGGJHeo87CZSYiXVUJbCjW7Qkh1bxRy9NAjjWYlgoce+gzXrr1OLneZanUJJclscOedP9zUf+ru3D5x4rHm9sAQvPjizwOCaDRDNDrMwsJbeF4Dy0o1G/lsmzX1pjQt9HLVj2CWxY0bMDwMN26QeXeJRr2kJMprNQiFaIQhs9x56PwR1ZQnhfIoglCU4UK84ndzr9VzoQf73D76MxwgDIaGjlEsXmdo6DjV6gKtL4jke9/7A65f/1bPcasXLvwhxeJNcrl3WVrK4nkO9XqRer3Q7PKW0kEIg3J5rvkTdYPe+tGeRS/aPYlcDpJJOOS7DY7D5Ndcpj8BdlgSbkgawsMzYHKm9RYdISjXFw8UEKkqI1HMoPSiXHChOWY1WgPHACeCjrvfLt0GQucy9jgehcINisUbSLlSF8fzJN/4xhf4R/9oucMTyGanWVi4gBAmQpi4rk2hcAXXdTGM1v2waUZwHBvXrTe36Qa99aM9i266PYlSCebnIe93is7NMTEbZupPDVIlQS0GqZJgahomLplgKpXAIAQlUF6F5wsIlofUTIt4EUIuYIDlQDIPJ2Zh4vu+odDcHqt5EtrD2OO4PQ1F4B7W6yv1cgJpD1UhJZrSHkLIZoIcAikf1xcKlSsVGTRroj2Lbtqn4gHEYircNDcH6bT6vxBM3Egw8doYvPceVP1ObVoXee6QCjV5hi9RHlQBSshn/DkWbaa6moDZFNp8azSr4hEKxVZszeUukUodpVC4iufhi3+q1yIRJdMTDscRIkQsNkIiMUqtttypyKDpizYW3Vy6pDyKgLExuHxZ5SakVElt21bbAVxXGRfXbT1H5S+KKTURz3TBkCokJaQqm5UmzbBIIFkeJL8124wOR+1pAqXnXjhOjfPnn+HRRz/X3BZIe6TTxymV5nDdOoYR4vDhH+Dxx5/tkOv5xCd+SxuHTXJw5T66KpyCQUacPKlCUInWaEbm5lQ46tAhZTCuX2/JfAihQk+B0fC9jKBzOz/kCwcGMtoeuLoZb2+g/wYDhoFhhDAMg09/+oWe0h6rzbbYIwz0FXcwjUWQl7AsNUu7UlHewtmz8K1vwec/D44D0ShkMmq/s2fVsU89pfIXuVzLmzh0SO1v2ypM5ZO9H154UnkXIUfNsyik1xh+FPwpBvqSGjD0Zz0QCBEiFLJwXQfPaxAKRbjrro80w0jtgp+RSAopwbYLe63xbqCvtoNpLHp5D+WyMgqlklr0czm18Jsm/Oqvwuc+p47LZlXC2zDUfoF3cccdsLiodKLaPtN2bahwQ2lDOX7OTcjWiFVAVQnqQT5bR/A5rnWJ6896T2KaUb9qSWKaEQxDzZZpr2QKheJEIimeeOJ3B8XLGOir7WCmUy9dUh5FO/G4mnRnWSof8QM/AB/8IJw4Aa++qryR11+Ha9eUFyElhP1VX0ooFiGV6jAUoDq6H/6fUE7A3FhnW4XsvnQMBvxy2kPoz3FgEcJkdPQHyWTuBUBKies2OgyF2u5SrS7x8sufbW4LKqNWm6Cn2TwHM8E9Pr7Ss6ioSXk9jciFCyr8FNRsex7U62Tvh5kfQ82wWC4xOVNiYrHz8Oz98MZfgURZqc4W0tv3a2l8ug3Fanp02qDsSSxrCADTDBOPj1KpzPfcT0qJEAZLS283t+Vyl4hGhzv20413W8PB9CxOn1YhpHJZeQLlsnr+wAMtowEqN/HWW3DrFty8CZFWA0QQXiomlREIxASz93f+qK98XO2zPKzKab3V8hWarUN2/YPeBkSzJ5GygW2XqdWW8bwe8v/N/Ryk7HTmLWuIhYW3mJ//LouLb1Or5XXj3RZxMI3F1JRKWB89CsvL6vHsWXj22ZYRyeVUyaxtK4/CddU2oVaZdulxgXo0XLU9IHs/3DrSEhDUC9Qu0W4wdE5oz5JMHiUcTtJo1EiljpJMjuE4dVb7o6nyWpeRkQcAla8ol+fxPBspBa5rk89fplZb3nTjXTY7zblzJ3nuuXHOnTtJNju9+V9wwNl1YyGE+KQQ4vtCiItCiM/2P2KLmJqCV16Bd99Vj1NTnUbkvfdUTuLIEXXrUq/TvI0RgtwhlbBuJ9ygQ0wwMCiByqxmF9Gf/x5HEA4nGBq6i3vu+XFOnXqFer2A5zV8KXK1T/cxsdgIjz/+LKDyFdFohkzmBKGQBagO7mRybFPJ7SBZXize6FCpPagGY1eNhRDCBL4ITAEPAj8nhHhwN8+paUTuuEMZjcXFpjfRREoyy9AId25uFxPM3g9X7wZHKQ+sPcdFoznQqPkUxeL1FXO0DSOMWiZWqkEaRohHHvmlpiHI5S4RDseJRNIcPvwAo6MPMTLyPur14qbOSifLO9ltz+KvAhellO9KKW3geeCJXT4nRZAEF0JVSFlWx8uTMyr/YIfVJWyH1fPJmVY+Q0gwg3ktsHr4QxsRzQFGKcJ6SOk1S1yz2WkqlUVct94zb5FIHOXQofuYnX21uS2TGafRqHTsdzv5isD4tHOQk+W7bSyOAVfbnl/zt3UghPhFIcS3hRDfvnXr1s6cWZAEDzAMZTDuuAPCYSYuwtRLkCoqGfJUUT2fuNgKP6WKgPA7uPvw0F/AoWWIlbftN9Jo9iSGYSGEYGTkgaahmJ5+Cte1SadPYJqtwhLTjHLo0P0MDd25YuGenDyN59nYdnlLhAK32vgMOrtdOtvrXnvFfbaU8neA3wHVlLfdJwWocNSDD6pqqEAoMBpV2lDvfz9897tMXJQ9527nDqkKKQGQV8OPbEHLNLd7GlJN0XvyReWRPP9z/kvdDXsazb5FiXUFVU3t4R+AWCzDrVtvAnDkyPubR3Uv3CocdbZDC+p2ureDQUu2TUeD30FVqd1tY3ENuLvt+V3A9V06l5X89E/Db/6mSnQH+k83b6rnsVhnmW0bgYig1YBoXf2zw7B4uCUqGOCJlsWcuKi+NmZD2ZXunIhGsx9xXRsQLC6+TTY73eyVqNfzFArv4Tg1gjusYnGOZHJ01YV7YmJqyzq1t9r4DDq7KvchhAgBbwOPA+8B3wI+I6X83mrHbJmQ4HoI5D3yeVUNFYkomfKFBYjFyB5eZmZSeRKZZZWvmLi4UuKj4eczPKA41DZeVah/qTwcyqv3KQz54raeNhZbivbS9jgCwwiTTt+FZSWpVpcoFK7TPUrSNC0ikUMkk2MIAfX6ntN/WouBvgp31bOQUjpCiKeA/waYwL9by1DsOJcuwehoS44c1Ep+4wbZu2pNg9DelIeft+AlVhgSgBefgHpE9V0YHoQa4Fjq+GgVqmGoJ6Ax0JfVLtNdebbWZ7ne/TTbjERKl3z+GiBXbcaTUg0xajRKHfO2p6efAvaE/tO+5WAKCXazXrnyfB6uXAHb5tzfaYWaAuywSmqf+n9X/1HZ+zuNSDWuyms9Q+U2XBPc9vyGZnP0EhHsNgZaAmQPElz83dPyWn+YUChCOn2imdMAsO0yqdRRTp16ZUfOcpMM9NWll6S2MarZHzA5975v8txXfopzv/Uw2b/3WGdH9+ysUpWFdTXlddNtKCZnlJfhGr50udE2+wII11d/L00fVpbldz7XY1f3KJLu0FNru/rjKA1PXdK602hj4Y9RzZ5wmP7gVYpxl2jDpHjzItO3zpH9v061OrotSyW6DaNvU1432ftVCOraXVBIqccXn4BIHUopVFWU9O2EH0Zp6FncW482BgPA6jfghhFmZOQBXdK6C2hj4cuVz9w3h+EJLNdEGCZW1VHdmvarrY7u971PVUEZxoqmvGJczda+NQrn/s5KQcGXH4dqTCW0DU89VmNQ85PfwX2Tpzu9NQecbsGEtlf42Mf+CY8//uym+ym01tPm6WsshBBDQoj7emz/4Pac0g4zPg6VCrl4nbDbJkEeiXS6tv5+QbK7vSmvkIJaAuJlSBV6K9AujagKqMB7MPyKqPIQHLmlnktDdXxHqgx4dHMAWO3z1Z/7rtM5f1v9QYQwicdHefTRzzExMcXU1FlSqaPUasukUkfXNdxokLSehBClPq+fEEL85Qbf898LIX5ms+e0ZjWUEOLTwL8C5oUQYeDvSim/5b/874EPbfYH7xlOn4anniJTClGMu1gNoYzF2Fina/vYY2rcqus2x6lOXFT/upPdVgNsVH4iaNpbLYQuhYrBeqaqrBIO1LtGami2ENH1f10NtYcxME0L8EgkRhkenmi+spl+iu5mP8tKYNtqu66i6k8/z+IfAz8spfwh4OeB/yCEeNJ/bX98tXyl2cm5+/Gkix0xkHffjR0LtVzb6Wk4dw5GRlSvRVcF2XqS3SOLgFBhpiDcJP1PvxJThsIJa0OxrfS6YrVs+Z7DMIJ7WOVhJBKjmKZ1253Tg6j1JIRICiG+IoT4H0KI7woh2rXzQkKIc0KI7wgh/rMQIu4f88NCiPNCiD8XQvw3IcTRrTiXfsbClFLeAJBS/nfgJ4B/IoT4B+ynyPrUFBN/+AZTf/+PSf3gj1CLeJ2urZ8Eb45bNTsnGK0n2f34yxCrqtCTZ7Qa88I1qCZ9aY/984nuPbQxGBhM08IwLII/2vDwxJbM0B5Qraca8DeklB9Crb//t2hptv8A8DtSyg8CBeCX/QjQvwZ+Rkr5w8C/A/7pVpxIv6a8ohDiPinlOwBSyhtCiMeAPwJ+cCtOYC+xqmt76RIMt41q7MrATc6oHIWN8ihKcZXDqMVUiCro7H7iRTU5b+EwzcXLVbL7GNLvr9BsD0p+SDMAOE4NIUwsK0E8fmTLeicGVOtJAP9MCPExlKt1DAi6hK9KKb/u//8/Av8A+FPgA8D/59sUE7ixFSfSz1j8El1fMSllUQjxSeDTW3ECA0H3zO42XajsuMvMJNiWPxEP8MIq2R2vdHZ2g9ovk1NGZWFEhZ4MF+1V7ATaYAwEgVy5EENbetc/oFpP/ytwBJUOaAghZoGo/1qvNKgAviel/MhWn8iaxkJK+cYq2xvA7wXPhRCvb8fJ7Rn8JDgA8bjSh6pWyf5wmukfWVRy5AUVespnlKFI+N5ue7IbWqNYAYaKsHSoJf2xLvSCp9nnCGH4Y1HLW37Xv5VCgztEGpj3DcVPAPe0vXZcCPERKeXrwM8BrwHfB44E2/2w1ANbIaO0VX0W0f67DDDdM7snJuDXf52ZH66tmMPtGap/op0g2d2dCI/UIVFCJb716FXNgaRVXSCEQSgUwzBChMMxIpH0oC3s28HvAR8WQnwb5WW81fbam8ApIcR3gGHgS/4QuZ8B/rkQ4g3gL4DJrTiRrRIS3P/LXDCjO+CZZ8hZZaJdKuUhxx+l2kZ7srtbTypah0RZpUHmRtf4+YFHob0Kzb6i3VU2kdIlnT6OECFSqS0p4hlIpJRJ/3EBWC1q03MEtZTyL4CP9dj+d2/nnA5WB/f0tBIHHB9Xj9ObaMaZnoaHH4bf+A0ySyuroKyqSlYHnd1lv7N7/ogqka1FV45i/fhX1LGG1zt/YXgbCFNpNANGMGfbMFT4qVC4vq7Ec79ubN2tvbWsy1gIIVZYML8qqvl0q05o22gTDGR4WD0+9dTGDMYzz8CTT8J3vgP0nsMd9uCjX1Wd3cUhqCQgWlb5Cdf340LuylGsC4cBf2Z32AHDafu5nh+m0twee/8qPTCoULr6g0jpYllDeJ6D59VxnArl8iIzM2dWXeD7dWMPUrf2oLDeJejLQoh/JBQxIcS/Bj7f9vrf3oZz21qCXolEQsV8Egn1/MyZ9R3/Ez8Bv/EbUKs1N602h/vRrymZ8iPzkM5BqtLKaURrEKvA08+pfYIO754RJt/D8MzuFzSawUZKh5aKrEu9vkxrXILAdassLGRXXeDbu7GFEFhWQmm5zZxZ1+uajbPenMWPAP8cmAFSqKTLjwUvSik3pFGyK3T3SoCqbJqd7X/s3/pb8OqrPV8KJD96Eczibmc1GfPhBVg4okpvJeAFfxl/RnczAa7vjve0CY8AACAASURBVDX7gl5pThVrNQzlddh2nkjkWE85jmD0ajvt3dj9XtdsnPV6Fg2gCsRQlU+XZKfa194nEAJsp1KBEyf6H/sHf7CpH7kRGfOPf6Wtw9ukaSSSJd/raJtzodkE+rPbo3T61IZhYZohpPSw7RLLy+9y7drrK7yLft3YA9qtvadZr7H4FspYPAJ8FPg5IcR/3raz2g5On24NMpJSPdq22t4Px+m/Tw965TQ8szViNXu/6vB+7mnVh/HIf4e7rvkyIA3VvJcqwVAewq6/3eYg1J5tLdpQ7FEEhmFiWUkMw0IIEyEMXNdpG6uqJud1h6MmJ0+vKVPe73XNxlnXWFUhxIellN/u2va3pZT/YdvObBVua6xqMD51dlZ5FMH41H6YplKi3QS9puMBvPQpPxwlwXRUHiPsqZzHzOTKEttyHKpRlST3DPQCuBn0Z7ZnUGKBBqGQhWlGqVYXeuwTRgjRUUrbLv2RzU6v2Y3d7/VdYM9cgb4Kx3MoOZB/I6V8tu8xegZ3H6an4Wd+ZmUIa5MEE/PKCVoS2UKFmRJFGFlqaU25foOfE/IVanXO4vbRn98uIzDNCInECJXKAqYZodGotnkSnetRKBQjlTpGJDJErbbM00+/u/OnvHXsiatPCGECbwN/DbiGihz9nJTywlrH6YLMfpw5A8ePQ2hr+hdnJpWH0LxsfIMhBdgx5W1MXIS7LqvS24bVpkjrHyOC4Ri9hmRo1kZ/XrtGLHaYUChCNJpmeHiCT3/6BY4e/RCGEcI0I4TDMV9tNkDgug2Wl99hbu4NarXlg1f6KsQnEeIrCPGu//jJLXjXvwpclFK+63d8Pw880ecYbSz64o9d5cQJVWobiaw197Ev80faKp0C/ByfE1Khquz98OYHlLcRarReX4FQDYAazd5GXcBSurz//T+NbZeYnf0z/tN/+klmZ8/junWEUEuRKqkNkEjpNh8bjSrPP/9TfOlLDx8Mo6EMwxeBo8CS//jFLTAYx4Crbc+v+dvWRBuLfgRVVOm08jAsC4zNf2xBpVMvDE+FoGYmW6Wy7moOjW88dLPeJtAGdoeRmKZFo2Hzl3/5PI1GWW2VHuAhpdsMQ7X3WgTHtrvhQpgsLV08KA12p4E6EMTAK/7z283S97r17Put0EvNagTSIBcuqIT43BwMDcGxY3DPPc1Z3BvF9PwwUg8++lUVgsodAuEqD6TXrnJPRD4HHG0wdhCB69Zx3YrvKQiEMJreBIDnNXD9ccUAhmFiGCFCoRhBRZRhGL4irXNQGuzGaRmKgIq//Xa4Btzd9vwu4Hq/g7Sx6EW7NMixY2qc6sICvPeeUp49e1bNtLjvvg2/9ZF51ekdChpY/WqosTnV+Q1+H0a7QZCr/B+CPiaNZg/Ta+xCJ0Hlk2GEMAyLZPIOQqEonucCHqapGpak9Pz8xoFosLsEdA9ajvvbb4dvARNCiHEhhAX8LPDH/Q7SxqIX3dIgY2MqZ/Hgg/DKK6rcdnwcwuEN5y8mZ5R3kc7BHTfh8KLSjXr85c59pKG8CwHNJHjUFyk03LY37PcXbEuCm5trF9FobgshDEwzsubrx49/lF/7tSq//us2P/uzf8Tw8AShUALDUMdKCa7bwHVtGo0KCwtvEYmkdvC32BXOABFaBiPuP78tl0qqxNBTwH9DyZx/eT3zLrSx6EWQ1G6nWxrk9GnI5VSD3wZYTU+qXTJk4iIcuaWMiuGBZcPwskqAG+4Gk9p+cjzUAFdrTLXQobwdQc2oiGIYJqpiUxFMwwOwrFRHs9zExBSnTr3C6dM3+fSnX2Bk5AeQsoHnOYAKT3leg1Jpbn/nLaT8U+BXUGNRh/3HX/G33+Zby5eklA9IKe+TUq5rRrfus+jFyZOtMar5vMpXVKuQTMLv/36rke/hh5sKtFtN0I9Rj7Qm6bkhiJegktxE3kL3aLTQn8O2YhgWUrp+fsLENEO4bgPwfIMh/bGpBkNDx/nrf/3/YWJiqq2J7hKZzHhHE92XvvQwS0sXkdLx+zTGMIyVjXp7nIG+8rZq+NH+IhijWiopQyGEqoBKJtX2s2fVfm++uam379XV3S1GeP2oasgLqp2EVE5MJbHJBPdAX6aaQcLzbIQwMAzLT1zXUdIeKu8Qiw3zxBO/u6Lbenr6KQzD6pAUh7NMTExh2wVGRt6HaAv7SikPQt5iz6CNRS+CMaqf+Yx6HomovEU6rTSlfuVXlOfRaKz9Pj3I3q+6sw1X5SCKKfWcl9TrM5P+oKSE32fh+MYhaN5rL73tVV2o6Y/+vLadROJOkslR5uf/Es9rYBgmoVC06RF0K8m2S4oDWFYC26a5XyYzTrF4o/k6aGHAnUYbi9WYmoJMBu69FwoF5WFcuaI8jFpNJcA3wcwkzbndoB5t4OXHoRFRrzkWIJRhkPgS5aBUXDw1PMkJ00peaykQzf/f3rkHx3We5/15z9kb9oI7AZCCSEImKInRxfJYlkRXIQu6saA40lQZ17GdWiNnxiMpalVPC1uuUsmy24lczjjVjBy7mUYOW0dOnFipMrY3jiValivKimXHjCPTJiiBpCReQBDEbQHs5Zyvf7znwzm72BsWl7OLfX8zGHB3z+5+WALnOd97ed46I5Hg0nLelIdgOL1JMzOnYRghLC7mWy9XshTfu3cEyeQDyGT4/mx2XowBNxhJcJdjYAA4dQp4/XVgdpZ3Enr4kVI1NedNdbCjrJdgFpjsdkXE8uwevGEoLQbK4B1HIMfJ7sgCsP/7jiOtIPiO+3fBVVA2crkFWFYWRAFYVgbp9HRecrqSpfjg4DCGh59EIrEVi4uXkEhsxfDwk34bAzYVsrMox7ZtwPe/7972FgPkcit2ox3dBSyGgZlWPtHH5oBImmdcKLgiYlrOACSTQ1AKbCoI5eiFcnYaBs/AuPNZznm89F7IDqMaSvWwyOdWJeQ0x1lFHw0Go7hw4RfI5Rbh/YCJ2PKDCIhEuvNCUdXsHAYHh0UcfER2FqVIJoFvfrP040qtaM6FzlWQDdgEZELApQ5gspVFofuiOygppgcecT8SbIN3Ey3zPOPCtJ3tvc3H6uR4Lsi9FIY06pWmlFAUuy2UgBNohhGG9xRCZCAUakMut4BcbgGFH6htZ2GaIbS2bkc83pOXnJadw8ZBRE8R0TgRrWjCqYhFKQ4eBCyLQ026GoqoZhPBI3uBrAFkIs7J3NkBpKPA9f/ITXl6UFI4DYQWAThvGcrwfW0zvBPpugj0jAPdE/x6mnCaQ1jiF1WEQjPGUsIgglEFXPpqGCaCwQhMM4xQKIHe3nfCttNLsyr0DoQhBAIt6OrajUikbVlyOr9sti5mT2xm/gzAis0IJQxVirExroJaXGTRUGrFDXhepjrYghwKMBXvDhT4xH5ywLH6cAYfXejhZHdsBojP845juh2YjwIxT1jXO6J1dJfT2V1KKBxbEctEaRfbzYyE59YQBSJzqVeCp9ulMTc3jlwuA+6nMEAUdHoscrDtDHK5BSilloWYKpXNNjOPPUa3gY0DB8A2HwcffXR1TXlKqReJaOdKnyfXoKUYGOBqKMPg3MQqmxfbnQ5sr4mgIs5PTHXw7cETwN3/m/2j2qaAxDyf30JZIJLictpiI1p1iCtT2lEBIBaatik070lTZoCsGbFYDwAbtm0v2XDMzb0F/eEqZcO2s8jlMp7BRsDs7NvLQkzeslkiQigUaxajwLI4QrHMoty5f8MRsfCinWYHBoCLF7n6qaeHxUJjmjWFovYecUpgDWdH4Vzdtyy4uwNNsYqp+DyHporZhOhy3Fw5Ow8FzMWB2bYVL725aFYhBVDtD09kIBxuQ2vr5Y5tRw7aPjwf5TzG9d3xeB9CodiyENPU1BiCwXx7nSYxCqzEelmU14SEoTTaaTYUAjo7eYYFEf97YoK9oXQoSikgnV7Ryw+eYAvy//frLBimxUJh2u5cbk37peUzuLNBoOcC7zwKmerg8tlyf+tmzhnN2tQnwwo0/WdT3ZZLKRup1Dnnlr0UktL24+68YMY0I2hr60c43IZMJoXnnnsIzz//ECYmjjtpQAO5XG6pNwOQhjuHAfCOwstaWJTXhIiFxus0C7jfu7qAvj4WjbNngXPnWCyIVhya2vdDYNvZylYfegZ3BrzDyHpCTsXQ4lLuZNc+DUy1s0hlausn3Jw0vUB4oSXDP7boKHMkBbCwMIniHvn5fxetrSwUADvHTk6OOrYdpnNfGtksj1OIx3uk4c5lDBx68jagrIVFeU2IWGjGxlgQvGin2YEBYHSUQ1OBAJfM1pjDGDyxXByKHaOT3eVERbNzjHcspSDlluemIiIWwnJ4noQBIhOJRB8WFiaxuDhV8viFhYsoMZpr2T0zM6cBbEck0oa5ubNQSsEwgktd3Zpcbg6Li0GphnI5CM5ZACwYa2JRTkRfB7AfQDcRvQXgUaXUn1Z6nm9iQUQHAfwW+AL6dQD3KKVK/3auNwMDrtOsZn6e51iMjAB33cX3BYNuM14ut6KmvJVQjagAnNw+egMQTDs5C32lbGMpItB6CRh2mmWTt0Mqg4Q8iNi3iZvoCMPDT+KZZz7ihJeWN96ZZrjizsPz6lAKSKXOOdbiGRB5S2o5DGXbFiKRDjz44Btr80NtAh59VP3dY4/R72Ptq6E+XMvzfLMoJ6LfAHBYKZUjoi8AgFLq05Wet24W5d6cRTTKQpHJsKHg8DCHolIpvs8wWCQy/vtrHPoYMNEJLMTBnd1OAx/AAnLVMU5q6x3KzjHgpX8BZGV3wTSlaLo5BZ2U5vBTAPF4L9rbd+LkyRdQKodBFIRS1ZlommYIhhFCLreIHTtuxcLCRVy8eMLp0+BfVNu2QGSiv/+mRrIbr4WG/m3zrRpKKfX3zsQmAPgReA6sf2in2a1bgUuX3PGpenbFnj1Afz+wfXvdCAWQ379hKDYZDGS5OqpzCnhrB+cztMPt0Rt4JoZhyeS8xv7TXTncLAd4RYDnTvBW07LSmJo6hYmJ0Tx310J4Bnaw6vdrbe3Hjh234u67D+PAgccRDrdCqRwsy4JtW1DKRjickBxFnVMvpbMfB1By5BURfYKIXiWiVy9cuLB+qxge5rGpb7zhjk/VjIywQJw5s+qei7WkXP/GRJdrTqj7NQwLSLUCrdNcidU0FApDkwkFAGzZsgeRSLsTAvL+6duwLD2DIoBMZhqxWJ/ncfIcTyAidHdfjUikA+VOIbxjCOQlqwcHh3HnnU9hy5Y9TjiKX6twvoVQf6xrGIqIngPQV+Shh5VSzzrHPAzg3QDuUlUsZkMm5ZUimeTchXaeXUOqGYhU6nnf+BD3bRi2YzxosAeVZXI1VWKOrUAAvqa8sIU9pkJZ4Hxvk9mDNKFIaCKRdidpbaClpd2pZtIfCAGwYRghEAE9PdfmmQGGQgnE470gcqfTjY4m8fzzD+H8ee+0SB3i4u+9vdfhwIHHRQiYhv7t83WsKhHdDeBeAAeUUvOVjgd8FgtgTUepaoG40MPjUyMp197DNpfP5i7FD251+jfIqXwy+DvZjnAAaJtmwcgEOVSVCfEuY7rNmY1RLRZARo3T+uqFRl77qsjvf1j+GD8eCsXR1bUbi4vTmJ4+BcMIorv7qqWS1kKDv6985XpcvHgCtp2BUgpEBMMIoatrF+699+gG/FwNQ0P/5vl2TUlEtwH4NIA7qhWKDcXbzT00xLfXEG3RMZtwrDuIk9TpsBsuOrK3utfa90Pg3/wlsP003zYtdqgFWHQsA5hOuBYhB55jIUrMOtbnVUDa/LDRhQIQu48SGAYPUgmF2pzy1gBaWjrR1bWrrBPsgQOPI5HoQ0fHLvT2Xo+Ojl1IJPpw4MDj/vwgwrrgZwDiSQAJAN8jop8R0Vd8XEs+ujLq7Fm3Ge+BB/j+mRlOfq8S78Q823SdaFNxfjyYdT2jqkH7SiXmgPgssBgFQIDhJLEtZ0ehdyv6+LZZwMyi+AnUGevaOs1hKwANfm0klIZDR6YZwcLCBCYmfgHTDOHGG+9HS0sXyims2Is3B76GoWphQ8JQQ0PLey5SKVckzp7lCXrpdM19Fk88yBVKBOBiF1/hk+KwUc847wISs8XtPcpx6GPAW/189W/oaXtOeKprknccUx0cklIKuNjtut8W5DzRe4F7N0JZXmMmiNJiYaN+yiWqYdOLnndAe/njiAJQyoJhsH9TNMpd1JzfUIhEOvIGEokQ1ExD/9Y10p/3xjE2xr0WXrJZ4OWXgdde467uWKxmU0GAk9h62FF8DgC5nlGZCvYe5dh7xJmwh3zDwvAiJ7ZnExxSutANTGxhwVomFODnTLXxDicdLiEUyvOd86ONQQP/ybrlr5Wo1l5X5QlFIBDB5OQopqbeQCo1jsXFKXGDFQCIWBRnYICb8jTT08Dp09yM198PdHcDc3MsGF6xWMFwJH1SzwSBUJon3hkKCGTyHWVXyuAJYMsFfi1lcHls2zSQjrhhr1TcSYArnpuR9yMoFizD5sa9+SgnwYvi5EtDWR7ItPMUN/11TXCvR965qp6swYutq17WVhZWdcNYaUdleQdJwzDR0tIF0wxjevr00qxsLqlNI52eXjpW3GCbF/GGKsbICOcoAN5h6N6Kbdv4zNrbC8TjbljqlVc4JGUVn0lcjEL/p85J4APfKi4Q5cpqiz124DlOnhtWvhFhq5N3sEwnYe38W3mqHQNOjsMGYAd4hkal82h8jt9nfAuQSAGXOt2EONmAMrGUHN8wyhX+aIqNVa3TXQeRgWuu+TCOHfsmLKu67mmXwi3f8h9yy5Y9eOutV8BiZOQdNzd3fskIUNxgmxfZWRSjsJvbtoEdO4A2zyW2NhkcGWG/qBUIhUYnmR98gr+XEgpdNaW7sJO38/2lHgPcaic9+6J73G3CMy2nH8Np3iPHR0o39tkGYDmXEZbTv7F08i24KjcU5z/mojziNeXsYCxnsibZHAKru5Nwg41VbWvbgWuv/Shisb6ifk0rg39IbcVBZGLv3hHYdgYcllKwbcuxADGWJtxlMqm6cYMdHU3i0KEhPPHEAA4dGsLo6NpWKwrLkQR3NZRLeB8+zL0Xv/wl5zXW+PM89LHlsy108hso/VhhYlwLi2GxGOjQUus05yQWYgBZTurBmWETnWUBsA234c9QnIzX95HNz8kF+bUVeYY8eTyqlryoNlI0Su0uKu066k3YAITD7QD0SNKV7ixKQ2Sip+fXcO+9R/HlL1+PyckTUCoH0wwjFutFNruIXG4OkUhH3bjBesewNljivQ5/s6pHwlDVUBiWGh/nr9FR14FW5yvWWCyWBht58JbVlnvMS2HYa8sELzUTAXrHgcQ0MHol70RIsX9U6xywmGVhsZ0fz3J+3VvmOQ9iB9hjihyBAC3v3ciG4M+fSSmhaDA6OnZhdvYMLGtxyUqcL/JW+7vG1Q9KAU88MYBwuBXhcCKv+ikYDOGOO56uq5OwdwwrAIRCMWQyfH89rXOzIWJRDTosdfAg8ItfcGiq0J5cT9BbY0pNzdOjWMs9Vkg1tufekl4AiKQBNQ3MtQKhBd51mDbPCV+IFpTWlsq7VnuCXs+cQTFvqDoXk0AgCqVyiETaMD19CkrZjn9TEJa1eiNLwwjCtnOwrAwikU5ks/PgXosQFhcv1c1OopCpqTFEIvmzZyTxvv6IWBQjmWRhGBvjyqiRERaM4WEOSb3yCovFOuwkCqk0NW8lE/WqoZg4mTbQ/+by0JYWFoCT3JMFs6PyKCUEyglbmZ5k+HpQ7P0LBaOOhALg8aWmGQbAJ3Z2iLWdTuuQk+jWDS7V1i2T83oB2HYOgUAk7wodAKLRLtx3X/3adLS3D2B29myeM64k3tcfSXAXUq57G2ABWYfcRCkGTyxPVnu7sEs9Vivekl6F8j0f3l6RcDrf+bYYZHPVl7ZQjywAsXkOiYUXfLIRIc9XnaFUDrFYDzKZFMLhVoRCCShlO7begGEYiMX60NNzTd4woUICgeiS6ADuVDzARjye70bQCFfoOhmfyaTqLvG+mZGdRSGlZnEfPMg7i4EBzldY1oYKRikBqHai3kreq5qRrqO7gPkWLpM1LBYqM8d26QCKXsXHHPfb7jSL24NP5L/e1z+8wih8HZe6rg5CMNiCzs5dSKdnkUhsxfvf/0UAwHPPPYTJyeMAgO7uPThw4HE8++w9UKr0ziIQCCEW64Nt55BKnYdlpUFkIhiMwjTzXSQb4Qqdw2JP4siRg5iaOlm34bLNhohFIeVmcQPA/v3ASy+t2zjVeqCSAHkrq1qngLkEMN3OSe+c6TQEklsNBcVhqsQc3yyWV9HvZ2YBGNz/kYfrer1UbbUiofBVVMgZRZpB5XARwTBMbNv2Hpw//zOk07NYXLyEM2d+jH37Hll2QhwdTeY1zS17NTJhmmHYdgaGEUJn5+BS9dD119+No0cPIZNBXlVRI1yhDw4OizhsMBKGKqSwextwZ3Enk8ChQ0BPDxCJ5B9jGPzVBHhNEFvSHEbqmOQTfHyWDQsNBYQyQEuKb4fTlcNa4TQLjKHPpx47EcPTaU0l+j6MHBAsNv3PB6EIBKLo63sXurquQkfHAD70oWewf/+jqPwnxz0Op069gMXFWRhGCJnMPF588fP4wQ8+t+zoI0cOIhLpLvm6kUg7tmzZU9Tob9++R8QAUKga2VkUUlgmq2dxj4y4IaqODu7iBni2RTgMXHkl3/7JT/xZ9wZSqpw3HWZ3WvJorQIw28phqkqDnW5+GXhxv1N66wkx6V3JVCu76SoDMLJOKW8ASwJi1k3LkAEiwszM2+jp2bMUItEn4Zde+oJTeeQcbbDxFjfFebPuNpSyYJoBWFYOP/rRH2Hfvkfy3mlqagzxeA+CwQguXXoD3p1LS0sXIpG2Ze/vRa7QhWoRsSjEWyZ78iTvKHQ11P33c1/F8eNs7xEO823v5LzOTmBy0q/VbwilynnDaf5eeP+W8ercc/f9kL//6Bbu4YACIvMsFpkg+1jFZoCER4wmW4F0lKuosoZrY5LHBuc2wuE2xOO9SyEd78l4375HsG/fI0uNZbOz56CUtWThYZohWFbaXbrKAQiAyEQ6PbfsvXRlUCTShq1bb8Di4jRmZ89AKRu9vddJLF9YM0QsiqHLZAtpbQWOHeNwk2nyjiObBQIB7uiORjlENTvL99cJtY5sLUWpct6bXwaO3rC6Ut59P3RFo3DdCxGeJKhZDAO5iGOKmOMGwXUrvV0BmcwMQqErijaKjY4mncTsGHg3YcG2WRyITBSqmnZYUMpCOBxf9l57944gmXxgKe9gGAEkEn0SThLWHBGLlaKU263t/dIJ8K1ba/KJWi+8yeg8/6hVlNiWq5jadrb2WeLFnud9bqH1SSrOm4Zgjl1vAbZez5lYRTbO27NgAvD+X1bjTsgn9osXjyMW68krQ/XaVBCZmJ4+BYAQi/Vhfn4cSlkgCqGlpdOZj83ui5aVA2Dj5ps/mfdaWnTC4VYohbxGOgA4dGgIU1NjaG8fkB2GsGpELFbCzAwbCo6PcxjKdC5jbRu49lrOb5w8WVeVUt5kNMDfM879q9ldlKqYqqWUt1pBK9zR5Jx8RdwTnbGNIqf0CiGoUCgBwEQ2OwPAQCDQgmh0C4gImUwK8/PjzpHuq/LAoGLZdCaTSSGTeR3BYAyjo0kMDg7n2VRcvHjcsQFXyGZTaG+/Ykk82tp2AjAcwQBCoShuvvmTS/kKr+jozmvbzuD227+EwcHhZY/Pzp5FMvkAANltCLUjYrESBga4SW/3br59/DjvIsLOUIhYjDu764hK3lJrRbWhrmLHHdnLSe1UjCuqTEc0vIKmn5cOAnaEu8pDaa64CmcI6Qgwk1BLJohLlBEKogCCwRZ0de1GJpNCInED7r77sLvW0SSeffbjRYUhEAghm7VRuhSW77es7NKJ2mtTwb0OLBaWlUYk0galtmNu7m0sLl5Cb++1JXcDlbyRxDtJWA9ELFZCYaXUwgLnL3RlFMAltXPLE5F+UclbaiWUEoRqdwaljlto4QQ2OVVNlgHMxbnSaXQX8Pz7eMqfYXFprm3ynA0jy9+zAYXFFqevo+wugkBkwDSDyOUyICJEoz0lO4Cff/4hLCxMONPpTKdPgrPllpWDYZgwzShyuXmnKc7d07DZH8G2s5iePo1nnvkI2tq2I5udRygUc/oueEKU7q42zSD6+2/JE6xiVPJGEu8kYT1ojsaAtaJwzkU8zglt75yL9nb/1leEldh3lKPcXA1vqIvA3w2L79fPPfQx4BsfAmbjfFL3HqcNCA3F9+vZ4Zkgv8fFTgCORfpMO7+3TQBMIJJyy2lLCwUhEmnH/v2fxY4dv454fCt6e69Bd/fVAOyS/QUTE8cBmDAMwxEGbZ+rljqflcphy5Zr0Nt7Hbq6roRphhAIRJ1yWAssUEFkMinMzZ3H4uIUMpkUYrEeKJWDUnZZwSpGe/tAXuktkN95XelxQagF2VmsFG+llPaR0pVQ8/Pch/HRjwLf+EZdVERVa99RiXK5j3Khrrw5Gs7JfLqNR72G03ycDU4l284AJkV8SrZMZ0aGnp0BIOu8hmnx44l5YD7BzwlYhFzRZguFdHoWR48eWlGVUOGEXK+Bnx5AZBgBzM9PQqksWBgCToVTZulxwEYgEEEk0gHTDCEa7cLU1Els2bKHbeIzbOlRbRK6sAKqsPO60uOCUAsiFrWinWlTKeD8eS6jNU3OZ3z0o/z985+vixzGWvhHlROEcqEur8gELLfhbs6pZJpNACBA2Z654RbboS/E+D1M53k2sLR7yDlzNBAOI2BnkSEbVplKJXZrDS3F7b3VRKWqhTo7d2Ni4pgzrsSAUvZS6Mi2bRhGAJaVxcLCOFpaehCNdiIe78GFC78AABAFQcQ7kVisF8FgFIuLl1bt6FrJG0m8k4T1QCbl1YLeUYRCwNQUi4VSnK/o7uYdBREwMcE5jczqZw/4TbmJfbpKqXDm9/B3gO/8cKmKHgAAFedJREFUpjsfYzEMzLQBcKbt6XBTyzzPxlDk7CRM3knEZ50RsAZwqd3jNeVg2EDbFJANAPNxZ+eidwNFQlLt7e8AYOP2279U1aQ1TnDfg3R6Fradc2y9s0vW4XrnYdsWTDOIaHQLHnzwDYyOJvHMMx9BJpNCIBBBLNaLSKTNSaJvrZiTEDYtDW17KTmLWtC2H7mcKxQAl9OOj7OAzMwALS11sbNYC8rlPspZpWsb88Ww0xdBzthWZ/xq+zRP5IvO82OWs2OIppznRRwrDwv5PlDO7PDZBBC0gF9/wRUfIP9YhjA3dxbt7TvzqoWICKFQbGnX4WVwcBh33vlV9PffhNbWy9DffxNaWjpBFEAut4hsdgG53CKUspHLpZdyAoODw7jrrqfR1nY5EonLEA63io220PBIGKoWtDPt6CicwQKuYBDxTiIQALZtA06soX+4j1TKfZQKde09Ajx7J1c8kQLgTKCFwaISdpwtMmFHAAjodhrsMkFHOBZ4sFIw6zrXzsWd5jsAw0nC4BsGTg4Ck90GZlqc7Y/H5oPIXDpZf+c791ddLVTonfTlL1+P+fkJ6NJYHpCYhWEE84RAQkHCZkPEohZ0v0U6nT8tj8h1ng06cwIMo66a9FZDLbmPwRNAPMUmg8rgvEV8DphJOJVVjlhYJouJ6WmYDmZ5p3Lf/ywIgxEhPAlkggoJlcDgqTQQMrH3pQySd9BS8tvFgGmG0NW1C4ODw6uatMYN+1zhpCfXAUBra/8yIWhWk75q8kFC4yFhqFoYGXF3D4GAKxbBIDfpBQLcqKctQJqcdBjongB6xtmWI5x2+yV0WIuc0tiYp0XF2w+yFAYLERSxUNgBwt7Xe/mztiwMnoli+Ofb0Z1qgaGAgGWgvf0d6OrajUSiDwcOPM6vtYpJa+n0DNradiAQiMAwAgiFEujoeMc6fGqNiWuQeDave3x0NOn30oRVImJRC8PDwN138791MjsQYKEwDOAP/gC4/HKujrLt5TWYTYZ3/KrGtIHucTfP0TUJtCxwHqNYP8hSXiQbxmLUQGKWMPxyJwbHW7nPxbaBtjYMjrfi3uQO/M7hXvR3XINifRSDg8M1z3Fobx+AYQTR1bUbPT3XoqtrNwwjKD0MDtXmg4TGQ8JQtaCHIG3fzmJx9izblLe0sOXHCy9w4jsYzA9NLS66ISlv+GqTU8ql9v3fXd7hXa4fZPAEMHjKAq7+NeC3fxs4/QLv3jo7gUAAoy1v48h1ZzC1xUR739XY+1uPlxSAWkNE0sNQHuke37yIWNRC4ZxuIuDUKT759/ezeExPs5CEw+7zdB7Dtl1n2iYQjGobAyvmRHTIb3wcuPFG4JFHlsqYR3fmkLyRYOQCiCwozGYm18U8TxLX5VlNPkiob6TPohYGBvhqVoeXjh/nHIZS7D4L8M7izBk+welwlFJsDzI4CFy8yJVShSNchdJow0bDAG66CTh8GBgaAs6exaEDb2M2kkHIcj7rYBCZnZdJX8MG43W8LdfD0qQ0dDxacha1UDinO+2U9Hh3ET09QCLBJzbda0EEnDvHoZNjx9znCeXRM0OyWf48s1m3eGBsDIhGMRVNI2g5v86GAaTTEv7wgdXkg4T6RsJQtVDoPhsIcD7CMICf/5xFo60NeNe7+Nj77+cwFcC7jLEx/9beyCi1tGvAzp18n1PG3D4fzt9ZhMMS/vCJZi0Z3uzIzqIWCt1ne3r4yteyWAzSad5B7N/Px587xye4YHDTdHRvKNz5BmccHIvxiJNQdsqY977WBttQyFAOyraQ6W7bkMTz6GgShw4N4YknBnDo0JCUiAqbFhGLWhke5pj5G2/w1e3WrUv1/giHgb4+roo6eJAFIhAQoVgLwmEuLtA4wj1Igxg+0oGEHcXitk4ktg6ue/hDegqEZkIS3GtBYcIb4KvgS5f4+8QEC0U63RTVT2tONMrWKXpuSCrF4nzY38T1oUNDyyp/xCxQKIMkuJuewoQ3wLcTCTYVXFgonsxu8ma9qggGgauuyh8wFY3WRXf81NQYgsFo3n2SVBc2KyIWa4G2/0ileOeQSrnW5fE45zG8ZoOawtveaqpmxTDc8ljT5K9iQqwT3D4iE+mEZkLEYi0oTHhv3cpJ744OnnGhzQT1STAU4q9Egru+TZMf2ySGgyuCiH9+gD+bSIQ/k2CQxXP37nwhPn+edxU//Sl/vn193GuR3Pg8wWo8pgSh0ZCcxXoxMMAnwTffdEXCtrlHoKeHZ3VHo9yNfO4ch1kuXWo+wejr4+T/228vd++Nx4EdO1ggdL4nnebO+dlZPlYpoLeXxffJJ92RtxuE67Aq3dxCRRo67ixisV4MDQGvvOKW0wIsBIYB7NoFdHXxFfLOnZy8/fa3OXTVbJgm7yZ0qEkpVywSCRYTPd/85EmeRDg97Y6x1X0Xl11WF0nvUohtt4AGFwvfw1BE9J+ISBFRt99rWVN0HkNj2/y1dStfFeuy25ER4OWX+Sra8P2/Y+OxLBYIHZoLhYB3vINDUJkM7yKI+Hsux4KaTrufldOtXS9J72JIia2wGfD17ERElwP4VwBO+7mOdWF4GNizh69+czm++t2+3e0+TiZ593HXXRyCOXfOzV00G/PzHJ4D+Oc/f55vF/alRCIsDOGwG65zurXrJeldDLHtFjYDfl/K/hGAT6HIxORNweOPcxjliis4T3HmDJsHjo0BH/84u9PaNl9B61kY3oazZkFXimkzxnSaRTVQ4EbT3s6C2tbGx+Vy/Lm1tfFzR+ozsSwltsJmwDexIKI7ALytlDpaxbGfIKJXiejVCxcubMDq1ghdJRUKud5QO3ZwUls36oXDbmI3m23O3UXQMxkpm2WRaG3lnIW3HDkUAj7zGXbt7ejg0FNnJ9/2IbldLVJiK2wG1jXBTUTPAegr8tDDAP4zgN9QSk0T0UkA71ZKTVR6zYZJcHtxbLSX5l/8/OcsCKEQV/K8/rorGNoypBprENN052I0Mt5qMaWA667jXRnAdim6EGBkpG4FoRxi2y04NPRVoC/VUER0LYDnAejLrX4AZwC8Ryl1rtxzG1IsKs2/OHOGcxYAC0p7OyfBp6f9W7NfxGLAX/1VQ4pCOaTEVoCIxRosotl2FtPTHJbSVhbz81zl09PDIqGvoj/4QQ6/NDre4U/lxskScQmsTz0TgrDONLRY+J3gbg4K7UACAd5p7Nrldnw/9RRw9CiX0x4+zCfKT32qdP6invMahsFCmEiwUGiRKBQK78+gO7l1g91BqRQShHqiLoYfKaV2+r2GdUUnur3x9y9+sfKV8yOPcMjq61/P7+wm4pNxNlt/LrbasiOXY6sOpbg0OJ3mHVUxkdPioivB6rhnQhCaFdlZbBTe+RcjIywcAwOVfY2+9jXgW9/iE3AoxFffgUD5cE4lvNVHa422NNGmiPPz3G/y9NMcZvOu3WvvoW079HPqtGdCEJoVEYuNJpkE7rmHrUDOnOHv99xTXjCGh/kq3bL4K5utfX53NOqaG64X2Sx3pKdSbv/D8DDw1a9yTkKLhGG4VVCJBJfLep8jCELdUBdhqKbioYeAyUnXgtuy+PZDD5UOSyWT3MxXqkxWh6SqIZPhk7NOOq81Osx09izP8fjkJ92fa3iYd1ORiJvsB7hje26O8zcNXCIrCJsZ2VlsNMePu0Khk7qGwfeX4uBBnj1dimBw5Qlvy3JDWzr5HI0Cjz3Gpbu1ohT/TKEQd68fOpS/axob4/fxou3cvcl9QRDqChGLRmBszN0FlAofaWuMSqKRy7mzMwwDuPxy4Mor+Yr/r/+ak+pPP+3mD2pBzyGPxZZXNpWaKig5CkGoa0QsNpqtWznfMD/Pu4Vslk+uu3eXfs7AQL5IFAqC7vbu6uKTdDBYPiehfagA7vco7GvQJoj9/fzYSvMb3mR1YWVTsamCkqMQhLpHxGIjSSZZJLyT8bJZTu5qe4tijIzwMYDb3KaJRl2PpOuuA555Bnj2WQ7tlIMIuPpqd7aGN/STTLJt+vi4O4tjpeiZ2YW7hmJTBaUBTxDqHklwbyQHD3JsPpFwp78FAhwKKneyHB7mXov77wdOn2bBiMW4ae+RR4o/5777gM9+tny3NMBi9fLLvHsZGAD27+c8g3aC1XM4KlFYyjs1xTucYruG4WERB0FoMOrC7mMlNKTdh6bQIwrgE+ylS5zcXUuGhoDR0fxxpRrT5N1Iby+HiEyTT+zpNIeoOjr46/RpFo1yyXWNTpLrUthQCLjlFqlsEgSXOrZdqIzsLDaSgYF8jyhg/ZK7Y2McirIsNikstNnQ8zX0XAiAdznZLJfydnTwsKbz512xCASKu+EGAiwOeheyYwd/r9MRp4IgrBzJWWwkG5nc1VVH27bxmNJEgncP0Sgnr/WJPRDgL51H0eGn8+c577B7Nz8nHAZaWtxjTZNLbHfuZKEoNg1QEIRNg4jFRrKS5K4eu+q1BCl2Xym8wtTayp3Tl1/O5bHasPCWW9wSWo1OaM/N8dyNY8f45N/Rwa9x3XU8bGjHDi6x/eM/dqcBDg6ymEh1kyBsOiRnUY8kk8ADD/AVezTqWpgr5U6Im5/nk3K5SqJksvzwoGSSZ4DrHYauzir0bCLiiqlUivMa8Th3ZuvkeqX3EQQBaPCchYhFPVI4/wLgK3yAy101qRTvTlaTG/jc54A//EMOI0Uibt9HTw+PfjUMd3KfafKOQlc5ScmrIKyEhhYLCUPVI8UsMbLZ5cnlWqy8C0NZN97IvRm33gps2cK7iB07eAiTNvrL5dw5HOPjxTuzV/Ke5cJngiDUJSIW9UgxS4xg0LX00Ky0kkqHt86e5RLes2f5NuDap99yi1tGq3MZSrnCod1uC4WqlCCUek8RDEFoKEQs6pFiVVOtrVzRtJpKqoMHeUcQi/EOotgOQb93IOBaonv7J7xzKrRQlROEat5TEIS6R8SiHilWNfXUUzwPYjU2GcXCW4U7BP3eu3axUOhRp0QcjurpyReqZBL4yEe4ge/tt4GZmXxBqOY9BUGoeyTB3UwUS5yXS5J7q5y0N9XsrFvxBPAO4s03eSei7UG2b+edkJ5PsZL3FITNS0MnuKWDu5kYGXFzFN7y21KhrEoeTkNDvIOIRPh1dI/G+fMsHlpUVvKegiDUJRKGaibW2vFVh5h6e3lXofMbCwv541TFZVYQGh4JQwm14w1rTU+7PlKxGHd3iyAIgpeGDkPJzkKonVKWIiIUgrDpELEQakdCTILQNEiCW1gdMshIEJoC2VkIgiAIFRGxEEojnk6CIDiIWDQLKz3xi6eTIAgeRCyagVpO/OLpJAiCBxGLZqCWE794OgmC4EHEohmo5cRfzCZ9pZbogiBsGkQsmoFaTvzFbNLF00kQmhYRi2aglhO/NNwJguBBvKGaBa/duHaDlRO/IGwkDe0NJR3czYJ0WguCsAokDCUIgiBURMRCEARBqIiIhSAIglAREQtBEAShIiIWgiAIQkVELARBEISK+CoWRPTviOhXRPQaEf13P9ciCIIglMa3Pgsi+pcA7gRwnVIqTUQ9fq1FEARBKI+fO4v7ADyulEoDgFJq3Me1COWQIUiC0PT4KRa7AdxKRK8Q0Q+I6MZSBxLRJ4joVSJ69cKFCxu4REGGIAmCAKyzNxQRPQegr8hDDwP4bwAOA3gQwI0A/hLAFarCgsQbaoMZGmKBiMXc+1IpNhY8fNi/dQlC4yHeUKVQSr2v1GNEdB+AZxxx+AcisgF0A5CtQz0xNsY7Ci8yBEkQmg4/w1D/F8AQABDRbgAhABM+rkcohgxBEgQB/orFUwCuIKJ/BvAXAO6uFIISfECGIAmCAB/FQimVUUr9rlLqGqXUu5RSEgCvR2QIkiAIkHkWQjXILAxBaHrE7kMQBEGoiIiFIAiCUBERC0EQBKEiIhaCIAhCRUQsBEEQhIqIWAiCIAgVEbEQBEEQKiJiIQiCIFRExEIQBEGoyLpalK8HRHQBwKk1ftluNI6JYSOtFWis9cpa1wdZKzOhlLptnV573Wk4sVgPiOhVpdS7/V5HNTTSWoHGWq+sdX2QtW4OJAwlCIIgVETEQhAEQaiIiAXzJ34vYAU00lqBxlqvrHV9kLVuAiRnIQiCIFREdhaCIAhCRUQsBEEQhIqIWDgQ0UEi+iUR/RMR/Q0Rtfu9plIQ0QeJ6DUisomoLsv8iOg2IvoVEZ0goof8Xk85iOgpIhp35sHXNUR0ORF9n4iOOb8DD/q9plIQUYSI/oGIjjprfczvNVWCiEwi+kci+pbfa6k3RCxcvgfgGqXUdQCOA/iMz+spxz8DuAvAi34vpBhEZAL4EoBhAHsAfJiI9vi7qrL8GYBGaZbKAfiPSqmrAdwM4Pfr+LNNAxhSSl0P4J0AbiOim31eUyUeBHDM70XUIyIWDkqpv1dK5ZybPwLQ7+d6yqGUOqaU+pXf6yjDewCcUEq9oZTKAPgLAHf6vKaSKKVeBDDp9zqqQSl1Vin1U+ffs+AT22X+rqo4iplzbgadr7qtqCGifgC/CeB/+b2WekTEojgfB5D0exENzGUA3vTcfgt1ekJrZIhoJ4AbALzi70pK44R1fgZgHMD3lFJ1u1YA/wPApwDYfi+kHgn4vYCNhIieA9BX5KGHlVLPOsc8DN7q//lGrq2QatZax1CR++r2irIRIaI4gG8C+A9KqRm/11MKpZQF4J1ODvBviOgapVTd5YaI6AMAxpVSPyGi/X6vpx5pKrFQSr2v3ONEdDeADwA4oHxuQKm01jrnLQCXe273Azjj01o2HUQUBAvFnyulnvF7PdWglJoiohfAuaG6EwsA7wVwBxHdDiACoJWIvqaU+l2f11U3SBjKgYhuA/BpAHcopeb9Xk+D82MAg0Q0QEQhAL8D4G99XtOmgIgIwJ8COKaU+qLf6ykHEW3RVYVE1ALgfQB+6e+qiqOU+oxSql8ptRP8+3pYhCIfEQuXJwEkAHyPiH5GRF/xe0GlIKJ/TURvAbgFwLeJ6Lt+r8mLUyjwAIDvghOw31BKvebvqkpDRF8H8DKAK4noLSL6Pb/XVIb3Avi3AIac39OfOVfD9chWAN8non8CX0B8TyklJakNith9CIIgCBWRnYUgCIJQERELQRAEoSIiFoIgCEJFRCwEQRCEiohYCIIgCBURsRAEQRAqImIhCEUgor8joimxqhYERsRCEIpzENz8JggCRCyEJoKIbnSGW0WIKOYM5Lmm2LFKqecBzG7wEgWhbmkqI0GhuVFK/ZiI/hbAfwXQAuBr9eiAKgj1iIiF0Gx8DuxTtAjg3/u8FkFoGCQMJTQbnQDiYNPIiM9rEYSGQcRCaDb+BMB/AQ+3+oLPaxGEhkHCUELTQEQfA5BTSj1NRCaAI0Q0pJQ6XOTYHwK4CkDcsYP/PaVUXVnBC8JGIhblgiAIQkUkDCUIgiBURMJQQtNCRNcC+D8Fd6eVUjf5sR5BqGckDCUIgiBURMJQgiAIQkVELARBEISKiFgIgiAIFRGxEARBECry/wF1WcETcTdvswAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 402.375x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "centers_true = numpy.array([[-0.5, -1],\n",
    "                           [1.0, 0.5]])\n",
    "covars_true = numpy.array([[[0.125, 0.0],\n",
    "                            [0.0, 1.0]],\n",
    "                           [[0.5, 0.0],\n",
    "                            [0.0, 1.0]]])\n",
    "\n",
    "points, labels = gen_mixture_of_gaussians(centers_true, covars_true, 100000)\n",
    "n, d = points.shape\n",
    "classes = numpy.unique(labels)\n",
    "k = len(classes)\n",
    "\n",
    "print(\"There are {} {}-dimensional points.\".format(n, d))\n",
    "print(\"There are {} classes: {}\".format(k, classes))\n",
    "print(\"\\nCoordinates of the first few points:\\n{}\\n...\".format(points[:5, :]))\n",
    "print(\"\\nClass labels of those first few points: {} ...\".format(labels[:5]))\n",
    "\n",
    "visualize_clusters(points, labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-cb3579f34a2e2239",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Background: The $k$-means clustering criterion ##\n",
    "\n",
    "Recall the $k$-means problem. You are giving a $m$ data points in a $d$-dimensional data points, denoted by the $m \\times d$ matrix $X$. That is, let each data point be a $d$-dimensional column vector $x_i$ organized as follows:\n",
    "\n",
    "$$\n",
    "  X\n",
    "  \\equiv \\left(\\begin{array}{c} \\hat{x}_0^T \\\\ \\vdots \\\\ \\hat{x}_{m}^T \\end{array}\\right)\n",
    "  = \\left(\\begin{array}{ccc} x_0 & \\cdots & x_{d-1} \\end{array}\\right).\n",
    "$$\n",
    "\n",
    "In the $k$-means problem, we hypothesize that these data may be divided into $k$ clusters.\n",
    "\n",
    "For each cluster $C$, consider its center $\\mu$ and measure the distance $\\|x-\\mu\\|$ of each observation $x \\in C$ to the center. Add these up for all points in the cluster; call this sum is the _within-cluster sum-of-squares (WCSS)_. Then, set as our goal to choose clusters that minimize the total WCSS over _all_ clusters.\n",
    "\n",
    "More formally, given a clustering $C = \\{C_0, C_1, \\ldots, C_{k-1}\\}$, let\n",
    "\n",
    "$$\n",
    "  \\mathrm{WCSS}(C) \\equiv \\sum_{i=0}^{k-1} \\sum_{x\\in C_i} \\|x - \\mu_i\\|^2,\n",
    "$$\n",
    "\n",
    "where $\\mu_i$ is the center of $C_i$. This center may be computed simply as the mean of all points in $C_i$, i.e.,\n",
    "\n",
    "$$\n",
    "  \\mu_i \\equiv \\dfrac{1}{|C_i|} \\sum_{x \\in C_i} x.\n",
    "$$\n",
    "\n",
    "Then, our objective is to find the \"best\" clustering, $C_*$, which is the one that has a minimum WCSS.\n",
    "\n",
    "$$\n",
    "  C_* = \\arg\\min_C \\mathrm{WCSS}(C).\n",
    "$$\n",
    "\n",
    "**Lloyd's algorithm.** The standard $k$-means algorithm (also known as _Lloyd's algorithm_) is a heuristic iterative method for finding a local minimum for this optimization problem. (Finding the global optimum is [NP-hard](https://en.wikipedia.org/wiki/NP-hardness).) The procedure alternates between two operations: _assignment_ and _update_.\n",
    "\n",
    "**Step 1: Assignment.** Given a fixed set of $k$ centers, assign each point to the nearest center:\n",
    "\n",
    "$$\n",
    "  C_i = \\{\\hat{x}: \\| \\hat{x} - \\mu_i \\| \\le \\| \\hat{x} - \\mu_j \\|, 1 \\le j \\le k \\}.\n",
    "$$\n",
    "\n",
    "**Step 2: Update.** Recompute the $k$ centers (\"centroids\") by averaging all the data points belonging to each cluster, i.e., taking their mean:\n",
    "\n",
    "$$\n",
    "  \\mu_i = \\dfrac{1}{|C_i|} \\sum_{\\hat{x} \\in C_i} \\hat{x}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-13840d179f2b5013",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Initial code implementation ##\n",
    "\n",
    "The following cells implement Lloyd's method. The data matrix $X$ is stored in the 2-D Numpy array, `points`, and the labels as a 1-D Numpy array, `labels`. Note that the k-means algorithm you will implement should **not** reference `labels` -- that's the solution this algorithm will try to predict given only the point coordinates (`points`) and target number of clusters (`k`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-7dee06f2303f7153",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**Step 0: Initialization.** To start the algorithm, you need an initial guess. Let's randomly choose $k$ observations from the data. The function, `init_centers(X, k)`, does so by randomly selecting $k$ of the given observations to serve as centers. It returns a Numpy array of size `k`-by-`d`, where `d` is the number of columns of `X`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "init_centers",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def init_centers(X, k, init_seed=1032742):\n",
    "    \"\"\"\n",
    "    Randomly samples k observations from X as centers.\n",
    "    Returns these centers as a (k x d) numpy array.\n",
    "    \"\"\"\n",
    "    from numpy.random import choice, seed\n",
    "    seed(init_seed)\n",
    "    samples = choice(len(X), size=k, replace=False)\n",
    "    return X[samples, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 1: Computing squared-distances.** Given the $m \\times d$ data `X` and a set of $k$ candidate centers `centers` (of size $k \\times d$), the function `compute_d2(X, centers)` returns a $m \\times k$ matrix $S$ whose $s_{ij}$ entry stores the _squared_ distance between each point $x_i$ and center $\\mu_j$.\n",
    "\n",
    "> The cell below first defines a function, `compute_d2__0()`, and then creates an alias `compute_d2`. You'll see why when you get to the exercise at the very bottom of this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "compute_d2",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def compute_d2__0(X, centers):\n",
    "    return numpy.linalg.norm(X[:, numpy.newaxis, :] - centers, ord=2, axis=2) ** 2\n",
    "\n",
    "compute_d2 = compute_d2__0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 2: Re-assignment.** Given the $m \\times k$ squared-distances as a 2-D Numpy array `S`, the function `assign_cluster_labels(S)` assigns each point to its closest center. The return is a 1-D array, `y`, whose `y[i]` entry indicates that point `X[i, :]` belongs to the cluster whose label is `y[i]`. By convention, the labels are integer values in the range of $[0, k-1]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "assign_cluster_labels",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def assign_cluster_labels(S):\n",
    "    return numpy.argmin(S, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 3: Update centers.** Given the $m \\times d$ data set `X` and $m$-vector of cluster assignments `y`, the function `update_centers__0(k, X, y)` determines the centroid $\\mu_i$ of each cluster $0 \\leq i < k$. The centroids are returned as an $k \\times d$ array.\n",
    "\n",
    "The code cell below further defines a variable, `update_centers`, which serves initially as an alias for `update_centers__0`. The reason is so that we can later redefine the alias to an Intrepydd version and easily compare different implementations by simply redefining `update_centers`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "update_centers",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def update_centers__0(k, X, y):\n",
    "    # X[:m, :d] == m points, each of dimension d\n",
    "    # y[:m] == cluster labels (max of k distinct labels, 0 through k-1)\n",
    "    m, d = X.shape\n",
    "    assert m == len(y), \"Number of points ({}) and labels ({}) don't match.\".format(m, len(y))\n",
    "    assert (min(y) >= 0), \"Labels must be positive.\"\n",
    "    assert (max(y) < k), \"Labels must be between 0 and {}, inclusive.\".format(0, k-1)\n",
    "    centers = numpy.zeros((k, d))\n",
    "    for j in range(k):\n",
    "        centers[j, :d] = numpy.mean(X[y == j, :], axis=0)\n",
    "    return centers\n",
    "\n",
    "update_centers = update_centers__0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 4: Calculating the within-cluster-sum-of-squares metric.** To verify the progress of the method, the following function will help compute the within-cluster-sum-of-squared distances, given the squared-distances matrix `S` as defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "WCSS",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def WCSS(S):\n",
    "    return numpy.sum(numpy.amin(S, axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 5: Convergence testing.** Lastly, to determine whether we can terminate the method, here is a simple function that tests whether two sets of center values are equal or not. Note that this implementation uses Python sets so that it does not depend on the ordering of the centroids."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-84eaf8eb026c9273",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "def has_converged(old_centers, centers):\n",
    "    return set([tuple(x) for x in old_centers]) == set([tuple(x) for x in centers])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Putting it all together.** From the six building blocks above, here is the entire implementation of Lloyd's algorithm. It accepts the dataset `X`, a target number of clusters `k`, and two optional arguments: the initial centers (`starting_centers`) and the maximum number of iteration steps (`max_steps`). Following the definition of this function is some code that executes the algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-7d55bc2221f76601",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def kmeans(X, k,\n",
    "           starting_centers=None,\n",
    "           max_steps=numpy.inf,\n",
    "           verbose=False):\n",
    "    if starting_centers is None:\n",
    "        centers = init_centers(X, k).copy()\n",
    "    else:\n",
    "        centers = starting_centers.copy()\n",
    "    old_centers = numpy.empty(centers.shape)\n",
    "    \n",
    "    converged = False\n",
    "    i = 1\n",
    "    while (not converged) and (i <= max_steps):\n",
    "        old_centers[:, :] = centers\n",
    "        S = compute_d2(X, centers)\n",
    "        clustering = assign_cluster_labels(S)\n",
    "        centers = update_centers(k, X, clustering)\n",
    "        converged = has_converged(old_centers, centers)\n",
    "        if verbose: print (\"iteration\", i, \"WCSS = \", WCSS (S))\n",
    "        i += 1\n",
    "    return clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iteration 1 WCSS =  565260.5748947266\n",
      "iteration 2 WCSS =  348579.85189060034\n",
      "iteration 3 WCSS =  233326.26372745118\n",
      "iteration 4 WCSS =  226161.20338560524\n",
      "iteration 5 WCSS =  225616.41204630243\n",
      "iteration 6 WCSS =  225564.38665114215\n",
      "iteration 7 WCSS =  225559.57379539014\n",
      "iteration 8 WCSS =  225559.0815381379\n",
      "iteration 9 WCSS =  225559.0249336476\n",
      "iteration 10 WCSS =  225559.02059976387\n",
      "iteration 11 WCSS =  225559.0201870468\n",
      "iteration 12 WCSS =  225559.01980372274\n",
      "iteration 13 WCSS =  225559.01936284977\n",
      "iteration 14 WCSS =  225559.01931455365\n",
      "iteration 15 WCSS =  225559.01927797496\n"
     ]
    }
   ],
   "source": [
    "update_centers = update_centers__0\n",
    "clustering = kmeans(points, k, starting_centers=points[[0, 187], :], max_steps=50, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Quick visual inspection.** Just to make sure the algorithm computed something \"reasonable,\" let's create a quick visualization of the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "kmeans_test",
     "locked": true,
     "points": 3,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "=== Estimated centers ===\n",
      "[[-0.37903808 -1.1583286 ]\n",
      " [ 0.93029433  0.72554208]]\n",
      "\n",
      "=== True centers ===\n",
      "[[-0.5 -1. ]\n",
      " [ 1.   0.5]]\n",
      "\n",
      "176215 matches out of 200000 possible (~ 88.1%)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAFgCAYAAABKY1XKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvX10XOd93/l57r1zZ4CZwYAgCZAiRRGyoET0i7quXce0YimUuzHcrHlWdr2xmh4mOa1PonDruBs2Tncb1XIdq8s2rlrZ8ea0TrCNbMeJdVZnayPd2iopqbRz7KSmZVOJQAngiwgSAIEZzAwwc+fe++wfz73zAsxg8DLAYIDncw4OMDN3Lh4MZp7f/b19f0JKiUaj0Wg0K2G0ewEajUaj2f5oY6HRaDSapmhjodFoNJqmaGOh0Wg0mqZoY6HRaDSapmhjodFoNJqmaGOh0Wg0mqZoY6HRaDSapmhjodFoNJqmWO1ewFp5//vfL//sz/6s3cvQaDSatSLavYCN0HGexczMTLuXoNFoNLuOjjMWGo1Go9l6tLHQaDQaTVO0sdBoNBpNU7Sx0Gg0Gk1TtLHQaDQaTVO0sdBoNBpNU7Sx0Gg0Gk1TtLHQaDQaTVO0sdBoNBpNU7Sx0Gg0Gk1TtLHQaDQaTVM6TkhQszMZGxvlwoWzpNPj9PYOcvz4GYaGhtu9LI1GE6CNhabtjI2NMjp6GsOwicX6yGYnGR09DTxd12Bow6LRbD06DKVpOxcunMUwbGw7jhAC245jGDYXLpxddmxoWLLZyRrDMjY22oaVazS7B20sNG0nnR4nEumuuS8S6Sadnlh27FoMi0ajaR3aWGjaTm/vIKXSQs19pdICvb1Hlx27FsOi0WhahzYWmrZz/PgZfN/BcfJIKXGcPL7vcPz4mWXHLjUshUKGmZm/IpudZGTkhA5HaTSbhJBStnsNa+Id73iH/P73v9/uZWjWwUqJ6cpjE/T2Hm2YtK5Ohnteifn5q4AklboLw4jg+w7Dw/UT45v5dwA66a5pRkePVdXGQrMlVG/ykUg3pdLCqjb2lTbm69e/Axj09NxBNJoCwHHyJJMHOXXq+S37OwqFNCCJxfas6W/T7Do62ljoMJRmS1hPYrpR5RPAqVPPk0gcYN++nywbCtj8/EW9v6NYnKdYzOqku2ZHo42FZktYT2K6mYFZS2K8VdT7O3y/hO+7NffppLtmp6GNhWZLWM/G3szArCUx3irq/R2GEcEwavtbN9toaTRbjTYWmk1lbGyUkZETTE39mExmgnz+1qo39mYGZmhomOHhp0kmD1IozJFMHtz0PEE9AxWN9hCNJldttMLX5KmnBnUFl6Zj0AluzaaxNBmcy01RKMwQjabYv/9Y04qh9SbFN4sw2T49fQnPK2IYUfr7jy2phlp9Ndd2+Js0W0pHJ7i1sdCsivXoMY2MnCCbncS24+X71lqttNqS2s2mVZt8vdckm72F6+aIxXp12e3OpqONhRYS1DRlrUJ/Ien0OLFYX819a038Dg0Nb9rGuRYDWJ1sB7DtOI6j7l/L+pa+JoVChsXFKaT06e29e9WvrUaz1eichaYp69VjWk1Su13x+7UKErZKZmTpa6JyOGBZXbrsVrOt0cZC05T1bpTNqpXaqSC7VgPYqjLdpa+J6xYASSIxUD5Gl91qtiM6DKVpSm/v4LI4+2o2ShVGebphzqFVoZ31UC9E5vslrl37Dv/qXx0IEtg2/f1v5vjxMxw/fobR0dM4DjU5i6UVT81CW0tfE9uOY9uJmsZCXXar2Y7oBLemKWtJ7q4lD/DUU4PEYn0IUcn7SSnJZt9g//5ja0qmrzUBX51oLhQyZLNv4LqLqBykwDAshICurn4iEZvh4aeBlSue1pME19VRu4qOTnBrY7HDadVUudVUJa3VqDz77KM4Tg7L6sK2kzhONgj1SJLJO+ju7t+0DXdsbJTnnvtlFhdn8X1n2eOGYSOEwDRtkslDNRVcjV7TkZETzMyM4TgZPK+IaUax7RT79g2tWP21XSq+NJtORxuLbRGGEkKYwPeBN6SUP9fu9ewU1lvFVI/VVCWtNqwUrsuyEpRKC5RKCzhOFiFMpPQBg3x+CsuKEY2mmoam1h/Okkjp1rlf4PslhBB4XpFs9jqLi3M1a6/3mk5N/ZhiMQ0IhDDxPIfFxSmmpkorvm6bWfEVokfRajbKdklwfxx4pd2L2Gls9VS51SbCw3UlkwP09BypeUwIA8uyAYNc7hZQySU0qpha+nuLxQzz89e5cuWFhhVWFy6cJRbbg2FYQSWSQeXCT6IMiQ8IXNfBcTLlDbfRa+r7DlKCYZgIITAMEynB94vLfv9WVoHpUbSaVtB2z0IIcRj4O8BngH/c5uXsKFrR57AWensHG4ZhGq0rFkthGCZgAx6mGcXznODKvEixmCGTubKid1SdgFfHX0VKiES6GnpT4RrC32cYETxvaTiqkrvo7t7Ht771SWZnX8X3fSwrSiSSoFTK4bpFMpkJLCsGSHzfRwgjMDYS04zWnHW1Hl+rvIF2FhJodg5tNxbAvwH+CZBsdIAQ4mPAxwCOHDnS6DDNEtZbxVSP1WxcR48+xNWrLwIGQpi4bhHXvcnRo/9wxXWpYwvBowZSekgpMU2b+fkbgCCZvKN8Jb90o6uuVMpmbyIlCAHx+IGGG2O4hnh8IBigJIIwmIfyLCivLZU6jO9LZmYuYRg2oEpeVejMQkVRDRwnD4DnKU/CNKMkEgP09dUay3qbdz6f40/+5O8GfzskEgdw3QKxWO+GQ4hbfdGg2Zm0NQwlhPg5YEpK+RcrHSel/H0p5TuklO/Yv3//Fq2u82mFKuvY2Chf/OL9fPWrH+T69T8HzIZhjImJcyQSB7CsKOBhWVESiQNMTJxbtq5CIc309CvcvPmDoApJAibhlTnIYDP16ek5QizWeGZFtaCg6y5iWXbNcyKRbqamLtWEfY4efQjfdzAMi2TyzmDD9wMFWRvDiGKaNuAjJeRyk4EnkML3HXxf5SGkdPH9YiBT7uH7JQzDxjRtpPQoFnMsLNyuCTdVh80KhQzT05eYn79OqZTHdR08rxiILt6kWJzfcAhxq6XctVDizqTdOYv3AB8UQkwAXwVOCCH+qL1L2jlsVJU1DJfcvn0ZUFfd2ew1fN+t2bjCzeHKlRcoFNLE4wP097+VvXvvpbu7f9kV7I0b32Nx8TauuxhcyQMYWFYEkFhWjP377+M3fuMmhw+/G9OM1Dy/3kY3NDTMqVPPc9dd7yWZPFxjXHK5KRwnUxOzv3hxhPvvP0UyeTCoyLIRwsQwLLq79yKEJCyjzWZv4PsO0WgPi4u3l8mRC2Hh+x7gI4QZ5FwkQhiUSnk8z6nxDqLRHkqlBbLZG8zNvRYYyxBlKCtrv0GxmAHW7w1spZS7zo/sXNoahpJS/hbwWwBCiIeA35BS/kI717TT2Eilzbe+9Umy2Zu47kKQALYQQpDP36Kvb4h0eqIm/m5ZXbiuE4R11JV9PXmPl176bBDPr8UwLPbuPYaUkkJBVR+tthkupN7xhcIM3d37lsXsJybOlY83jD04zmtI6bG4eJuurr04Ti4Ijwn27z8WGE2BaVplz0JRu8Hv3XsvADMzf43nFZf9XimhUEiTz9+qeV4jZmcvY9tJYrFUTUhrbbPAGzdHthKdH9m5bIechWYbMjY2yszMpSA8I5DSR8oShmHhecWyEQg3ByldpPTKPQtzc5cDAwOWZTM2NsrQ0DAXLpzF80p1jIVfvsIOzx1uho6Tx/NmayTBG2084cb4rW99kpmZS0gJUnoYRm2SObxKr97cLCuG5ykj4Dg54vEBstkb5bV6XgEhIixtTaruVapuMPS8QpD0rv29hcIc8Xh/YCxWh+sWyeVu8va3q/xPvST5c8/9MuEs8Hp5jo1u1kuN09GjDzExca7GMOn8yM6l3WGoMlLKc7rHYvsQbqLqSlp9V/kENwi7OOXNwfdLZDJXl40WDUtPFxdny6GIdHq8vLmqjbWyufq+Xw6RHD36UDmcYdsJXNehULjNwsJtoHlcvFTKkUodZf/+YxhGhPn5qxQKmarHF7DtJNeufYe5ude4fftVbDuJlC6e5+A4WebmXsPzivT0HMJ1nSCZrfIUQphlQwoqma3CWHY51GMYFrbdu2RdyhA6zjy2Hcc0o8uqpephWVG6uyv5n62eBb40vDQzM8YLL3ya2dmxGsNk2z1bPupWszVsG2Oh2V6k0+MkkwdRIZZag9HXd0859xFWFYVVTNWbvxAGhhGhUMiUN63e3sHy1bcyGtXHU86rTEycK3ss8/PXkFLlA2ZnL/Pcc7/Ec8/9csO4+Le/rcJn6fTrwWbWC8jAS1AbeaEwRz4/hRBGkFsokMtNBgYv9BRUr0WptIhtx+nu7kcIQW/v3fT2Hg16Kazy7e7ufSQSA0xPX2Ju7jV8X5LL3WB6+hUKhXRNrqC3dzAwJM3DUJbVxd6995JIVPI/q5kFrmRMrnPlyosbTjQvNU6OkwEMCoVMjWESgi0fdavZGnQYSlOX3t5BZmfHqspa1cZuGBG6u/eWjwtDEfU2PSlVktfziuVQxAc+8HmuX/9OVShKGaOurj4GBt5WlsX45jcfIxbrY3Z2DDAwDAMpVcf1wsJtpHTJ52+iDFkE3/f46lc/SCSiei2EsDBNm1KpgOPkEMLAdQvMz79Bf/8xLMtmcXEWKf1yqWsjcrlJAIrFeVy3SD5/E9OMsn//MQqFeXK5NwKZ8Sj5/NSy18J1F5mbG2dg4C08/PDvlsNBKkTUj+Okgx6PauNZOUcyeSjY+FVIbGTkRDlJXl0WbRiVQoBCIcP8/NVgXbENz8lYGl7yvGK5FyYkDLF94AOf1/IlOxCtDaWpy/nzT/DCC59GeQwSKVUsPx4/WC4fvf/+U1y8OFI3BFWNqhDqxvMWAwNi1fQzSCmxLJsHHvgtHnzwt4GK0F86/Xo5/KMqjmSd5rmGvxmW9ExUqq/C+yLlv2015wgN19K/1zCsFV4DEXSmR4nH++skoyeIRpPkcrcQIoLjpCmVFpHSJxrtwzQjLCzcrPpdSgPLsmLEYr1Vifw0Yc4im70ehM4glToSyKasbUphNUsn/N2+/SquW8SyouWE/kbOv0voaG0oHYba5TSK/U9MnKO7+0AgWRFupoJiMV0OOXz3u5+r2rhr8w/VqFBOFt/3yk1tS0NWnlfipZc+y/nzTwCVcs+wLLViKFbWWVrym5esw1t+xIqGYvk5wK9rFFYyluHfWyot1ITNbtz4Xvnxrq693H33+3DdHKVSgWg0yVvf+iip1KEaQ6F0qxyKxXksK1ZTFn3y5Jc4efIPSCYPUioVME2bVEo1sd6+/Spzc69x7dp31hWOWlp+a9spwCcWS+lw0y5BexY7mGZd1yuptX7zm48hhMn8/LXAIFTCI31992DbPdy69UMsK1aO8/t+bZWTSgKL8kaqrq5VdVBtEtTANCNIKTEMg/vu+xBjY98ox8OV59FFInGQdPq1zX3RNplQi0qFyPKkUkeJRLpZWJgil7tJd7dqalT5EwfTjFIq5Rue76GHPlX2xqoJPQEp3cDzq4gmWla0xotbLUvVcSvVUDrctEo62rPQxmKHshrZ7qWhBaiEEgCuX/9zpPQDgbzwfaKSmcnkIfL5mzhOrpzYBBUqUlpLFr29dzE//8aSpjOBcmi9qttQidfXvh9VxRH09NyJ6y6uqdx0O6IS6ma5T8O2kyQSA+Ryt3DdYlBurPI4jUQIqzEMi/e+958t2/gr1Us3cV0HcIPj7eC7wUc+8uwG5eq1gu0a6WhjocNQO5SlFUFSusvKKFdSiQ3DDqoTWclwgMQwLFx3Ed93+Kmf+kRV38UipdJCWQRQCIN0emKJoSA4T3U4qNpA1EuSe0jpBfIXnW0ogMD4VkJfjpNldvYyjpPF9x08r1AOtSmBxWYYfPe7n1sWTgQYHn46+F1hiEwE51X5lfWU1J4//wRf+9ojXL36Ivn8DLOzY7pDe5egjcUOZGxslOnpS/i+hxAWnhf2QZRqmqNW0gwaGhpm375j5VCSacawrC6UZlOC4eGnefDB3+a++z5UVldVbycj2BDDPEMzlndyaySeV2xapRX2dhQK83UlNm7c+N4yg+P7JTyvhGnG1twoF3bf+76PYdj4vks+P4XnOZsme6/ZPmhjsQMJa+KFUL0LhqE28Wx2sqY5qplm0Pve9yTJ5AF6e+9m//5jpFJ3kUrdySOPfLkcdshmb9DTcwjbTmKaFrbdTU/PHTWdzPUJezc0jVhaBLAUZVAKwf94+YyNF174HVy32uCEczpcLCtGoTC3JrE/NbPDC+Z1VN5Xi4tp3aG9C9B9FjuQdHqcROJgIPrnBXFwuaxapZlmUPXj09OXcJwsnlfij//4Efbtu5eHH36SdHqc7u5+4vGB8nmllOTzM5imscLVcVjZtDxPoamm2WsjkZJgs/bLsujFYm6FfIcMJEcOrkn+PJ0ex7KieF6pnEtSfTQF3aG9C9DGYgcSdlX39Bwhn79VnmsthFkOF1QbhJU2iPCx5577pWCTUH0X09OX+PrXHw2kPq6Vj49EuonFUiSTB8lkrjZZqTYSrUGFrSyrqxxybF5iLEgmlYFfrdhfb+8gruuyuDhVvghRnoalS2Z3AboaagdSXQnleaVABVaSSt2FYURqqqLGxkbLE+CkpOwxVG8aIyMngsoorxwDd91SkDQ3G/YYCGFhGMYamug09Vmd9yWEiWnaQb7IqduEWE0icZBk8g5AeYNhd3uzUutSycFx0rhuEcMw11WGu0vp6GoobSw6kNWULobHXL/+HcCgp+cOolE14yEsjz169CFeeOFfBBUzAjAxDOjq2sfJk18qn/PJJ3vLMxVAGQE1FMin+UZmBLFtsUTWW7M1rPz/2bPnHmKxFNnsLQqFmXLfR71Sa1jea1HPoOiy2oZoY7GV7HZjsZr+ieoPbDZ7k56eQ4GYniK8inScTJAArWgSKU2lCIcPv4tTp57n2Wd/gZdffqYNf6lmKzDNKL29g2QyE3R37yMeH6BYzAR9H4vYdqKmoCGk0SyNZu/NXU5HGwuds2gza/ES0ulxCoU0tp0gFtsDLI83L51zkMtNkU5PlMd8mmY00HYqVoWPKlefUpZw3RJTU5cYGxvlRz/6yha+GpqtRuU6bKLRFN3d/RSLmSDXZAQ6Vfllye96szRGR08TiST04KMdjC6dbSOrGUG59BjHyZPL3SqHhQqFDJnMBBMT/5XPfKaLP/mTD1MqOeUyymg0hZReoBxrlofoAJhmOJhnuXe5sDDFN77xWFWoSbMzEThOjkRigFJpgVzuFpXQoY9lxZY1c9abpWEYNrOzrzZs8tR0PtqzaCOrGUG59BjLigUb/i2khExmIvAQlIaS6y7gukUikVgw1jRXpfLqYVlRbDuF6+bo6updQXdIkslMlH/W7FQk2exNbDuB48yXx8hKqUboxuMDyzb8qakf47qqW980o8TjA0SjPUjJMtl0Pfho56CNRRtZzQjKpceoD3UOx1HT3ELUlDYTz1MjUPP5W8RiqfLcgUikqywlrXIWJTyvetDPRtC9Ep2M5xVZWFggGt0bzBnx8H2XROIAsZiSNg83/LGxURxnPujiNvE8NXO9q6ufffvuDd6bq5uXrukstLFoI2E/xEpXYtXHFIsZFhdv1z2XlB6e55ZLY123UDU7wsG2E9y+/WpgPCzi8X5cdzHom6iW3Kje+MOf6xmDUPBu6XM0nUZYXut5C+Vkt+975HKTFAoZLCuKado89dQghUKaSER5Ier9pUqnC4UZPvjB3wdYsVpK07loY7EOWlUeePz4GZ577pfKw4MMwyIaTfKzP/u7NceMjp7GcSCbvUlYvGYYkUBkr7LRq76HCKYZJRLpolCYo6/vHjKZqywszADh3AWHTOYKXV19xOMHyOVuVK1KLvm5kSHwG/ys6VRcd5FCYbamb8b3SxQKBSwrTiLRz/z8GwB0d+8POvqLmKYKk9Z2/mt2GtpYrJFGlSDrHVe5fIxmbTK5WnJjdnaMSKQL1/UxDAspq+U0ZDDRzqerq4+TJ/+gvJ4vfvF+Zmb+OjjWwDRtPK/I4uJtDMMKvJFGPRDaY9g9CBYXZ6m8B9UsEiEMHCeDEAPlnJnjZJdNyKuH7rvYOWhjsUZWk5Rey7lisV56eg6V73Oc/LJzhZIc4fyJbPaNYGaEiZQVD0MIg3377uN976vtwC4W54P5zCKocgEvaOwNNwOlIlqi1jjo8NLuoloqXgRS5qWaOeXx+ADz81dxXTUid6W8RKMLqxs3TgVDk7QB6SS0sVgjq0lKb+Rcvl/i2rXv8NRTg+UKE8eZp7d3kKNHH+LixRFsO8Xi4lSQoIZk8g5M027Y/NTbO8j8/PXy4JvweSFS+khZT5JDG4rdi/rfq0mFSnJ+auplTDOKbSeR0qNQmFsxL1HvwiqbzfHSS58llTraIs9cs1VoY7FGVpOUXu+5VEPUlUBe3GR6+hIg6Ok5QjY7ycWLIxw+/G7Gxr4RJLMNbDtBX99Qww/ss8/+AleunA8+8C4qMb0anaGliW/NbkMNv1KoXJnq03HdxbrT+ZZS72LIcdL4vqsb9zoQ3ZS3RprNgNjIuebnb6DUQO8gn59CCAshDBYWprDtOKWSw49//Ce4rgpBGYaNacZWNBQvv/zMkk0/HFRUzfKmO20oNNWo8KXq04lEUnz3u59rOguj3nAt1y1WNYMqdONeZ6A9izXSbAbERs4FPj09R4jFUszPXw0E+2Q5XlwozARqokZQDeWzuDjDt7/9SYByIjEMX01N/XCVK9HhJs3K+L7Lnj1vAmB+/ipS+vT23r0sjFSd0I5GeygU5oBK34VhmHR19dacWzfudQZaSHAbESawbTse9ESohLNp2uzdey+Tk38JQCTSVX6O53nBB/fOcpI6k7mCVnnVtBolea48TsuKsX//MaBSDRWWeVcLCRYKaeLxfhwnS2/v0XLebZeKDXa0bo72LLYR1T0V8Xh/edPv7u7HcZQsh2HU/suEIJgzYZeNTOiRaDStpHo2hu+XKBYzRKOpchipXkIboLt7L7/6qxfLz73jjnfqxr0ORHsW24zqeQHRaDKohlJXZXNz44EiaG15qxACy+rC952qCWmd9X/VdAIiGNGrJGNABgrIKfr6hsoJ7er561KqEa4f//jrDc+6i3oxtGehaR0rjTk9f/4Jzp9/ovxBrfRBqGl06kOqjYRmc7CsGKYZo1icK98Xqhi//e3/kIkJ1lwp2PomV81moY1Fm1jP1dTExDmSyYMUCplAZiGK5zmBRg9NxpdqQ6LZGFL6NYYCwDBMTDPJd7/7OQzDplCYK89prydfs5RWNrlqNhdtLNrAWq6m6k29SyQGyOVu4XnFsoEQwmRlY6ANhWZj+H7tPG8hVIe356URQtDTcxeLi7fx/VLwfrRoFnlpZZOrZnPRxqINrHQ1FT6eTo9j2z3k81PEYr01U++EEEESu/JBVJUqXr1fp9G0BCkrnf+qx8eiVFoEIBJJsLAwhXpPGoEKcpF8fopnn3207mhWqN+YqvqNfEZGTuzk/EXHoZvy2kA6PV53otj09KWaqXizs5dZXJxBShchVLNeOGvA8wpVIoK13bYazWbj+25gKMJEdxLHyQYGRZXXqosXn0Ihzde+9gjnzz9Rfv7Y2CgjIyeYnr5EJjNBNnuLQiFNOj2B75dIJA7VnRypaR/as2gDjSRDPK+IYewp368+eCa53C2i0VSbVqvR1KO2wz+Xm1z5aN/npZc+yx13vBOgHIZNJg9hmhEWFmbI5z1M0yaZvINYTL3fdf5i+6CNRRuo7qfw/RLZ7GTgGQhisUojXZjADj2IfP5Wm1as0TSjeU7M971yqHXpqGAVuipgGNGa5+j8xfZBh6HawNDQMMPDT2NZNpnMFTzPwfc9fL9EOv0as7OqJj0SSZST2JOTf1kzRlWj6SR838H3S0xMnOPq1ZfK6gJKPPNqkDwXuK4a01ooZAAtBbKd0J7FFtCoTPbChbNEIimKxdpRqcXiHJOTP2D51DqNpvPxfZe5uXEsK4rrFlCzM0wsKxbMZoF8/mYwT0PP8N4uaGOxyaxUJptOj+M46QbPVJVNQphYlnLNS6UCeoSpppNRBkAlvpWhUAly3/fp7t4XzLy4SalUKOtN6XzF9kAbi01mpTLZ3t7BpvFYKf1lMs8aTaeiKqSWTmNUPRuOkyOZvAMhLJLJg5w69Xxb1qipj85ZbDKNymTT6QmOHz9T1tppjA4/aXYOSrVWSdUIYZTf/0KA6xY2NB9Gs7m01VgIIe4UQvxXIcQrQogfCyE+3s71bAb1BsDkclMUCnN885uP0d29v00r02jaiUQIA9OMljvBfd8lnZ7AthPtXpymDu32LFzgf5NS3gf8FPBrQohjbV5TS1k6DS+bvcXCws1ArbOPWGwPsdie8nxsjWanE44NVl5FKIqpkNLj9u3LPPfcL+tmvG1GW42FlHJSSvmXwc9Z4BXgUDvX1GrCMlklADiH6+ZIJA4Qjw8ghMC248TjAxw58h66u/vbvVyNZhMRWFYMIQy6u/djmna5yzuk3gRIzfag3Z5FGSHEUeB/AP68vStpLWHZ7PT0JQqFOYrFeRYX0xSLmfIxYQ5jaW5Do9lZSGy7FyEMFhdnSCQG6O9/K+FMHaV5Fs77NpmZebW9y9XUsC2qoYQQCeDrwK9LKefrPP4x4GMAR44c2eLVrZ+wbNbzHBYXZwGBlJJSKcfs7GUsq4uenkMsLqYpFGbLIys1mp3K4uI08Xg/CwszzM/fYN++nvJjhhEp/xxOgBwZObEbhiJ1BG2flCeEiAD/CfjPUsrGwvcBnTQpL5ypnc2+UZ6nXS3+p8oGBbp3QrM7USODS6UsrlvEMKxgFoaP77sIYbBnz5t20qzujp6U1+5qKAH8B+CV1RiKTiOdHsf3SzhObplKrEKiDYVm9yLxfYf3vOc36e7ehxAmvu+Vk9+JxAC+7zI7O0Y6/TrZ7E2+9S2dx2gX7Q5DvQf4+8DLQogfBPf9UynlN9u4pjWx0sS7aLSH6elLbV6hRrN9KRTmeOWVr3Py5B+UZ8/39h5levqcAPXJAAAgAElEQVQShhFlfv4qqidDzWuZmbnE2NjomryLXTTje1Npq7GQUr5EB7tmzSbeSUnVvGyNRlOPW7d+yIULZ2s28ZGRE1y//ueEulEAUgoMw16TZLme8d062u1ZdDSNpDy+/e1PBhVQP9bT6zSaVTA7O8bXvvYItt1Df/+bOXr0Ia5efREwg4suH/BJJg+tSbJcz/huHdumdLYTqSfl4fslpqcvkc1OEol0tWllGs32wDAiwQjgxpimTT4/he/7uO4C2ewkFy+O0NNzJ4ZhIqWLaUZIpY5gGJE1SZavJLejWRvaWGyAelIe2exk+UomHj/QppVpNO1HCKssO9442mwEY4JL+H4J1y1g23EMww5CuD5S+kgpcd3CmnWj6n1G9YyM9aGNxQZYKuURiqAlEgcBiMVSNbXjGs1uQggjGI/aOG9nGEZVf5HE9z0KhUwwQfI6sdg+LCuG5xVZWJjh/vtPrSl81OgzqoUK1442FhtgqZRHMnmQffuOYZoVA5FK3YUQZhtXqdG0ByEsFhZm6txf+Tyo2RZQUaI1yedvlT30ZHKAvXt/goGBt5FKHWVi4tya1lDvM9rhvRpto+1NeWtlK5vy1lNyd/78E7z00mfxfQ/LimLbvUhZwrK6yGSuoCujNLsFw7Do7j5ALvcGjd/3goqhsDAMCyldQBKL7cN183heEdOMBtppPh//+Otb9Se0mo6t/ATtWTQkLLnLZidrSu5WUsIcGxvl4sURIpE4Uvo4To58/gZ33/0+fv3Xx0ml7trCv0CjaS++77KwcJOVL5AkQkSwrBiRSAwpS9h2nGTyMIXCbTzPQQgTz1OzuW07uVXL1yxBG4sGVJfcheqwYY33Ss8plRwcJxvILwuk9Hn55Wf47Gd7yWSubt0foNFsA3zfbXqMlC6uW8B1CxiGxU/91CeIxXoIPY7Kl0RUXZuPjY0yMnKCp54aZGTkhJY032S0sWjAekrupqcvkc9P4nkOvq+0oEIcZx4t7aHR1COcw+0RiSS4eHGEXO4WqdRdmGakqnT2LorFLLA+z1+zMbSxaMBaS+7GxkYpFjMrNOHpXIVG0wjDiGCaNo4zj+c5eF4Rw4iwd++99Pe/lb17763psViP56/ZGNpYNGCtJXff+tYn8X3tOWg068H33WCOhcHiYhrDiK74+dPNdluPNhYNWG3J3djYKL/3e/czNfXDIPSk0WjWjsTzXHzfoVTK4ftF7r//VMPPXyPP37aTOo+xSejS2Q1QiZvexHUXg3s76/XUaLYjPT2HMU27YU9EtUBgOO+iUJgDBLFY73adgaFLZ3cFo6Nw4gQMDqrvo6PluKmUru7U1mhagqCrq49CIUMmc5Vnn320rndQz/NPJAaIxXp1HmOT0Kqzq2F0FE6fBtuGvj6YnITTp0n/gzyxvkOYZhTPczBNG89z0N6FRrNeJIuLc8HUvAiOk28oKT40NFxz31NPDRKL9dUco/MYrUN7Fqvh7FllKOJxNRw4HgfbpvdWkVJpgXh8gNBAGIaNflk1mo0g8f0Svu8gpU8mc62hh1GNFg3cXPSuFlIdZrr/fvUVhpx+/GPorq28oLub499XFRuGYZFM3hlo3ngIofRvdGhKo1k/UnpIKRHCwnFyTfsotGjg5qKNBVTCTJOTYJpw6RK88or6eWwMZmbg4kV49VXIZNRzFhYYso+V46Ys5Dg8bfPzf9YHno+QQldHaTQbxkdKF8vqapp/0KKBm4vOWUBtmOnVV8GyQEp44w3wfRV68jwoFuHqVejvV8efUcKCQ5eBJ06DvQe6u4mWZijYbofXPmg02wPfL1Eq+czPXw8qnhSNhD7rlbfrGdwbR5fOggo39fUpo/DyyxVjUSgoo2Ca4DjQ1QWLi5BIwJe/DMPBG+7ECeWVxNXoxvP33ODcT06qx5YajOqXWxsTjWaNCAYG3sp9932IixdHakpnl5bJjo2N8u1vf5Lp6UsYhk0icRDTjLSznLajP/E6DAXKWCwEibFoVHkTYTe2Yaifu7rg3nvhbW+DPXsqhmJ0FL7zHXjttXKY6sHLd9SXgVpql+WSL41G0wTJ1NSPeOGFf0Gp5DQskw37MG7fvoya4+2RzV4LOsV1Oe160MYC4MwZ5Tnk8yrE5LrKQESjKvzk+zAwoI5dWICjR9XPTzwBjzyiwlOuqzyRq1chk2FgPlprAFZjDLTB0GiaIqWP75dYXJyuub+6TLa2B8rEMExAkM/f0uW060QbC1BewtNPw8GDyjAcOwb33QeplPIs+vuhp0cZE8dRxmV0FD77WXV8JKLCVqWSun3jBg//RS9xI4HhsTYjoL0NjWZVeF6RQiFTvl1dJhtqR5lmtDy2VQgDzyvqctp1ohPcIcPDldBSNaOjKgE+MaE8ijNn1HEnTihvwrbVcb6vbrsueB5DryY4ObXAhfcYpPsMonmXbBwWEmtcV2gwOjraqdFsDpnMFaQ8Us5FhGWyvb2DZLOTJBIDZDJXg6iyKsPV5bTrQye418vgoCqpDUNWjrP8GCNw3KrUaMfugQvH4cpdINfr12nDodEEGBiGyZEjD9RUOVVrR/l+iWx2Et932LfvGO9735Ptqobq6E+u9izWy+CgMhRTUyr8VI86kuVDl9XX2D3w9Ueg2LWO3629Dc2uIZzRvRylliAAb1k5rPr56aBkdoLDh9+lS2Y3iDYW6+Whh1QVlOepfMUaGboMH3oWvjEMmb7mx9dFl+FqdjShoRAIYVQNFhPBsCQL3/cQQlU3NdOOAt1zsRF0GGo9hB3fmQyk08pgbIAwNHXjADjr8TSWog2HZgdgGFZZMsf3S5hmtDynWxkPH5Akk3cCPh/4wOdXNAT1ZM23uOeioz+Z2lishxMnlAzI1JTKSxQKLTv1+Z+Gcz9Da95WHf3W1OxWhDCJRnuw7TiLi2k8TxmIBx74LS5d+jqzs5eDudxR4vEBDMPCsmwcJ7fMENx//ykmJs6RTo9TKKSxrATJ5ED5dzlOnmTyIKdOPb8lf9pW/JLNQoehVkNYETU+rnIVf/mXkMspj8JobfXxgy/CHZPK05juh3w363+L6dyGpsMQwuDBB38bgBde+EyVvprB9773e7zznb9at3O7UCiQz0/VGJFSyeGllz5LKnWUWKyP+fnrlEoLRCIxYrEUoCXM14I2Fs1YOstibKwiJijEuvIVzQiT4CFf+BhMH9zACZcuURsPzTYlFtvDK698nenpv8b3wwpDgZSSxcUZXnnl6wwPVxLXvb1HOXr0IV544dMIYSKEiec5zM9fDZ7nYdtKhseyunDdIvn8rbKx0D0Xq0cbi3pUexLptNKC2rNHPZbJKK2oDeYp1sJjv6++n/9pOPdeNv5f04lxzTZlcXGWxcXbS+6VSFlCSouZmVeXJa5HRk4E3do+QgiEMPF9D89ziETi5eMSiQHS6Su4bgEpZdkr0T0Xq0N3cC+lWq68r0+Fm6amKt5EsVgRF4TlnoXYvN33wRfh8c/AQ89TX3tqPXRWykqz41npDenV/Xil0+NqTAA+vu8jZfixlHR19ZaPi0ZTJBID2HZcS5ivA+1ZLKVarhyUgGCxCLduKfmPaFTd7u6GZFIZFSlV7qLyLt1UHnxRfY3dA1/5CMiNzljSnoamI5D4vs/Y2GjNBh92a6dSR8jlbuF5RQzDoqdnECEEjpMv5zdM0+aRR76sDcQ60J7FUsbHa6fiDQxU5MqlVAbD99X3+XmlC2Xbla8tZOgy/PbvwOOfgrtea9FJtS6Vpu00vmKR0uO5536pZmJeOCFPCIu+viF6e+8mmTzABz7weT0MqYXs3tLZpRVOoebT/ffD5cuqOzsaVcaiUFDhqD17lDcxPw83b6r7YzE4fFipzZqmmnfRRs7/NJx7iNZeBmhvQ7MNEMJCCJW0Nk2bw4ffXe6lqDTbqaT3Nm226+hP0u40FtUVTt3dSnbcceDUKfjCF2B2tlIS63mwbx986Uvq9unT6th0WhkQIeDAAfWz46gQ1Tbhc/8rzK+3O7weHf1W13QuqmNbCIHnqc+XYdjs2XN33aa6bdyl3dGfoN1pLJZMtgOU/PjNm2rjd12VoygW1dS8e+6BJ5+ERx+FbFaFoSxLGZRiURmMgQG4fVs9d5u9ps+ehJfvp7Vv1Y5+22u2O3v23EM2ewPXLWCadjDDolas07aT2HaKffuGyk1126BLeyU6+lOzOxPc4+Oq0qma7m5lCO6+W23+KVWHXZ7Fffq0MiihOGAoT27bSkgwnYZDh9SxjYQF28Qjz6mvliXEQSfFNZuGEBFisRRSSjKZCTyvvkJCqbSI5xWZmqp83sKhR2FvhW3HcRzqakdp1sbuNBaDg8s9i4UFlY9YWFh+f7Go8hWxmDIoIWHISQgVgrp6ta7S7HYhTIiDMhxf/iityW3oTnFNCzFNK/geoatrLwsLUw2O9JHSwPcrod90epxYrPZCUHdpt4bdWQ1VPUZVysoEvE98ovb+W7fU0KP5ebh+XYWe6iFl7dzuDmDoMjz+6RZWUYGupNK0AAPf93CcPIXCHMViuuGRoZCgaUbL99l2DzMzf8XU1Mvcvv0qhUJGd2m3iLZ7FkKI9wNPASbw76WUT276Lw3HqNabgPfOd6r7L11SjXj79qkQU2hENknio1384h9Vfm5pbkN7G5o1kkgcDMak+iSTB7Esm8XF24HC7PKrECFMEokB+vqGAJWvyOengtyGkv3IZK7Q1dXHz/7s765rTds4Wb7ltDXBLYQwgVeBvw1cB74HfFRKeanRc7ZMdTZMgruu8iqqlWV3mMGox6c/CX60+XFrQhsOTUMEicQhIhG7nIx+6qlB8vkZPK+I77ssNRZdXXuJxVLl40dGTpDNTiKlW27OE8Ji7957+JVfubjmFW1CsryjPwHt9iz+FnBZSvk6gBDiq8BJoKGx2DLGx1XfxLVryjhEItsucb2Z/LPAv2tp+a32NjRLEMIMNJ0MXDfHBz9Y6a7u7R3E913y+aklw48ABKVSnne96x+Vjw/zFUIIolFVoCKlpFCYW9fadLK8lnbnLA4B16puXw/uaz9hElwIZTQiEfUlxKbqP203PvHvVIf445+itXpUOrehQXVkW1aM3t6jxGJ7amZoLyzcJpe7iVKPDd98gnj8IAcPvp1U6igTE+fK5+rtHaRUWqg5/0byFen0OJFId819uzlZ3m5jUW/XXbaFCCE+JoT4vhDi+9PT01uwLCpJcM9TIaiFBRWSCktqdyGPf7rF0iKgjYYmKJG9QjSaBCrhH89z6Om5CxWtllhWjD173kRPzx3A8o07lP1wnDxSShwnvyFV2VYbn06n3cbiOnBn1e3DwI2lB0kpf19K+Q4p5Tv279+/NSsbHlYyHp6nqpyEUNVQ+bwSF2xUGbUL+MU/Ukaj5eq32mjsUtSc7TANGIZ/pHRZWFADjUAlucM5FLB84x4aGm6pFlSrjU+n0+4d73vAkBBiEHgD+Hng0fYuqYqeHhWCCru1fV95F44Dd90Fr7XyErvzCNVvoYWVVLrZb9fheQ4guH37VcbGRkmnxxHCZH7+GlLKIFfh43lF0ukJUqm7Gs6iWDrrYiOo8zzdCZpTW0Lb5T6EEB8A/g2qdPZLUsrPrHT8ls7gHhxUxmJqSjXgRaPQ36+6tBMJmJnZmnV0EH/4C3DlTS0+qTYauwCl/5RKHca2E9y+fRnXLbDcdTXo7t5PIjGAEFAszndSSWtHv5PbbizWypYai0YaUrOzyoBoVuTZk/Dy32jxSTv646ZZiXAsqppb0ajy0KC//y2USrntqv+0Eh397m13zmJ7MDqqDMOBA0rWY2BA3X7ooeUd3a+9pg3FKnnkOZXbeOsPaF0+Quc2dixSevh+CX8FJQQhBLOzr5ZLWoUQ2HYcw7C5cOHsFq5299HunEX7CeXKHUd5DEKoyqexMbhyRcmWnztX6eg2tH1dK6GQIbQwTKVzGzsUCaw8315KdElrG9DGIhyjOj2t8hOGoSqgMhmlInvuHDz/fCUk9frrlWS3Zs2E8iJaWkTTGEEj99E0I+zdey+l0kK5WQ5WX9Kq5TvWj75MDseoFosVryGcU9HdrbSjqo+LRrV30QIeeQ4efyJo9nNbdFLd7LfjeeCB3+Lhh59cV0lr2L+RzU4Si/WRzU4yOnq6ZkSrpjFNdz0hRI8QYlngQAjxts1Z0hYzOKjCTtFoxVvwfXV7YUGJDIIqo/2rv1JjU72V3WTN2nj8M1W5Dd23oakTYxTCoL//bTz44G+vu5+iWr5ju+c6hBC5Jo8fFUL8aI3n/EMhxIfXu6YVw1BCiI+gylqnhBAR4BellN8LHv5D4O3r/cXbhjNnVM4ilVIJ7NBgpFIqj3HmjMprTE2p26apei00Lad6SNM3PgCZPS04qc5tdDACw7ABj66ufbzvfRVB6vX0U+hZFxujmWfxT4G/KaX8G8AvAf9RCPFI8NjO+OiFcuVDQ2p6Xne3qogaGlL3Dw+rvEZvr/IybLvdK97xDF2GX/+3lS5x0cowlWbbU5lPIQGf/fuPcfLklzacW+hE+Q4hREII8W0hxF8KIV4WQpysetgSQowIIX4ohPhTIUR38Jy/KYQ4L4T4CyHEfxZCHGzJWlbqsxBCvCylfGvV7YPAfwJGUF7GlnsWW9pnETI4qAxJKCD4gx/oUFQb+MLHYLolb/uAnXG5s+MQwkAIC5AcOfJAeb72RtkG87lX/Y4TQuSklAmhXohuKeW8EGIf8F1gCLgLGAcekFL+NyHEl1Bq3U8B54GTUsppIcT/AvyslPKXhRB/CPwnKeWfrmfxzaqhskKIN0kpXwOQUk4KIR4C/h/gzev5hR3J0jGsodHYBXMtthOP/b763rKRsDpEtS1REh8OlhVraYioQ+U7BPA7Qoj3ojJ6h4CB4LFrUsr/Fvz8R8A/Av4MeAvwX4Tap0xgshULaWYsfpUlHyMpZTaYbveRViygIwjzGqDCVLYdFnvXDkXSbAnhSFhoobehDce2w/f9loeIWqkdtUX8PWA/Kh1QEkJMALHgsaVXqkqREX4spXx3qxey4rWZlPKilPJynftLUspnwttCiO+0emHbijCvcfAgzM3BPfeocav79rV7Zbuex35f5TYefQbMYotOqiup2oq6IrYAf9cqvFaRAqYCQ/EzqPBTyBEhRGgUPgq8BPw1sD+8XwgREUK0JArUqqa8WPNDOpzhYfUVMjoKf/fvtm89mhqGLsP/ERTL/J//GBaTLTipbvbbIsImPAOQCGFgmjZ7997TaV7AZvAM8P8KIb4P/AD4q6rHXgFOCSH+L2AM+D0ppROUx/5bIUQKtcf/G+DHG11Iq4zF7rsOe+YZpRml2Xb8k99V31sum66NxiZgYFlRPM/BMCykBMMwSSYP8PDDT674zJ3cjS2lTATfZ4BGIaVjDZ77A+C9de7/xY2saXe1IoeCgYOD6vvoOjo3R0fV8595pvmxmrZS3SXekul+ukO8hQiEMDFNi2TyEL29g5im6qlw3QLFYo4LF8427K5eTTf22NgoIyMneOqpQUZGTuhO7Q2yKmMhhFhmwYKqqPLNVi1o0wgFAycnVRns5KS6vRaD8cQT8MgjFQkQTccQTvdrmQKuNhrrxjAs4vEDgBp8lE6PUyrl8X0X33cRAmw7uaIcR7NubC3t0XpW61l8TQjxm0LRJYT4d8Bnqx7/+5uwttYSCgbG46rkNR5Xt8+ustX/Z34GHn9cVz91OKG38egzMHCrBSfU3saa8X2XfP4mYY+XlB653CSuuxgcYZDNXsP33YZyHOn0+IrKs50k7dEprNZYvAs1K/sCahTqDeA94YNSyjVplLSFUAiwmmqhwJX4hV9Q6rOaHcPQZfiVL2pvo32EL1btFmQYNpZlA4J8/lZDOY5m3djNjIlm7azWWJSARaALVfk0LqXsLI3uUDCwmmqhwJX44z/elCVptgfVuY2e2RacUHsbq0S9SIZhYZpdqDyG2pKk9HGcHFNTL1MozC0LHx0/fmZF5dlOlPbY7qzWWHwPZSzeCTwAfFQIsa6W8bZx5kzt1Lt8viIU2AwtHLhr+MS/U0bj8U9BV7YFJ9SGYwUEhmFiWV1I6WFZMcDHdUs1Y1VtO7Es39BMebaZMdGsnVXN4BZCvENK+f0l9/19KeV/3LSVNWBD2lCjoypHMTGhPIozZ2p7JxphmnrY0S7l/E/Di+8Fr5VjwrZ/Ocgmo4yEKpu1Mc0Yi4sz5cdCy2pZXfT0HCIaTeE4eZLJg2vSiaqU1m4baY9t858PVDieQsmB/Hsp5cp1yqzSWGwntlxIcHQUPvzh5SEsza7j/E/Df3sPlKLNj10V22br2Fri8YOYpsnCwgymGaVUWqzyJCr7kWFYCGEhhMD3PQzD4CMfebbdG/5G2Bb/cSGECbwK/G3gOipy9FEp5aWVnre7+izWw9mzcOQIWHoC7W7nwRfhnz5ZkRfpXnE8zSrYZSGqVOooAwNvIxKJ0dc3xEc+8iwHD749yFlEiUS6gvkVCt938bwCrruISpEau6/8VYj3I8S3EeL14Pv7W3DWvwVcllK+LqV0gK8CJ5s8RxuLpoRVVOEsi2i0ojqr2bUMXYYz/1oZjfhGjQbsaKNhmlEsq5s9ewa5774Pkc9PMTHxX/nKV36OiYnzeF6xKrG9ND8oAIGUHiDJZK7y7LOP7g6DoQzD54GDwGzw/fMtMBiHgGtVt68H962IvlxuRihPnkopD+PWLZXw1vMsNCij8Rv/Wv08dg989SPgRzZwwh0oLeJ5RUwzyrVrf86VKy8EGz+EBZVSgu+XMM0IlbB4mLuQ5Z+llAgRwXHyjI6eBrZsDkW7OAMUgTAGvlB1/59t4Lz13l1NL1W0Z9GIUBrk0iWVEL91S83hPnQI7rpLGZFoq4LXmp3A0GX4Z7/TolniOypEJfC8Ip63EBgKEQw4qmw/vl/Cq7oAC6uk1NQ8ARgYhgH4WFZstzTYDVIxECELwf0b4Tqqby7kMKp3bkW0Z1GPUBrEtpVxiERgZgZKJTh2TFVRPfYYvPnN8N//ux6ApKkhnCUOLRrU1PHeRqOxCxUMI4IQAsOwAEEi0Y9hRJmfvwr4GIaN55XKkiDz89cpFOa2aP1tYxwVeqo2GN3B/Rvhe8CQEGIQeAP4eeDRZk/SxqIe1dIgAAMDkEioeRbPB6V7g4MwNqYNhWZFwkFNY/fA1x+BYtcGTtYxRqNS/gpqTKphRPC8+gNHhDBqxqdWl7zu23cfc3Ov4boFpPQRwsI0I3ieQ7HoMjY2upNDUWdROQtQBqMbiAb3rxsppSuEOA38Z1Tp7JeklE0lzHXpbD2WztwGZRTm5uD119Xt0VElKqi1ojRrYOwe+NbDMNXPxoPA295oKI/BNFUSx3WL5XxFNdFoig996CsAdSXHx8ZG+drXHsH3fQzDDHIdkq6ufvbtG2rZjO4tYO3/MZXMPoMKPY0DZ5FyI/mKdaM9i3pUz9zOZFS+YnFReRejo5VBSNGoNhaaNTF0WX2FjN0Df/o/g9Pd+DkN2aajYIUwgkS1xPfd4LsP+KgSfxl4CQY9PUf4O3/nCwCMjp7GMOwaldgwiW3bPbjuAp7nYJpR4vEBotGena/1pAxDW4zDUrSxqEc4czuXU4ZCCDAMZSxOn1YjVkEPP9JsmKHL8FtnVcPfufey/k/kNjIcUvoYRoRIJE6xmA4MhvIyALq6+jh58g9qwkcjIyfKKrEAth3HcZSnMTQ0TH//m8lmJ8uPAzhOXms9bSG6Gqoe4cztXFBAH42qstmBAZXL+LVfUyEorRmlaREPvgiPf6bS8JeaY/2VUNugiiqVuou+vjdhmjHCXgnLipFK3UUstmdZJVMzlVit9dR+tGfRiOFh6O2Fu++G+XnlYVy9qjyMQkEZDcPQmlGaljN0GX7936qfN5QYX150tGXEYqngJx/TjCGlkvOYn7+KYdjLKpl6eweXeQ7VKrHKC3l6u2k97Sp0gnslTpyAH/4QZmdVgluISvVTJKJuO87WrEWz6xm7B547CfnEGp+4xQZDCIs9e46Sy93CcXKEVss0owhhBEONanWewsl2hmETiXRTKi3g+06NkuwOYBtlltaODkOtxB13wO3bFQNRbVhdV8t+aLaUsFs8lFC3F5s/B6ht8NtwiEqFlFSiejmGEcWyoszOvobjVGu8q6S2lB5CQCy2ryYU1UxyXNN+dBiqEaOj8PWvN35cSijWrxvXaLaCDz8LX/swuEuFBFYTflp3QjycbicxjCjgB3OzjaBKqZ9stn4zsO+XsO1Ew0qmoaFhbRy2ACHEl4CfA6aklG9Z7fO0sWjE2bNK/8kwloegOix0p9mZDF2Gj/wpXDgO6T3QOwfHL6jHvvrz4Ne/+F/OmnMbSqfJNE2EMANDYdPXN8TMzCWUJEdX2fEulZQLZJoR9u69F1heyVRpxKvtsdBsCn8IPA3832t5kjYWjRgfr/RReJ4yENpIaLYZS/s2Qt57Hs4/BHK1geY1GQyJEGa5V0IIA88rkstN4boOqq+iEMyjqHRze55LoZDGMCI1lUzV+Yp6PRa7mU99anlT3uOPb6wpT0r5ghDi6Fqfp3MW1YTigYODkE4rY2EYalKeNhSaDuLBF+GjX4GBW2xc1LAO8Xg/KgTl43klPM8hl7tOxer4+L5TI/FhmhEymStYll2Tj7hw4Wy5x0IIgW3Hd4tQ4IoEhmKZRHlw/5ajPYuQavHAvj6VwL55U0mT53KVnopwxKo2HpptTrXXMXYPfOMDkNkTPBgmuptdLjbIbUSjKSKRbjKZa8EMChF0bqtO7VoEicQBksk7cJw8XV17azyGdHqcWKyv5hnVPRa7mM2SKF8X2liE1BMPBGUohFADkMJQlE5uazqMsHfj/E/DS+8FX4DpgW+oL3yQEWUN1KbvrSgWOz9/LdB88sshqVB+fKmQ4J49byr3XUQi3UxNXeKLX7yfmZlXEUL9Ptd1SdTdU7wAAB+dSURBVCYHys+p7rHYxQyiPIpqWiFRvi60sQgZH1ceRTX9/aqfQkr12OSk8jaWJrw1mg7hwRfhjslKUjxahFwcYkVBBIvSoQH8bpv7X8xx7r6p5Z5HYDB83w0S1/ViXI0/F7ncFIXCHIXCbZTgqRqOVCqpCqpEor/cY6G7szdNonxdtM1YCCHOAv8T4ACvAb8kpUy3az014oEhCwtqnCooOfLbt9UsbtfVhkLTsdQTM7zwvhjphEvvzRzHf+3LDP3mI0wkBBN3SnVdFIagJAjDaCg3Xs9QZDJXkPIIphmhUJgJkt5WMMyoguvmKBQiuju7wqZIlAshvgI8BOwTQlwHHpdS/oemz2tXB7cQ4n8Eng+01f8lgJTyN5s9b9M6uKtzFt3dylA4TkU08JFHVK7CstQQJG0wNDsB1SGnhnpVy/B3dTF21OVrH3LxRVCGG3gVhhnB90urOr1p2ghhAT6HD7+b6elLLC7OBsOO1DFSSnzfo6fnEB//+Oub9ZduB9bcxbsZ1VDrpW2ehZTy/6u6+V3gw+1aC1ARDzx7Vo1RPXpUqc8OB1c3qZRSmV1cVEZDh6E0nY4QqmAD4OWX1YXQwICqCPQ8hv7K5YEXghxH8FY3PBBGZS52DXWS4UKY7Nv3kxQKc5w69TwjIye4fv3Py2W3EKrUWjpHUYfAMGwLifLtUjr7y8BouxfB8LCahPf66+r7cJUbfOwYHD6svA7b1gKCms6jt1cZg3i80mzquspTFkIVbVy5okKuR44AKsfxkT+G/ikwpPrqyVuYpl11YlE3GW6KCFJ6Ncnq48fPEI32IKWL53n4voeUPtFoUucotjmbaiyEEN8SQvyoztfJqmP+d8AFnlnhPB8TQnxfCPH96enpzVxyY86cUWGpQkF7FJrO5MgR9d71PKV7FlHzJXBd9d42DOVdZDLKsFgWCMHQa4JffSbJz587yOF0Aum57N37E8RiewCj9vMQXkNJkK6LEFZNsnpoaJiTJ7/E/v3HEEIghGDfvvuWzbfQbD/aqjorhDgF/ArwsJRyodnxsMWqs0sZHYVHH1UNexpNp9HdrcKooAyD51XEMIVQ3rIdeAxvfSu8+qoyIlKq26BCscEs+rGxUb797U8yPflDhBT4yKAdD0wEUkj2H3wbDz/8pDYEio5WHm1nNdT7gd8EHlytoWg7w8PwiU/A44+3eyUazdqp9oq9YBZ29cWiYaiQVCLQQB8YUGGpsHw8LPo4U/EShoaGGfvQ/VwYuMxUqohvgOkL9mdsjt+6h6FPXdzCP1CzmbSzz+JpVBnYf1GldHxXSvkrbVxPLaOjKtk9Pq7KasNk97lzSgZEN+VpOo1meTbTVMYilVLGwbJUf9HAgKqSWlr0ETD0D55kqG4l4ZOb97dotpy2JbillPdIKe+UUv6N4Gt7GYrTp1XfRdiMd/q0un98XCW6NZqdRNhoGovBzAxcuqQ2/8ceg717V87ThZWEBw8qo3LwoLo9rENPOwk9Ka8eJ04sb9C7dUtJfxQKyl13XfWlq6I025lmJd5CKA8ilOM/cEApFywsqNyclLBnz/LeI20I1kNH5yy2S+ns9mJ8XH04QjIZmJpSxuLQIeWqh8nBpdPy9PQ8zWYSiahNfbXvs2YXg2F1VGgoYjFVOvv66+o9n06riyYh1HfbVuFZza5DG4t6DA6qq6iQW7fU964uVVJ4110V+fLqD61u1NNsJnv2bM4ALtNUoaZoFK5eVRdDlqW85mJRXSyFdHerplXNrkMbi3qEPRX5vPpQhlUkoRJtKgU/+ZMqd/Hgg+oDpGdeaDaT3l545pmKLMdGCD3ias/42DEVejUM9VV9XHixBLV6aW1kbGyUkZETPPXUICMjJxgba39P705HG4t6LE3YxePKUKRSlWPCD82ZMyo0EJYiajSbQdBRDVSS0esNeYbGxrbVl2lWLpBCCX7PU/cbhurNkFJdPFWVzraLcLJeNjtZM1lPG4zNRSe4V0M9kcF0WsV3b95UrnqHvY6aLcQwNlYIYRjQ01PpdWiliKVpwpvfDBcvwv33w+XL6vzRqLpAKhRUrm7Pnoals1vNyMgJstlJbLtSgOI4eZLJg5w69XwbV9aUjk5o6nkWq2GpyKCUMD1dSXJrQ6FZiY0YCiHg7rvhxg21cUci6qIl9AI2SuidDA4qg5RM1lY/2TZ8+cttNxDV6Ml67UGHoepRPYv7xAl1OxQZ/PznlTcRbgDaUGg2i+5uZRxSKZV09v3KTHjbXl5c0Yylx4Qej+OofiLHUcfY9rbul+jtHaRUqhV90JP1Nh9tLJayUkMeKO8inMety2Q1m4nvq3AQVEpmw4sU01SPRaMqHNrsoiX0gMNcRzSqzhWL1ZbG9vaqyqh6ysvbhOPHz+D7Do6TR0qJ4+T1ZL0tQBuLpVTP4q5XWz4+rj5g2lBoNhvXVQ1y+XwlROT7KvzpeerxPXtUf8RKhOOBQ0KtJ1DeQzUdUBo7NDTM8PDTJJMHKRTmSCYPMjz8tBYr3GR0zmIp9WZxV3+ABgfVh7RQ0BVQms0hDAVZlnrfJZNKwPKd74RPflKpwYIqd/3Qh+DTn175fAcOKHFAx1GFGcWi8kyOHq01IrBtSmObEYoYarYO7VksZWlDHtR+gB56SGnn1HP7w6ljGk01oeaSba/u2DAvIaX67jjwhS+oxy9eVKWsi4vq53PnVj6vaUI2q3IPQ0Owfz888AA8+6w6Z3U/0TYpjdVsT3Tp7FKazeI+fbr2Cs3zKhpSxaJ6TKPp7laNm0vfPx/9KMzPr64wwrLUlb/vq6/77lMGoprBQWUQxsfrV13t3Qtve5vKP9QjVFeuN0pY02o6OnatjUU9Gn2A6gkM/vCHKln4Ez+hbv/FX2zu2jSdQTyuQj/HjtVuwKOjKpT0ox9VDEYwka5cjRTeL4QyFqHQnxCV4UUh4XvSdVVSutpg7N2rKqm2YUXTLkUbi62krZPywqu4qSnlRUSjqqTRdVVDE6grvNnZ9qxPsz2wbdUbsZJCa+jB3rypDEFYGls9KyWsXIrFGhuLpZ7w1JQKk6ZSyw2Vpt10tLHQOYu10NOjhNYcpxJLDmcXh3Hf/v6Kto5md1IqNVZoDXt4HnusIqtRLCpDYZq1Mh6h9EZYAXXvvct/11JpmqEhlY+4eXPblr5qOhO9q62V6lr16q+JCTUwJvzgazqHauG8MCQUSnZXsxZZ8FdfVUajugy1uofHNOHaNfVeOXCgdib2wYOVdYTr27cPnnyy9lxh4+jZs8qDqO6NqNdYqtFsAG0s1sL8vJInj0RU6Cn8MPu+cvmPHlUbge7BaA/rfd3jcbVhR6Nw553w1rfCW96iJOmrxSOrQ7ZWk6rzfB5ee02FhaobOsMenqkpdY7QK737bvW+EkKtJVxPfz+8613wpS/V5j1Wahxt9rhGsw60sVgLg4PqA33vvWpDsSx1hRg26cXjrRV506yO6rDNWrAs1cNw772qc/nd71ZX5xcvqq9f/EW1kdczQrZd6VGo93joXZZKtSN5w6FaxWJFDrxYVEbpyBH1vGbhpGaNo80e12jWgTYWa2HpnIsw2RjOuQBlODRby2qNhBAV4x4Orgo7pOv1F3zuc+q4ri51lV9tFELJl3BokGFUvkOlV6JUUnmuRx9VOa+whyeU21gq6REarJXyDUsnOUJt42izxzWadaCNxVpYmkxMJNRmUx2q6O2t/KzDUduHeBz++T+H975X/f/e8hbVt+D7jQXzstlKo+VS8b7Qq8jnlWfylrcoz9OyKgKAYQ4iElHH3bql+nPyefW+CWe4r2Sw6tGscbTZ4xrNOtDGYq1Uq88eOaLCBK+8UtkEbBv+3t+r1d/RbA0rGef/v737j427Pu8A/n7u7LNjx0mcgO38oIkJBEgp69ZSGrVILGMT3iqqIaFupVtVKlWi+8GiLRmMrurQpmaLtO0PqKZqS4eWQjcxuiEmr4O5rJPG0jIGpSgFAh4hJLGTkh+OTS7n87M/Hj/6fn25u+/57Lvv9+77fklWfOez/YmdfJ/v5/N5nuczPQ088kiwEexLTdXu4vv6FrZ0yWbt99rVFTzf0WGp0h4IvCbCizP9eNLubuvjNDBgwcn3uaICVjmlM9zSQBP1caI6MFjUwzcQL160De9i0TYz33jDZht33QV86UvRm6C0dLlcsPQTFZzD6/a1ZAvt2mUXct+Hmp21t2LRnu/osGWmyUkLIgMDwFVXBYHEN6z9SN6eHputjI3VHrDKKZ3hlgaaqI8T1YFFefUIV3IfO2azC+/js369XUBErDgqk2ELkEar9QCqrVvtIv/ww5VbupReUB980PYuzp+3G4H+fkt5DfdjKhYtMFx+uV34R0dtj2J62mYUfiTv9LT9+6jUeoPaXUuvSzNY1GN42FISz52z2UT4Z+jr2tmsvU1Ps+4iCTo77e2mm+xxaduWWi/kw8P2uYVCUHPj+xof/3jw+dV6jPEOP61aOlhwGaoevoE4MWEXDM+s8T8vXrSLiZ8/QI3lmUd+IFAp/934uv1SsoVWrQo2pr3CulCwj4X3BNK8FMSCwLbEYFEP30C8cGHhEoifZgbYn6dOMSOqEcItMfr6bO/Al4X8bj/8Wg8m27fbxXqp2UKe4ZTJBL/7TZsuDQSeDJHgU+eWHQsC2xaDRT38rrG3NwgOXqBXLNr7fhHhGRfLK1yAJ2L7AZ6d1NNjGWpe65LJWGX01VdbRbS3y1hKtpBX8Xd3B0V9W7c25u/ailgQ2LYYLJbife8LTjXzxoKZjGVCrVlj74fbgtDSeeGbiN25rlplz/nF/8QJK5QbGrLah3JpqUtZIiqt4t+2zR6zhsGwILBtMbezHuHNy82bbap94YJV+vb22ullQ0OWQePVwoC9xvcwas3goYWKReD977fjRJ991po3nj9vAWJ62t4/f94u3nv3Vg4AIyP1LQvt3m2/e2DhxjVrGIwnAISTB1gQ2BY4s6hHeKq9Zg2wYUNwDOamTfafZWLCev74JqgvTw0NBfn3nHEsji/vTU7aedRjY0EDx/7+oPitq8te04i18jRvXNeCBYFti6mz9fDUWb/Yv/aa/YdQtaUJwP6THDli78/O2gVscNAueOvXAz/5CXD48KUbrVSZ92fKZCwFdmws+F28/rplJfmmc7Fo+xWsa2g+HtVaSUvfHXIZqh6lU+183i5i4bTNnh7b/PQNvkLBCvguXrTXHz3KtNpaeVAuFGzmED4nwn8X+XxQMe/N+bhWHo96l/go0bgMVY/Sqba3fZidBV5+2WYak5O2RPLQQ7ZX8eabQZfa8XF7fbjvEEXzU+PCG8r+u/CeTMVi0F6Da+VEy4bBoh6l69YDA3b3Wyza3kU+b1k5t9xirz9xIqgg9tbWVDvf91G1wNvVFayB++/CezJls3aAUUdHc9bKWYBGKcFgUa9wwdXwsAUOz/fv6rKN7GeftbXb2Vm7eDFQLF1X18K+TID9Ll56CXjySdvLWGwX13qxAI1ShBvcy6F0wxuwu+DTp+3PU6csUHh2FC1OT49lnPm5IUlpyBduKOmSMjZKopbe4ObMYjlUah/R12dnHLz3ngWKUkydjdbZCVx77cIDppKycc0CNEoRBovlUC63/MwZq7VYudLW0cN9hFzp43JN8NLGq7M9uyybTe6pbzyRjlKEwWI5lCvUGhiwQjGv4J6bC2oEcjl76+uzqu9s1j6WxlRakaCXkweIXC44H2LbtoWBeGLC7txfeMF+vkND8W0sswCNUoR7Fo0yPGxB4O23gyAxN2cpswMDVvnd02MptidO2DLL6dPpCxhDQ/azOHLEluv8EKk1aywpYHBwYTV8Pm97BFNT9vmeJpvLxVNJzQI0ql1LrztzZtEoXizmh+P40kp3d3AO8+nT1hH1U58KjupMm5MnLaBeuGAXfK/SnpqywHDxIrBxowWVfB647DJ7TiTo9Hv2bHydTWttQ84UW2pxsQcLEfk9EVERuSzusSwrX6Jwc3NBSmf4HObdu4HnnrO9jUzsv47m8yK6zk57rGp36F4nEW51PTtre0H5/MJzQ/L5ZG8sM8WW2kCsVycRuQLAzwM4Euc4GmJkxCq4s1m7yHV2Wktzrz72O8077gDeeceWonwGkjYzM8FpcyK27OQV8WHd3RYYurqCWZi39kjyxjLPeKA2EPet7F8A2AOgtTZOarV3ry2fXHmlLT0dO2bNA8fHgbvvtjvMuTm7gy4Wg83vtPFMMW/GmM9bUO0oaV22Zo0F1NWr7XWzs/ZzW7062RvLTLGlNhBbsBCR2wG8o6ov1fDaL4jI8yLy/MmTJ5swumXiWVK5HPDWW/bc5s22qe2Fel1dQQptoZDO2YUvQQH2M+josEON+voWZhrlcsD999s+T3+/XXDXrrXHSW4TzhRbagMNzYYSkWcADJX50AMA/gDAL6jqWRH5PwAfVtVTUV+zZbKhwkorfV9+OThhb3AQeOONIGB4y5BaWoP4Ma6trvQgqC1bgK99zd5vh0yj8GFZ4QOTkhzgqBFa+i6woS3KVfXWcs+LyAcADAN4SewuehOAF0TkI6p6opFjisX4uN0BOz8GNJ+3JZShIduzAIL6gqkpy/Kpph0CBbAwUPixqUD7tLr2GWY7BD5KrUTUWaRuZnH2rC1LeSuLmRnL8hkYsCDhF5M777Tll1bn6bDe8qTSvzkRYOvW4IAo9lei9tLSM4u4N7jTodz5F2vXWlttr/jev986p4bz9ffsqbx/keR9jUzGMpc+9CE7ZtbbtwMLA0X47+D1KKtXc/OXKIESESxUdUsts4qWVa4dyDe+cWlwKPXlLwOf/vSl9RciwTneSeOBwlNbBwYsAPisKnz2uAcO379ZscIec/OXKHESESxSIVzpu3u3rV/XUs174ADw1FN2Ac7l7O67o+PSTeHFCGcfLTdvaeJNEWdmrN7k0UctcFRqqDg7ax9nfyWiRGKwaLbRUeBznwMOHrS6i4MH7XG1gDEyYg31/NjQQqF8y/NahLu6uuWeoRQKVpEevvCPjNhsyvcvfCwdHfacF9o149AiIlq0hmZDURn33Qe8+65dKD319d137flqfYUOH66c/dTZGVRAR1m5Ejh37tKDmpw3PKyXf93jx60x4K5dwd9rZMRmD5UOinrzzfq/LxE1FGcWzfbaa0Gg8E3dTMaer2TfPmu0V8li9i8KBdtAXrHCNthXrAhmGj09tuR1112L+zuFedfYXM5Sgh95ZOGsiQVqRC2JwaIVjI8Hd/uVmg16a4yooPHee7Y0NDNjy2AbNgDXXGMX8ccft7v/AweWFjC8Mr1cDySeAUHUkhgsmm39ettvmJmx2UKhYMtL27ZV/pzh4YVBojQgeLX3unV2ke7srN7B1hsbAlbvUe4siGPHLO01l6uvG+7goP1ZmgZbLjOMexREicc9i2YaHbUgkc0GLcvn5qzP0d69lT9v9247Ge7s2Uv3E3p6gnMgtm8P7tDvvjuoCi+nowO47jq7s1+3buHFenTU2qbPzQVdcxfLz8wut8TULpXZRCnCYNFM+/ZZYOjrC05/6+gArrii+sVzZAR47DHgi1+0E+Xm5myJZ88eq8Uo5557gK98pXx6bbhIrlCwwDA8bG+33GL7DJ7i6gEtSmkq75kzNnvhEhNRW0hEu4/FaMl2H254uHmZQDt3Aq+/bmdlqC7McurstLqNwUFbIspm7bl83oJIf7+9HTlin1dtc935JrmfYCcC7NjBHkhEgQRW0daOM4tm8qNWvZoZaFwm0Pi4pakWi7YcFZ4dqAbna3hBHGAX+ULBUnn7++2wpomJIFh0dCxs25HNBs+r2lLY4KC1Fz99mr2diNoIg0Uz7d5traqBha2qG7FM44FpwwYLTn7R7+y0lNmpKQsgnkXlF36fgUxM2Kb76tXAoUMWJHI5y6by2WhfnwUUP/7UTU8zFZaozTAbqpkWkwnkx66GW4KUe66ScIrqqlXAxo22N/L440FPqh07LDCEs508aJw/b+duHDpkAaa/377GDTfYYUObN1sLj717mQpLlALcs0iicoflnDljF2M/Ia6WA3RGR6ufoTA6ameA+wzD+zqFG/yp2p/r1lkgyOetCnzXrmBzPer7EBHQ4nsWDBZJVHr+BWB3+IClu7rp6aWf+/Dgg8BXv2r7Ft3dQd3HwIAd/ZrJBCf3ZbM2o/AsJ9ZHEC1GSwcLLkMl0fi4zR7CCoVL6x3qOfehdCnrxhuBJ54Abr7ZTugTsYAwNWWBIpOx7+vncExOlq/MXsz3rLZ8RkSJxGCRROX6J3V2BpvRbrGZVL68dfy4pfAePx5suHv79B07gjRa38vw1NtMJuh2WxqoKgWESt+TAYOopTBYJFG5/kmrVln20VI2kvftsxlBb6/NIKr1bvI02WLRXutZUuFzKjxQVQsItXxPIko8BoskKpc1tX+/nQexlJ5K5Za3KvVuuuoqCxTZrNVOiJQ/oGh01E7zO3LECgDPnVsYEGr5nkSUeNzgTpNyG+fVNsnDWU59ffbc1FSQ8QTYDOLtt4PCvLk5q73wwrwtWxb3PYnaV0tvcLMoL00WWxQY1fBv506bQXR329fxGo2JCQseHlSaVYhIRA3DZag0We724L7ENDhoswrf3/AzM7zegi3JiVoel6GofuFlrbNng5Yivb1W3c2AQBTW0stQnFlQ/Sq1FGGgIGo7DBZUPy4xEaUGN7hpaXjqHVEqcGZBRESRGCyoMvZ0IqJ5DBZpsdgLP3s6EVEIg0Ua1HPhZ08nIgphsEiDei787OlERCEMFmlQz4W/XJv0xbZEJ6K2wWCRBvVc+Mu1SWdPJ6LUYrBIg3ou/Cy4I6IQ9oZKi3C7ce8Gyws/UTO1dG8oVnCnBSutiWgJuAxFRESRGCyIiCgSgwUREUVisCAiokgMFkREFInBgoiIIsUaLETkt0TkVRF5RUT+LM6xEBFRZbHVWYjIzwL4JIAbVDUvIgNxjYWIiKqLc2ZxD4C9qpoHAFWdjHEsVA0PQSJKvTiDxTYAN4vIQRH5DxG5sdILReQLIvK8iDx/8uTJJg6ReAgSEQEN7g0lIs8AGCrzoQcA/AmAMQD3ArgRwN8DuFIjBsTeUE22c6cFiN7e4LnpaWssODYW37iIWg97Q1WiqrdW+piI3APgifng8H0RmQNwGQBOHZJkfNxmFGE8BIkodeJchvonADsBQES2AcgBOBXjeKgcHoJERIg3WOwHcKWI/AjAtwB8NmoJimLAQ5CICDEGC1W9qKqfUdXrVfVnVJUL4EnEQ5CICDzPgmrBszCIUo/tPoiIKBKDBRERRWKwICKiSAwWREQUicGCiIgiMVgQEVEkBgsiIorEYEFERJEYLIiIKFJDW5Q3goicBPDWMn/Zy9A6TQxbaaxAa42XY20MjtWcUtXbGvS1G67lgkUjiMjzqvrhuMdRi1YaK9Ba4+VYG4NjbQ9chiIiokgMFkREFInBwnw97gEsQiuNFWit8XKsjcGxtgHuWRARUSTOLIiIKBKDBRERRWKwmCci+0TkxyLyQxH5toisiXtMlYjInSLyiojMiUgi0/xE5DYReVVEDovIfXGPpxoR2S8ik/PnwSeaiFwhIt8VkUPz/wbujXtMlYhIt4h8X0Remh/rH8U9pigikhWR/xWRp+IeS9IwWASeBnC9qt4A4DUA98c8nmp+BOAOAN+LeyDliEgWwMMARgBsB/CrIrI93lFV9bcAWqVYahbA76rqdQA+CuA3EvyzzQPYqao/BeCDAG4TkY/GPKYo9wI4FPcgkojBYp6q/puqzs4//G8Am+IcTzWqekhVX417HFV8BMBhVX1TVS8C+BaAT8Y8popU9XsA3o17HLVQ1eOq+sL8+1OwC9vGeEdVnprz8w87598Sm1EjIpsA/BKAv457LEnEYFHe3QBG4x5EC9sI4O3Q46NI6AWtlYnIFgA/DeBgvCOpbH5Z50UAkwCeVtXEjhXAXwLYA2Au7oEkUUfcA2gmEXkGwFCZDz2gqv88/5oHYFP9bzZzbKVqGWuCSZnnEntH2YpEZCWAfwTwO6p6Lu7xVKKqRQAfnN8D/LaIXK+qidsbEpFPAJhU1f8RkVviHk8SpSpYqOqt1T4uIp8F8AkAP6cxF6BEjTXhjgK4IvR4E4BjMY2l7YhIJyxQfFNVn4h7PLVQ1TMi8ixsbyhxwQLAxwDcLiK/CKAbwCoROaCqn4l5XInBZah5InIbgN8HcLuqzsQ9nhb3AwBXi8iwiOQA/AqAJ2MeU1sQEQHwNwAOqeqfxz2eakTkcs8qFJEVAG4F8ON4R1Weqt6vqptUdQvs3+sYA8VCDBaBhwD0AXhaRF4Ukb+Ke0CViMgvi8hRADsA/IuIfCfuMYXNJwr8JoDvwDZg/0FVX4l3VJWJyGMAngNwjYgcFZHPxz2mKj4G4NcA7Jz/d/ri/N1wEq0H8F0R+SHsBuJpVWVKaotiuw8iIorEmQUREUVisCAiokgMFkREFInBgoiIIjFYEBFRJAYLIiKKxGBBVIaI/KuInGGraiLDYEFU3j5Y8RsRgcGCUkREbpw/3KpbRHrnD+S5vtxrVfXfAUw1eYhEiZWqRoKUbqr6AxF5EsAfA1gB4EASO6ASJRGDBaXNg7A+RRcA/HbMYyFqGVyGorRZC2AlrGlkd8xjIWoZDBaUNl8H8Ieww63+NOaxELUMLkNRaojIrwOYVdVHRSQL4L9EZKeqjpV57X8CuBbAyvl28J9X1US1gidqJrYoJyKiSFyGIiKiSFyGotQSkQ8A+LuSp/OqelMc4yFKMi5DERFRJC5DERFRJAYLIiKKxGBBRESRGCyIiCjS/wMLvyFEt1exgwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 402.375x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Test cell: `kmeans_test`\n",
    "\n",
    "# Calculate the estimated centers\n",
    "centers_est = update_centers(k, points, clustering)\n",
    "print(\"\\n=== Estimated centers ===\\n{}\".format(centers_est))\n",
    "print(\"\\n=== True centers ===\\n{}\".format(centers_true))\n",
    "\n",
    "n_matches = count_matches(labels, clustering)\n",
    "print(\"\\n{} matches out of {} possible (~ {:.1f}%)\".format(n_matches, n, 100.0*n_matches/n))\n",
    "\n",
    "visualize_clusters(points, clustering)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Profiling\n",
    "\n",
    "What is (or are) the bottleneck(s) in this code? Let's use Python's [`line_profiler`](https://jakevdp.github.io/PythonDataScienceHandbook/01.07-timing-and-profiling.html) to find out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext line_profiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "%lprun -f kmeans kmeans(points, k, starting_centers=points[[0, 187], :], max_steps=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Observation.** You should have seen that there are two bottlenecks: computing the sum-of-squared distances and updating the centers.\n",
    "\n",
    "Let's use Intrepydd to try to improve the performance of the latter. Then, as an exercise for you, try optimizing the former.\n",
    "\n",
    "> Before we begin, estimate how much faster you can go, ideally, by improving `update_centers__0`. A quick-and-dirty method is to use [Amdahl's Law](https://en.wikipedia.org/wiki/Amdahl%27s_law). That is, suppose that the profiler shows that `update_centers__0` takes 60% of the total execution time. Further suppose that, magically, we can shrink that time to zero, while the time to perform all other steps remains the same. Then, the new time will be 100% - 60% = 40% of the original time. The **_speedup_**, or ratio of the time before to the time after, is 100% / 40% = 2.5. In other words, the improved time will never be more than two-and-a-half times faster than the original time, all other parts of the program remaining as they are."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Intrepydd reimplementation ##\n",
    "\n",
    "In Intrepydd v0.1, you'll need to write a \"scalar\" version of the original array-based code. Here is one such implementation. Notice that it has an identical signature as `update_centers__0`, meaning we should be able to \"drop it in\" as a direct substitute for the original algorithm.\n",
    "\n",
    "Take a moment to study the code to see that you follow it's steps, then read the explanation below for what it does to try to speed up the original code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting opt.pydd\n"
     ]
    }
   ],
   "source": [
    "%%writefile opt.pydd\n",
    "# opt.pydd\n",
    "\n",
    "def update_centers(k: int64, X: Array(float64, 2), y: Array(int64)) -> Array(float64, 2):\n",
    "    m = shape(X, 0) # type: int64\n",
    "    d = shape(X, 1) # type: int64\n",
    "    centers = zeros((k, d), float64())\n",
    "    counts = zeros(k, int64())\n",
    "    \n",
    "    # Sum each coordinate for each cluster\n",
    "    # and count the number of points per cluster\n",
    "    for i in range(m):\n",
    "        c = y[i] # type: int64\n",
    "        counts[c] += 1\n",
    "        for j in range(d):\n",
    "            centers[c, j] += X[i, j]\n",
    "            \n",
    "    # Divide the sums by the number of points\n",
    "    # to get the average\n",
    "    for c in range(k):\n",
    "        n_c = counts[c]\n",
    "        for j in range(d):\n",
    "            centers[c, j] /= n_c\n",
    "    return centers\n",
    "\n",
    "# eof"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Performance improvement techniques.** This Intrepydd-based implementation tries to do several things to improve the original code.\n",
    "\n",
    "First, this new code is explicit about types. This information allows Intrepydd to generate more specialized (and presumably faster) code.\n",
    "\n",
    "Secondly, it uses a better \"algorithm\" or procedure, in particular, by making fewer passes over the data. Recall the original code:\n",
    "\n",
    "```python\n",
    "    centers = numpy.zeros((k, d))\n",
    "    for j in range(k):\n",
    "        centers[j, :d] = numpy.mean(X[y == j, :], axis=0)\n",
    "    return centers\n",
    "```\n",
    "\n",
    "Notice that it loops over the entire dataset, `X` and `y`, **once per cluster**, or `k=2` times in total. The expression, `X[y == j, :]` extracts all rows of `X` that match the label `j`, which requires one pass over `X`; since `j` ranges from 0 to `k-1=2-1=1`, or two times, then the original loop must make two passes over the data.\n",
    "\n",
    "By contrast, the new version only makes **one** pass over `X` and `y`, in the first loop nest (lines 12-16). The original code is more \"Pythonic,\" but the Intrepydd code effectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Compile using Intrepydd:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!pyddc opt.pydd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Initial test.** The following cell runs the Intrepydd version. You should observe identical (or nearly identical) convergence behavior of the Intrepydd version is correct.\n",
    "\n",
    "> The phrase \"_nearly identical_\" refers to the fact that, due to floating-point roundoff errors, there may be slight differences in the calculation. That's because the Intrepydd version might reordering operations and parallelize as part of trying to speed up your code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iteration 1 WCSS =  565260.5748947266\n",
      "iteration 2 WCSS =  348579.85189060034\n",
      "iteration 3 WCSS =  233326.26372745118\n",
      "iteration 4 WCSS =  226161.20338560524\n",
      "iteration 5 WCSS =  225616.41204630243\n",
      "iteration 6 WCSS =  225564.38665114215\n",
      "iteration 7 WCSS =  225559.57379539014\n",
      "iteration 8 WCSS =  225559.0815381379\n",
      "iteration 9 WCSS =  225559.0249336476\n",
      "iteration 10 WCSS =  225559.02059976387\n",
      "iteration 11 WCSS =  225559.0201870468\n",
      "iteration 12 WCSS =  225559.01980372274\n",
      "iteration 13 WCSS =  225559.01936284977\n",
      "iteration 14 WCSS =  225559.01931455365\n",
      "iteration 15 WCSS =  225559.01927797496\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([0, 1, 0, ..., 1, 1, 1])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import opt\n",
    "update_centers = opt.update_centers\n",
    "kmeans(points, k, starting_centers=points[[0, 187], :], max_steps=50, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Profile run.** Let's re-run the profiler to see if the bottleneck sped up at all."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%lprun -f kmeans kmeans(points, k, starting_centers=points[[0, 187], :], max_steps=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.29 s ± 77.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "777 ms ± 87.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "\n",
      "==> Speedup is ~1.65x\n"
     ]
    }
   ],
   "source": [
    "update_centers, compute_d2 = update_centers__0, compute_d2__0\n",
    "t_0 = %timeit -o kmeans(points, k, starting_centers=points[[0, 187], :], max_steps=50)\n",
    "\n",
    "update_centers, compute_d2 = opt.update_centers, compute_d2__0\n",
    "t_opt = %timeit -o kmeans(points, k, starting_centers=points[[0, 187], :], max_steps=50)\n",
    "\n",
    "print(\"\\n==> Speedup is ~{:.2f}x\".format(t_0.average / t_opt.average))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does this speedup compare to the ideal speedup you estimated above?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Improve the squared-distance computations ##\n",
    "\n",
    "Recall that a second bottleneck is the computation of squared distances. Try replicating the basic development procedure above to see if you can obtain a speedup using Intrepydd. To help you get started, here is the original implementation as a reminder, as well as a placeholder for your Intrepydd code and some timing code. Good luck!\n",
    "\n",
    "```python\n",
    "    # Original implementation, for your reference\n",
    "    def compute_d2__0(X, centers):\n",
    "        return numpy.linalg.norm(X[:, numpy.newaxis, :] - centers, ord=2, axis=2) ** 2\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile opt_d2.pydd\n",
    "# opt_d2.pydd\n",
    "\n",
    "def compute_d2(X: Array(float64, 2), centers: Array(float64, 2)) -> Array(float64, 2):\n",
    "    m = shape(X, 0) # type: int64\n",
    "    d = shape(X, 1) # type: int64\n",
    "    k = shape(centers, 0) # type: int64\n",
    "    \n",
    "    # Place your Intrepydd implementation here\n",
    "    D2 = zeros((m, k))\n",
    "    ###\n",
    "    ### YOUR CODE HERE\n",
    "    ###\n",
    "    return D2\n",
    "\n",
    "# eof"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hint.** If you get stuck, there is [a sample solution](https://github.com/hpcgarage/intrepyddguide/blob/master/notebooks/opt_d2_soln.pydd)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pyddc opt_d2.pydd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test run. Make sure the convergence of `WCSS` looks the same.\n",
    "import opt_d2\n",
    "compute_d2 = opt_d2.compute_d2\n",
    "kmeans(points, k, starting_centers=points[[0, 187], :], max_steps=50, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Timing run using `line_profiler`\n",
    "update_centers, compute_d2 = opt.update_centers, opt_d2.compute_d2\n",
    "%lprun -f kmeans kmeans(points, k, starting_centers=points[[0, 187], :], max_steps=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Testing run using %timeit\n",
    "update_centers, compute_d2 = opt.update_centers, opt_d2.compute_d2\n",
    "t_opt_d2 = %timeit -o kmeans(points, k, starting_centers=points[[0, 187], :], max_steps=50)\n",
    "\n",
    "print(\"\\n==> Speedup is ~{:.2f}x\".format(t_0.average / t_opt_d2.average))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary ##\n",
    "\n",
    "This example is designed to show you how to construct and optimize a larger workflow using a canonical algorithm, the Lloyd's k-means algorithm. Some good practices to observe:\n",
    "\n",
    "- Get a working solution and use profiling to help triage performance bottlenecks.\n",
    "- Use aliasing to help swap different implementations.\n",
    "- When trying to improve performance, look for opportunities to specialize types and reduce the number of passes through the data"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Create Assignment",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}