{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TP1 KMEANS\n",
    "\n",
    "On nous propose de coder l'algorithme des kmeans afin de faire du clustering sur 2 classes puis plus de 2 classes.\n",
    "Plus tard, on utilisera notre algorithme pour segmenter une image sur l'information de couleur."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 379,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.spatial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 380,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mean = [1,2,3,4]\n",
    "# sd = [0.25, 0.25, 0.1, 0.2]\n",
    "clusters = 2\n",
    "mean = np.random.randint(5, size=clusters)\n",
    "sd = [0.25, 0.25, 0.3]\n",
    "dim = 2\n",
    "nb = 50\n",
    "K= clusters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fonctions à utiliser pour le clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 381,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_points(mean=1,sd=0.5, nb=100, dim=2, clusters=2):\n",
    "    size = []\n",
    "    # for i in range(0,dim):\n",
    "    size.append(nb)\n",
    "    size.append(dim)\n",
    "    points = np.random.normal(mean[0],sd[0],size=size)\n",
    "    for i in range(1,clusters):\n",
    "        points = np.concatenate((points,np.random.normal(mean[i],sd[i],size=size)),axis=0)\n",
    "    \n",
    "    return points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 382,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance(points,Pc):   \n",
    "    return scipy.spatial.distance.cdist(points[:,:], Pc[:,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 383,
   "metadata": {},
   "outputs": [],
   "source": [
    "def kmeans(points = [0,0], K = 1, nb=1, dim=2):\n",
    "    # Initialisation K prototypes\n",
    "    Pc_index = []\n",
    "    Pc_save = np.zeros([K,dim])\n",
    "    clusters = []\n",
    "    iter = 0\n",
    "    eps = 0.1\n",
    "\n",
    "    for i in range(0,K):\n",
    "        Pc_index.append(np.random.randint(0,nb*dim))\n",
    "    Pc = points[Pc_index,:]\n",
    "\n",
    "    while (np.mean(distance(Pc,Pc_save)) > eps and iter < 10):\n",
    "        iter += 1\n",
    "        Pc_save = Pc\n",
    "        # print(Pc)\n",
    "        # print(points[:,:Pc.shape[0]])\n",
    "        dist = distance(points=points[:,:Pc.shape[0]],Pc=Pc)\n",
    "        clust = np.argmin(dist, axis=1)\n",
    "        clust = np.expand_dims(clust, axis=0)\n",
    "        points = np.append(points[:,:Pc.shape[0]], clust.T, axis=1)\n",
    "        # print(points)\n",
    "        Pc = np.zeros([K,dim])\n",
    "        index = np.array([])\n",
    "\n",
    "        for n in range(0,2*nb):\n",
    "            for k in range(0,K):\n",
    "                index = np.append(index, (clust==k).sum())\n",
    "                if points[n,-1] == k:\n",
    "                    # print(points)\n",
    "                    # print(Pc)\n",
    "                    Pc[k,:] = np.add(Pc[k,:], points[n,:-1])\n",
    "\n",
    "        for k in range(0,K):\n",
    "            Pc[k,:] = np.divide(Pc[k,:],index[k])\n",
    "\n",
    "        # print(Pc)\n",
    "    return Pc, points\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 384,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors=['red', 'green','yellow','blue','purple', 'orange']\n",
    "def visualisation(points, Pc=[0,0], dim=2, K=1):\n",
    "    if(dim==2):\n",
    "        for k in range(0,K):\n",
    "            for n in range(0,len(points)):\n",
    "                plt.plot(points[n,0], points[n,1], 'o', color=colors[int(points[n,-1])])\n",
    "        plt.plot(Pc[:,0],Pc[:,1],'r+')\n",
    "        plt.grid(True)\n",
    "        plt.axis([min(mean)-1,max(mean)+1,min(mean)-1,max(mean)+1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 385,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = gen_points(mean,sd,nb,dim,clusters)\n",
    "# print(points.shape)\n",
    "# print(points.mean(axis=0))\n",
    "# print(points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 386,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist = distance(points,points)\n",
    "# print(dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 387,
   "metadata": {},
   "outputs": [],
   "source": [
    "Pc, clusters = kmeans(points,K=K,nb=nb,dim=dim)\n",
    "# print(Pc)\n",
    "# print(clusters)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 388,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2.00659379  2.0037594 ]\n",
      " [-0.05586229 -0.02372516]]\n",
      "[0 2]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGiCAYAAADulWxzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAABGdklEQVR4nO3dfVhUdcI//vcwMoOkaJgCKj5vPpXPiuheKquE6Tfj/u66rfVdtTVbVzFwVMK92tSyMFRATdf6eif3tnnX1pXUtzWNEPQ2EfOBX+qWJZq0BViakJAzI3N+f9BMDMzDOcM5M2fOvF/XxWVz+JyZz2cGlvd+HnWCIAggIiIi0qCwQFeAiIiISCkMOkRERKRZDDpERESkWQw6REREpFkMOkRERKRZDDpERESkWQw6REREpFkMOkRERKRZDDpERESkWQw6REREpFmKBp2//vWvGDFiBKKiohAVFYXExES8//77Hu958803MWTIEERERODee+/F/v37lawiERERaZiiQad3797YuHEjTp06hZMnT+JXv/oVHnzwQZw/f95l+WPHjmHevHlYtGgRzpw5g9TUVKSmpuLcuXNKVpOIiIg0SufvQz2jo6OxadMmLFq0qM33HnroITQ0NOC9995zXJs4cSJGjRqFXbt2+bOaREREpAEd/PVCTU1NePPNN9HQ0IDExESXZcrKymAymZyupaSkoLCw0O3zms1mmM1mx2ObzYbr16+jW7du0Ol0stSdiIiIlCUIAn744Qf07NkTYWHyDTgpHnTOnj2LxMRE3Lp1C506dcK+ffswbNgwl2VramoQExPjdC0mJgY1NTVunz87Oxvr16+Xtc5EREQUGF999RV69+4t2/MpHnQGDx6MiooK1NXV4a233sKCBQtw+PBht2FHqjVr1jj1AtXV1aFPnz74/PPPER0dLctrBAOr1YqSkhIkJSUhPDw80NXxG7ab7Q4FbDfbHQquX7+Ou+++G507d5b1eRUPOgaDAYMGDQIAjB07Fh9//DG2bt2Kl156qU3Z2NhY1NbWOl2rra1FbGys2+c3Go0wGo1trkdHR6Nbt27trH3wsFqtiIyMRLdu3ULqF4PtZrtDAdvNdocSuaed+H0fHZvN5jSnpqXExEQUFxc7XSsqKnI7p4eIiIjIE0V7dNasWYP7778fffr0wQ8//IC9e/eitLQUBw8eBADMnz8fvXr1QnZ2NgAgPT0dU6dOxZYtWzB79my8/vrrOHnyJF5++WUlq0lEREQapWjQuXr1KubPn4/q6mp06dIFI0aMwMGDB5GcnAwAqKqqcppZPWnSJOzduxdPPfUU/vznP+MXv/gFCgsLcc899yhZTSIiItIoRYPOf/7nf3r8fmlpaZtrc+fOxdy5cxWqEREREYUSnnVFREREmsWgQ0RERJrFoENERESaxaBDREREmsWgQ0RERJrFoENERESaxaBDREREmsWgQ0RERJrFoENERESaxaBDREREmsWgQ0RERJrFoENERESaxaBDREREmsWgQ0RERJrFoENERESa1SHQFSAiIj+zWICdO4HKSmDgQGDpUsBgCHStiBTBoENEFEoyM4HcXKCp6edrq1YBJhOQkxO4ehEphEGHiChUZGYCmza1vd7U9PN1hh3SGM7RISIKBRZLc0+OJ7m5zeWINIRBh4goFOzc6Txc5UpTU3M5Ig1h0CEiCgWVlfKWIwoSDDpERKFg4EB5yxEFCQYdIqJQsHQpoNd7LqPXN5cj0hAGHSKiUGAwNC8h98Rk4n46pDlcXk5EFCrsS8db76Oj13MfHdIsBh0iolCSkwNs2MCdkSlkMOgQEYUagwHIyAh0LYj8gnN0iIiISLMYdIiIiEizGHSIiIhIsxh0iIiISLMYdIiIiEizGHSIiIhIsxh0iIiISLMYdIiIiEizFA062dnZGD9+PDp37owePXogNTUVFy5c8HhPQUEBdDqd01dERISS1SQiIiKNUjToHD58GMuWLcPx48dRVFQEq9WK++67Dw0NDR7vi4qKQnV1tePrypUrSlaTiIiINErRIyAOHDjg9LigoAA9evTAqVOnMGXKFLf36XQ6xMbGKlk1IiIiCgF+Peuqrq4OABAdHe2x3M2bN9G3b1/YbDaMGTMGzz//PIYPH+6yrNlshtlsdjyur68HAFitVlitVplqrn72toZSmwG2m+0ODWw32x0KlGqvThAEQZFnbsVms2HOnDm4ceMGjh496rZcWVkZvvjiC4wYMQJ1dXXYvHkzjhw5gvPnz6N3795tyq9btw7r169vc33v3r2IjIyUtQ1ERESkjMbGRjz88MOoq6tDVFSUbM/rt6Dzpz/9Ce+//z6OHj3qMrC4Y7VaMXToUMybNw/PPvtsm++76tGJj49HdXU1unXrJkvdg4HVakVRURGSk5MRHh4e6Or4DdvNdocCtpvtDgXXrl1DXFyc7EHHL0NXaWlpeO+993DkyBFJIQcAwsPDMXr0aFy8eNHl941GI4xGo8v7QukHxI7tDi1sd2hhu0NLqLVbqbYquupKEASkpaVh3759OHToEPr37y/5OZqamnD27FnExcUpUEMiIiLSMkV7dJYtW4a9e/finXfeQefOnVFTUwMA6NKlCzp27AgAmD9/Pnr16oXs7GwAwDPPPIOJEydi0KBBuHHjBjZt2oQrV67gscceU7KqREREpEGKBp2//vWvAIBp06Y5Xd+zZw8WLlwIAKiqqkJY2M8dS99//z0WL16Mmpoa3HnnnRg7diyOHTuGYcOGKVlVIiIi0iBFg46Yec6lpaVOj/Py8pCXl6dQjYiISJMsFmDnTqCyEhg4EFi6FDAYAl0rUgG/7qNDREQku8xMIDcXaGr6+dqqVYDJBOTkBK5epAoMOkREFLwyM4FNm9peb2r6+TrDTkjj6eVERBScLJbmnhxPcnOby1HIYtAhIiJ1s1iA/Hxg+fLmf+3BZedO5+EqV5qamstRyOLQFRERqdawggJ0+N//G7DZfr5on3/z44/inqSyUpnKUVBg0CEiIlUKy8rCoMLCtt+wz79JShL3RAMHylovB670CgoMOkREpD4WC8Ly8wEAOndlDh8GwsKce3ta0+ubA0iL55UlnHClV9DgHB0iIlKfnTuhs9nchxygOeBMner5eUymn4NMZiYQGQmsWAG8+GLzv5GRzdelsK/0aj0/yN7TJPX5SFEMOkREwcbd5FwtETuvZvhwYPXq5p6blvT65uv23hW5wglXegUdDl0REQUTLQ+ZtBxW+ve/xd0zcCCQkQFs2OB+SEpsONmwwfUwVut6iV3plZEhrg2kKAYdIqIgEZaV5foPtlo3x5MyH8ZFgLMfIuR2+Krl/BuDwX2wkLIMvfVzuAqWYnCll2pw6IqIKBi0mJzrlpqGTKTMh3EzrGQPOG5PTZwyRVxdxIaO1uXcDXeJodRKL5KMQYeIKAgMOHAAOk+riwD1bI4nZT6MmGEld0pKxE0mFhs6WpZrT71ar/SigGLQISIKAnfU1IgrGOghE6mTdb0MK+ngYegKEDeZeOnStpOVW2sdTsQMd7nTcqUXBRyDDhFREGiIjRVXMNBDJlKPZZArmHkatjMYmsOHJ63DiS/10umcV3qRKjDoEBEFgUszZ0II8/I/2WoYMpE6H0auYOZt2C4nR9wydDtf6qXTNa/cIlVh0CEiCgYGA2zeliurYchE6nyYpUubA4IcvIWsnBygsRHIywPS0pr/bWx03QMjZrirNZtNHXOkyAmDDhFRkLBt3CitVyIQpM6HMRiaV2R54HbVVWtiQpZ9Gfr27c3/uguGYoa7XAn0HClqg/voEBEFk5wcz5vjBZo9INj39XElPd25/tnZwP/8D/Dxx22KOkKOTgcIHiKPEsN29uAoZR+dQM+RojYYdIiIgo2nzfHUwF1A0OuBMWOArVtd7+z8y182H2nRItAIYWEQMjKg1+s9hyelhu3swXLrVu/L2NUwR4raYNAhIiL5uep5qqpqnhfTmn2J+OrVwK1bjnua+vXDe336YFZqKvTh4c1lXYUnpY6/aL2z84oVrutvp4Y5UtQGgw4RESmjZc+TxdK8uZ8n9vOmfrrHZrUC+/f//H1/Dtu5OvpBrwfGjwdOn/Zf2KJ2Y9AhIiJpxJxh1bqM1SrPYZj+GLaz7+zcWlNT8zwikwmIj1fnHClqg0GHiIjEH8Ap5vR0V2XELiEP9KolMTs75+UB9fVAp07Sn9ufk8j9/XoqxaBDRBTqxIQXezl3PR0tr7sq42nFVEuBXrUkZmdnQQCioprfI7HDVWLfY7n4+/VUjEGHiCiUiQkvOTniz7DydvCoJ2pYtfT55+LKCYLz++OJ2PdYLv5+PZXjhoFERKFKygGcYs+wEttz48qYMYEfWqmullbe0xlbgPRDTtvLYgG2bPHf6wUBBh0iolAl5QBOf8yd+fhj73vVKE3s4al23s7YknrIaXvNnOm9V03O1wsCDDpERKFKygGc/po7E+jehsGDpd/j6X2Ueshpe2RmAiUl/nu9IMGgQ0QUqqQcwCn2DKv2HtAZ6N4GXw7z9PQ+Sj3k1FdihsjkfL0gwqBDRBSqpBzAKeaQS5OpeWVPewWyt0HqYZ7eJlBLPeTUV2KGyOR8vSDCoENEFKrEhhf7BOGcHO+np3sqk5Qkrl6B7m1w1wZXvB37IPU99pWUcBhiR1Uw6BARhTIx4aV1+cbG5k3z0tKa/21sdC7nrsyBA/7p3ZBDyzb06OG6zPjx4pZpS32PfdGnj7hyU6eG1NJygPvoEBGR1DOkxBzD0LqMfYn6iBHAmTPu71NTb4PBAHzzDXD1quvv21eJiQ07/jqny5PZs/37eirAoENERMqeIeVql97W1Hgwptg9cDZsEBdYlHyPq6rkLachig5dZWdnY/z48ejcuTN69OiB1NRUXLhwwet9b775JoYMGYKIiAjce++92N/y9FoiIgoe9l163YWcUaNcD3/9xFhbC/3s2cAf/gDk5/t36bm/98BpD3+t7gpCigadw4cPY9myZTh+/DiKiopgtVpx3333oaGhwe09x44dw7x587Bo0SKcOXMGqampSE1Nxblz55SsKhERyU1Mj8jZs26HcMKysnDfkiUIKyoC9uwBVqwAIiOlbSposTQHpOXLpQclf+6B017+Wt0VhBQNOgcOHMDChQsxfPhwjBw5EgUFBaiqqsKpU6fc3rN161bMnDkTq1evxtChQ/Hss89izJgxePHFF5WsKhERya09PSKZmQjLzYWu9ZES9vOaxISdzMzmYLRiBfDii9KDUjD1kvhrdVcQ8uscnbq6OgBAdHS02zJlZWUwtfqwUlJSUFhY6LK82WyG2Wx2PK6vrwcAWK1WWK3WdtY4eNjbGkptBthutjs0BGu7wz7/HGK23mv6/HPYWratqgodNm8GALjaflAAgNxc3F671u0f7rCsLIT91JvU8jmEn4KSrakJto0bPVds8WJ0WLkSsNnc10Ovx+3FiwEZPxufP+/nnkNYUxPC8vOha3EMhKDXw5aeDttzz8laT7kp9fPtt6Bjs9mQkZGByZMn45577nFbrqamBjExMU7XYmJiUFNT47J8dnY21q9f3+Z6SUkJIiMj21fpIFRUVBToKgQE2x1a2O7gMMBsxr0iyv3LbMalFnMxJ65bhxgPh4PqAKCpCZ898QQuzZnTtoDFggfy838u2+peAYAuPx/7J0702sMxbM4cDCosbL6nxXV77S4+8AD+9eGHHp/DVz593lOmABMnYsCBA7ijpgYNsbG4NHNmcztVPt+1sbFRkef1W9BZtmwZzp07h6NHj8r6vGvWrHHqAaqvr0d8fDySkpLQrVs3WV9LzaxWK4qKipCcnIzw8PBAV8dv2G62OxQEbbtnzIBQUOC1R2TItm0Y0iJwhL39NlBR4fXphxmNGDJrVpvrYdu2IczDwZY6ADqbDf+rqgq2J57w/CKzZsGWlYWw/HznwzJ/6iXpt3Ej+nmtqTSyfN6pqY7/HCJPtRR37do1RZ7XL0EnLS0N7733Ho4cOYLevXt7LBsbG4va2lqna7W1tYh1c6Ks0WiE0Whscz08PDy4/gdBJmx3aGG7Q0vQtTs8HFi5snlOjQs6ADCZEH7HHc7fGD0a+NvfvD69/u67oXf1fnz5pajq6b/80vX9rW3ZAmRnO+2Bo1u6FHqDQdTQnK+C7vNuJ6XaqmjQEQQBy5cvx759+1BaWor+/ft7vScxMRHFxcXIaLHXQFFRERITExWsKRERSWLfANDb5nf2JeOt99HxtG/O0qXAqlUQmppc9gQ57ne3gkiJScRK7oFDilI06Cxbtgx79+7FO++8g86dOzvm2XTp0gUdO3YEAMyfPx+9evVCdnY2ACA9PR1Tp07Fli1bMHv2bLz++us4efIkXn75ZSWrSkQUvMSGDqll3XG1AeCqVe6Diy87L5tMwKZNbebGOHhaQSRmU7wQXWodkgQFoXkIts3Xnj17HGWmTp0qLFiwwOm+f/zjH8Ldd98tGAwGYfjw4cI///lP0a9ZV1cnABC+++47mVoRHCwWi1BYWChYLJZAV8Wv2G62OxR4bPfq1YKg1wsC8POXXt98vT1l3Vm92vn+1l9SnsuL2yaT0BQWJq2+3uqnQD3lFqo/5999950AQKirq5P1eRUfuvKmtLS0zbW5c+di7ty5CtSIiEhD7LsOt2bfawb4uYdFSll35D4SwQvbxo3YP3Ei/ldVFfRffimut8pb/QD1HTVBiuJZV0REwcBiQdi2bc0TbQcOBB57THzosP+3mLKeAorYDQAfegjo3VuegysNBtieeELcpGEx9QOA+Hjf60NBh0GHiEjNLBaE3Xcf5hw+7DxXxWRqHoTxpOWuw2J3KPY04VbsUQctN3j1NHdHbsF0ZAP5DYMOEZFaZWYCmzdD7yrQiJgaAEDaH3VvZX056sDb0Jgck6Ol1s9dOTnrQqqh6FlXRETkI/ucGrGBxp2BA+Vbbv3YY77XIze37YGa7T2LqrX2HGwpd11INRh0iIjURuykWm/sf9TlONnaYgF+/3vf69L68E57kGs9pCbl0M7WfD3YUom6kGow6BARqY3YSbXe2P+ot/dka3tvh5vDlUWzD42JXb3VugdIjJwcYPXqtsFOr2++3nr4TMm6kCow6BARqY3UybK6VlvqufqjLjUA2Lnr7fCFfWhM5OqtsF27fHudnBygsRHIywPS0pr/bWx03UaxK8la9kZRUOFkZCIitZEy6VevB27cAHbvFnccg5QdiuUaQrPXc+nS5ud86y1x91RWAoMG+fZ6Yo9s4EotzWPQISJSm5/OehLVi2IyAZ06iT+HScqZTXINoQHN9XzqqbZHR3jiyyovqZQ4F4tUhUNXRERqI2ZOjU7necgJaO49yc8Hli9v/lfqPBMpvRh6PTB+vPuhMUDaEJheD9uSJeJf31dyTNQmVWOPDhGRGrk79VunA6ZNAw4c8LzHi9SDN10R24uRmgq88UZzfW7ebF6dZR8ae/XV5uuRkeKey87T5Gg5tThANOB1IUWwR4eISK1+mlTbtHkzLs2ahabNm4Fbt4BDh7yHHDmWS4vt7bCHnMxMoGvX5tVZZ882/9u1KzBzpqSeHK89VXLzdaI2BQUGHSIiNfvprKezjz8O2xNPeO9ZkHO5tJRl6Z7CVUmJ99cCgMmT3a+OUpqUlVoUVDh0RUSkJVKWS4uZlOxuCE2v/3kYTK7VWb/5TWCHiKRM1KagwaBDRKQlSiyX9rYsXY7VWZzwSwph0CEi0hKllkt76u2QY48ZTvglhXCODhGRlgRiubTY0JSUxAm/5Hfs0SEi0pJALJcWs8GhXt+8JB4QvzOzP1gs6qoPyY5Bh4hIa8RMIJaT1HCllgm/cuw1RKrHoENEpEVSz7WSg04HCELba6tWqS842JfDt2bfawgQV2f2CKkegw4RkVb5a7m0u9AAtA0+gdIykPTpA2ze7Ll8bm5zUFR692lSHCcjExGR7+TcoFApmZnNR1CsWAG8+GLzY28BzL7XkKfnlGP3aVIcgw4REflOygaFgeAukIjhbtl8MIQ7cmDQISIi3ymxQaFc2rtjs7tl82oPd+SEQYeIiHyn1AaFcmjPjs2e9hpSc7ijNhh0iIjId75sUGixAPn5wPLlzf8qNcTTnqDhaa8hNYc7aoNBh4iIfCflhHOg7cTgFSuaHysxedeXoCFmp+ZA7D5NPmPQISKi9snJaQ4H3o538PdKJbGBJCcHSEsD8vKAxkbvS8OlhjsKKO6jQ0RE7edtg0KxK5W87V0jhdgdm1evlv7c/t59mnzGoENERPLwtEGhlJVKcm5yKGcgab0L8oYN/t99miRj0CEiCmX+OsIgkCuV5DgOg7sgBy0GHSKiULVyZfO8lJa7BCv1xzvQK5XacxyGXOdiUUBwMjIRUSiaMKG5h6L1UQiBnhistpVK3AU56DHoEBGFGpMJ+Phjz2Xk/uMdrCuVuAty0GPQISIKJfbN+rxR4o+32GXoasJdkIMe5+gQEYWSnTu9n9xtp9aJwf4U6LlF1G6K9ugcOXIEDzzwAHr27AmdTofCwkKP5UtLS6HT6dp81dTUKFlNIqLQISW8KD0xePv25n/VGnKA4J1bRA6KBp2GhgaMHDkSO3bskHTfhQsXUF1d7fjq0aOHQjUkIgoxYsOLTsc/3kDwzi0iB0WHru6//37cf//9ku/r0aMHunbtKqqs2WyG2Wx2PK6vrwcAWK1WWK1Wya8drOxtDaU2A2w32x0aZG334sXosHIlYLNB5+Lb9kEtW3o6bDodEMD3WjWf93PPIaypCWH5+dDZbI7Lgl7f/D4995ys75Nq2u1nSrVXJwhiB2vb+UI6Hfbt24fU1FS3ZUpLS5GUlIS+ffvCbDbjnnvuwbp16zB58mS396xbtw7r169vc33v3r2IjIyUo+pERJoyrKAAg36aStA67AgAvh80CP+zebO/q6V+FgsGHDiAO2pq0BAbi0szZ7InR0aNjY14+OGHUVdXh6ioKNmeV1VB58KFCygtLcW4ceNgNpuxe/duvPrqqygvL8eYMWNc3uOqRyc+Ph7V1dXo1q2b3M1QLavViqKiIiQnJyM8PDzQ1fEbtpvtDgVKtDssK6ttD4VOB9sTT8Dm6WwoP+LnHVrtvnbtGuLi4mQPOqpadTV48GAMHjzY8XjSpEmorKxEXl4eXn31VZf3GI1GGI3GNtfDw8ND6gfEju0OLWx3aJG13Vu2ANnZTqufdEuXQm8wwMvUW7/j5x0alGqrqoKOKxMmTMDRo0cDXQ0iIu1pz7EIREFC9RsGVlRUIC4uLtDVICIioiCkaI/OzZs3cfHiRcfjy5cvo6KiAtHR0ejTpw/WrFmDr7/+Gn/7298AAPn5+ejfvz+GDx+OW7duYffu3Th06BA++OADJatJREREGqVo0Dl58iSSkpIcj00/7UWwYMECFBQUoLq6GlVVVY7vWywWrFy5El9//TUiIyMxYsQIfPjhh07PQURERCSWokFn2rRp8LSoq6CgwOlxZmYmMuU+MZeIiIhClurn6BARERH5ikGHiIiINItBh4iIiDSLQYeIiIg0i0GHiIiINItBh4iIiDSLQYeIiIg0i0GHiIiINItBh4iIiDSLQYeIiIg0i0GHiIiINItBh4iIiDSLQYeIiIg0i0GHiIiINItBh4iIiDSLQYeIiIg0i0GHiIiINItBh4iIiDSLQYeIiIg0i0GHiIiINItBh4iIiDSLQYeIiIg0i0GHiIiINItBh4iIiDSLQYeIiIg0i0GHiIiINItBh4iIiDSLQYeIiIg0i0GHiIiINItBh4iIiDSLQYeIiIg0i0GHiIiINItBh4iIiDSLQYeIiIg0i0GHiIiINEvRoHPkyBE88MAD6NmzJ3Q6HQoLC73eU1paijFjxsBoNGLQoEEoKChQsopERESkYYoGnYaGBowcORI7duwQVf7y5cuYPXs2kpKSUFFRgYyMDDz22GM4ePCgktUkIiIijeqg5JPff//9uP/++0WX37VrF/r3748tW7YAAIYOHYqjR48iLy8PKSkpLu8xm80wm82Ox/X19QAAq9UKq9XajtoHF3tbQ6nNANvNdocGtpvtDgVKtVcnCIKgyDO3fiGdDvv27UNqaqrbMlOmTMGYMWOQn5/vuLZnzx5kZGSgrq7O5T3r1q3D+vXr21zfu3cvIiMj21ttIiIi8oPGxkY8/PDDqKurQ1RUlGzPq2iPjlQ1NTWIiYlxuhYTE4P6+nr8+OOP6NixY5t71qxZA5PJ5HhcX1+P+Ph4JCUloVu3borXWS2sViuKioqQnJyM8PDwQFfHb9hutjsUsN1sdyi4du2aIs+rqqDjC6PRCKPR2OZ6eHh4SP2A2LHdoYXtDi1sd2gJtXYr1VZVLS+PjY1FbW2t07Xa2lpERUW57M0hIiIi8kRVQScxMRHFxcVO14qKipCYmBigGhEREVEwUzTo3Lx5ExUVFaioqADQvHy8oqICVVVVAJrn18yfP99RfsmSJbh06RIyMzPx2WefYefOnfjHP/6BFStWKFlNIiIi0ihFg87JkycxevRojB49GgBgMpkwevRoPP300wCA6upqR+gBgP79++Of//wnioqKMHLkSGzZsgW7d+92u7SciIiIyBNFJyNPmzYNnlavu9r1eNq0aThz5oyCtSIiIqJQoao5OkRERERyYtAhIiIizWLQISIiIs1i0CEiIiLNYtAhIiIizWLQISIiIs1i0CEiIiLNYtAhIiIizWLQISIiIs1i0CEiIiLNYtAhIiIizWLQISIiIs1i0CEiIiLNYtAhIiIizWLQISIiIs1i0CEiIiLNYtAhIiIizeoQ6ApQcLHctmDnyZ2ovF6JgdEDsXTcUhg6GAJdLSIiIpcYdEi0zKJM5Jblokloclxb9cEqmBJNyEnOCWDNiIiIXGPQIVEyizKx6dimNtebhCbHdYYdIiJSG87RIa8sty3ILcv1WCa3LBeW2xY/1YiIiEgcBh3yaufJnU7DVa40CU3YeXKnn2pEREQkDoMOeVV5vVLWckRERP7CoENeDYweKGs5IiIif2HQIa+WjlsKvU7vsYxep8fScUv9VCMiIiJxGHTIK0MHA0yJJo9lTIkm7qdDRESqw+XlJIp96XjrfXT0Oj330SEiItVi0CHRcpJzsCFpA3dGJiKioMGgQ5IYOhiQMTEj0NUgIiIShXN0iIiISLMYdIiIiEizGHSIiIhIsxh0iIiISLMYdIiIiEizGHSIiIhIs/wSdHbs2IF+/fohIiICCQkJOHHihNuyBQUF0Ol0Tl8RERH+qCYRERFpjOJB54033oDJZMLatWtx+vRpjBw5EikpKbh69arbe6KiolBdXe34unLlitLVJCIiIg1SPOjk5uZi8eLFePTRRzFs2DDs2rULkZGReOWVV9zeo9PpEBsb6/iKiYlRuppERESkQYrujGyxWHDq1CmsWbPGcS0sLAwzZsxAWVmZ2/tu3ryJvn37wmazYcyYMXj++ecxfPhwl2XNZjPMZrPjcX19PQDAarXCarXK1BL1s7c1lNoMyNtuy20Ldp3ehcrvKzHwzoFYMmaJao+34OfNdocCtjs02y03nSAIgiLPDOCbb75Br169cOzYMSQmJjquZ2Zm4vDhwygvL29zT1lZGb744guMGDECdXV12Lx5M44cOYLz58+jd+/ebcqvW7cO69evb3N97969iIyMlLdBpFkFXxfg3W/fhQ02x7UwhGFO9zlY2Gth4CpGRBQiGhsb8fDDD6Ourg5RUVGyPa/qzrpKTEx0CkWTJk3C0KFD8dJLL+HZZ59tU37NmjUwmUyOx/X19YiPj0dSUhK6devmlzqrgdVqRVFREZKTkxEeHh7o6viNHO3OKs5C4beFba7bYEPht4UYMGAANk7f2M6ayoufN9sdCtju0Gr3tWvXFHleRYPOXXfdBb1ej9raWqfrtbW1iI2NFfUc4eHhGD16NC5evOjy+0ajEUaj0eV9ofQDYsd2S2O5bcHWE1s9ltl6YiuyZ2SrchiLn3doYbtDS6i1W6m2KjoZ2WAwYOzYsSguLnZcs9lsKC4uduq18aSpqQlnz55FXFycUtWkELbz5E40CU0eyzQJTdh5cqefakRERHJSfOjKZDJhwYIFGDduHCZMmID8/Hw0NDTg0UcfBQDMnz8fvXr1QnZ2NgDgmWeewcSJEzFo0CDcuHEDmzZtwpUrV/DYY48pXVUKQZXXK2UtR0RE6qJ40HnooYfw7bff4umnn0ZNTQ1GjRqFAwcOOJaMV1VVISzs546l77//HosXL0ZNTQ3uvPNOjB07FseOHcOwYcOUriqFoIHRA2UtR0RE6uKXychpaWlIS0tz+b3S0lKnx3l5ecjLy/NDrYiApeOWYtUHqzwOX+l1eiwdt9SPtSIiIrnwrCsKCZbbFuQfz8fy/cuRfzwfltsWAIChgwGmRJPHe02JJlVORCYiIu9Ut7ycSG6ZRZnILct16rVZ9cEqmBJNyEnOQU5yDgC0KaPX6R1liIgoODHokKZlFmVi07FNba43CU2O6/awsyFpA3ae3InK65UYGD0QS8ctZU8OEVGQY9AhzbLctiC3LNdjmdyyXGxI2gBDBwMMHQzImJjhn8oREZFfcI4OOXE3lyUYcY8cIiJij04Qsdy2KDq04m0uS7DhHjlERMSgEySUDiFi57IEE+6RQ0REHLoKAvYQ0noYxh5CMosy2/X8YueyBNsw1tJxS6HX6T2W4R45RETaxqCjcv4IIYGey6LUvCDukUNERBy6UjmxIWTX6V0YhEE+vUYg57K4GpIzHTRhVOwozB85v93zkLhHDhFRaGPQUTnRIeT7Sp+DTqDmsribFyRAwJmaMzhTc0aWeUjcI4eIKHQx6Kic6BBy50DgO99eQ8x5Tzro8Ngo+U6QFzMkB8g3Gdrfe+QovUKOiIjE4RwdlRM7oXbJmCU+v4aYuSwCBHTN6druic92YobkWgqmydCZRZmIfD4SKw6uwIsfv4gVB1cg8vlI2d47IiISj0FH5fw1oTYnOQerJ632GKrkWuUFSJ/vEywb+5kOmhRdIUdERNIw6ASBDUkbkNQvqc11vU6P1ZNWyzahNic5Bzcyb0AHncdycvSu+DLfR+0b+608uBJ5x/M8lgmmnikiIi1g0FE5+zBIyZcljms66JDUNwmNf26UfdXQ7ordECB4LCNH74qYIbnW1LyxX2ZRJnKPi5tzFAw9U0REWsGgo2LuNgoUIKDkSgmeKnnK7b2+7k3jr6XmYobkWlLzxn5iJ1bbqb1niohISxh0VKo9GwW2ZzKsP5eai5kXZKfmjf2kTqxWc88UEZHWMOiolK+7FWcVZ7VrMqy/j03ISc5B458bkZeSh1Exo9rMD5J7HpJUYnrGpPTQqLlniohIi7iPjkr5MoRksVmQfyLfY/ncslxsSNrgtnfEPqTkaiM/u5a9K3LsF2Pf4yZjYoaq9p8Re5CqlB4aNfdMERFpEYOOSvkyhHTguwOwCTaP5ZuEJjz01kPoHdXbbZAQe2yCEieq+3tjP3eknOYuZsNFAFgxcQWPnCAi8jMOXamUL0NINZYaUc9deKHQ69ydlkNKaePTkJeS57TKS+kT1QNJ6vwoUXsdTTQhN0X8hGUiIpIHg45K+bJRYKwhVvLreAom9t6V7bO2I2NihtNwldInqsvFl9VnvsyPcjex2j7HaEvKFt8aQERE7cKgo2Le/ni2HgaZeddMhOl8+0ilBBNfJ0r7m6+rz3xdYu+tF4yIiPyPc3RUTsrJ24YwAzImZCC3XPoQiT2YiJkf46+9dtpDyhyb1tqzxF4tc4yIiKgZe3SCgLshJFc2Tt8oem+a1sQGE3/uteOL9g6t+XuJPRERKYdBR4NaD6GkDk4VdZ/YYKL2INDeoTV/HaRKRETK49CVRrUcQrHctiDy+UiPf/ylBBOpe+34m9ieqbf+9Zbb4UCxS+yJiEjdGHRCgBLBRM1BQGzP1EdffYSPvvoIgOv9f6TMjyIiInVi0AliLXcR7telH/rY+rgtq0QwUWsQELuBX0v2SconvzmJ4d2HO7WFk4uJiIIXg06QcrUrcRjCkNExA1tmut6zRYlgosYgIKYHy52SL0tQ8mUJgPbv8uyOmo65ICLSOgadIORu6bQNNuSW5+JM7RkceOSA6yXoKgwmSnDXgyWFmKXoUilxbAYREbnHVVdBRszS6ZIvS9Dx+Y6ijmHwZefgYNF69dnk+Mk+PY9cuzxr+dgMIiK1YtAJMmKWTgOATbB5/ePp687BwaTlHkS/GfYbn55Djl2eg+nYDCIiLWHQCTJSdxt298czFHsXxOz/4057d3kOlmMziIi0xi9BZ8eOHejXrx8iIiKQkJCAEydOeCz/5ptvYsiQIYiIiMC9996L/fv3+6OaQUHqbsOu/niGau+CmI0A3WnvLs/BcGwGEZEWKR503njjDZhMJqxduxanT5/GyJEjkZKSgqtXr7osf+zYMcybNw+LFi3CmTNnkJqaitTUVJw7d07pqgYFX3olWv/xDNXeBcttC3p27omenXpKuk+OXZ7VfmwGEZFWKR50cnNzsXjxYjz66KMYNmwYdu3ahcjISLzyyisuy2/duhUzZ87E6tWrMXToUDz77LMYM2YMXnzxRaWrGhR86ZVo/ccz6HoXqquBdeua//2JfRJ1+sF0vHv1Xa+9Ty3nI31z8xtJLy/HLs9qPzaDiEirFF1ebrFYcOrUKaxZs8ZxLSwsDDNmzEBZWZnLe8rKymAyOf8hT0lJQWFhocvyZrMZZrPZ8bi+vh4AYLVaYbVa29kCdXpu2nNoampCXnkeBAgey+p1eiweudjpvejXpZ+o1+nXpZ863sOvvkL4+vWwzpoF3HUXsoqzkH8iHzbB5ihSsKkAGRMysHH6xja3ZxVn+XSiuw46rEhYgeemPdfu90EHHdInpHusR/qEdOgEnajXspdRxefjR2w32x0KQr3dclM06Hz33XdoampCTEyM0/WYmBh89tlnLu+pqalxWb6mpsZl+ezsbKxfv77N9ZKSEkRGRvpYc/WbgimYOGIinql8Buca3A/rPXDXA/jwgw+drvWx9UEYwmCDzc1dzZsP9rnaRxXzo7pUVmIagI+OHsXWk9ko/LawTRmb0LyH0KVLl7Cw10LHdYvNgvxP8n163dfueQ2R5kjZ3oMpmIJL3S/h3W/fdXrvwxCGOd3nYIp5iuTXKioqkqVuwYbtDi1sd2hobGxU5HmDfsPANWvWOPUA1dfXIz4+HklJSejWrVsAa+YfqUh12cOh1+mRPiHdZQ8HAGR0zPDYu5CRkIHU6alyV1e86mrgp3Cr++nfxAgjVn3yDkYLQHVnoKZz29v+33f/D39/9O+OoaZtJ7Z5DHTumBJM+M105+XoltsW7Dq9C5XfV2LgnQOxZMwSyUNaszBLluexWq0oKipCcnIywsPDJd0bzNhutjsUhGq7r127psjzKhp07rrrLuj1etTW1jpdr62tRWxsrMt7YmNjJZU3Go0wGo1troeHh4fMD8iWmVvwzNRn8MTfn4Axzoi777rb67ECW2ZugV6vV+WhnACAV14BWvXUGZel4eRP/71uKrA+qe1tTUIT/u//938duz9/WfelpJd1135XOxo/WfykT+9VeHg4Vk5eKekeT88VKj/nLbHdoYXtDg1KtVXRoGMwGDB27FgUFxcjNTUVAGCz2VBcXIy0tDSX9yQmJqK4uBgZGRmOa0VFRUhMTFSyqkHP0MGAOT3mYFbKLNE/LIE+lNPjmU9//CMwZ07zf58+DSxejL1PJGHzreZzqKpd9ObYtZxELXYV0+T4yfjNsN+4bL+7IzeUOCKCiIjkpfjQlclkwoIFCzBu3DhMmDAB+fn5aGhowKOPPgoAmD9/Pnr16oXs7GwAQHp6OqZOnYotW7Zg9uzZeP3113Hy5Em8/PLLSlc1JAXq7CuvZz7FxTV/tSCMGY0zPx246UnLcLN03FKYDpo8TtrWQYdD8w+5DHhi9xzakLSBB3MSEamQ4svLH3roIWzevBlPP/00Ro0ahYqKChw4cMAx4biqqgrVLZYNT5o0CXv37sXLL7+MkSNH4q233kJhYSHuuecepatKfuLrrsy/HfpbRZZoewpBobrnEBGRVvhlZ+S0tDRcuXIFZrMZ5eXlSEhIcHyvtLQUBQUFTuXnzp2LCxcuwGw249y5c5g1a5Y/qkl+4NOuzHFxwNq1CI/v43UPodZ73uw8udPrEnwAmPnaTJfXg27PISIichL0q660xOOcFY2Q0kPiGFKLi2veMBBATlzzXBixk6jFBpCSL0tw89ZN7K7Y7fT+c0djIqLgxqCjEl7nrGiE2ODx+bXPkX8832XoazmJ+vPvPoe52oxt/2cb7uh4R5vnkRJAOr/gPMPZdNCEFRNXQK/Tewxn3NGYiEi9GHRUIJRW9YgNHrtO7nIacmod+uyTqK1WK/bv3+92IrG1yfedNgUIyD2ei/E9x+Pjbz52W06OIyKIiEgZfpmjQ+6F2kniYg8lbT2vxttE5dbsZ1tlfiiuvCcnvznp6NlpSa/TY/Wk1ZoJoUREWsSgE2D+XtVjPwxz+f7lyD+e7/cA5cuhpC21Dn2W2xa8e/VdpB9Md7TH3aouXwkQENcpDo1/bkReSh7SxqchLyUPjX9uZMghIlI5Dl0FmD9X9ahlHpD9tVrXRQed1xVSLScqZxZlYkvZluajL346kHzlByshCN5XWUn1zoV3sHry6oDsOURERL5j0Akwf63qUds8IFe7Ml/47gJ2ndrl9d7K65Vu29PyvC8iIiIOXQWYmDkr7V3Vo9Z5QPYJxdtnbUfGxAwMvmuwqPv6dOnjtT1ye3Dwg359PSIikgeDToCJmbPS3lU9wbK7r9jQB0C2+Tdi6KBDekK6316PiIjkw6CjAjnJOVg9abViq3qCZXdfsaGvqq7KTzVqtmrSKkfQDPRkbiIikoZzdFRCyZPEg2l3X3cTlVvufJx/PN8vdWm927JaJnMTEZF4DDoqotRJ4kvHLcWqD1YFze6+3kKfmPZIpdfpkZ6Qjvgu8S5fU22TuYmISBwGnRBgHxJy9YfaTm27+3oKfWLaI9bk+Mn4zbDfeOw9EzuZe0PSBlW9h0RExKATMsQMCfmDXAeX2uu7+dhmUaeTuzM6drTXXjSfDiKVKBQOdCUiCgQGnRCi5DwgMeSe45KTnIO6W3V4+fTLPtep5bwkd2FD6cncnPtDRKQcBp0Qo9Q8IG+UmuNyd/TdPtep5bwkT2FDycncnPtDRKQsLi8nxSm5YeGSMUsQ5uOPsX1ekruzsexho6quSpFNHcW8L5uPbcbNWzclPS8REf2MQYcUp+SGhYYOBszpPkfSPS33JxITNraVb/O6YaAvk7nFvC8CBES9ECX61HYiInLGoEOKU3qOy8JeC2FKMHntdQGApH5JTqeOiw1h8V3iZd/UUWx7BQjYdGwTww4RkQ8YdEhx/tiwcOP0jbiReQM66DyWO3LliNNjKSEsJzkHjX9uRF5KHtLGpyEvJc8pNEkltb2BOI+MiCjYMeiQ4vxxcCkA7K7Y7XWpeeshMqkhrPVBpO1ZsSbmfWlJDeeREREFGwYdUpxcB5d6O2fKlyEyf4UwV8S8L60F+jwyIqJgw+Xl5Bft3bDQ3fLv9AnpmIIpAHwbIgv0rtFSNz5Uw3lkRETBhD065De+znHxtPw7tzwXBV8XAPC9d0bp0+O9yUnOQf2T9V7nF6npPDIiomDBHh3yK6kbFt68dRObj232WObdb9+F5bYFd3S8w+femUDvGt0pohNWTVoVVOeREREFAwYdUq3MokxRQzo22LDr9C6snLyyXUNkgdo12k4t55EREWkJgw6pkrujEdyp/P7nSbqB7p1pj2CuOxGRGjHokOqI2a24tYF3Ok/SDXTvTHsEc92JiNSGk5FJdcTsVtxSGMKwZMwSBWtERETBikGHVEfqXjFzus/h0A4REbnEoSuSleW2pd3zS8TuFaODDisSVmCKeYovVSUiohDAoEOycbepn9QVQ0vHLcWqD1Z5HL7SQYf6J+th1Buxf//+dtWbiIi0i0NXJAtPm/pJPXlbzNEIqyatQqeITj7VlYiIQgeDDrWbmFVSUk/eDvRuxUREpA0cugohcsyfcUXMKin7ydtSlk1zTxkiImovRXt0rl+/jkceeQRRUVHo2rUrFi1ahJs3b3q8Z9q0adDpdE5fS5Zw6XB7ZRZlIvL5SKw4uAIvfvwiVhxcgcjnIyUNKbnjy6nhYtn3lNk+azsyJmYw5BARkSSK9ug88sgjqK6uRlFREaxWKx599FE8/vjj2Lt3r8f7Fi9ejGeeecbxODIyUslqap67XYbt82cAtGsoyJdTw4mIiPxBsR6dTz/9FAcOHMDu3buRkJCAX/7yl9i+fTtef/11fPPNNx7vjYyMRGxsrOMrKipKqWpqnhLzZ1rz9dRwIiIipSnWo1NWVoauXbti3LhxjmszZsxAWFgYysvL8R//8R9u733ttdfw97//HbGxsXjggQfwl7/8xW2vjtlshtlsdjyur68HAFitVlitVplao372trZu8/YT20XNn9levh1PTHjCp9fWQYf0CenILXcfqNInpEMn6Bz1s9y2YNfpXaj8vhID7xyIJWOW+DQs5a7dWsd2s92hgO0OzXbLTbGgU1NTgx49eji/WIcOiI6ORk1Njdv7Hn74YfTt2xc9e/bEJ598gieffBIXLlzA22+/7bJ8dnY21q9f3+Z6SUlJSA55FRUVOT0+9O9Dou47VHEIg74b5PPrTsEUXOp+Ce9++y5ssDmuhyEMc7rPwRTzFMd+NwVfF7Qpl/lhJuZ0n4OFvRb69Pqt2x0q2O7QwnaHllBrd2NjoyLPKznoZGVl4YUXXvBY5tNPP/W5Qo8//rjjv++9917ExcVh+vTpqKysxMCBbed4rFmzBibTz3uu1NfXIz4+HklJSejWrZvP9Qg2VqsVRUVFSE5ORnh4uOP6xRMXsf9D7xvq/WrUrzBrwqx21WEWZnntqckqzkLht4Vt7rXBhsJvCzFgwABsnL5R9Gu6a7fWsd1sdyhgu0Or3deuXVPkeSUHnZUrV2LhwoUeywwYMACxsbG4evWq0/Xbt2/j+vXriI2NFf16CQkJAICLFy+6DDpGoxFGo7HN9fDw8JD6AbFr3e7lCcvxZPGTHoev9Do9licsR3iH9r9f4eHhWDl5pcvvWW5bsPXEVo/3bz2xFdkzsiUPY/HzDi1sd2hhu0ODUm2VHHS6d++O7t27ey2XmJiIGzdu4NSpUxg7diwA4NChQ7DZbI7wIkZFRQUAIC4uTmpVCT/vMuxq1ZWdKdHkl2XbSu23Q0RE5I5iq66GDh2KmTNnYvHixThx4gQ++ugjpKWl4Xe/+x169uwJAPj6668xZMgQnDhxAgBQWVmJZ599FqdOncKXX36Jd999F/Pnz8eUKVMwYsQIpaqqeWrZZVjJ/XaIiIhcUXQfnddeew1paWmYPn06wsLC8Otf/xrbtm1zfN9qteLChQuOCUgGgwEffvgh8vPz0dDQgPj4ePz617/GU089pWQ1Q4IadhnmfjtERORvigad6Ohoj5sD9uvXD4IgOB7Hx8fj8OHDSlYppNl3GQ4UMaeSc78dIiKSEw/1JL8Rcyq5v+YLERFRaOChnuRX9vlAuWW5Tj07ep0epkQTTyUnIiJZMeiQ36lhvhAREYUGBh0KiEDPFyIiotDAOTpERESkWQw6REREpFkMOkRERKRZDDpERESkWQw6REREpFkMOkRERKRZDDpERESkWQw6REREpFkMOkRERKRZDDpERESkWQw6REREpFkMOkRERKRZDDpERESkWQw6REREpFkMOkRERKRZDDpERESkWQw6REREpFkMOkRERKRZDDpERESkWQw6REREpFkMOkRERKRZDDpERESkWQw6REREpFkMOkRERKRZDDpERESkWQw6REREpFkMOkRERKRZDDpERESkWQw6REREpFkMOkRERKRZDDpERESkWQw6REREpFmKBZ3nnnsOkyZNQmRkJLp27SrqHkEQ8PTTTyMuLg4dO3bEjBkz8MUXXyhVRSIiItI4xYKOxWLB3Llz8ac//Un0PTk5Odi2bRt27dqF8vJy3HHHHUhJScGtW7eUqiYRERFpWAelnnj9+vUAgIKCAlHlBUFAfn4+nnrqKTz44IMAgL/97W+IiYlBYWEhfve737m8z2w2w2w2Ox7X1dUBAK5fv96O2gcfq9WKxsZGXLt2DeHh4YGujt+w3Wx3KGC72e5QYP+7LQiCrM+rWNCR6vLly6ipqcGMGTMc17p06YKEhASUlZW5DTrZ2dmOUNXS3XffrVhdiYiISBnXrl1Dly5dZHs+1QSdmpoaAEBMTIzT9ZiYGMf3XFmzZg1MJpPj8Y0bN9C3b19UVVXJ+kapXX19PeLj4/HVV18hKioq0NXxG7ab7Q4FbDfbHQrq6urQp08fREdHy/q8koJOVlYWXnjhBY9lPv30UwwZMqRdlZLCaDTCaDS2ud6lS5eQ+gGxi4qKYrtDCNsdWtju0BKq7Q4Lk3f6sKSgs3LlSixcuNBjmQEDBvhUkdjYWABAbW0t4uLiHNdra2sxatQon56TiIiIQpukoNO9e3d0795dkYr0798fsbGxKC4udgSb+vp6lJeXS1q5RURERGSn2PLyqqoqVFRUoKqqCk1NTaioqEBFRQVu3rzpKDNkyBDs27cPAKDT6ZCRkYENGzbg3XffxdmzZzF//nz07NkTqampol/XaDRi7dq1LoeztIztZrtDAdvNdocCtlvedusEuddx/WThwoX4r//6rzbXS0pKMG3atOYX1+mwZ88ex3CYIAhYu3YtXn75Zdy4cQO//OUvsXPnTq6gIiIiIp8oFnSIiIiIAo1nXREREZFmMegQERGRZjHoEBERkWYx6BAREZFmaSLoPPfcc5g0aRIiIyPRtWtXUfcIgoCnn34acXFx6NixI2bMmIEvvvhC2YrK7Pr163jkkUcQFRWFrl27YtGiRU7L912ZNm0adDqd09eSJUv8VGPf7NixA/369UNERAQSEhJw4sQJj+XffPNNDBkyBBEREbj33nuxf/9+P9VUXlLaXVBQ0OZzjYiI8GNt5XHkyBE88MAD6NmzJ3Q6HQoLC73eU1paijFjxsBoNGLQoEGiDxJWE6ntLi0tbfN563Q6j8flqE12djbGjx+Pzp07o0ePHkhNTcWFCxe83hfsv9++tFsLv99//etfMWLECMduz4mJiXj//fc93iPXZ62JoGOxWDB37lxJGwvm5ORg27Zt2LVrF8rLy3HHHXcgJSUFt27dUrCm8nrkkUdw/vx5FBUV4b333sORI0fw+OOPe71v8eLFqK6udnzl5OT4oba+eeONN2AymbB27VqcPn0aI0eOREpKCq5eveqy/LFjxzBv3jwsWrQIZ86cQWpqKlJTU3Hu3Dk/17x9pLYbaN4uvuXneuXKFT/WWB4NDQ0YOXIkduzYIar85cuXMXv2bCQlJaGiogIZGRl47LHHcPDgQYVrKi+p7ba7cOGC02feo0cPhWoov8OHD2PZsmU4fvw4ioqKYLVacd9996GhocHtPVr4/fal3UDw/3737t0bGzduxKlTp3Dy5En86le/woMPPojz58+7LC/rZy1oyJ49e4QuXbp4LWez2YTY2Fhh06ZNjms3btwQjEaj8N///d8K1lA+//rXvwQAwscff+y49v777ws6nU74+uuv3d43depUIT093Q81lMeECROEZcuWOR43NTUJPXv2FLKzs12W/+1vfyvMnj3b6VpCQoLwxz/+UdF6yk1qu8X+7AcTAMK+ffs8lsnMzBSGDx/udO2hhx4SUlJSFKyZssS0u6SkRAAgfP/9936pkz9cvXpVACAcPnzYbRmt/H63JKbdWvz9FgRBuPPOO4Xdu3e7/J6cn7UmenSkunz5MmpqajBjxgzHtS5duiAhIQFlZWUBrJl4ZWVl6Nq1K8aNG+e4NmPGDISFhaG8vNzjva+99hruuusu3HPPPVizZg0aGxuVrq5PLBYLTp065fQ5hYWFYcaMGW4/p7KyMqfyAJCSkhI0nyvgW7sB4ObNm+jbty/i4+M9/j8lLdHC590eo0aNQlxcHJKTk/HRRx8FujrtUldXBwAeT67W4uctpt2Atn6/m5qa8Prrr6OhoQGJiYkuy8j5WUs660or7OPYMTExTtdjYmKCZoy7pqamTTd1hw4dEB0d7bENDz/8MPr27YuePXvik08+wZNPPokLFy7g7bffVrrKkn333Xdoampy+Tl99tlnLu+pqakJ6s8V8K3dgwcPxiuvvIIRI0agrq4OmzdvxqRJk3D+/Hn07t3bH9UOCHefd319PX788Ud07NgxQDVTVlxcHHbt2oVx48bBbDZj9+7dmDZtGsrLyzFmzJhAV08ym82GjIwMTJ48Gffcc4/bclr4/W5JbLu18vt99uxZJCYm4tatW+jUqRP27duHYcOGuSwr52et2qCTlZWFF154wWOZTz/9FEOGDPFTjfxDbLt91XIOz7333ou4uDhMnz4dlZWVGDhwoM/PS4GVmJjo9P+MJk2ahKFDh+Kll17Cs88+G8CakRIGDx6MwYMHOx5PmjQJlZWVyMvLw6uvvhrAmvlm2bJlOHfuHI4ePRroqviV2HZr5fd78ODBqKioQF1dHd566y0sWLAAhw8fdht25KLaoLNy5UrHGVjuDBgwwKfnjo2NBQDU1tYiLi7Ocb22ttZxcnqgiG13bGxsm4mpt2/fxvXr1x3tEyMhIQEAcPHiRdUFnbvuugt6vR61tbVO12tra922MTY2VlJ5NfKl3a2Fh4dj9OjRuHjxohJVVA13n3dUVJRme3PcmTBhQlAGhbS0NMdiCm+9E1r4/baT0u7WgvX322AwYNCgQQCAsWPH4uOPP8bWrVvx0ksvtSkr52et2jk63bt3x5AhQzx+GQwGn567f//+iI2NRXFxseNafX09ysvL3Y4X+ovYdicmJuLGjRs4deqU495Dhw7BZrM5wosYFRUVAOAU+NTCYDBg7NixTp+TzWZDcXGx288pMTHRqTwAFBUVBfxzlcKXdrfW1NSEs2fPqvJzlZMWPm+5VFRUBNXnLQgC0tLSsG/fPhw6dAj9+/f3eo8WPm9f2t2aVn6/bTYbzGazy+/J+ln7MFFada5cuSKcOXNGWL9+vdCpUyfhzJkzwpkzZ4QffvjBUWbw4MHC22+/7Xi8ceNGoWvXrsI777wjfPLJJ8KDDz4o9O/fX/jxxx8D0QSfzJw5Uxg9erRQXl4uHD16VPjFL34hzJs3z/H9f//738LgwYOF8vJyQRAE4eLFi8IzzzwjnDx5Urh8+bLwzjvvCAMGDBCmTJkSqCZ49frrrwtGo1EoKCgQ/vWvfwmPP/640LVrV6GmpkYQBEH4/e9/L2RlZTnKf/TRR0KHDh2EzZs3C59++qmwdu1aITw8XDh79mygmuATqe1ev369cPDgQaGyslI4deqU8Lvf/U6IiIgQzp8/H6gm+OSHH35w/P4CEHJzc4UzZ84IV65cEQRBELKysoTf//73jvKXLl0SIiMjhdWrVwuffvqpsGPHDkGv1wsHDhwIVBN8IrXdeXl5QmFhofDFF18IZ8+eFdLT04WwsDDhww8/DFQTJPvTn/4kdOnSRSgtLRWqq6sdX42NjY4yWvz99qXdWvj9zsrKEg4fPixcvnxZ+OSTT4SsrCxBp9MJH3zwgSAIyn7Wmgg6CxYsEAC0+SopKXGUASDs2bPH8dhmswl/+ctfhJiYGMFoNArTp08XLly44P/Kt8O1a9eEefPmCZ06dRKioqKERx991CncXb582el9qKqqEqZMmSJER0cLRqNRGDRokLB69Wqhrq4uQC0QZ/v27UKfPn0Eg8EgTJgwQTh+/Ljje1OnThUWLFjgVP4f//iHcPfddwsGg0EYPny48M9//tPPNZaHlHZnZGQ4ysbExAizZs0STp8+HYBat4992XTrL3tbFyxYIEydOrXNPaNGjRIMBoMwYMAAp9/zYCG13S+88IIwcOBAISIiQoiOjhamTZsmHDp0KDCV95Gr9rb+32kt/n770m4t/H7/4Q9/EPr27SsYDAahe/fuwvTp0x0hRxCU/ax1giAI0vuBiIiIiNRPtXN0iIiIiNqLQYeIiIg0i0GHiIiINItBh4iIiDSLQYeIiIg0i0GHiIiINItBh4iIiDSLQYeIiIg0i0GHiIiINItBh4iIiDSLQYeIiIg06/8HD3ZTHB4gITkAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "visualisation(clusters, Pc, dim=dim, K=K)\n",
    "print(Pc)\n",
    "print(mean)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "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.8.10"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}