{ "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 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": [ "
" ] }, "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.sub( 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( Expr.sub(di, dijAux),Domain.greaterThan(0)), # di >= dijAux\n", " DJC.term( label , Domain.equalsTo(clusterInd) )) # Y_i = j\n", " \n", " M.disjunction( [ ANDs[i] for i in range(0, numberOfClusters) ] )\n", "\n", "\n", " # symmetry breaking constraints\n", " M.constraint( Expr.sub(centroids.slice([0, 0], [numberOfClusters-1, 1]), \n", " centroids.slice([1, 0], [numberOfClusters, 1])), Domain.lessThan(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": [ "
" ] }, "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": [ "
" ] }, "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.sub( 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( Expr.sub(di, dijAux),Domain.greaterThan(0)), # di >= dijAux\n", " DJC.term( label , Domain.equalsTo(clusterInd) )) # Y_i = j\n", "\n", " M.disjunction( [ ANDs[i] for i in range(0, numberOfClusters) ] )\n", "\n", " # symmetry breaking constraints\n", " M.constraint( Expr.sub(centroids.slice([0, 0], [numberOfClusters-1, 1]), \n", " centroids.slice([1, 0], [numberOfClusters, 1])), Domain.lessThan(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": [ "
" ] }, "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": [ "
" ] }, "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": [ "\"Creative
This work is licensed under a Creative Commons Attribution 4.0 International License. The **MOSEK** logo and name are trademarks of Mosek ApS. 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 }