{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](img/banner.png)\n",
    "%%HTML\n",
    "<script src=\"require.js\"></script>\n",
    "# k-means Clustering\n",
    "\n",
    "**Mahmood Amintoosi, Fall 2024**\n",
    "\n",
    "Computer Science Dept, Ferdowsi University of Mashhad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Definition of Clustering**\n",
    "\n",
    "A clustering of a set of datapoints $\\{x_1, x_2, \\ldots, x_n\\}$ is a partition of the datapoints into k disjoint subsets (or clusters) $\\mathcal{C} = \\{C_1, C_2, \\ldots, C_k\\}$ in such a way that datapoints in the same group are more similar (in some specific sense defined by the analyst) to each other than to those in other groups (clusters).\n",
    "\n",
    "**k-means Clustering**\n",
    "\n",
    "k-means clustering is a type of clustering that partitions the datapoints into k clusters based on their similarity to the centroid of each cluster. The goal of k-means clustering is to find the partition $\\mathcal{C}$ that minimizes the sum of squared distances between each datapoint and the centroid of its assigned cluster.\n",
    "\n",
    "**Lemma 1: Sum of Squared Distances to Any Point**\n",
    "\n",
    "Let $\\{x_1, x_2, \\ldots, x_n\\}$ be a set of datapoints and let $\\mu$ be the centroid of the datapoints. The sum of squared distances of the datapoints to any arbitrary point $z$ equals the sum of squared distances to the centroid plus the number of datapoints times the squared distance from the point $z$ to the centroid. That is,\n",
    "\n",
    "$$\\sum_{i=1}^n |x_i - z|^2 = \\sum_{i=1}^n |x_i - \\mu|^2 + n |\\mu - z|^2$$\n",
    "\n",
    "Proof: \n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\sum_{i=1}^n |x_i - z|^2 &= \\sum_{i=1}^n |x_i - \\mu + \\mu - z|^2 \\\\\n",
    "&= \\sum_{i=1}^n |x_i - \\mu|^2 + 2(\\mu - z) \\cdot \\sum_{i=1}^n (x_i - \\mu) + n |\\mu - z|^2 \\\\\n",
    "&= \\sum_{i=1}^n |x_i - \\mu|^2 + n |\\mu - z|^2\n",
    "\\end{align*}\n",
    "$$\n",
    "since $\\sum_{i=1}^n (x_i - \\mu) = 0$.\n",
    "\n",
    "**Corollary 1: Centroid Minimizes Sum of Squared Distances**\n",
    "\n",
    "The centroid minimizes the sum of squared distances since the second term, $n |\\mu - z|^2$, is always positive.\n",
    "\n",
    "Proof:\n",
    "\n",
    "This follows directly from Lemma 1, since the second term, $n |\\mu - z|^2$, is always positive.\n",
    "\n",
    "**Lemma 2: Sum of Squared Distances Between All Pairs of Points**\n",
    "\n",
    "Let $\\{x_1, x_2, \\ldots, x_n\\}$ be a set of datapoints and let $\\mu$ be the centroid of the datapoints. The sum of squared distances between all pairs of points equals the number of points times the sum of squared distances of the points to the centroid of the points. That is,\n",
    "\n",
    "$$\\sum_{i=1}^n \\sum_{j>i} |x_i - x_j|^2 = n \\sum_{i=1}^n |x_i - \\mu|^2$$\n",
    "\n",
    "Proof:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\sum_{i=1}^n \\sum_{j>i} |x_i - x_j|^2 &= \\frac{1}{2} \\sum_{i=1}^n \\sum_{j=1}^n |x_i - x_j|^2 \\\\\n",
    "&= \\frac{1}{2} \\sum_{j=1}^n (\\sum_{i=1}^n |x_i - x_j|^2) \\\\\n",
    "&= \\frac{1}{2} \\sum_{j=1}^n (\\sum_{i=1}^n |x_i - \\mu|^2 + n |\\mu - x_j|^2) \\tag{Lemma 1}\\\\\n",
    "&= \\frac{1}{2} \\left(\\sum_{j=1}^n (\\sum_{i=1}^n |x_i - \\mu|^2) + n \\sum_{j=1}^n(|\\mu - x_j|^2)\\right)\\\\\n",
    "&= \\frac{1}{2} \\left(n\\sum_{i=1}^n |x_i - \\mu|^2) + n\\sum_{i=1}^n |x_i - \\mu|^2)\\right)\\\\\n",
    "&= n \\sum_{i=1}^n |x_i - \\mu|^2\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "**k-means Clustering Algorithm**\n",
    "\n",
    "The k-means clustering algorithm is a natural algorithm for k-means clustering. The algorithm starts with k initial centroids and iteratively updates the centroids and assigns each datapoint to its nearest centroid until convergence.\n",
    "\n",
    "1. Initialize k initial centroids $\\mu_1, \\ldots, \\mu_k$ randomly.\n",
    "2. For each iteration, perform the following steps:\n",
    "   1. Assign each datapoint $x_i$ to the cluster $C_j$ with the nearest centroid $\\mu_j$.\n",
    "   2. Update the centroids $\\mu_1, \\ldots, \\mu_k$ as the mean of all datapoints assigned to each cluster.\n",
    "3. Repeat step 2 until convergence.\n",
    "\n",
    "**Convergence of k-means Algorithm**\n",
    "\n",
    "The k-means algorithm always converges, but possibly to a local minimum. To show convergence, we argue that the cost of the clustering, the sum of the squares of the distances of each datapoint to its cluster centroid, always improves.\n",
    "\n",
    "Here is a Python implementation of the k-means clustering algorithm:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def kmeans(data, k):\n",
    "    # Initialize centroids randomly\n",
    "    centroids = data[np.random.choice(range(len(data)), size=k, replace=False)]\n",
    "\n",
    "    while True:\n",
    "        # Assign each data point to the closest centroid\n",
    "        labels = np.argmin(np.linalg.norm(data[:, np.newaxis] - centroids, axis=2), axis=1)\n",
    "\n",
    "        # Compute new centroids as the mean of all data points assigned to each centroid\n",
    "        new_centroids = np.array([data[labels == i].mean(axis=0) for i in range(k)])\n",
    "\n",
    "        # Check for convergence\n",
    "        if np.all(centroids == new_centroids):\n",
    "            break\n",
    "\n",
    "        centroids = new_centroids\n",
    "\n",
    "    return centroids, labels\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using OOP for implememntation of k-means\n",
    "\n",
    "\n",
    "This notebook first generates some sample data and defines a KMeansClustering class that implements the K-Means algorithm. The fit method initializes the centroids randomly and then iteratively updates them until convergence. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.cluster import KMeans\n",
    "\n",
    "# Generate some sample data\n",
    "np.random.seed(0)\n",
    "# n_samples = 100\n",
    "# data = np.random.rand(n_samples, 2)\n",
    "\n",
    "# Define the KMeans class\n",
    "class KMeansClustering:\n",
    "    def __init__(self, k):\n",
    "        self.data = None\n",
    "        self.k = k\n",
    "        self.centroids = None\n",
    "        self.labels = None\n",
    "\n",
    "    def fit(self, data):\n",
    "\n",
    "        self.data = data\n",
    "        n_samples = len(self.data)    \n",
    "        # Initialize centroids randomly\n",
    "        self.centroids = self.data[np.random.choice(range(n_samples), size=self.k, replace=False)]\n",
    "\n",
    "        while True:\n",
    "            # Assign each data point to the closest centroid\n",
    "            self.labels = np.argmin(np.linalg.norm(self.data[:, np.newaxis] - self.centroids, axis=2), axis=1)\n",
    "\n",
    "            # Compute new centroids as the mean of all data points assigned to each centroid\n",
    "            new_centroids = np.array([self.data[self.labels == i].mean(axis=0) for i in range(self.k)])\n",
    "\n",
    "            # Check for convergence\n",
    "            if np.all(self.centroids == new_centroids):\n",
    "                break\n",
    "\n",
    "            self.centroids = new_centroids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following script creates an instance of the KMeansClustering class, fits the model to the data, and plots the data with centroids and labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAGdCAYAAABO2DpVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABz9ElEQVR4nO3dd3hUVfrA8e+5d9J76L13pBcpKiiCiCAg2BXL6qrYd/2pu2tdFd21rcqqa28IIlJEEZEqVXrv0juk98zc8/tjkkAkM5kkUzLJ+3mePMncOXPvO5Nk7jvnnvMepbXWCCGEEEJ4gRHoAIQQQghRdUhiIYQQQgivkcRCCCGEEF4jiYUQQgghvEYSCyGEEEJ4jSQWQgghhPAaSSyEEEII4TWSWAghhBDCa2z+PqBlWRw9epSYmBiUUv4+vBBCCCHKQWtNeno69evXxzBc90v4PbE4evQojRo18vdhhRBCCOEFhw4domHDhi7v93tiERMTAzgDi42N9ffhhRBCCFEOaWlpNGrUqOg87orfE4vCyx+xsbGSWAghhBBBprRhDDJ4UwghhBBeI4mFEEIIIbxGEgshhBBCeI0kFkIIIYTwGkkshBBCCOE1klgIIYQQwmsksRBCCCGE10hiIYQQQgiv8XuBLCEqM0trVh85zOH0NOLDwunfuAlhNvk3EUIIT8k7phAFlh48wN8W/MzhtLSibbFhYTx6YT9u6dRFFs0TQggPSGIhBPDbkcPcPnMali6+PS03l2cXLyDfsriza/fABCeEEEFExlgIAby0dDEa0OgS739txVIy8vL8G5QQQgQhSSxEtbc/JZlNJ45j6ZKTCoAcu52f9+72Y1RCCBGcJLEQ1d6prMxS25hKedROCCGqO0ksRLVXJyq61DYOrakTFeOHaIQQIrhJYiGqvcZx8XSrWx/DzayPyJAQBrdo6ceohBAiOEliIQTw94suwVTKZXLxRL+LiQwJ8XNUQggRfCSxEALoWq8+X42+lpaJNYptrxUZyb8GDeHmTl0CE5gQQgQZpbWbofA+kJaWRlxcHKmpqcTGxvrz0EKUSmvN1lMnOZSWSkJ4BD3qN8BmSP4thBCenr+lQJYQ51BK0bF2HTrWrhPoUIQQIijJRzEhhBBCeI0kFkIIIYTwmjInFkeOHOHmm2+mRo0aREREcMEFF7BmzRpfxCaEEEKIIFOmMRbJycn069ePgQMHMmfOHGrVqsXu3btJSEjwVXxCCCGECCJlSixeeeUVGjVqxCeffFK0rVmzZl4PSgghhBDBqUyXQmbNmkWPHj0YO3YstWvXpmvXrnzwwQduH5Obm0taWlqxLyGEEEJUTWVKLH7//XfeffddWrVqxdy5c7n33nt58MEH+eyzz1w+ZsKECcTFxRV9NWrUqMJBCyGEEKJyKlOBrNDQUHr06MHy5cuLtj344IOsXr2aFStWlPiY3NxccnNzi26npaXRqFEjKZAlhBBCBBFPC2SVqceiXr16tG/fvti2du3acfDgQZePCQsLIzY2ttiXEEIIIaqmMiUW/fr1Y+fOncW27dq1iyZNmng1KCGEEEIEpzIlFo888ggrV67kpZdeYs+ePUyaNIn//e9/jB8/3lfxCSGEECKIlCmx6NmzJ9OnT+frr7+mY8eO/POf/+TNN9/kpptu8lV8QgghhAgisrqpEKLayLXbOZSWSohh0jguDqVUoEOqUs5kZTFpy0a+276NtNwcGsfFc+MFnbm6TTtCTTPQ4YkKktVNhRCiQFZ+Pv9ZtZyvt2wiIy8PgMZxcdzbozfXtu8oCYYX/J6cxPXTppCUnY1V8Hk1JecEG0/MZfr2bXx89SjCbSEBjlL4gyxCJoSo0nLtdm6ZPpWP1q8tSioADqam8uT8n3ltxbIARlc1aK2594dZJJ+TVABonD//dvQwb6xc7urhooqRxEIIUaV9uXkjG44fK3bCO9d/16xi15nTfo6qall99Ai7k87gcPEaW1ozafNGsvPz/RyZCARJLIQQVdqXmza4vd9UislbN/snmCpq7bEjmKVcTsrMz2dPcpKfIhKBJImFEKJKO5SWirsR6g6t2Zec7Ld4qiJTGW5f40Jywqke5PcshKjSIkPcDxg0lCImLNRP0VRNfRs1dnmpqVBCeDitatT0U0QikCSxEEJUaSPatHPbTW9pzVWt2vgxoqqnY+069KhX3+XrrIA7unaXKafVhCQWQogq7U9duxNm2jBKOOmZStG+Vi0ubdYiAJFVLROvHEGT+HiAote6MNEY3rot93TvFajQhJ9JgSwhRJW3/thR7v1hFiezMrEZBlprHFrTo34D3r1yBDUiIwMdYpWQY89n9q6dzNy5neScbJrFJ3Bdh070a9RYaoVUAZ6evyWxEEJUC3bLYv6+vWw+cYIQ02BA0+Z0rlM30GEJETSk8qYQQpzDZhgMadGKIS1aBToUIao0GWMhhBBCCK+RHgshhKjGlh86yOcb17P++FFCTJNBzVpwa+euNE9IDHRoIkjJGAshhKimXluxlImrV2EqVVSO21QKQyneHXY1lzZrHuAIRWXi6flbLoUIIUQ1NP/3vUxcvQqg2BofDq2xWxbjf5zF6aysQIUngpgkFkIIUQ19tH5tibU9ADSQb1l8I2uoiHKQxEIIIaqhdcePui3DbWnNmqNH/BiRqCoksRBCiGpIUXrBKsOQolai7CSxEEKIaqhPw0Zu11BRQN+Gjf0XkKgyJLEQQohq6E/dehQbtHkuA4gKDWVM+w7+DUpUCZJYCCFENdS3UWOeunggQLGeC0MpwkNC+GjEKGLDwgMVnghiUiBLCCGqqdu7dOPCho34avNG1h07QqhpY1CzFlzX4QJqRUUFOjwRpCSxEEKIaqxdzVq8MHBQoMMQVYhcChFCCCGE10hiIYQQQgivkcRCCCGEEF4jiYUQQgghvEYSCyGEEEJ4jSQWQgghhPAaSSyEEEII4TWSWAghhBDCaySxEEIIIYTXSOVNISqB7Px8DqalEmqaNI2LR7lZdVIIISozSSyECKCMvDxeX7mMKVs2k23PB6BxXBz39ejN2PYdJcEQQgQdSSyECJCs/HxumDaF7adPYZ2zfPXB1FSemP8zR9LTeOTCfgGMUAghyk7GWAgRIJ9tXHdeUnGut39bye/JSX6OSgghKkYSCyEC5KtNG10mFQCmUnyzdbMfIxJCiIorU2Lx7LPPopQq9tW2bVtfxSZEleWwLI5mpLttY2nN/pQU/wQkhBBeUuYxFh06dOCXX345uwObDNOojizL4vThM1iWplajGpimGeiQgoqhFOE2Gzl2u9s2MWFhfoxKCCEqrsxZgc1mo27dur6IRQQBrTWz/juXqa/O4sSBUwAk1ktg1INXMvYvwzFtkmB4QinFiNZtmbZ9Kw4Xl0McWnNVqzZ+jkwIISqmzGMsdu/eTf369WnevDk33XQTBw8e9EVcohLSWvPWfR/wzgMfFSUVAEnHkvn4b1/xwvVvYFlWACMMLnd370mIaWKUMKXUVIoudepyUZOm/g9MCCEqoEyJRe/evfn000/56aefePfdd9m3bx8XXXQR6emurxXn5uaSlpZW7EsEp02LtzH7/Xkl3qc1LP1uFb9OW+XnqIJX84REvhg1hhoRkQDYDAOzIMm4sGEjPr56dIlJhxBCVGZKazfD0kuRkpJCkyZNeP3117nzzjtLbPPss8/y3HPPnbc9NTWV2NjY8h5aBMCLN77Br9+uxGEvuVfCMA0uuKgdry541r+BBTm7ZbFg3162njpJqGkysGlz2teqHeiwhBCimLS0NOLi4ko9f1do5GV8fDytW7dmz549Lts8+eSTPProo8UCa9SoUUUOKwLkwNbDLpMKAMthcXD7YT9GVDXYDIPBLVoxuEWrQIcihBAVVqHEIiMjg71793LLLbe4bBMWFkaYjGz3mbQz6cz9ZCErvl9DXk4ebXu14qp7BtO0g/eTt6i4SFCAmz6uiJgIrx9XCCFE8ChTYvHXv/6V4cOH06RJE44ePcozzzyDaZrccMMNvopPuLFzzV6eGPJPMlOz0JbzbL97/T5mTvyJe1+/jdEPD/Pq8S65ti9blu1web9hGlx6Q3+vHlMIIURwKdPgzcOHD3PDDTfQpk0brr32WmrUqMHKlSupVauWr+ITLmRnZPPkFS+QlXY2qQCwCi5VvPvop6z7ZZNXj3n5rZdQs0ENDPP8PxvDNIiKjWD4vYO9ekwhhBDBpUw9FpMnT/ZVHKKM5n+1lPTkDJeXJQzTYOqrs+g2qJPXjhkVG8lrC5/lH8MncGjHUWfNCgWOfAc1GyTyz1lPkFg3wWvHE0IIEXykbGaQWjd/E0opXE3qsRwW6xdsRmvt1aW367eoy4db3mDdL5vZsGAz2tJ06NeW3ld1q5LVN/McDn7YtZNvtm3mWEY6daKiGdu+I8NbtyVMqs4KIcR55J0xSFkOy1k8wl0bq9wzid0yDIMegzvTY3Bnn+y/skjPzWXcjG/ZcOI4hlJYWnM4LY3VR4/w2cb1fDV6LLFh4YEOs9rLzs/naHoaESEh1IuO8WoiLYQoO0ksglT7C1uzfOZql8mFYSja9Gopb7IV8Oyi+Ww6eQKgaBXSwu87Tp/ib/Pn8c6VwwMWX3V0OiuLLSdPYCpFs/gEPli/hqnbthStudKuZi0e7N2HITJ1V4iAkcQiSA25fSCfPfsN+Tl5JeYWlqUZ/ZB3Z4VUJ6eyMpm1a4fLZc0dWjNnzy6OpqdRP6ZyF3rLdzg4lpFOiGFSNzo6KJPN1Jwcnl08n9m7dhatrVL4LM79De04fZp7f5jFs5dcyq2du/o9TiGEJBZBK65mLE9/8yjPjv43WuuiwlWGaWA5LEY+MJRLru0b4CiD1/pjR10uDlZIA2uPHa20iUWu3c5/16zii00bSMnJAaBlQiL39ujNqHbtAxyd57Ly87l+2hT2JJ0p9jsp6bejC7b+c8lChrZqTa3IqLP3ac2G48eYum0LR9LTqRkZycg27ejXuImUThfCiySxCGK9h3XnvQ2vMvPtOSybuZr83Hza9GjByAeG0uvKbkH5ybSy8HR0SgUq4vtUrt3ObTOnsfrokWK9LnuTk/jLvDkcSE3h4QuDI/H8Zutmdp057fHvBJy/v+nbt3F3954AOCyL//tlLtN3bMNUCofWmEoxfcc2+jZsxP+GjyIyJMQn8QtR3UhiEeSatGvIg/+9iwf/e5ffjpmWlM7BbYexhdpo2bUZtpCq92fUrW79ogGbriige/0G/guqDKZs3cxvRw6fdzIuvP3WbysY1qoNrWrU8HdoZTZ5a9nrsRhKsS8luej2W7+tYMaObQBFvR6F31ceOczfF8zjjSFXlukYG08cZ8G+veTa7bStWZuhLVvJTCEhkMRClEFaUjrv//VzFnz1K/Z8BwBxNWO49rGrGfOX4RhGmeqtVWq1oqK4qlUbfti9s8RLIqZSDGrekgaV9DLIF5s2uL3fVIrJWzfx1MUD/RNQBRxLTy9Tb0Wh6NBQ7JZFrt3OJxvWudyHpTXf79rB//W9iHoxMaXuNzk7m/t+nMWqI4cxlUIphd2yeG5xOG8PvYr+jZuUI1ohqo6qcyYQPpWZlsWjFz/NL18sKUoqAFJPp/PB418y8cGPAxidbzw/cFDRKqOqYKhg4cWl1jVqMuGyywMUWen2pyS7PRk7tGZvcpLf4qmIwmXly8JuWfy8dzet33mDTu+9TUZentv2ltb8enB/qfu1tObO76ez5ugRwPk62i3n+Ka03FzunPUd20+dLHO8QlQlklgIj8x4ew4Hdxxx1s8owaz/zmXPhn1+jsq3YsPC+GbM9fxr0BC61atH/ZgYutStx4RLL2fatTcQH155F1yLKGW8gKEUMaGhfoqmYsa071iuwZWH0tIAz8fL5FuuV+4ttOzgATYcP1ZiL5ZGY2nN+2tXlyVMIaocSSyER354f16xNUn+yLQZzPlwvh8j8o8wm40x7TsydewNLL39bqZdeyPXdexEuK1yD/Qb3rotppuTsaU1Q1u28WNE5XfTBZ2pGx3t9vnYDMPt/Z64oHadUtv8uGeX2+M4tObHPbvcjs0RoqqTxEKUSmvN6SNn3LZx2C1O7D/lp4hEae7s2p0Q0yzxk76pFK0Sa3B58xYBiKzs4sLDmTrmhvMGyipgSPOWvDjwcv7UtQcP9u5Dy4REDMqWYJhKcUHtOnSqU7fUtpl5eaX2gNgti3yHo5RWQlRdMnhTlEopRVRcFBkpmS7bmDaD2JqlD3wT/tE8IZHPRl7DvT/MIik7G5thoDU4tEX7WrX5cPgoQoJobZd6MTFMvuY6dp05zfrjx7AZBhc2bHTe4Nm3f1uJVYahnqZSxIWHezwjpFlC6Yvs1Y6MktkholqTv37hkUG3XMys/851OcbCYbe49MaL/ByVcKdn/YYsv+PPzN27m80njhNimgxo2owe9RoEbY2T1jVq0rpGTZf3hxhG0WBKV8JNk3zLIi48nDHtO3JHl27Ujor26PjXtr+Ad35b6fJ+Qylu6dzFo30JUVVJYiE8MubR4fzyxRKy0rPPSy4M06BDvzZ0G3RBgKITroSaJsNbt2V467aBDsUvBjVvyY8upggXen7gIMa071iu/TeIjeWJ/hczYekSFMUHhhpK0b5mLW7v0r1c+xaiqpDEIsBysnL58YNf+OF/v3Dy4Clia8QweNwARoy/goTacYEOr0idJrV4bdFz/PO61zm88yiGYaC1RmtNnxE9+L9P769SdSyE9/2enMTnG9czd+8e8i0HF9Suw7jO3bikSVOv9aD8qVsPfti987yTPjgve9SKjOKq1hUbtHpXt57Ui47hndWr2HXmNOCsmXFDx0482KuPVPAU1Z7Sfq5JnJaWRlxcHKmpqcTGVs7iQv6SmZrJXwY+y+8bDzjXOCj4TRimQVytWN5Y8jwNWtYLbJB/oLVm05Jt7F77O7ZQG72GdqV+i9IHvYnqbfH+fdw9ewaW1kW9CYWlte/o0o2/XzTAa8nFnD27eGTuj+Q7LJRyDvJ0aE296Bi+GDWG5gmJXjmO1prjGRnkOuzUi44p97iKU5mZzNi5jaPp6YSbNrrUrUebmjVpGl/6eA4h/MnT87ckFgH02p/e5efPFpU4bsG0GTTv1IT/rvlXACITwnuSs7Pp98n/yLXbXQ6rfGfoVVzZynvTX5Oys5i6bQubT5woGltyRYvKVXJba81/1/zGGyuXOXv//nB/25q1+L++FzGgabOAxCfEH3l6/q48/2XVTHpyBr98ucTtYMjd6/axc/Ue2vRs6efohPCeadu3uk0qDKX4aP1aryYWiRGR/Ll7L6/tzxcmbdnEayuWurx/x+lT3DnrO9664iqGVfDyjRD+JBfFA+T3TQew59ndtlFKsX3lbj9FJIRvrDl2xO39VsFy5pV1pVhfsFsW/1m13KO2f18wj1y7+/cKISoT6bEIENNWeg0Bjca0Se4nyud4Rjpfb9nEgn2/Y7csutdvwE0XdKZdzVp+jcNUpf8Nl6dkdzDbeOIYp7OySm2ngbS8XOb9voerqsnMHhH8JLEIkFbdmhEZG0FWWrbrRhq6Xd7Jf0GJKmPFoYPc+f108hyOovLSe5LO8PXmjTxzyaXc2rmr32Lp26gxc/bscnm/qRR9GjYK2toa5ZGVl+9xW1MpDqam+jAaIbxLPg4HSFhEGKMevNLlm6lhGvQZ3qPSzQoRld+ZrCz+9P10cu2OYmtWOAoGCD67eAG/HTnst3hGtmlHXFi4y14Jh9bc1a2n3+KpDDyp4FnI0pq48HAfRiOEd0liEUC3PD2WS67tA1B0ycMwnd9bdWvO/312f8BiE8Hrm22bybE7OH+egZOpFB+vX+u3eKJCQ/l05DVEh4QWLT9fGAfA3/pfwkVNmp7/wF27oG5d5/cqpmFsHP0bNfboDdhQiiEtWvk8JiG8RS6FBJBpM/nbpIe56p7B/PTRAo7uPU5C3XgG3XwxfYb38GgchhB/tPzQQZdJBTh7CJYdOuDHiKBznbosGHcHU7dtYd7eveQ67HSuW8/9mI9Jk+DECfj6a3jmGb/G6w/PDxzE6G8mkZaTg7si5Hd07U7NyEi/xSVERUliEWBKKTpf0oHOl3QIdCiiivBkckUg5l/EhYXTMCaOUJvJ4fRUfjtyiHrRMdSJiiIxooQT55QpZ79XwcSiaXwCM6+7mf+sWs6MnduLLlsV9ukYSnFH1+78X19Zg0cEFymQJUQV859Vy52rfLr41zaVon/jJnxy9TV+iynf4WD8j9/zy769GEoVxWYoRWJEBJOvua54RcydO6Ft2+K3W7f2W7z+lpWfz9G0NNYcO8LprCziwsO5omUrakVGBTo0IYp4ev6WMRZCVDHXd+iEqQxczbFwltH270JZ7675jfn79gIUS3gsrUnOzubu72cUr2MxbRoULutuGM7bVVhkSAgta9Tg+o6duL/XhdzSqYskFSJoSWIhRBVTJzqaiVdehWkYRQMk4exgyUcv7FfyYEkfyXM4+GzjOpeXXxxa83tKMssOHTy7ccoUKFz+3LLOXhYRQlR6MsZCiCpoUPOW/HTTOD7fuJ75BQWyutWrz7jOXenVoKFfYzmYmkJyTk6xbaH5+XQ8fARVuCCZYXBkzo/QsRMkJcGmTcV3snEjfP89JJawgJhS0K0byJRMISoFGWMhhPCpvUlnuPzLT4ttu3XJUp6dNsP1gwzjbI9FSbf/6K234IEHKhRnMNFas2Df73y+aT1bT50kzLQxpGUrxnXqSpP4+ECHJ6ooGWMhhKgUmsQnnDddcnKf3nx6UT+Akqda/jGJKCmpKLzM88ADcNddxe5Kyclmb9IZkrJLL5sdbLTWPL1oPnfNnsHyQwdJys7mWEY6X2xczxVffeb3qcRC/JH0WAghfO79tb/xyrJfz9t+2eatvPbV10Tm5mFz1yPxRzYbREfD55/D8OFFm/ckneHV5Uv5Zd9eLK1RwCVNmvGXPv3oULuOF55J4H23fSt/nfdTifcpFJEhNpbf8WdiwsL8GteprEx2nDpFiGnSuU5dIkJC/Hp84Xuenr8lsRCiFDn2fKZu28rkLZs4mp5OYkQEY9p34MaOnaXUsocclsVf581h5s4dmErhKDjpA1xgtzNl2kzClq/wfIcXX+wsoNWgQdGm7adOMvbbyeTa7TjOeVszlcJmGHw1+lq61avvpWcUOMMmfc7O06exXAyHVeDX9WCSsrN4dtECftyzq2jGT3RoKLd36caDvfpgGtIxXlVIYiGEF6Tl5nLzd9+w9dRJ4GxhKUMp6kXHMGXMddSPkb9jT2itWXJgP19v2cTe5CTiwsIZ0aYto9t1INo04ZVX4B//cF/hSyl44QV4/PGz01ELjJryFZtPniixfoehFM3iE/j55tuCerGzPIeDthPfdNvGUIphrdrwnyuG+TyetNxcRn/zFQdSUoolc+BMcEa2acerg4cG9WsuzvL0/C2zQnzozLFk0k6nUaN+IrE1YgIdjiiHF5YsZPvpU+d9NrS05nhGOo/M/ZEpY64PSGzBRinFJU2bcUnTZiU3+NOfnIlFae6667ykYteZ02w8cdzlQyyt2ZucxPrjx4K618KT07OCYtOMfemLTevZn5JSYjKngek7t3Nzpy50DeLXXJSd9FH5wNblO/nLwGe4vsHd3N35r4ypcyfPXfNvDu86GujQRBkkZWcxY+f28z6JFXJozeqjR9hx+pSfI6uiZswod7vfk5M9eui+FM/aVVYhpkmPeg1crhQLzr/Lfo2b+CWerzdvclnhFZwJztRtW/wSi6g8JLHwsrXzNvKXAc+w5dftRdu0pVk+aw33936SA9v9t1y1qJitp05i92BA4bpjkjC6k5qTw+cb1/PSr4uYuHoVB1NTSm44derZmR7gHKB57ndwTjv95pvzHhoTFupRLNGhnrWrzO7u3sNtufYaEZEMa+Wf8ucnMjPc3u/QmsPpaX6JRVQeFUosXn75ZZRSPPzww14KJ7g5HA5eveO/WJaFZRX/x7ccFtkZObzzwEcBik6Ulak8+/fwV7dzMPpi0wZ6f/Qezy1ewGcb1/PGymUM/Owjnpz/M/kOx9mGSUmwcOHZaaWGAe3awezZ0KaN8zaAw+Fs94ceip71G5JQykDaCFsIFzVu6sVnFxiDmrfkL336A8X/9hQQExbGZyOvIdzmnxkZ8aW85qZS1CxpgTlRpZU7sVi9ejXvv/8+nTp18mY8QW3dvE2cPpKEtkr+NGE5LDYs2MKxfSf8HJkoj0516hJhcz8MSQF9GjX2T0BB5vtdO3hm0XzyHA40kG9ZWNq5oPs3Wzfzwq+LzjaeNcuZNJxbm2L1ahg2DNasOVv8Silnu1mzih0r1DR5sHcft/Hc17MXkVVkCuT4nr2ZfcMtjG3fkQ61atO9Xn2e7H8JC269g/a1avstjmvad3SbWDu0ZlTb9n6LR1QO5Rq8mZGRwU033cQHH3zACy+84O2YgtaR3cdRSlHaRJtje09Qr1nVmFNflUWHhnLTBZ35aP3aEif2mUpxabPmNI6L93doPpeVn8+snduZu3c3qTk5GIZB7agoGsXGcWXL1nSqU9ftSH+tNa+tWIqi5CXaNfDV5o2M79mb2lHRzssgAPHx8MUXzoSiUHg4vPkmDBoEt97q7K2YOhXGjSu2z1s7dSUjL4//rFqBZWlMQxWNj7m7W0/u69G7Ii9JpdO+Vm1eumxwQGO4o0s3pm3bSkpO9nljkQyluLBBQ7+N9xCVR7kSi/HjxzNs2DAGDRpUamKRm5tLbm5u0e20tKp7vS0qLrLUpKKwnQgOf+17EftSkpm/7/ei+guFy363r1Wbfw26ItAhet3+lGRu+m4qxzLSz7tPofhg3RouatyEiVeOcDlmYeeZ0xxMTXV7HK01837fy00XdIYtW2DAAGdtinr1Sn7AVVc52914o/P7H2NTivE9L+S6Dp34ftcOTmSkUzMyiuGt21InOrrU5y3KrnZUNN+MvZ6Hf/qBzSfP9sQq4KpWbXjpssFuB5qKqqnMicXkyZNZt24dq1ev9qj9hAkTeO6558ocWDC6cHh3QsJs5OfaXbap3bgmrbo392NUoiJCTZP3rxrJkgP7+WbrZg6lpVIzMopr2rVnSItWhPxh2mOwc1gWt82cxkkXg/J0Qf/DskMHeein2Xw0YnSJ7TLz80o9lqEUGXkFHzq2boWoqOKDN0tSv75zjEVmpssmNSMjub1LN5f3p+Xm8N32baw6chiNplf9hlzTroMUOyunZvEJzLz+ZjadOM7mkycIMQz6N24i9V2qsTIlFocOHeKhhx5i3rx5hHv4T/jkk0/y6KOPFt1OS0ujUaNGZYsySMQkRDPm0eF8/fL0kvt/gdv+eT2GVKILKoZSDGjajAGu6i9UIQv3/15qTwM460Is3L+P7adP0a5mrfPubxwbX9Sz44pDa5rFJzhvlKVHQamytT/Hb0cO86fvp5OZdzbxmbd3D6+vXMb/rhpJXxkvU26d6tSlU526gQ5DVAJlSizWrl3LyZMn6dbt7KcBh8PBkiVLeOedd8jNzcX8wye4sLAwwvxcsz6Qbvvn9djz7Hz7xmwATNPAYbewhdn4879v5fJbLglwhEK4tuTAfmyG4dE0W1Mp5uzexc7Tp/hi0wZ2J50hKiSEYa3acluXrlzWrDkL9v1eYh0QBSRGRDKwqWe9d3uTzvDVlk1sOn6cUJvJoGYtGNO+A7FhnvcyHM9I546Z35HjsJ+X92fn27lz1nR+ueV2GkhFYCEqpEwlvdPT0zlwoPjKebfffjtt27bl8ccfp2PHjqXuo7qU9D59NIkl36wg5VQqtRvXYsB1fYmOjwp0WEK49fcF85i6bYtHiYVNKRrExnIgNbVY74SpFGE2G69efgVPLfyFlJycYsmFoRQK+GD4KI96gT7buI7nFy/EUGcHYyogLjycL0aO8XhxsTdWLmPi6lVua0Dc2a0HT/S72KP9CVHd+KSkd0xMzHnJQ1RUFDVq1PAoqahOatZPZPTDvq/VL4Q3dalbj6+3bPKorV1rDhRcNjn3ZO3Qmhy7nacW/sI3Y67n7d9WMnv3zqJk5cKGjXjkwr50r9egxP2ea9mhAzy3eGHRfgtpnOtUjJs5jSW33eXRNNKf9+4p9dLMvL17JLEQooJkrRAhRJGrWrXhxV8XkZGX5/YkDLicSgrORONMdjabT57g9SFX8vzAQZzMzCAuLJwakZ7Pivpg7Zqi2TglHSMpO5tZO7dzfcfS6+nknVuQy2Ub1wOvhRCeqfAowkWLFvHmm296IRQhRKBFhITw3rCrCTFMl4WPDNwnFYVshsH648cAZ02Q5gmJZUoqtNYsO3TA5VothbEsPXjA5f3n6lK3nttiTqZSdKnrYqqrEMJjMj0B5yqknz49mdvaPMh19e/i/y5/nl+/W4XlwXVmIaqaCxs24sebbuWmCzqTEB5BiGEUq0VwUZOmvH/VSI/25WlZdFdK6zXRgF179n96S6cubpMUh9bc0qlrWcITQpSg2l8K2bN+H49d9hxZ6dlYDucbVMqpNNbP38zA6/vx+BcPnDfTRVQ9ezbsY+Y7P7F+/maUUnQbdAFX3z+U5p2qZ9XAZvEJPDvgMp4dcBkA+Q4HKTk5RISEEB0aitaalgmJ7E1OctlzYbcs+leg6qJSigvq1GXLyRNuEgxFVw97GbrUrcejF/bj9ZXLil1eKfz5/p4X0qtBw3LHK4RwqtY9FvZ8O/8YPqFYUgEU/bxwyjJmvDUnUOEJP5nz0Xzu6/448z5fxIkDpzi+/yRzP13IPV0fY+6nCwMdXqUQYprUiooqqrSplOLPPXq5TCpMpWgen8DFTZpW6Lh3dOnmMqlQQKhpMKad5wPH7+91IR8OH0WvBg0xlcJUiu71GvC/q67m0T79KhSrEMKpWvdYrJi1hjNHk1030DDtzdmMeuhKKWpVRf2+6QBv3P0+Wmsc9nNmNtidyeVrf3qXNj1b0rRD1SzqVhGj27ZnT9IZ3l+7GlMZOLRVNPaiTnQ0H189usLlnIe3bsvqo0f4avPG83oZlFK8M3R4mcZtAFzarDmXNmteVH7f3ZonQoiyq9aJxZalOzBDTBz5rkeLnzp0hqRjydRsUMOPkQl/mTXxJwxTFUsqzmUYilkTf+LB/97l58gqP6UUj/e7mCtbtWHS5o3sPHOa6JBQhrVqzfA27byykqhSiucHXMZFjZvw+cb1bD55glDTZHCLVozr3JXWNWpWaN9CCO+r1omFx+8r8gZUZa1bsLmod6IkDrvF+oXnL3glzrqgdh0m+HCVTaUUg1u0YnCLVj47hhDCe6p1/36XSy9w21uhlKJ+y7rUqJfgx6iEP3nyqVXySiGE8Fy1Tix6Du1C/RZ1MMySXwatNWP/MkK6TKuw7pd3xrS5/jcwbAbdB3X2Y0RCCBHcqnViYZomL8x+kvjaccWSh8ITzYj7hjDs7kHnPS47M4dpb8zm9rYPcVXUTVzf8G4+/vskko67GQgqKqWrxw/BbakEDSPGD/FbPEIIEezKtAiZN1TGRcgyUzOZ++kiFn2znKy0bJpd0Iir/jyYThe3P6+3IjM1k78MeIbfNx1Eo4vKDxqmQWyNGN5Y8jwNW9cPwLMQ5bVg0q+8Mu4dUGAVjLcwbAYKeOKLBxlwXfmmIeba7czZs5slB/aRbznoVKcuY9p1JCEiwovRi+oiLTeX7Px8EiMiCJHaOiIAPD1/S2JRRq/f9S5zP11UrO5FIcM0aNqxEe+t+7dcPilBZmom875YwuZft6MUdLq4A5fdfBFRsWWbLugLh3YeYdZ/57J+gXOgZrfLLmDEfUPKnSTuTTrDrTO+5VhGBqZSaJyX1kJNkzevGMYQGYgoPPTbkcO8/dsKlh06CDjLo1/f4QLG97yQuHDPl40XoqIksfCB9OQMrq13F/Y89wsV/Wf5i7S/sLWfogoOm5Zs46kRL5OVnl2UdGmtiYqN5IXvn6Bj/3YBjtB7svPzufTzjzidlXVeCWmFc9nwmdffTPtatQMToAgac/bs4oE5s1EUX93VVIqm8Ql8O/YGSS6E33h6/q7WYyzKav+WQ6UmFcpQ7Fi1208RBYeTh07ztytfIicjBzRoS6Mt52WkrPRsnhz6IqePnAl0mF4za9cOTmRmlrguReGWD9et8W9QIuhk5uXx2LyfnMXb/vC35NCa/SnJ/GfV8gBFJ4RrkliUgbvZA0U0mDa5/nmu79/9mfzcfCyrhBOtpcnLyWf2+/MCEJlvzP99L+4uhDm0Zt7ve/wWjwhOs3fvJCs/32XZdIfWTNm6hRx7vl/jEqI0kliUQcuuzYiOj3LbRmtN98EyPfFcy2f+VuKYlEKWw2LZjN/8GJFv5TrspS4pnudwXT9FCIA9SWewlbKUQLY9nxMZmX6KSAjPSGJRBqHhoYx+eBiuPo4apkGf4T1o2Mqz1Rari7yc0j9R5WVXnU9dHWrVwXQzeNdQirY1a/kxIlFWJzMzmLh6FQ/P/YG/zf+ZRfv3lbqEu7dFhoTgyRA4b5ROF8KbJLEooxv/PppBN10MnL00Ulhgq3WP5vzfZ/cHLLbKqk3Plm4vI5k2gza9WvoxIt+6vuMFbk9Cltbc1rmbHyMSZTF5yyb6ffw/3li5jNm7djJ12xbumPUdIyZ/waks//UODGnRqsRxOoUMpehcpy61otz3ogrhb5JYlJFpmvzfZ/fzxpLnuezmi7ngonb0G9mTZ797jDd/faHUSyXV0dXjryh1PY6rx1/hx4h8q3FcPM8PdBZWO7fnovCnEa3bcnXbqjMLpipZcmA/f1swD4fWWAVfhSf3nadP86dZ0z3qRfCG9rVqM6BJU5crxFpa82DvPn6JRYiykOmmwi8++cfXTHrpOwzTKBpvUfjzzU+NYdxz1wU4Qu9bdugA/1u7mqUHD6CBVok1uKNLN8Z2uKDCy4kL37ju28msPXbUbY/Tl6PG0rdRY7/Ek5GXx/gfZ/HrwQOYykApcFgWNsPkhUsHMbZ9R7/EIQRUszoWWmtOHT6DPc9OrUY1CAmVa46V0dLpq/j29e/ZtnwnKEWHvm0Y8+hw+o3sFejQfMphWVhaS7XESi49N5fO77/jto3NMLj5gs48fcmlforKaeOJ48zZs4vMvDyaJyQyqm074sOlgqvwL0/P30G/bPovXy7h6wnfcXD7EQCiE6IYfs9gbvrHNYRFhAU4OnGu/qN6039UbyyroMeilBHvVYVpGEhKUfnlejhTx9N23tS5Tl0616nr9+MKUR5B/c7+xXNTeeXWtzm040jRtozkTKa8MoMnhrxAXk5eAKMTrhiGUW2SChE8EsLDqRHhvry8w7JoU6OmnyISIjgF7bv7wR1H+Py5bwDOW53SsjRbl+2sUkWXhBC+ZRoGt3Tq4nL8iwLCbDZGtm3v38CECDJBm1jM+XA+hpspjBrNrHfn+jEiIUSwu7t7D7rWrYf6Q7EaUymUUrw++Epiw+QSqxDuBG1icXDHkaIlrkuk4eie436bGiaECH7hthC+GDWGx/r2p250NOCsF3Fps+ZMHXM9V7SUVWmFKE3QDt6MiosoNnWxJOFRYbJ8ufC6I3uOceLAaWJrRNOic9Ny/Y1l5+dzKC2VUNOkSVy8/J1WIuG2EO7p0Ys/d+9Jjt1OiGmWWlpbCHFW0CYWF13Th4VfL3N5v2kzGHhdPz9GJKq6Pev3MfGhj9mydEfRtgat63HXyzd7PGU2PTeXN1YuY8rWzWTbnSvlNo2P5/6eFzK6XQefxC3KRylFhJTLFqLMgjYN7zuiB80uaFxiqWjDUJg2G9c8OjwAkYmqaM/6fTzc/x9sW7Gr2Paju4/x7Oh/s2DSr6XuIzMvj+unTeGLTRuKkgqAAykp/HXeT7zz20qvxy2EEP4WtImFaTN55eenaNW9RcFtA1uIs1pATGIML8/9B43bNijxsacOn+HTpyfz2KDneGLIP5n66izSzqT7LXYRfP77yCfk59nPu/RWOITnnQc+Ii/X/UJqH29Yy84zp89b/6Hw1hsrl3EwNcVLEQshRGAEfeVNrTXbVuxi1Q9rsefZadW9Bf1H93JZfXPx1BVMuPk/aEsXnSSUoQiPCuelH56kY39Zw0EUd2zfCW5tUfrick998ygXj3G9dkPvD99zu4iVqRR/7t6Lv/btX644hRDCl6pN5U1VUBq6Q982pbbdu3E/L9345vmfOi1NTmYOfxv2Ep/veYf4WnG+ClcEoZMHT5faxjANThxw3S7Xbi91ZUxLw/6U5DLHJ4QQlUnQJxZlMf2tH3E1+N6ZXOTy00cLuP6JUf4NTFRqcTVL71mzHBZxNWNc3h9qmoQYJvmW63LQhoKYIKyRsPnkCaZs3cy+5CTiw8MZ1qotlzdvIWujCFFNBe0Yi/JY9cM6t8t3a0uz6sd1foxIBIMm7RvSpEMjt1NCQ8ND6Deyp8v7lVJc1bpNsWXU/8ihNcNald7zVllorXl+8QKunvwlU7ZsYsXhQ8zdu4f753zP1VO+4kxWVqBDFEIEQLVKLBz20hcPsuf7f4EhUbkppbjrlZsLbpTc5qZ/jCEqLsrtfv7cvSc2wyixZLSpFN3r1aefn5bj9oYvNm3g043rAYoGpBYuN777zGnGz/k+YLGJkmmt0VYK2kqW4oHCZ6pVYtGuT+sSp6cWMkzDo7EaovrpfWU3nvrmUeJqOC93KMOZHIRFhHLHizdyw5OlXz5rXaMmn40cQ0LBctc2wyjqwejbqDEfjRgVNIWyHJbFe2t/c32/1vx25DBbTp7wY1TCFa01Outb9Omh6JO90Cd7o08PQWdNlgRDeF21GmMx6oEr+e0H15c6tNZcdc9gP0YkgslF11xInxE9+O3H9Rzff5LYGjH0GdGDqFj3K2Keq1eDhiy/425+2beX7adOEWYzubRpc9rVqu3DyL1vX0oyxzMy3LYxlWLxgf10rF3HT1GJkmit0Wn/hOwvKdbl5jiATnsa8rdA7D+DJqkVlV+1Six6DO7MDU+O4usJ04uVAzdtBg6HxSPv/5mGreoFOEpRmdlCbPS9+uxYivJ82gsxTYa2bM3Qlq29GZpf2S036/QUUEphdzNYVfhJ3sqCpALOVk055+fsbyB8MIRd7O/IRBVVrRILgDtevJEO/doy/T8/sHX5TgzToMeQzlzz8FW07yOXQUTpTh85w7Q3fmDupwvJSM4gsV4CV/5pEKMeupKYhOhAh+cXTePjiQoJITPfdVEwu2XRqU5dP0YlSqKzJgEm4CrJM9FZk1CSWAgvKVOBrHfffZd3332X/fv3A9ChQweefvpphg4d6vEBvV0gy1+yM3PYtWYvlsOiRZemxCa6nlooqq5DO4/wyEVPkZ6cWaweimEa1G1aizeXvkBCnfjABehHL/66iE82rCsasHkuUynqRsewaNydmLKAV0BZpy4HxwH3jYwGGLUX+icgEbR8UiCrYcOGvPzyy7Rq1QqtNZ999hlXX30169evp0OHqrmAkj3fzmfPfMPMd+aQnZEDgC3UxqCbL+ae124tdSaAqFom3PzWeUkFOOtYnDhwirfGf8gz3/41QNH51yMX9mPdsaNsOH4MONvJbipFZEgI7w4bIUlFZaA8eI9SEb6PQ1QbFS7pnZiYyL///W/uvPNOj9oHU4+F1prnx77Gsum/nXct3TANmndqwhu//pPwyMpR1CgvJ49fp61ix2+7MW0mPYZ0odugCzDkzd0rdq3dy/ieT7htowzFpIPvUbN+op+iCqxcu53JWzfx5aaNHEpLJTo0lJFt2nN7l240qOT/39WFzngPnfEm4GpcjIGKvh8VXXrZelG9+bykt8PhYOrUqWRmZtKnj+v1EXJzc8nNzS0WWLBY98smln63qsT7LIfFng37mPvJQq4ef4WfIzvflqXbeWbUv0k7k44ZYoKGaW/MpkmHRrz0w5PUblyr9J3s2gUXXwxLlkDr4B1Y6Cu71/5eahttafZtOlBtEoswm41xnbsxrnO3QIciXIm8FjI/Bp3O+eMsDFDREHFdICITVVSZP8pu3ryZ6OhowsLCuOeee5g+fTrt27d32X7ChAnExcUVfTVq1KhCAfvTnI8XuK17AfDD+/P8FI1rR/Yc44khL5Ce7Jz+58h3FBUDO7zzCI9d9hx5OXml72jSJDhxAr7+2pfhBi1bqGd5eEhYyQvg+UJWfj5ZbgZQCqGMRFTi52DUKNhio+gzpVEDlfgZyvTgg4cQHipzj0WbNm3YsGEDqampfPvtt4wbN47Fixe7TC6efPJJHn300aLbaWlpQZNcHP/9hNsS4Gg4cfCU/wJy4bs3fyA/3462zr+q5bBbHN17gsVTV3D5LZe439GUKWe/P/OMDyINbt0Hd0YZqsTXuVBkbATtLmzl0zi01szYsZ0P1q9mx2nnwmdta9bkrq49Gdm2ndQjEOdRIW2h1kLImYfO+w3QqNCeED4YpUIDHZ6oYsqcWISGhtKyZUsAunfvzurVq/nPf/7D+++/X2L7sLAwwoJwYSWAhDrxxepdlMSTBap8bdE3y7HcJEDKUCz5tpTEYudO2LHD+fP27c7LInI5pJia9RMZdPPFzP9yCVZJyYWC0Q8NIyzCd3/vWmte+NU5G+Pc9GHn6dP8Zd4cNp88zlMXD6zyyYXWuso/R29TKgQirkRFXBnoUCo1nb/DWTRMhUBoX+nNKYcK17GwLKvYGIqqZNAtl7By9lqX9ytDMeT2gX6MqGS5me5ff21pstKy3e9k2jQwTXA4wDCct5980otRVg0P/vcuko+nsObnjc7CanYL02bisDsYPG4ANz89xqfHX3n4EJ9scFaPLaHUEZ9uXM+g5i3pG0RrjnjqTFYWn2xYxzfbNnMmK4uEiAjGtu/IHV27UytSZmeJitH2g+jUv0L+hnO2GujwUai4Z1AqPFChBZ0yJRZPPvkkQ4cOpXHjxqSnpzNp0iQWLVrE3LlzfRVfQPUb2ZM2PVuwe92+83otDJtBYp14hleCEuCN2jZg78b9LrvoTZtB0w6lXH6aMgUKqylalvO2JBbnCY8M46U5f2fT4m3M+3wxKSdTqd24JkNuH0ibni19fvwvN2/AVKpo0a8/MpXiy00bqlxicTQ9jbFTJ3MiM6OobkZSdjYfrlvDjB3bmDr2BhrGxpV5v1pr1h0/yqHUNOLDw+nTsBFhtmpXN7Da045T6KTrwUr+wz0W5ExHWych4UPpJfNQmf6DTp48ya233sqxY8eIi4ujU6dOzJ07l8svv9xX8QWULcTGhJ/+wb9vm8iK79cAoBRoDW26t+BvXz9MbI3AF8oacd8QXr/rPZf3O+wWw2+7CJYvdwb/R0lJsGlT8W0bN8L330NiCbMblIJu3SC8embwSik6D+hA5wH+r92y5eRJl0kFOBf/2nb6pB8j8o8n5//MyXOSikIOrTmdlcX/zZvLpGuuLdM+Vxw6yN8XzmN/SkrRtviwcB7t04+bO3XxQtQiWOisTwuSipKqk1qQ96uzNHqY6xmQ5+1TZ0P2bHTeCkCjQrpDxEiUUfWr85Ypsfjoo498FUelFZMQzfMzH+fInmNsWLAFh92ifZ/WtOzaLNChFRk8bgBLv1vF6rkbivVaFA40vPmpMTRd8TM8+KDrnRjG2R6LwtsjRrhu/9Zb8MADXohelEVkSOkzTiJs/puV4g8HU1P49aDrypEOrVl55BB7k87QIrGGy3bnWn30MONmTjsvUUnJzeHpRfPJtyxu7yJTaKuNrGm4LnkOYKJzZqI8TCx0/iZ00l2gkymcfKlzfoSM1yD+XVTYhRUOuTKTPj8PNWhZjwYtK+cCZabN5LkZ/8e3r33Pd2/9SPLxFACatG/IDU+M4tIbL4KcHOeAzHfeOdvtcq4/LipV0iJThY974AG46y7fPJkKOHHgFFuX70QpRcf+banV0LOTTDAZ2rI1u5POlFhGG8BQKqgXNyvJlpOe9cBsOXXS48RiwtIlWFq7fB1fXf4rY9t3JDrUuzMmtOOo85OvtiC0C8rm+8tnwgM6pZQGDnB4NgNQW0nopNtBZxZsOee9VGejk++Gmj+gbMExO7I8JLGoImwhNq5/YhRjHxtB8vEUTJtJfO24s9cEw8Ph7bfh8sth3DjIyAC7vQwHsEF0NHz+OQwf7psnUU5pZ9J59c53WfH96qJRjMpQXDzmQh55/89Vquz6DR078dH6NWTm5593UjSUIiokhBs6dgpQdL4RanpWbifEMD1qdyAlpagMuSvZdjvz9u5hVDvXNXrKQlsZ6NS/Q+5PnDvsVof0QsW/ijJlsbaAMmqC5S6BNcH08INl1pSCpKKkmXoWkO9c9C328bLHGSSk1nMVY5omNRvUIKFOfMkDjUaMgC1bwE211BL17et8XCVLKnKycvnLgGdY9cPaYtMktKX5ddoq/u/y58nLrToFpGpFRfHFqLEkFIxvMZXCLPg9x4eH8/mosdSKqjqJFECvBo0INd0nDSGGQZ+Gnn0CPJmVUWobUymP2nlCazs6+U7InUvxuTxA/lr0mRvQVqpXjiXKKeJa3J8OHaiIazzalc75Bdfl0537IvfnMgQXfCSxqCYcdgeZqZk4HA5o0AAWLoQXX3Re3nBHKWe7BQucj6tk5n22iP3bDpVYa8RyWOxa8ztLpq4IQGS+06lOXZbefjevXn4Fo9t1YHS7Dvz78itYdvvddK6Cy5THhoVxS6cuuPpLNVBc37ETCRGeLaRVJ6r0wXMOrT1q55Hc+ZC/npJPNg6wjkLWZO8cqwTafhidswCduwytc3x2nGCmom4FswHO5eXPuxfCh0NIFw/35sFrrKtmiYZCklhUccd+P8Hrd73LiNhbGJlwGyPjxzHxwY85fSIV/vQnz3Zy113OGheV0E+fLHR5wgEwDMVPnyzwWzz+EmazMbpdB14ZNIRXBg3hmnYdqvQ0yf/rexFDWznHjpjKKPZ9UPMW/K1/KVVlz9E4Lp5udetjuEmqI2whXN7cO+MfdPZ03L/VanT2t145VrG9Oo5gJd2JPn0ZOuUedPLt6JN90RnvorW7T9TVjzLiUYmTIWwgnPuOoiIg6m5U3CueTzW1XUDJCUohE0I6ViDayq/qvhMJ9m05yCMXPUVOZk5RafKczFxmvTeXxVOX88EDbfFo5v+MGZVysCbAmaPJJc6gLWRZmjNH/jg3XQSbENPk7Suu4o4ux/h22xZOZGZQOyqa0e3a06NegzLXF/jbRZdww7QpoHWJ/QiP97uIKG8N3HScwn3XOGAleedYBbTjJPrMtQX7PfcaYQY64w1wnELFPe3VYwY7ZdZCJfwX7TgO+dudlTdDuqKMsl1aVFE3onO+c9PCgYq8qWLBVnKSWFRRWmteueUtsjNyzrtMYNkt0s6kc/K1d4k7d4aIzeYc0Fn4HZzTTr/5ptImFrUaJpJ0PNllcTDDUNRqVPVmh1RHSim61atPt3r1K7yvbvXq88WosfxjwTz2JJ89qdeIiOSvffpxnTcHwJoNwL4N19MZFRh1vHc8QGd+UJBUuDhm9pfoqBtlVkoJlFkXKjCYVoV0gugH0Rlv4eypKnz/Lfg58lYI7e+FSCsvSSyqqF1rf2fvRtdz/yPzc2ietI+iTzOGAe3awYQJ8PjjzvVCLMtZ4nvhQkhOhoQE/wRfBkPvvIwdv+1xeb9laYbeeZkfIxLBoleDhsy9+TY2nTzBodQU4sMj6N2gISFevuynIq9B5/5USpvrvXY8rS3I/pZS6zJkT0fFPOa144qzVPT9YGuLzvwI8guWhbB1QEXdDuHDqnwFTxljUUXt2+Q6qQDow1FMNLrwD/yBB2D1ahg2DNasOVv8SilncjFrlo8jLp9Bt1xMy67NMEqYkmiYBu37tuGia3oHILLKTWvNiYwMDqWmkudwdwKq2pRSdK5Tl6tat6V/4yZeTyoACL0IQgdAiaOBTLC1gggvrjGjs86poeCyETiOe++Y4jwqfBBGja9Rdbah6mzFqDkNFXFVlU8qQBKLKis0wv314Ys5DIAjKhpmz4Y334TCVWjDw523v/8e4uOd26ZO9VmsFREaHsq/5z/DJdf2LZZcmDaTQbdczMs//R1biHTMneuHXTsZ9vUX9Pn4fS757EN6ffgu/17+K1n5VWdabmWilIFKeAcixwHnrnxrQviVqMSvUEakFw8Y8YfjlNgIDLlE6A9K2Zwry1YjSmt3Q9+8Ly0tjbi4OFJTU4mNDfyS41VV2pl0rqt/F/b8kj+Nfql/5IQZQ9vdqwht1sT1jo4ehRtvhP37nV+V2JljyWxfuQulFO37tiGhdtkXparq/rd2NS8vW4KieEUFQyk61anLpNFjCa9iJcErE22lQ/5GwA62jiizpk+OY6X+A7Ldl6lWNaajQvy/3o0IXp6ev6XHooqKrRHD8HuHuOx2+xOD2fzsO+6TCoD69Z1jLLZs8UGU3lWjXgL9R/Wm38heklSU4HBaKq8sWwKcV6YJS2s2Hj/OZxvX+z+wakQZMaiw/qiwAT5LKgBU1J9BReK6LsMwSSqEz0hiUYXd/e9buPRG5+hj02ZgmAaGzfkrv+LB4dzwt9Ge7UgpZzlvEdSmbN3stnaDRvPlpo1+jEj4irI1ctZlsP1x3RgbRNyMivtXQOIS1YNcfHYjKz2b1T9tIDs9m4Zt6tOhb5ugGnhjC7HxxBcPMvavI/jliyUkn0yhVoMaXD5uAI3bVr4qmsK39iYluVx0q9CR9DTsloXNkM8cwU6FtIIaM8C+BfJ3ggqDsP4oo/LN7hJVS7VMLA7vPkbqqTRqNUykduNa591vWRZfPv8t3/x7JrnZeUXbG7apz2Mf30f7Pm38GW6FtejclBadmwY6DBFgUaGhGErhcJNchBhm0dojIvgppSDkAueXEH5SrRKL9Qs288HjX7J77e9F2zoP6MCfX72VVt2aF2378PEvmfra9+c9/uju4/z1sud4a/mLtOzSzC8xC+EtQ1u2Ztr2rS7vN5Xiylatg6pXTghR+VSb/s7f5qzniSEvsGf9vmLbN/+6nYcveoqda/YCcOrwGb59fXaJ+7AsC0e+g8+enuLzeIXwtkuaNKVDrdol9kgoFIZS3N29ZwAiE0JUJdUisXA4HLxx93toS59X+tlyWNjz7Lxz/4cALPx6Kcpw/YnNclis+mEdaUnpPo1ZCG8zDYNPr76GznXqAWBTRtFYipiwUD4cMYp2Nc+/NCiEEGVRLS6FrJ+/hdNHXC/yYzksdvy2hwPbD5N8IhXDUFhuihFqrUk7nU5sYowPohXCd2pERjJ17PWsP36MBft+J9dhp0Ot2gxt2bpKr44qRHWg7QcgbyWgnQuohQRmPGC1eCc5vu+kZ+1+P0HNBok4HO5XIjRMg3ipkyCClDcX8xJCBJ62UtCpj0PuwuLbQ3qg4l9DmfX8Gk+1uBQSk+hZDYaYGjFcemN/DDeXQgybQf/RvYiOL9tSukIIIYS3aZ2HThoHuUvOvzN/PfrMjWgrza8xVYvEotfQLoRHua+dX7txTdr2aklCnXhueebaEtsYpkF4ZBi3Pe+9lQiF8CWH3cHS6av457Wv8dig53jznv8VDVQWQlQBOT+BfTsll293gHUUsr/xa0jVIrGIiI7g5qfGum1z54SbMAoGst34t9GMf+sO4moWH0PR7sJWvLn0BRq1keJSovJLO5POAxf+jeeueZWl039jw4It/PTxfO7v9QRv3/8hfl4mSAjhAzp7Ou5P5RqdNc1f4QDVZIwFwLWPjcBhd/DlC9+Sn5uPYRpYdouImHDuee02Lr2hf1FbpRQj7x/KVX++nK3LdpKZlkXD1vWlWmU1prUOuvoOL974Jns37gecA5QBHHbn91n/nUv9FnW55pGrAhWeEMIbrNOA+3GBWK4nL/hCtVvdNCMlk6XfrSLlZCq1GtWk36hehEeWtsSwADiw/TA7Vu3GtJl0GdiBmg2q9rLLSceTmfbGD/z0yQLSz6QTXyeeoXdcyuiHhxFXs3KvzLtv8wHu7vxXt20S68Uz6cB7mLaSFqoSQgQDK/m+gkGbrqYyKrC1x6g5vcLH8vT8XW16LApFx0dxxR2XBjoMsjOySU/OJLZGTKVPbE4eOs2/xr3NxkXbirYpQ3HpDf156N27iIiOCGB0vnHs9xM83P8fpJxKK/q0n3w8hcmvzOCXL5fwn2UvVOrEas3cjRiGgWW5/iSTdCyFA9sO07xTKSvcCiEqLRUxFp37i5sWGhV5nd/igWoyxqIyObD9MC9c/zojE27jpib3MiphHC/f8hZH9hwLdGglSktK55GLnmLz0h3FtmtLs3DyMv4x/GUcDjdFP4LUv257h5TTZ5OKQpbD4vTRJN64538Biswz+Xl28ODKjT3f7vtgRNDTjtNY6f/BOnUZ1oleWGfGorO/Q+v8QIcmwi6B0AGUfDo3wNYJIkb5NSRJLPxo97rfub/XE/w6bVXRCcue72DRlGWM7/UEB7YfDnCE55s1cS6nD5/Bsp//yddyWGxavI3Vczb4PzAf2r/1EFuW7ijxOQNYdovfflzHiQOn/ByZ59r0bHFeUvRHYZFhNGojtSyEe9q+B316GGS+C45DoFMgfzM69Ql08l1onVfqPoTvKGWgEt6ByNuAc3uPQyDiGlTipyjl315xSSz8RGvNv2+fSF5O/nlv+A67RXZ6Dm/c/V6AonNt7icLsSzXw3AM0+Dnzxf5LyA/OHeROpc05607U5l0vewC6reog2GW/C9umAZD77i0Sl7GEt6jtYVOvg90GsUHCBb8nLcSnfFOIEIT51AqFCP2CVTtZaiEz1EJn6JqL8OIexFleFbHyZsksfCTnav3sG/zQZefIi2HxdZlOzm444ifI3Mv5VSq2/sth0XS0WQ/ReMftlDPhh552i4QDMPg6W//SmRMRPHkQjnHx7Ts2ozbX7whcAGK4JC3Ahz7cT0w0IKsr6TXopJQRjQq7EJUWF+UER+wOCSx8JOD2z1LGA5WssshNeoluL3fsBnUblLTT9H4R9fLOpY6UyIsIpQLLmrnp4jKp0Xnpry/8VVGPXglcbViCQmz0bBVfe55dRyvL36OyBjprRClyF8HlDJrSKeDvfL23gn/q7wfuaqYiOhwr7bzlyvvGsSHT3513qqwhSy7xRW3B36WjTfF14rjijsv5ccPfinxeSuluPr+oUFxYq7dqCb3vDaOe14bF+hQRAHtOIbOmgS5C0DnQUgXVOTNqNDOgQ6tBB5ORVYyZVmcJT0WftLt8k6ERYS6bRMaEUpert3tFEF/u+qewTRsXb/Ea/XKUPQZ0YOul10QgMh86743bqP3sG4AmDaj4LvzzfOSa/tw+wtS1l2Unc5diT41BDI/APtucByAnNnopLHojMo3xorQPri+DFLAqAlmU39EI4JEtSuQFUifPjWZr16aBi5ecaVAa6jdpBZPfvEAHftXjq721NNpvHXfB/z63aqiT/Ch4SEMv2cwd758EyGhIQGO0De01mxZuoOfP1tE8vEUajZIZPDtA2nXu1XQVeEUgaetFPSpAaCzcfUmoBI+RIVd7Ne43NFao8+MAfs2XCUYKuZxVNSd/g1MBISn529JLPzIsizefeRTZrwzB8Mw0JYucb0GZShsoTbeXvESLTo39X+gLpw+msSedfswbQbt+7YhKjYy0CEJETR05sfo9Fdw+ckCE0IvxEj8xJ9hlUo7jqOTbi0YxKlwxm8CDggfjYp7CaWk87s6kMSiEju27wTTXp/NzIk/uWxjmAZ9r+7JM9+6L8sshAgOVvI9znEVboVg1N3ql3jKQutsyP4BnTMbrBSwNUNFXAehvaX3rhqRkt6VWL1mdYiMjcC0GUWLQv2R5bBYPuM3sjNziIiqXAM6hRDl4NFnuMq54qxSERA5BhU5JtChiCAg/VcBknYmwzmowg3L0mSmZvkpIiGEL6mwXrivs25CaC9/hSOEz5QpsZgwYQI9e/YkJiaG2rVrM3LkSHbu3Omr2Kq0Ok1quZzCWSg0PJTYGjF+ikgI4VMRo4FwXCcXDlTU7X4MSAjfKFNisXjxYsaPH8/KlSuZN28e+fn5DB48mMzMTF/FV2VdfuvFJQ7cLGTaDC6/9RJCw6rmjAshvEnrPHT2bKyUx7FS/orO/AxtpQU6rGKUkYBKeBcIpXh9COfPKvohVNglgQhNCK+q0ODNU6dOUbt2bRYvXszFF3s2RUoGb5711YvT+PSpyedtN0yD+FqxTFzzCjXrJwYgMiGCh7bvRSfdAdYxnCdpXfAVhor/Dyp8YGAD/ANtP4zO/hpyfgHyCwpk3YQK7R7o0IRwyy+DN1NTnetIJCa6Pvnl5uaSm5tbLDDhdOPfRpNQO44vnp/K6SNJQEHRqeE9uO/N2ySpEKIU2spCJ40D60zBlnNrLeSiU8ZDjRmokNaBCK9EytYQFfMYxDwW6FCE8Ily91hYlsWIESNISUlh6dKlLts9++yzPPfcc+dtlx6LsxwOB3vW7SMnM5cGretJQiGEh3TWN+i0f7hpYUL4KIz4l/wWkxBVlc/rWNx7773MmTOHpUuX0rBhQ5ftSuqxaNSokSQWQgiXtH0f5M5HW1mokFYQdhlKnV8S30q6G/IW43aaporFqLPGd8FWIVrnOUuMZ00Fx1EwajqnmEaMdE45FdWaTy+F3H///cyePZslS5a4TSoAwsLCCAsLK89hhKgyLMsiPSkDW6hNKpa6oa0sdOrjkDsX59hyA40dVALEv44K6/eHR7guj312p7nu7xcAaCsDnXwH5G/A+dpbYB1Hp22GrC8g8UuUIb2ponRlmhWiteb+++9n+vTpLFiwgGbNmvkqLiGqBHu+nSn/mslNTe5lTO07GRk/jgcufJLlM1cHOrRKSac+DLnzCm5ZgL3gjhR08l3o/M3FH2DrgPsVOA0IaeP1OKsinfYi5G8quFVYuK8gabPvQ6c84f1jWsno/C1o+z63s+REcCnTpZD77ruPSZMmMXPmTNq0OfvPGhcXR0SEZ91kMitEVBf2fDtPX/0Ka+ZuLPamqQyFtjR/fvVWxjw6PIARVi46f5NzwSuXTAgbiIp/GxyHAI3WdjhzFe56LVTcv1ERV3s73CpFW0nok/0pSuRcUDV/QdkaV/x4jmPotFcKeqYKBtyaLVExD6HCh1R4/9WF1jmQuxgcp8CsDWEDSrxk6C0+GWPhqib8J598wm233ebVwIQIdt+/9zNvjf/A9TlPwac736JBy3p+jauystJehqzPKHWZblUL9Cnnz0YNsF0AeYso6r53NgI0hA9Dxb0mi2SVQucuQiffXWo7FfsyKnJ0xY7lOI4+MxqsZIr/rp2/MxX7Airy2godozrQWZPR6f8CnUHR37uKQ8U+iYqo2O/IFU/P32W+FFLSl6dJhRDVycyJc1BuSjgbhsGPH8z3Y0SVnPZwKnphUgHOaaZ5iyGkK9h6nN1uNkPFPi9Jhcc8XEjMC+uN6fQ3S0gqoDAD12n/RFvpFT9QFaazpqLTni5IKqDo04tORac+gc7+PmCxgSxCJspp6/KdzH7/Z/ZvOURkbAQDru3LoFsuJiJaRo4XOrTjqNvrxpbDYv/Wg36MqHJTZhN0uRbh0pC/HhX/FoR9AlgoJQPGyySkE87TgbtLIQpCeri5v3TayoCc73HfK+WcmULkDRU6VlWldR46/VX3bdL/BeFXopS78Ue+I6m8KBOtNf99+BMe7v8PFn69lD3r97F5yTbeuv9D7mj/CEf3Hg90iJVGWIT7a53KUITLyrVnRYyqwIMNdNbXKBUiSUU5KCOhYC0TV6eEgvEtFR1fYZ0E8ktpZKIdhyp2nKosbwXoZPdtrBOQv9Y/8ZRAEgtRJnM+nM/0t34EKFryXRdUUE46lszfr5qAZZW8FHx10390b0yb638xbWn6jZTVLAspszYq5vHCW2V8tAX2370dUkBp+x509kx09o9oK8nnx1Mxf4OQbgW3Cv9uC34PtpaouAleOIgn4+o0Ssniiy55+rdQVI3W/ySxEB7TWjPl3zNdvudbDovDO4+y9ueN/g2skhr7l+Eow0AZ579gps2gfsu69B/dOwCRVV4q6nZU3OtgNj1naygoD+onqGhfheVX2n4Q68wN6NNXolMfQ6c+jD7ZHyv1OWcBKx9RRiQq8TPn6x/SC8zGENINFfsSqsZUZ69GRY9h1nTu2+2pxwHhwyp8rCrLqO/ddj4giYXw2JljyRzdc9xtPSLTZrLul82uG1QjzS5owvMzH3de7lDO18a0Oa95NmhVj3//8rSsXlsCFXEVquZPqJo/o2rMRNVegYq+D/e9GAYqYoS/QvQZ7TiFTrq+oEjVueyQ/TU65WGf1ntQKgQVcRVGjc8xav2CUeNrVOQYlPLeJTsV82DhTyXdC+GjvTKltcoK7QlGPVz/PygwmxWMmwkMGbwpPGY5PLjEoTxsV030HNKFKUfeZ8GkpexasxdbqI1eV3aj5xVdMAzJ611RSoGtadFtHTEaMj8E6xTnD/wzwYiDyOv8GaJP6KxPXMyYALAg9xfIXwdBvBKqCu0F8RPRqU+ATsV5Gip4vhFjULHPBDK8Sk8pA+KeQyffU7Dl3ERTAQYq9hmX5SH8QRIL4bEa9ROoUT+BM0ddDxxy5Dvo2L+tH6Oq/CKiIxh29+UMu/vyQIcStJQRDYlfopPvBcduzr512cFsiEp4t2qUm876FvczJkx09vSgX2JdhV8GYcucS8c79jkvY4UPRplS08UTKmwAJHyITnsJHHvO3mFrjYr5OyrswoDFBpJYiDIwTZPRDw3jwye+pKTeWMM0iK8dR58RFZuSJkRJlK0x1JwNeSvReasAjQrtAaH9qkStCq016JRSWjnAOu2PcNzSVhbkfI/OXQo4UCGdIXJsmZI7pUIh4krfBVnFqbD+UPMHsO9w/k0YdcDWKqA9FYUksRBlcs0jV7F91W6WfrcKw1BYljPDMEyDiOhw/jnrcWwh8mclfEMpBWF9UGF9Ah2Ks5y4fSfofLC1QBkVm8mglEKrRNDuRv2bzhNIAOn8beik2wumPBqARucugIy3If5NVPiggMZXnSilIKRdoMM4T7mXTS8vKekd/BwOB0umrmTWu3M5uO0wEdHhDLy+HyPGX0GthjUCHZ4QPqW1hqxP0ZkfnNN7EAoRo1Axj6GM8r+vWemvOceSuLkcohKnokI7l/sYFaGtNPSpQQVVUv84lkoBJqrGDFRI6wBEJ3zNJ2uFeIMkFv6ntWb9gi388uViUk+mUatRTa64YyBte7UKdGhCBB0r9Z+Q/UUJ95jOeg+Jk1FGVLn2ra0k9OlRBYWk/phcqHPWPglMd7fO/Ayd/hKup4aZEDEaI+5Ff4Yl/EQSCwFATlYuz476F2vnbcK0GTjsVtH3y8ddwl8+vBfTDEzZVyGCjc7fhj4z0k0LAxX9CCr6z+U/huM4OvUfkPcrZ0/g4RB1Cyr6YZQK3BRlK2mcs/KjOyoBo84q/wQk/MrT87dcDK/i/nPv/1i/YAtwtlJm4fd5ny+mbpPa3PqsrCQohCd09lTAxPWlCgudNblCiYUy66ISP0TbD4J9OxAKoT2dM2MCTed60Ki0kt2iqgv+odTCpVOHzzD/q19d15XQMO3N2eRme/JmIYTAfohSl3W3jnrlUMrWGBU+BBU+sHIkFVBQdMldD6cBIR39FY2opCSxqMLWztuEttxf6cpKy2b7yt1+ikiIIGfE4/7ECgT7Ohe7dkHdus7vf6Air+f8QZvnslCRt/gsNBEcJLGowux57pZAPivfw3ZCVBfacQqd+QlW+r/QmZ+gHacAZ7nx0gpYETHSHyH6zqRJcOIEfP31eXcpW3NUzFMFt85NsApOJRE3QphMN63uJLGowlp1b15qG8M0aN6piR+iEaLy01pjpf8HfepidPorkPkpOv0V9KmLndtD+kNIV0rutTBBRaIib/Nz1F42ZUrx73+gom5GJXwGof1xDtMzIOQCVNzrAS8lLSoHGbxZhbXp0YJW3Zqxd+OBEsdZGDaD/iN7UaNexVctFKJKyPoYMiees+Gc/5vMic6xDgkfolP/CrkLcX42U4DDWVo8/i2UraF/Y/amnTthxw7nz9u3Oy+HtD6/JoUqKFJWOKlQkglxLumxqCCtNRsXb2XG23P48cP5nDnmeh2NQHjyq4eISYjCMIv/qg1DUbdpbe5/508BikyIykXrHHTGRPdtMiaCCsFIeN+5AmvME6joR1EJnzlXY62EVRDLZNo0KJx+bhjO224opSSpEOeROhYVsGvtXibc9B8O7zrmLMerNYZpMOS2Adz/9p2EhocGOkQATh9N4rs3ZjP300WkJ2eQWC+BYXcNYuQDQ4lJqCSjzYUIMJ27CJ18d6ntVML/nItAVUWdO8PmzRQtBtS5M2zYENCQROUhBbJ87PCuo9zX43Fys/POu8ygDEX/Ub14eupfAxSdEKKsdPZsdOqjpbZTca8XDOIMQjk5sG4dJa4imJQEI0acv33WLEgsYXExpaBbNwgP936colKSAlk+9vWE6eTlnJ9UAGhL8+u0Vexau5fW3VsEIDpxrkM7j/Dzp4s4deQM8bXiGHTLxbTs0izQYYnKxtbUw3ZB/LfzwQfw4IOu7zcMsKzit0tKNgq99RY88IDXwtNWpnMdEiMepSK8tl/hX5JYlIPD7mDB10uLKliWxLQZLPjqV0ksAsiyLN579DOmv/UjylBFNT2mvTGbll2b8sq8p4lNDPKaA8J7bB3A1gbsuym5VoMBttZga+/vyLznrrucAzLfeQetFOqPPReW5f42OHsqtHYmFHfd5ZWwdP4udMZbkPsLztfehg4fjop+ILgHw1ZTMnizHHIyc0qtEaE1pJxO81NEoiRTXpnJ9Ld+BDivUNie9fu5teX9JJ9MDURoohJSSqHiXgJCOX86qQmEoeJeDO7BiuHh8Pbb6O/egxgDXdZlgmw2iItzXh556y2vXAbReRvRZ8ZA7nzOJnR2yJmFPjMabd9X4WMI/5LEohzCo8OJiCn9H6pO41p+iEaUJC8njyn/muG2TWZKFq/c+rZ/AhJBQYVcgKoxFcIuwTmNFOf3sEtQNb5BhVwQyPC8QmsL3X8qelFT6BHucp3SEvXtC1u2wPDhXopFo1MfA/I4v/CYA3Q6Ou0ZrxxL+I8kFuVgmiZX3nnZeVM4z2U5LAbfNsB/QYlitizbSWZqVqnt1v68kWP7TvghIhEsVEgbjIT3ULVXomrOQdVe5bwd0ibQoXlH3q/gOAj1DPS0hugnaqBL64RRCl58ERYsgAYNvBdL/lpw7Md1mXAH5K10LsgmgoYkFuV03eMjSayXgGEr+SW89rGradCynp+jEoXysvM8brtj1R4fRiKClTISULYWKCM+0KF4lc5bTdHwOlPBTR7OzrvrrrM1LrzFvte77USlIIlFOSXUieftFS/SZ3gPlHE23Y+rFcs9r43jTy/fFMDoRNOOjTxu667nSYgqb06mZ+1mzPD+sVWUZ+0MD9uJSkFmhVRAzQY1eHbaY5w5lszB7YcJiwildY8W2ELkZQ20uk1r02VgRzYs3OK2nTIVnS4O8mqJQpSBCu2Jzvzf2dvfpzuHkxQMttAmKAdomw1lLxikbhjwzTdemwVSJOwinINl3fQwqgQI6eLd4wqfko9qXlCjXgJdL72A9n3aSFJRifzlo3sJjXBd/VQZikE3XUxCnXj/BSVEoIVeBGYTwIRkByzPRhUMcdAKaBWK/vZBVJs2zoQCwOGAhQsh2btLFigjDiJv5exA2RLaRI9HqcpRxVh4RhILUWXVbVqb/216jZoNa5R4f8d+bXngnTv9HJUQgaWUgUp4D4wE+DnT2TtReF6/Mx696F4Y/TqsWXO2+JVSzuRi1izvxxPzF4i4vuCWSdGKqRgQNR4ib/H6MYVvSUlvUeU5HA5WzFzDjx/9QtKxFGo3qsng2wbQZ3gPTJuXB6MJESS0lQJXXoyauxkdZ0O/Pwh19f9B2IDitTpmz4Zbb3X2Vgwb5rzti3jsByFnJto6gzLqQcTVKLOuT44lykfWChFCCOFekybQvDlMmgT13MxiO3oUbrwR9u93folqSdYKEUII4d7WrRAV5bzU4U79+s4xFpkeziAR1ZokFkIIUV1FR3veVqmytRfVlgzeFEIIIYTXSGIhhBBCCK8pc2KxZMkShg8fTv369VFKMcMX1diEEEIIEZTKnFhkZmbSuXNnJk6c6It4hBCiWtNa4+fJekJ4VZkHbw4dOpShQ4f6IhYhhKi2dO4ydOZHkLcS0OiQbqio21Dhlwc6NCHKRGaFCCGEB7TWkLccnTUZHL+DikdFDIfwESgjsmL7zvwInf4KzsqTDufG/LXolNXoqLsxYv5a4fiF8BefJxa5ubnk5uYW3U5LS/P1IYUQwqu0dqBTn4CcmZw9+St0/mrI/AASv0SZbgpMudt3/raCpAKKkgoAChbwyPwfOrQvKqxv+Z+AEH7k81khEyZMIC4uruirUSPPl7MWQohKIfOjgqQCzp78C8ZBOI6ik8eXe1yEzvoaZ7LiionO+rJc+xYiEHyeWDz55JOkpqYWfR06dMjXhxRCCK/R2o7O+sRNCwfYt0D++vIdIH89xXsqSth//sby7VuIAPD5pZCwsDDCwsJ8fRghhPANx36wzpTSyHQOugztVo4DePL+GFKO/QoRGGVOLDIyMtizZ0/R7X379rFhwwYSExNp3LixV4MTQoiA05aHDT1tV5wKvwydsdXN402QmSEiiJT5UsiaNWvo2rUrXbt2BeDRRx+la9euPP30014PTgghAs7WFFRpKzE7IKQ8vRVAxLWgIij57VgBJirypvLtW4gAKHOPxYABA6R4ixCi2lAqFB15M2S+S9GAzWJMMBtDaJ/y7d+sCQkfoZPvAp3xh3tDUQkTUbam5dq3EIEgdSyEEKIUKvo+dP4WyFuCs2eh8LKFAUY8KuFdVGlLj7vbf2g3qLUIsmeg81YAFiq0J0SMRhkJFX8CQviRJBZCCFEa6wzYmkP+OtCZgAkqDiJvREXdjDISK3wIZcRA1C2oqFsqHm81pXUu5MxB5y4EnQe2dqjIa1Fm3UCHVq1IYiGEEG7o/O3opJtBZ3F2WqgFOhnyfoWoOwMZniig7fvQSbeBdYyiXqXchejMdyH2RVTk6NL3kb8d8paBtkNIFwjtXaGeqOpKEotKbt+Wg+zfcoiwyFC6DOxIZExEoEMSotrQ2kKnPPCHpAKKxlrkb0ZnvImK/Xsgwgta2soAnQ1GAkpV/DSkdV5BUnGyYItV7LtOexJsjZyXl0qMJwmd/BDkr8KZlCjAAWYziH8HFdKqwjFWJ5JYVFIHdxzh1dsnsn3V7qJtYZFhjHnkKm55diym6a5SnxDCK/KWg+OgmwYWZH+Djn6kwuuFVAc6bzU6Y6LzdQVQMejI61BR9zovBZVXztyCngpXDHTmRyUmFs6kZBzYC8sonDPt13EQnXQT1JyNMmuXP75qxueVN0XZHdt3gof7/Z2da/YW256blctXL03jnQc+DlBkQlQz+ZtxX24b5ydvx35/RBPUdM5P6KRbIG/VORvTIfNjdNJ1aKv860jp3IW4P505IHdxyTMac34G+05Krn7qAJ0uJdXLSBKLSujrl74jKz0by1FCwRwNs9/7mUM7j/g/MCGqHZOSp5j+kXT+uqOtDHTK4zhfyz+ewC2w73P2ZJT7APmU/ntyUFIRMp3zPaUmJdkzyh1adSSJRSVjz7fzy1e/4rC7ruJn2gx+/myxH6MSopoK60+pFTWN2mBr4ZdwglbObCAb1yd/h/OSks4r1+5VSAec4yJctgBbK5QqoffJSqHU37GWVbnLQhKLSiYrLZv8nPxS2yUdT/ZDNEJUbyqkPYRciLvLISrqTyWfsEQRbd9Lqb06OhMcJ923cSViDO4vWWlU5K0l32U2LeWxCkxZlbssJLGoZCJjIwgJL33BocS6UjRHCH9QCW+CrXXBrcK3zIITUcRNEDkuAFEFGRWBR5eUyjkAVpk1UXGv4Pz9nJskFPRihA0tSD5KeGzktbhfXRZU5PXliqu6ksSikrGF2Bh000WYNte/GofdYvC4S/wYlRDVlzISUTW+RcW/DWGXQ2hviBiDqvEtRtwzUufAAyp8MO5P3gaEdK1QoTEVcRUqcYrzd0QIzssfbVCxL6HiX3fdqxTSDcJLTjqccXV2mZSIksmIo0roxr9fw9LvVpGZVsIATgXD7xlMozYNAhOcENWQUiEQPgQVPiTQoQQlFdIRHdoPCsqVn0+josdX/DihnVGhbzn3qLVHSZ9SCuJeAFszdNbHziqr4OxlibgWFf0ISnmytL0opLSfVxRLS0sjLi6O1NRUYmNLWzGw+jq44wiv3vlftq/YVbQtPCqMax65iluekToWQojgoq10dMr9BcmFSVERKmyo2OdRkdcENkBA63yw78ZZHKuF1Cb5A0/P35JYVHL7thzkwNZDhEWG0WVgByKipfKmECI4aa0hfxM65yfQmShbC4i4GmXEBzo04QFPz99yKaSSa9axMc06Ng50GEIIUWFKKQjtjArtHOhQhA/J4E0hhBBCeI0kFkIIIYTwGkkshBBCCOE1klgIIYQQwmsksRBCCCGE10hiIYQQQgivkcRCCFGtaZ1X7lU1hRDnkzoWQohqR2sNOT+iMz8B+ybnNlsnVNQdED5U1v8QogIksRBCVDs649+Q+SHFOm3tW9CpD4N9KyrmsUCFJkTQk0shQohqRef9VpBUQPEFsQp+zvwAnbfa32EJUWVIYiGEqFZ05lc4F8FyxSxoI4QoD0kshBDVi30zzlU1XXGAfYu/ohGiypExFkKIaia89CYqzPdhVBE6fzc660vI+xW0BaE9UVG3oEI6BTo0ESDSYyGEqF7Ch+D+rc+AsCH+iiao6ezZ6DPDIfsbcBwG6yjkzEafGYPO/DzQ4YkAkcRCCFGtqMjrQYVT8tufASoCFXmdv8MKOtp+EJ36GM5Br+deWnL+rNNfQOdtCEBkItAksRBCVCvKrINK+BhUdMEWg6K3QhWDSvgYZdbxybG1zkNnz8BKuhXr1BVYSXegc+agtd0nx/MlnTWplBYmOkt6LaojGWMhhKh2VGg3qLUEcr4vmlqqQntB+FUoI9Inx9RWKjppHNi34UxkLHDsR+cthZBekPgBSkX45Ng+kbeKUgfB5q30VzSiEpHEQghRLSkjEiKv89tlD536d7DvLLhlFf+evwad9hIq7p9+icU7PKlOKhVMqyO5FCKEED6mHUcgdx6uP+FbkP0d2krxY1QVFNYP96cQE0L7+ysaUYlIYiGEEL6W9xugS2mUD/kb/BCMd6jI63EWGnPVK2Ghom71Y0SispDEQgghfK60pKKwmYftKgFlNkDFv43zivq5lUydyYaKfREV0iEwwYmAkjEWQgjhayHdPWhkg9DOPg/Fm1T4pVDzJ3T215D7K1BQICvyZpStZaDDEwEiiYUQQviYsjVBh14MecsoeZyFUTAjJdHfoVWYsjVCxfwfxPxfoEMRlUS5LoVMnDiRpk2bEh4eTu/evfntt9+8HZcQQlQpKu5lMBvjHJNQOC6h4LutPSr2qQBFJoR3lTmxmDJlCo8++ijPPPMM69ato3PnzgwZMoSTJ0/6Ij4hhKgSlFkTVeM7VMxTYGsHRm0I6eQci1Dja5QRE+gQhfAKpXXZRgv17t2bnj178s477wBgWRaNGjXigQce4Iknnij18WlpacTFxZGamkpsbGz5ohZCCCGEX3l6/i5Tj0VeXh5r165l0KBBZ3dgGAwaNIgVK1aU+Jjc3FzS0tKKfQkhhBCiaipTYnH69GkcDgd16hSvo1+nTh2OHz9e4mMmTJhAXFxc0VejRo3KH60QQgghKjWf17F48sknSU1NLfo6dOiQrw8phBBCiAAp03TTmjVrYpomJ06cKLb9xIkT1K1bt8THhIWFERYWVv4IhRBCCBE0ytRjERoaSvfu3Zk/f37RNsuymD9/Pn369PF6cEIIIYQILmUukPXoo48ybtw4evToQa9evXjzzTfJzMzk9ttv90V8QgghhAgiZU4srrvuOk6dOsXTTz/N8ePH6dKlCz/99NN5AzqFEEJUXVprdN4yyHjbuRy8zgYVAxGjUFG3ocwGgQ5RBEiZ61hUlNSxEEKI4KatTHTyPZC/qoR7FagoVOKXqJD2fo9N+I5P6lgIIYQQOvVJF0kFgAadiU55AK0tv8YlKgdJLIQQQnhM2w9C7k+ltQLHIchb7peYROUiiYUQQgjP5S7ysKEB+Zt9GYmopCSxEEII4Tmdi8enDhXi01BE5SSJhRBCCM+FtAM8GTthQehFvo5GVEKSWAghhPBcaF8wGwDKTSMFof1QIW38FZWoRCSxEEII4TGlDFTcm4CbpRpsbVDxb/grJFHJSGIhhBCiTFRoZ1TNmRA+BigcR2EDsyXEvYqq8R3KiA9ghCKQylx5UwghhFC2Zqj4l9D6RSAPCEUpd5dHRHUhiYUQQohycyYTsoK1OEsSCyGEECJIacdxyJ6Ftk6jzNoQPhxlBnbtLkkshBBCiCCjtUZnvAqZHxVsMdBYkP4qOuoeVPRDAbs0JYM3hRBCiGCT+R5kfoCzpogF2M/+nPlfyPo4YKFJYiGEEEIEEW1loTPfd98m479oneuniIqTxEIIIYQIJnkrQGe5b6PTIc/VCrS+JYmFEEIIEUx0pnfbeZkkFkIIIUQwsTXzrJ3Z3LdxuCCJhRBCCBFMbB3B1hrXp3ADbJ0CtlaLJBZCCCFEEFFKoeJexlmYzPzDvSaocFTcPwMQmZMkFkIIIbxKW0no7B/Q2TPQ9j2BDqdKUiEdUTWmQtilnD2VGxB2OarGNFRIu4DFJgWyhBBCeIXWuei0FyB7Gs66CgXbQ3qh4l9BmQ0CF1wVpEJaoxImoq10sJLBSEQZ0YEOS3oshBBCVJzWGp38AGRP5dykAoD8tegz16EdZwISW1WnjBiUrXGlSCpAEgshhBDekLcK8hbhrP74Rw6wTqOzPvdzUCIQJLEQQghRYTpnBucPJDyXVdCbIao6SSyEEEJUnOMk4HDfxkr2SygisCSxEEIIUXFmXdz3WABGTb+EIgJLEgshhBAVpiJG4b7HwkBFXuuvcEQASWIhhBCi4kJ6QNgVgCrhThPMehB5i7+jEgEgiYUQQogKU0qh4l+DyNtxVoQsugfCLkYlTkEZ8QGKTviTFMgSQgjhFUqFoGKfQEffD3mrgXwI6SCFsaoZSSyEEEJ4lTKiIXxgoMMQASKXQoQQQgjhNZJYCCGEEMJrJLEQQgghhNdIYiGEEEIIr5HEQgghhBBeI4mFEEIIIbxGEgshhBBCeI0kFkIIIYTwGkkshBBCCOE1fq+8qbUGIC0tzd+HFkIIIUQ5FZ63C8/jrvg9sUhPTwegUaNG/j60EEIIISooPT2duLg4l/crXVrq4WWWZXH06FFiYmJQqqTldcsvLS2NRo0acejQIWJjY72672BR3V+D6v78QV6D6v78QV4Def6+ef5aa9LT06lfvz6G4Xokhd97LAzDoGHDhj49RmxsbLX8YzpXdX8NqvvzB3kNqvvzB3kN5Pl7//m766koJIM3hRBCCOE1klgIIYQQwmuqVGIRFhbGM888Q1hYWKBDCZjq/hpU9+cP8hpU9+cP8hrI8w/s8/f74E0hhBBCVF1VqsdCCCGEEIEliYUQQgghvEYSCyGEEEJ4jSQWQgghhPCaKp1YjBgxgsaNGxMeHk69evW45ZZbOHr0aKDD8ov9+/dz55130qxZMyIiImjRogXPPPMMeXl5gQ7Nb1588UX69u1LZGQk8fHxgQ7HLyZOnEjTpk0JDw+nd+/e/Pbbb4EOyW+WLFnC8OHDqV+/PkopZsyYEeiQ/GrChAn07NmTmJgYateuzciRI9m5c2egw/Krd999l06dOhUVhurTpw9z5swJdFgB8/LLL6OU4uGHH/brcat0YjFw4EC++eYbdu7cybRp09i7dy9jxowJdFh+sWPHDizL4v3332fr1q288cYbvPfee/ztb38LdGh+k5eXx9ixY7n33nsDHYpfTJkyhUcffZRnnnmGdevW0blzZ4YMGcLJkycDHZpfZGZm0rlzZyZOnBjoUAJi8eLFjB8/npUrVzJv3jzy8/MZPHgwmZmZgQ7Nbxo2bMjLL7/M2rVrWbNmDZdeeilXX301W7duDXRofrd69Wref/99OnXq5P+D62pk5syZWiml8/LyAh1KQPzrX//SzZo1C3QYfvfJJ5/ouLi4QIfhc7169dLjx48vuu1wOHT9+vX1hAkTAhhVYAB6+vTpgQ4joE6ePKkBvXjx4kCHElAJCQn6ww8/DHQYfpWenq5btWql582bpy+55BL90EMP+fX4VbrH4lxJSUl89dVX9O3bl5CQkECHExCpqakkJiYGOgzhA3l5eaxdu5ZBgwYVbTMMg0GDBrFixYoARiYCJTU1FaDa/s87HA4mT55MZmYmffr0CXQ4fjV+/HiGDRtW7P3An6p8YvH4448TFRVFjRo1OHjwIDNnzgx0SAGxZ88e3n77bf785z8HOhThA6dPn8bhcFCnTp1i2+vUqcPx48cDFJUIFMuyePjhh+nXrx8dO3YMdDh+tXnzZqKjowkLC+Oee+5h+vTptG/fPtBh+c3kyZNZt24dEyZMCFgMQZdYPPHEEyil3H7t2LGjqP1jjz3G+vXr+fnnnzFNk1tvvRUdxMVGy/r8AY4cOcIVV1zB2LFjueuuuwIUuXeU5/kLUd2MHz+eLVu2MHny5ECH4ndt2rRhw4YNrFq1invvvZdx48axbdu2QIflF4cOHeKhhx7iq6++Ijw8PGBxBF1J71OnTnHmzBm3bZo3b05oaOh52w8fPkyjRo1Yvnx50HaNlfX5Hz16lAEDBnDhhRfy6aefYhhBl0sWU57f/6effsrDDz9MSkqKj6MLnLy8PCIjI/n2228ZOXJk0fZx48aRkpJS7XrqlFJMnz692GtRXdx///3MnDmTJUuW0KxZs0CHE3CDBg2iRYsWvP/++4EOxedmzJjBqFGjME2zaJvD4UAphWEY5ObmFrvPV2w+P4KX1apVi1q1apXrsZZlAZCbm+vNkPyqLM//yJEjDBw4kO7du/PJJ58EfVIBFfv9V2WhoaF0796d+fPnF51MLcti/vz53H///YENTviF1poHHniA6dOns2jRIkkqCliWFdTv+WVx2WWXsXnz5mLbbr/9dtq2bcvjjz/ul6QCgjCx8NSqVatYvXo1/fv3JyEhgb179/LUU0/RokWLoO2tKIsjR44wYMAAmjRpwquvvsqpU6eK7qtbt24AI/OfgwcPkpSUxMGDB3E4HGzYsAGAli1bEh0dHdjgfODRRx9l3Lhx9OjRg169evHmm2+SmZnJ7bffHujQ/CIjI4M9e/YU3d63bx8bNmwgMTGRxo0bBzAy/xg/fjyTJk1i5syZxMTEFI2tiYuLIyIiIsDR+ceTTz7J0KFDady4Menp6UyaNIlFixYxd+7cQIfmFzExMeeNqSkcY+jXsTZ+nYPiR5s2bdIDBw7UiYmJOiwsTDdt2lTfc889+vDhw4EOzS8++eQTDZT4VV2MGzeuxOe/cOHCQIfmM2+//bZu3LixDg0N1b169dIrV64MdEh+s3DhwhJ/3+PGjQt0aH7h6v/9k08+CXRofnPHHXfoJk2a6NDQUF2rVi192WWX6Z9//jnQYQVUIKabBt0YCyGEEEJUXsF/0V0IIYQQlYYkFkIIIYTwGkkshBBCCOE1klgIIYQQwmsksRBCCCGE10hiIYQQQgivkcRCCCGEEF4jiYUQQgghvEYSCyGEEEJ4jSQWQgghhPAaSSyEEEII4TWSWAghhBDCa/4f/C0RNp0teCEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.datasets import make_blobs\n",
    "from scipy.spatial import distance_matrix\n",
    "\n",
    "# Create an instance of the KMeansClustering class\n",
    "kmeans = KMeansClustering(k=3)\n",
    "\n",
    "n_samples = 100\n",
    "X, y = make_blobs(n_samples=n_samples, centers=3, n_features=2,\n",
    "                  random_state=0)\n",
    "# Fit the model to the data\n",
    "kmeans.fit(X)\n",
    "\n",
    "# Plot the data with centroids and labels\n",
    "plt.scatter(X[:, 0], X[:, 1], c=kmeans.labels)\n",
    "plt.scatter(kmeans.centroids[:, 0], kmeans.centroids[:, 1], marker='*', s=200, c='red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's use the Lemma 2 to prove that the centroid minimizes the sum of squared distances\n",
    "\n",
    "$$\\sum_{i=1}^n \\sum_{j>i} |x_i - x_j|^2 = n \\sum_{i=1}^n |x_i - \\mu|^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "function `lemma_2` implements the Lemma 2. The lemma states that the sum of squared distances of all data points is equal to the number of data points times the squared distance from the centroid to the mean of the data points. The notebook tests the lemma on a subset of the data (cluster 0).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1812.8736866847603 1812.8736866847605\n"
     ]
    }
   ],
   "source": [
    "def lemma_2(data, centroid):\n",
    "    n = len(data)\n",
    "    sum_squared_distances_to_centroid = np.sum((data - centroid) ** 2)\n",
    "    # Compute pairwise distances\n",
    "    distances = distance_matrix(data, data)\n",
    "    sum_squared_distances_of_all_points = np.sum(distances**2)/2\n",
    "    print(sum_squared_distances_of_all_points, n * sum_squared_distances_to_centroid)\n",
    "\n",
    "# Test the lemma for  cluster 0\n",
    "centroid = kmeans.centroids[0]\n",
    "lemma_2(X[kmeans.labels == 0], centroid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optional materials related to k-means\n",
    "\n",
    "* [SE Question: Minimizing the sum of distances between points and a point on the plane?](https://math.stackexchange.com/questions/3037564/minimizing-the-sum-of-distances-between-points-and-a-point-on-the-plane)\n",
    "* [Geometric Median](https://en.wikipedia.org/wiki/Geometric_median)\n",
    "* [Weiszfeld's algorithm (Andrew Vázsonyi)](https://en.wikipedia.org/wiki/Andrew_V%C3%A1zsonyi)\n",
    "  - [M. Amintoosi, (2023) Challenges and Requirements of Persian Language Support in Citation Software (In Persian)](https://www.dropbox.com/scl/fi/meqz3eah19a0kxsswxpje/persian-ref-man-challenges-final-draft.pdf?rlkey=odeveabt22we2axrm26hb7utz&st=ofdoo1y0&dl=0)\n",
    "* [Data clustering: 50 years beyond K-means - ScienceDirect](https://www.sciencedirect.com/science/article/abs/pii/S0167865509002323)\n",
    "* [K-Means Factorizaton](https://www.dropbox.com/s/2fsmpj7i7x3rngy/K-Means%20Factorizaton1512.07548.pdf?dl=1)  \n",
    "* [Balanced K-Means for Clustering | SpringerLink](https://link.springer.com/chapter/10.1007/978-3-662-44415-3_4)\n",
    "* [M. Amintoosi, (2002) Student’s sectioning using fuzzy inference system (In Persian)](https://www.dropbox.com/s/7ohrff9c6im8bye/1382-StudentSectioning.pdf?dl=1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pth",
   "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.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}