kmeans-clustering/kmeans.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clustering using Disjunctive Constraints\n",
    "\n",
    "[K-Means clustering](https://en.wikipedia.org/wiki/K-means_clustering) is one of the most used clustering problems in unsupervised learning. Typically, a heuristic algorithm is used to solve the K-Means clustering problem. Such algorithms however do not guarantee a global optimum. In this notebook, we show how K-Means can be expressed as a Generalized Disjunctive Program (GDP) with a Quadratic Rotated Cone, and we demonstrate how to implement it in MOSEK using Disjunctive Constraints (DJC). Such problem was for example studied by [Papageorgiou, Trespalacios](https://link.springer.com/article/10.1007/s13675-017-0088-0) and [Kronqvist, Misener, Tsay](https://link.springer.com/chapter/10.1007/978-3-030-78230-6_19). We further show a modification called Euclidean clustering by changing only a few lines of code. \n",
    "\n",
    "We assume a set of points $\\textbf{p}_1, \\ldots, \\textbf{p}_n \\in \\mathbb{R}^\\mathcal{D}$ and a natural number $\\mathcal{K} \\in \\{1, 2, ..., n\\}$ which specifies number of centroids. We want to find positions of the centroids such that the overall squared Euclidean distance from each point to the closest centroid is minimized. The formulation using disjunctions can look as follows\n",
    "\n",
    "$$\\begin{array}{rll}\n",
    "\\text{minimize} & \\sum_{i=1}^n d_i & \\\\\n",
    "\\text{subject to} & \\bigvee_{j \\in \\{1, .., \\mathcal{K}\\}} \\Bigl[ d_i \\geq || \\textbf{c}_j - \\textbf{p}_i ||_2^2  \\wedge Y_i = j \\Bigr], & \\forall i \\in \\{1, ..., n\\},\\\\\n",
    "& c_{j-1, 1} \\leq c_{j, 1}, & \\forall j \\in \\{2, .., \\mathcal{K}\\},\\\\\n",
    "\\\\\n",
    "& d_1, ..., d_n \\in \\mathbb{R}, & \\\\\n",
    "& \\textbf{c}_1, ..., \\textbf{c}_{\\mathcal{K}} \\in \\mathbb{R}^\\mathcal{D}, \\\\\n",
    "& Y_1, ..., Y_n \\in \\{1, .., \\mathcal{K}\\}, &\n",
    "\\end{array}$$\n",
    "\n",
    "where $\\textbf{c}_1, ..., \\textbf{c}_{\\mathcal{K}}$ are the positions of the centroids, $d_1, ..., d_n$ are auxiliary variables representing the shortest squared distance to the nearest centroid for each point and $Y_1, ..., Y_n$ are classification labels for each point, indicating the index of the nearest centroid. The first constraint is a disjunctive constraint representing the choice of a centroid for each point. Exactly one of the clauses in each of $n$ disjunctions is \"active\" and determines the index of the nearest centroid and the distance to is. The second constraint in the formulation is a symmetry-breaking constraint. \n",
    "\n",
    "**MOSEK** only supports Disjunctive Normal Form (DNF) of affine constraints. Formally, this means that each Disjunctive Constraint (DJC) is of the form $ \\bigvee_i \\bigwedge_j T_{i, j}$, where $T_{i,j}$ is an affine constraint. Such constraint is satisfied if and only if there exists at least one term $i$, such that all affine constraints $T_{i,j}$ are satisfied. We therefore need to move the non-linearity out of the disjunction. This can be tackled by a using new auxiliary variables $dAux_{i, j}$ and constraining them outside of the dicjunctions. The program then looks in the following way: \n",
    "\n",
    "$$\\begin{array}{rll}\n",
    "\\text{minimize} & \\sum_{i=1}^n d_i & \\\\\n",
    "\\text{subject to} & \\bigvee_{j \\in \\{1, .., \\mathcal{K}\\}} \\Bigl[ d_i \\geq dAux_{i, j} \\wedge \\hspace{0.2cm} Y_i = j \\Bigr], & \\forall i \\in \\{1, ..., n\\},\\\\\n",
    "& dAux_{i, j} \\geq || \\textbf{c}_j - \\textbf{p}_i ||_2^2, & \\forall j \\in \\{1, .., \\mathcal{K}\\}, \\forall i \\in \\{1, ..., n\\} \\\\\n",
    "& c_{j-1, 1} \\leq c_{j, 1}, & \\forall j \\in \\{2, .., \\mathcal{K}\\}, \\\\\n",
    "\\\\\n",
    "& dAux \\in \\mathbb{R}^{n \\times \\mathcal{K}}, &  \\\\\n",
    "& d_1, ..., d_n \\in \\mathbb{R}, & \\\\\n",
    "& \\textbf{c}_1, ..., \\textbf{c}_{\\mathcal{K}} \\in \\mathbb{R}^\\mathcal{D}, & \\\\\n",
    "& Y_1, ..., Y_n \\in \\{1, .., \\mathcal{K}\\}. &\n",
    "\\end{array}$$\n",
    "\n",
    "### Preparing synthetic data\n",
    "\n",
    "To prepare the synthetic data, we generate 3 clusters with the same number of points. These clusters are generated randomly according to the Multivariate normal distribution each with an appropriate mean vector and a covariance matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from mosek.fusion import *\n",
    "import mosek.fusion.pythonic\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "\n",
    "# make the randomness deterministic\n",
    "np.random.seed(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUCElEQVR4nO3de7BdZXnH8e9jiBIDEpWjgRCMVEURYmBOUczUdrxEFBW0aLH1gpdh1NLGKSpQZjTqqFhaMY5YS71OxToUASlREUWqVlHDRRRjHKTVJCRwKEYxBrn49I+1D+wcT86FvfZe693n+5nZk7PWXnvtZyc5v7PO86613shMJEnlelDTBUiSemOQS1LhDHJJKpxBLkmFM8glqXAGuSQVziCXehQRJ0bEt5quQ3OXQa6+iIgTIuK7EbEjIm7tfP2miIima5soIq6MiNf3ad/LIiIj4jedxy0RcWlEPGcW+/AHhaZkkKt2EXEKsBY4C1gMPBp4A7ASePCAa9ljkO83hUWZuRfwFOBy4KKIOLHZkjQ0MtOHj9oewD7ADuDPp9nuIcA/Ar8AbgE+CizoPPdnwGbgFOBWYCvwmlm+9lRgG/BvwMOBS4Ex4Jedrw/obP8e4F7gTuA3wIc7659IFbi3AxuBl3W9/yOBS4BfA98D3g18azefcxmQwB4T1r+lU/uDOsunAT8D7gB+DLy4s/5Jndru7dS3vbP+GODaTg2bgDVN/9v7aO7hEbnqdhRV0H5hmu3eDzwBWAE8DlgCvL3r+cVUPxSWAK8DzomIh8/itY8AHgOcRPWb5yc7ywcCO4EPA2TmGcA3gZMzc6/MPDkiFlKF+GeBRwEvBz4SEU/u7P8cqnDdD3ht5zFbF3b2fXBn+WfAn3Q+8zuBz0TEfpm5geq3me906lvU2X4H8CpgEVWovzEijnsAdWgYNP2TxMdwPYBXANsmrPs2sJ0qQJ8BBFUQ/VHXNkcB/9P5+s862+7R9fytwNNm+Nq7gD2nqHEF8Muu5SuB13ct/wXwzQmv+RfgHcA84G7giV3PvZfZH5Hv2Vm/cjevuw44tvP1ibvbf9f2HwTObvrf30czj7b0DzU8/g/YNyL2yMx7ADLz6QARsZnq6HgEeChwddfYZ1CF5H37GX99x2+BvWb42rHMvPO+JyMeCpwNHE3VZgHYOyLmZea9k3yGxwBPjYjtXev2oGrTjHS+3tT13M8n/ZuY2pLOn7d3anwV8HdUwQ/VZ913dy+OiKcCZwKHUo07PAT4jwdQh4aArRXV7TvA74Bjp9jmNqoj7idn5qLOY5+sBgOnM5PXTryl5ylULYynZubDqH4rgOoHwGTbbwL+q2v/i7Jqa7yRqs9+D7C0a/sDZ1D3RC+m+i1jY0Q8BvhX4GTgkVm1T340RX1QtX0uAZZm5j5U4wStOyNIg2GQq1aZuZ2qx/uRiDg+IvaKiAdFxApgYWeb31MF19kR8SiAiFgSEc+dwf4fyGv3pgr/7RHxCKoWSbdbgIO6li8FnhARr4yI+Z3HH0fEkzpH8BcCayLioRFxCPDq6eoeFxGPjoiTOzWc3vk8C6nCeqyzzWuojrS76zsgIrrP+NkbuD0z74yII4G/nGkNGj4GuWqXmf9A1SZ4G9VR5y1UPeZTqfrldL6+EbgqIn4NfJX7B/6mM9vXfhBYQHU0fxXw5QnPrwWOj4hfRsSHMvMOYBVwAnAz1dkv76dqX0B15LxXZ/2nqAZSp7M9InYAPwSeD7w0Mz8BkJk/Bv6J6reZW4DDgP/ueu0VwA3Atoi4rbPuTcC7IuIOqoHe82dQg4ZUZDqxhCSVzCNySSqcQS5JhTPIJalwBrkkFa6RC4L23XffXLZsWRNvLUnFuvrqq2/LzJGJ6xsJ8mXLlrF+/fom3lqSihURk15FbGtFkgpnkEtS4QxySSqcQS5JhTPIJalwBrmKt+6mday6YBXLP72cVResYt1N65ouSRooJ5ZQ0dbdtI41317DnfdW80hs3bGVNd9eA8AxBx3TYGXS4HhErqKtvWbtfSE+7s5772TtNWsbqkgaPINcRdu2Y9us1kvDyCBXa82k97144eJJX7u79dIwMsjVSuO97607tpLkfb3viWG++ojV7Dlvz13W7TlvT1YfsXqQ5UqNMsjVSjPtfR9z0DGsefoa9lu4H0Gw38L9WPP0NQ50ak7xrBW10mx638ccdIzBrTnNI3K1kr1vaeYMcrWSvW9p5mytqJXGWyVrr1nLth3bWLxwMauPWG0LRZqEQa7WsvctzYytFUkqnEEuSYUzyCWpcAa5JBXOIJekwhnkklQ4g1y7cLYdqTyeR677ONuOVCaPyHUfZ9uRymSQD5Fe2yLOtiOVySAfEjOdiGEq3nFQKlMtQR4RiyLigoj4SURsiIij6tivZq6Otoh3HJTKVNdg51rgy5l5fEQ8GHhoTfvVDNXRFvGOg1KZeg7yiHgY8AzgRIDMvAu4q9f9anYWL1zM1h1bJ10/G95xUCpPHa2Vg4Ax4JMRcW1EfCwiFtawX82CbRFp7qojyPcAjgD+OTMPB3YAp03cKCJOioj1EbF+bGyshrdVNychluauOoJ8M7A5M7/bWb6AKth3kZnnZuZoZo6OjIzU8Lbqtu6mdfa2pTmq5yDPzG3Apog4uLPqWcCPe92vZq6OUw8llauu88j/BjgvIq4HVgDvrWm/mgGvyJTmtlpOP8zM64DROval2fOKTGlu88rOIeAVmdLcZpAPAU89lOY2b2M7BLwiU5rbDPIh4RWZ0txla0WSCmeQS1LhDHJJKpxBLkmFM8glqXAGuSQVziCXpMIZ5JJUOINckgpnkEtS4QzyAVt30zpWXbCK5Z9ezqoLVjn5g6Seea+VARqfyWd8EojxmXwA75Mi6QHziHyAnMlHUj8Y5APkTD6S+sEgHyBn8pHUDwb5ADmTj6R+cLBzgJzJR1I/GOQD5kw+kupma0WSCmeQt4AXCUnqha2VhnmRkKReeUTeMC8SktQrg7xhXiQkqVcGecO8SEhSrwzyhnmRkKRe1TbYGRHzgPXAlsx8QV37HXZeJCSpV3WetbIa2AA8rMZ9zgleJCSpF7W0ViLiAOAY4GN17E+SNHN1HZF/EHgbsPfuNoiIk4CTAA488MCa3laSZufia7dw1mUbuXn7TvZftIC3Pvdgjjt8SdNl9aTnI/KIeAFwa2ZePdV2mXluZo5m5ujIyEivbytJs3bxtVs4/cIfsmX7ThLYsn0np1/4Qy6+dkvTpfWkjtbKSuBFEfG/wOeAZ0bEZ2rYryTV6qzLNrLz7nt3Wbfz7ns567KNDVVUj56DPDNPz8wDMnMZcAJwRWa+oufKJKlmN2/fOav1pfA8cklzxv6LFsxqfSlqDfLMvNJzyCW11VufezAL5s/bZd2C+fN463MPbqiienj3Q0lzxvjZKcN21opBLmlOOe7wJcUH90T2yCWpcAa5JBXOIJekwhnkklQ4BzslNW4Y738ySAa5pEaN3/9k/NL58fufAIb5DNlakdSoYb3/ySAZ5JIaNaz3Pxkkg1xSo4b1/ieDZJBLatSw3v9kkBzslNSoYb3/ySAZ5JIaN4z3PxkkWyuSVDiDXJIKZ5BLUuEMckkqnEEuSYUzyCWpcAa5JBXOIJekwhnkklQ4r+yU1FdzddKIQX5ug1xS38zVSSMG/bltrUjqm7k6acSgP7dBLqlv5uqkEYP+3LZWJPXN/osWsGWS8Or3pBFN9+UH/bk9IpfUN01MGjHen96yfSfJ/f3pi6/d0rf3nGjQn7vnII+IpRHx9YjYEBE3RMTqOgqTVL7jDl/C+15yGEsWLSCAJYsW8L6XHNbXo+M29OUH/bnraK3cA5ySmddExN7A1RFxeWb+uIZ9SyrcoCeNaEtffpCfu+cj8szcmpnXdL6+A9gADO95RZJabS5O5lxrjzwilgGHA9+d5LmTImJ9RKwfGxur820l6T5T9acvvnYLK8+8gseeto6VZ14x0L55P9UW5BGxF/B54M2Z+euJz2fmuZk5mpmjIyMjdb2tJO1id/1poPFB0H6p5fTDiJhPFeLnZeaFdexTkh6oyfrTK8+8YreDoKVfZVrHWSsBfBzYkJkf6L0kSapfWwZB+6GO1spK4JXAMyPius7j+TXsV5JqM8yDoHWctfKtzIzMXJ6ZKzqPL9ZRnCTVpYmLkwbFS/QlzQnjffBhvKWuQS5pzhj0xUmD4r1WJKlwBrkkFc4gl6TCGeSSVDiDXJIK51krkjQA/Zy1yCCXpD4bn7Vo/F4v4zfsAmoJc1srktRn/Z61yCCXpD7r9w27DHJJ6rN+37DLIJekPuv3Dbsc7JSkPuv3DbsMckkagH7esMvWiiQVziCXpMIZ5JJUOINckgpnkEtS4QxySSqcQS5JhTPIJalwBrkkFc4gl6TCGeSSVDiDXJIKZ5BLUuEMckkqXC1BHhFHR8TGiLgxIk6rY5+SpJnpOcgjYh5wDvA84BDg5RFxSK/7lSTNTB1H5EcCN2bmTZl5F/A54Nga9itJmoE6gnwJsKlreXNn3S4i4qSIWB8R68fGxmp4W0kS1BPkMcm6/IMVmedm5mhmjo6MjNTwtpIkqCfINwNLu5YPAG6uYb+SpBmoI8i/Dzw+Ih4bEQ8GTgAuqWG/kqQZ2KPXHWTmPRFxMnAZMA/4RGbe0HNlkqQZ6TnIATLzi8AX69iXJGl2vLJTkgpnkEtS4QxySSqcQS5JhTPIJalwBrkkFc4gl6TCGeSSVDiDXJIKZ5BLUuEMckkqnEEuSYUzyCWpcAa5JBXOIJekwhnkklQ4g1ySCmeQS1LhDHJJKlw5QX79+XD2obBmUfXn9ec3XZEktUItky/33fXnw3/+Ldy9s1r+1aZqGWD5y5qrS5JaoIwj8q+96/4QH3f3zmq9JM1xZQT5rzbPbr0kzSFlBPk+B8xuvSTNIWUE+bPeDvMX7Lpu/oJqvQbPgWepVcoY7Bwf0Pzau6p2yj4HVCHuQOfgOfAstU4ZQQ5VSBgUzZtq4Nl/H6kRZbRW1B4OPEutY5Brdhx4llqnpyCPiLMi4icRcX1EXBQRi2qqS23lwLPUOr0ekV8OHJqZy4GfAqf3XpJabfnL4IUfgn2WAlH9+cIP2R+XGtTTYGdmfqVr8Srg+N7KUREceJZapc4e+WuBL+3uyYg4KSLWR8T6sbGxGt9WPSn9nPDS65dqMO0ReUR8FVg8yVNnZOYXOtucAdwDnLe7/WTmucC5AKOjo/mAqlW9Sj8nvPT6pZpMG+SZ+eypno+IVwMvAJ6VmQZ0SUo/J7z0+qWa9NQjj4ijgVOBP83M39ZTkgam9HPCS69fqkmvPfIPA3sDl0fEdRHx0Rpq0qCUfk546fVLNekpyDPzcZm5NDNXdB5vqKswDUDp54SXXr9UE6/snMtKPye89PqlmkQT45Ojo6O5fv36gb+vJJUsIq7OzNGJ6z0il6TCGeQaHC/ekfqinPuRqznXn9/7pB5evCP1jUfkmtp4AP9qE5D3B/Bsj6anunhHUk8Mck2trgD24h2pbwxyTa2uAPbiHalvDHJNra4A9uIdqW8Mck2trgD24h2pbzxrRVMbD9pez1oZ35fBLdXOINf0DGCp1WytSFLhDHJJKpxBLs2EtxdQi9kjl6bj7QXUch6RS9Px9gJqOYNcmo63F1DLGeTSdLy9gFrOIJem4+0F1HIGuTQdby+glvOsFWkmvLpVLeYRuSQVziCXpMIZ5JJUOINckgpnkEtS4SIzB/+mEWPAz2ve7b7AbTXvsy5trg3aXV+ba4N219fm2qDd9bW1tsdk5sjElY0EeT9ExPrMHG26jsm0uTZod31trg3aXV+ba4N219fm2iZja0WSCmeQS1LhhinIz226gCm0uTZod31trg3aXV+ba4N219fm2v7A0PTIJWmuGqYjckmakwxySSrc0AV5RLwlIjIi9m26lm4R8e6IuD4irouIr0TE/k3XNC4izoqIn3TquygiFjVdU7eIeGlE3BARv4+IVpwSFhFHR8TGiLgxIk5rup5uEfGJiLg1In7UdC0TRcTSiPh6RGzo/JuubrqmbhGxZ0R8LyJ+0KnvnU3XNBNDFeQRsRR4DvCLpmuZxFmZuTwzVwCXAm2aleBy4NDMXA78FDi94Xom+hHwEuAbTRcCEBHzgHOA5wGHAC+PiEOarWoXnwKObrqI3bgHOCUznwQ8Dfjrlv3d/Q54ZmY+BVgBHB0RT2u2pOkNVZADZwNvA1o3gpuZv+5aXEiLaszMr2TmPZ3Fq4BWzWGWmRsyc2PTdXQ5ErgxM2/KzLuAzwHHNlzTfTLzG8DtTdcxmczcmpnXdL6+A9gALGm2qvtl5TedxfmdR2u+V3dnaII8Il4EbMnMHzRdy+5ExHsiYhPwV7TriLzba4EvNV1Eyy0BNnUtb6ZFYVSKiFgGHA58t+FSdhER8yLiOuBW4PLMbFV9kylqhqCI+CqweJKnzgD+Hlg12Ip2NVV9mfmFzDwDOCMiTgdOBt7Rlto625xB9avveYOqa9xM6muRmGRd64/a2iQi9gI+D7x5wm+rjcvMe4EVnbGiiyLi0Mxs3XhDt6KCPDOfPdn6iDgMeCzwg4iAqjVwTUQcmZnbmq5vEp8F1jHAIJ+utoh4NfAC4FnZwMUFs/i7a4PNwNKu5QOAmxuqpTgRMZ8qxM/LzAubrmd3MnN7RFxJNd7Q6iAfitZKZv4wMx+VmcsycxnVN9oRgwzx6UTE47sWXwT8pKlaJoqIo4FTgRdl5m+brqcA3wceHxGPjYgHAycAlzRcUxGiOtL6OLAhMz/QdD0TRcTI+FlbEbEAeDYt+l7dnaEI8kKcGRE/iojrqVpAbTrt6sPA3sDlndMjP9p0Qd0i4sURsRk4ClgXEZc1WU9nYPhk4DKqwbrzM/OGJmvqFhH/DnwHODgiNkfE65quqctK4JXAMzv/166LiOc3XVSX/YCvd75Pv0/VI7+04Zqm5SX6klQ4j8glqXAGuSQVziCXpMIZ5JJUOINckgpnkEtS4QxySSrc/wMvhZM2L2Ns9AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# specify the number of points\n",
    "nData = 21\n",
    "\n",
    "# generate data\n",
    "numberOfClusters = 3\n",
    "size = nData // numberOfClusters\n",
    "\n",
    "pointsClass1 = np.random.multivariate_normal(np.array([1, 1])*2, np.array([[1, 0], [0, 1]])*0.7, size=size)\n",
    "pointsClass2 = np.random.multivariate_normal(np.array([-1, -1])*2, np.array([[1, 0], [0, 1]])*0.7, size=size)\n",
    "pointsClass3 = np.random.multivariate_normal(np.array([-1, 3])*2, np.array([[1, 0], [0, 1]])*0.7, size=size)\n",
    "\n",
    "points = np.vstack((pointsClass1, pointsClass2, pointsClass3))\n",
    "\n",
    "# plot the generated data\n",
    "\n",
    "labels = np.zeros(3*size)\n",
    "labels[size:2*size] = 1\n",
    "labels[2*size:] = 2\n",
    "\n",
    "plt.title(\"Generated Data\")\n",
    "plt.scatter(pointsClass1[:, 0], pointsClass1[:, 1])\n",
    "plt.scatter(pointsClass2[:, 0], pointsClass2[:, 1])\n",
    "plt.scatter(pointsClass3[:, 0], pointsClass3[:, 1])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fusion Model\n",
    "In the following block, we show the K-Means model in the **Mosek Fusion for Python** in a vectorized fashion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Problem\n",
      "  Name                   :                 \n",
      "  Objective sense        : minimize        \n",
      "  Type                   : CONIC (conic optimization problem)\n",
      "  Constraints            : 2               \n",
      "  Affine conic cons.     : 63              \n",
      "  Disjunctive cons.      : 21              \n",
      "  Cones                  : 0               \n",
      "  Scalar variables       : 112             \n",
      "  Matrix variables       : 0               \n",
      "  Integer variables      : 21              \n",
      "\n",
      "Optimizer started.\n",
      "Mixed integer optimizer started.\n",
      "Threads used: 64\n",
      "Presolve started.\n",
      "Presolve terminated. Time = 0.00, probing time =  0.00\n",
      "Presolved problem: 420 variables, 290 constraints, 664 non-zeros\n",
      "Presolved problem: 21 general integer, 63 binary, 336 continuous\n",
      "Presolved problem: 63 cones\n",
      "Presolved problem: 63 disjunctions\n",
      "Clique table size: 63\n",
      "BRANCHES RELAXS   ACT_NDS  DEPTH    BEST_INT_OBJ         BEST_RELAX_OBJ       REL_GAP(%)  TIME  \n",
      "0        1        1        0        NA                   0.0000000000e+00     NA          0.0   \n",
      "0        1        1        0        3.6682139754e+02     0.0000000000e+00     100.00      0.0   \n",
      "Cut generation started.\n",
      "0        1        1        0        3.6682139754e+02     0.0000000000e+00     100.00      0.0   \n",
      "Cut generation terminated. Time = 0.00\n",
      "9        13       10       2        3.6682139754e+02     0.0000000000e+00     100.00      0.1   \n",
      "23       27       24       3        3.6682139754e+02     0.0000000000e+00     100.00      0.1   \n",
      "41       45       42       5        3.6682139754e+02     0.0000000000e+00     100.00      0.1   \n",
      "54       58       55       6        3.6682139754e+02     0.0000000000e+00     100.00      0.1   \n",
      "72       76       73       7        3.6682139754e+02     0.0000000000e+00     100.00      0.1   \n",
      "95       99       96       8        3.6682139754e+02     0.0000000000e+00     100.00      0.1   \n",
      "128      132      129      9        3.6682139754e+02     0.0000000000e+00     100.00      0.1   \n",
      "206      210      207      11       3.6682139754e+02     0.0000000000e+00     100.00      0.2   \n",
      "398      402      399      14       3.6682139754e+02     5.7818891965e-09     100.00      0.2   \n",
      "782      786      783      20       3.6682139754e+02     5.7818891965e-09     100.00      0.2   \n",
      "1550     1554     1489     11       2.9209669476e+02     5.7987916610e-09     100.00      0.3   \n",
      "3022     2964     2649     15       2.9209669476e+02     1.3458281073e-08     100.00      0.4   \n",
      "4885     4827     4054     23       2.5587728623e+02     1.0159788076e-06     100.00      0.5   \n",
      "6751     6686     5698     31       6.1363276101e+01     2.4285911784e+00     96.04       0.6   \n",
      "9534     8501     5095     17       3.0246866008e+01     4.8619330908e+00     83.93       0.7   \n",
      "14484    8860     389      29       2.9900620501e+01     1.1705980010e+01     60.85       0.8   \n",
      "14868    8886     31       23       2.9900620501e+01     1.1904381561e+01     60.19       0.8   \n",
      "14899    8901     12       24       2.9900620501e+01     1.1904381561e+01     60.19       0.8   \n",
      "14917    8917     12       26       2.9900620501e+01     1.6935494637e+01     43.36       0.9   \n",
      "14931    8930     6        25       2.9900620501e+01     1.6935494640e+01     43.36       0.9   \n",
      "An optimal solution satisfying the relative gap tolerance of 1.00e-02(%) has been located.\n",
      "The relative gap is 0.00e+00(%).\n",
      "An optimal solution satisfying the absolute gap tolerance of 0.00e+00 has been located.\n",
      "The absolute gap is 0.00e+00.\n",
      "\n",
      "Objective of best integer solution : 2.990062050092e+01      \n",
      "Best objective bound               : 2.990062050092e+01      \n",
      "Initial feasible solution objective: Undefined\n",
      "Construct solution objective       : Not employed\n",
      "User objective cut value           : Not employed\n",
      "Number of cuts generated           : 0\n",
      "Number of branches                 : 14937\n",
      "Number of relaxations solved       : 8936\n",
      "Number of interior point iterations: 125523\n",
      "Number of simplex iterations       : 0\n",
      "Time spend presolving the root     : 0.00\n",
      "Time spend optimizing the root     : 0.00\n",
      "Mixed integer optimizer terminated. Time: 0.88\n",
      "\n",
      "Optimizer terminated. Time: 0.96    \n",
      "\n",
      "\n",
      "Integer solution solution summary\n",
      "  Problem status  : PRIMAL_FEASIBLE\n",
      "  Solution status : INTEGER_OPTIMAL\n",
      "  Primal.  obj: 2.9900620501e+01    nrm: 2e+02    Viol.  con: 0e+00    var: 0e+00    acc: 0e+00    djc: 0e+00    itg: 0e+00  \n"
     ]
    }
   ],
   "source": [
    "# get shape\n",
    "n, d = points.shape\n",
    "\n",
    "xs = np.repeat( points, numberOfClusters, axis=0)\n",
    "\n",
    "# create model\n",
    "with Model() as M: \n",
    "\n",
    "    ########## create variables ##########\n",
    "\n",
    "    centroids = M.variable( [numberOfClusters, d], Domain.unbounded() )\n",
    "\n",
    "    distances = M.variable( n, Domain.greaterThan(0) )\n",
    "    distanceAux = M.variable( [n, numberOfClusters] , Domain.greaterThan(0) )\n",
    "\n",
    "    # create classification labels\n",
    "    Y = M.variable( n, Domain.integral( Domain.inRange(0, numberOfClusters-1) ) )\n",
    "\n",
    "\n",
    "    ########## create constraints ##########\n",
    "\n",
    "    # quadratic cone constraints\n",
    "    a1 = Expr.flatten( distanceAux )\n",
    "    a2 = Expr.constTerm([n*numberOfClusters, 1], 0.5)\n",
    "    a3 = Expr.repeat(centroids, n, 0) - xs\n",
    "    \n",
    "    hstack = Expr.hstack([a1, a2, a3])                      # ( d_{ij}Aux, 1/2, x_i - c_j ) in Qr\n",
    "    M.constraint( hstack, Domain.inRotatedQCone() )\n",
    "\n",
    "    # create disjunctive constraints\n",
    "    for pointInd in range(0, n):\n",
    "        label = Y.index(pointInd)\n",
    "        di = distances.index(pointInd)\n",
    "\n",
    "        ANDs = {}\n",
    "\n",
    "        # create AND constraints\n",
    "        for clusterInd in range(0, numberOfClusters):\n",
    "            dijAux = distanceAux.index([pointInd, clusterInd])\n",
    "\n",
    "            ANDs[clusterInd] = DJC.AND( DJC.term( di >= dijAux),\n",
    "                                        DJC.term( label == clusterInd) )\n",
    "        \n",
    "        M.disjunction( [ ANDs[i] for i in range(0, numberOfClusters) ] )\n",
    "\n",
    "\n",
    "    # symmetry breaking constraints\n",
    "    M.constraint( centroids[:-1,0] <= centroids[1:,0] )\n",
    "\n",
    "    ########## solve ##########\n",
    "\n",
    "    M.objective( ObjectiveSense.Minimize, Expr.sum(distances) )        \n",
    "\n",
    "    M.setLogHandler(sys.stdout)     # Enable log output\n",
    "\n",
    "    # Solve \n",
    "    M.solve()\n",
    "\n",
    "    # get centroids and clusters\n",
    "    labels_KMEANS = Y.level()\n",
    "    centroids_KMEANS = np.array(centroids.level()).reshape(numberOfClusters, d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA070lEQVR4nO3deXxU9b3w8c93JpNkskMSAtlYBILsIItWBVcWodStrbW11t6r3Wy9PpW2PtpKF6/28bFWb3uvtU9vta2t9iqlFkTUClWpgCwqSwiLbAlbSDLZt5n5PX+cM8MkZINMMjPh+369eJE563cmZ775ne/5nd8RYwxKKaVilyPSASillOodTeRKKRXjNJErpVSM00SulFIxThO5UkrFOE3kSikV4zSRh4mIjBARIyJx57j+50Xk9XDH1RdE5FkR+Umk44hFIvIlEXk30nEAiMhOEbki0nF0R0TWici/djKvV9+7gUITeTsicpmI/FNEqkWkUkTWi8jMMO/jjIPPGPO8MWZeOPdj7+sKESntxfp9mnjaf0nteKtE5JYOlj0oIi0iktVu+gf25zmir+KMNvb7Hd2bbRhjJhhj1oUpJBVBmshDiEgasBL4D2AwkAf8EGiOZFznCxGZB6wAvmyMeaGTxQ4AnwtZZxLg7vvoBo5oar2KRfNQL+kH2NZYAGPMn4wxPmNMozHmdWPMRwAi4hCRB0XkkIicFJHfiUh6RxuyW4/XhLxeJiJ/sF++bf/vEZE6EbmkfctXRD4hIu/bZwbvi8gnQuatE5Ef22cLtSLyevtWqr1cMrAayLX3UyciuSKSICI/F5Gj9r+fi0hCB+tfCDwNXGKv6wmZPUhEVtn73ygiF4SsN05E3rDPaEpE5DPdffAishj4M3CrMeYvXSz6e+CLIa9vB37XblsJIvJ/ReSwiJwQkadFxG3PGyQiK0Wk3G75rxSR/JB1O/1sRSRRRP4gIhUi4rF/LzmdvJ8CEVlu76dCRH7RwTJnnJmFnqGIyGgR+Yd9DJwSkRft6YHj50P79/LZwGdon514xDqrnByy3YMi8l0R+QioF5G40GPUPj7/bB/TtWKVXWaErD9dRLbZ8/5HRF6UTspr9rG8XkT+w459t4hc3e49Piwi64EGYFRXx7vtAhHZZM//q4gM7mTf6SLyGxE5JiJlIvITEXG2i+sJ+zP62N7vl0TkiFjf6ds72m7UM8boP/sfkAZUAM8BC4FB7eZ/GdgHjAJSgOXA7+15IwADxNmvDwLXhKy7DPhDR8va074EvGv/PBioAm4D4rBaoFVApj1/HbAf6w+P2379aCfv6QqgtN20HwEbgCFANvBP4MedrB+MK2Tas0AlMMuO73ngBXteMnAEuMOeNx04BUzoZPvrgL8CntDPq5NlDwLXACXAhYDT3tdw+/McYS/3c+AV+3NMBf4GPGLPywRuApLsef8DrGgXT4efLfAVe1tJ9r4vAtI6iNMJfAg8YX8eicBlHfyeOzoO1gH/av/8J+ABrAZXcBv2PAOMDnk9HTgJzLb3f7v9eSWEfHYfAAWAu/0xinV8NgHX2es/Amyw58UDh4B7ABdwI9AC/KSLY8YL3Gsv/1mgGhgc8h4PAxOwjpEcuj/ey4CJ9uf5Mp18l7DO6H5lLzcE2AR8pV1cd9jv8Sd2HL8EEoB5QC2QEulcdLb/tEUewhhTA1yGdWD8GigXkVdCWl2fB35mjPnYGFMH3A/cIuE/VV0E7DXG/N4Y4zXG/AnYDXwyZJnfGmP2GGMasVqyU89i+58HfmSMOWmMKccqH912ljEuN8ZsMsZ4sRJ5YP+LgYPGmN/asW/F+uLd3MW2rgT2AOt7uO9Aq/xarM+lLDBDRAS4E7jXGFNpjKkF/h24BcAYU2GMedkY02DPexiY2277nX22rVh/CEYb64xti33MtDcLyAWWGmPqjTFNxphzuc7QivVHKrcH27gT+JUxZqMd23NYJcGLQ5Z5yhhzxH5fHXnXGPOqMcaH9RlPsadfjJVgnzLGtBpjlmMlyK6cBH5uL/8i1h/fRSHznzXG7LSPn3l0f7z/3hizwxhTD3wf+EygpR1gf08XAv9mf+4nsf6Yhl5vOWAfmz7gRaw/bD8yxjQbY17H+gPVq2sPkaCJvB1jTLEx5kvGmHysFkAuVgsP++dDIYsf4nSLIpza7yewr7yQ18dDfm7AOkM41+0fsqedjc72PxyYbZ+6euxyzOeBoV1s6/tYSWeF2CUeEVktp8tBn2+3/O+BW7FaWL9rNy8bq8W8JWT/r9nTEZEkEfmVWOWxGqwyV0a7pNDZe/s9sAZ4QayS1P8REVcH76cAOGQnqd74DiDAJrvU8eUulh0OfLvd515A29/rkW721/59J9qNlFygzNjN2h5uq/3y7Y+x0PV7crwfaTfPBbQvJw63px8L+Qx+hdUyDzgR8nMjgDGm/bSz+S5FBU3kXTDG7MYqI0y0Jx3FOlgCCrFO1U5wpnqshBIQmsi6G3Ky/X4C+yrrYNnudLSvjt7H0bNYvytHgH8YYzJC/qUYY77WxTr1WKf06cBLIuIyxiy010sxxjzfJiBjDmFd9LwOq7wV6hTWl3FCyP7TjTGBL+e3gSJgtjEmDZhjT5fu3pjduvyhMWY88Amss48vdrDoEaCwB2dq9fb/HR4nxpjjxpg7jTG5WGWd/5TOe6ocAR5u97kn2a3b4Ca7iaczx4A8+2wnoKCbddov3/4YC42lJ8d7Qbt5rVi/61BHsBoEWSGfQZoxZkI3scY8TeQhxLpI922xL36JSAFWvW6DvcifgHtFZKSIpGCdsr/YScvrA6yyi8u+aBRaWigH/Fi19o68CowVkVvti1KfBcZj9ag5WyeATGl7UfZPwIMiki3WhbwfAH/ocG1r/XwRie/h/lbasd9mv3eXiMwU68Jpp+wyxwKs1tkf2582d+BfgKvsU+3Q7fixymJPiMgQABHJE5H59iKpWIneY18we6iH7wsRuVJEJtmx1WAlE18Hi27CSn6PikiyWBdJL22/kF3WKgO+ICJOu8UdetH403L6QmwVVvIL7O8EbY+fXwNfFZHZYkkWkUUiktrT99eF9+z93m0fj5/CKh91ZQjwLfv3/2msaxqvdrJsT473L4jIeBFJwrrG85JdHgkyxhwDXgceF5E0sTonXCAi7UtnA44m8rZqsS4WbRSReqwEvgOrFQfw31in129jtQibgG92sq3vY30pq7Bq0H8MzDDGNGDVZtfbp4ChdUyMMRVYrb1vY118/Q6w2BjTvgXSLfus4k/Ax/a+crEu8mwGPgK2A1vtaR15C9gJHBeRbvdvJ+R5WHXJo1in6z/FupjU3boerLr3WOB30kW3NGPMfmPM5k5mfxfrovQGu3zyJlYrHKwymRurNbcBq+zSU0OBl7CSeDHwDzr4A2gnmE9i1VoPA6VYF/w6ciewFOv3PAHrwnPATKxjsQ7r4u09xpgD9rxlwHP27/Qz9mdxJ/ALrGNuH1bpqdeMMS1YFzj/Beui9BewkmxX3XI3AmOwPueHgZvt47qj7ffkeP891tnxcawLv9/qZL9fxLo4uwvrc3gJGNbNW4x50raMpZRS3RORjcDTxpjfdjDvS1g9by7r98DOU9oiV0p1S0TmishQu/RxOzCZszubUX0oau7wUkpFtSKsrpgpWP3sb7Zr0ioKaGlFKaVinJZWlFIqxkWktJKVlWVGjBgRiV0rpVTM2rJlyyljTHb76RFJ5CNGjGDz5s56jimllOqIiLS/AxbQ0opSSsU8TeRKKRXjNJErpVSM037kSqmwam1tpbS0lKampkiHErMSExPJz8/H5epocM0zaSJXSoVVaWkpqampjBgxgrYDIKqeMMZQUVFBaWkpI0eO7NE6mshVzFuxrYzH1pRw1NNIboabpfOLuH5aXvcrqj7R1NSkSbwXRITMzEzKy8t7vI4mchXTVmwr4/7l22lstUY0LfM0cv/y7QCazCNIk3jvnO3npxc7VUx7bE1JMIkHNLb6eGxNSYQiUqr/aSJXMe2I5xg+qs+YftTT2WMpleqex+PhP//zP896vc2bN/Otb3U8VPqIESM4deqsHynQI1paUVGrq9r3iboTbDm2BVfqFurqB+H2T2+zbm6GOxIhqwEikMi//vWvnzHP5/PhdHb8AKsZM2YwY8aMvg7vDNoiV1EpUPsu8zRiOF37fva9D1i1ZxV/LfkrpxpOcfelCxjknNZmXbfLydL5RR1vWJ0Xfve73zF58mSmTJnCbbfdRnl5OTfddBMzZ85k5syZrF+/HoBly5bx5S9/mSuuuIJRo0bx1FNPAfC9732P/fv3M3XqVJYuXcq6deu48sorufXWW5k0aRJNTU3ccccdTJo0iWnTprF27VoA1q1bx+LFiwGoqKhg3rx5TJs2ja985SsERpqtr69n0aJFTJkyhYkTJ/Liiy/2+v1qi1xFpUDt2+AFnAhCY6uPn699n28vSuTi/IsZnz2eOEccI9O110o0+1vJ386YdsHgCxifPR6v38vqvavPmF+UVcTYzLE0eZt4Y/8bbeZ9suiTXe5v586dPPzww6xfv56srCwqKyu5++67uffee7nssss4fPgw8+fPp7i4GIDdu3ezdu1aamtrKSoq4mtf+xqPPvooO3bs4IMPPgCsBL1p0yZ27NjByJEjefzxxwHYvn07u3fvZt68eezZs6dNHD/84Q+57LLL+MEPfsCqVat45plnAHjttdfIzc1l1apVAFRXn1kaPFuayFVUKvPU0SKHaHJsR0wCiWYSLpOHpyaHz026jjjH6UP3+ml5mrhV0FtvvcXNN99MVlYWAIMHD+bNN99k165dwWVqamqora0FYNGiRSQkJJCQkMCQIUM4ceJEh9udNWtWsF/3u+++yze/aT2ud9y4cQwfPvyMRP7222+zfPny4D4GDRoEwKRJk7jvvvv47ne/y+LFi7n88st7/Z41kauo4jd+dpXvQlLWUNt0AHDiogCMVQXMy0hpk8RV9OuqBR3niOtyfmJcYrct8PaMMWd03/P7/bz33nu43WdeO0lIOP1ccKfTidfr7XC7ycnJbfbREx11Ixw7dixbtmzh1Vdf5f7772fevHn84Ac/6NH2OqM1chVVBOHPO//MhQX1JDlGkuZbRKpvAS4zTGvfqkeuvvpq/vznP1NRUQFAZWUl8+bN4xe/+EVwmUDJpDOpqanBFntH5syZw/PPPw/Anj17OHz4MEVFRZ0us3r1aqqqqgA4evQoSUlJfOELX+C+++5j69atZ/0e29OmjYoov/FTcqqEfxz6B1+c/EWS4pO4bfJt3DrJx67DyTz++l6tfauzMmHCBB544AHmzp2L0+lk2rRpPPXUU3zjG99g8uTJeL1e5syZw9NPP93pNjIzM7n00kuZOHEiCxcuZNGiRW3mf/3rX+erX/0qkyZNIi4ujmeffbZNyx7goYce4nOf+xzTp09n7ty5FBYWAlZdfenSpTgcDlwuF//1X//V6/cckWd2zpgxw+iDJc5vfuNnb8VeVu9bTfGpYpzi5PYptzM7f3akQ1O9VFxczIUXXhjpMGJeR5+jiGwxxpzRv1Fb5KrfNXubeXrz01YtXIQpOVOYd8E8Rg0aFenQlIpJmshVvzDGUNFYQVZSFglxCTT7mpmcM5kFoxcwatAoHZtDqV7QRK763MdVH/O3kr9RfKqYH17xQ3JScvjW7G+R4EzQBK5UGGgiV33moOcgr5S8wkcnPkJEmJozFadYtzYnxiVGODqlBg5N5KpPeBo9PPz2wxgM04ZNY9GYRYzIGBHpsJQakDSRq7A5XH2YD45/wJKiJWS4M/j0hE8zLmschemFkQ5NqQFNE7lq41yetlNaU8qK3SvYemwrgjA7bzY5KTnMu2BeP0WtVOeWLVtGSkoK9913X4fzV6xYwdixYxk/fnw/RxY+mshV0Nk+baeqsYrntz/P1mPWnWkXDbuITxV9ipyUnP4LWqleWrFiBYsXL47pRK636KugwIiDPjwYWoCOn7bj9VtjUTjFyY6TO7ho2EX86Mof8Y1Z3yA/Pb/f41aqvYcffpiioiKuueYaSkqs4/fXv/41M2fOZMqUKdx00000NDTwz3/+k1deeYWlS5cydepU9u/f3+Fy0U5b5ANIbx9CHHiqjpOMDqcfrzvOy8Uvc7TmKA9f/TBpiWk8Pu9xkuOT229KqaD+HsZ2y5YtvPDCC2zbtg2v18v06dO56KKLuPHGG7nzzjsBePDBB/nNb37DN7/5TZYsWcLixYu5+eabAcjIyOhwuWimiXyACMdDiLPTWzhW3XBGIs9Ka+CX7/+S98veB2BG7gyavE0kxiVqEldR55133uGGG24gKSkJgCVLlgCwY8cOHnzwQTweD3V1dcyfP7/D9Xu6XDQJSyIXkQzg/wETAQN82RjzXji2rXom9CHELbIfwY1pzeKxNSU9TuRXT6nkuQ0lOFouR+yqm7gOkZKzlk1lSczMnckN424gL00HrlI919/D2ELHw8d+6UtfYsWKFUyZMoVnn32WdevWdbhuT5eLJuGqkT8JvGaMGQdMAYrDtF3VQ4Hyh8FPi+MwTY6PqHO+xd6at9h8dDMNrZ3U+exB0/ZX7mdkjo8HrllEdloLLXxMXoabx65fxOenzeeRqx/h7ll3axJXUW/OnDn85S9/obGxkdraWv72N6u0U1tby7Bhw2htbQ0OLwtnDlnb2XLRrNctchFJA+YAXwIwxrSAfaVM9ZvcDDdlnkYEB8m+K/BTg1fKyUipYduxbYzLGgfA0dqj1LXUUZBWgPvhn4LHg+/x/8umsk04xYF35TImJlVxyZzLeXLBN3A4HMC/RvbNKXUWpk+fzmc/+1mmTp3K8OHDg0/g+fGPf8zs2bMZPnw4kyZNCibvW265hTvvvJOnnnqKl156qdPlolmvh7EVkanAM8AurNb4FuAeY0x9Z+voMLbh175GDtZDiB+5cRLXTc4m3hkPwLqD69hTsQeMIevFv1Hwx1Uc+txC/nJJBi1bNuEu3s+sMVdw44N/IFdb3+oc6DC24dHfw9jGAdOBbxpjNorIk8D3gO+3C+Au4C4gOMC6Cp9AHby7Xitzh8+lKLOIbce38dGSS/lN/E4OVq0mdwUsKYGbPvEv5P7016CDWSkVM8KRyEuBUmPMRvv1S1iJvA1jzDNYLXdmzJjR/0+zGOC66nroN34cYl0O2VC6gZ3lOznVcIpd5cU0jBnORX89yLK1MKoa2KhJXKlY0+tEbow5LiJHRKTIGFMCXI1VZlH9pH1ZpdTTwNLl73KwejBFeV5O1J/g1km3Eu+MJ9mVTH1LPXUtdbhdbj6zy8EXV0Fyq72xe++FJ57QZK56paMHIKueO9uSd7h6rXwTeF5EPgKmAv8epu2qHgjteuiVk9Q736LS/y7/sX41tS21jBk8Jng35rjscSTEJeByxHHl24e5/edrSf76PeD3wz33wJNPWsk8Ao8AVANDYmIiFRUVZ52MlMUYQ0VFBYmJPR/qOSz9yI0xHwBnFOBV/wh0PQRwmCTiTDZOk4m3NpPPTLiJg56DvHPoHa694FrinfF8cuwn2Vi2kQXJrSR9veh0C/yJJ6yNZGRoi1yds/z8fEpLSykvL490KDErMTGR/PyeD3ehD18eAC599C3KQpJ5QE66j4dujqO0ppTB7sHMGT6HIclDgJBTX2PaJu32r5VSUaOzXis6aNYAsHR+EW6XM/ja4APXPmZcWMLJ+pN8ouATzM6bzco9K9lXuQ8IufOtfdLWJK5UzNGxVgaAM7oepidx2eRErpswhdl5s6lsrOS1fa+RnphOXqr2DVdqoNFEPkDMHZdEarrh8uFWHdzrv4I4h1VWWbNvDemJ6Sweu1iflanUAKSJPMa1+FrYemwrO07uIM4RR2VjJUNThhLniKOupY41+9aQkZjBorGLNIkrNUBpIo9h+yr3saF0Aw2tDYzLGsesvFltknVKfAqXD7+cwvRCTeJKDWCayGPY3oq9JLuSmXfBvGBvFLAegpzgTCAnJYexmWMjGKFSqj9oIo8hLb4WthzdwsQhE0lNSOWqkVcR74xvcwfdIc8h3vj4DYalDGPR2EURjFYp1V80kccAYwx7K/eysXQjTd4mBrkHMS7BukMz1EHPQd78+E0y3Zlce8G1EYpWKdXfNJFHuVMNp1h/eD0n6k+Qk5zDwjELyUrKOmO5QBLPSsriujHXBYetVUoNfJrIo9zuU7upaa7hihFXMGbwmE4HIvq46mOyk7JZOGahJnGlzjOayKOMMYbdp3aTmZTJkOQhzMqbxay8WZ0m58AQtVeMuAKf34fL6erniJVSkaaJPIqcrD/J+sPrKW8oZ3z2eIYkD+mydb2/cj/bjm9j0ZhFuF1uHE4dcUGp85Em8ijQ5G1iU9kmdp/aTZIriatGXsXowaO7XGdf5T7WHlgbvPlHKXX+0gzQzzp6ks+o3Ar2VOxhSs4Upg+b3m15ZG/FXtYdXMew1GEsGL1AE7lS5znNAP0o9Ek+Pqo4VF3O/ctbePiGCdw8/mYyEjO63caBqgOsPbiW3NRcTeJKKUCHse1XgSf5+KiiwbmRFtlHY6uPx1/f26MkDpCTksOFWRdqEldKBWki70eBJ/k4GUSCfzxu/6w207tSVlOG3/hJciVx+fDLNYkrpYI0kfej3Ax38Od4U4jgPGN6R3af2s2qvavYfmJ7n8anlIpNmsj7Ufsn+QC4XU6Wzi/qdJ3i8mLePvQ2BWkFTBwysa9DVErFID0/70dnPMnH7rUSmN7ervJdvHv4XQrTC7l21LU4Hc4Ol1NKnd80kfez66fldZq4QwX6lmsSV0p1RxN5lEqMS2RJ0RLSE9I1iSuluqSJPAqE3iSUnnacOy4dzj1zr2Gwe3CkQ1NKxQC92BlhgZuEyjyNNMtBjtR/wON/38SKbWWRDk0pFSM0kUdY4CYhQwstjv3EmaFIyyQeW1MS6dCUUjFCSysRFrgZSIgnyXcxghvB0aObhJRSCrRFHnGhNwM5SEbsX0l3NwkppVSAJvIIO5ebhJRSKlTYSisi4gQ2A2XGmMXh2u5Ad7Y3CSmlVHvhrJHfAxQDaWHc5nmhpzcJKaVUR8JSWhGRfGAR8P/CsT2llFI9F64W+c+B7wCpnS0gIncBdwEUFhaGabdKKXV2OnpKV6yfEfe6RS4ii4GTxpgtXS1njHnGGDPDGDMjOzu7t7tVSqmzFnoDngHKPI3cv3x7zN+AF47SyqXAEhE5CLwAXCUifwjDdpVSKqxOP6WrGh9VADS2+mL+BrxeJ3JjzP3GmHxjzAjgFuAtY8wXeh2ZUkqFWeBGuxbHPpocO8+YHqu0H7lS6rwRuNHO0IqQcMb0WBXWRG6MWad9yJVS0SpwA56RZgQXMDBuwNOxVpRS541A75Rlq49wqtZB3gDptaKJXCl1XrFuwLsz0mGEldbIlVLnFb/x09DagN/4Ix1K2GgiV0qdVyobK/nDR3/gcPXhSIcSNprIlVLnlSZvE2A9F3eg0ESulDqvNLZafcbdcbHd5TCUXuxUSkVcf45/MhBb5JrIlVIRFRj/pLHVB5we/wTok2Te6G3EIQ4S4hK6XzhGaGlFKRVRoQ8g91MP9O34J4XphVycf3GfbDtStEWulIqowDgnzY79tMph4v0jiTcX9Nn4J0NThjI0ZWifbDtStEWulIqowDgn8f6RuMwwWhwf0+B8h8HpVX2yP0+Th4bWhj7ZdqRoIldKRVRg/BMHiST6J5Pkm028M57Z48rYemxr2Pf32r7X2FC6IezbjSQtrSilIqr9A8gLM3L59ry5jM2vZkTGCAAaWhtwOVy4nK5e76+xtXFA9VgBTeRKqSjQ8QPIC4I/rTu4jqrGKi4puIRRg0ad8358fh+t/tYB1YcctLSilIoBFw27CLfLzZsfv8nKPSupajy3+nmj174ZyKWJXCml+lVOSg43jLuBywov41TDKV4ufplDnkNnvZ2BeDMQaGlFKRUjRITx2eMZNWgU245tY1jqMMCqefe0hZ0Sn8KVI64kO2lgPQBeE7lSKqYkxiVyScElgDUk7co9K0mMS+TSwksZ7B7c7bpjMsf0R5j9SksrSqmYJQgTh0yksrGSl3e9zHtH3qPF19Lp8jXNNRyvO44xph+j7HuayJVSMUtEuDD7Qj478bOMyxrH9pPbeXHHi3iaPB0uv/vUblbuWYmI9G+gfUxLK0qpmJcYl8jlwy9nXNY4dpXvIj0hHYBWX+vpvufGWPX0QNdDY2CAJHRN5EqpASM7OZu5yXMBaPY28z+7/oeRGSOZ+d+vEe+ppfEb860eK8bAvfdCRgYsWxbRmMNBSytKqQFJRBiZMZJd5Tt5oWETJX94ksZnfok7LtFK4k8+CR6PldRjnLbIlVJ9qj8eGuH1e/H6vfiNnyRXEgC1zbUUpheSEJfAppu9/Emq2V+yillfWMWejyH+3jtwPXQf8Q2nSEtIC/v45P35sAxN5EqpPhP60AiDlyOequBDI64an0pNc00wCbf6WvEZHxOHTARgb8VejtUdo9XXas33t+IUJwvHLASs2/YPVB3A6/disFrVKfEp3DrpVgA2lm2ktKYUAKfDSfUVF/NBzUbK0qAiCXJvnQ37XgXgyhFXMiZzDMfrjrNyz0pcDhfxznhcTuv/WXmzGJoyFE+Th5JTJcHpgX9DU4aSGJeI1++lxdfC6o/KeXBFcb89LEMTuVKqzwQeGgHQ5NiBV44jrdfw2JoScrLcbD+5/Yx1JmRPQEQ41XCKI9VHiHPEEeeIw+V04XKcHjRraMpQEpwJwXlxjrg2d2xenH8xXr+XOEccNU3V+P/8IrIfErwwuBEyX3yFGf/7l4jDEex/nuRKYkrOFFp8LbT6W2nxtdDia8EhVhW6uqmaHSd34DO+NjEvKVrC0JShHKg6wNqDa/npa7vx+L2IM44k3ydw4A4+LKMvErlEoj/ljBkzzObNm/t9v0qp/jXye6sIZBgvpzDSSJzJxYGTD394eXBUw9BkHO+MD2sMjS0N/OWBGzFr1jBo3hI8X/wME557lffX/ZGkq+ZzzUPPMSQl56y2GRh8K5Do0xPScTldVDdVU1ZbxjU/+zuGVox4SfAXIXabWYADjy465/ciIluMMTPaT9cWuVKqz+RmuCmzn/QTRxaBrJ6b4SYtIY20hLQ+2W+gPl3mqced+iFXJw5i2dV3cfDef6Xs+AdMeux3DPt2Mn9PrqCmpZYhnF0idzqcOB3OM8ZsSU9MJz0xnZHpx6333a6dHHiIRrhprxWlVJ8JPDQilNvlZOn8oj7bZ6AuX+ZpxOCjst7Laq7l3du+T2pCGgZDXWs9Q372Kz79gxcZPXg0AIerD9PYGp7Hy/X3++51IheRAhFZKyLFIrJTRO4JR2BKqdh3/bQ8HrlxEnkZbgTIy3DzyI2T+qz3BrStywsu3P6Z+FpzeOz1PWQkZpCdlI3X7wUR4hxWUaLF18JbB97i5eKXOVp7tNcx9Pf77nWNXESGAcOMMVtFJBXYAlxvjNnV2TpaI1dK9ZXQunyo7urTlY2VvPnxm3iaPEwfNp3pw6YHL3JGi85q5L2O0hhzzBiz1f65FigG+u7PrVJKdaGzOnR39enB7sHceOGNjM0cy9ZjW1m1ZxV+4++LEMMurBc7RWQEMA3Y2MG8u4C7AAoLC8O5W6WUClo6vyjYdz0gUJ9esa2MB1b/kco6P6PSZp5xk06cI44rRlxBXmoeNc01Udci70zYErmIpAAvA/9mjKlpP98Y8wzwDFillXDtVymlQrV/mHPgrkqA+5dvp9JXj6G1y5t0QscsP153nANVB5iVNwuno+0FzGgRlkQuIi6sJP68MWZ5OLaplFLnqqOHOV/66Fs0tvqIk6H4aQbo0U06R2uPsv3kdo7VHeOaUdf0WZfJ3ghHrxUBfgMUG2N+1vuQlFIq/I7a/dldpoAEM/qM6Z2ZPmw68y6YR21zLS/vepl9lfv6NM5zEY4C0KXAbcBVIvKB/e+6MGxXKaXC5lwvggKMyBjBTeNvIjMpk7cOvMXh6sPhDq9Xel1aMca8i9WzRymlolZXF0F7IiU+hcVjF7O3Yi8FaQWAdat+NNTNY+OSrFJK9VI4btJxiIOirCJEhPqWel7Y8QLF5cV9F3QP6VgrSqnzRkcXQc+VQxwMcg/incPvUFZbxpzhc8I+4FePY4nIXpVSKsa5XW4Wjl7IrLxZHKg6wPLi5ZTXl0ckFk3kSil1jkSEqUOnsqRoCX7j73B89f6gpRWllOqlnJQcbrrwJqze2FDTXEO8M/6MYW77iiZypZQKg9Bnfr514C3qWuq4auRV5Kbm9vm+tbSilFJhdnnh5bgcLlbtWcWWo1vo6yexaYtcKaXCLDMpkxsvvJF3D7/LlmNbOFp7lMba8Tz55qE247+EqweNJnKllOoDLqeLK0deSV5aHr/d+C5/eW83Ta3WvK4G7DoXWlpRSqk+NDZzLG9/NCKYxAMCA3aFgyZypZTqY50NzNXdgF09pYlcKaX6WG8G7OoJTeRKKdXHls4vwu1qO7jW2QzY1R292KmUUn2ss6cWaa8VpZSKIeEcsKs9La0opVSM00SulFIxThO5UkrFOE3kSikV4zSRK6VUjNNErpRSMU4TuVJKxThN5EopFeM0kSulVIzTRK6UUjFOE7lSSsU4TeRKKRXjNJErpVSM00SulFIxLiyJXEQWiEiJiOwTke+FY5tKKaV6pteJXEScwC+BhcB44HMiMr6321VKKdUz4WiRzwL2GWM+Nsa0AC8AnwrDdpVSSvVAOBJ5HnAk5HWpPa0NEblLRDaLyOby8vIw7FYppRSEJ5FLB9PMGROMecYYM8MYMyM7OzsMu1VKKQXhSeSlQEHI63zgaBi2q5RSqgfCkcjfB8aIyEgRiQduAV4Jw3aVUkr1QFxvN2CM8YrI3cAawAn8tzFmZ68jU0op1SO9TuQAxphXgVfDsS2llFJnR+/sVEqpGKeJXCmlYpwmcqWUinGayJVSKsZpIldKqRiniVwppWKcJnKllIpxmsiVUirGaSJXSqkYp4lcKaVinCZypZSKcZrIlVIqxmkiV0qpGKeJXCmlYpwmcqWUinGayJVSKsZpIldKqRiniVwppWKcJnKllIpxYXlmZ39Ysa2Mx9aUcNTTSG6Gm6Xzi7h+Wl6kw1JKqYiLiUS+YlsZ9y/fTmOrDz9NlHr83L98O4Amc6XUeS8mSiuPrSmhsdUHQKNjM3XO1znlX8cDq5+n2dsMQJO3Ca/fG8kwlVIqImKiRX7U0xj8Od6Mxm9q8Esd5XVVuJwuALYc3cLO8p2kxqcyyD2IjMQMBiUOoiirKFJhK6VUv4iJRJ6b4abMTuYuMxQYCgbyMtw4xDqpGDVoFG6Xm6rGKqqaqiirKcPtcgcT+duH3qa2ubZNkh/kHkRiXGKk3pZSSoVFTCTypfOLgjXyALfLydL5p1vbw1KHMSx1WPC1MYZG7+mWfGJcIhUNFZScKqHV3wpATnIOnxr3KQC2HtuKy+FikHsQgxIHkRyf3NdvK2bphWeloktMJPJAkjib5CEiJLmSgq9n5c0Ce/G6ljo8TZ42y+8+tZu6lrrga5fDxYQhE6z1gCPVR0hPTCc1PhURCdM7iz2hF54ByjyNeuFZqQgTY0y/73TGjBlm8+bN/b7f7jS2NuJp8lDVVEVVYxXZydmMzRxLs7eZ5z58DgCnOMlIzCAjMYOirCLy0/IxxmAwwTLPQHbpo29R5mnEhwe/1BNnchGEvAw36793VaTDU2pAE5EtxpgZ7afHRIu8v7hdbtwud5sSDYDL6eL6cdcH6++eJg8n60+Sl2a1QKuaqlhevJy0hLQ29ffc1Nw2ZwUDQeDCs1eO0+o4jMuX12a6Uqr/aSLvAYc4GJI8hCHJQzqcH++MZ0rOlGCSP+Q5hMGwcPRCktKTOFZ7jA9PfNgmyWckZhDvjO/nd9J7gQvPfmlGSGgzXSkVGb1K5CLyGPBJoAXYD9xhjPGEIa6YkhKfwsy8mcHXPr+PmuYaUuJTAGj1t1LXUkdZTRk+c/qC7afHf5pB7kGcqDtBRWNFTPSkCVx4bvVmYUwacOaFZ6VU/+pti/wN4H5jjFdEfgrcD3y392HFNqfDySD3oODrwvRCCtMLMcZQ21IbLNGkJViJ8KDnIB+e+DC4fGJcIhmJGSwaswinw0lNcw0OcQT/METS6QvP8dprRako0atEbox5PeTlBuDm3oUzsIkIaQlppCWkMZzhwemz82czYcgE60KrneQbWhtwOpwAbCrbxMdVHwe7R2YkZpCVlMXEIRMj8j6un5bHtRMycLvc58UFXqWiXThr5F8GXuxspojcBdwFUFhYGMbdDgwp8SmkxKeQn5Z/xrwpOVPITc0NJvnSmlIqGyuDiXzlnpU0e5utGrzdDz4zKTPY4u/KufQJb/W18vz255mdN5spQ6ec2xsOE+3TrlQPErmIvAkM7WDWA8aYv9rLPAB4gec7244x5hngGbC6H55TtOep7ORsspOz20zz+U/X2oelDKO8oZyT9SfZX7UfgOHpw5k/ej5g3dXqjnMHW/MZiRnEOeLOrU+4MTS0NgBYPXKMgQj1q9c+7UpZuk3kxphrupovIrcDi4GrTSQ6pZ+nAmUXgItyLwr+7PV729zs5Dd+jtcdp7qpGsPpX8+0odN4bE0tDa2teOUoLXIQB8m0eDP599caWTL1U2eWTZYtA4+Hhh99B4CkODfcey9kZFjz+oAxhhZfCy2+FuIccbhdbrx+L/sq97Fs9Tt4vHUYh5cE/4UIQmOrj8fWlGgiV+eV3vZaWYB1cXOuMaYhPCGp3ohzxJGVlBV87RAHn5nwGfzGT3VTdbCL5JDkIRz1nMTQSJNjB61yGL80AMLu+gRe25fAdWOuo6m1iQOeAxSk5ZPi8cCTT1Lv8MCtl5C87GF48mm4554uW+Zev5cWXwtAsF/9/sr9NHmbaPY10+JrodnbTE5KDuOyxmGM4YUdLwTnBUzOmczF+RdjjOHtQ29TWr8d4wDBRQJjCRzO2qddnW96WyP/BZAAvGHftr7BGPPVXkelws4hDqt+HtKbJjdjD6UeQ7JvDn6pw08dXqkgJclDQVoBADvLd/LEhidwOV1kLEgnP+Ey/O88x9Htz3HwI5D/dSvN936OliPrSYlPYerQqYBVt69qrKLZ14zf+AEYkTGCeRfMA2D9kfU0eZsA649PgjMh2O1SRMhLyyPOEUe8M54EZwIJcQkMdg8GrBu0bp10K799/Z8cq249471qn3Z1vultr5XR4QpE9b/Tg5El4TBJwBDS48bwyKJJTMqxShNjMsdw5/Q7Oeg5yJGaI+yfM5GSY+9Sngxb86DgEzBq7yry0/IZNWhUcNtZSVmkJ6RbiTgugQRnAumJ6cH5N4y7AZfTRbwzvsOeL3OGz+ky9pT4FL67YEK3g6kpdT7QOzvPYz0ZjCwtIY25I+Yyl7lgDM33fpOnd0N9HLS44PDmbRya5aSxtZHC9NO9kWblzeqya2JqQmq/xK/U+UAT+Xnu+ml5PUt8xsC997LnT7/EfftVfOGhF8l84CeUP/UkG782mi2LRwVLI54mDy/teon8tHwK0wvJT8vvUVfIPo1fqQFME7nqGRHIyKD4iwsZcsedZCZnwRNPkA0sTsng2svuJ85hHU67T+3moOcgx+uOs6diD4lxiaQnpHPlyCs7Ha9GKXXuNJGrHjv67bvwlAzjiiETrAki8MQTIKHDZ1kXNecMn8MhzyEaWhvwGz9N3ibW7qrhqb/v4FD1XtJSKvnKZbO44+KZfdZaV+p8oYlcdStw9+QRz1FS0jxkzk9gbKY9s4Muh0NThjI0ZSj1LfUUnypm96nd7D7q5+X1+2hs9eEXH+UNHh79+18pqdzMnDHDKUwv5JL8S87rh3Yoda40kasuhd496WAQDTWDeHBFMU5HXLe16eT4ZGbkzmD6sOl8YtVqK4nThMOkWF0efQ28V1zLkkkZVDRUBJP45qObSXAmUJBeQEZiRj+8S6VimyZy1aXH1pTQ2OrDRxWCGweJZ333pEMcnKg+/dqJ1R/cQRKemiQWjF4QnGeM4ZDnEBWNFbxX+h6p8akUpBcwevBohqZ0NFKEUkoTuerS6bskHXjlKPFmVLvpPRN4IIWDxDOmhxIRbhp/E7XNtRypOcKR6iPsqdhDkiuJoSlDafW1UnyqmML0Qm2tK2XTRK66FEjATtJxmvQ208/G6ZuPenbzTmpCKuOzxzM+ezw+vy/4QI6T9SfZULqBDaUbSIlPoSCtgIL0AvLT8oO9ZpQ63+hg0qpLS+cX4XY520w7l7snr5+WxyM3TiIvw40AeRluHrlxUo/KM06HM/hYvLy0PG6ddCuXF15OVlIW+yr38fr+16lusmo31U3VVDZWnlVsSsU6icSAhTNmzDCbN2/u9/2qcxPNY377jZ8TdSeCD8x++9Db7D61m2RXMgXpBRSkFZCXlheTz0dVqj0R2WKMmXHGdE3kaiCpb6kP1tbLasto8bWQlpDGLRNvCc5Pjk+OcJRKnZvOErkWFdWAkhyfzLiscYzLGhdsrTf7mgGrR8zy4uWISJvaurbWVazTRK4GLIc4giUXAINhVt4sjtQc4YDnACUVJQjC7PzZTM6ZbC1jjN6UpGKOJnJ13nCIg6KsIoqyivAbPyfrT3Kk+khw/JeT9SdZs29NsLaen5ZPQpw1+EA0XydQShO5Oi85xBEcSiDAKU5yU3M55DnEnoo9CMKQ5CE01I7jR68c0GeDqqiliVwpW2ZSJlePuhpjDOUN5RyuPszR2qP8x5uHaWz14eUUcViP0dNng6pooolcqXZErJZ4oORyT/UqazquNsvps0FVtNAbgpTqRuAuVifpHU5XKtI0kSvVjXDd3apUX9HSilLd0GeDqminiVypHtBng6popqUVpZSKcZrIlVIqxmkiV0qpGKeJXCmlYpwmcqWUinERGY9cRMqBQ2HebBZwKszbDJdojg2iO75ojg2iO75ojg2iO75ojW24MSa7/cSIJPK+ICKbOxpwPRpEc2wQ3fFFc2wQ3fFFc2wQ3fFFc2wd0dKKUkrFOE3kSikV4wZSIn8m0gF0IZpjg+iOL5pjg+iOL5pjg+iOL5pjO8OAqZErpdT5aiC1yJVS6rykiVwppWLcgEvkInKfiBgRyYp0LKFE5Mci8pGIfCAir4tIbqRjChCRx0Rktx3fX0QkI9IxhRKRT4vIThHxi0hUdAkTkQUiUiIi+0Tke5GOJ5SI/LeInBSRHZGOpT0RKRCRtSJSbP9O74l0TKFEJFFENonIh3Z8P4x0TD0xoBK5iBQA1wKHIx1LBx4zxkw2xkwFVgI/iHA8od4AJhpjJgN7gPsjHE97O4AbgbcjHQiAiDiBXwILgfHA50RkfGSjauNZYEGkg+iEF/i2MeZC4GLgG1H22TUDVxljpgBTgQUicnFkQ+regErkwBPAd4Cou4JrjKkJeZlMFMVojHndGOO1X24A8iMZT3vGmGJjTEmk4wgxC9hnjPnYGNMCvAB8KsIxBRlj3gYqIx1HR4wxx4wxW+2fa4FiIGoGejeWOvuly/4XNd/VzgyYRC4iS4AyY8yHkY6lMyLysIgcAT5PdLXIQ30ZWB3pIKJcHnAk5HUpUZSMYoWIjACmARsjHEobIuIUkQ+Ak8Abxpioiq8jMfWEIBF5ExjawawHgP8NzOvfiNrqKj5jzF+NMQ8AD4jI/cDdwEPREpu9zANYp77P91dcAT2JL4pIB9OivtUWTUQkBXgZ+Ld2Z6sRZ4zxAVPta0V/EZGJxpiou94QKqYSuTHmmo6mi8gkYCTwoYiAVRrYKiKzjDHHIx1fB/4IrKIfE3l3sYnI7cBi4GoTgZsLzuKziwalQEHI63zgaIRiiTki4sJK4s8bY5ZHOp7OGGM8IrIO63pDVCfyAVFaMcZsN8YMMcaMMMaMwPqiTe/PJN4dERkT8nIJsDtSsbQnIguA7wJLjDENkY4nBrwPjBGRkSISD9wCvBLhmGKCWC2t3wDFxpifRTqe9kQkO9BrS0TcwDVE0Xe1MwMikceIR0Vkh4h8hFUCiqZuV78AUoE37O6RT0c6oFAicoOIlAKXAKtEZE0k47EvDN8NrMG6WPdnY8zOSMYUSkT+BLwHFIlIqYj8S6RjCnEpcBtwlX2sfSAi10U6qBDDgLX29/R9rBr5ygjH1C29RV8ppWKctsiVUirGaSJXSqkYp4lcKaVinCZypZSKcZrIlVIqxmkiV0qpGKeJXCmlYtz/BxmvIXZanvYLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plotTheResults(centroids, points, labels, kmeans=True):\n",
    "    # scatter centroids and data\n",
    "    plt.scatter(centroids[:, 0], centroids[:, 1], marker='x', color='r')\n",
    "    plt.scatter(points[:, 0], points[:, 1], marker='o')\n",
    "\n",
    "    # draw lines to centroids\n",
    "    for i in range(0, numberOfClusters):\n",
    "\n",
    "        centroid = centroids[i, :]\n",
    "        cluster = points[labels == i, :]\n",
    "        nPointsInCluster = cluster.shape[0]\n",
    "\n",
    "        centroidCopied = np.tile(centroid, (nPointsInCluster, 1))\n",
    "\n",
    "        xx = np.vstack( [ centroidCopied[:, 0] , cluster[:, 0] ])\n",
    "        yy = np.vstack( [ centroidCopied[:, 1] , cluster[:, 1] ])\n",
    "\n",
    "        # plot the lines\n",
    "        plt.plot(xx, yy, '--', color='green', alpha=0.4)\n",
    "\n",
    "\n",
    "    # show the points\n",
    "    if kmeans:\n",
    "        plt.title(\"Solution to the K-Means clustering problem\")\n",
    "    else: \n",
    "        plt.title(\"Solution to the Euclidean clustering problem\")\n",
    "    \n",
    "    plt.legend([\"centroids\", \"data\"])\n",
    "    plt.show()\n",
    "\n",
    "# call the function\n",
    "\n",
    "plotTheResults(centroids_KMEANS, points, labels_KMEANS, kmeans=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Euclidean Clustering\n",
    "\n",
    "In order to change the proposed model to Euclidean Clustering, we need to only change the squared norm of the distances to standard Euclidean norm. This approach is more robust to outliers compared to the K-Means algorithm. Euclidean clustering problem has the following form:\n",
    "\n",
    "$$\\begin{array}{rll}\n",
    "\\text{minimize} & \\sum_{i=1}^n d_i & \\\\\n",
    "\\text{subject to} & \\bigvee_{j \\in \\{1, .., \\mathcal{K}\\}} \\Bigl[ \\hspace{0.2cm} d_i \\geq dAux_{i, j} \\wedge Y_i = j \\Bigr], & \\forall i \\in \\{1, ..., n\\},\\\\\n",
    "& dAux_{i, j} \\geq || \\textbf{c}_j - \\textbf{p}_i ||_2, & \\forall j \\in \\{1, .., \\mathcal{K}\\}, & \\forall i \\in \\{1, ..., n\\}, \\\\\n",
    "& c_{j-1, 1} \\leq c_{j, 1} , & \\forall j \\in \\{2, .., \\mathcal{K}\\}, \\\\\n",
    "\\\\\n",
    "& dAux \\in \\mathbb{R}^{n \\times \\mathcal{K}}_+  & \\\\\n",
    "& d_1, ..., d_n \\in \\mathbb{R}_{+}, & \\\\\n",
    "& \\textbf{c}_1, ..., \\textbf{c}_{\\mathcal{K}} \\in \\mathbb{R}^\\mathcal{D}, & \\\\\n",
    "& Y_1, ..., Y_n \\in \\{1, .., \\mathcal{K}\\}. &\n",
    "\\end{array}$$\n",
    "\n",
    "Heuristic algorithms usually solve such a problem by finding a geometric median. In the **Fusion API** model, we need to only swap the Rotated Quadratic Cone for Quadratic Cone.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUCElEQVR4nO3de7BdZXnH8e9jiBIDEpWjgRCMVEURYmBOUczUdrxEFBW0aLH1gpdh1NLGKSpQZjTqqFhaMY5YS71OxToUASlREUWqVlHDRRRjHKTVJCRwKEYxBrn49I+1D+wcT86FvfZe693n+5nZk7PWXnvtZyc5v7PO86613shMJEnlelDTBUiSemOQS1LhDHJJKpxBLkmFM8glqXAGuSQVziCXehQRJ0bEt5quQ3OXQa6+iIgTIuK7EbEjIm7tfP2miIima5soIq6MiNf3ad/LIiIj4jedxy0RcWlEPGcW+/AHhaZkkKt2EXEKsBY4C1gMPBp4A7ASePCAa9ljkO83hUWZuRfwFOBy4KKIOLHZkjQ0MtOHj9oewD7ADuDPp9nuIcA/Ar8AbgE+CizoPPdnwGbgFOBWYCvwmlm+9lRgG/BvwMOBS4Ex4Jedrw/obP8e4F7gTuA3wIc7659IFbi3AxuBl3W9/yOBS4BfA98D3g18azefcxmQwB4T1r+lU/uDOsunAT8D7gB+DLy4s/5Jndru7dS3vbP+GODaTg2bgDVN/9v7aO7hEbnqdhRV0H5hmu3eDzwBWAE8DlgCvL3r+cVUPxSWAK8DzomIh8/itY8AHgOcRPWb5yc7ywcCO4EPA2TmGcA3gZMzc6/MPDkiFlKF+GeBRwEvBz4SEU/u7P8cqnDdD3ht5zFbF3b2fXBn+WfAn3Q+8zuBz0TEfpm5geq3me906lvU2X4H8CpgEVWovzEijnsAdWgYNP2TxMdwPYBXANsmrPs2sJ0qQJ8BBFUQ/VHXNkcB/9P5+s862+7R9fytwNNm+Nq7gD2nqHEF8Muu5SuB13ct/wXwzQmv+RfgHcA84G7giV3PvZfZH5Hv2Vm/cjevuw44tvP1ibvbf9f2HwTObvrf30czj7b0DzU8/g/YNyL2yMx7ADLz6QARsZnq6HgEeChwddfYZ1CF5H37GX99x2+BvWb42rHMvPO+JyMeCpwNHE3VZgHYOyLmZea9k3yGxwBPjYjtXev2oGrTjHS+3tT13M8n/ZuY2pLOn7d3anwV8HdUwQ/VZ913dy+OiKcCZwKHUo07PAT4jwdQh4aArRXV7TvA74Bjp9jmNqoj7idn5qLOY5+sBgOnM5PXTryl5ylULYynZubDqH4rgOoHwGTbbwL+q2v/i7Jqa7yRqs9+D7C0a/sDZ1D3RC+m+i1jY0Q8BvhX4GTgkVm1T340RX1QtX0uAZZm5j5U4wStOyNIg2GQq1aZuZ2qx/uRiDg+IvaKiAdFxApgYWeb31MF19kR8SiAiFgSEc+dwf4fyGv3pgr/7RHxCKoWSbdbgIO6li8FnhARr4yI+Z3HH0fEkzpH8BcCayLioRFxCPDq6eoeFxGPjoiTOzWc3vk8C6nCeqyzzWuojrS76zsgIrrP+NkbuD0z74yII4G/nGkNGj4GuWqXmf9A1SZ4G9VR5y1UPeZTqfrldL6+EbgqIn4NfJX7B/6mM9vXfhBYQHU0fxXw5QnPrwWOj4hfRsSHMvMOYBVwAnAz1dkv76dqX0B15LxXZ/2nqAZSp7M9InYAPwSeD7w0Mz8BkJk/Bv6J6reZW4DDgP/ueu0VwA3Atoi4rbPuTcC7IuIOqoHe82dQg4ZUZDqxhCSVzCNySSqcQS5JhTPIJalwBrkkFa6RC4L23XffXLZsWRNvLUnFuvrqq2/LzJGJ6xsJ8mXLlrF+/fom3lqSihURk15FbGtFkgpnkEtS4QxySSqcQS5JhTPIJalwBrmKt+6mday6YBXLP72cVResYt1N65ouSRooJ5ZQ0dbdtI41317DnfdW80hs3bGVNd9eA8AxBx3TYGXS4HhErqKtvWbtfSE+7s5772TtNWsbqkgaPINcRdu2Y9us1kvDyCBXa82k97144eJJX7u79dIwMsjVSuO97607tpLkfb3viWG++ojV7Dlvz13W7TlvT1YfsXqQ5UqNMsjVSjPtfR9z0DGsefoa9lu4H0Gw38L9WPP0NQ50ak7xrBW10mx638ccdIzBrTnNI3K1kr1vaeYMcrWSvW9p5mytqJXGWyVrr1nLth3bWLxwMauPWG0LRZqEQa7WsvctzYytFUkqnEEuSYUzyCWpcAa5JBXOIJekwhnkklQ4g1y7cLYdqTyeR677ONuOVCaPyHUfZ9uRymSQD5Fe2yLOtiOVySAfEjOdiGEq3nFQKlMtQR4RiyLigoj4SURsiIij6tivZq6Otoh3HJTKVNdg51rgy5l5fEQ8GHhoTfvVDNXRFvGOg1KZeg7yiHgY8AzgRIDMvAu4q9f9anYWL1zM1h1bJ10/G95xUCpPHa2Vg4Ax4JMRcW1EfCwiFtawX82CbRFp7qojyPcAjgD+OTMPB3YAp03cKCJOioj1EbF+bGyshrdVNychluauOoJ8M7A5M7/bWb6AKth3kZnnZuZoZo6OjIzU8Lbqtu6mdfa2pTmq5yDPzG3Apog4uLPqWcCPe92vZq6OUw8llauu88j/BjgvIq4HVgDvrWm/mgGvyJTmtlpOP8zM64DROval2fOKTGlu88rOIeAVmdLcZpAPAU89lOY2b2M7BLwiU5rbDPIh4RWZ0txla0WSCmeQS1LhDHJJKpxBLkmFM8glqXAGuSQVziCXpMIZ5JJUOINckgpnkEtS4QzyAVt30zpWXbCK5Z9ezqoLVjn5g6Seea+VARqfyWd8EojxmXwA75Mi6QHziHyAnMlHUj8Y5APkTD6S+sEgHyBn8pHUDwb5ADmTj6R+cLBzgJzJR1I/GOQD5kw+kupma0WSCmeQt4AXCUnqha2VhnmRkKReeUTeMC8SktQrg7xhXiQkqVcGecO8SEhSrwzyhnmRkKRe1TbYGRHzgPXAlsx8QV37HXZeJCSpV3WetbIa2AA8rMZ9zgleJCSpF7W0ViLiAOAY4GN17E+SNHN1HZF/EHgbsPfuNoiIk4CTAA488MCa3laSZufia7dw1mUbuXn7TvZftIC3Pvdgjjt8SdNl9aTnI/KIeAFwa2ZePdV2mXluZo5m5ujIyEivbytJs3bxtVs4/cIfsmX7ThLYsn0np1/4Qy6+dkvTpfWkjtbKSuBFEfG/wOeAZ0bEZ2rYryTV6qzLNrLz7nt3Wbfz7ns567KNDVVUj56DPDNPz8wDMnMZcAJwRWa+oufKJKlmN2/fOav1pfA8cklzxv6LFsxqfSlqDfLMvNJzyCW11VufezAL5s/bZd2C+fN463MPbqiienj3Q0lzxvjZKcN21opBLmlOOe7wJcUH90T2yCWpcAa5JBXOIJekwhnkklQ4BzslNW4Y738ySAa5pEaN3/9k/NL58fufAIb5DNlakdSoYb3/ySAZ5JIaNaz3Pxkkg1xSo4b1/ieDZJBLatSw3v9kkBzslNSoYb3/ySAZ5JIaN4z3PxkkWyuSVDiDXJIKZ5BLUuEMckkqnEEuSYUzyCWpcAa5JBXOIJekwhnkklQ4r+yU1FdzddKIQX5ug1xS38zVSSMG/bltrUjqm7k6acSgP7dBLqlv5uqkEYP+3LZWJPXN/osWsGWS8Or3pBFN9+UH/bk9IpfUN01MGjHen96yfSfJ/f3pi6/d0rf3nGjQn7vnII+IpRHx9YjYEBE3RMTqOgqTVL7jDl/C+15yGEsWLSCAJYsW8L6XHNbXo+M29OUH/bnraK3cA5ySmddExN7A1RFxeWb+uIZ9SyrcoCeNaEtffpCfu+cj8szcmpnXdL6+A9gADO95RZJabS5O5lxrjzwilgGHA9+d5LmTImJ9RKwfGxur820l6T5T9acvvnYLK8+8gseeto6VZ14x0L55P9UW5BGxF/B54M2Z+euJz2fmuZk5mpmjIyMjdb2tJO1id/1poPFB0H6p5fTDiJhPFeLnZeaFdexTkh6oyfrTK8+8YreDoKVfZVrHWSsBfBzYkJkf6L0kSapfWwZB+6GO1spK4JXAMyPius7j+TXsV5JqM8yDoHWctfKtzIzMXJ6ZKzqPL9ZRnCTVpYmLkwbFS/QlzQnjffBhvKWuQS5pzhj0xUmD4r1WJKlwBrkkFc4gl6TCGeSSVDiDXJIK51krkjQA/Zy1yCCXpD4bn7Vo/F4v4zfsAmoJc1srktRn/Z61yCCXpD7r9w27DHJJ6rN+37DLIJekPuv3Dbsc7JSkPuv3DbsMckkagH7esMvWiiQVziCXpMIZ5JJUOINckgpnkEtS4QxySSqcQS5JhTPIJalwBrkkFc4gl6TCGeSSVDiDXJIKZ5BLUuEMckkqXC1BHhFHR8TGiLgxIk6rY5+SpJnpOcgjYh5wDvA84BDg5RFxSK/7lSTNTB1H5EcCN2bmTZl5F/A54Nga9itJmoE6gnwJsKlreXNn3S4i4qSIWB8R68fGxmp4W0kS1BPkMcm6/IMVmedm5mhmjo6MjNTwtpIkqCfINwNLu5YPAG6uYb+SpBmoI8i/Dzw+Ih4bEQ8GTgAuqWG/kqQZ2KPXHWTmPRFxMnAZMA/4RGbe0HNlkqQZ6TnIATLzi8AX69iXJGl2vLJTkgpnkEtS4QxySSqcQS5JhTPIJalwBrkkFc4gl6TCGeSSVDiDXJIKZ5BLUuEMckkqnEEuSYUzyCWpcAa5JBXOIJekwhnkklQ4g1ySCmeQS1LhDHJJKlw5QX79+XD2obBmUfXn9ec3XZEktUItky/33fXnw3/+Ldy9s1r+1aZqGWD5y5qrS5JaoIwj8q+96/4QH3f3zmq9JM1xZQT5rzbPbr0kzSFlBPk+B8xuvSTNIWUE+bPeDvMX7Lpu/oJqvQbPgWepVcoY7Bwf0Pzau6p2yj4HVCHuQOfgOfAstU4ZQQ5VSBgUzZtq4Nl/H6kRZbRW1B4OPEutY5Brdhx4llqnpyCPiLMi4icRcX1EXBQRi2qqS23lwLPUOr0ekV8OHJqZy4GfAqf3XpJabfnL4IUfgn2WAlH9+cIP2R+XGtTTYGdmfqVr8Srg+N7KUREceJZapc4e+WuBL+3uyYg4KSLWR8T6sbGxGt9WPSn9nPDS65dqMO0ReUR8FVg8yVNnZOYXOtucAdwDnLe7/WTmucC5AKOjo/mAqlW9Sj8nvPT6pZpMG+SZ+eypno+IVwMvAJ6VmQZ0SUo/J7z0+qWa9NQjj4ijgVOBP83M39ZTkgam9HPCS69fqkmvPfIPA3sDl0fEdRHx0Rpq0qCUfk546fVLNekpyDPzcZm5NDNXdB5vqKswDUDp54SXXr9UE6/snMtKPye89PqlmkQT45Ojo6O5fv36gb+vJJUsIq7OzNGJ6z0il6TCGeQaHC/ekfqinPuRqznXn9/7pB5evCP1jUfkmtp4AP9qE5D3B/Bsj6anunhHUk8Mck2trgD24h2pbwxyTa2uAPbiHalvDHJNra4A9uIdqW8Mck2trgD24h2pbzxrRVMbD9pez1oZ35fBLdXOINf0DGCp1WytSFLhDHJJKpxBLs2EtxdQi9kjl6bj7QXUch6RS9Px9gJqOYNcmo63F1DLGeTSdLy9gFrOIJem4+0F1HIGuTQdby+glvOsFWkmvLpVLeYRuSQVziCXpMIZ5JJUOINckgpnkEtS4SIzB/+mEWPAz2ve7b7AbTXvsy5trg3aXV+ba4N219fm2qDd9bW1tsdk5sjElY0EeT9ExPrMHG26jsm0uTZod31trg3aXV+ba4N219fm2iZja0WSCmeQS1LhhinIz226gCm0uTZod31trg3aXV+ba4N219fm2v7A0PTIJWmuGqYjckmakwxySSrc0AV5RLwlIjIi9m26lm4R8e6IuD4irouIr0TE/k3XNC4izoqIn3TquygiFjVdU7eIeGlE3BARv4+IVpwSFhFHR8TGiLgxIk5rup5uEfGJiLg1In7UdC0TRcTSiPh6RGzo/JuubrqmbhGxZ0R8LyJ+0KnvnU3XNBNDFeQRsRR4DvCLpmuZxFmZuTwzVwCXAm2aleBy4NDMXA78FDi94Xom+hHwEuAbTRcCEBHzgHOA5wGHAC+PiEOarWoXnwKObrqI3bgHOCUznwQ8Dfjrlv3d/Q54ZmY+BVgBHB0RT2u2pOkNVZADZwNvA1o3gpuZv+5aXEiLaszMr2TmPZ3Fq4BWzWGWmRsyc2PTdXQ5ErgxM2/KzLuAzwHHNlzTfTLzG8DtTdcxmczcmpnXdL6+A9gALGm2qvtl5TedxfmdR2u+V3dnaII8Il4EbMnMHzRdy+5ExHsiYhPwV7TriLzba4EvNV1Eyy0BNnUtb6ZFYVSKiFgGHA58t+FSdhER8yLiOuBW4PLMbFV9kylqhqCI+CqweJKnzgD+Hlg12Ip2NVV9mfmFzDwDOCMiTgdOBt7Rlto625xB9avveYOqa9xM6muRmGRd64/a2iQi9gI+D7x5wm+rjcvMe4EVnbGiiyLi0Mxs3XhDt6KCPDOfPdn6iDgMeCzwg4iAqjVwTUQcmZnbmq5vEp8F1jHAIJ+utoh4NfAC4FnZwMUFs/i7a4PNwNKu5QOAmxuqpTgRMZ8qxM/LzAubrmd3MnN7RFxJNd7Q6iAfitZKZv4wMx+VmcsycxnVN9oRgwzx6UTE47sWXwT8pKlaJoqIo4FTgRdl5m+brqcA3wceHxGPjYgHAycAlzRcUxGiOtL6OLAhMz/QdD0TRcTI+FlbEbEAeDYt+l7dnaEI8kKcGRE/iojrqVpAbTrt6sPA3sDlndMjP9p0Qd0i4sURsRk4ClgXEZc1WU9nYPhk4DKqwbrzM/OGJmvqFhH/DnwHODgiNkfE65quqctK4JXAMzv/166LiOc3XVSX/YCvd75Pv0/VI7+04Zqm5SX6klQ4j8glqXAGuSQVziCXpMIZ5JJUOINckgpnkEtS4QxySSrc/wMvhZM2L2Ns9AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# use the same data\n",
    "\n",
    "plt.title(\"Generated Data\")\n",
    "plt.scatter(pointsClass1[:, 0], pointsClass1[:, 1])\n",
    "plt.scatter(pointsClass2[:, 0], pointsClass2[:, 1])\n",
    "plt.scatter(pointsClass3[:, 0], pointsClass3[:, 1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fusion Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Problem\n",
      "  Name                   :                 \n",
      "  Objective sense        : minimize        \n",
      "  Type                   : CONIC (conic optimization problem)\n",
      "  Constraints            : 2               \n",
      "  Affine conic cons.     : 63              \n",
      "  Disjunctive cons.      : 21              \n",
      "  Cones                  : 0               \n",
      "  Scalar variables       : 112             \n",
      "  Matrix variables       : 0               \n",
      "  Integer variables      : 21              \n",
      "\n",
      "Optimizer started.\n",
      "Mixed integer optimizer started.\n",
      "Threads used: 64\n",
      "Presolve started.\n",
      "Presolve terminated. Time = 0.00, probing time =  0.00\n",
      "Presolved problem: 357 variables, 290 constraints, 664 non-zeros\n",
      "Presolved problem: 21 general integer, 63 binary, 273 continuous\n",
      "Presolved problem: 63 cones\n",
      "Presolved problem: 63 disjunctions\n",
      "Clique table size: 63\n",
      "BRANCHES RELAXS   ACT_NDS  DEPTH    BEST_INT_OBJ         BEST_RELAX_OBJ       REL_GAP(%)  TIME  \n",
      "0        1        1        0        NA                   0.0000000000e+00     NA          0.0   \n",
      "0        1        1        0        8.5213968665e+01     0.0000000000e+00     100.00      0.0   \n",
      "Cut generation started.\n",
      "0        1        1        0        8.5213968665e+01     0.0000000000e+00     100.00      0.0   \n",
      "Cut generation terminated. Time = 0.00\n",
      "9        13       10       2        8.5213968665e+01     0.0000000000e+00     100.00      0.1   \n",
      "24       28       25       3        8.5213968665e+01     0.0000000000e+00     100.00      0.1   \n",
      "43       47       44       5        8.5213968665e+01     0.0000000000e+00     100.00      0.1   \n",
      "55       59       56       6        8.5213968665e+01     0.0000000000e+00     100.00      0.1   \n",
      "75       79       76       7        8.5213968665e+01     0.0000000000e+00     100.00      0.1   \n",
      "97       101      98       8        8.5213968665e+01     0.0000000000e+00     100.00      0.1   \n",
      "131      135      132      9        8.5213968665e+01     0.0000000000e+00     100.00      0.1   \n",
      "207      211      208      11       8.5213968665e+01     0.0000000000e+00     100.00      0.1   \n",
      "399      403      400      14       8.5213968665e+01     0.0000000000e+00     100.00      0.2   \n",
      "783      787      784      20       8.5213968665e+01     0.0000000000e+00     100.00      0.2   \n",
      "1551     1555     1506     7        8.0007253389e+01     0.0000000000e+00     100.00      0.2   \n",
      "3023     3019     2752     20       7.7542742673e+01     1.2384027469e-09     100.00      0.3   \n",
      "4943     4939     4282     10       6.4230909292e+01     1.2384027469e-09     100.00      0.4   \n",
      "6863     6859     6190     26       6.4230909292e+01     1.5939326100e-09     100.00      0.5   \n",
      "8778     8774     7569     21       6.4230909292e+01     3.3543150697e+00     94.78       0.5   \n",
      "10672    10668    9177     24       5.3741376919e+01     6.1996500725e+00     88.46       0.6   \n",
      "12586    12571    10627    29       3.0022302617e+01     7.8071468376e+00     74.00       0.7   \n",
      "14549    14469    11228    24       2.4972155345e+01     8.9100468605e+00     64.32       0.8   \n",
      "17060    16350    11243    32       2.4972155345e+01     9.4691749685e+00     62.08       0.9   \n",
      "19230    18233    11387    14       2.0852857330e+01     9.9792349266e+00     52.14       1.0   \n",
      "22397    20093    10512    13       2.0852857330e+01     1.0770938019e+01     48.35       1.1   \n",
      "24852    21978    10029    16       2.0852857330e+01     1.1656163183e+01     44.10       1.2   \n",
      "27019    23877    9694     20       2.0852857330e+01     1.2027249317e+01     42.32       1.3   \n",
      "29093    25778    9428     24       2.0852857330e+01     1.2350676113e+01     40.77       1.4   \n",
      "31121    27685    9144     18       2.0852857330e+01     1.2509919303e+01     40.01       1.4   \n",
      "33180    29593    8753     20       2.0852857330e+01     1.2792794764e+01     38.65       1.5   \n",
      "35230    31500    8335     16       2.0852857330e+01     1.3012038572e+01     37.60       1.6   \n",
      "37275    33413    7852     14       2.0852857330e+01     1.3338809996e+01     36.03       1.7   \n",
      "39308    35325    7359     11       2.0852857330e+01     1.3661633587e+01     34.49       1.8   \n",
      "41351    37235    6848     18       2.0852857330e+01     1.3874268700e+01     33.47       1.9   \n",
      "43417    39148    6230     17       2.0852857330e+01     1.4262195538e+01     31.61       2.0   \n",
      "45434    41013    5585     14       2.0852857330e+01     1.4539497525e+01     30.28       2.1   \n",
      "47571    42925    4798     13       2.0852857330e+01     1.4558626153e+01     30.18       2.1   \n",
      "49732    44834    3933     18       2.0852857330e+01     1.5158292439e+01     27.31       2.2   \n",
      "51896    46751    2989     20       2.0852857330e+01     1.5629371465e+01     25.05       2.3   \n",
      "54255    48662    1682     18       2.0852857330e+01     1.6574615587e+01     20.52       2.4   \n",
      "55919    49807    514      14       2.0852857330e+01     1.7679256490e+01     15.22       2.4   \n",
      "56423    49986    62       21       2.0852857330e+01     1.9995164234e+01     4.11        2.5   \n",
      "56487    49997    0        19       2.0852857330e+01     2.0852857330e+01     0.00e+00    2.5   \n",
      "An optimal solution satisfying the relative gap tolerance of 1.00e-02(%) has been located.\n",
      "The relative gap is 0.00e+00(%).\n",
      "An optimal solution satisfying the absolute gap tolerance of 0.00e+00 has been located.\n",
      "The absolute gap is 0.00e+00.\n",
      "\n",
      "Objective of best integer solution : 2.085285733044e+01      \n",
      "Best objective bound               : 2.085285733044e+01      \n",
      "Initial feasible solution objective: Undefined\n",
      "Construct solution objective       : Not employed\n",
      "User objective cut value           : Not employed\n",
      "Number of cuts generated           : 0\n",
      "Number of branches                 : 56487\n",
      "Number of relaxations solved       : 49997\n",
      "Number of interior point iterations: 655050\n",
      "Number of simplex iterations       : 0\n",
      "Time spend presolving the root     : 0.00\n",
      "Time spend optimizing the root     : 0.00\n",
      "Mixed integer optimizer terminated. Time: 2.47\n",
      "\n",
      "Optimizer terminated. Time: 2.51    \n",
      "\n",
      "\n",
      "Integer solution solution summary\n",
      "  Problem status  : PRIMAL_FEASIBLE\n",
      "  Solution status : INTEGER_OPTIMAL\n",
      "  Primal.  obj: 2.0852857330e+01    nrm: 4e+01    Viol.  con: 0e+00    var: 0e+00    acc: 0e+00    djc: 0e+00    itg: 0e+00  \n"
     ]
    }
   ],
   "source": [
    "# get shape\n",
    "n, d = points.shape\n",
    "\n",
    "xs = np.repeat( points, numberOfClusters, axis=0)\n",
    "\n",
    "# create model\n",
    "with Model() as M: \n",
    "\n",
    "    ########## create variables ##########\n",
    "\n",
    "    centroids = M.variable( [numberOfClusters, d], Domain.unbounded() )\n",
    "\n",
    "    distances = M.variable( n, Domain.greaterThan(0) )\n",
    "    distanceAux = M.variable( [n, numberOfClusters] , Domain.greaterThan(0) )\n",
    "\n",
    "    # create classification labels\n",
    "    Y = M.variable( n, Domain.integral( Domain.inRange(0, numberOfClusters-1) ) )\n",
    "\n",
    "\n",
    "    ########## create constraints ##########\n",
    "\n",
    "    # quadratic cone constraints - THIS IS THE ONLY DIFFERENCE TO K-MEANS MODEL\n",
    "    a1 = Expr.flatten( distanceAux )\n",
    "    a2 = Expr.repeat(centroids, n, 0) - xs\n",
    "    \n",
    "    hstack = Expr.hstack([a1, a2])                      # ( d_{ij}Aux, x_i - c_j )\n",
    "    M.constraint( hstack, Domain.inQCone() )            # now a quadratic instead of rotated quadratic cone\n",
    "\n",
    "    # create disjunctive constraints\n",
    "    for pointInd in range(0, n):\n",
    "        label = Y.index(pointInd)\n",
    "        di = distances.index(pointInd)\n",
    "\n",
    "        ANDs = {}\n",
    "\n",
    "        # create AND constraints\n",
    "        for clusterInd in range(0, numberOfClusters):\n",
    "            dijAux = distanceAux.index([pointInd, clusterInd])\n",
    "\n",
    "            ANDs[clusterInd] = DJC.AND( DJC.term( di >= dijAux ),\n",
    "                                        DJC.term( label == clusterInd ) )\n",
    "\n",
    "        M.disjunction( [ ANDs[i] for i in range(0, numberOfClusters) ] )\n",
    "\n",
    "    # symmetry breaking constraints\n",
    "    M.constraint( centroids[:-1,0] <= centroids[1:,0] )\n",
    "\n",
    "    ########## solve ##########\n",
    "\n",
    "    M.objective( ObjectiveSense.Minimize, Expr.sum(distances) )        \n",
    "\n",
    "    M.setLogHandler(sys.stdout)     # Enable log output\n",
    "\n",
    "    # Solve \n",
    "    M.solve()\n",
    "\n",
    "    # get centroids and clusters\n",
    "    labels_EUCL = Y.level()\n",
    "    centroids_EUCL = np.array(centroids.level()).reshape(numberOfClusters, d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA070lEQVR4nO3deXxU9b3w8c93JpNkskMSAtlYBILsIItWBVcWodStrbW11t6r3Wy9PpW2PtpKF6/28bFWb3uvtU9vta2t9iqlFkTUClWpgCwqSwiLbAlbSDLZt5n5PX+cM8MkZINMMjPh+369eJE563cmZ775ne/5nd8RYwxKKaVilyPSASillOodTeRKKRXjNJErpVSM00SulFIxThO5UkrFOE3kSikV4zSRh4mIjBARIyJx57j+50Xk9XDH1RdE5FkR+Umk44hFIvIlEXk30nEAiMhOEbki0nF0R0TWici/djKvV9+7gUITeTsicpmI/FNEqkWkUkTWi8jMMO/jjIPPGPO8MWZeOPdj7+sKESntxfp9mnjaf0nteKtE5JYOlj0oIi0iktVu+gf25zmir+KMNvb7Hd2bbRhjJhhj1oUpJBVBmshDiEgasBL4D2AwkAf8EGiOZFznCxGZB6wAvmyMeaGTxQ4AnwtZZxLg7vvoBo5oar2KRfNQL+kH2NZYAGPMn4wxPmNMozHmdWPMRwAi4hCRB0XkkIicFJHfiUh6RxuyW4/XhLxeJiJ/sF++bf/vEZE6EbmkfctXRD4hIu/bZwbvi8gnQuatE5Ef22cLtSLyevtWqr1cMrAayLX3UyciuSKSICI/F5Gj9r+fi0hCB+tfCDwNXGKv6wmZPUhEVtn73ygiF4SsN05E3rDPaEpE5DPdffAishj4M3CrMeYvXSz6e+CLIa9vB37XblsJIvJ/ReSwiJwQkadFxG3PGyQiK0Wk3G75rxSR/JB1O/1sRSRRRP4gIhUi4rF/LzmdvJ8CEVlu76dCRH7RwTJnnJmFnqGIyGgR+Yd9DJwSkRft6YHj50P79/LZwGdon514xDqrnByy3YMi8l0R+QioF5G40GPUPj7/bB/TtWKVXWaErD9dRLbZ8/5HRF6UTspr9rG8XkT+w459t4hc3e49Piwi64EGYFRXx7vtAhHZZM//q4gM7mTf6SLyGxE5JiJlIvITEXG2i+sJ+zP62N7vl0TkiFjf6ds72m7UM8boP/sfkAZUAM8BC4FB7eZ/GdgHjAJSgOXA7+15IwADxNmvDwLXhKy7DPhDR8va074EvGv/PBioAm4D4rBaoFVApj1/HbAf6w+P2379aCfv6QqgtN20HwEbgCFANvBP4MedrB+MK2Tas0AlMMuO73ngBXteMnAEuMOeNx04BUzoZPvrgL8CntDPq5NlDwLXACXAhYDT3tdw+/McYS/3c+AV+3NMBf4GPGLPywRuApLsef8DrGgXT4efLfAVe1tJ9r4vAtI6iNMJfAg8YX8eicBlHfyeOzoO1gH/av/8J+ABrAZXcBv2PAOMDnk9HTgJzLb3f7v9eSWEfHYfAAWAu/0xinV8NgHX2es/Amyw58UDh4B7ABdwI9AC/KSLY8YL3Gsv/1mgGhgc8h4PAxOwjpEcuj/ey4CJ9uf5Mp18l7DO6H5lLzcE2AR8pV1cd9jv8Sd2HL8EEoB5QC2QEulcdLb/tEUewhhTA1yGdWD8GigXkVdCWl2fB35mjPnYGFMH3A/cIuE/VV0E7DXG/N4Y4zXG/AnYDXwyZJnfGmP2GGMasVqyU89i+58HfmSMOWmMKccqH912ljEuN8ZsMsZ4sRJ5YP+LgYPGmN/asW/F+uLd3MW2rgT2AOt7uO9Aq/xarM+lLDBDRAS4E7jXGFNpjKkF/h24BcAYU2GMedkY02DPexiY2277nX22rVh/CEYb64xti33MtDcLyAWWGmPqjTFNxphzuc7QivVHKrcH27gT+JUxZqMd23NYJcGLQ5Z5yhhzxH5fHXnXGPOqMcaH9RlPsadfjJVgnzLGtBpjlmMlyK6cBH5uL/8i1h/fRSHznzXG7LSPn3l0f7z/3hizwxhTD3wf+EygpR1gf08XAv9mf+4nsf6Yhl5vOWAfmz7gRaw/bD8yxjQbY17H+gPVq2sPkaCJvB1jTLEx5kvGmHysFkAuVgsP++dDIYsf4nSLIpza7yewr7yQ18dDfm7AOkM41+0fsqedjc72PxyYbZ+6euxyzOeBoV1s6/tYSWeF2CUeEVktp8tBn2+3/O+BW7FaWL9rNy8bq8W8JWT/r9nTEZEkEfmVWOWxGqwyV0a7pNDZe/s9sAZ4QayS1P8REVcH76cAOGQnqd74DiDAJrvU8eUulh0OfLvd515A29/rkW721/59J9qNlFygzNjN2h5uq/3y7Y+x0PV7crwfaTfPBbQvJw63px8L+Qx+hdUyDzgR8nMjgDGm/bSz+S5FBU3kXTDG7MYqI0y0Jx3FOlgCCrFO1U5wpnqshBIQmsi6G3Ky/X4C+yrrYNnudLSvjt7H0bNYvytHgH8YYzJC/qUYY77WxTr1WKf06cBLIuIyxiy010sxxjzfJiBjDmFd9LwOq7wV6hTWl3FCyP7TjTGBL+e3gSJgtjEmDZhjT5fu3pjduvyhMWY88Amss48vdrDoEaCwB2dq9fb/HR4nxpjjxpg7jTG5WGWd/5TOe6ocAR5u97kn2a3b4Ca7iaczx4A8+2wnoKCbddov3/4YC42lJ8d7Qbt5rVi/61BHsBoEWSGfQZoxZkI3scY8TeQhxLpI922xL36JSAFWvW6DvcifgHtFZKSIpGCdsr/YScvrA6yyi8u+aBRaWigH/Fi19o68CowVkVvti1KfBcZj9ag5WyeATGl7UfZPwIMiki3WhbwfAH/ocG1r/XwRie/h/lbasd9mv3eXiMwU68Jpp+wyxwKs1tkf2582d+BfgKvsU+3Q7fixymJPiMgQABHJE5H59iKpWIneY18we6iH7wsRuVJEJtmx1WAlE18Hi27CSn6PikiyWBdJL22/kF3WKgO+ICJOu8UdetH403L6QmwVVvIL7O8EbY+fXwNfFZHZYkkWkUUiktrT99eF9+z93m0fj5/CKh91ZQjwLfv3/2msaxqvdrJsT473L4jIeBFJwrrG85JdHgkyxhwDXgceF5E0sTonXCAi7UtnA44m8rZqsS4WbRSReqwEvgOrFQfw31in129jtQibgG92sq3vY30pq7Bq0H8MzDDGNGDVZtfbp4ChdUyMMRVYrb1vY118/Q6w2BjTvgXSLfus4k/Ax/a+crEu8mwGPgK2A1vtaR15C9gJHBeRbvdvJ+R5WHXJo1in6z/FupjU3boerLr3WOB30kW3NGPMfmPM5k5mfxfrovQGu3zyJlYrHKwymRurNbcBq+zSU0OBl7CSeDHwDzr4A2gnmE9i1VoPA6VYF/w6ciewFOv3PAHrwnPATKxjsQ7r4u09xpgD9rxlwHP27/Qz9mdxJ/ALrGNuH1bpqdeMMS1YFzj/Beui9BewkmxX3XI3AmOwPueHgZvt47qj7ffkeP891tnxcawLv9/qZL9fxLo4uwvrc3gJGNbNW4x50raMpZRS3RORjcDTxpjfdjDvS1g9by7r98DOU9oiV0p1S0TmishQu/RxOzCZszubUX0oau7wUkpFtSKsrpgpWP3sb7Zr0ioKaGlFKaVinJZWlFIqxkWktJKVlWVGjBgRiV0rpVTM2rJlyyljTHb76RFJ5CNGjGDz5s56jimllOqIiLS/AxbQ0opSSsU8TeRKKRXjNJErpVSM037kSqmwam1tpbS0lKampkiHErMSExPJz8/H5epocM0zaSJXSoVVaWkpqampjBgxgrYDIKqeMMZQUVFBaWkpI0eO7NE6mshVzFuxrYzH1pRw1NNIboabpfOLuH5aXvcrqj7R1NSkSbwXRITMzEzKy8t7vI4mchXTVmwr4/7l22lstUY0LfM0cv/y7QCazCNIk3jvnO3npxc7VUx7bE1JMIkHNLb6eGxNSYQiUqr/aSJXMe2I5xg+qs+YftTT2WMpleqex+PhP//zP896vc2bN/Otb3U8VPqIESM4deqsHynQI1paUVGrq9r3iboTbDm2BVfqFurqB+H2T2+zbm6GOxIhqwEikMi//vWvnzHP5/PhdHb8AKsZM2YwY8aMvg7vDNoiV1EpUPsu8zRiOF37fva9D1i1ZxV/LfkrpxpOcfelCxjknNZmXbfLydL5RR1vWJ0Xfve73zF58mSmTJnCbbfdRnl5OTfddBMzZ85k5syZrF+/HoBly5bx5S9/mSuuuIJRo0bx1FNPAfC9732P/fv3M3XqVJYuXcq6deu48sorufXWW5k0aRJNTU3ccccdTJo0iWnTprF27VoA1q1bx+LFiwGoqKhg3rx5TJs2ja985SsERpqtr69n0aJFTJkyhYkTJ/Liiy/2+v1qi1xFpUDt2+AFnAhCY6uPn699n28vSuTi/IsZnz2eOEccI9O110o0+1vJ386YdsHgCxifPR6v38vqvavPmF+UVcTYzLE0eZt4Y/8bbeZ9suiTXe5v586dPPzww6xfv56srCwqKyu5++67uffee7nssss4fPgw8+fPp7i4GIDdu3ezdu1aamtrKSoq4mtf+xqPPvooO3bs4IMPPgCsBL1p0yZ27NjByJEjefzxxwHYvn07u3fvZt68eezZs6dNHD/84Q+57LLL+MEPfsCqVat45plnAHjttdfIzc1l1apVAFRXn1kaPFuayFVUKvPU0SKHaHJsR0wCiWYSLpOHpyaHz026jjjH6UP3+ml5mrhV0FtvvcXNN99MVlYWAIMHD+bNN99k165dwWVqamqora0FYNGiRSQkJJCQkMCQIUM4ceJEh9udNWtWsF/3u+++yze/aT2ud9y4cQwfPvyMRP7222+zfPny4D4GDRoEwKRJk7jvvvv47ne/y+LFi7n88st7/Z41kauo4jd+dpXvQlLWUNt0AHDiogCMVQXMy0hpk8RV9OuqBR3niOtyfmJcYrct8PaMMWd03/P7/bz33nu43WdeO0lIOP1ccKfTidfr7XC7ycnJbfbREx11Ixw7dixbtmzh1Vdf5f7772fevHn84Ac/6NH2OqM1chVVBOHPO//MhQX1JDlGkuZbRKpvAS4zTGvfqkeuvvpq/vznP1NRUQFAZWUl8+bN4xe/+EVwmUDJpDOpqanBFntH5syZw/PPPw/Anj17OHz4MEVFRZ0us3r1aqqqqgA4evQoSUlJfOELX+C+++5j69atZ/0e29OmjYoov/FTcqqEfxz6B1+c/EWS4pO4bfJt3DrJx67DyTz++l6tfauzMmHCBB544AHmzp2L0+lk2rRpPPXUU3zjG99g8uTJeL1e5syZw9NPP93pNjIzM7n00kuZOHEiCxcuZNGiRW3mf/3rX+erX/0qkyZNIi4ujmeffbZNyx7goYce4nOf+xzTp09n7ty5FBYWAlZdfenSpTgcDlwuF//1X//V6/cckWd2zpgxw+iDJc5vfuNnb8VeVu9bTfGpYpzi5PYptzM7f3akQ1O9VFxczIUXXhjpMGJeR5+jiGwxxpzRv1Fb5KrfNXubeXrz01YtXIQpOVOYd8E8Rg0aFenQlIpJmshVvzDGUNFYQVZSFglxCTT7mpmcM5kFoxcwatAoHZtDqV7QRK763MdVH/O3kr9RfKqYH17xQ3JScvjW7G+R4EzQBK5UGGgiV33moOcgr5S8wkcnPkJEmJozFadYtzYnxiVGODqlBg5N5KpPeBo9PPz2wxgM04ZNY9GYRYzIGBHpsJQakDSRq7A5XH2YD45/wJKiJWS4M/j0hE8zLmschemFkQ5NqQFNE7lq41yetlNaU8qK3SvYemwrgjA7bzY5KTnMu2BeP0WtVOeWLVtGSkoK9913X4fzV6xYwdixYxk/fnw/RxY+mshV0Nk+baeqsYrntz/P1mPWnWkXDbuITxV9ipyUnP4LWqleWrFiBYsXL47pRK636KugwIiDPjwYWoCOn7bj9VtjUTjFyY6TO7ho2EX86Mof8Y1Z3yA/Pb/f41aqvYcffpiioiKuueYaSkqs4/fXv/41M2fOZMqUKdx00000NDTwz3/+k1deeYWlS5cydepU9u/f3+Fy0U5b5ANIbx9CHHiqjpOMDqcfrzvOy8Uvc7TmKA9f/TBpiWk8Pu9xkuOT229KqaD+HsZ2y5YtvPDCC2zbtg2v18v06dO56KKLuPHGG7nzzjsBePDBB/nNb37DN7/5TZYsWcLixYu5+eabAcjIyOhwuWimiXyACMdDiLPTWzhW3XBGIs9Ka+CX7/+S98veB2BG7gyavE0kxiVqEldR55133uGGG24gKSkJgCVLlgCwY8cOHnzwQTweD3V1dcyfP7/D9Xu6XDQJSyIXkQzg/wETAQN82RjzXji2rXom9CHELbIfwY1pzeKxNSU9TuRXT6nkuQ0lOFouR+yqm7gOkZKzlk1lSczMnckN424gL00HrlI919/D2ELHw8d+6UtfYsWKFUyZMoVnn32WdevWdbhuT5eLJuGqkT8JvGaMGQdMAYrDtF3VQ4Hyh8FPi+MwTY6PqHO+xd6at9h8dDMNrZ3U+exB0/ZX7mdkjo8HrllEdloLLXxMXoabx65fxOenzeeRqx/h7ll3axJXUW/OnDn85S9/obGxkdraWv72N6u0U1tby7Bhw2htbQ0OLwtnDlnb2XLRrNctchFJA+YAXwIwxrSAfaVM9ZvcDDdlnkYEB8m+K/BTg1fKyUipYduxbYzLGgfA0dqj1LXUUZBWgPvhn4LHg+/x/8umsk04xYF35TImJlVxyZzLeXLBN3A4HMC/RvbNKXUWpk+fzmc/+1mmTp3K8OHDg0/g+fGPf8zs2bMZPnw4kyZNCibvW265hTvvvJOnnnqKl156qdPlolmvh7EVkanAM8AurNb4FuAeY0x9Z+voMLbh175GDtZDiB+5cRLXTc4m3hkPwLqD69hTsQeMIevFv1Hwx1Uc+txC/nJJBi1bNuEu3s+sMVdw44N/IFdb3+oc6DC24dHfw9jGAdOBbxpjNorIk8D3gO+3C+Au4C4gOMC6Cp9AHby7Xitzh8+lKLOIbce38dGSS/lN/E4OVq0mdwUsKYGbPvEv5P7016CDWSkVM8KRyEuBUmPMRvv1S1iJvA1jzDNYLXdmzJjR/0+zGOC66nroN34cYl0O2VC6gZ3lOznVcIpd5cU0jBnORX89yLK1MKoa2KhJXKlY0+tEbow5LiJHRKTIGFMCXI1VZlH9pH1ZpdTTwNLl73KwejBFeV5O1J/g1km3Eu+MJ9mVTH1LPXUtdbhdbj6zy8EXV0Fyq72xe++FJ57QZK56paMHIKueO9uSd7h6rXwTeF5EPgKmAv8epu2qHgjteuiVk9Q736LS/y7/sX41tS21jBk8Jng35rjscSTEJeByxHHl24e5/edrSf76PeD3wz33wJNPWsk8Ao8AVANDYmIiFRUVZ52MlMUYQ0VFBYmJPR/qOSz9yI0xHwBnFOBV/wh0PQRwmCTiTDZOk4m3NpPPTLiJg56DvHPoHa694FrinfF8cuwn2Vi2kQXJrSR9veh0C/yJJ6yNZGRoi1yds/z8fEpLSykvL490KDErMTGR/PyeD3ehD18eAC599C3KQpJ5QE66j4dujqO0ppTB7sHMGT6HIclDgJBTX2PaJu32r5VSUaOzXis6aNYAsHR+EW6XM/ja4APXPmZcWMLJ+pN8ouATzM6bzco9K9lXuQ8IufOtfdLWJK5UzNGxVgaAM7oepidx2eRErpswhdl5s6lsrOS1fa+RnphOXqr2DVdqoNFEPkDMHZdEarrh8uFWHdzrv4I4h1VWWbNvDemJ6Sweu1iflanUAKSJPMa1+FrYemwrO07uIM4RR2VjJUNThhLniKOupY41+9aQkZjBorGLNIkrNUBpIo9h+yr3saF0Aw2tDYzLGsesvFltknVKfAqXD7+cwvRCTeJKDWCayGPY3oq9JLuSmXfBvGBvFLAegpzgTCAnJYexmWMjGKFSqj9oIo8hLb4WthzdwsQhE0lNSOWqkVcR74xvcwfdIc8h3vj4DYalDGPR2EURjFYp1V80kccAYwx7K/eysXQjTd4mBrkHMS7BukMz1EHPQd78+E0y3Zlce8G1EYpWKdXfNJFHuVMNp1h/eD0n6k+Qk5zDwjELyUrKOmO5QBLPSsriujHXBYetVUoNfJrIo9zuU7upaa7hihFXMGbwmE4HIvq46mOyk7JZOGahJnGlzjOayKOMMYbdp3aTmZTJkOQhzMqbxay8WZ0m58AQtVeMuAKf34fL6erniJVSkaaJPIqcrD/J+sPrKW8oZ3z2eIYkD+mydb2/cj/bjm9j0ZhFuF1uHE4dcUGp85Em8ijQ5G1iU9kmdp/aTZIriatGXsXowaO7XGdf5T7WHlgbvPlHKXX+0gzQzzp6ks+o3Ar2VOxhSs4Upg+b3m15ZG/FXtYdXMew1GEsGL1AE7lS5znNAP0o9Ek+Pqo4VF3O/ctbePiGCdw8/mYyEjO63caBqgOsPbiW3NRcTeJKKUCHse1XgSf5+KiiwbmRFtlHY6uPx1/f26MkDpCTksOFWRdqEldKBWki70eBJ/k4GUSCfzxu/6w207tSVlOG3/hJciVx+fDLNYkrpYI0kfej3Ax38Od4U4jgPGN6R3af2s2qvavYfmJ7n8anlIpNmsj7Ufsn+QC4XU6Wzi/qdJ3i8mLePvQ2BWkFTBwysa9DVErFID0/70dnPMnH7rUSmN7ervJdvHv4XQrTC7l21LU4Hc4Ol1NKnd80kfez66fldZq4QwX6lmsSV0p1RxN5lEqMS2RJ0RLSE9I1iSuluqSJPAqE3iSUnnacOy4dzj1zr2Gwe3CkQ1NKxQC92BlhgZuEyjyNNMtBjtR/wON/38SKbWWRDk0pFSM0kUdY4CYhQwstjv3EmaFIyyQeW1MS6dCUUjFCSysRFrgZSIgnyXcxghvB0aObhJRSCrRFHnGhNwM5SEbsX0l3NwkppVSAJvIIO5ebhJRSKlTYSisi4gQ2A2XGmMXh2u5Ad7Y3CSmlVHvhrJHfAxQDaWHc5nmhpzcJKaVUR8JSWhGRfGAR8P/CsT2llFI9F64W+c+B7wCpnS0gIncBdwEUFhaGabdKKXV2OnpKV6yfEfe6RS4ii4GTxpgtXS1njHnGGDPDGDMjOzu7t7tVSqmzFnoDngHKPI3cv3x7zN+AF47SyqXAEhE5CLwAXCUifwjDdpVSKqxOP6WrGh9VADS2+mL+BrxeJ3JjzP3GmHxjzAjgFuAtY8wXeh2ZUkqFWeBGuxbHPpocO8+YHqu0H7lS6rwRuNHO0IqQcMb0WBXWRG6MWad9yJVS0SpwA56RZgQXMDBuwNOxVpRS541A75Rlq49wqtZB3gDptaKJXCl1XrFuwLsz0mGEldbIlVLnFb/x09DagN/4Ix1K2GgiV0qdVyobK/nDR3/gcPXhSIcSNprIlVLnlSZvE2A9F3eg0ESulDqvNLZafcbdcbHd5TCUXuxUSkVcf45/MhBb5JrIlVIRFRj/pLHVB5we/wTok2Te6G3EIQ4S4hK6XzhGaGlFKRVRoQ8g91MP9O34J4XphVycf3GfbDtStEWulIqowDgnzY79tMph4v0jiTcX9Nn4J0NThjI0ZWifbDtStEWulIqowDgn8f6RuMwwWhwf0+B8h8HpVX2yP0+Th4bWhj7ZdqRoIldKRVRg/BMHiST6J5Pkm028M57Z48rYemxr2Pf32r7X2FC6IezbjSQtrSilIqr9A8gLM3L59ry5jM2vZkTGCAAaWhtwOVy4nK5e76+xtXFA9VgBTeRKqSjQ8QPIC4I/rTu4jqrGKi4puIRRg0ad8358fh+t/tYB1YcctLSilIoBFw27CLfLzZsfv8nKPSupajy3+nmj174ZyKWJXCml+lVOSg43jLuBywov41TDKV4ufplDnkNnvZ2BeDMQaGlFKRUjRITx2eMZNWgU245tY1jqMMCqefe0hZ0Sn8KVI64kO2lgPQBeE7lSKqYkxiVyScElgDUk7co9K0mMS+TSwksZ7B7c7bpjMsf0R5j9SksrSqmYJQgTh0yksrGSl3e9zHtH3qPF19Lp8jXNNRyvO44xph+j7HuayJVSMUtEuDD7Qj478bOMyxrH9pPbeXHHi3iaPB0uv/vUblbuWYmI9G+gfUxLK0qpmJcYl8jlwy9nXNY4dpXvIj0hHYBWX+vpvufGWPX0QNdDY2CAJHRN5EqpASM7OZu5yXMBaPY28z+7/oeRGSOZ+d+vEe+ppfEb860eK8bAvfdCRgYsWxbRmMNBSytKqQFJRBiZMZJd5Tt5oWETJX94ksZnfok7LtFK4k8+CR6PldRjnLbIlVJ9qj8eGuH1e/H6vfiNnyRXEgC1zbUUpheSEJfAppu9/Emq2V+yillfWMWejyH+3jtwPXQf8Q2nSEtIC/v45P35sAxN5EqpPhP60AiDlyOequBDI64an0pNc00wCbf6WvEZHxOHTARgb8VejtUdo9XXas33t+IUJwvHLASs2/YPVB3A6/disFrVKfEp3DrpVgA2lm2ktKYUAKfDSfUVF/NBzUbK0qAiCXJvnQ37XgXgyhFXMiZzDMfrjrNyz0pcDhfxznhcTuv/WXmzGJoyFE+Th5JTJcHpgX9DU4aSGJeI1++lxdfC6o/KeXBFcb89LEMTuVKqzwQeGgHQ5NiBV44jrdfw2JoScrLcbD+5/Yx1JmRPQEQ41XCKI9VHiHPEEeeIw+V04XKcHjRraMpQEpwJwXlxjrg2d2xenH8xXr+XOEccNU3V+P/8IrIfErwwuBEyX3yFGf/7l4jDEex/nuRKYkrOFFp8LbT6W2nxtdDia8EhVhW6uqmaHSd34DO+NjEvKVrC0JShHKg6wNqDa/npa7vx+L2IM44k3ydw4A4+LKMvErlEoj/ljBkzzObNm/t9v0qp/jXye6sIZBgvpzDSSJzJxYGTD394eXBUw9BkHO+MD2sMjS0N/OWBGzFr1jBo3hI8X/wME557lffX/ZGkq+ZzzUPPMSQl56y2GRh8K5Do0xPScTldVDdVU1ZbxjU/+zuGVox4SfAXIXabWYADjy465/ciIluMMTPaT9cWuVKqz+RmuCmzn/QTRxaBrJ6b4SYtIY20hLQ+2W+gPl3mqced+iFXJw5i2dV3cfDef6Xs+AdMeux3DPt2Mn9PrqCmpZYhnF0idzqcOB3OM8ZsSU9MJz0xnZHpx6333a6dHHiIRrhprxWlVJ8JPDQilNvlZOn8oj7bZ6AuX+ZpxOCjst7Laq7l3du+T2pCGgZDXWs9Q372Kz79gxcZPXg0AIerD9PYGp7Hy/X3++51IheRAhFZKyLFIrJTRO4JR2BKqdh3/bQ8HrlxEnkZbgTIy3DzyI2T+qz3BrStywsu3P6Z+FpzeOz1PWQkZpCdlI3X7wUR4hxWUaLF18JbB97i5eKXOVp7tNcx9Pf77nWNXESGAcOMMVtFJBXYAlxvjNnV2TpaI1dK9ZXQunyo7urTlY2VvPnxm3iaPEwfNp3pw6YHL3JGi85q5L2O0hhzzBiz1f65FigG+u7PrVJKdaGzOnR39enB7sHceOGNjM0cy9ZjW1m1ZxV+4++LEMMurBc7RWQEMA3Y2MG8u4C7AAoLC8O5W6WUClo6vyjYdz0gUJ9esa2MB1b/kco6P6PSZp5xk06cI44rRlxBXmoeNc01Udci70zYErmIpAAvA/9mjKlpP98Y8wzwDFillXDtVymlQrV/mHPgrkqA+5dvp9JXj6G1y5t0QscsP153nANVB5iVNwuno+0FzGgRlkQuIi6sJP68MWZ5OLaplFLnqqOHOV/66Fs0tvqIk6H4aQbo0U06R2uPsv3kdo7VHeOaUdf0WZfJ3ghHrxUBfgMUG2N+1vuQlFIq/I7a/dldpoAEM/qM6Z2ZPmw68y6YR21zLS/vepl9lfv6NM5zEY4C0KXAbcBVIvKB/e+6MGxXKaXC5lwvggKMyBjBTeNvIjMpk7cOvMXh6sPhDq9Xel1aMca8i9WzRymlolZXF0F7IiU+hcVjF7O3Yi8FaQWAdat+NNTNY+OSrFJK9VI4btJxiIOirCJEhPqWel7Y8QLF5cV9F3QP6VgrSqnzRkcXQc+VQxwMcg/incPvUFZbxpzhc8I+4FePY4nIXpVSKsa5XW4Wjl7IrLxZHKg6wPLi5ZTXl0ckFk3kSil1jkSEqUOnsqRoCX7j73B89f6gpRWllOqlnJQcbrrwJqze2FDTXEO8M/6MYW77iiZypZQKg9Bnfr514C3qWuq4auRV5Kbm9vm+tbSilFJhdnnh5bgcLlbtWcWWo1vo6yexaYtcKaXCLDMpkxsvvJF3D7/LlmNbOFp7lMba8Tz55qE247+EqweNJnKllOoDLqeLK0deSV5aHr/d+C5/eW83Ta3WvK4G7DoXWlpRSqk+NDZzLG9/NCKYxAMCA3aFgyZypZTqY50NzNXdgF09pYlcKaX6WG8G7OoJTeRKKdXHls4vwu1qO7jW2QzY1R292KmUUn2ss6cWaa8VpZSKIeEcsKs9La0opVSM00SulFIxThO5UkrFOE3kSikV4zSRK6VUjNNErpRSMU4TuVJKxThN5EopFeM0kSulVIzTRK6UUjFOE7lSSsU4TeRKKRXjNJErpVSM00SulFIxLiyJXEQWiEiJiOwTke+FY5tKKaV6pteJXEScwC+BhcB44HMiMr6321VKKdUz4WiRzwL2GWM+Nsa0AC8AnwrDdpVSSvVAOBJ5HnAk5HWpPa0NEblLRDaLyOby8vIw7FYppRSEJ5FLB9PMGROMecYYM8MYMyM7OzsMu1VKKQXhSeSlQEHI63zgaBi2q5RSqgfCkcjfB8aIyEgRiQduAV4Jw3aVUkr1QFxvN2CM8YrI3cAawAn8tzFmZ68jU0op1SO9TuQAxphXgVfDsS2llFJnR+/sVEqpGKeJXCmlYpwmcqWUinGayJVSKsZpIldKqRiniVwppWKcJnKllIpxmsiVUirGaSJXSqkYp4lcKaVinCZypZSKcZrIlVIqxmkiV0qpGKeJXCmlYpwmcqWUinGayJVSKsZpIldKqRiniVwppWKcJnKllIpxYXlmZ39Ysa2Mx9aUcNTTSG6Gm6Xzi7h+Wl6kw1JKqYiLiUS+YlsZ9y/fTmOrDz9NlHr83L98O4Amc6XUeS8mSiuPrSmhsdUHQKNjM3XO1znlX8cDq5+n2dsMQJO3Ca/fG8kwlVIqImKiRX7U0xj8Od6Mxm9q8Esd5XVVuJwuALYc3cLO8p2kxqcyyD2IjMQMBiUOoiirKFJhK6VUv4iJRJ6b4abMTuYuMxQYCgbyMtw4xDqpGDVoFG6Xm6rGKqqaqiirKcPtcgcT+duH3qa2ubZNkh/kHkRiXGKk3pZSSoVFTCTypfOLgjXyALfLydL5p1vbw1KHMSx1WPC1MYZG7+mWfGJcIhUNFZScKqHV3wpATnIOnxr3KQC2HtuKy+FikHsQgxIHkRyf3NdvK2bphWeloktMJPJAkjib5CEiJLmSgq9n5c0Ce/G6ljo8TZ42y+8+tZu6lrrga5fDxYQhE6z1gCPVR0hPTCc1PhURCdM7iz2hF54ByjyNeuFZqQgTY0y/73TGjBlm8+bN/b7f7jS2NuJp8lDVVEVVYxXZydmMzRxLs7eZ5z58DgCnOMlIzCAjMYOirCLy0/IxxmAwwTLPQHbpo29R5mnEhwe/1BNnchGEvAw36793VaTDU2pAE5EtxpgZ7afHRIu8v7hdbtwud5sSDYDL6eL6cdcH6++eJg8n60+Sl2a1QKuaqlhevJy0hLQ29ffc1Nw2ZwUDQeDCs1eO0+o4jMuX12a6Uqr/aSLvAYc4GJI8hCHJQzqcH++MZ0rOlGCSP+Q5hMGwcPRCktKTOFZ7jA9PfNgmyWckZhDvjO/nd9J7gQvPfmlGSGgzXSkVGb1K5CLyGPBJoAXYD9xhjPGEIa6YkhKfwsy8mcHXPr+PmuYaUuJTAGj1t1LXUkdZTRk+c/qC7afHf5pB7kGcqDtBRWNFTPSkCVx4bvVmYUwacOaFZ6VU/+pti/wN4H5jjFdEfgrcD3y392HFNqfDySD3oODrwvRCCtMLMcZQ21IbLNGkJViJ8KDnIB+e+DC4fGJcIhmJGSwaswinw0lNcw0OcQT/METS6QvP8dprRako0atEbox5PeTlBuDm3oUzsIkIaQlppCWkMZzhwemz82czYcgE60KrneQbWhtwOpwAbCrbxMdVHwe7R2YkZpCVlMXEIRMj8j6un5bHtRMycLvc58UFXqWiXThr5F8GXuxspojcBdwFUFhYGMbdDgwp8SmkxKeQn5Z/xrwpOVPITc0NJvnSmlIqGyuDiXzlnpU0e5utGrzdDz4zKTPY4u/KufQJb/W18vz255mdN5spQ6ec2xsOE+3TrlQPErmIvAkM7WDWA8aYv9rLPAB4gec7244x5hngGbC6H55TtOep7ORsspOz20zz+U/X2oelDKO8oZyT9SfZX7UfgOHpw5k/ej5g3dXqjnMHW/MZiRnEOeLOrU+4MTS0NgBYPXKMgQj1q9c+7UpZuk3kxphrupovIrcDi4GrTSQ6pZ+nAmUXgItyLwr+7PV729zs5Dd+jtcdp7qpGsPpX8+0odN4bE0tDa2teOUoLXIQB8m0eDP599caWTL1U2eWTZYtA4+Hhh99B4CkODfcey9kZFjz+oAxhhZfCy2+FuIccbhdbrx+L/sq97Fs9Tt4vHUYh5cE/4UIQmOrj8fWlGgiV+eV3vZaWYB1cXOuMaYhPCGp3ohzxJGVlBV87RAHn5nwGfzGT3VTdbCL5JDkIRz1nMTQSJNjB61yGL80AMLu+gRe25fAdWOuo6m1iQOeAxSk5ZPi8cCTT1Lv8MCtl5C87GF48mm4554uW+Zev5cWXwtAsF/9/sr9NHmbaPY10+JrodnbTE5KDuOyxmGM4YUdLwTnBUzOmczF+RdjjOHtQ29TWr8d4wDBRQJjCRzO2qddnW96WyP/BZAAvGHftr7BGPPVXkelws4hDqt+HtKbJjdjD6UeQ7JvDn6pw08dXqkgJclDQVoBADvLd/LEhidwOV1kLEgnP+Ey/O88x9Htz3HwI5D/dSvN936OliPrSYlPYerQqYBVt69qrKLZ14zf+AEYkTGCeRfMA2D9kfU0eZsA649PgjMh2O1SRMhLyyPOEUe8M54EZwIJcQkMdg8GrBu0bp10K799/Z8cq249471qn3Z1vultr5XR4QpE9b/Tg5El4TBJwBDS48bwyKJJTMqxShNjMsdw5/Q7Oeg5yJGaI+yfM5GSY+9Sngxb86DgEzBq7yry0/IZNWhUcNtZSVmkJ6RbiTgugQRnAumJ6cH5N4y7AZfTRbwzvsOeL3OGz+ky9pT4FL67YEK3g6kpdT7QOzvPYz0ZjCwtIY25I+Yyl7lgDM33fpOnd0N9HLS44PDmbRya5aSxtZHC9NO9kWblzeqya2JqQmq/xK/U+UAT+Xnu+ml5PUt8xsC997LnT7/EfftVfOGhF8l84CeUP/UkG782mi2LRwVLI54mDy/teon8tHwK0wvJT8vvUVfIPo1fqQFME7nqGRHIyKD4iwsZcsedZCZnwRNPkA0sTsng2svuJ85hHU67T+3moOcgx+uOs6diD4lxiaQnpHPlyCs7Ha9GKXXuNJGrHjv67bvwlAzjiiETrAki8MQTIKHDZ1kXNecMn8MhzyEaWhvwGz9N3ibW7qrhqb/v4FD1XtJSKvnKZbO44+KZfdZaV+p8oYlcdStw9+QRz1FS0jxkzk9gbKY9s4Muh0NThjI0ZSj1LfUUnypm96nd7D7q5+X1+2hs9eEXH+UNHh79+18pqdzMnDHDKUwv5JL8S87rh3Yoda40kasuhd496WAQDTWDeHBFMU5HXLe16eT4ZGbkzmD6sOl8YtVqK4nThMOkWF0efQ28V1zLkkkZVDRUBJP45qObSXAmUJBeQEZiRj+8S6VimyZy1aXH1pTQ2OrDRxWCGweJZ333pEMcnKg+/dqJ1R/cQRKemiQWjF4QnGeM4ZDnEBWNFbxX+h6p8akUpBcwevBohqZ0NFKEUkoTuerS6bskHXjlKPFmVLvpPRN4IIWDxDOmhxIRbhp/E7XNtRypOcKR6iPsqdhDkiuJoSlDafW1UnyqmML0Qm2tK2XTRK66FEjATtJxmvQ208/G6ZuPenbzTmpCKuOzxzM+ezw+vy/4QI6T9SfZULqBDaUbSIlPoSCtgIL0AvLT8oO9ZpQ63+hg0qpLS+cX4XY520w7l7snr5+WxyM3TiIvw40AeRluHrlxUo/KM06HM/hYvLy0PG6ddCuXF15OVlIW+yr38fr+16lusmo31U3VVDZWnlVsSsU6icSAhTNmzDCbN2/u9/2qcxPNY377jZ8TdSeCD8x++9Db7D61m2RXMgXpBRSkFZCXlheTz0dVqj0R2WKMmXHGdE3kaiCpb6kP1tbLasto8bWQlpDGLRNvCc5Pjk+OcJRKnZvOErkWFdWAkhyfzLiscYzLGhdsrTf7mgGrR8zy4uWISJvaurbWVazTRK4GLIc4giUXAINhVt4sjtQc4YDnACUVJQjC7PzZTM6ZbC1jjN6UpGKOJnJ13nCIg6KsIoqyivAbPyfrT3Kk+khw/JeT9SdZs29NsLaen5ZPQpw1+EA0XydQShO5Oi85xBEcSiDAKU5yU3M55DnEnoo9CMKQ5CE01I7jR68c0GeDqqiliVwpW2ZSJlePuhpjDOUN5RyuPszR2qP8x5uHaWz14eUUcViP0dNng6pooolcqXZErJZ4oORyT/UqazquNsvps0FVtNAbgpTqRuAuVifpHU5XKtI0kSvVjXDd3apUX9HSilLd0GeDqminiVypHtBng6popqUVpZSKcZrIlVIqxmkiV0qpGKeJXCmlYpwmcqWUinERGY9cRMqBQ2HebBZwKszbDJdojg2iO75ojg2iO75ojg2iO75ojW24MSa7/cSIJPK+ICKbOxpwPRpEc2wQ3fFFc2wQ3fFFc2wQ3fFFc2wd0dKKUkrFOE3kSikV4wZSIn8m0gF0IZpjg+iOL5pjg+iOL5pjg+iOL5pjO8OAqZErpdT5aiC1yJVS6rykiVwppWLcgEvkInKfiBgRyYp0LKFE5Mci8pGIfCAir4tIbqRjChCRx0Rktx3fX0QkI9IxhRKRT4vIThHxi0hUdAkTkQUiUiIi+0Tke5GOJ5SI/LeInBSRHZGOpT0RKRCRtSJSbP9O74l0TKFEJFFENonIh3Z8P4x0TD0xoBK5iBQA1wKHIx1LBx4zxkw2xkwFVgI/iHA8od4AJhpjJgN7gPsjHE97O4AbgbcjHQiAiDiBXwILgfHA50RkfGSjauNZYEGkg+iEF/i2MeZC4GLgG1H22TUDVxljpgBTgQUicnFkQ+regErkwBPAd4Cou4JrjKkJeZlMFMVojHndGOO1X24A8iMZT3vGmGJjTEmk4wgxC9hnjPnYGNMCvAB8KsIxBRlj3gYqIx1HR4wxx4wxW+2fa4FiIGoGejeWOvuly/4XNd/VzgyYRC4iS4AyY8yHkY6lMyLysIgcAT5PdLXIQ30ZWB3pIKJcHnAk5HUpUZSMYoWIjACmARsjHEobIuIUkQ+Ak8Abxpioiq8jMfWEIBF5ExjawawHgP8NzOvfiNrqKj5jzF+NMQ8AD4jI/cDdwEPREpu9zANYp77P91dcAT2JL4pIB9OivtUWTUQkBXgZ+Ld2Z6sRZ4zxAVPta0V/EZGJxpiou94QKqYSuTHmmo6mi8gkYCTwoYiAVRrYKiKzjDHHIx1fB/4IrKIfE3l3sYnI7cBi4GoTgZsLzuKziwalQEHI63zgaIRiiTki4sJK4s8bY5ZHOp7OGGM8IrIO63pDVCfyAVFaMcZsN8YMMcaMMMaMwPqiTe/PJN4dERkT8nIJsDtSsbQnIguA7wJLjDENkY4nBrwPjBGRkSISD9wCvBLhmGKCWC2t3wDFxpifRTqe9kQkO9BrS0TcwDVE0Xe1MwMikceIR0Vkh4h8hFUCiqZuV78AUoE37O6RT0c6oFAicoOIlAKXAKtEZE0k47EvDN8NrMG6WPdnY8zOSMYUSkT+BLwHFIlIqYj8S6RjCnEpcBtwlX2sfSAi10U6qBDDgLX29/R9rBr5ygjH1C29RV8ppWKctsiVUirGaSJXSqkYp4lcKaVinCZypZSKcZrIlVIqxmkiV0qpGKeJXCmlYtz/BxmvIXZanvYLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzaElEQVR4nO3deXxU5b348c93liSTjexAFggIhB2CuLQKaFGQgpTaXmtbbdVftbZarbdyldrb2l/1V+/1tlZvb+ulm12s1rZKtaioRXCpiGwKCsgiSxKWkJ1kkszy/P44M+NkT8gkM5N8368XLzJn/Z4zJ9885znPeR4xxqCUUip+2aIdgFJKqf7RRK6UUnFOE7lSSsU5TeRKKRXnNJErpVSc00SulFJxThP5GRCRYhExIuI4w/W/KCIvRjqugSAij4rIvdGOoyvtvwsReV5EvtybZWOBiBwSkUtiII64uCZF5CIRKetmfkxfrwNlWCdyEblQRP4pInUiUi0ib4jIORHeR4fkYYx5zBizKJL7Ceyr24u8F+tfKyKvRzKmdtvfICLNInI67N+zkdyHMWaJMea3kdxmrBORe0TkD/3ZxkBdk2pwxEzJZLCJSDrwd+BrwJNAAjAPaIlmXMPALcaYX0Y7CPUREXEYY7zRjgNiK5Z4MpxL5JMAjDGPG2N8xhi3MeZFY8y7ACJiE5HviMhhETkpIr8TkRGdbaj97XG7EtKrgf9rAyXQj7Uv+YrIx0Xk7cCdwdsi8vGweRtE5AeBu4UGEXlRRHI6iSEFeB7IDyvt5otIooj8REQqAv9+IiKJnaw/BXgE+Fhg3dqw2Zkisjaw/7dE5Kyw9SaLyEuBO5q9InJlTye+i3PY4W4gcCczIfCzS0R+FPg+6kTkdRFxdbKdDSLylcDPdhH5LxE5JSIHgaXtlh0hIr8SkWMiUi4i94qIPTDvLBFZLyJVgfUfE5GMsHUPicgdIvJuIJ4/iUhSN8d3g4jsDpzD90VkTifLtKkWaH+HJSJ3BuJsCJzrhSJyGfBt4HOB7+2dXhzbtYHr6UERqQbu6eSaNCJyk4jsE5EaEfkfEZGw8/qjwHn5UERukW6qrALnalXguGtE5DfBcxU8xsCxHQd+05trVkS+Hdj/IRH5YjfnfZmI7BCRWrHuvme2i2tl4DtsDJyvkWJVzzWIyMsiktnVtmPJcE7kHwA+EfmtiCzp5Au7NvDvYmA8kAr89Az2Mz/wf4YxJtUY82b4TBHJAtYCDwPZwI+BtSKSHbbYF4DrgDysO4c72u/EGNMILAEqAvtJNcZUAHcD5wOzgVnAucB3Oll/N3AT8GZg3Yyw2Z8Hvg9kAvuB+wKxpwAvAX8MxPZ54GciMq1XZ6Zv/gs4G/g4kAX8G+DvYZ0bgGVAKTAX+Gy7+b8FvMCEwDKLgK8E5gnwQyAfmAIUAfe0W/9K4DJgHDAT63rpQET+JbDul4B0YDlQ1UPs7bdRAtwCnGOMSQMWA4eMMS8A/w/4U+B7m9WLYwM4DziI9b3d18VulwHnYF03Vwb2CdZ5XYJ1Tc0BVvTiEL4YWP8srEJU+DU4Cus7HQvcSM/X7CggBygAvgysDpyfNgJ/LH8NfBXrd+t/gWfa/VH4DHBpIKbLsQpD3w5s3wbc2otji7phm8iNMfXAhYABfgFUisgzIjIysMgXgR8bYw4aY04Dq4Cruip19MNSYJ8x5vfGGK8x5nFgD9ZFFfQbY8wHxhg3VjXQ7D5s/4vA/zXGnDTGVGIl5Gv6GONTxpjNgVvex8L2vwwrmfwmEPs24K90TJjhHg6UjoL/ftDTzkXEBlwP3GaMKQ/cQf3TGNNTNdiVwE+MMUeNMdVYiTm4zZFYyeibxphGY8xJ4EHgKgBjzH5jzEvGmJbAefsxsKD9sRhjKgLbfpauv5evAP9pjHnbWPYbYw73dNzt+IBEYKqIOI0xh4wxBzpbsKdjC6gwxvx34Htzd7HP+40xtcaYI8ArYcd3JfCQMabMGFMD3N+L+H8a9j3ch/VHP8gPfC9wrt307pr998DyG7EKQp3dCd4A/K8x5q3ANfNbrKrT88OW+W9jzAljTDnwGvCWMWZ74Np6GuuPYMwbtnXkECqFXgtWFQHwB+AnWBdZPhD+y3YY63yNJLLa7ye4r4Kwz8fDfm7Cujs40+0fDkzri672PxY4T9pWwziA33ezrVvPoI48B0gCOk1c3cgHjoZ9Dj8PYwEncCxQYwBWweYogIjkYd0lzQPSAvNq2m2//Xnp6rwWnUHsbRhj9ovIN7FK9tNEZB3wr4G7rva6PbaAo+1X6kRX33v789qbbbX/HsLPVaUxpjnsc0/XbE3gDrSr+UFjgS+LyDfCpiW0W/ZE2M/uTj735XctaoZtibw9Y8we4FFgemBSBdaFEDQG61b1BB01Aslhn0eFb7qHXbffT3Bf5T2s15nO9tXZcXT2y9/V+t05Cmw0xmSE/Us1xnytj9uBdudQRMLP4SmgGeu2vC+OYSXRoDFhPx/FKp3lhMWebowJVgv9EOt8zDTGpANXY1W3nImj9C727q4jjDF/NMZciPV9GuA/grM62V93x9bZOn1xDCgM+1zU1YJdLNP+GmwfS0/XbGagWq+r+UFHgfvaXZ/JgbveIWXYJnKxHtJ9S0QKA5+LsErimwKLPA7cLiLjRCSVj+ohO3uivgOr2sUpIu3rYiuxbh3HdxHKc8AkEfmCiDhE5HPAVKwWNX11AsiWtg9lHwe+IyK5Yj0k/S7WnUdX6xeKSEIv9/f3QOzXBI7dKSLniPXgtK/ewSppzg48CLsnOMMY48eq6/yxWA9w7WI9NO7w0LadJ4FbRaQw8AzkrrBtHgNeBH4kIuliPdw+S0SC1SdpwGmsh9QFwMozOKagXwJ3iMjZYpkgIu3/eIN1HX1SRLICf8i+GZwhIiUi8onAMTdjlRZ9gdkngOJAFVRvjq2/ngRuE5ECsR4A39mLdW4OfA9ZWHXQf+pm2d5cs98XkQQRmYdVxffnTrbzC+AmETkvcN5TRGSpiKT1It64MmwTOdCA9cDnLRFpxErgu4BvBeb/GquK4FXgQ6xfnm90sh2Af8cqcdVg1ef9MTjDGNOEVSf4RqBOOLx+DmNMFdaF+C2sB2D/Biwzxpzq6wEF7ioeBw4G9pUP3AtsAd4FdgLbAtM6sx54DzguIj3u3xjTgPUQ7SqsEtFxrFJidwn2p9K2HfnWwLY+AP4v8DKwD2jfnv2OQPxvA9WB/fR0/f4CWIf1R2Ib8FS7+V/CutV+H+u7+wswOjDv+1gP8uqw6mDbr9trxpg/Y10Df8S67tZgPdxr7/eBWA9hJeLwZJeIVRd9Cus852ElRPgoiVWJyLZeHFt//SIQ37vAdqzCiJeP/rB05o+BdQ4G/nX30k5P1+xxrGOqwHpmc1Pg2m/DGLMFq578p4Hl99PFA+l4J0YHllBK9YOILAEeMcZ0dpeBiBwCvmKMeXlQAxtGhnOJXCl1BsRq0//JQFVgAfA9rBYeKko0kSul+kqwqp5qsKpWdmPVY6so0aoVpZSKc1oiV0qpOBeVF4JycnJMcXFxNHatlFJxa+vWraeMMbntp0clkRcXF7Nly5Zo7FoppeKWiHTatYNWrSilVJzTRK6UUnFOE7lSSsW5Yd37oVIq8jweD2VlZTQ3N/e8sOpUUlIShYWFOJ3OXi2viVwpFVFlZWWkpaVRXFxMWDe6qpeMMVRVVVFWVsa4ceN6tY4mchX31mwv54F1e6modZOf4WLl4hJWlBb0vKIaEM3NzZrE+0FEyM7OprKystfraCJXcW3N9nJWPbUTt8fqeK+81s2qp3YCaDKPIk3i/dPX86cPO1Vce2DdXtweHz5qMHgAcHt8PLBub5QjU2rwaCJXcavV18qhut002l6nyf4WHjkWmldR29UwlEr1rLa2lp/97Gd9Xm/Lli3cemvn4zUXFxdz6lSfhxnoFa1aUTGru7rvHcd3sO3YNpJSD+BpTCbJPx2H+WjchPwMV7TCVkNAMJF//etf7zDP5/Nht9s7XW/u3LnMnTt3oMPrQEvkKiYF677La90YoKy2gX996jn+vPVDAFKcKUzImsD3Fl5Htv1CnKYQwfrlcjntrFxcEsXoVbT97ne/Y+bMmcyaNYtrrrmGyspKPvOZz3DOOedwzjnn8MYbbwBwzz33cP3113PRRRcxfvx4Hn74YQDuuusuDhw4wOzZs1m5ciUbNmzg4osv5gtf+AIzZsygubmZ6667jhkzZlBaWsorr7wCwIYNG1i2bBkAVVVVLFq0iNLSUr761a8S7Gm2sbGRpUuXMmvWLKZPn86f/tTdqHe9oyVyFZOCdd9eavDYyvDJcYzfxw/XJfAvZ49jYvZEJmZPZP5YSE/M0lYrMezZvc92mHZW1llMzZ2K1+/l+X3Pd5hfklPCpOxJNHubeenAS23mXV5yebf7e++997jvvvt44403yMnJobq6mltuuYXbb7+dCy+8kCNHjrB48WJ2794NwJ49e3jllVdoaGigpKSEr33ta9x///3s2rWLHTt2AFaC3rx5M7t27WLcuHH86Ec/AmDnzp3s2bOHRYsW8cEHH7SJ4/vf/z4XXngh3/3ud1m7di2rV68G4IUXXiA/P5+1a9cCUFdX14uz2D1N5ComldXWU2d/Cq+tHIe/iCRTgtNfRE19ZodlV5QWaOJWIevXr+ezn/0sOTk5AGRlZfHyyy/z/vvvh5apr6+noaEBgKVLl5KYmEhiYiJ5eXmcOHGi0+2ee+65oXbdr7/+Ot/4hjWE7+TJkxk7dmyHRP7qq6/y1FNPhfaRmWlduzNmzOCOO+7gzjvvZNmyZcybN6/fx6yJXMWMysZKdp/azfyx8ynMSKf+tI1E3wyS/WdjD4xVrHXf8ae7ErTD5uh2fpIjqccSeHvGmA7N9/x+P2+++SYuV8frJzHxo7HC7XY7Xq+30+2mpKS02UdvdNaMcNKkSWzdupXnnnuOVatWsWjRIr773f4NsKR15Cqq/MbPloot3P/6/dzx4h38evuvaWhpYOXiEkbJVaT5Lw0lca37Vr2xcOFCnnzySaqqqgCorq5m0aJF/PSnPw0tE6wy6UpaWlqoxN6Z+fPn89hjjwHwwQcfcOTIEUpKSrpc5vnnn6empgaAiooKkpOTufrqq7njjjvYtm1bn4+xPS2Rq6h55/g7/Gr7r6h2V5PiTGHe2HlcMv4S0hLTWFGaBqB136rPpk2bxt13382CBQuw2+2Ulpby8MMPc/PNNzNz5ky8Xi/z58/nkUce6XIb2dnZXHDBBUyfPp0lS5awdOnSNvO//vWvc9NNNzFjxgwcDgePPvpom5I9wPe+9z0+//nPM2fOHBYsWMCYMWMAq1595cqV2Gw2nE4nP//5z/t9zFEZs3Pu3LlGB5YYfvzGz7vH3yXJmcTknMlU1Ffwsy0/46Lii5g3Zh6JjsSeN6Ji3u7du5kyZUq0w4h7nZ1HEdlqjOnQvlFL5GrA1TXX8Y8P/8Frh1+jsqmSmSNnMjlnMvnp+dz7iXujHZ5ScU8TuRpQj+54lI2HNuLxeyhKL+Kamdcwf+z8aIel1JCiiVxF1OnW02w8tJHFExbjsDlIT0xnbv5cFp21iInZE6MdnlJDkiZyFRF7Kvfw0sGX2H58O62+VrJd2ZxfdD5XTLki2qEpNeRpIlf9Uu2u5v7X76eioQKn3UnpqFIuHX8pU3L1YZdSg0UTueqzfVX7KG8o56Lii8hyZZHlyuLCMRdycfHFpCWmRTs8pYYdTeSqja56HGz2NLPx8EY2HNrA0fqjpCakMm/MPOw2O3ddeFe0w1aqS/fccw+pqanccccdnc5fs2YNkyZNYurUqYMcWeRoIlchwR4HmzytCPbQaDs7K9/gcPOLNHubyUnO4dNTPs3CcQux2zrvylOpeLJmzRqWLVumiVwNDf/vhbep9u3Eb2vFRhIO/0jwjOapLU18YUEJl551KbNGztJhvFTMu++++/jd735HUVERubm5nH322fziF79g9erVtLa2MmHCBH7/+9+zY8cOnnnmGTZu3Mi9997LX//6V9avX99hueTk5GgfUrc0kQ8hZzoIcZOniU1lm9h3+mk89nL8NCEkksQMnP7R1NbnsvKCawf+ANSQNNjd2G7dupUnnniC7du34/V6mTNnDmeffTZXXHEFN9xwAwDf+c53+NWvfsU3vvENli9fzrJly/jsZz8LQEZGRqfLxTJN5ENEt4MQz86H8FK0MSAS6iXOLnZWb1mNSTyI8aSQ4J9Ekn8WicZq9609Dqp48tprr/HpT386VIpevnw5ALt27eI73/kOtbW1nD59msWLF3e6fm+XiyURSeQikgH8EpgOGOB6Y8ybkdi26p3gQAwABg8+6mjyZPPAn95ixW9fhQcftJK5Mbhv/wYvZ5zi1Y8X8Z+X/ieJjkQuL7mc/KRant2cjcf3UcsT7XFQ9ddgd2MLnXcfe+2117JmzRpmzZrFo48+yoYNGzpdt7fLxZJIdWP7EPCCMWYyMAvYHaHtql4KH2zYKydw27fQaN/Ah/Ih1asfgttvp7HlNKvvvIQvHP8f7vNu4J0T77D3lDXa/BdnfpEfX34zD1wxj4IMFwIUZLj44RUztMdBFVfmz5/P008/jdvtpqGhgWeftap2GhoaGD16NB6PJ9S9LHTssrar5WJZv0vkIpIOzAeuBTDGtAKt/d2u6pv8DBflgWTuMKNJ8jvwSjnJqRU8ectF7HnnIXZ89SHcDhh51jiu+9SdrJiygpGpI9tsR0fbUfFuzpw5fO5zn2P27NmMHTs2NALPD37wA8477zzGjh3LjBkzQsn7qquu4oYbbuDhhx/mL3/5S5fLxbJ+d2MrIrOB1cD7WKXxrcBtxpjGrtbRbmwjr30dOQSqRZbk8eqJh9n88m8pOQk37IBFe1px2J3RC1YNadqNbWT0pRvbSFStOIA5wM+NMaVAI9DhDRERuVFEtojIlsrKygjsVoVbUVrAD6+Y0aZa5MaFNt489XP2b32R5Xvgz0/BJ/eD41srrQeeSqkhIRKJvAwoM8a8Ffj8F6zE3oYxZrUxZq4xZm5ubm4EdqvCtW96eMGMfWyv+QPH3/oHSzYe466ZN+Hw+uG22+Ahq85ck7lSQ0O/68iNMcdF5KiIlBhj9gILsapZ1CDprOnhHzYdZPYEO4udk/jkOfPI+NHPrFYrDz5orZSR0bZJolIR1NkAyKr3+lrlHal25N8AHhORBOAgcF2Etqt6Idj00CNl+E0riYzH6ZlHWdlhSq+bwcQx8z5K2sFkrr9kaoAkJSVRVVVFdna2JvMzYIyhqqqKpKSkXq8TkURujNkBdKiAV4OjotaNjzqabe/jlVMkeIsRbDQ0FHc+Go/+cqkBVFhYSFlZGfos7MwlJSVRWFjY6+X1zc4hICO9kiON27CRSJr3kzTb3sXpL2JsRu8vBKUixel0Mm7cuGiHMaxE6oUgFQV+4+fVw6/ysSknSLbnkOz7OD6pxCvHSXS26huZSg0TmsjjmE1stPpaufbci3l4xfXkjmjFYzvA6JRx/OiKJfpij1LDhFatxKGKhgpSE1JJT0xn4biFiAhNeU2ctreS6JjHFVOuwGHTr1ap4UJL5HHm3RPvsvaDtWwu3wx81DnQnlN78Pg9XDL+Ek3iSg0z+hsfJzw+DxsPb+RgzUHGZYxjwdgFbebPGT2H4oxislxZUYpQKRUtmsjjwOnW0zy/73lqm2s5r+A8Zo2aFZp3svEkLoeLtMQ0TeJKDVOayONAkiOJZGcyHy/6OAXpHz3AdHvcvHjgRVITUlkxeUX0AlRKRZXWkccoYwy7Tu7C4/PgsDlYOmlpmyRujOGVQ6/Q4m1h3ph5UYxUKRVtWiKPQS3eFtZ/uJ6j9Uexi50puR27BN1xfAdl9WXMGzOP7OTsKESplIoVmshjTFVTFS8eeJFGTyPzxszrNImfOH2CLRVbOCvzrE7nK6WGF03kMeRI3RFeOvASSY4klpcsJy8lr9PlslxZzBw5k9LRpYMcoVIqFmkijyFZrizGjBjDhWMuxOXsOHK9MQaf8eG0Ozmv8LwoRKiUikX6sDPKmjxNbKnYgjGG1IRULj3r0k6TOFgvAz29+2mavc2DHKVSKpZpiXyQhY/kkzPCzQXTKphZlMr4zPHdtgM/cfoEm8s3My5zHEmO3vdTrJQa+jSRD6LwkXxa5TAHT++hYksyZ4+8utsk3uxt5h8f/oPUhNTO+xdXSg1rWrUyiIIj+XjlJC223ThMDo7W83nklZPdrrfx0EaaPE1cMv4SEuwJgxStUipeaIl8EFXUugFwmDyS/LNxmJEIEpremRZvC42eRs4vPJ/cFB20WinVkSbyQZSf4aI8kLSdZlSb6V1JdCSyYvIKbKI3T0qpzml2GEQrF5fgctrbTHM57Z2O5NPibeH1I6/T4m3RJK6U6pZmiEG0orSAH14xg4IMFwIUZLj44RUzOh3JZ+Phjew5tYf6lvrBD1QpFVe0amWQrSgt6HEItl0nd3Go9pDWiyulekVL5DGmsrGSTWWbGDtiLDNHzox2OEqpOKAl8hgQ/pJQUto2lszM5JqZF0U7LKVUnNASeZQFXxIqr3VjgMaGqfx9cx7P7zwV7dCUUnFCE3mUBV8SCrKRRKsnhQfW7Y1iVEqpeKKJPMq6ehmou5eElFIqnCbyKOvqZaDuXhJSSqlwmsijrC8vCSmlVGci1mpFROzAFqDcGLMsUtsd6oJtyoOtVvIzXKxcXNJjW3OllAqKZPPD24DdQHoEtzks9OYlIaWU6kpEqlZEpBBYCvwyEttTSinVe5Eqkf8E+DcgrasFRORG4EaAMWPGRGi3SinVN+Ev4A2Vqsx+l8hFZBlw0hiztbvljDGrjTFzjTFzc3O1/xCl1OBr/wJeea2bVU/tZM328miH1i+RqFq5AFguIoeAJ4BPiMgfIrBdpZSKqOALeD5q8FEHgNvji/sX8PqdyI0xq4wxhcaYYuAqYL0x5up+R6aUUhEWfNGuxfYBLbY9HabHK21HrpQaNoIv2rn8c0jyz+gwPV5FNJEbYzZoG3KlVKwKvoAnOLGRDAyNF/C0RK6UGjZWlBawamkBI9KPYGjtdpSueKL9kSulhpXZ4ww3L7bzpVlLSHIkRTuciNASuVJqWKlx1+ByuIZMEgdN5EqpYaamuYZMV2a0w4goTeRKqWHDGEO1u5osV1a0Q4koTeRKqWHD7XUjyJBL5PqwUykVdYPV/0myM5lrZ1+LwUR829GkiVwpFVXB/k+CY9cG+z8BBiSZiwiCRHy70aRVK0qpqGo/ADkMXP8n245t462ytyK+3WjTRK6UiqpgPyd+mmiyvYWfpjbTI+lQ7SGq3dUR3260aSJXSkVVsJ8Tgwe/nKbJ/hY+GiLe/4kxhhr30Gt6CJrIlVJRFuz/xM4IXL7zAPAmvM0NCyLbsqS+pR6f8Q25FiugDzuVUlHWdgBymJB6ERfMOII/aRvHT49iVOqoiOynprkGgMykoVci10SulIq69gOQuz1uNpdvjmjp2W/8ZCRlkJGUEbFtxgpN5EqpmONyulhQvAAAr9/LkbojjM8c369tjs8c3+9txCqtI1dKxbSdJ3by8sGX2VrR7bDAw5omcqVUTJs1ahaTsiex9dhW3jjyBsb0/a1Mv/Hz+M7H2Xsqvsfm7IpWrSilYppNbCwYu4AkRxLvnniXFl8LFxVfhE16Xw6ta66jobUBu80+gJFGjyZypVTMExHOLzyfJEcSO0/spMnTRGpCaq/XH8otVkATuVIqjsweNZspOVNIdCQC4PF5cNqdPa5X7a5GkCHZYgW0jlwpFWeCSfzt8rf5296/0eRp6nGdGncN6YnpQ7ZqRRO5UiouFaQX0NDSwN/2/I36lvpul81NyWVC1oRBimzwaSJXSsWl/LR8lk1ahsfv4W97/sapplNdLjt71GzOzj97EKMbXJrIlVIDas32ci64fz3j7lrLBfevZ8328ohtOzcll+Uly7Hb7Dy37zk8Pk+HZXx+H37jj9g+e2sgj7s9fdiplBowgzFoREZSBstLllPVVNXpg8+DNQfZeHgjV067kvTE9IjssyeDPViGJnKl1IBpP2hEi3xAk+8Udz+/k9G580lLTCMtIY2iEUX92k9qQmqoOeKB6gP4jI9J2ZPAmFDTw9SEVDAGJLKjA/n8Pho9jZxuPc3p1tM0tjby/ee34/aMbLNccLAMTeRKqbjSfnAIIRHBSWXjKXad3IXP+Eh2JnP1zKsB2HBoAzXumlCCT0tMIyMpg/y0/F7vc1/1Po7UHaHlvx9kRm0i1TcvJiMpAxsCt98OGRlwzz293p7b46a+pb5NsraJjfMLzwfg2Q+e5WTjybbHfbocFyM7bGsgBssATeRKqQGUn+GiPCx5JZixYMZSkOHi+tKLcXvdNHubQ/NHJI6gydNEVVMVh2sP4zM+8lLyWDF5BQDP7XuOFm9Lm0Sf5cpq09XtpeMv5d6XnmRFqx+3awypv17H1RedD6tvh4cegttua1Myr2+pp9pdTWOrlagbPY20+lq5bMJlALxx9A0O1hwMbd9pc5KdnB36PHPkTLx+L6kJqaQ4U0hJSOG3L73a5rjDz8dA0ESulBowKxeXtKkrBnA57axcXIKIkOxMJtmZHJpXOrqUUkoBa0Qft9fd5gFmliuLand1m0Q/PnN8KJH/+b0/s/VwHY9tOkkTTrz2zdRSyc/X20l++zXmfuszNN5wCaffe5LPTP0MDpuD906+x86TVv21TWyhahq/8WMTGzPyZjApe1JoeoI9oc0xdtajYnfHPRD6nchFpAj4HTAK8AOrjTEP9Xe7Sqn413bQCDf5GS5WLi7pVT1xMNET9vwyWJ0BHyX6YIsUYwyjUkfxzPb9NPvq8IsHjxzCSxU1CdU8/LFi7rt6OWnGS3ZyNl6/F4fNwbS8aUzImkBqQipJjiSkXR36yNSOVSQDedxnQs6kJ7E2GxAZDYw2xmwTkTRgK7DCGPN+V+vMnTvXbNmypV/7VUqpzoy7ay0G8OPFbduElyo8tkPYzQh+2drE3FX/zbis+OyXXES2GmPmtp/e73bkxphjxphtgZ8bgN3AwPzZUUqpHgTroW04SPFfyAj/p0jxX0yGcXD09ed46f4beO3Qq1GOMrIiWkcuIsVAKfBWJ/NuBG4EGDNmTCR3q5RSIZ3VT2faZ3HlBYXs8zXzqn02DY8coShjPbcuHMP8ySmMThsdxYj7L2KJXERSgb8C3zTGdOj4wBizGlgNVtVKpParlFLhuqqf9vm93Ja4H+MZiZNMymvd3PnMWj5Z1syy6TM4J/8cclNyoxz9mel3HTmAiDiBvwPrjDE/7ml5rSNXSg22C+5fT1ltE8JHDzMNPjLTT/CvSxNp9jZTnFHM3Py5ER30OZIGrI5crEe8vwJ29yaJK6VUNFTUutskcQDBTl19PldNv4q5+XMpry9nx/Ed0QmwHyJRtXIBcA2wU0R2BKZ92xjzXAS2rZRSEdH+5aTw6Qn2BOaMnsPU3Kmh5oxVTVW8V/kec0bP6dNoRNHQ70RujHkdiGznBUopFWG9eUknyZEU+rmyqZIPqj5gX9U+puZOZfao2bicA/NmZn/pm51KqWGhry/pTM6ZTEFaAduObWPXyV3sPrWb0lGllI4uHcywe0UTuVJq2FhRWtCntyvTEtNYULyAWaNmsbViK16/NzQv+GZoLIiNKJRSKoZlJGWwcPzC0OcjdUfYcGgDpaNKmZo7NepjgeoIQUop1UfJzmSyXdm8WfYmT+x6gt2Vu6MyClGQJnKllOqjnOQclk5ayrJJy0hJSOG1I6+x9oO1UYtHq1aUUuoM5afls2LyCo7UHcHnt1rD+I2fo3VHGZsxdtDi0ESulFL9NGbER/1H7a/ez4ZDG8hNzuWcgnMoTC8c8P1rIldKqQiakDUBYwxbj23luX3PMTp1NOcUnMOmfb4B659cE7lSSkWQTWyU5JQwIWsCe07tYfvx7Ty48Vme2VQUehmpvNbNqqesUYkikcz1YadSSg0Au83OtLxpXDX9Kja8O6rNG6UAbo9VQo8ETeRKKTWAHDYHJ+s6T7UVnfT9ciY0kSul1AALjlrU2+l9pYlcKaUG2MrFJbicbd/+bN9hV3/ow06llBpgfe2wq680kSul1CDoa4ddfaFVK0opFec0kSulVJzTRK6UUnFOE7lSSsU5TeRKKRXnNJErpVSc00SulFJxThO5UkrFOU3kSikV5zSRK6VUnNNErpRScU4TuVJKxTlN5EopFec0kSulVJyLSCIXkctEZK+I7BeRuyKxTaWUUr3T70QuInbgf4AlwFTg8yIytb/bVUop1TuRKJGfC+w3xhw0xrQCTwCfisB2lVJK9UIkEnkBcDTsc1lgWhsicqOIbBGRLZWVlRHYrVJKKYhMIpdOppkOE4xZbYyZa4yZm5ubG4HdKqWUgsgk8jKgKOxzIVARge0qpZTqhUgk8reBiSIyTkQSgKuAZyKwXaWUUr3g6O8GjDFeEbkFWAfYgV8bY97rd2RKKaV6pd+JHMAY8xzwXCS2pZRSqm/0zU6llIpzmsiVUirOaSJXSqk4p4lcKaXinCZypZSKc5rIlVIqzmkiV0qpOKeJXCml4pwmcqWUinOayJVSKs5pIldKqTiniVwppeKcJnKllIpzmsiVUirOaSJXSqk4p4lcKaXinCZypZSKc5rIlVIqzmkiV0qpOBeRMTsHw5rt5Tywbi8VtW7yM1ysXFzCitKCaIellFJRFxeJfM32clY9tRO3x4efZspq/ax6aieAJnOl1LAXF4n8gXV7cXt8ALhtW7CRiHjO4YF1e5lUWEeyM5nMpEzSE9MRkShHq5RSgysuEnlFrTv0c4I5C8EOQHltI/88+k/8xg+AXexkujKZkjOFKblTAGhoaSA1IVUTvFJqyIqLRJ6f4aI8kMydZnRoekFGCl+e9WVqmmuocddQ01xDtbs6NL+xtZHHdz2O0+Yk05VJZlImma5Mxo4Yy4ikEYN+HEopNRDiIpGvXFwSqiMPcjntrFxcgtPuJC8lj7yUvA7rOWwO5o2ZF0rwR+qOsLdqL8nOZEYkjeBU0yn+efSfZCZlkuXKItNl/Z/kSBrMw4s7+uBZqdgSF4k8mCT6mjwSHYmhKpagZm8zdrGqZnx+6w/DgZoD7D61O7TM5ZMuZ3TaaKqaqjjZeDKU4BPsCZE8rLgU/uAZoLzWrQ+elYoyMcYM+k7nzp1rtmzZMuj77U6Tp4lqdzU17homZU8i0ZHI9mPbebvi7dAyKc4UMl2ZLBy3kERHYuiPgtPujGLkg+uC+9dTXuvGTwteKcdpChESKMhw8cZdn4h2eEoNaSKy1Rgzt/30uCiRD4ZkZzLJzmQK0wtD00pHlzIxe2Iowdc011DbXBsqmb9d/ja7T+0mLSEtVGrPcmUxIWtCtA5jwAUfPBuaabF9gPidOE1RmwfSSqnBpYm8B6kJqaQmpDJmxJgO8yZmTyQlISWU5Mvqy0hxpoQS+etHXqfJ02TVvwfq4UckjcAm8ftCbfDBs50R2HDhkeM4TRH5Ga5oh6bUsNWvRC4iDwCXA63AAeA6Y0xtBOKKC6NSRzEqdVTos9/4afI0hT4LQm1zLYdrD2OwqrBGp47m8pLLAdhduRuX0xVXbeDDHzw7/Pm02g6Q6PSycnFJtENTatjqb4n8JWCVMcYrIv8BrALu7H9Y8ckmNlITUkOfLxhzAWA9VK1trqWmuQaHzTrlxhjeLHsTr98LWG3gM5IymJwzmWl504DYbAMf/uD5aO0okpOPcPOFafqgU6ko6lciN8a8GPZxE/DZ/oUzNNltdrKTs8lOzg5NExG+NOtLbdq/17hrQvPdHjeP73och83RpnnkmBFjyEjKiMJRfGRFaUEoca//MI/ijOKoxqPUcBfJOvLrgT91NVNEbgRuBBgzpmN983DksDnITcklNyW3wzy7zc78sfNDCT7YBj7RnkhGUgbV7mpeO/xam/bvmUmZuJx9q6vub5vwT4yLbksVbdOuVC8SuYi8DIzqZNbdxpi/BZa5G/ACj3W1HWPMamA1WM0PzyjaYSTBnsDknMltpjV7m0MPSr1+LyLCwZqDtJxqCS2zdOJSCtILqHZXc/z08VCCT3QkdthHpNqEt/paafI0DfqdgrZpV8rSYyI3xlzS3XwR+TKwDFhootEofRgJf+M0LyWP5SXLAasNfI3bqp4JVt+U1ZexqWxTaPlkZzJZriwuLr4Yl9NFs7eZ/3jhvbDOyN5BcOD3jOeBdXv7lAif2/ccACsmr+jvIfaKz++jobWBe194lQZPKw5yQvPcHl+f41cq3vW31cplWA83FxhjmnpaXg2MYBv4gvSPktfMkTMZnzk+lOCDbeCDJfNtx7ax7/SziD0JMal4pRwjLST4z6Ki1s2+qn0cO32M0amjGZ02us1D3PbGZYzjrbJNNLQ0kJaYZk00Bs7wIa0xhiZPEw2tDXj93lDb/vUfrqeioSLUMujQ6Z3YbFk4/Dlt1tc27Wq46W8d+U+BROClQMuKTcaYm/odlYqIYBv4ohFFHeadlXkWBSnTOX66Cj8NCEmIScNGEvkZLt44+ga7Tu4iwZ5AijOFvJQ8JmZP5KLiizpsa/wjT/BW02YOFJzL7NGlVhK//XbIyIB77uk0tlZfK/Ut9TR7m0OJenP5Zj6s+ZCG1oZQj5bpielcNf0qgNALW+mJ6aQlpDEuNYsTdR3b5GubdjXc9LfVytB9hXGIG5k6ku8vWRaqYzb4MbSEOiMbme3C5XBR3lBOXXMd71W+x5H6I6FE/uiOR3HanEzIOouJdTXkvbCeg+bfmf2fz8Ltt+N/+CFO33oT9XVHKQz8IXm/8n32nNpDQ0sDLT6rXt8udq4vvR4RwWFzkJ2cTXFGMWmJaaQnppOemB6K+fzC89scw92XJXfZmZpSw4m+2TmMdeyMLCus1UcBHyv6GD6/j7qWOqrd1aHOxgD2nNrDycaTPL//eWzzbbQkTaTu0FrevdjGmDrw3bEQc+Vs2P8818y8JtSaJsmRRF5KHmkJVqIOVcUAc0bP6Wf82mpFDU/aaZY6I16/N/TW6qayTaz/8B+cfP5p7EBaC0y5biVn55/NxKyJzBw5E4ddywxK9Zd2mqUiKvii0j7fPpw2B5/f0sLFT8K+LNhYDGWvrGP9x2v5sPZDZo+eDcDx08dJdia3qS5RSvWfJnJ1Rpo8Tbyw/wVONVYy/fF/cN6Dz2G/9TayH3yQ827/Jh/+78Ns9Rbg+srlobbvD7zyV9a8e4D6JgejUgq5dcH5XPexOaFuC5RSZ0Z/g1SPOnt78vJZo0iwJ7BowmKK047BrbfBgw+CCPLgTxiPMD4lA+/EJQD85s3t/OHt3fi8yQiJHGs8wHfX7edY4yG+u+hKAGqba6Pe/YBS8UjryFW3wt+eNHhplf2McEzi/ivmtH2o2L7deLvP5/3wGQ417MQrx3H55mInEx815I9IY/OqK6hvqeeJXU+Q4kyhaEQRRelFFKQX6KhMSoXROnJ1Rh5Ytxe3x4ef07jtW/Hj5rQ3s+Pbk+1f/mn3+WSdHRez8dOEjWQAHORQWWfNT7QnMn/sfI7WHeVgzUH2nNqDICyZuITC9EJ8fh82scVUT5BKxQpN5KpbwbckhQTEJJDsn4mdzD6/PRkckCKYxMOngzW+6uScyUzOmYzf+DnZeJKy+jJykq23Nt+rfI8dx3dQlF5E0YgiCtMLdZBspQI0katuBROwkECy/2NtpvdF+IAUQV29vGMTW4dBO7Jd2RSmF3K0/ij7qvcBVn8znyr5lJbS1bCniVx1qy8JuDv9fXmnIL2AgvQCjDGcajpFWX0ZTZ6mUBJ/8cCL2MRGYXohRelFpCSk9Ck+peKZJnLVrUi+PRk+IMWZEpFO+3BPdiZzuPYwB2sOApDlymJ63vQOXQErNRRpqxU1pFS7qymrL+No3VHGZoxlet50WrwtvHLolVD9ur6QpOKVtlpRw0KWK4ssVxYzR84MTWtobaC2uZYjdUfgqNWjYlF6ETNGztCkroYETeRqyMtJzuGq6VdR31LP0bqjHK0/yt6qvUzPmw5AeX05Ve4qitKLyHRlRjlapfpOE7kaNtIT05mWN41pedPw+X3YbVZvjmX1Zbxz4h02sanNC0njMseF1tWxQVUs00SuhqVgEgc4r/A8puVNC9WtH6w5yPHTx0OJ/H9ff5P/WneEVk8KgujYoCrmaCJXCms0pfAXkhpbGwFr2Ln/2riGan8TdlsWyf5zAR0bVMUWTeRKtWMTW2jACxGhteE8kuQUYG+znI4NqmKFJnKlelCYkUF5bWKH6To2qIoVHUeuVUq1sXJxCS5n29K4jg2qYomWyJXqgY4NqmKdJnKleiES3QsoNVC0akUppeKcJnKllIpzmsiVUirOaSJXSqk4p4lcKaXiXFT6IxeRSuBwhDebA5yK8DYjJZZjg9iOL5Zjg9iOL5Zjg9iOL1ZjG2uMyW0/MSqJfCCIyJbOOlyPBbEcG8R2fLEcG8R2fLEcG8R2fLEcW2e0akUppeKcJnKllIpzQymRr452AN2I5dggtuOL5dggtuOL5dggtuOL5dg6GDJ15EopNVwNpRK5UkoNS5rIlVIqzg25RC4id4iIEZGcaMcSTkR+ICLvisgOEXlRRPKjHVOQiDwgInsC8T0tIhnRjimciPyLiLwnIn4RiYkmYSJymYjsFZH9InJXtOMJJyK/FpGTIrIr2rG0JyJFIvKKiOwOfKe3RTumcCKSJCKbReSdQHzfj3ZMvTGkErmIFAGXAkeiHUsnHjDGzDTGzAb+Dnw3yvGEewmYboyZCXwArIpyPO3tAq4AXo12IAAiYgf+B1gCTAU+LyJToxtVG48Cl0U7iC54gW8ZY6YA5wM3x9i5awE+YYyZBcwGLhOR86MbUs+GVCIHHgT+DYi5J7jGmPqwjynEUIzGmBeNMd7Ax01AYTTjac8Ys9sYszfacYQ5F9hvjDlojGkFngA+FeWYQowxrwLV0Y6jM8aYY8aYbYGfG4DdQMx09G4spwMfnYF/MfO72pUhk8hFZDlQbox5J9qxdEVE7hORo8AXia0SebjrgeejHUSMKwCOhn0uI4aSUbwQkWKgFHgryqG0ISJ2EdkBnAReMsbEVHydiasRgkTkZWBUJ7PuBr4NLBrciNrqLj5jzN+MMXcDd4vIKuAW4HuxEltgmbuxbn0fG6y4gnoTXwyRTqbFfKktlohIKvBX4Jvt7lajzhjjA2YHnhU9LSLTjTEx97whXFwlcmPMJZ1NF5EZwDjgHREBq2pgm4ica4w5Hu34OvFHYC2DmMh7ik1EvgwsAxaaKLxc0IdzFwvKgKKwz4VARZRiiTsi4sRK4o8ZY56KdjxdMcbUisgGrOcNMZ3Ih0TVijFmpzEmzxhTbIwpxvpFmzOYSbwnIjIx7ONyYE+0YmlPRC4D7gSWG2Oaoh1PHHgbmCgi40QkAbgKeCbKMcUFsUpavwJ2G2N+HO142hOR3GCrLRFxAZcQQ7+rXRkSiTxO3C8iu0TkXawqoFhqdvVTIA14KdA88pFoBxRORD4tImXAx4C1IrIumvEEHgzfAqzDelj3pDHmvWjGFE5EHgfeBEpEpExE/k+0YwpzAXAN8InAtbZDRD4Z7aDCjAZeCfyevo1VR/73KMfUI31FXyml4pyWyJVSKs5pIldKqTiniVwppeKcJnKllIpzmsiVUirOaSJXSqk4p4lcKaXi3P8HvA9amGFhGLAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# call the plotting functions\n",
    "plotTheResults(centroids_KMEANS, points, labels_KMEANS, kmeans=True)\n",
    "plotTheResults(centroids_EUCL, points, labels_EUCL, kmeans=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\"><img alt=\"Creative Commons License\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by/4.0/80x15.png\" /></a><br />This work is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\">Creative Commons Attribution 4.0 International License</a>. The **MOSEK** logo and name are trademarks of <a href=\"http://mosek.com\">Mosek ApS</a>. The code is provided as-is. Compatibility with future release of **MOSEK** or the `Fusion API` are not guaranteed. For more information contact our [support](mailto:support@mosek.com). "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  },
  "vscode": {
   "interpreter": {
    "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}