Clustering_Kmean_class.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/OptimizationExpert/Pyomo/blob/main/Clustering_Kmean_class.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oPrjdRysqXo5"
      },
      "source": [
        "Developed by Alireza Soroudi\n",
        "alireza.soroudi@gmail.com\n",
        "\n",
        "Linkedin News letter : https://www.linkedin.com/newsletters/optimization-in-open-source-6874020019009859585/\n",
        "\n",
        "Git hub: https://github.com/OptimizationExpert/Pyomo\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "HtMwjJMo_kE6"
      },
      "outputs": [],
      "source": [
        "!pip install -q pyomo\n",
        "from pyomo.environ import *\n",
        "import copy\n",
        "\n",
        "\n",
        "from pyomo import environ as pe\n",
        "import os\n",
        "\n",
        "os.environ['NEOS_EMAIL'] = 'XXXXXX@gmail.com'\n",
        "opt = pe.SolverManagerFactory('neos')\n",
        "\n",
        "#!wget -N -q \"https://github.com/thomasfork/ipopt_linux/raw/main/ipopt-linux64.zip\"\n",
        "#!unzip -o -q ipopt-linux64\n",
        "#opt = SolverFactory('ipopt', executable='/content/ipopt')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pQC338vq_XDn"
      },
      "outputs": [],
      "source": [
        "import networkx as nx\n",
        "import matplotlib.pyplot as plt # Data visualization\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "import time\n",
        "import random"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Us3Rpb-VqNlc"
      },
      "source": [
        "$$\\forall_{i,k} \\ \\ \\ \\ \\  d_{i,k} \\geq \\sum_p (X_{i,p}-\\mu_{k,p})^2 - M(1-U_{i,k})$$\n",
        "\n",
        "$$\\forall_i \\ \\ \\ \\sum_k U_{i,k}=1$$\n",
        "\n",
        "$$\\min_{U_,\\mu} OF= \\sum_{i,k}\\forall_{i,k}$$\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "nT5VbbHRqEZr"
      },
      "outputs": [],
      "source": [
        "df= pd.DataFrame()\n",
        "N = 30\n",
        "Kcluster = 3\n",
        "df['x'] = [random.randint(0,3)+(0.5-random.random()) for i in range(N)]\n",
        "df['y'] = [random.randint(0,3)+(0.5-random.random()) for i in range(N)]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DtsOAEUfxQZD",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 540
        },
        "outputId": "359f0413-5f14-41f9-865c-6396a4937ce2"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.collections.PathCollection at 0x7bba8f679cf0>"
            ]
          },
          "metadata": {},
          "execution_count": 38
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 600x600 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhAAAAH6CAYAAABWN3XwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA6LklEQVR4nO3dfXBU9b3H8c/ZJezWMYlmhARKmkDxEkPkoQQhFJRWHqqOY2buqKVS1Eto6cAduThtSadilbmTPmgfppcr9uYqvVDHlrbEe6kiAcqTJJogTElA5iJJoJZN6I1mQ9qNYffcP2JWQx5P2LO7Z/f9mtk/9uR3ku9x2d2Pv/N7MEzTNAUAAGCBK9YFAAAA5yFAAAAAywgQAADAMgIEAACwjAABAAAsI0AAAADLCBAAAMAyAgQAALBsVKwLiLRQKKS//OUvSk1NlWEYsS4HAADHME1T7e3tGj9+vFyuwfsYEi5A/OUvf1F2dnasywAAwLEuXLigCRMmDNrG1gDx3HPP6bnnnlNjY6MkaerUqdq4caPuuuuufttv3bpVjz76aK9jHo9HgUBg2H8zNTVVUvfFp6WljaxwAACSkN/vV3Z2dvi7dDC2BogJEybo+9//vm6++WaZpqlf/vKXuu+++3T8+HFNnTq133PS0tJ05syZ8HOrtyF62qelpREgAAAYgeF899oaIO69995ez//1X/9Vzz33nKqrqwcMEIZhKCsry86yAADANYraLIxgMKiXX35ZHR0dKioqGrDd5cuXlZOTo+zsbN13332qr68f9Pd2dnbK7/f3egAAAHvZHiBOnjyp66+/Xh6PR6tXr9bOnTuVn5/fb9spU6bohRde0CuvvKLt27crFApp3rx5+vOf/zzg7y8rK1N6enr4wQBKAADsZ5imadr5Bz788EOdP39ebW1t+u1vf6vy8nIdPHhwwBDxSV1dXbrlllu0bNkybdq0qd82nZ2d6uzsDD/vGQDS1tbGGAgAACzw+/1KT08f1neo7dM4R48ercmTJ0uSZs2apZqaGv3sZz/T888/P+S5KSkpmjlzps6ePTtgG4/HI4/HE7F6AQDA0KK+EmUoFOrVYzCYYDCokydPaty4cTZXBQAArLC1B6K0tFR33XWXPvOZz6i9vV0vvfSSDhw4oNdff12StGLFCn36059WWVmZJOnpp5/W3LlzNXnyZH3wwQf60Y9+pKamJpWUlNhZJgAAsMjWANHS0qIVK1bo4sWLSk9P17Rp0/T6669r8eLFkqTz58/3Wirz/fff16pVq+Tz+XTjjTdq1qxZOnr06LDGSwAAgOixfRBltFkZAAIAAD5m5TuU3TgBAIBlBAgAAGAZAQIAAFhGgADgeIGuoC61dyrQFYx1KUDSsH0hKQCwS01jq8oPn1PlqWaFTMllSIvzM7VqwSQV5mbEujwgodEDAcCRtlU36YEtVdp7ukWhj+aShUxp7+kW3b+lSturm2JbIJDgCBAAHKemsVUbK+pkSgqGes9ED4ZMmZKeqKhTbWNrTOoDkgEBAoDjlB8+J5fLGLSNy2Wo/EhDlCoCkg8BAoCjBLqCqjzV3Kfn4WrBkKk99T4GVgI2IUAAcJT2wBUNkR3CQmZ3ewCRR4AA4Cip3lEa4u5FmMvobg8g8ggQABzFm+LW4vxMuYdIEW6XoSVTs+RNcUepMiC5ECAAOE7JgkkKDXEfIxQyVTJ/YpQqApIPAQKA48zOzdCm4gIZUp+eCLfLkCFpU3EBi0kBNuLmIABHWj43R3lZqSo/0qA99b5eK1GWzJ9IeABsRoAA4FiFuRkqzM1QoCuo9sAVpXpHMeYBiBICBADH86a4CQ5AlDEGAgAAWEaAAAAAlhEgAACAZQQIAABgGQECAABYRoAAAACWESAAAIBlBAgAAGAZAQIAAFhGgAAAAJYRIAAAgGUECAAAYBkBAgAAWEaAAAAAlhEgAACAZQQIAABgGQECAABYRoAAAACWESBsEugK6lJ7pwJdwViXAgBAxI2KdQGJpqaxVeWHz6nyVLNCpuQypMX5mVq1YJIKczNiXR4AABFBD0QEbatu0gNbqrT3dItCZvexkCntPd2i+7dUaXt1U2wLBAAgQggQEVLT2KqNFXUyJQV70sNHgiFTpqQnKupU29gak/oAAIgkAkSElB8+J5fLGLSNy2Wo/EhDlCoCAMA+BIgICHQFVXmquU/Pw9WCIVN76n0MrAQAOB4BIgLaA1c0RHYIC5nd7QEAcDICRASkekdpiLsXYS6juz0AAE5GgIgAb4pbi/Mz5R4iRbhdhpZMzZI3xR2lygAAsIetAeK5557TtGnTlJaWprS0NBUVFem1114b9JwdO3YoLy9PXq9Xt956q1599VU7S4yYkgWTFBriPkYoZKpk/sQoVQQAgH1sDRATJkzQ97//fR07dky1tbX64he/qPvuu0/19fX9tj969KiWLVumlStX6vjx4youLlZxcbHq6ursLDMiZudmaFNxgQypT0+E22XIkLSpuIDFpAAACcEwTXOYw/8iIyMjQz/60Y+0cuXKPj978MEH1dHRoV27doWPzZ07VzNmzNCWLVuG9fv9fr/S09PV1tamtLS0iNU9XLWNrSo/0qA99b7wSpRLpmapZP5EwgMAIK5Z+Q6N2mi+YDCoHTt2qKOjQ0VFRf22qaqq0vr163sdW7p0qSoqKgb8vZ2dners7Aw/9/v9Eal3pApzM1SYm6FAV1DtgStK9Y5izAMAIOHYHiBOnjypoqIiBQIBXX/99dq5c6fy8/P7bevz+ZSZmdnrWGZmpnw+34C/v6ysTE899VREa77aSMKAN8VNcACuASEciG+2B4gpU6boxIkTamtr029/+1s9/PDDOnjw4IAhwqrS0tJevRZ+v1/Z2dkR+d1sjAVEH+87wBlsDxCjR4/W5MmTJUmzZs1STU2Nfvazn+n555/v0zYrK0vNzc29jjU3NysrK2vA3+/xeOTxeCJbtLo3xtpYUSeXy+izMdae+mZtKi7Q8rk5Ef+7QDLjfQc4R9TXgQiFQr3GLHxSUVGR9u3b1+tYZWXlgGMm7MLGWED08b4DnMXWAFFaWqpDhw6psbFRJ0+eVGlpqQ4cOKCHHnpIkrRixQqVlpaG2z/22GPavXu3nn32Wb3zzjv63ve+p9raWq1du9bOMvtgYywg+njfAc5i6y2MlpYWrVixQhcvXlR6erqmTZum119/XYsXL5YknT9/Xi7Xxxlm3rx5eumll/Td735X3/nOd3TzzTeroqJCBQUFdpbZS8/GWEPtbfHJjbEY4AVcG953gPPYGiD+8z//c9CfHzhwoM+x+++/X/fff79NFQ1tJBtj8UEGXBved4DzsBfGVdgYC4g+3neA8xAgrsLGWED08b4DnIcA0Q82xgKij/cd4CwEiH6wMRYQfbzvAGfhRuIAls/NUV5Wap+NsRbnZ7IxFmAT3neAc0R9N0672bEbJ2vyA9HH+w6IvrjcjdPJ2BgLiD7ed0B8YwwEAACwjAABAAAsI0AAAADLCBAAAMAyAgQAALCMAAEAACwjQAAAAMsIEAAAwDICBAAAsIwAAQAALCNAAAAAywgQAADAMgIEAACwjAABAAAsI0AAAADLCBAAAMAyAgQAALCMAAF8JNAV1KX2TgW6grEuBQDi3qhYFwDEWk1jq8oPn1PlqWaFTMllSIvzM7VqwSQV5mbEujwAiEv0QCCpbatu0gNbqrT3dItCZvexkCntPd2i+7dUaXt1U2wLBIA4RYBA0qppbNXGijqZkoI96eEjwZApU9ITFXWqbWyNSX0AEM8IEEha5YfPyeUyBm3jchkqP9IQpYoAwDkIEEhKga6gKk819+l5uFowZGpPvY+BlQBwFQIEklJ74IqGyA5hIbO7PQDgYwQIJKVU7ygNcfcizGV0twcAfIwAgaTkTXFrcX6m3EOkCLfL0JKpWfKmuKNUGQA4AwECSatkwSSFhriPEQqZKpk/MUoVAYBzECCQtGbnZmhTcYEMqU9PhNtlyJC0qbiAxaQAoB/c2EVSWz43R3lZqSo/0qA99b5eK1GWzJ9IeACAARAgkPQKczNUmJuhQFdQ7YErSvWOYswDAAyBAAF8xJviJjgAwDAxBgIAAFhGgAAAAJYRIAAAgGUECAAAYBkBAgAAWEaAAAAAltkaIMrKyjR79mylpqZq7NixKi4u1pkzZwY9Z+vWrTIMo9fD6/XaWSYAALDI1gBx8OBBrVmzRtXV1aqsrFRXV5eWLFmijo6OQc9LS0vTxYsXw4+mpiY7ywQAABbZupDU7t27ez3funWrxo4dq2PHjun2228f8DzDMJSVlWVnaQAA4BpEdQxEW1ubJCkjY/D9BS5fvqycnBxlZ2frvvvuU319/YBtOzs75ff7ez0AAIC9ohYgQqGQ1q1bp89//vMqKCgYsN2UKVP0wgsv6JVXXtH27dsVCoU0b948/fnPf+63fVlZmdLT08OP7Oxsuy4BAAB8xDBN04zGH/rGN76h1157TUeOHNGECROGfV5XV5duueUWLVu2TJs2berz887OTnV2doaf+/1+ZWdnq62tTWlpaRGpHUB8YyM0IDL8fr/S09OH9R0alc201q5dq127dunQoUOWwoMkpaSkaObMmTp79my/P/d4PPJ4PJEoE4DD1DS2qvzwOVWeau61FfuqBZPYih2wma23MEzT1Nq1a7Vz507t379fEydOtPw7gsGgTp48qXHjxtlQIQCn2lbdpAe2VGnv6RaFPupHDZnS3tMtun9LlbZXM3sLsJOtAWLNmjXavn27XnrpJaWmpsrn88nn8+nvf/97uM2KFStUWloafv70009rz549OnfunN5++20tX75cTU1NKikpsbNUAA5S09iqjRV1MiUFQ73vwgZDpkxJT1TUqbaxNSb1AcnA1lsYzz33nCRp4cKFvY6/+OKLeuSRRyRJ58+fl8v1cY55//33tWrVKvl8Pt14442aNWuWjh49qvz8fDtLBeAg5YfPyeUy+oSHT3K5DJUfaeBWBmCTqA2ijBYrA0AAOE+gK6j8jbs1SHYIcxnSqae/xMBKYJisfIeyFwYAR2kPXBlWeJC6x0S0B67YWxCQpAgQABwl1TtKLmN4bV1Gd3sAkUeAAOAo3hS3Fudnyj1EinC7DC2ZmsXtC8AmBAgAjlOyYJJCQ9zHCIVMlcy3PnUcwPAQIAA4zuzcDG0qLpAh9emJcLsMGZI2FRcwAwOwETcHATjS8rk5ystKVfmRBu2p9/VaibJk/kTCA2AzAgQAxyrMzVBhbgZ7YQAxQIAA4HjeFDfBAYgyxkAAAADLCBAAAMAyAgQAALCMAAEADhDoCupSe6cCXcFYlwJIYhAlAMS1msZWlR8+p8pTzb2mqq5aMImpqogpeiAAIE5tq27SA1uqtPd0S3gDsZAp7T3dovu3VGl7dVNsC0RSI0AAQByqaWzVxoo6mZKCVy3bHQyZMiU9UVGn2sbWmNQHECAAIA6VHz4n1xAbhrlchsqPNESpIqA3AgQAxJlAV1CVp5r79DxcLRgytafex8BKxAQBAgDiTHvgiobIDmEhs7s9EG0ECACIM6neURri7kWYy+huD0QbAQIA4ow3xa3F+Zl9tiq/mttlaMnULPYBQUwQIAAgDpUsmKTQEPcxQiFTJfMnRqkioDcCBJIWK/shns3OzdCm4gIZUp+eCLfLkCFpU3EBi0khZrhxhqTDyn5wiuVzc5SXlaryIw3aU+/r9e+1ZP5E/r0ipgzTNIc51tcZ/H6/0tPT1dbWprS0tFiXgzizrbpJGyvq5HIZvabIuV2GQiFTm4oLtHxuTgwrBPoX6AqqPXBFqd5RjHmAbax8h3ILA0mDlf3gZN4Ut8akeggPiBsECCQNVvYDgMghQCApsLIfAEQWAQJJgZX9ACCyCBBICqzsBwCRRYBAUmBlPwCILAIEkgYr+wFA5BAgkDRY2Q8AIocbvUgqrOwHAJFBgEDSKczNUGFuBiv7AcA1IEAgaXlT3AQHABghxkAAAADLCBAAAMAyAgQAALCMAAEAACwjQAAAAMsIEAAAwDICBAAAsIwAAQAALCNAAAAAywgQAADAMlsDRFlZmWbPnq3U1FSNHTtWxcXFOnPmzJDn7dixQ3l5efJ6vbr11lv16quv2lkmAACwyNYAcfDgQa1Zs0bV1dWqrKxUV1eXlixZoo6OjgHPOXr0qJYtW6aVK1fq+PHjKi4uVnFxserq6uwsFQAAWGCYpmlG649dunRJY8eO1cGDB3X77bf32+bBBx9UR0eHdu3aFT42d+5czZgxQ1u2bOnTvrOzU52dneHnfr9f2dnZamtrU1paWuQvAgCABOX3+5Wenj6s79CojoFoa2uTJGVkZAzYpqqqSosWLep1bOnSpaqqquq3fVlZmdLT08OP7OzsyBUMAAD6FbUAEQqFtG7dOn3+859XQUHBgO18Pp8yMzN7HcvMzJTP5+u3fWlpqdra2sKPCxcuRLRuAADQ16ho/aE1a9aorq5OR44ciejv9Xg88ng8Ef2dAABgcFEJEGvXrtWuXbt06NAhTZgwYdC2WVlZam5u7nWsublZWVlZdpYIAAAssPUWhmmaWrt2rXbu3Kn9+/dr4sSJQ55TVFSkffv29TpWWVmpoqIiu8oEAGBIga6gLrV3KtAVjHUpccHWHog1a9bopZde0iuvvKLU1NTwOIb09HR96lOfkiStWLFCn/70p1VWViZJeuyxx3THHXfo2Wef1T333KOXX35ZtbW1+sUvfmFnqQAA9KumsVXlh8+p8lSzQqbkMqTF+ZlatWCSCnMHnhSQ6GydxmkYRr/HX3zxRT3yyCOSpIULFyo3N1dbt24N/3zHjh367ne/q8bGRt1888364Q9/qLvvvntYf9PKFBQAAAazrbpJGyvq5HIZCoY+/rp0uwyFQqY2FRdo+dycGFYYWVa+Q6O6DkQ0ECAAAJFQ09iqB7ZUabAvSUPSjtVFCdMTEbfrQAAA4BTlh8/J5eq/J72Hy2Wo/EhDlCqKLwQIAACuEugKqvJUc6/bFv0Jhkztqfcl5cBKAgQAAFdpD1zRENkhLGR2t082BAgAAK6S6h2lIe5ehLmM7vbJhgABAMBVvCluLc7PlHuIFOF2GVoyNUveFHeUKosfBAgAAPpRsmCSQkPcxwiFTJXMH3qRxEREgAAQN1jpD/Fkdm6GNhUXyJD69ES4XYYMSZuKCxJmCqdVyXfTBkDcYaU/xKvlc3OUl5Wq8iMN2lPv6/Xvs2T+xKT+98lCUgBiKtlW+oNzBbqCag9cUap3VMKOeWAhKQCOUNPYqo0VdTKlPvPtgyFTpqQnKupU29gak/qAT/KmuDUm1ZOw4cEqAgSAmGGlP8C5CBAAYoKV/gBnI0AAiAlW+gOcjQABICZY6Q9wNgIEgJhgpT/A2QgQAGKGlf4A5yJAAIgZVvoDnIubigBiipX+AGciQACIucLcDBXmZiTFSn9AoiBAAIgb3hQ3wQFwCMZAAAAAywgQAADAMgIEAACwjAABAAAsI0AADhboCupSeycbTQGIOmZhAA5U09iq8sPnVHmqude6CasWTGLdBABRQQ8E4DDbqpv0wJYq7T3dEt7NMmRKe0+36P4tVdpe3RTbAgEkBQIE4CA1ja3aWFEnU1Lwqj0kgiFTpqQnKupU29gak/oAJA8CBOAg5YfPyTXE7pUul6HyIw1RqghAsiJAAA4R6Aqq8lRzn56HqwVDpvbU+xhYCcBWBAjAIdoDVzREdggLmd3tI41ZHwB6MAsDcIhU7yi5DA0rRLiM7vaRwqwPAFejBwJwCG+KW4vzM+UeYgyE22VoydSsiG1KxawPAP0hQAAOUrJgkkJDdEGEQqZK5k+MyN9j1geAgRAgAAeZnZuhTcUFMqQ+PRFulyFD0qbigojdVmDWB4CBMAYCcJjlc3OUl5Wq8iMN2lPv6zUmoWT+xIiFh55ZH0ONufjkrI9I3TYBEP8IEIADFeZmqDA3Q4GuoNoDV5TqHRXxL++RzPogQADJgwABRJidX+pX86a4bfsbsZz1ASD+8Y4HIiTRpjr2zPrYe7pl0MWr3C5Di/Mz6X0AkgyDKIEISNSpjtGe9QHAOQgQwDVK5KmO0Z71AcA5CBDANUr0qY7L5+Zox+oiLc7PVM9l9tye2bG6SMvn5sS2QAAxwRgI4Boky1THaMz6AOAsBAjgGiTbVEc7Z30AcBZbb2EcOnRI9957r8aPHy/DMFRRUTFo+wMHDsgwjD4Pn89nZ5nAiPVMdRwOpjoCSCS2BoiOjg5Nnz5dmzdvtnTemTNndPHixfBj7NixNlUIXJtYbXAFALFm6/8O3XXXXbrrrrssnzd27FjdcMMNw2rb2dmpzs7O8HO/32/57wHXomTBJO2pbx60DVMdASSauJyFMWPGDI0bN06LFy/WG2+8MWjbsrIypaenhx/Z2dlRqhLoxlRHAMkorgLEuHHjtGXLFv3ud7/T7373O2VnZ2vhwoV6++23BzyntLRUbW1t4ceFCxeiWDHQjamOAJJNXI3omjJliqZMmRJ+Pm/ePL377rv6yU9+om3btvV7jsfjkcfjiVaJwICY6oh4xb9J2CGuAkR/brvtNh05ciTWZQDDxlRHxItE258F8SWubmH058SJExo3blysywAAR0nU/VkQP2ztgbh8+bLOnj0bft7Q0KATJ04oIyNDn/nMZ1RaWqr33ntP//Vf/yVJ+ulPf6qJEydq6tSpCgQCKi8v1/79+7Vnzx47ywSAhDLU/ixS9/4seVmp9ERgxGwNELW1tfrCF74Qfr5+/XpJ0sMPP6ytW7fq4sWLOn/+fPjnH374oR5//HG99957uu666zRt2jTt3bu31+8AAAyuZ3+WwbZh79mfhQCBkTJM0xzmQrzO4Pf7lZ6erra2NqWlpcW6HACIqkBXUPkbdw9riXWXIZ16+kuM2UGYle/QuB8DAQAYvpHsz5KMAl1BXWrvVKArGOtSHCvuZ2EAAIavZ3+W4fZAJNv+LIk0MyXW03OT618OACS4nv1Z9p5uGXQMhNtlaHF+ZlLdvthW3aSNFXVyuYw+M1P21DdrU3GBIxZ9i5cQxC0MAEgwJQsmKTREF0Sy7c8y1MwUU90zU2obW2NS33DF0/RcAgQAJBj2Z+mrZ2bKYHpmpsSreAtBBAgASEDsz/KxQFdQlaeaB72lI3V/Ce+p98XtwMp4C0GMgQCABMX+LN1GMjMl3v479YSgoa7jkyHI7msgQABAgkv2/VkSYWZKPIYgbmEAABJaz8yUq8eDXM3tMrRkalZchq2eEDQc0QpBBIgkxAIqAJKN02emxGMIir9+GtgmXuYOA0C09cxMeeKjdSA+OaDS7TIUCplxPzOlZMEk7alvHrRNNEMQPRBJIp7mDgNALDh9Zkq8Tc9lM60kUNPYqge2VGmwF9qQtGN1UVynbwCIFCfPTKltbFX5kQbtqfeFe5OXTM1SyfyJ1/wZbuU7lFsYSYCtfQGgNyfPTImX6bkEiAQXj3OHAQDXLtYhiDEQCY6tfQEAdiBAJLh4nDsMAHA+AkSCi8e5wwAA5yNAJAGnL6ACAIg/BIgkEG9zhwEAzscN7ySxfG6O8rJS+8wdXpyfGZG5wwCA5EKASCLxMncY9uM1BmA3AkQSivXcYdiH/U4ARAtjIIAEwX4nAKKJAAEkgJrGVm2sqJMp9VmyPBgyZUp6oqJOtY2tMakPQOIhQAAJoGe/k8H07HcCAJFAgAAcrme/k8E2S5N673eSSAJdQV1q70y46wLiHYMoAYcbyX4niTCIlgGjQGzRAwE4XDLud8KAUSD2CBCAwyXbficMGAXiAwECSADJtN8JA0aB+ECAABJAsux3kuwDRoF44vyboQAkJcd+J8k6YBSIRwQIIIEk+n4nPQNGhxMiEmXAKBCvuIUBJCBviltjUj0JFR6k5BswCsQzAgRijoWAYEUyDRgF4hn9e4gZFgLCSPQMGH2iok4ul9FrQKXbZSgUMhNiwCgQ7wzTNIc5JMkZ/H6/0tPT1dbWprS0tFiXgwFsq27SxiG+AJbPzYlhhYh3tY2tfQaMLpmalTADRoFYsPIdSg8Eom6ohYCk7oWA8rJS+SLAgBJ9wCgQ7xgDgahjISBEUqIOGAXiHQECUcVCQACQGAgQiKqRLAQEAIg/BAhEVTLuHAkAiYgAgahiISAASAy2BohDhw7p3nvv1fjx42UYhioqKoY858CBA/rc5z4nj8ejyZMna+vWrXaWiBhgISAAcD5bA0RHR4emT5+uzZs3D6t9Q0OD7rnnHn3hC1/QiRMntG7dOpWUlOj111+3s0xEWbLsHAkAiSxqC0kZhqGdO3equLh4wDbf/va39Yc//EF1dXXhY1/+8pf1wQcfaPfu3f2e09nZqc7OzvBzv9+v7OxsFpJyABYCAoD44tiFpKqqqrRo0aJex5YuXap169YNeE5ZWZmeeuopmyuDHVgICACcK64GUfp8PmVmZvY6lpmZKb/fr7///e/9nlNaWqq2trbw48KFC9EoFRHEQkAA4Dxx1QMxEh6PRx6PJ9ZlAACQVOKqByIrK0vNzc29jjU3NystLU2f+tSnYlQVAAC4WlwFiKKiIu3bt6/XscrKShUVFcWoIgAA0B9bA8Tly5d14sQJnThxQlL3NM0TJ07o/PnzkrrHL6xYsSLcfvXq1Tp37py+9a1v6Z133tG///u/6ze/+Y3+5V/+xc4yAQCARbYGiNraWs2cOVMzZ86UJK1fv14zZ87Uxo0bJUkXL14MhwlJmjhxov7whz+osrJS06dP17PPPqvy8nItXbrUzjIBAIBFUVsHIlqszGEFAAAfs/IdGldjIAAAgDMQIAAAgGUECAAAYBkBAgAAWEaAAAAAlhEgAACAZQQIAABgGQECAABYRoAAAACWESAAAIBlBAgAAGAZAQIAAFhGgACAGAp0BXWpvVOBrmCsSwEsGRXrAgAgGdU0tqr88DlVnmpWyJRchrQ4P1OrFkxSYW5GrMsDhkQPBABE2bbqJj2wpUp7T7coZHYfC5nS3tMtun9LlbZXN8W2QGAYCBAAEEU1ja3aWFEnU1KwJz18JBgyZUp6oqJOtY2tMakPGC4CBABEUfnhc3K5jEHbuFyGyo80RKkiYGQIEAAQJYGuoCpPNffpebhaMGRqT72PgZUJLBEGzzKIEgCipD1wRUNkh7CQ2d3em+K2tyhEVSINnqUHAgCiJNU7SkPcvQhzGd3tkTgSbfAsAQIAosSb4tbi/Ey5h0gRbpehJVOz6H1IIIk4eJYAAQBRVLJgkkJD3McIhUyVzJ8YpYoQDYk4eJYAAQBRNDs3Q5uKC2RIfXoi3C5DhqRNxQWOux+OgSXq4FlusAFAlC2fm6O8rFSVH2nQnnpfr8F0JfMnEh4STKIOniVAAEAMFOZmqDA3Q4GuoNoDV5TqHeWILw1Y1zN4djghwkmDZ7mFAQAx5E1xa0yqh/CQwBJ18CwBAgAAmyXi4FkCBAAANkvEwbPOuNECAIDDJdrgWQIEgIhiUCAwsEQaPEuAABARibTGP2A3b4rbscGhB2MgAFyzRFvjH9cuEXabxODogQBwTYZa41/qXuM/LyuVnogkQE9U8qAHAsA1ScQ1/jEy9EQlFwIEgBFL1DX+YV0i7jaJwREgAIzYSNb4R2KiJyr5ECAAjFjPGv/D4aQ1/mENPVHJiQABYMQSdY1/WENPVHIiQAC4Jom4xj+soScqOREgAFyTRFzjH9bQE5WcCBAArtnyuTnasbpIi/Mzw/8n2jP/f8fqIi2fmxPbAmE7eqKSD/1IACIikdb4h3U9PVFPVNTJ5TJ6Dah0uwyFQiY9UQmGAAEgohJhjX+MTKLtNonBESAAABFDT1TyiMoYiM2bNys3N1der1dz5szRW2+9NWDbrVu3yjCMXg+v1xuNMgEAEeJNcWtMqofwkMBsDxC//vWvtX79ej355JN6++23NX36dC1dulQtLS0DnpOWlqaLFy+GH01NrJ8OAEA8sT1A/PjHP9aqVav06KOPKj8/X1u2bNF1112nF154YcBzDMNQVlZW+JGZmWl3mQAAwAJbA8SHH36oY8eOadGiRR//QZdLixYtUlVV1YDnXb58WTk5OcrOztZ9992n+vr6Adt2dnbK7/f3eiDxBbqCutTeyZK4ABAjtg6i/Otf/6pgMNinByEzM1PvvPNOv+dMmTJFL7zwgqZNm6a2tjY988wzmjdvnurr6zVhwoQ+7cvKyvTUU0/ZUj/iT01jq8oPn1PlqeZeI7xXLZjECG8AiKK4W0iqqKhIK1as0IwZM3THHXfo97//vcaMGaPnn3++3/alpaVqa2sLPy5cuBDlihEt26qb9MCWKu093RJedz9kSntPt+j+LVXaXs1YGQCIFlt7IG666Sa53W41Nzf3Ot7c3KysrKxh/Y6UlBTNnDlTZ8+e7ffnHo9HHo/nmmtFfKtpbNXGijqZUp8d/3qeP1FRp7ysVHoiACAKbO2BGD16tGbNmqV9+/aFj4VCIe3bt09FRUXD+h3BYFAnT57UuHHj7CoTDlB++JxcQ6yz73IZKj/SEKWKACC52b6Q1Pr16/Xwww+rsLBQt912m37605+qo6NDjz76qCRpxYoV+vSnP62ysjJJ0tNPP625c+dq8uTJ+uCDD/SjH/1ITU1NKikpsbtUxKlAVzA85mEwwZCpPfU+BbqCzD0HAJvZHiAefPBBXbp0SRs3bpTP59OMGTO0e/fu8MDK8+fPy+X6uCPk/fff16pVq+Tz+XTjjTdq1qxZOnr0qPLz8+0uFXGqPXBlyPDQI2R2tydAAIC9DNM0h/nR7Ax+v1/p6elqa2tTWlparMtBBAS6gsrfuHtYIcJlSKee/hIBAgBGwMp3aNzNwgCu5k1xa3F+ptxDjIFwuwwtmZpFeACAKCBAwBFKFkxSaIguiFDIVMn8iVGqCACSGwECjjA7N0ObigtkSH16ItwuQ4akTcUFTOEEgChhO284xvK5OcrLSlX5kQbtqff1WomyZP5EwgMARBEBAo5SmJuhwtwMBbqCag9cUap3FGMeACAGCBBwJG+Km+AAADHEGAgAAGAZAQIAAFhGgAAAAJYRIAAAgGUECAAAYBkBAgAAWEaAAAAAlhEgAACAZQQIAABgGQECAABYRoAAAACWESAAAIBlBAgAAGAZAQIAAFhGgAAAAJYRIAAAgGUECAAAYBkBAgAAWEaAAAAAlhEgAACAZQQIAABgGQECAABYRoAAAACWESAAAIBlBAgAAGAZAQIAAFhGgAAAAJYRIAAAgGUECAAAYBkBAgAAWEaAAAAAlhEgAACAZQQIAABgGQECAABYRoAAAACWESAAAIBlBAgAAGAZAQIAAFgWlQCxefNm5ebmyuv1as6cOXrrrbcGbb9jxw7l5eXJ6/Xq1ltv1auvvhqNMgHbBLqCutTeqUBXMNalAEBE2B4gfv3rX2v9+vV68skn9fbbb2v69OlaunSpWlpa+m1/9OhRLVu2TCtXrtTx48dVXFys4uJi1dXV2V0qEHE1ja36+rZa5W/crdn/ulf5G3fr69tqVdvYGuvSAOCaGKZpmnb+gTlz5mj27Nn6t3/7N0lSKBRSdna2/vmf/1kbNmzo0/7BBx9UR0eHdu3aFT42d+5czZgxQ1u2bBny7/n9fqWnp6utrU1paWmRuxDAom3VTdpYUSeXy1Aw9PHbzO0yFAqZ2lRcoOVzc2JYIQD0ZuU71NYeiA8//FDHjh3TokWLPv6DLpcWLVqkqqqqfs+pqqrq1V6Sli5dOmD7zs5O+f3+Xg8g1moaW7Wxok6m1Cs86KPnpqQnKuroiQDgWLYGiL/+9a8KBoPKzMzsdTwzM1M+n6/fc3w+n6X2ZWVlSk9PDz+ys7MjUzxwDcoPn5PLZQzaxuUyVH6kIUoVAUBkOX4WRmlpqdra2sKPCxcuxLokJLlAV1CVp5r79DxcLRgytafex8BKAI40ys5fftNNN8ntdqu5ubnX8ebmZmVlZfV7TlZWlqX2Ho9HHo8nMgUDEdAeuKIhskNYyOxu701x21sUAESYrT0Qo0eP1qxZs7Rv377wsVAopH379qmoqKjfc4qKinq1l6TKysoB2wPxJtU7SkPcvQhzGd3tAcBpbL+FsX79ev3Hf/yHfvnLX+r06dP6xje+oY6ODj366KOSpBUrVqi0tDTc/rHHHtPu3bv17LPP6p133tH3vvc91dbWau3atXaXCkSEN8WtxfmZcg+RItwuQ0umZtH7AMCRbP9fnwcffFCXLl3Sxo0b5fP5NGPGDO3evTs8UPL8+fNyuT7OMfPmzdNLL72k7373u/rOd76jm2++WRUVFSooKLC7VCBiShZM0p765kHbhEKmSuZPjFJFABBZtq8DEW2sA4F4sb26SU+wDgQAB7HyHcrNV8Amy+fmKC8rVeVHGrSn3qeQ2T3mYXF+pkrmT1RhbkasSwSAESNAADYqzM1QYW6GAl1BtQeuKNU7ijEPABICAQKIAm+Km+AAIKE4fiEpAAAQfQQIAABgGQECAABYRoAAAACWESAAAIBlBAgAAGAZAQIAEBGBrqAutXeyRX2SYB0IAMA1qWlsVfnhc6o81dxrxdVVCyax4moCowcCADBi26qb9MCWKu093aKeLV9CprT3dIvu31Kl7dVNsS0QtiFAAABGpKaxVRsr6mRKvTaM00fPTUlPVNSptrE1JvXBXgQIAMCIlB8+J5fLGLSNy2Wo/EhDlCpCNBEgAACWBbqCqjzV3Kfn4WrBkKk99T4GViYgAgQAwLL2wBUNkR3CQmZ3eyQWAgQAwLJU7ygNcfcizGV0t0diIUAAACzzpri1OD9T7iFShNtlaMnULLazT0AECADAiJQsmKTQEPcxQiFTJfMnRqkiRBMBAgAwIrNzM7SpuECG1Kcnwu0yZEjaVFzAYlIJiptSAIARWz43R3lZqSo/0qA99b5eK1GWzJ9IeEhgBAgAwDUpzM1QYW6GAl1BtQeuKNU7ijEPSYAAAQCICG+Km+CQRBgDAQAALCNAAAAAywgQAADAMgIEAACwjAABAAAsI0AAAADLCBAAAMAyAgQAALCMAAEAACwjQAAAAMsIEAAAwDICBAYU6ArqUnunAl3BWJcCAIgzbKaFPmoaW1V++JwqTzX32pp31YJJbM0LAJBEDwSusq26SQ9sqdLe0y0Kmd3HQqa093SL7t9Spe3VTbEtEAAQFwgQMRRvtwhqGlu1saJOpqRgT3r4SDBkypT0REWdahtbY1IfAMSjePssjxZuYcRAvN4iKD98Ti6X0Sc8fJLLZaj8SAO3MgAkvXj9LI8WeiCiLF5vEQS6gqo81TxoeJC6eyL21PuSLmkDwCfF62d5NBEgoiiebxG0B65oiOwQFjK72wNAMornz/JoIkBEUc8tgsH03CKItlTvKA1RWpjL6G4PAMkonj/Lo4kAESXxfovAm+LW4vxMuYd4U7hdhpZMzZI3xR2lygAgfsT7Z3k0ESCixAm3CEoWTFJoiCJDIVMl8ydGqSIAiC9O+CyPFtsCRGtrqx566CGlpaXphhtu0MqVK3X58uVBz1m4cKEMw+j1WL16tV0lRpUTbhHMzs3QpuICGVKfngi3y5AhaVNxQVKMLgaA/jjhszxabAsQDz30kOrr61VZWaldu3bp0KFD+trXvjbkeatWrdLFixfDjx/+8Id2lRhVTrlFsHxujnasLtLi/Mzwm6RnatKO1UVaPjcnJnUBQDxwymd5NNgSjU6fPq3du3erpqZGhYWFkqSf//znuvvuu/XMM89o/PjxA5573XXXKSsry46yYq5kwSTtqW8etE083CIozM1QYW6GAl1BtQeuKNU7KqHfBABghVM+y+1mSw9EVVWVbrjhhnB4kKRFixbJ5XLpzTffHPTcX/3qV7rppptUUFCg0tJS/e1vfxu0fWdnp/x+f69HvHLaLQJviltjUj2EBwD4BKd9ltvFlh4In8+nsWPH9v5Do0YpIyNDPp9vwPO+8pWvKCcnR+PHj9ef/vQnffvb39aZM2f0+9//fsBzysrK9NRTT0Wsdrstn5ujvKxUlR9p0J56X6/Vy0rmT0z4f3AAkAj4LLcYIDZs2KAf/OAHg7Y5ffr0iIv55BiJW2+9VePGjdOdd96pd999V5/97Gf7Pae0tFTr168PP/f7/crOzh5xDdHALQIAcL5k/yy3FCAef/xxPfLII4O2mTRpkrKystTS0tLr+JUrV9Ta2mppfMOcOXMkSWfPnh0wQHg8Hnk8nmH/znjiTXEn1T82AEhEyfpZbilAjBkzRmPGjBmyXVFRkT744AMdO3ZMs2bNkiTt379foVAoHAqG48SJE5KkcePGWSkTAADYzJZBlLfccou+9KUvadWqVXrrrbf0xhtvaO3atfryl78cnoHx3nvvKS8vT2+99ZYk6d1339WmTZt07NgxNTY26r//+7+1YsUK3X777Zo2bZodZQIAgBGybR2IX/3qV8rLy9Odd96pu+++W/Pnz9cvfvGL8M+7urp05syZ8CyL0aNHa+/evVqyZIny8vL0+OOP6x//8R/1P//zP3aVCAAARsgwTXOYi3I6g9/vV3p6utra2pSWlhbrcgAAcAwr36HshQEAACwjQAAAAMsIEAAAwDICBAAAsIwAAQAALCNAAAAAywgQAADAMgIEAACwzJbtvGOpZ10sv98f40oAAHCWnu/O4awxmXABor29XZLifktvAADiVXt7u9LT0wdtk3BLWYdCIf3lL39RamqqDMPo83O/36/s7GxduHAh4Za6TtRrS9Trkrg2p0rUa0vU65K4tuEyTVPt7e0aP368XK7BRzkkXA+Ey+XShAkThmyXlpaWcP+IeiTqtSXqdUlcm1Ml6rUl6nVJXNtwDNXz0INBlAAAwDICBAAAsCzpAoTH49GTTz4pj8cT61IiLlGvLVGvS+LanCpRry1Rr0vi2uyQcIMoAQCA/ZKuBwIAAFw7AgQAALCMAAEAACwjQAAAAMsIEAAAwLKkCBCtra166KGHlJaWphtuuEErV67U5cuXBz1n4cKFMgyj12P16tVRqnhgmzdvVm5urrxer+bMmaO33npr0PY7duxQXl6evF6vbr31Vr366qtRqtQaK9e1devWPq+N1+uNYrXDd+jQId17770aP368DMNQRUXFkOccOHBAn/vc5+TxeDR58mRt3brV9jpHwuq1HThwoM/rZhiGfD5fdAoeprKyMs2ePVupqakaO3asiouLdebMmSHPc8J7bSTX5oT323PPPadp06aFV2IsKirSa6+9Nug5Tni9JOvXFs3XKykCxEMPPaT6+npVVlZq165dOnTokL72ta8Ned6qVat08eLF8OOHP/xhFKod2K9//WutX79eTz75pN5++21Nnz5dS5cuVUtLS7/tjx49qmXLlmnlypU6fvy4iouLVVxcrLq6uihXPjir1yV1L9n6ydemqakpihUPX0dHh6ZPn67NmzcPq31DQ4PuuecefeELX9CJEye0bt06lZSU6PXXX7e5UuusXluPM2fO9Hrtxo4da1OFI3Pw4EGtWbNG1dXVqqysVFdXl5YsWaKOjo4Bz3HKe20k1ybF//ttwoQJ+v73v69jx46ptrZWX/ziF3Xfffepvr6+3/ZOeb0k69cmRfH1MhPcqVOnTElmTU1N+Nhrr71mGoZhvvfeewOed8cdd5iPPfZYFCocvttuu81cs2ZN+HkwGDTHjx9vlpWV9dv+gQceMO+5555ex+bMmWN+/etft7VOq6xe14svvmimp6dHqbrIkWTu3Llz0Dbf+ta3zKlTp/Y69uCDD5pLly61sbJrN5xr++Mf/2hKMt9///2o1BQpLS0tpiTz4MGDA7ZxynvtasO5Nqe+32688UazvLy835859fXqMdi1RfP1SvgeiKqqKt1www0qLCwMH1u0aJFcLpfefPPNQc/91a9+pZtuukkFBQUqLS3V3/72N7vLHdCHH36oY8eOadGiReFjLpdLixYtUlVVVb/nVFVV9WovSUuXLh2wfSyM5Lok6fLly8rJyVF2dvaQadxJnPCaXasZM2Zo3LhxWrx4sd54441YlzOktrY2SVJGRsaAbZz6ug3n2iRnvd+CwaBefvlldXR0qKioqN82Tn29hnNtUvRer4TbjfNqPp+vTxfpqFGjlJGRMei916985SvKycnR+PHj9ac//Unf/va3debMGf3+97+3u+R+/fWvf1UwGFRmZmav45mZmXrnnXf6Pcfn8/XbPp7uOY/kuqZMmaIXXnhB06ZNU1tbm5555hnNmzdP9fX1w9qJNZ4N9Jr5/X79/e9/16c+9akYVXbtxo0bpy1btqiwsFCdnZ0qLy/XwoUL9eabb+pzn/tcrMvrVygU0rp16/T5z39eBQUFA7ZzwnvtasO9Nqe8306ePKmioiIFAgFdf/312rlzp/Lz8/tt67TXy8q1RfP1cmyA2LBhg37wgx8M2ub06dMj/v2fHCNx6623aty4cbrzzjv17rvv6rOf/eyIfy+uXVFRUa/0PW/ePN1yyy16/vnntWnTphhWhsFMmTJFU6ZMCT+fN2+e3n33Xf3kJz/Rtm3bYljZwNasWaO6ujodOXIk1qVE3HCvzSnvtylTpujEiRNqa2vTb3/7Wz388MM6ePDggF+0TmLl2qL5ejk2QDz++ON65JFHBm0zadIkZWVl9RmMd+XKFbW2tiorK2vYf2/OnDmSpLNnz8YkQNx0001yu91qbm7udby5uXnA68jKyrLUPhZGcl1XS0lJ0cyZM3X27Fk7SoyqgV6ztLQ0R/c+DOS2226L2y/ntWvXhgddD/V/bk54r32SlWu7Wry+30aPHq3JkydLkmbNmqWamhr97Gc/0/PPP9+nrdNeLyvXdjU7Xy/HjoEYM2aM8vLyBn2MHj1aRUVF+uCDD3Ts2LHwufv371coFAqHguE4ceKEpO5u2FgYPXq0Zs2apX379oWPhUIh7du3b8B7YUVFRb3aS1JlZeWg986ibSTXdbVgMKiTJ0/G7LWJJCe8ZpF04sSJuHvdTNPU2rVrtXPnTu3fv18TJ04c8hynvG4jubarOeX9FgqF1NnZ2e/PnPJ6DWSwa7uara9XVIZqxtiXvvQlc+bMmeabb75pHjlyxLz55pvNZcuWhX/+5z//2ZwyZYr55ptvmqZpmmfPnjWffvpps7a21mxoaDBfeeUVc9KkSebtt98eq0swTdM0X375ZdPj8Zhbt241T506ZX7ta18zb7jhBtPn85mmaZpf/epXzQ0bNoTbv/HGG+aoUaPMZ555xjx9+rT55JNPmikpKebJkydjdQn9snpdTz31lPn666+b7777rnns2DHzy1/+sun1es36+vpYXcKA2tvbzePHj5vHjx83JZk//vGPzePHj5tNTU2maZrmhg0bzK9+9avh9ufOnTOvu+4685vf/KZ5+vRpc/Pmzabb7TZ3794dq0sYkNVr+8lPfmJWVFSY//u//2uePHnSfOyxx0yXy2Xu3bs3VpfQr2984xtmenq6eeDAAfPixYvhx9/+9rdwG6e+10ZybU54v23YsME8ePCg2dDQYP7pT38yN2zYYBqGYe7Zs8c0Tee+XqZp/dqi+XolRYD4v//7P3PZsmXm9ddfb6alpZmPPvqo2d7eHv55Q0ODKcn84x//aJqmaZ4/f968/fbbzYyMDNPj8ZiTJ082v/nNb5ptbW0xuoKP/fznPzc/85nPmKNHjzZvu+02s7q6OvyzO+64w3z44Yd7tf/Nb35j/sM//IM5evRoc+rUqeYf/vCHKFc8PFaua926deG2mZmZ5t13322+/fbbMah6aD1TF69+9FzPww8/bN5xxx19zpkxY4Y5evRoc9KkSeaLL74Y9bqHw+q1/eAHPzA/+9nPml6v18zIyDAXLlxo7t+/PzbFD6K/a5LU63Vw6nttJNfmhPfbP/3TP5k5OTnm6NGjzTFjxph33nln+AvWNJ37epmm9WuL5utlmKZpRr5fAwAAJDLHjoEAAACxQ4AAAACWESAAAIBlBAgAAGAZAQIAAFhGgAAAAJYRIAAAgGUECAAAYBkBAgAAWEaAAAAAlhEgAACAZf8PkUZ5b65XoWwAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ],
      "source": [
        "plt.figure(figsize=(6,6))\n",
        "plt.scatter(df[\"x\"],df[\"y\"], s=  50)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "n-p84G1k8Uep"
      },
      "source": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vWCtZjqL8Ubk"
      },
      "source": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "MiBZog3dx-qQ"
      },
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "L7gUFyvdqEeX",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 990
        },
        "outputId": "27f0ff3d-5e97-4085-a62c-6fde5f9b4c15"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "           x         y\n",
              "0   0.187026 -0.422885\n",
              "1   0.420190  1.057513\n",
              "2   3.482497  1.365564\n",
              "3   3.205376 -0.456993\n",
              "4   1.199737  1.185130\n",
              "5   2.399816  1.426568\n",
              "6  -0.260194  3.057673\n",
              "7  -0.377020  2.942529\n",
              "8   0.939371 -0.023406\n",
              "9   1.614755  1.816510\n",
              "10  3.030900  1.372515\n",
              "11  0.735238  2.662303\n",
              "12  3.149366  1.908935\n",
              "13  1.644236  2.942909\n",
              "14  2.177196  0.766498\n",
              "15  0.539502  0.451881\n",
              "16  2.265196  2.286545\n",
              "17  1.483374  0.616614\n",
              "18  1.400418  1.602659\n",
              "19 -0.179868 -0.452139\n",
              "20  0.902008  1.540868\n",
              "21  0.697423  2.192745\n",
              "22  2.664715  0.595353\n",
              "23  1.746641  3.348775\n",
              "24  1.733469  1.142788\n",
              "25  3.322407  0.722252\n",
              "26  1.746695  2.434644\n",
              "27  1.141096  1.657486\n",
              "28  2.637402 -0.201436\n",
              "29  0.095249  1.297269"
            ],
            "text/html": [
              "\n",
              "  <div id=\"df-f1a558e1-f525-4593-ac77-4e6ec98e2ad9\" class=\"colab-df-container\">\n",
              "    <div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>x</th>\n",
              "      <th>y</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>0.187026</td>\n",
              "      <td>-0.422885</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>0.420190</td>\n",
              "      <td>1.057513</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>3.482497</td>\n",
              "      <td>1.365564</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>3.205376</td>\n",
              "      <td>-0.456993</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>1.199737</td>\n",
              "      <td>1.185130</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>5</th>\n",
              "      <td>2.399816</td>\n",
              "      <td>1.426568</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>6</th>\n",
              "      <td>-0.260194</td>\n",
              "      <td>3.057673</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>7</th>\n",
              "      <td>-0.377020</td>\n",
              "      <td>2.942529</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>8</th>\n",
              "      <td>0.939371</td>\n",
              "      <td>-0.023406</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>9</th>\n",
              "      <td>1.614755</td>\n",
              "      <td>1.816510</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>10</th>\n",
              "      <td>3.030900</td>\n",
              "      <td>1.372515</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>11</th>\n",
              "      <td>0.735238</td>\n",
              "      <td>2.662303</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>12</th>\n",
              "      <td>3.149366</td>\n",
              "      <td>1.908935</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>13</th>\n",
              "      <td>1.644236</td>\n",
              "      <td>2.942909</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>14</th>\n",
              "      <td>2.177196</td>\n",
              "      <td>0.766498</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>15</th>\n",
              "      <td>0.539502</td>\n",
              "      <td>0.451881</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>16</th>\n",
              "      <td>2.265196</td>\n",
              "      <td>2.286545</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>17</th>\n",
              "      <td>1.483374</td>\n",
              "      <td>0.616614</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>18</th>\n",
              "      <td>1.400418</td>\n",
              "      <td>1.602659</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>19</th>\n",
              "      <td>-0.179868</td>\n",
              "      <td>-0.452139</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>20</th>\n",
              "      <td>0.902008</td>\n",
              "      <td>1.540868</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>21</th>\n",
              "      <td>0.697423</td>\n",
              "      <td>2.192745</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>22</th>\n",
              "      <td>2.664715</td>\n",
              "      <td>0.595353</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>23</th>\n",
              "      <td>1.746641</td>\n",
              "      <td>3.348775</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>24</th>\n",
              "      <td>1.733469</td>\n",
              "      <td>1.142788</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>25</th>\n",
              "      <td>3.322407</td>\n",
              "      <td>0.722252</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>26</th>\n",
              "      <td>1.746695</td>\n",
              "      <td>2.434644</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>27</th>\n",
              "      <td>1.141096</td>\n",
              "      <td>1.657486</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>28</th>\n",
              "      <td>2.637402</td>\n",
              "      <td>-0.201436</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>29</th>\n",
              "      <td>0.095249</td>\n",
              "      <td>1.297269</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>\n",
              "    <div class=\"colab-df-buttons\">\n",
              "\n",
              "  <div class=\"colab-df-container\">\n",
              "    <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-f1a558e1-f525-4593-ac77-4e6ec98e2ad9')\"\n",
              "            title=\"Convert this dataframe to an interactive table.\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\" viewBox=\"0 -960 960 960\">\n",
              "    <path d=\"M120-120v-720h720v720H120Zm60-500h600v-160H180v160Zm220 220h160v-160H400v160Zm0 220h160v-160H400v160ZM180-400h160v-160H180v160Zm440 0h160v-160H620v160ZM180-180h160v-160H180v160Zm440 0h160v-160H620v160Z\"/>\n",
              "  </svg>\n",
              "    </button>\n",
              "\n",
              "  <style>\n",
              "    .colab-df-container {\n",
              "      display:flex;\n",
              "      gap: 12px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert {\n",
              "      background-color: #E8F0FE;\n",
              "      border: none;\n",
              "      border-radius: 50%;\n",
              "      cursor: pointer;\n",
              "      display: none;\n",
              "      fill: #1967D2;\n",
              "      height: 32px;\n",
              "      padding: 0 0 0 0;\n",
              "      width: 32px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert:hover {\n",
              "      background-color: #E2EBFA;\n",
              "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "      fill: #174EA6;\n",
              "    }\n",
              "\n",
              "    .colab-df-buttons div {\n",
              "      margin-bottom: 4px;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert {\n",
              "      background-color: #3B4455;\n",
              "      fill: #D2E3FC;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert:hover {\n",
              "      background-color: #434B5C;\n",
              "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "      fill: #FFFFFF;\n",
              "    }\n",
              "  </style>\n",
              "\n",
              "    <script>\n",
              "      const buttonEl =\n",
              "        document.querySelector('#df-f1a558e1-f525-4593-ac77-4e6ec98e2ad9 button.colab-df-convert');\n",
              "      buttonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "      async function convertToInteractive(key) {\n",
              "        const element = document.querySelector('#df-f1a558e1-f525-4593-ac77-4e6ec98e2ad9');\n",
              "        const dataTable =\n",
              "          await google.colab.kernel.invokeFunction('convertToInteractive',\n",
              "                                                    [key], {});\n",
              "        if (!dataTable) return;\n",
              "\n",
              "        const docLinkHtml = 'Like what you see? Visit the ' +\n",
              "          '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
              "          + ' to learn more about interactive tables.';\n",
              "        element.innerHTML = '';\n",
              "        dataTable['output_type'] = 'display_data';\n",
              "        await google.colab.output.renderOutput(dataTable, element);\n",
              "        const docLink = document.createElement('div');\n",
              "        docLink.innerHTML = docLinkHtml;\n",
              "        element.appendChild(docLink);\n",
              "      }\n",
              "    </script>\n",
              "  </div>\n",
              "\n",
              "\n",
              "<div id=\"df-50cadbbb-9a57-4b4d-9ff5-5466c5b46776\">\n",
              "  <button class=\"colab-df-quickchart\" onclick=\"quickchart('df-50cadbbb-9a57-4b4d-9ff5-5466c5b46776')\"\n",
              "            title=\"Suggest charts\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "<svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "     width=\"24px\">\n",
              "    <g>\n",
              "        <path d=\"M19 3H5c-1.1 0-2 .9-2 2v14c0 1.1.9 2 2 2h14c1.1 0 2-.9 2-2V5c0-1.1-.9-2-2-2zM9 17H7v-7h2v7zm4 0h-2V7h2v10zm4 0h-2v-4h2v4z\"/>\n",
              "    </g>\n",
              "</svg>\n",
              "  </button>\n",
              "\n",
              "<style>\n",
              "  .colab-df-quickchart {\n",
              "      --bg-color: #E8F0FE;\n",
              "      --fill-color: #1967D2;\n",
              "      --hover-bg-color: #E2EBFA;\n",
              "      --hover-fill-color: #174EA6;\n",
              "      --disabled-fill-color: #AAA;\n",
              "      --disabled-bg-color: #DDD;\n",
              "  }\n",
              "\n",
              "  [theme=dark] .colab-df-quickchart {\n",
              "      --bg-color: #3B4455;\n",
              "      --fill-color: #D2E3FC;\n",
              "      --hover-bg-color: #434B5C;\n",
              "      --hover-fill-color: #FFFFFF;\n",
              "      --disabled-bg-color: #3B4455;\n",
              "      --disabled-fill-color: #666;\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart {\n",
              "    background-color: var(--bg-color);\n",
              "    border: none;\n",
              "    border-radius: 50%;\n",
              "    cursor: pointer;\n",
              "    display: none;\n",
              "    fill: var(--fill-color);\n",
              "    height: 32px;\n",
              "    padding: 0;\n",
              "    width: 32px;\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart:hover {\n",
              "    background-color: var(--hover-bg-color);\n",
              "    box-shadow: 0 1px 2px rgba(60, 64, 67, 0.3), 0 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "    fill: var(--button-hover-fill-color);\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart-complete:disabled,\n",
              "  .colab-df-quickchart-complete:disabled:hover {\n",
              "    background-color: var(--disabled-bg-color);\n",
              "    fill: var(--disabled-fill-color);\n",
              "    box-shadow: none;\n",
              "  }\n",
              "\n",
              "  .colab-df-spinner {\n",
              "    border: 2px solid var(--fill-color);\n",
              "    border-color: transparent;\n",
              "    border-bottom-color: var(--fill-color);\n",
              "    animation:\n",
              "      spin 1s steps(1) infinite;\n",
              "  }\n",
              "\n",
              "  @keyframes spin {\n",
              "    0% {\n",
              "      border-color: transparent;\n",
              "      border-bottom-color: var(--fill-color);\n",
              "      border-left-color: var(--fill-color);\n",
              "    }\n",
              "    20% {\n",
              "      border-color: transparent;\n",
              "      border-left-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "    }\n",
              "    30% {\n",
              "      border-color: transparent;\n",
              "      border-left-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "      border-right-color: var(--fill-color);\n",
              "    }\n",
              "    40% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "    }\n",
              "    60% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "    }\n",
              "    80% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "      border-bottom-color: var(--fill-color);\n",
              "    }\n",
              "    90% {\n",
              "      border-color: transparent;\n",
              "      border-bottom-color: var(--fill-color);\n",
              "    }\n",
              "  }\n",
              "</style>\n",
              "\n",
              "  <script>\n",
              "    async function quickchart(key) {\n",
              "      const quickchartButtonEl =\n",
              "        document.querySelector('#' + key + ' button');\n",
              "      quickchartButtonEl.disabled = true;  // To prevent multiple clicks.\n",
              "      quickchartButtonEl.classList.add('colab-df-spinner');\n",
              "      try {\n",
              "        const charts = await google.colab.kernel.invokeFunction(\n",
              "            'suggestCharts', [key], {});\n",
              "      } catch (error) {\n",
              "        console.error('Error during call to suggestCharts:', error);\n",
              "      }\n",
              "      quickchartButtonEl.classList.remove('colab-df-spinner');\n",
              "      quickchartButtonEl.classList.add('colab-df-quickchart-complete');\n",
              "    }\n",
              "    (() => {\n",
              "      let quickchartButtonEl =\n",
              "        document.querySelector('#df-50cadbbb-9a57-4b4d-9ff5-5466c5b46776 button');\n",
              "      quickchartButtonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "    })();\n",
              "  </script>\n",
              "</div>\n",
              "\n",
              "  <div id=\"id_c5dfdc73-2920-415b-885d-e4f9c9b837af\">\n",
              "    <style>\n",
              "      .colab-df-generate {\n",
              "        background-color: #E8F0FE;\n",
              "        border: none;\n",
              "        border-radius: 50%;\n",
              "        cursor: pointer;\n",
              "        display: none;\n",
              "        fill: #1967D2;\n",
              "        height: 32px;\n",
              "        padding: 0 0 0 0;\n",
              "        width: 32px;\n",
              "      }\n",
              "\n",
              "      .colab-df-generate:hover {\n",
              "        background-color: #E2EBFA;\n",
              "        box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "        fill: #174EA6;\n",
              "      }\n",
              "\n",
              "      [theme=dark] .colab-df-generate {\n",
              "        background-color: #3B4455;\n",
              "        fill: #D2E3FC;\n",
              "      }\n",
              "\n",
              "      [theme=dark] .colab-df-generate:hover {\n",
              "        background-color: #434B5C;\n",
              "        box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "        filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "        fill: #FFFFFF;\n",
              "      }\n",
              "    </style>\n",
              "    <button class=\"colab-df-generate\" onclick=\"generateWithVariable('df')\"\n",
              "            title=\"Generate code using this dataframe.\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "       width=\"24px\">\n",
              "    <path d=\"M7,19H8.4L18.45,9,17,7.55,7,17.6ZM5,21V16.75L18.45,3.32a2,2,0,0,1,2.83,0l1.4,1.43a1.91,1.91,0,0,1,.58,1.4,1.91,1.91,0,0,1-.58,1.4L9.25,21ZM18.45,9,17,7.55Zm-12,3A5.31,5.31,0,0,0,4.9,8.1,5.31,5.31,0,0,0,1,6.5,5.31,5.31,0,0,0,4.9,4.9,5.31,5.31,0,0,0,6.5,1,5.31,5.31,0,0,0,8.1,4.9,5.31,5.31,0,0,0,12,6.5,5.46,5.46,0,0,0,6.5,12Z\"/>\n",
              "  </svg>\n",
              "    </button>\n",
              "    <script>\n",
              "      (() => {\n",
              "      const buttonEl =\n",
              "        document.querySelector('#id_c5dfdc73-2920-415b-885d-e4f9c9b837af button.colab-df-generate');\n",
              "      buttonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "      buttonEl.onclick = () => {\n",
              "        google.colab.notebook.generateWithVariable('df');\n",
              "      }\n",
              "      })();\n",
              "    </script>\n",
              "  </div>\n",
              "\n",
              "    </div>\n",
              "  </div>\n"
            ],
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "dataframe",
              "variable_name": "df",
              "summary": "{\n  \"name\": \"df\",\n  \"rows\": 30,\n  \"fields\": [\n    {\n      \"column\": \"x\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 1.1337888708647177,\n        \"min\": -0.37701964996148585,\n        \"max\": 3.4824970449019275,\n        \"num_unique_values\": 30,\n        \"samples\": [\n          1.1410962432101135,\n          0.5395019652374079,\n          1.7466413479503187\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"y\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 1.0820566550875095,\n        \"min\": -0.45699307498294683,\n        \"max\": 3.348775369957172,\n        \"num_unique_values\": 30,\n        \"samples\": [\n          1.6574860493823633,\n          0.4518814147120972,\n          3.348775369957172\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    }\n  ]\n}"
            }
          },
          "metadata": {},
          "execution_count": 39
        }
      ],
      "source": [
        "df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "D0m6Hv6MqEge",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "54fe2e09-5292-4747-9c3d-6bd0a1a96dc4"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "<ipython-input-40-82040cfc532e>:3: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
            "  dist =[ sum( (df.iloc[ii-1][j-1]-df.iloc[i-1][j-1])**2 for j in [1,2]) for i in range(1,N+1) for ii in range(1,N+1)]\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "24.390311990225158\n"
          ]
        }
      ],
      "source": [
        "M = 100\n",
        "\n",
        "dist =[ sum( (df.iloc[ii-1][j-1]-df.iloc[i-1][j-1])**2 for j in [1,2]) for i in range(1,N+1) for ii in range(1,N+1)]\n",
        "print(max(dist))\n",
        "M = max(dist)"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "df.columns"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Hzds1m_Yzc9u",
        "outputId": "704d7c90-de55-401e-c152-139aa624bba6"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Index(['x', 'y'], dtype='object')"
            ]
          },
          "metadata": {},
          "execution_count": 41
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "89e0aBt4klX5",
        "outputId": "6952eeeb-8e7a-4e68-d2e4-44a7ed08443c"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "12.768760442733765\n"
          ]
        }
      ],
      "source": [
        "model = AbstractModel()\n",
        "model.i = RangeSet(N)\n",
        "model.j = Set(initialize=[c for c in df.columns])\n",
        "model.k = RangeSet(Kcluster)\n",
        "\n",
        "def rule_mu(model,k,j):\n",
        "  return (df[j].min(), df[j].max())\n",
        "model.mu = Var(model.k,model.j, bounds=rule_mu, within=Reals)\n",
        "\n",
        "model.U = Var(model.i,model.k, within=Binary)\n",
        "def rule_initx(model,i,j):\n",
        "  return df.loc[i-1,j]\n",
        "model.x = Param(model.i,model.j, initialize=rule_initx)\n",
        "model.d = Var(model.i,model.k, bounds=(0,M), within= Reals)\n",
        "\n",
        "def rule_c1(model,i,k):\n",
        "  return model.d[i,k]>= sum( (model.x[i,j]-model.mu[k,j])**2 for j in model.j) - M*(1-model.U[i,k])\n",
        "model.c1 = Constraint(model.i,model.k, rule=rule_c1)\n",
        "\n",
        "def rule_c2(model,i):\n",
        "  return sum(model.U[i,k] for k in model.k)==1\n",
        "model.c2 = Constraint(model.i, rule=rule_c2)\n",
        "\n",
        "def rule_c3(model,k1,k2):\n",
        "  if k1<k2:\n",
        "    return model.mu[k1,'x']<= 0.01+ model.mu[k2,'x']\n",
        "  else:\n",
        "    return Constraint.Skip\n",
        "model.c3 = Constraint(model.k,model.k, rule=rule_c3)\n",
        "\n",
        "def rule_c4(model,i1,i2,k):\n",
        "  if i2>i1 and sum( (model.x[i1,j]-model.x[i2,j])**2 for j in model.j)>M/2:\n",
        "    return model.U[i1,k]+model.U[i2,k]<=1\n",
        "  else:\n",
        "    return Constraint.Skip\n",
        "model.c4 = Constraint(model.i,model.i,model.k, rule=rule_c4)\n",
        "\n",
        "def rule_of(model):\n",
        "  return sum(model.d[i,k] for k in model.k for i in model.i)\n",
        "model.obj = Objective(rule=rule_of, sense=minimize)\n",
        "\n",
        "\n",
        "instance = model.create_instance()\n",
        "st= time.time()\n",
        "results = opt.solve(instance, solver = \"cplex\")\n",
        "print(time.time()-st)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 407
        },
        "id": "EcKCwdwHn5GO",
        "outputId": "7b10748d-6242-4957-c817-f2813b49c52b"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "30\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 600x400 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiYAAAF1CAYAAADcPCGCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABE10lEQVR4nO3deXxU1f3/8dedLJOFJEBICIGQkIQQUFFBNCAWNLYi4gZSUVS0RaXtt/60tt+vsSrar4r2Wy1VW6nVolUrViMqglDABYtAlU2LIUpYEgiQECATskyWOb8/QkYi2UlmSd5PH/fhzL1zZj6HSWbeuffccy1jjEFERETEB9i8XYCIiIhIAwUTERER8RkKJiIiIuIzFExERETEZyiYiIiIiM9QMBERERGfoWAiIiIiPkPBRERERHyGgomIiIj4DAUTERER8RmB3i6gPebPn8/ChQvZs2cPlZWVxMTEMHbsWO6//35GjhzZbLubb76Zl1566aT1AwcOZO/evW1+fZfLRWFhIREREViW1aE+iIiI9ETGGMrKyoiPj8dma36/iF8Fk48//pji4mKSk5OpqqoiNzeXN998kw8++ID8/HzCw8NbbD9w4EAGDRrkvh8bG9uu1y8sLCQhIaFDtYuIiAgUFBQ0+i7+Lr8KJq+99hohISHu+/fffz8PP/wwhw8fZvv27YwePbrF9rNnz+bBBx/s8OtHREQA9f+okZGRHX4eERGRnsbhcJCQkOD+Lm2OXwWTkJAQFi9ezOOPP47D4SA3NxeAmJgY0tLSWm0/f/585s2bR2xsLOeffz6PPPIIKSkpzT7e6XTidDrd98vKygCIjIxUMBEREemA1oZC+N3g14MHD7JhwwZycnJwuVwMGTKEDz/8sNUEFhwczIABAxg0aBB79+7l9ddfZ8yYMezbt6/ZNvPmzSMqKsq96DCOiIhI1/K7YDJnzhxcLhd79uzh2muvZdeuXVx77bXuvRlN+eUvf0lJSQk5OTnk5eWxYMECAI4cOcLChQubbZeVlUVpaal7KSgo6PT+iIiIyLf8LphA/W6gwYMHc++99wKwbds2XnvttWYff/rpp9OrVy/3/ZkzZ7pv5+fnN9vObre7D9vo8I2IiEjX85tgUlJSwssvv0x1dbV73bJly9y3y8vLAUhPTyc9PZ1nnnnGvW3u3LkUFxe77y9atMh9OykpqQurFhERkfawjDHG20W0xe7duxkyZAihoaGkpKQ0OrQSERHBl19+SWJiontQzdy5c91n4FiWhc1mIzk5GWMMeXl5AMTFxbF169Y2nzbscDiIioqitLRUe09ERETaoa3foX6zx6R3797MmDGDAQMGkJeXx/79+0lISOCGG25gw4YNJCYmNtv2kUceYdy4cTgcDvbt20dqaipz5szh888/b/dcJiIiItJ1/GaPiS/QHhORnqW6ro6q2hoigu2a7VnkFLX1O9Sv5jEREfGE0qoq5v3rYxZv/4oal4uk3r25O2M8l6UN83ZpIt2egomIyAnqXC5uXPwGOYeKqTu+Q3nP0aP8fPl7WBZMHqpwItKV/GaMiYiIJ3y0exf/KS5yhxIAA1jAk+vWoqPfIl1LwURE5ARbDu4nsIkrnxpg59EjVNTUeL4okR5EwURE5AS9Q0JxNbNXxB4QgD1QR8BFupKCiYjICa5ISyfAsvjuOTgBlsXU4ac1uTdFRDqPfsNERE4QEx7O7y+5jECbDQsItOo/Jk+P7c//nP897xYn0gNon6SIyHdMHprGmIEDWfp1LkeqKhkVF88FiUnYNJeJSJdTMBERaUJMWDg3nzXK22WI9Dg6lCMiIiI+Q8FEREREfIaCiYiIiPgMBRMRERHxGQomIiIi4jMUTERERMRnKJiIiIiIz1AwEREREZ+hYCIiIiI+Q8FEREREfIaCiYiIiPgMBRMRERHxGbqInx/Z6yjl9W1fkl9aSkqfvlx72hn079XL22WJiIh0GgUTP/Hx7l3c9t7buIzBHF/33KbP+NtV1zBqQLxXaxMREeksOpTjB5y1tdz1z2XUulzUGYPr+FJVW8tdK5bhMqb1JxEREfEDCiZ+YMO+vRytquK78cNlDAWOUr4qLvJKXSIiIp1NwcQPVNbWnNJ2ERERf+FXwWT+/PmceeaZ9O7dG7vdzqBBg5g+fTpffPFFq22ffvppRowYgd1uJzY2lh/96EccPHjQA1WfujHxAwm0Nf1W9QoO5vSY/h6uSEREpGv4VTD5+OOPKS4uJjk5mZSUFPbv38+bb77JhRdeSHl5ebPt7r//fu644w5ycnJITEzk2LFjLFy4kIkTJ1JRUeHBHnRM39AwfjbmPACs4+tsVv2tX44dT2hQkJcqExER6Vx+FUxee+01CgsL2bRpE1999RX33nsvAIcPH2b79u1Ntjl48CCPP/44AHfffTdff/0169evx7Istm/fzoIFC5p9PafTicPhaLR4yx3njuV335/E8H4x9AoO5ozY/vxp8hXcdObZXqtJRESks/nV6cIhISEsXryYxx9/HIfDQW5uLgAxMTGkpaU12WbVqlXU1NSPwZg2bRoAI0eOJDU1lW+++Ybly5fzi1/8osm28+bN46GHHuqCntSrrKlhT+lRosPCiAkLb/GxlmUxdfhpTB1+WpfVI9LdGWP45nAJdS4XadH9CGjmEKmIeI9fBROo3wOyYcMG9/0hQ4awZMkSIiIimnx8QUGB+3ZsbKz7dv/+/fnmm2/Iz89v9rWysrIahRaHw0FCQsKplA/Un03zzL/X89ymz6ioqcECLhqSzLzMS+gXFnbKzy8iJ1u/t4B7Vq8gv7QUgNjwcH4zMZMfpAz1cmUiciK/+3Nhzpw5uFwu9uzZw7XXXsuuXbu49tprKSsra9fzmDbM/WG324mMjGy0dIYFn/+b+Rs+peL4nhwDfLR7F7PeflNzkoh0gZ1HDnPzO9kUHA8lAMXl5fx02RI27S/0YmUi8l1+F0yg/rDG4MGD3WNMtm3bxmuvvdbkY0/cw1FUVHTS7cGDB3dhpSerrqvjzxs/O2l9nTHkHCpmbcEej9Yj0hO88sUW6lyuRnMBGeoHk7+w+XMvVSUiTfGbYFJSUsLLL79MdXW1e92yZcvctxvOyklPTyc9PZ1nnnkGgMzMTAID649YZWdnA/DFF1+wY8cOACZNmuSR+hsUl5dTVu1scpvNssg9dMij9Yj0BDmHiqlrYm9knTF8VVzshYpEpDl+E0zKysq46aab6N27N2eccQaDBw8mKysLgIiICKZOnQpAbm4uubm5HDr+BR8XF8evfvUrAJ544gmGDRtGRkYGxhiGDh3K7bff7tF+9A4JaXZOEpcxxOmifCKdLiEyigDLOmm9zbJIiIryQkUi0hy/CSa9e/dmxowZDBgwgLy8PPbv309CQgI33HADGzZsIDExsdm2jzzyCPPnzyc9PZ1du3YRHh7OrFmzWLNmDeHhLZ8N09nCg4OZmj7CPQ9JA5tl0Tc0lO8np3q0HpGeYOYZZzY5fstlDDeNPMvzBYlIsyzTllGgAtSflRMVFUVpaekpDYQ9Vl3Nz5a9yyf5344niQkL5/krruaMWM3iKtIV/rHtS+Z+tBpnXR0AgTYbv8g4nznnnOvlykR6hrZ+hyqYtENnBZMG/yk6yH+KDhITFs73EpMICgjohCpFpDkOp5N/5e+m1uXi/IREonV6vojHKJh0gc4OJiIiIj1FW79D/WaMiYiIiHR/CiYiIiLiMxRMRERExGcomIiIiIjPUDARERERn6FgIiIiIj5DwURERER8hoKJiIiI+AwFExEREfEZCiYiIiLiMxRMRERExGcEersAke6stKqKd3Jz2FN6lKTefbhy2HAi7XZvlyUi4rMUTES6yBcHD3DT4jcpq3YSYLNR53Lx5Lq1vDp1OiNiYr1dnoiIT9KhHJEu4DKG/1q2hGM11Rig1uXCAGXVTv7r/SXoot4iIk1TMBHpApv2F7K3zIHrOwHEZQy7jx7ly6KDXqpMRMS3KZiIdAGH09ni9lJnlYcqERHxLwomIl1gZP84AiyryW1BNhtnxPb3cEUiIv5BwUSkC/QLC+PHZ49ucttto8fQOyTUwxWJiPgHnZUj0kX++/zvERPei+c3fc7B8mPE94rg1tHncNPIs71dmoiIz7KMTg9oM4fDQVRUFKWlpURGRnq7HPEjtS4XgTbtoBSRnqut36H6pBTxAIUSEZG20aeliIiI+AwFExEREfEZCiYiIiLiM/wmmDzxxBNMnDiRAQMGYLfbSUxMZNasWezcubPFdg8++CCWZTW51NbWeqh6ERERaQu/OV346aefJj8/n2HDhhEaGsquXbv429/+xj//+U9yc3NbPUumX79+pKSkNFpnNTMBloiIiHiH3+wxufXWW9m9ezc5OTns3LmTO++8E4ADBw6wevXqVttfdtllrF+/vtESEBDQYhun04nD4Wi0iIiISNfxm2Dy61//msGDB7vvX3DBBe7bdru91fbZ2dmEhoYyYMAApkyZwubNm1ttM2/ePKKiotxLQkJCx4oXERGRNvGbYHKiuro6nnvuOQCSk5PJzMxs8fEBAQHExcWRlJTEgQMHWLp0KWPHjm01nGRlZVFaWupeCgoKOq0PIj7DGCgvgSN76v+vORcxxvD29hyuf+sfTHr1Re7/cBW7jx7xdlkiPYLfzfxaXl7Oddddx5IlS4iLi2P16tWMGDGi2cd//fXX9OvXj759+wKwYsUKJk2aBMCPf/xjnn/++Ta/tmZ+lW6l8ihsfQ02/BmO7Pp2fZ8hcN7tcOZ1ENrbW9V51X0fruLvX27FZlm4jCHAsrAHBPL6Nddymi7AKNIh3XLm1wMHDjBhwgSWLFlCWloaa9eubTGUAKSlpblDCcAll1xCdHQ0APn5+V1ar4jP2rEKnhwBy7PgyO7G247srl//5Ij6x/Uw24oO8vcvtwLgOv53W50xVNXV8ui/PvZmaSI9gt8Ek23btpGRkcHGjRu54IILWLduHcnJyY0ek5mZSXp6OllZWe51jz/+eKMAsnLlSkpKSgBISkrySO0iPmXHKnj1h1BTCZjjy4mOr6uprH9cDwsnH+zeSUATZ+y5jGHd3gIqa2q8UJVIz+E3wWTq1Kns2bMHgLKyMiZPnkxGRgYZGRnuwzF5eXnk5uayf/9+d7tnn32WpKQkEhMTGTFiBJdccgkA4eHh7jN7RHqMyqPw+k3Hx5G4Wnmwq/5xr99U366HsGh+GgEL0CwDIl3Lb4KJ0+l0396yZQsbNmxwL3v37m223b333ktmZiY1NTXs3LmTxMREZs6cycaNG1s9DCTS7Wx9DWoqaD2UNHDVP37roq6syqd8PyWVuiaG3gVYFuMHJxISGOSFqkR6Dr8b/OpNGvwqfs0YeOrs42NK2vNrb0GfJLhjc4/ZXTDvXx/zl02fNxr8Gh4UzBvTr2Po8TFqItI+bf0O9ZuZX0XkFFUcbnz2TZuZ+naVRyCsb+sP7wbuOf97jIkfyD+2/YeSygrGxA/kpjPPJj5Cf5CIdDUFE5GeovrYqbV3lvWYYGJZFhcnp3Jxcqq3SxHpcfxmjImInKLgXqfW3h7ROXWIiLRAwUSkpwjrWz95WgtnnTTNqm8X2qcrqhIRaUSHckR6Csuqn9F1eVbrj/2u8+b0mIGvPcWRykr++NkGln6znVqXi8whKfz83LEM1MB+8TLtMRHpSc68DoLCaPOvvmWrf/yZM7q0LPGsMqeTa954jZe2buJgeTkllZVk52zjytdfYX9ZmbfLkx5OwUSkJwntDdf+7fjej9Z+/W2ABde+3GOvmdNd/eOr/7D76JFG87XUGUNpVRV/2fSZFysTUTAR6XlSL4aZ/4CgUI7PZfqdBxxfFxQKM9+A1Jav3i3+55P83U3OZFNnDB/u7sgp5SKdR2NMRHqi1IvhF1/Vz+i6YcF3ri6cVD+m5KzrICTKayVK1wkNDHRPHtfUNhFv0k+gSE8V2hsy5tQPiK08Uj9PiT2i/uwbDXTt1qYMTWdF3o6T1lvAVem6VId4lw7liHQyYwx1rrZei8YHWNbxU4kT6/+vUNLtXTo0jSlDhwH11wBquJryOfEDuenMs7xYmYj2mIh0GmdtLU//ez2vfrmVUmcVqX2j+X/njuWytGHeLk2kEZtlMX/SZVyZPpzlO76h1uXiwqQhXJqaRlBAgLfLkx5OF/FrB13ET1py25K3+WD3Tvdxe4v6S+X93/cnMW34aV6tTUTE29r6HapDOSKd4Muig6zalddoMGHDrf/79BP/OrQjIuJFCiYineDzwn3NTvReVF5OoSatEhFpEwUTkU4Qabc3OS8E1B/S6RUc7MlyRET8loKJSCf4fnIqIYGBJ+01CbAsJiQOoU9oqFfqEhHxNwomIp0g0m7nqUmXEWizYbMsAm31v1rxEZE8fNHFXq5ORMR/6HRhkU5ycXIqa26+lbdzv6KovJzTY/ozeWgads2kKSLSZvrEFOlE/Xv14vbR53q7DBERv6VDOSIiIuIzFExERETEZyiYiIiIiM9QMBERERGfoWAiIiIiPsNvgskTTzzBxIkTGTBgAHa7ncTERGbNmsXOnTtbbbto0SJGjRpFaGgoffv25ZprriEvL88DVYuIiEh7+M3VhZOSksjPz2fYsGE4nU527doFQFxcHLm5uc1eqfCFF15g9uzZAAwZMoSSkhIcDgexsbFs3bqVuLi4NtegqwuLiIh0TLe7uvCtt97K7t27ycnJYefOndx5550AHDhwgNWrVzfZprq6mnvuuQeAadOmsXPnTnJycoiIiKCoqIhHH33UU+WLiIhIG/hNMPn1r3/N4MGD3fcvuOAC92273d5km88++4xDhw4B9cEEID4+noyMDACWL1/e4ms6nU4cDkejRURERLqO3wSTE9XV1fHcc88BkJycTGZmZpOPKygocN+OjY113+7fvz8A+fn5Lb7OvHnziIqKci8JCQmnWrqIiIi0wO+CSXl5OVdffTUrVqwgLi6OJUuWNLvHpDltHVaTlZVFaWmpezkx6IiIiEjn86tr5Rw4cIApU6awceNG0tLSeP/990lOTm728Sfu4SgqKjrp9omHhppit9vbHXpERETayhgDrv1AMFZAP2+X4xP8Zo/Jtm3byMjIYOPGjVxwwQWsW7fupFCSmZlJeno6WVlZAIwZM4bo6GgAsrOzASgsLGT9+vUATJo0yYM9EBER+ZZxfow59ANM8URM8ThcJddianK9XZbX+U0wmTp1Knv27AGgrKyMyZMnk5GRQUZGBs8//zwAeXl55Obmsn//fgCCg4PdZ95kZ2eTnJzM8OHDKSsro1+/fu4zdkRERDzJVG/CHLkd6k4Y61izFXP4ekxdsfcK8wF+cyjH6XS6b2/ZsqXRtpb2fNx2222Eh4fzu9/9jpycHEJCQpg6dSqPPfYY8fHxXVWuiIhIs0z5c4AFuE5Y6wJTDpWvQ6//8lJl3uc3E6z5Ak2wJiIincFVNA5ch5reaM/E1udZzxbkAd1ugjUREZFuw9af+j0m3xUAttgm1vccCiYiIiIeZoVdBzR1wKIOK+waT5fjUxRMRKRbKy4v5/0dX/PJnt1U19V5uxyReqHXQOgNx+9Yx5cArMgHsYLO8GJh3uc3g19FRNrDGMPjn37C85s+x3V8KF10aCh/mDSFcQktz2Ek0tUsy4YV9QAm/EZwfgJWMNgv1lwmaI+JiHRTr3y5lec2fuYOJQCHK6v48buLOXjsmBcrE/mWFTgEK/wmrLAZCiXHKZiISLe0cPPGk4YWGgw1rjrezNnmlZpEpHUKJiLSLe0rczQ5tNBmWex1lHq8HhFpGwUTEemWhvTp0+TJmHUuFyl9+nq8HhFpGwUTEemW5ow+96Q9JjbLoldwMNOGn+aVmkSkdQomItItXZU+gnvHTyAsKMi9LimqN69O/SF9QkO9WJmItERT0reDpqQX8T8VNTV8VVxEeHAw6dH9sKymDvCISFdr63eo5jERkW4tLCiIc+IHersMEWkjHcoRERERn6FgIiIiIj5Dh3JExK3W5WLVzjw+3L2TQJuNSalDGZ+QqHEZIuIxCiYiAoCztpYfL1nMpwX5BBwPIq/95wuuTh/B/31/EjaFExHxAB3KEREAXv5iC+sK8gGoM4a64yfsLd7+FSvyvvFmaSLSgyiYiAhQH0Cam8L93dztHq9HRHomHcoREaB+vo+muIyhoqa6y163oLSUN3P+w6GKCk6LieXKYcMJDw7ustcTEd+mYCIiAExITOLVL7e6D+E0sIDzExK75DWXfp3LnSuW1r+OZfGay8Uz/17P69fMICEqqkteU0R8mw7liAgAt44eQ4Td7h74ChBgWQyKjGLG6Wd0+usdrark7pXvu8ez1LpcABRXlHPfhys7/fVExD8omIgIAAMjIln8w5lcMWw4EcF2+oSEct3pI3nzh9cRaQ/p9Nf7Z94OquvqTlpfZwyf5O/hcGVFp7+miPg+HcoR8RPGGCpqaggOCCAoIKBLXiOxd2+e+MGlXfLc31VeU4MFTQ64BaisqQVda0+kx1EwEfED/8z7ht99+i92HDlMkC2Aq9OHc8/479E7xH+/uTMGDmo2lMRHRDAgIsKj9YiIb9ChHBEftzJvB3OWvkvekcMA1LjqyM7Zxg2L33CPy/BHw2NimTJ0GCdO29Ywidv/nP89Tegm0kMpmIj4uCfWrz3pkEedMXxVXMwHu/K8VVaneOIHl/KLsefTP7wXgTYbp8fE8pcpV3F5Wrq3SxMRL/GrYLJmzRomT55MTEwMlmVhWRYLFixotd2DDz7ofvx3l9raWg9ULtIxztpavi451OQhj0CbjU37Cz1eU2cKCgjgZ2MyWPfj2/n6v+7i7Rk3kJmc4u2yRMSL/GqMyaZNm1i5ciXJyckcOnSo3e379etHSkrjDz1dnEx8WVBAACGBgVQ1EaBdxtAn1H/HmIiINKVde0yOHTvGO++8w/Lly6n5ziyRf/3rX/nNb37TqcV914033ojD4WDFihUdan/ZZZexfv36RktAF53dINIZbJbF9BGnNznewgKuSBvu+aJERLpQm4PJrl27GDZsGFOnTuWyyy7jtNNOIzc31739L3/5Cw899FCXFNkgOjqa0FP4CzE7O5vQ0FAGDBjAlClT2Lx5c4uPdzqdOByORouIp/1y7HhG9o8D6g/f2CyLAMvidz+4VGeuiEi30+ZDOQ899BD79+9339+xYweZmZmsXbuWxMSuma66MwUEBBAXF0dgYCDbt29n6dKlrFq1inXr1nH22Wc32WbevHldHrZEWhNht/Pm9OtYs2c3G/fvI8oewuVp6fTv1cvbpYmIdLo27zH56KOPsCyLO++8k+zsbH7wgx9QWFjIpZdeytGjR7uwxFN3/fXXU1RUxDfffENOTg7Lly8H6veI/PGPf2y2XVZWFqWlpe6loKDAUyWLNGKzLCYmDeHuseOZPeochRLxKmOqMBVv4Dr6S1yOBzHVLe99FmmPNu8xOXDgAEOGDOHJJ58E4Morr+Sqq65i6dKlXHXVVTidzi4r8lSlpaU1un/JJZcQHR1NSUkJ+fn5zbaz2+3Y7fauLk9ExG8Y1xFMyfVQlwfUj9EzFX/HhP8MW8T/825x0i20eY9J3759G31J22w2Fi1axNlnn80nn3zCli1buqK+dsvMzCQ9PZ2srCz3uscff7xRAFm5ciUlJSUAJCUlebpEERG/Zcr+AHW7j9+rO74A5X/E1PzHS1VJd9LmYDJs2DByc3PZu3eve11YWBjvvfceCQkJXVLcd7311lukpqYyceJE97oHHniA1NRUZs6cCUBeXh65ubmNxsM8++yzJCUlkZiYyIgRI7jkkksACA8P58477/RI7SIi3ULVu7jDSCMBmMqlnq5GuqE2B5MpU6aQkJDAyy+/3Gh9XFwcS5cuZejQoV0+CNbhcJCXl8eePXvc64qLi8nLy2Pfvn3Ntrv33nvJzMykpqaGnTt3kpiYyMyZM9m4cSMjRozo0ppFRLoVU9XMBguo9GQlPsWYakzlO7iO/gLX0f/GVH2AMf55yQjjKsXU5mGafa+7lmWMae46Wk368MMPufDCC5vc9uKLL3LzzTd3Rl0+yeFwEBUVRWlpKZGRkd4uR0TE41yHfwzVn9LUXhOr9zNYIT/wfFFeZkwl5vDNULOZ+r/3LaAOQi7DinoCy/KPSdaNqwzjeBCqlgIusMIg7BasXv+FZZ36nF9t/Q5t97/W97//fX71q181mmCtpKSEadOmMXv27I5VKyIifsGKuIv68yZO/KKyQdAosF/kpaq8rPwlqNl6/I4Ld2irWgrOf3qrqnYxxmCO/PTbUAJgKurHDh172qO1tDuYuFwunnzySc4991y2bdvGsmXLOOOMM1i8eDGBgX41w72IiLSTFXQ6VvTr9SHEigBbHITPwerzVyyrZ34HmMp3cX+ZN2Lzn3E3tV9CzQaa7Ef5XzGuCo+V0u6fonfffZfbb7+drVu3Mnr0aGpqajDGMHLkSF566aWuqFFERHyIFTQCq0/zc0D1PM1Nl2FaGJPjY2pyWthYBXUFYBvmkVLavcdkypQpbNy4kf79+7tDydlnn82///1vzjzzzK6oUfxERVkl+3cdpNpZ0/qDRUS6C/tEGh/a+pZlv8CjpXSYLbaFjRbY+nmulPY2WLNmDRMmTKCoqMh9Zd4tW7Zw8cUXk5eX1+kFiu+rKKvkdz/+E1Ojb+GmlP/ih3GzefXhbFwu/xyRLiLSHlb4bLAiaRxOAiBgCIRO81ZZ7WMfD7b+nBywAsB+MVZAtMdKafdZOQEBARhjSElJ4ZVXXmHPnj385Cc/4ciRI4SHh1NWVtZVtXqdzspp2n//4Dds/XAbrrrGQeSmuT/kxrnTvVSViIjnmLp9mGMLwLkaCITQKVjht2HZenu7tDYzNdsxR34MrmLqzywyEDQSq89fsGx9Tvn52/od2u5gYrPZuPXWW3nyyScJDw8HoLCwkFmzZvHBBx9QV9fUxDvdg4LJyb7emMfPxtzT5LaQcDtvHHyBkDBN6y8i4g+MqQHnx1B3AIKGQdA57qMjp6qt36HtHvz69ttvc8UVVzRaFx8fz8qVK3nqqafaX6n4tR2bdjW7rarcyf68Aww5w/evPi0iImBZQRBysVdraPcYk++GkhPdcccdp1SM+J8+cb2b32hB79goj9UiIiL+zz+moxOfNWbSWfSN640toPGPki3Axtgp59Cnf2/vFCYiIn5JwUROSWBQIA+/l0Vk314AWLb6Y5HJIxP5xfNzvFmaiIj4oZ45TZ90qqGjknk1fwHr3v2cQ3tLGDIykbMuPA2bTblXRETaR8FEOkWwPYgJ08d6uwzpIsYYVrz4EUueXcGhfSWkjU7h2v+5itPPT/d2aSLSzSiYiEirFtz9Em/NX4plWRhj+Pf7m9mwbBMPLf5vxl5+jrfLE5FuRPvaRaRFhXkHeOsP9Rcia5j2yFXnAmN49q6FmuFXRDqV9piISIs2rvyiyfXGwP6dRRzYVUR8SpyHq+o6xnUYU/4KVH8KVjhW6FUQchmWpb/jRDxBwUREWhQUHAgtzA8dZA/yXDFdzNTtx5RcA64S6i//bsNUfwLOjyDqd502A6aINE9/AohIizIuH01g0MlXTrUF2Bg2JpWYQZ67uFdXM2V/ANdh6kMJ3/6/aglUr/NWWSI9ioKJiLSod0wUP39mNlgQEFj/kWELsBHaK4Q7/3ybl6vrZM4VQFPX+wrEOFd6uhqRHkmHckSkVZNvvZjUUUN4//nVlOw/QsqZSVx2+/fpF9/X26V1shaOWbXveqci0kEKJiLSJmmjU0gbneLtMrqW/WKoWsrJe01qsUIu8kZFIj2Ogol0S7u3FbD+vY1YlsW4K88hYdhAb5ckfsDqdQfG+QkYB/XhxAJMfWAJHu/l6kR6BssY7Z9sK4fDQVRUFKWlpURGRnq7HGmCMYZnf/Eii/+wrP7CggZcLhcz7rmaHz1ync6qkFaZugOY8hehei1YvepPFw6dhmXp7ziRU9HW71D9pkm3suaNdSz+wzLg+CRgxy16bDGnjRtGxpTR3ipN/IQVEIcVeY+3yxDpsXRWjnQrS/+yCpvt5L0itgAb77+w2gsViYhIeyiYSLdytKgUl+vko5OuOheHDxz1fEEiItIufhVM1qxZw+TJk4mJicGyLCzLYsGCBW1qu2jRIkaNGkVoaCh9+/blmmuuIS8vr4srFk8bMXaYe66NE9kCbJw2bpgXKhIRkfbwq2CyadMmVq5cSd++7Zs74YUXXuC6665j8+bNDBgwgLq6OrKzsxk3bhwHDhzoomrFG675xRQCgwLrB74eZwuwYQ8L5qqfX+rFykREpC38KpjceOONOBwOVqxY0eY21dXV3HNP/UC2adOmsXPnTnJycoiIiKCoqIhHH320q8oVLxiUFs8TH/+G4RlD3etOH5/O79f8L3FJsV6sTERE2sKvzsqJjm7/NTk+++wzDh06BNQHE4D4+HgyMjJYuXIly5cvb7at0+nE6XS67zscjna/vnjesHNSmP/Jw5SXloNlER4Z5u2SRESkjfxqj0lHFBQUuG/Hxn77F3P//v0ByM/Pb7btvHnziIqKci8JCQldV6h0uvCocIUSERE/0+2DSXPaMq9cVlYWpaWl7uXEkCMiIiKdz68O5XTEiXs5ioqKTro9ePDgZtva7XbsdnvXFSciIl3GGMNR51EqaisICwyjt723Zn/2A90umGRmZrJv3z6uvvpq5s2bx5gxY4iOjqakpITs7Gyuu+46CgsLWb9+PQCTJk3ycsUiItKZHNUO3t3xLn/f/ncKyr7d050QkcD16ddzReoVRAbrsiK+yq8O5bz11lukpqYyceJE97oHHniA1NRUZs6cCUBeXh65ubns378fgODgYPeZN9nZ2SQnJzN8+HDKysro16+f+4wdERHxf2v3reXiNy7mt5/9lr1lextt21u2l99+9lsufuNi1u5b66UKpTV+FUwcDgd5eXns2bPHva64uJi8vDz27dvXbLvbbruNV155hbPOOovCwkIsy2Lq1Kl8+umnxMfHe6J0ERHpYmv3reWnq39KVW0V5vh/J2pYV1VbxU9X/1ThxEfp6sLtoKsLi4j4Jke1g4vfuNgdSlpjYRESGMKq6at0WMdD2vod6ld7TERERJry7o532xxKAPeekyV5S7q4MmkvBRMREfFrxhj+vv3vHWr7as6rbZo+QjxHwURERPzaUedRCsoK2ry3pIHBUFBWQKmztIsqk45QMBEREb9WUVtxSu3La8s7qRLpDAomIiLi18ICT+3SE+GB4Z1UiXQGBRMRkW7ImBqM6zDG1Hm7lC7X296bhIgELNo3q6uFRUJEAlH2qC6qTDpCwUREpBsxphqX43FM0TmYogxM0ThM+QvdeoCnZVlcn359h9rOHD5T09T7GAUTEZFuxJT+N1QsBFN5fMURTNnjUP6MdwvrYlekXkFIYEib95rYsBESGMLlKZd3cWXSXgomIiLdhKndCVXLANfJ2479BePqvoM8I4Mj+f3E32NZVqvhxMICC+ZPnK/J1XyQgomISHdR82ULG6ugLs9jpXjD+QPP50+Zf3LvOfluQGlYFxIYwrOZzzJu4DgvVdr5TN1BXI55uIon4zo0DVP+N4yp9nZZHdLtri4sItJj2fq2vN1qZXs3cP7A81k1fRVL8pbwas6rja4uPChiEDOHz+SKlCuICI7wYpWdy9Ttw5RMA1cpUAdYmLL/gPND6PM8lhXg7RLbRcFERKS7CB4LthhwldD4cE4ABI3CChzkrco8KjI4kpnDZ3J9+vWUOkspry0nPDCcKHtUtxzoasr+eEIoARommqteC85VEHKJt0rrEB3KERHpJiwrEKvPs2A17A04/pdyQDxW7996rS5vsSyL3iG9GdhrIL1DenfLUALUhw+aOi08AOP80NPVnDLtMRER6UasoJEQ8zFULQdXIQQOBftFWFaQt0uTLtPSPgb/2/+gYCIi0s1YtjAIm+rtMsRTQi+FikWcvNekDivkB96o6JT4X5QSERERNyv8ZxAQx7df6cf/b78Ugr/nrbI6THtMRMQv1NXV4ayoJrRXSPcdKyDSAVZAP4h+GyoWYao/AcKwQq+AkEuxLP/b/6BgIiI+rbamlpcfeoN3/ric8tIK+g7ow4z/uYqrfn6pAorIcZYtCnrdjsXt3i7llCmYiIhPe2L2s6x+5RP3tV4O7z/Cn+5cSIWjkpn3TfNydeIpxhio2QzVG8AKh5BJWAGx3i5LuoD/7eMRkR6jMO8Aq15e0+QF6F57bDGVxyq9UJV4mjHVmKNzMIdnYI49hSl7BFM8AVOR7e3SpAsomIiIz9q+4ZtmtzkrnOzetteD1YjXlP8ZnB8dv1NH/QRidRjHvZjaXd6rS7qEgomI+KyIvr1a3B4Z3fJ26R5MxSLcs5k2YsNULvZ0OdLFFExExGednXkGfeN6Y7M1/qiyBdgYnpHGwNQBXqpMPMpV2swG6/j0+9KdKJiIiM8KDApk7lu/IjQiBCwICKqfYj06vg/3vPxzL1cnHhN0Ok1/XdViBZ3p6Wqki1mmqVFl0iSHw0FUVBSlpaVERkZ6uxzxAGMMRypqKHfWEm4PpE9YkE5R9YJyRwUfLVrLwT3FJI5I4IJp5xEcEuztssRDjHMt5siPGu4d/38ABMRhRS+tn+lWfF5bv0MVTNpBwaTnKK2sIXvjXl76dDd7Dle41yf2DWPWuCSmjR5EVKiuPSLiKabqQ0zZb6EuD7CB/WKsyPuwAuK8XZq0UVu/Q/3uUM6iRYsYNWoUoaGh9O3bl2uuuYa8vLwW29x8881YlnXSMmhQz7gEuLTPx18XM3beav73va/IPyGUAOQfruB/3/uKsfNW8/HXxV6qUKTnsUIuxOq3DCt2A1bsJmx9nlEo6ab8aoK1F154gdmzZwMwZMgQSkpKyM7O5pNPPmHr1q3ExbX8Qzpw4MBGYSQ2VpPzSGMff13MLQv/jaHpcwAa1lXW1HHLwn+z8JZzmZAW48EKRXouy7LA6uPtMqSL+c0ek+rqau655x4Apk2bxs6dO8nJySEiIoKioiIeffTRVp9j9uzZrF+/3r28++67LT7e6XTicDgaLdJ9lVbW8JNXNtaHklYOcBpTH1J+8spGSitrPFGeiEiP4DfB5LPPPuPQoUNAfTABiI+PJyMjA4Dly5e3+hzz58/HbreTkJDAjBkzWj0ENG/ePKKiotxLQkLCKfZCfNWRolLuvPdNKpy1rYaSBsZAZXUdb23SJF8iIp3Fb4JJQUGB+/aJh2D69+8PQH5+fovtg4ODGTBgAIMGDWLv3r28/vrrjBkzhn379jXbJisri9LSUvdyYg3SfRw7Ws7Px97Lvypaf2xTXly7u8kp00VEpP38Jpg0py1fCL/85S8pKSkhJyeHvLw8FixYAMCRI0dYuHBhs+3sdjuRkZGNFul+lj63iv1FDmr6REI7TwU2wJ7DFRyt0OEcEZHO4DfB5MTDKEVFRSfdHjx4cLNtTz/9dHr1+nbq6pkzZ7pvt7anRbq/jSu3UhdwauPAjzlrO6kaEZGezW+CyZgxY4iOjgYgO7v+ipKFhYWsX78egEmTJgGQnp5Oeno6zzzzjLvt3LlzKS7+9tTORYsWuW8nJSV1deni40J7hRBQd2rBopfdr05wExHxWX4TTIKDg91n3mRnZ5OcnMzw4cMpKyujX79+7jN2cnNzyc3NdQ+UBfjNb35DXFwcQ4cOJTU1lVtvvRWAuLg49+nH0nNdOON8rPIqgo44wNW+sSIW9ZOu9Q7TZGsiIp3Bb4IJwG233cYrr7zCWWedRWFhIZZlMXXqVD799FPi4+ObbffII48wbtw4HA4H+/btIzU1lTlz5vD5559rLhPhe9PHctF144nalFOfNNrp5vOTNE29iEgn0ZT07aAp6bsvYwwfv7eJW/9VSA1WmwbB2iwICQpgXVampqcXEWlFt52SXqQrWJbFxMtH8/yPM7DZrFZzScP2BTeMVigREelECiYiJ5iQFsPCW84lNCgAi5OP7DSsCw0K4MVbzuV7mo5eRKRT6VQCke+YkBbDuqxM3tq0lxfXNr668OC+Ydx8fv3VhSNDtKdERKSzaYxJO2iMSc9jjOFoRQ3HnLX0sgfSOyxIA11FRDqgrd+h2mMi0gLLsugTHkyf8GBvlyIi0iNojImIiIj4DAUTERER8RkKJiIiIuIzFExERETEZyiYiIiIiM9QMBERERGfoWAiIiIiPkPBRERERHyGgomIiIj4DAUTERER8RkKJiIiIuIzFExERETEZyiYiIiIiM9QMBERERGfoWAiIiIiPkPBRERERHyGgomIiIj4DAUTERER8RkKJiIiIuIzFExERETEZyiYiIiIiM/wu2CyaNEiRo0aRWhoKH379uWaa64hLy+v1XZPP/00I0aMwG63Exsby49+9CMOHjzogYpFRESkrSxjjPF2EW31wgsvMHv2bACGDBlCSUkJDoeD2NhYtm7dSlxcXJPt7r//fh5++GEAhg4dyt69e6msrCQ9PZ2NGzcSFhbWptd3OBxERUVRWlpKZGRk53RKpAN2bN7Fiw8sYvPqL7GH2cmceQE3PfhDIvr08nZpIiJNaut3qN/sMamuruaee+4BYNq0aezcuZOcnBwiIiIoKiri0UcfbbLdwYMHefzxxwG4++67+frrr1m/fj2WZbF9+3YWLFjQ7Gs6nU4cDkejRcTbdmzZxf87/9d8tnwL1VU1lB0+xrt/WsHdE+dSXVXt7fJERE6J3wSTzz77jEOHDgH1wQQgPj6ejIwMAJYvX95ku1WrVlFTU9Oo3ciRI0lNTW2xHcC8efOIiopyLwkJCZ3TGZFT8PJDb1BbU4erzuVe56pzsevLfD56/VMvViYicur8JpgUFBS4b8fGxrpv9+/fH4D8/PxObQeQlZVFaWmpeznxuUS8ZcsH/2kUShoEBNrY8tF/vFCRiEjnCfR2Aaeqo0Nk2tLObrdjt9s79PwiXSUk3E5FWWUTWyzCeoV6vB4Rkc7kN3tMTjyMUlRUdNLtwYMHd2o7EV/1g1kTsQWc/KtbV1vHRdeP90JFIiKdx2+CyZgxY4iOjgYgOzsbgMLCQtavXw/ApEmTAEhPTyc9PZ1nnnkGgMzMTAIDAxu1++KLL9ixY0ejdiL+4rp7pzJ0VDIAAYEBBAQGADDz19MYMXaYN0sTETllfnW68HPPPcftt98OND5duF+/fmzdupX4+HgsywJg7ty5PPjggwDce++9zJs3D4C0tDQKCgqorKxk6NChbN68mfDw8Da9vk4XFl9RW1PL2sX/ZvPqLwkJt3PhdeMZNibV22WJiDSrrd+hfjXG5LbbbiM8PJzf/e535OTkEBISwtSpU3nssceIj49vtt0jjzxC//79WbBgAXl5eURFRfHDH/6Qxx57rM2hRMSXBAYFMuGH45jww3HeLkVEpFP51R4Tb9MeExERkY7pdhOsiYiISPenYCIiIiI+Q8FEREREfIZfDX4VEZGew9TmY8qfA+e/wBaOFXo1hN2IZWniy+5MwURERHyOqd2JKZkOpgKoAxeYsv8D5xrosxDLCvB2idJFdChHRER8jin7w7eh5Nu1UL0enB96qyzxAAUTERHxPdUf0TiUNAjEOD/2cDHiSQomIiLig4KaWW/Aam6bdAcKJiIi4ntCpgBNjSOpwwq51NPViAcpmIiIiM+xev0cAuIB6/ia4yEldAYEneOtssQDdFaOiIj4HCsgGqLfgco3MM51YOuFFXI52Ce6L9Yq3ZOCiYiI+CTL1gvCb8EKv8XbpYgH6VCOiIiI+AwFExEREfEZCiYiIiLiMxRMRERExGcomIiIiIjPUDARERERn6FgIiIiIj5D85iIRx0pKmXVy2s4uLuIxBGDuOj68YRHhXu7LBER8REKJuIxX6z5insnP0p1VTW2ABt1tXW8NPd1fvfhQySdluDt8kRExAfoUE43U1lexYeL1vLOH5ezY/Mub5fjVltTy2+mP0F1VTXGZairqQMDZUfKmTfzDxhjvF2iiIhPMaYS4yr3dhkepz0m3cjGlVv5zfQnqHBUYllgDIy9/Bzue/0ugkOCvVrb5g/+Q2mx46T1rjoXO7/YQ/72fSQOH+SFykREfIup3YFxPArVawGDCToHKzILK+gMb5fmEdpj0k2UHnIw96rfUllWBdSHEoD1Szfy0gOve7GyehWlFS1uLz/a8/4qEBH5LlN3AFMyA6rXAcc/yGs2YUpmYmp3erU2T1Ew6SY+fG0t1c6akw6JGJdhyZ//icvl8lJl9UaMG9bsFUFDwu0MGZno4YpERHyPqXgFzDGg7oS1LqAGU/5XL1XlWX4TTMrKyrjrrrsYNGgQwcHBpKSk8NBDD1FbW9tqW8uymlzuu+8+D1TuGSWFh7EFNP12VpZV4ays9nBFjcUMiuaKn10CJ2aT47dvfGA6oeEhXqlLRMSnVG+mPoh8Vx1Uf+7parzCL8aYuFwuLr/8cj7++GOCgoJITk7mm2++4cEHHyQvL4+//e1vbXqes846C7vd7r6fkNB9zgRJOSupfkBpE+KGxBISZm9ymyf9dP4t9B8cw1tPLeXQ3sMMTI1jxv9czSW3XOjt0kREfIOtLxBA4z0mALbj27o/y/jB6RBvvfUW06ZNA2DJkiVMmTKFp59+mjvuuAOAjRs3MmrUqGbbNxxC2LVrF0lJSR2uw+FwEBUVRWlpKZGRkR1+nq5QU13D7NPu4sDuYlx1jdP23S/8lEk+9uVvjGn20I6ISE9lnB9hjtzW5DYr6nGs0Ks9XFHnaet3qF8cynn//fcBCA0NZfLkyQDuoAKwfPnyNj3POeecQ1hYGKeddhqPPfYYTqezxcc7nU4cDkejxVcFBQfxxEcPMfr7I92HSKL6RfDzZ2b7XCgBFEpERJoSPAHCbz9+x0b93hMg9FoIudJbVXmUXxzKKSgoACA6OhqbrT5L9e/f3709Pz+/1efo06cPgwYNYs+ePXz11VdkZWWxefNmXn+9+TNW5s2bx0MPPXSK1XtOv4HRPLrs1xwtLuXYkXLihsQSGOQXb7GIiHB8TGTE3ZjQK6FqJeAC+4VYQSO8XZrHeHWPyT333NPswNSGZfv27U22bc8RqPXr11NSUsKWLVvYt28fF110EQD/+Mc/3KGnKVlZWZSWlrqXlh7rS3rHRDEoLV6hRETET1mBqVi9foLV62c9KpSAl/eYjB49mlmzZrX4mKioKPcg1UOHDuFyubDZbBQVFbkfM3jw4Baf47zzznPfDgsL4+qrr+aDDz4A6vfGNDcI1m63NxosKyIiIl3Lq8Fk+vTpTJ8+vdXHTZo0ieeff56qqiqWLVvGlClTyM7ObrQdYPHixWRlZQGwevVqBg4cyJo1aygqKuLqq68mICCAqqoq3nnnHXfbxETNnyEiIuIr/GJf/1VXXcX48eP517/+xdSpU0lJSeHrr78G4Prrr3efkVNaWkpubi4ANTU1AOzcuZNbbrmF8PBwkpOT2bt3L0eOHAHglltuYeDAgV7okYiIiDTFL87KCQgIYOnSpdxxxx3ExMSQl5fH4MGDeeCBB3jxxRdbbDt+/HjmzJnD4MGD2bVrFy6Xi9GjR7NgwQKee+45z3RARERE2sQv5jHxFb48j4mIiIgv61bzmIiIiEjPoGAiIiIiPkPBRERERHyGX5yV4ysahuP48tT0IiIivqjhu7O1oa0KJu1QVlYGdK+rEouIiHhSWVkZUVFRzW7XWTnt4HK5KCwsJCIi4qSL0DkcDhISEigoKOh2Z+yob/6nu/YL1Dd/pb75n87ulzGGsrIy4uPj3de9a4r2mLSDzWZj0KBBLT4mMjKyW/1gnkh98z/dtV+gvvkr9c3/dGa/WtpT0kCDX0VERMRnKJiIiIiIz1Aw6SR2u525c+d2y6sRq2/+p7v2C9Q3f6W++R9v9UuDX0VERMRnaI+JiIiI+AwFExEREfEZCiYiIiLiMxRMRERExGcomHRQWVkZd911F4MGDSI4OJiUlBQeeughamtrW21rWVaTy3333eeByhtbtGgRo0aNIjQ0lL59+3LNNdeQl5fXarunn36aESNGYLfbiY2N5Uc/+hEHDx70QMVt15G+3XzzzU2+N61NrOdJa9asYfLkycTExLjrW7BgQZvadvT99oSO9uvBBx9s9neqLb+PnvDEE08wceJEBgwYgN1uJzExkVmzZrFz585W2/ryewYd75s/vG/z58/nzDPPpHfv3tjtdgYNGsT06dP54osvWm3r65+RHe2bRz4jjbRbXV2dmTBhggFMUFCQGTZsmLHZbAYwN954Y6vtAQOYs846y5x33nnuZcGCBR6o/lvPP/+8u5YhQ4aYyMhIA5jY2Fizf//+Ztvdd9997nZDhw41oaGhBjDp6emmvLzcgz1oXkf7NmvWLAOYgQMHNnpvLr/8cg9W37Lf//73JjAw0KSlpbn7+Oyzz7barqP/Jp7S0X7NnTvXAKZfv36N3rPzzjvP1NbWeqDy1iUmJhrLskx6eroZMmSIu39xcXGmtLS02Xa+/p4Z0/G++cP7dtVVV5kBAwaYs88+2wwfPtz9Od+3b19z7NixZtv5w2dkR/vmic9IBZMOyM7Odv/QLVmyxBhjzFNPPeVet3HjxhbbNzxu165dHqi2aU6n0/Tr188AZtq0acYYY/bt22ciIiIMYH7+85832e7AgQMmKCjIAObuu+82xhizdetWY1mWAcwTTzzhsT40p6N9M+bbX7q5c+d6qNr2O3TokKmoqDC7du1q8xf4qfybeEpH+mXMt19ws2bN6voiO+jhhx82e/bscd+/88473X186623mmzjD++ZMR3rmzH+8b5VVlY2un9i4Pj888+bbOMPn5HGdKxvxnjmM1KHcjrg/fffByA0NJTJkycDMG3aNPf25cuXt+l5zjnnHMLCwjjttNN47LHHcDqdnV9sMz777DMOHToEfFt7fHw8GRkZQPN9WLVqFTU1NY3ajRw5ktTU1BbbeVJH+3ai+fPnY7fbSUhIYMaMGT616zw6OprQ0NB2temMf5Ou1pF+nSg7O5vQ0FAGDBjAlClT2Lx5cydWd2p+/etfM3jwYPf9Cy64wH27ucmr/OE9g4717US+/L6FhISwePFiMjIyGDFiBI8++igAMTExpKWlNdnGHz4joWN9O1FXfkYqmHRAQUEBUP9B2nCFxP79+7u35+fnt/ocffr0YdCgQdjtdr766iuysrK46aabuqbgJjT0ASA2NtZ9u6EfzfWho+086VRrDA4OZsCAAQwaNIi9e/fy+uuvM2bMGPbt29c1BXuAP7xvpyIgIIC4uDiSkpI4cOAAS5cuZezYsT71Jdegrq6O5557DoDk5GQyMzObfJw/vmdt7VsDf3jfDh48yIYNG8jJycHlcjFkyBA+/PBDIiIimny8P71v7e1bg67+jFQwOcE999zT7GCshmX79u1NtjXtmEB3/fr1lJSUsGXLFvbt28dFF10EwD/+8Y9GP9Te0J5+dEY7T2pLjb/85S8pKSkhJyeHvLw89+DLI0eOsHDhwq4u0eP84X1rzfXXX09RURHffPMNOTk57r9InU4nf/zjH71cXWPl5eVcffXVrFixgri4OJYsWdLu6b599T1rb9/85X2bM2cOLpeLPXv2cO2117Jr1y6uvfZaysrK2vU8vvi+daRvnviMVDA5wejRo5k1a1aLS1RUFAkJCQAcOnQIl8sFQFFRkft5Ttyt2ZTzzjsPy7IACAsL4+qrr3Zv81QwaegDNK694XZzfehoO086lRpPP/10evXq5b4/c+ZM921f+kunvfzhfeuotLQ0+vbt675/ySWXEB0dDfjWe3bgwAEmTJjAkiVLSEtLY+3atYwYMaLZx/vTe9bevoH/vG9Qfybl4MGDuffeewHYtm0br732WpOP9af3DdrXN/DMZ6SCyQmmT5/Oiy++2OIyYMAAJk2aBEBVVRXLli0D6o+TNmjYvnjxYtLT00lPT3fv4lqzZg1vvvkmdXV17ud455133G0TExM90tcxY8a4PwQaai8sLGT9+vWN+tBQ/zPPPANAZmYmgYGBjdp98cUX7Nixo1E7b+po3wDmzp1LcXGx+/6iRYvct5OSkrq69E6TmZlJeno6WVlZQNv/TXzdd/sF8Pjjjzf6QFy5ciUlJSWA77xn27ZtIyMjg40bN3LBBRewbt06kpOTGz3GX9+zjvQNfP99Kykp4eWXX6a6utq9ruHzHur3EIF/fkZ2tG/goc/ILhtW243V1taa8ePHu08XTk9Pd59qdf3117sft3DhwpPOwGlYFx4ebs444wzTp08f92NuueUWj/bjz3/+c5OnIvbr18/s27fPGPPtGUQnjsDOyspyr09LS3OfCjd06NAWTzPzpI72DTA2m82kpqaalJSURqc+Hjx40Eu9aSw7O9ukpKSYxMREd30xMTEmJSXF/fPXsO3EMx7a8m/iTR3tV8PpqoMHDzbDhw93n/0QHh5utm3b5qXeNHbiKdDfnSbgL3/5izHGP98zYzreN19/3xrODgsNDTWnn366SUhIcPczIiLC7N692xjjn5+Rp9I3T3xGKph0UGlpqbnjjjtMfHy8CQoKMklJSeaBBx4w1dXV7sc0FUy++eYbM2fOHDN8+HDTq1cvExUVZUaPHm0WLFhgampqPN6PV155xZx11lnGbrebqKgoM3XqVPP111+7tzf1g+lyucz8+fNNenq6CQoKMv369TOzZs3ymXkVGnSkb4888ogZP368iY2NNSEhISY1NdXMmTPH7N271ws9aNqJP1ffXSZMmGCMafqLwJjW/028qaP9+vOf/2wuvvhiM2DAAGO3201SUpKZOXOm2b59u3c60oQTw9Z3l4afP398z4zpeN98/X07cuSImTFjhklOTjahoaEmMDDQJCQkmBtuuMF89dVX7sf542fkqfTNE5+R1vEXFxEREfE6jTERERERn6FgIiIiIj5DwURERER8hoKJiIiI+AwFExEREfEZCiYiIiLiMxRMRERExGcomIiIiIjPUDARERERn6FgIiLdxpo1a5g8eTIxMTFYloVlWe7LsouIf1AwEZFuY9OmTaxcuZK+fft6uxQR6SAFExHxC++88w6WZWGz2fjoo48AeO+999zrVq9ezY033ojD4WDFihXeLVZEOkzBRET8wpVXXsns2bMxxnDrrbeyf/9+br/9dgDuvPNOMjMziY6OJjQ01MuVisipUDAREb/x+9//npSUFHbs2MHZZ59NYWEhZ5xxBvPmzfN2aSLSSRRMRMRv9OrVi1deeQWbzcbBgwcJCgri1VdfxW63e7s0EekkCiYi4lfy8/NxuVwA1NTUsGfPHi9XJCKdScFERPzGvn37mDNnDgBnnXUWALNnz6a4uNiLVYlIZ1IwERG/YIzh5ptv5siRI4wbN45169YxcuRIDh48yG233QbAW2+9RWpqKhMnTnS3e+CBB0hNTWXmzJleqlxE2iPQ2wWIiLTFU089xapVqwgNDWXhwoWEhITw0ksvce655/L222/z17/+FZvNRl5eXqN2xcXFFBcXM2jQIC9VLiLtYRljjLeLEBEREQEdyhEREREfomAiIiIiPkPBRERERHyGgomIiIj4DAUTERER8RkKJiIiIuIzFExERETEZyiYiIiIiM9QMBERERGfoWAiIiIiPkPBRERERHzG/wfVzCQw9f4mbQAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ],
      "source": [
        "plt.figure(figsize=(6,4))\n",
        "\n",
        "all= [sum([value(k*instance.U[i,k]) for k in instance.k]) for i in instance.i]\n",
        "print(len(all))\n",
        "for i in instance.i:\n",
        "  pass\n",
        "  #print(\"TAG300\", i, [value(k*instance.U[i,k]) for k in instance.k] )\n",
        "\n",
        "plt.scatter(df[\"x\"],df[\"y\"], c=all, s=  20)\n",
        "\n",
        "for k in instance.k:\n",
        "  xk,yk = value(instance.mu[k,'x']), value(instance.mu[k,'y'])\n",
        "  plt.scatter(xk,yk, s=  100)\n",
        "plt.xlabel(\"x1\", fontweight='bold')\n",
        "plt.ylabel(\"x2\", fontweight='bold')\n",
        "plt.xticks(fontweight='bold')\n",
        "plt.yticks(fontweight='bold')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Silhouette Score"
      ],
      "metadata": {
        "id": "AT5ACMgRsUmb"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "1. For each data point, the average distance (a_i) to other data points within\n",
        "the same cluster is calculated. This value represents the similarity level of the data point to others in its cluster.\n",
        "\n",
        "2. For each data point, the average distance (b_i) to all other clusters it doesn’t belong to is computed. This value indicates how different the data point is from data points in other clusters."
      ],
      "metadata": {
        "id": "4CTwanSas6VJ"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Silhouette Score = (b_i — a_i) / max(a_i, b_i)\n",
        "\n"
      ],
      "metadata": {
        "id": "jBUzbWFqsXwL"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def dist(i,n):\n",
        "  return sum( (df.loc[i-1,j]-df.loc[n-1,j])**2 for j in df.columns)"
      ],
      "metadata": {
        "id": "x4Ao5mXOtXJ-"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "8PhFEynbbsVH",
        "outputId": "8a166f97-e84b-412c-d67a-099ba45da519"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0.6979704055226915\n"
          ]
        }
      ],
      "source": [
        "from re import S\n",
        "result = {}\n",
        "for i in instance.i:\n",
        "  cluster = int(sum([value(k*instance.U[i,k]) for k in instance.k]))\n",
        "  result[i] = cluster\n",
        "result\n",
        "\n",
        "silhouette = {}\n",
        "for i in instance.i:\n",
        "  k = result[i]\n",
        "  ai = sum([dist(i,n) for n in instance.i if result[n]==k])/ len([n for n in instance.i if result[n]==k])\n",
        "  bi = sum([dist(i,n) for n in instance.i if result[n]!=k])/ len([n for n in instance.i if result[n]!=k])\n",
        "  #print(i,k,ai,bi, (bi-ai)/max(ai,bi))\n",
        "  si = (bi-ai)/max(ai,bi)\n",
        "  silhouette[i] = si\n",
        "S = sum(silhouette.values())/len(silhouette)\n",
        "print(S)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "lrxX4hxvbsXt",
        "outputId": "f162eae3-2cc2-474e-fe11-d6ad9af715cf"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(           x         y\n",
              " 0   0.187026 -0.422885\n",
              " 1   0.420190  1.057513\n",
              " 2   3.482497  1.365564\n",
              " 3   3.205376 -0.456993\n",
              " 4   1.199737  1.185130\n",
              " 5   2.399816  1.426568\n",
              " 6  -0.260194  3.057673\n",
              " 7  -0.377020  2.942529\n",
              " 8   0.939371 -0.023406\n",
              " 9   1.614755  1.816510\n",
              " 10  3.030900  1.372515\n",
              " 11  0.735238  2.662303\n",
              " 12  3.149366  1.908935\n",
              " 13  1.644236  2.942909\n",
              " 14  2.177196  0.766498\n",
              " 15  0.539502  0.451881\n",
              " 16  2.265196  2.286545\n",
              " 17  1.483374  0.616614\n",
              " 18  1.400418  1.602659\n",
              " 19 -0.179868 -0.452139\n",
              " 20  0.902008  1.540868\n",
              " 21  0.697423  2.192745\n",
              " 22  2.664715  0.595353\n",
              " 23  1.746641  3.348775\n",
              " 24  1.733469  1.142788\n",
              " 25  3.322407  0.722252\n",
              " 26  1.746695  2.434644\n",
              " 27  1.141096  1.657486\n",
              " 28  2.637402 -0.201436\n",
              " 29  0.095249  1.297269,\n",
              " 30)"
            ]
          },
          "metadata": {},
          "execution_count": 46
        }
      ],
      "source": [
        "df,N"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# K-Memoids"
      ],
      "metadata": {
        "id": "Gtuo0_n07ijl"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "TL-0w1h5bsaL",
        "outputId": "3d76b33f-b74a-41cd-86b4-a8795d87ec8e"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "309.7168815135956\n"
          ]
        }
      ],
      "source": [
        "model = AbstractModel()\n",
        "model.i = RangeSet(N)\n",
        "model.j = Set(initialize=[c for c in df.columns])\n",
        "model.k = RangeSet(Kcluster)\n",
        "\n",
        "def rule_mu(model,k,j):\n",
        "  return (df[j].min(), df[j].max())\n",
        "\n",
        "model.mu = Var(model.i,model.k, within=Binary)\n",
        "model.U = Var(model.i,model.k, within=Binary)\n",
        "model.y = Var(model.i,model.i,model.k, bounds=(0,1), within=Reals)\n",
        "\n",
        "def rule_initx(model,i,j):\n",
        "  return df.loc[i-1,j]\n",
        "model.x = Param(model.i,model.j, initialize=rule_initx)\n",
        "\n",
        "def rule_c2(model,i):\n",
        "  return sum(model.U[i,k] for k in model.k)==1\n",
        "model.c2 = Constraint(model.i, rule=rule_c2)\n",
        "\n",
        "def rule_c3(model,k):\n",
        "  return sum(model.mu[i,k] for i in model.i)==1\n",
        "model.c3 = Constraint(model.k, rule=rule_c3)\n",
        "\n",
        "def rule_c4(model,i,k):\n",
        "  return model.mu[i,k] <= model.U[i,k]\n",
        "model.c4 = Constraint(model.i, model.k, rule=rule_c4)\n",
        "\n",
        "\n",
        "def rule_c5A(model,i1,i2,k):\n",
        "  if i2!=i1:\n",
        "    return model.y[i1,i2,k] <= model.mu[i1,k]\n",
        "  else:\n",
        "    return Constraint.Skip\n",
        "model.c5A = Constraint(model.i,model.i, model.k, rule=rule_c5A)\n",
        "\n",
        "def rule_c5B(model,i1,i2,k):\n",
        "  if i2!=i1:\n",
        "    return model.y[i1,i2,k] <= model.U[i2,k]\n",
        "  else:\n",
        "    return Constraint.Skip\n",
        "model.c5B = Constraint(model.i,model.i, model.k, rule=rule_c5B)\n",
        "\n",
        "def rule_c5C(model,i1,i2,k):\n",
        "  if i2!=i1:\n",
        "    return model.y[i1,i2,k] >= model.mu[i1,k]+ model.U[i2,k]-1\n",
        "  else:\n",
        "    return Constraint.Skip\n",
        "model.c5C = Constraint(model.i,model.i, model.k, rule=rule_c5C)\n",
        "\n",
        "\n",
        "def rule_c5D(model,i1,i2,k):\n",
        "  if i2!=i1:\n",
        "    return model.y[i1,i2,k] <= model.U[i1,k]\n",
        "  else:\n",
        "    return Constraint.Skip\n",
        "model.c5D = Constraint(model.i,model.i, model.k, rule=rule_c5D)\n",
        "\n",
        "\n",
        "def rule_c6(model,k1,k2):\n",
        "  if k1<k2:\n",
        "    return sum(model.x[i,'x']*model.mu[i,k1] for i in model.i) <= 0.1+ sum(model.x[i,'x']*model.mu[i,k2] for i in model.i)\n",
        "  else:\n",
        "    return Constraint.Skip\n",
        "model.c6 = Constraint(model.k,model.k, rule=rule_c6)\n",
        "\n",
        "def rule_of(model):\n",
        "  return sum( model.y[i1,i2,k]*(model.x[i1,j]-model.x[i2,j])**2 for j in model.j for i1 in model.i for i2 in model.i for k in model.k)\n",
        "model.obj = Objective(rule=rule_of, sense=minimize)\n",
        "\n",
        "\n",
        "instance = model.create_instance()\n",
        "st= time.time()\n",
        "results = opt.solve(instance, solver = \"cplex\", options={\n",
        "    'timelimit': 600,\n",
        "    'mipgap': 0.1\n",
        "})\n",
        "print(time.time()-st)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "jsZoDZ5Dbsc1",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 963
        },
        "outputId": "b92fbdd9-d720-4fac-c987-2fbd867c0f56"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "TAG300 1 1\n",
            "TAG300 2 1\n",
            "TAG300 3 3\n",
            "TAG300 4 3\n",
            "TAG300 5 1\n",
            "TAG300 6 3\n",
            "TAG300 7 2\n",
            "TAG300 8 2\n",
            "TAG300 9 1\n",
            "TAG300 10 2\n",
            "TAG300 11 3\n",
            "TAG300 12 2\n",
            "TAG300 13 3\n",
            "TAG300 14 2\n",
            "TAG300 15 3\n",
            "TAG300 16 1\n",
            "TAG300 17 2\n",
            "TAG300 18 1\n",
            "TAG300 19 2\n",
            "TAG300 20 1\n",
            "TAG300 21 2\n",
            "TAG300 22 2\n",
            "TAG300 23 3\n",
            "TAG300 24 2\n",
            "TAG300 25 3\n",
            "TAG300 26 3\n",
            "TAG300 27 2\n",
            "TAG300 28 2\n",
            "TAG300 29 3\n",
            "TAG300 30 1\n",
            "1 0.5395019652374079 0.4518814147120972\n",
            "2 0.6974232674574891 2.1927449116706126\n",
            "3 2.6647145217672126 0.5953532423897839\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 600x400 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiYAAAF1CAYAAADcPCGCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABPPUlEQVR4nO3deXxU5d3//9eZyWQy2UlICEsSsgABNwSpqFhpEaSoVUGtggu2Cty/+65r21uoC7Su31tbqvYWW1u1KmJvES0FhbhULKK1uKDIGpawJwRIQpbJZOb8/hgYCclkI5mZM7yfj0cezFmuM59PDpl8cs51XccwTdNEREREJALYwh2AiIiIyFEqTERERCRiqDARERGRiKHCRERERCKGChMRERGJGCpMREREJGKoMBEREZGIocJEREREIoYKExEREYkYKkxEREQkYsSEO4COmDt3Ls899xzbt2+nrq6OjIwMzjnnHO69915OP/30oO2mTp3KCy+80Gx937592blzZ7vf3+fzsXv3bpKSkjAMo1M5iIiInIxM06S6upo+ffpgswW/LmKpwuSDDz6gvLyc/Px86uvr2bBhA6+99hrvvfcepaWlJCQktNq+b9++9OvXL7CcmZnZofffvXs32dnZnYpdREREYMeOHU1+Fx/PUoXJK6+8QlxcXGD53nvv5YEHHuDAgQOsX7+e4cOHt9r+5ptvZvbs2Z1+/6SkJMD/TU1OTm6yzePxsHz5csaNG4fD4ej0e0Qi5WY90ZoXKDerUm7W09V5VVVVkZ2dHfhdGoylCpO4uDgWLVrEo48+SlVVFRs2bAAgIyODgQMHttl+7ty5PPzww2RmZnLeeefx4IMPUlBQEHR/t9uN2+0OLFdXVwPgcrlwuVxN9o2JiSE+Ph6XyxVV/zFBuVlRtOYFys2qlJv1dHVeHo8HoM2uEIZpmuYJv1sIzZs3j//4j/8ILOfl5bF48WJOOeWUoG2mTp3KK6+8Qn5+Pg0NDWzZsgWAHj168NVXX9G3b98W282ePZs5c+Y0Wz9//nzi4+NPMBMREZGTR21tLZMnT6aysrLZXYdjWa4wAX8Hmh07dvCLX/yCV199lVNOOYVVq1YFvTz09ddf079/fxITEwF45plnmDFjBgC//vWvueeee1psd/wVk6OXofbv39/irZzi4mLGjh0bVRUzKDcrita8QLlZlXKznq7Oq6qqip49e7ZZmFjqVs5RhmGQk5PDrFmzePXVV1m7di2vvPIK06ZNa3H/U089tcnylClTAoVJaWlp0PdxOp04nc5m6x0OR9CT1No2q1Nu1hOteYFysyrlZj1dlVd7j2GZeUwqKip48cUXaWhoCKxbunRp4HVNTQ0ARUVFFBUV8dRTTwW23X///ZSXlweWFyxYEHjdv3//boxaREREOsIyV0yqq6u54YYbmD59OgUFBVRWVrJjxw7AP1pm4sSJAIEOsfv37w+0/dWvfsUDDzxAfn4+pmlSUlICQFZWFjfffHOIMxEREZFgLHPFJDU1lWuuuYbevXtTUlLCnj17yM7O5rrrruOTTz4hNzc3aNsHH3yQc889l6qqKnbt2kVhYSEzZszg3//+d4fnMhEREZHuY5krJqmpqbzyyitt7tdSX95Zs2Yxa9as7ghLRKJYg9dLfaOHpFinZnsWCRHLFCYiIqFSWV/Pw//8gEXrv8Hj89E/NZW7Ro7i4oGDwh2aSNRTYSIicgyvz8f1i/6PdfvL8R65Arv90CF++vbfMQyYMEDFiUh3skwfExGRUPjHtq18XV4WKEoATMAAfrNqZYu3i0Wk66gwERE5xhf79hDTwpNPTWDLoYPUHplWW0S6hwoTEZFjpMa58AW5KuK023HG6A64SHdSYSIicowfDizCbhgcPwbHbhhMHHxKi1dTRKTr6CdMROQYGQkJ/Paii4mx2TCAGMP/MXlqZi/++7zvhjc4kZOArkmKiBxnwoCBjOjblyUbN3Cwvo5hWX04P7c/Ns1lItLtVJiIiLQgIz6BqUOHhTsMkZOObuWIiIhIxFBhIiIiIhFDhYmIiIhEDBUmIiIiEjFUmIiIiEjEUGEiIiIiEUOFiYiIiEQMFSYiIiISMVSYiIiISMRQYSIiIiIRQ4WJiIiIRAwVJiIiIhIx9BA/C9lZVcmra7+itLKSgh5p/OiU0+iVmBjusERERLqMChOL+GDbVqb9/Q18pol5ZN0fPvuUv1x+JcN69wlrbCIiIl1Ft3IswN3YyB3Ll9Lo8+E1TXxHvuobG7lj2VJ8ptn2QURERCxAhYkFfLJrJ4fq6zm+/PCZJjuqKvmmvCwscYmIiHQ1FSYWUNfoOaHtIiIiVmGpwmTu3LmcccYZpKam4nQ66devH1dddRVr1qxps+2TTz7JkCFDcDqdZGZm8uMf/5h9+/aFIOoTN6JPX2JsLZ+qxNhYTs3oFeKIREREuoelCpMPPviA8vJy8vPzKSgoYM+ePbz22mt873vfo6amJmi7e++9l1tvvZV169aRm5vL4cOHee655xg9ejS1tbUhzKBz0lzx/OeIswEwjqyzGf5XPztnFC6HI0yRiYiIdC1Ljcp55ZVXiIuLCyzfe++9PPDAAxw4cID169czfPjwZm327dvHo48+CsBdd93FY489xpo1axg6dCjr169n3rx53HnnnS2+n9vtxu12B5arqqoA8Hg8eDxNb58cXT5+fVf5jzPPIjshiRfXfM7O6iryUntw09DhXJhf0G3veVR35xZO0ZpbtOYFys2qlJv1dHVe7T2OYZrWGtKxaNEiHn30UaqqqtiwYQM+n4+MjAxKSkpISkpqtv/LL7/MddddB8BHH33EOeecA8DAgQPZtGkTY8eOZfny5S2+1+zZs5kzZ06z9fPnzyc+Pr4LsxIREYlutbW1TJ48mcrKSpKTk4PuZ6krJuC/AvLJJ58ElvPy8li8eHGLRQnAjh07Aq8zMzMDr3v16sWmTZsoLS0N+l4zZ85scjWlqqqK7Oxsxo0b1+yb6vF4KC4uZuzYsTjauLXiM03+sPpTnv/yM2o9Hgzggtw87h89hnSXq9W24dCR3KwmWnOL1ryg87l9unsXs//xLjuqKgHISEjgl6Mu4Pt5Bd0VaofpvFlTtObW1XkdvevQFssVJjNmzGD69Ons2LGDX/ziF7z66qv86Ec/YtWqVUGLk5a050KR0+nE6XQ2W+9wOIKepNa2HfW/n37Cbz/9uMm6d7ZvpfTvb7D42usD/UciTXtys6pozS1a84KO5bbl4AF+/Pc38Hi9gWH3uw4f5j+XLeWvV14TcZMU6rxZU7Tm1lV5tfcYlur8epRhGOTk5DBr1iwA1q5dyyuvvNLivtnZ2YHXZWVlzV7n5OR0Y6TNNXi9PLP602brvabJuv3lrNyxPaTxiJwMXlrzBV6fr8lcQCb+zuR/+vzfYYpKRFpimcKkoqKCF198kYaGhsC6pUuXBl4fHZVTVFREUVERTz31FABjxowhJsZ/YWjhwoUArFmzhs2bNwMwfvz4kMR/VHlNDdUN7ha32QyDDfv3hzQekZPBuv3leFu4Suo1Tb4pLw9DRCISjGUKk+rqam644QZSU1M57bTTyMnJYebMmQAkJSUxceJEADZs2MCGDRvYf+QXfFZWFj//+c8BePzxxxk0aBAjR47ENE0GDBjA9OnTQ5pHalxc0DlJfKZJlh7KJ9LlspNTsLdwi9RmGGSnpIQhIhEJxjKFSWpqKtdccw29e/empKSEPXv2kJ2dzXXXXccnn3xCbm5u0LYPPvggc+fOpaioiK1bt5KQkMCNN97IihUrSEhICGEWkBAby8SiIc36kdgMgzSXi7H5hSGNR+RkMOW0M1p8ppTPNLnh9KGhD0hEgrJM59fU1NSg/UiO1VKnVsMwuO2227jtttu6I7QOu+e732PP4Wo+LP22P0m6K55nf3gFzhjLnBIRyzgjqzcPjxnH/f94F7fXC0CMzcadI8/jQv0xIBJR9FswDBJjY3nh8iv5umwfX5ftIyM+ge/m9sdht4c7NJGodfUppzG+cCD/LN1Go8/Hedm5pGs+IpGIo8IkjE7N7MWpmXrOjUioJDudTBgwKNxhiEgrLNPHRERERKKfChMRERGJGCpMREREJGKoMBEREZGIocJEREREIoYKExEREYkYKkxEREQkYqgwERERkYihwkREREQihgoTERERiRgqTERERCRi6Fk5It2osr6eNzesY3vlIfqn9uCyQYNJdjrDHZaISMRSYSLSTdbs28sNi16jusGN3WbD6/Pxm1UreXniVQzJyAx3eCIiEUm3ckS6gc80+a+liznsacAEGn0+TKC6wc1/vbUY0zTDHaKISERSYSLSDT7bs5ud1VX4jitAfKbJtkOH+KpsX5giExGJbCpMRLpBldvd6vZKd32IIhERsRYVJiLd4PReWdgNo8VtDpuN0zJ7hTgiERFrUGEi0g16xsfzkzOHt7ht2vARpMa5QhyRiIg1aFSOSDf5xXnfJSMhkWc/+zf7ag7TJzGJW4afxQ2nnxnu0EREIpYKE5FuYjMMfnLmcH5y5nAafT5ibLpAKSLSFn1SioSAihIRkfbRp6WIiIhEDBUmIiIiEjFUmIiIiEjEsExh8vjjjzN69Gh69+6N0+kkNzeXG2+8kS1btrTabvbs2RiG0eJXY2NjiKIXERGR9rDMqJwnn3yS0tJSBg0ahMvlYuvWrfzlL39h+fLlbNiwgeTk5Fbb9+zZk4KCgibrjCATYImIiEh4WKYwueWWW7j++uvJyckB4I477mDu3Lns3buXd999lyuuuKLV9hdffDHPP/98h97T7XbjPmZq8aqqKgA8Hg8ej6fJvkeXj18fDZSb9URrXqDcrEq5WU9X59Xe4ximRR9z+vrrrzNp0iQAlixZwoQJE1rcb/bs2cyZM4fExEQaGxtJTU1l+PDh/PrXv+bMM1uf6Opo2+PNnz+f+Pj4E09CRETkJFFbW8vkyZOprKxs9S6HJQsTr9fLxRdfzLJly8jPz+ebb77B6XS2uO/s2bN54IEHyMvLIyYmhvXr1wPgdDpZtWpVq8VJS1dMsrOz2b9/f7Nvqsfjobi4mLFjx+JwOLogy8ih3KwnWvOC0ORmmiZLN23k9Q3fcLCulmG9+3LD6UPJSUntlvc7SufNmqI1t67Oq6qqip49e7ZZmFjmVs5RNTU1XHvttSxbtoysrCwWL14ctCgBmDx5MrfeeitpaWkALFu2jPHjx+N2u/n973/Ps88+G7St0+ls8dgOhyPoSWptm9UpN+uJ1ryge3O75/13mP/Vl9gMA59psv7AARauX8erV/6IU0LwAEadN2uK1ty6Kq/2HsMyo3IA9u7dywUXXMDixYsZOHAgK1euZMiQIa22GThwYKAoAbjoootIT08HoLS0tFvjFRHrWVu2j/lffQmA78gFZa9pUu9t5KF/fhDO0EROCpYpTNauXcvIkSNZvXo1559/PqtWrSI/P7/JPmPGjKGoqIiZM2cG1j366KNNCpDi4mIqKioA6N+/f0hiFxHreG/bFuwtjNjzmSardu6gLso6OIpEGssUJhMnTmT79u0AVFdXM2HCBEaOHMnIkSMDt2NKSkrYsGEDe/bsCbR7+umn6d+/P7m5uQwZMoSLLroIgISEBG6//faQ5yEikc0g+DQCBqBZBkS6l2UKk2M7oX7xxRd88sknga+dO3cGbTdr1izGjBmDx+Nhy5Yt5ObmMmXKFFavXt3mbSAROfmMLSjE28KYALthMConl7iY6OtDIBJJLNP5ddu2bZ3aZ9q0aUybNq3rAxKRqDQovSe3DDuLP37270DnV7thkOCI5Z7zvxfu8ESinmUKExGRULn7vO8yok9f/rr2ayrqahnRpy83nHEmfZJan2FaRE6cChMRkeMYhsGF+YVcmF8Y7lBETjoqTES6ktcD7mrwecFmB2cS2NUnQUSkvVSYiJyouoNQvhEOlEDVLvDUgWn6h284XJDcF9IKIGMguHqEO1oRkYimwkSks9zVsG0l7P4c6g5ATBzEJkBcKhg2MH3gqYWydf59XGnQ50zof57/SopIGB2sq+P3n37Ckk3rafT5GJNXwE+/cw5923hSu0h3U2Ei0hnlG2HTcji4DRIzIaPIX4wc72gBYvqgtgI2vwMVm2HAOP8VFJEwqHa7ufL/XqG08lBgaPTCdWt5Z2sJi6+5nt5JKpwlfCwzj4lIxNjzJaz5K1Tv9RckCRktFyXHMmz+/TKK/O3W/NV/HJEw+Os3X7Pt0MEm87V4TZPK+nr++NmnYYxMRIWJSMeUb4RvFgM+SC/wd3DtCJvd3w6f/zjlG7sjSpFWfVi6jZYeK+81Td7ftjXk8YgcS4WJSHu5q/23bxrrISX7xI6Vku0/zqbl/uOKhJArJgZbkLn1XTG6wy/hpcJEpL22rfT3KenRv2uO1yPXf7xtK7vmeCLtdMmAosCTk49lAJcX6VEdEl4qTETao+6gf2RNYmabt29M02zxWSvN2GL8x9v9OdQd6po4RdrhBwMGcsmAQYD/GUBHn6Z8Vp++3HDG0DBGJqJROSLtU77RPyQ4oyjoLm6fyZM76nh5n5vKRpNCl43bsl1c3NMZ/Ljx6VC+Hso3QM7Z3RC4SHM2w2Du+Iu5rGgwb2/eRKPPx/f65/GDwoE47B3sNyXSxVSYiLTHgRL/PCWtjL756YbDvHfQg+/Ickmdj59urKHeB5MygxQnhs1/3AMlKkwkpGyGwZi8AsbkFYQ7FJEmdCtHpC1ej39G19iEoLt8dbiRd44pSoDAqIf/2V7b+q2d2ASo2u1/HxGRk5wKE5G2uKv908w74oPu8u+qRloe4wBlHpPdbl+QrfiP66nV6BwREVSYiLTN5z3y7JvgPy7JMUaL80KAf6RDoj1Y2cKR6etN//uIiJzkVJiItMVm9z+Qzwx+1WNsmoM4G82umtiBC1Id9HC08qNm+vzH7+hkbSIiUUiFiUhbnEn+pwR7aoPukhxj44mBicQY/h+qmCMVSh+njQcKgt8CAvzHdcTrwX4iImhUjkjb7A5I7ut/SnArLkyLZcXwVN4od1PW4OPUhBgm9IzFaWvlNg5AQw1kDva/j4jISU6FiUh7pBX4J0Izfa32NekVa2N6X1f7j2v6/FPTp2nIpogI6FaOSPtkDARXGtRWdO1xayv8x80Y1LXHFRGxKBUmIu3h6gF9zoTDZV03esbX6D9enzPBldo1xxQRsTgVJiLt1f88/wP8Dm7rmuMd3O4/Xv/zuuZ4IiJRQIWJSHs5k2DAOP8U8pU7TuxYlTv8xxkwTqNxRESOocJEpCMyBsKQSwEbVJR0/LaOr9HfDpv/OBkDuyNKERHLskxh8vjjjzN69Gh69+6N0+kkNzeXG2+8kS1btrTZdsGCBQwbNgyXy0VaWhpXXnklJSUlIYhaolLvM+D0qyEpy/9k4JryVidfA/zba8r9TxFOyvK3731GaOIVEbEQywwXfvLJJyktLWXQoEG4XC62bt3KX/7yF5YvX86GDRtITk5usd2f/vQnbr75ZgDy8vKoqKhg4cKFfPjhh3z55ZdkZWWFMg2JFhkDIbk3bFvpH0Zcvt5/ayY2wT9ZmmHzFyOeWv88JY31/tE3hRf6+5To9o2ISIssc8XklltuYdu2baxbt44tW7Zw++23A7B3717efffdFts0NDRw9913AzBp0iS2bNnCunXrSEpKoqysjIceeihU4Us0cibBoPEwcgaceqV/kjRbDNQf8g8Drj/kX84c7N8+8j/8+6soEREJyjJXTH75y182WT7//POZO3cuAE6ns8U2n376Kfv37wf8hQlAnz59GDlyJMXFxbz99tutvqfb7cbtdgeWq6qqAPB4PHg8TR9Rf3T5+PXRQLm1ISYReg/zf3k90HDY3/fEZofYxKYzuoboe6hzZk3KzZqiNbeuzqu9xzFM0wz2UNSI5fV6ufjii1m2bBn5+fl88803LRYnCxYs4NprrwXgnXfeYcyYMQBcf/31vPTSSzidTurr64O+z+zZs5kzZ06z9fPnzyc+vo3nn4iIiEhAbW0tkydPprKyMmj3C7DQFZOjampquPbaa1m2bBlZWVksXrw46BWTYNpbi82cOZM777wzsFxVVUV2djbjxo1r9k31eDwUFxczduxYHI7oeuaJcrOeaM0LlJtVKTfr6eq8jt51aIulCpO9e/dyySWXsHr1agYOHMhbb71Ffn5+0P2zs7MDr8vKypq9zsnJafX9nE5ni0WPw+EIepJa22Z1ys16ojUvUG5WpdyaMk0TfHuAWAx7z+4J7AR11Tlr7zEs0/l17dq1jBw5ktWrV3P++eezatWqZkXJmDFjKCoqYubMmQCMGDGC9PR0ABYuXAjA7t27+fjjjwEYP358CDMQERH5lun+AHP/OMzy0Zjl5+Kr+BGmZ0O4wwo7yxQmEydOZPv27QBUV1czYcIERo4cyciRI3n22WcBKCkpYcOGDezZsweA2NjYwMibhQsXkp+fz+DBg6murqZnz56BETsiIiKhZDZ8hnlwOnhLv13p+RLzwGRMb3n4AosAlrmVc+zomC+++KLJttaufEybNo2EhAQee+wx1q1bR1xcHBMnTuSRRx6hT58+3RWuiIhIUGbNHwADOHZyRh+YNVD3KiT+V5giCz/LFCbbtm3r9D5TpkxhypQpXRuQiIhIZ3nWAC090sKH6VmLEep4IohlbuWIiIhEDVsvaLH8sIMtM9TRRBQVJiIiIiFmxF8LtDR1hRcj/spQhxNRVJiISFQrr6nhrc0b+XD7Nhq8HXwatEh3cV0JruuOLBhHvuwYybMxHKeFMbDws0wfExGRjjBNk0c/+pBnP/s3viOTKqa7XPxu/CWcm936HEYi3c0wbBgp92EmXA/uD8GIBeeFETuXSSjpiomIRKWXvvqSP6z+NFCUAByoq+cnf1vEvsOHwxiZyLeMmDyMhBsw4q9RUXKEChMRiUrPfb66WddCExOPz8tr69aGJSYRaZsKExGJSruqq1rsWmgzDHZWVYY8HhFpHxUmIhKV8nr0aHEwptfno6BHWsjjEZH2UWEiIlFpxvDvNLtiYjMMEmNjmTT4lLDEJCJtU2EiIlHp8qIhzBp1AfHHPNG0f0oqL0+8mh4uVxgjE5HWaLiwiEStm4edxeTTzuCb8jISYmMpSu+JYZzMk32LRD4VJiIS1eIdDs7q0zfcYYhIO+lWjoiIiEQMFSYiIiISMXQrR0QCGn0+3tlSwvvbthBjszG+cACjsnPVL0NEQkaFiYgA4G5s5CeLF/HRjlLsRwqRV75ewxVFQ/ifseOxqTgRkRDQrRwRAeDFNV+wakcpAF7TxHvkGTOL1n/DspJN4QxNRE4iKkxEBPAXIMGmcP/bhvUhj0dETk66lSMiANR6PC2u95kmtZ6GbnvfHZWVvLbua/bX1nJKRiaXDRpMQmxst72fiEQ2FSYiAsAFuf15+asvA7dwjjKA87Jzu+U9l2zcwO3LlvjfxzB4xefjqX99zKtXXkN2Skq3vKeIRDbdyhERAG4ZPoIkpzPQ8RXAbhj0S07hmlNP6/L3O1Rfx13FbwX6szT6fACU19Zwz/vFXf5+ImINKkxEBIC+ScksunoKPxw0mKRYJz3iXFx76um8dvW1JDvjuvz9lpdspsHrbbbea5p8WLqdA3W1Xf6eIhL5dCtHxCJM06TW4yHWbsdht3fLe+SmpvL4uB90y7GPV+PxYECLHW4B6jyNoGftiZx0VJiIWMDykk089tE/2XzwAA6bnSuKBnP3qO+SGmfd39wj+/YLWpT0SUqid1JSSOMRkcigWzkiEa64ZDMzlvyNkoMHAPD4vCxct5brFv1foF+GFQ3OyOSSAYM4dtq2o5O4/fd539WEbiInKRUmIhHu8Y9XNrvl4TVNvikv572tJeEKq0s8Pu4H3HnOefRKSCTGZuPUjEz+eMnlXDqwKNyhiUiYWKowWbFiBRMmTCAjIwPDMDAMg3nz5rXZbvbs2YH9j/9qbGwMQeQineNubGRjxf4Wb3nE2Gx8tmd3yGPqSg67nf8cMZJVP5nOxv+6gzeuuY4x+QXhDktEwshSfUw+++wziouLyc/PZ//+/R1u37NnTwoKmn7o6eFkEskcdjtxMTHUt1BA+0yTHi7r9jEREWlJh66YHD58mDfffJO3334bz3GzRP75z3/mV7/6VZcGd7zrr7+eqqoqli1b1qn2F198MR9//HGTL3s3jW4Q6Qo2w+CqIae22N/CAH44cHDogxIR6UbtvmKydetWRo0axd69ewEoKChg8eLFDBo0CIA//vGP/Otf/+K+++7rnkiB9PT0E2q/cOFCXn31VVJTUxk+fDi//vWvOfPMM4Pu73a7cbvdgeWqqioAPB5Ps8Ls6PLx66OBcguv20eMZH3ZPtbs20uMzYbPNDGAB78/lp5xcS3GboW8Oku5WZNys56uzqu9xzFM0ww2Yq+JqVOn8pe//KXJuj59+rBy5Upyc3M555xz+Ne//oW3hQmTutq2bdvIy8sD4Omnn2bGjBmt7j979mweeOAB8vLyiImJYf16/wPJnE4nq1atClqczJ49mzlz5jRbP3/+fOLj408wCxERkZNHbW0tkydPprKykuTk5KD7tbsw6d+/Pzt27OC2227j/PPP55lnnmH58uUUFRXx0Ucf8YMf/CBiC5ONGzfSs2dP0tLSAFi2bBnjx48H4Cc/+QnPPvtsi+1aumKSnZ3N/v37m31TPR4PxcXFjB07FofD0encIpFys55ozQuUWyQwTTfUv4XZ8G+wJWA4L8aIPb3VNlbJrTOiNbeuzquqqoqePXu2WZi0+1bO3r17ycvL4ze/+Q0Al112GZdffjlLlizh8ssvb/ILPNIMHDiwyfJFF11Eeno6FRUVlJaWBm3ndDpxOp3N1jscjqAnqbVtVqfcrCda8wLlFi6m7yBmxWTwlgBH+ug1vAQJ/4kt6bY220dybicqWnPrqrzae4x2d35NS0tr8kvaZrOxYMECzjzzTD788EO++OKLDgfZHcaMGUNRUREzZ84MrHv00UebFCDFxcVUVFQA/itBIiLSPmb178C77ciS98gXUPN7TM/XYYpKokm7C5NBgwaxYcMGdu7cGVgXHx/P3//+d7Kzs7sluOO9/vrrFBYWMnr06MC6++67j8LCQqZMmQJASUkJGzZsYM+ePYF9nn76afr3709ubi5DhgzhoosuAiAhIYHbb789JLGLiESF+r8RKEaasGPWLQl1NBKF2l2YXHLJJWRnZ/Piiy82WZ+VlcWSJUsYMGAAubm5XR7gsaqqqigpKWH79u2BdeXl5ZSUlLBr166g7WbNmsWYMWPweDxs2bKF3NxcpkyZwurVqxkyZEi3xiwiElXM+iAbDKAulJFEHLN+Kb5Dd+I79AvM+vcwTWs+MsL0VWI2lvj7EoVBu/uY3HXXXdx11128//77zbadcsopbNiwgeeff74rY2tm6tSpTJ06tdV9tm3b1mzdtGnTmDZtWvcEJSJyMok9Bxo+ovlVk0aM2HPDEVHYmUeKNbPyHojxAAZm/RsQdzGkPI5hWGOSddNXjVk1G+qXAD5Mbw9gFqbpBULXd6bD362xY8fy85//vMl45IqKCiZNmsTNN9/cpcGJiEhkMZLuwP837bGTU9rAMQyc3w9TVGFWO/+YBR+Boq1+CbiXhyOiDjNNE/Pg/xcoSvwra/3/1DwT0lg6XJj4fD5+85vf8J3vfIe1a9eydOlSTjvtNBYtWkRMjKVmuBcRkQ4yHKdipL/qL0KMJLBlQcIMjB5/xjBOzt8BZt3SIFts1ul30/gVeD4hUJQcq+YlTF9tyELp8P+iv/3tb0yfPp0vv/yS4cOH4/F4ME2T008/nRdeeKE7YhQRkQhiOIZg9Ph9uMOIIA1B1put9MmJMJ51rWysB+8OsA0KSSgdvmJyySWXsHr1anr16hUoSs4880z+9a9/ccYZZ3RHjGIRtdV17Nm6jwZ3dE3LLCLSKueooJsM5/khDOQE2DJb2WiArWfoQulogxUrVnDBBRdQVlYWeDLvF198wYUXXkhJSUmXByiRr7a6jsd+8r9MTL+JGwr+i6uzbublBxbi81mzR7qISEcY8TceeXVsvxs72PPANSkcIXWccxTYetE0h6PbRmPYT+xZdR3R4cLke9/7Hps2bSI/P5+VK1eyYMECUlNT+ec//8nQoUO7IUSJdLMn/Q/Ff/kAb6O/w1dNZS3P37eAl3+9MMyRiYh0P8Pey//CdZn/yoItCxJuwkhfgGFLCG9w7WQYDowefwRb2tE1325LvjeksXS4MDFNk1tuuYUvvviCs88+m6uvvpqvvvqKMWPGUFsbus4xEhk2ri7h83e+wudtfnXkr4+9SX1t5D6qQESkK9mS78GW+RG2zBXYkn6BYUsNd0gdYjiKMDL+gZH6vxhJ92L08D9HLtR5dLjz6xtvvMEPf/jDJuv69OlDcXExTzzxRJcFJtaw+bOtQbfV17jZU7KXvNO6d+I9ERHpGobhgLgL/a89HiDYiKPu0+ErJscXJce69dZbTygYsZ4eWanBNxqQmpkSslhERMT6rDEdnUSsEeOHkpaVis3e9L+SzW7jnEvOokev1PAEJiIilqTCRE5IjCOGB/4+k+S0RAAMm7/DVP7pudz57IxwhiYiIhZ0ck7TJ11qwLB8Xi6dx6q//Zv9OyvIOz2Xod87BZtNda+IiHSMChPpErFOBxdcdU64w5BuYpomy57/B4ufXsb+XRUMHF7Aj/77ck49ryjcoYlIlFFhIiJtmnfXC7w+dwmGYWCaJv9663M+WfoZcxb9gnMuPSvc4YlIFNG1dhFp1e6Svbz+O/+DyEzTBPDPW2OaPH3Hc5rhV0S6lK6YiEirVhevaXG9acKeLWXs3VpGn4KsEEfVfUzfAcyal6DhIzASMFyXQ9zFGIb+jhMJBRUmItIqR2wMmK1sdzpCF0w3M717MCuuBF8F/se/2zAbPgT3PyDlscDzwUSk++hPABFp1chLhxPjaP5gL5vdxqARhWT0C93DvbqbWf078B3AX5Tw7b/1i6FhVbjCEjmpqDARkValZqTw06duBgPsMf6PDJvdhisxjtufmRbm6LqYexngbWFDDKa7ONTRiJyUdCtHRNo04ZYLKRyWx1vPvkvFnoMUnNGfi6ePpWeftLYbW0or96zMVraJSJdRYSIi7TJweAEDhxeEO4zu5bwQ6pfQ/KpJI0bc98MRkchJR4WJRKVta3fw8d9XYxgG5152FtmD+oY7JLEAI/FWTPeHYFbhL04MwPQXLLGjwhydyMlBhYlEFdM0efrO51n0u6X+Bwua8OzdL3HN3Vfw4wev1agKaZURkwM938SseR4aVoKR6B8u7Jqk4cIiIaLCRKLKiv9bxaLfLQWOTAJ2xIJHFnHKuYMYecnwcIUmFmHYszCS7w53GCInLf0JIFFlyR/fwWZrflXEZrfx1p/eDUNEIiLSESpMJKocKqvE52s+esLn9XFg76HQByQiIh1iqcJkxYoVTJgwgYyMDAzDwDAM5s2b1662CxYsYNiwYbhcLtLS0rjyyispKSnp5ogl1IacMygw18axbHYbp5w7KAwRiYhIR1iqMPnss88oLi4mLa1jcyf86U9/4tprr+Xzzz+nd+/eeL1eFi5cyLnnnsvevXu7KVoJhyvvvIQYR4y/4+sRNrsNZ3wsl//0B2GMTERE2sNShcn1119PVVUVy5Yta3ebhoYG7r7b35Ft0qRJbNmyhXXr1pGUlERZWRkPPfRQd4UrYdBvYB8e/+BXDB45ILDu1FFF/HbFr8nqnxnGyEREpD0sNSonPb3jz+T49NNP2b9/P+AvTAD69OnDyJEjKS4u5u233w7a1u1243a7A8tVVVUAeDwePB5Pk32PLh+/PhpYLbf8M3L4n/fup7aqFgyD+CQX0HL8VsutvaI1L1BuVqXcrKer82rvcSxVmHTGjh07Aq8zM7/9i7lXr14AlJaWBm378MMPM2fOnGbrly9fTnx8fIttiouj93kays16ojUvUG5Wpdysp6vyqq2tbdd+UV+YBGO247kXM2fO5M477wwsV1VVkZ2dzbhx40hOTm6yr8fjobi4mLFjx+JwRM9j4EG5WVG05gXKzaqUm/V0dV5H7zq0JeoLk+zs7MDrsrKyZq9zcnKCtnU6nTidzmbrHQ5H0JPU2jarU27WE615gXKzKuVmPV2VV3uPYanOr+0xZswYioqKmDlzJgAjRowI9E1ZuHAhALt37+bjjz8GYPz48eEJVERERJqxVGHy+uuvU1hYyOjRowPr7rvvPgoLC5kyZQoAJSUlbNiwgT179gAQGxsbGHmzcOFC8vPzGTx4MNXV1fTs2TMwYkdERETCz1KFSVVVFSUlJWzfvj2wrry8nJKSEnbt2hW03bRp03jppZcYOnQou3fvxjAMJk6cyEcffUSfPn1CEbqIiIi0g6X6mEydOpWpU6e2us+2bdtaXD9lypTAVRURERGJTJa6YiIiIiLRTYWJiIiIRAwVJiIiIhIxVJiIiIhIxFBhIiIiIhFDhYmISBQyTQ+m7wCm6Q13KCIdYqnhwiIi0jrTbMCs/i3UzQezDowekDgN4n+MYRjhDk+kTSpMRESiiFn5C6h/G/AdWXEQs/pRDLMWEn8a1thE2kO3ckREooTZuAXqlxIoSo7ddviPmL6a0Acl0kG6YiIiEi08X7WysR68JWA7PWThdDWPz0NNQw1e04vdsJMQm4DDFn1P8+0M07sPs+bP4P4QDBeG6zKIvwbDiA13aB2mwkREJFrY0lrfbrSxPQJVuivZWrmV7VXb2Ve7j/rGekxMDAziYuLoFd+L3ORc8lLySHGmhDvcsDC9uzArJoGvEvACBmb11+B+H3o8i2HYwx1ih6gwERGJFrHngC0DfBU0vZ1jB8cwjJh+4Yqsw2o8NXy691O+qfiGQ+5DOO1O4h3xJDuTsWHDh4+6xjo2H9rM2oq1pDpTGZI+hBFZI0hwJIQ7/JAyq39/TFECYPr/aVgJ7ncg7qJwhdYpKkxERKKEYcRAj6cxD/wEzErADnjB3gcj9f+FO7x221K5hQ93fsiuw7tIj0unMLWwxRFFCY4EcIFpmhysP8hHuz+itKqUUf1GkZ+SH4bIw8T9Dt8WJceyY7rfx1BhIiIi4WI4ToeMD/wjc3y7IWYAOL+PYVijL8a6inW8U/oObq+b/JR87La2b0MYhkGaK40UZwo7D+9kyZYlXJhzIYPTB4cg4kjQ2jgW641xUWEiIhJlDFs8xE8MdxgdtqVyC++UvoNpmuQm53a4vd1mJzc5lz2H9/BO6Ts4Y5wnx5UT1w+gdgHNr5p4MeLGhSOiE2K9UkpERKJOjaeGD3d+iNvrpndi7xM6Vu/E3ri9bv6585/UeKJ/iLSR8J9gz+LbX+lH/nX+AGK/G66wOk2FiYhYgtfrpba6DtM0wx2KdINP937KrsO76JfYNR10+yX2Y+fhnXy699MuOV4kM+w9MdLfwEi8A2K/A7HfxUj5DUbqbzAM6/2a160cEYlojZ5GXpzzf7z5+7epqawlMy+Dq39zkQqUKFLpruSbim9Ij0tvV5+S9rDb7KTHpfNNxTcM6zWM5NjkLjlupDJsKZA4HYPp4Q7lhFmvlBKRk8rjNz/NKw8voqayFoCDew8B8H+PLQ5jVNKVtlZu5ZD7ED3ierSyl0mafTuDYt+nIPYj4oyqNo+bGpfKIfchthza0nXBSrdTYSIiEWt3yV7eeXFFi1dHXvvtYuoO14UhKulq26u247Q7gz5k0EYj57pe5HsJf2CI813OcC7lB4n/Q65jdavHtRk2nHYn26u2d0fY0k1UmIgc4fH6OFjTwP7Dbg7WNODxNn/eiITW+k82Bd3mrm1g29qdIYxGuoPH52Ff7T7iHfFB9xkU+wFZMRsBsBkmhmFi4GN43CISbftbPX68I559tfvw+DxdGrd0H/UxkZNaVZ3/w+rNL3azt9pDvceLzzSxGQZxDjtZKXHk9YynMCOJlHhrzAMRTZLSElvdnpze+naJfDUNNdQ31pPsDN4HJC/2Uwyj6VUzwwCfCbmOz1jrDj4k1hXjospdRU1DDalxqV0VtnQjFSZyUjrsbuSTkgq+3nmA/sDmfdW44pwkuxzYjnzg1TV42bi3mq92HqJHfCyn9U3h7IJ0Ep36sQmVM8ecRlpWKofKqvD5ml7BGjSigL6FJzasVMLPa3oxMbG1cgE/1qgNssXAabQ+HNiGDRMTr9nSzKgSiXQrR046m8sO88q/SvlgUzl2m/+edk56AhlJThKdMcTHxpDojCEjyUn/ngkM6JWE3WbwwaZyFvyrlM1lh8OcwckjxhHD/a//HFdSHBhgd3w7YuPOP1h/9IGA3bBjYOAj+K3Tg96++Mzm/U8MfBzwtj682IcPAwO7xR5kdzLTn35yUvl6VyVvr92L2+NjQGYSdnzQxvxLNsMgPdFJanwspQdqeeOLXYw/JYtT+56cTzINtSEjB/Ly9qf5x4KV7NteTs6QftRwkKy8XuEOTbpAQmwCcTFx1DXWBX343jr3GEbFP49p+m/hAPhMG3VmMjs8Z7R6/LrGOuJi4kiIPbke7GdlKkzkpLG57DBvr92LaUJezyMfUh3o32q3GeT1TGDXoTreXruXOIedwkz1cQiFhOR4Lp42FgCPx8PSpUvDHJF0FYfNQa/4Xmw+tBlcLe9T5i3ko7rrOM35Nsn2ckzTYHfjYNbUX4yX2FaPX+uppTC1EIdNfcSswnK3chYsWMCwYcNwuVykpaVx5ZVXUlJS0mqbqVOnYhhGs69+/azzCHA5MYfdjby/oQy3x0ff1CCffu3UN9WF2+PjHxvKOOxu7KIIRU5eucm5uL3uVifN29tYRHHNbSyu/iVvVt/HJ3WTqTNbv2rpM324ve5OPXdHwsdSV0z+9Kc/cfPNNwOQl5dHRUUFCxcu5MMPP+TLL78kKyur1fZ9+/ZtUoxkZmZ2a7wSOT4pqWDHgVoGZCZ1yfGy01xsLjvMJyUVjBmiWwoiJyIvJY9UZyoH6w+S5kprZU+DBjP4sOLjHao/RKozlfzUk+BBflHEMoVJQ0MDd999NwCTJk3itddeY/fu3RQVFVFWVsZDDz3EE0880eoxbr75ZmbPnt3u93S73bjd7sByVZV/pkGPx4PH03RM/NHl49dHA6vnVlXn4eudB8hMiPH3KTn29o3P2/TfdooBMhNi+HrnAYZmJ5EcF1mXia1+zlqj3KyptdzibfEMTh3MJ3s+ISUmpUumpff6vBysPcjZvc/GZbi69Xsareetq/Nq73EM0yIPnFi5ciWjRo0CYP78+Vx77bUAjBs3juLiYgYMGMDGjRtbbDt16lReeOEFUlJSqKurIzMzk/POO48HH3yQgoKCoO85e/Zs5syZ02z9/PnziY9vf9UuIiJysqutrWXy5MlUVlaSnBx83hrLXDHZsWNH4PWxt2B69fJfRi8tLW21fWxsLL1796ahoYEtW7bw6quvsnz5cr766iv69u3bYpuZM2dy5513BparqqrIzs5m3Lhxzb6pHo+H4uJixo4di8MRWX89nyir5/bmF7vZvK+anPQWeuX7vCTsX8urs96lsrQC0+ev0w27QWyii4vn/RfxGcHvY5dW1FDYK4nLhvbprvA7xernrDXKzZrak9u2qm0s27YM0zTJSmj91nxr9tbsxTAMLup/Ef2T+3f6OO0Vreetq/M6etehLZYpTIJpzwWfn/3sZzz11FMkJvpHUDzzzDPMmDGDgwcP8txzz3HPPfe02M7pdOJ0OputdzgcQU9Sa9uszoq5ebw+9lZ7cMU5oZXLw4e278dTe9ztubpGvnrlQ86+/bKg7VxxTvZVe8Bmx2GPvL7kVjxn7aXcrKm13AakD6CRRt4pfYftNdvpl9ivQ7d1vD4vOw/vxGl3MiZnDAPSB3RV2O0Sreetq/Jq7zEi75M0iOzs7MDrsrKyZq9zcnKCtj311FMDRQnAlClTAq/butIi1na4vpF6jxdXbBsfbi3Ut6bXx85V61pt5oq1U+fxcrheo3NEusLg9MFcnH8xGa4MtlRu4UDdgTb/APWZPg7UHWBL5RYyXZlcnH8xg9MHhyhi6WqWKUxGjBhBeno6AAsXLgRg9+7dfPzxxwCMHz8egKKiIoqKinjqqacCbe+//37Ky8sDywsWLAi87t+/f3eHLmHkNc0jz75pY8cgPwmO+OZXzJo0M8Bnmnit0VVLxBLyU/K5cuCVnNvnXLyml82HNlNaVcr+uv3UeGqo89RR46lhf91+SqtKKTlUgtf0cm6fc5k0cBL5KRqFY2WWuZUTGxvLQw89xPTp01m4cCH5+flUVFRQXV1Nz549AyN2NmzYAMD+/d8+cfJXv/oVDzzwAPn5+ZimGZj3JCsrKzD8WKKT3TCwGQa+tuqGliZaM6DgouGtNzP9M8PagzyuXUQ6J8GRwOjs0ZyZeSZbK7eyvWo7+2r3UeWuwsTEwCAuJo7C1EJyk3PJT80nOTZ4h0qxDssUJgDTpk0jISGBxx57jHXr1hEXF8fEiRN55JFH6NMneOfDBx98kLfeeouNGzdSVVVFYWEhF154Iffcc4/mMolyiXExxDns1DV4W334Xt73Tmfj0tUYR/qJmF4fvYcPoGjiua0ev67Bi8thJzHOUj9KIpaR4kxhaOZQhmYOxePzUNNQg9f0YjfsJMQmaEbXKGS5T9MpU6Y06SNyvJbuRc6aNYtZs2Z1Z1gSoRx2G1kpcWzcW01GUvDbMufOuprcC4ex/R9r8Hl9ZI8aQs6oU7DFtN43pcbdyMCspIjs+CoSbRw2B6lxqeEOQ7qZ5QoTkY7K6xnPVzsPHelr0vItF8Mw6HdOEf3OKWr3cX2mibvRS15PzWkjItJV9GeeRL3CjCR6xMdysKahS497sKaBHvGxFHbRNPciIqLCRE4CKfEOTuubQvlhN942e8G2T6PPR/lhN6f1TSHFpXvcIiJdRYWJnBTOLkgnJy2e0gO1XXK8HQfqyEmL5+yC9C45noiI+KkwkZNCojOG0YMycTps7DpUd0LH2nWoDqfDxuhBma2O9BERkY5TYSInjcLMRMafkoVhwNb9NR2+rdPo87F1fw2GAeNPyaIwM7HtRiIi0iH6c09OKqf2TSHOYef9DWVsKqsmMyGGFh7t14TPNDlY00D5YTc5afGMHpSpokREpJuoMJGTTmFmIlkpcXxSUsHXOw8A/qcEu+KcuGLtR6aZ90+eVuNuxN3opUd8LBcMyODsgnTdvhER6Ub6hJWTUqIzhjFDejG0XxL/fH8jhb2S2FftoarOE5jvxOWwMzAribye8RRmJmn0jYhICKgwkZNa8pFi47KhfcBm53B9I17TxG4YJMbFaEZXEZEQU2EicoTDbqNHQmy4wxAROanpz0ERERGJGCpMREREJGKoMBEREZGIocJEREREIoYKExEREYkYKkxEREQkYqgwERERkYihwkREREQihgoTERERiRgqTERERCRiqDARERGRiKHCRERERCKGChMRERGJGCpMREREJGKoMBEREZGIYbnCZMGCBQwbNgyXy0VaWhpXXnklJSUlbbZ78sknGTJkCE6nk8zMTH784x+zb9++EEQsIiIi7WWpwuRPf/oT1157LZ9//jm9e/fG6/WycOFCzj33XPbu3Ru03b333sutt97KunXryM3N5fDhwzz33HOMHj2a2traEGYg0jU2f76Vey59mIvjJzOx5038/rY/U33wcLjDEhE5YTHhDqC9GhoauPvuuwGYNGkSr732Grt376aoqIiysjIeeughnnjiiWbt9u3bx6OPPgrAXXfdxWOPPcaaNWsYOnQo69evZ968edx5550tvqfb7cbtdgeWq6qqAPB4PHg8nib7Hl0+fn00UG6RZctX2/nvsb+i0ePFNKC+zs1bz73L2lXr+Z/37ifWGWvJvNpLuVmTcrOers6rvccxTNM0u+Qdu9nKlSsZNWoUAPPnz+faa68FYNy4cRQXFzNgwAA2btzYrN3LL7/MddddB8BHH33EOeecA8DAgQPZtGkTY8eOZfny5S2+5+zZs5kzZ06z9fPnzyc+Pr5L8hIRETkZ1NbWMnnyZCorK0lOTg66n2WumOzYsSPwOjMzM/C6V69eAJSWlna43aZNm4K2A5g5c2aTqylVVVVkZ2czbty4Zt9Uj8dDcXExY8eOxeFwtDMra1BukeXaftOpPVzfbL09xuCCq8/ltv+dZsm82ku5WZNys56uzuvoXYe2WKYwCaazF3za087pdOJ0OputdzgcQU9Sa9usTrlFBrvdjqeu+SVRX4yduLi4JnlYKa+OUm7WpNysp6vyau8xLNP5NTs7O/C6rKys2eucnJwubScSqcbdOBqbvfmPrrfRy/cnjwpDRCIiXccyhcmIESNIT08HYOHChQDs3r2bjz/+GIDx48cDUFRURFFREU899RQAY8aMISYmpkm7NWvWsHnz5ibtRKzi2lkTGTAsHwB7jB17jB2AKb+cxJBzBoUzNBGRE2aZWzmxsbE89NBDTJ8+nYULF5Kfn09FRQXV1dX07NkzMGJnw4YNAOzfvx+ArKwsfv7zn/Pwww/z+OOPs3jxYnbs2IFpmgwYMIDp06eHLSeRzohPcjH3n79m5aJ/8fm7XxGX4OR7145i0IjCcIcmInLCLFOYAEybNo2EhAQee+wx1q1bR1xcHBMnTuSRRx6hT58+Qds9+OCD9OrVi3nz5lFSUkJKSgpXX301jzzyCAkJCSHMQKRrxDhiuODqc7ng6nPDHYqISJeyVGECMGXKFKZMmRJ0e0udWg3D4LbbbuO2227rztBERETkBFmmj4mIiIhEPxUmIiIiEjFUmIiIiEjEsFwfExEROTmYjaWYNX8A9z/BloDhugLir8cwmk98KdFDhYmIiEQcs3ELZsVVYNYCXvCBWf0/4F4BPZ7DMOzhDlG6iW7liIhIxDGrf/dtUfLtWmj4GNzvhyssCQEVJiIiEnka/kHTouSoGEz3ByEORkJJhYmIiESgYA98M8GIvgflybdUmIiISOSJuwRoqR+JFyPuB6GORkJIhYmIiEQcI/GnYO8DGEfWHClSXNeA46xwhSUhoFE5IiIScQx7OqS/CXX/h+leBbZEjLhLwTkawzDaPoBYlgoTERGJSIYtERJuwki4KdyhSAjpVo6IiIhEDBUmIiIiEjFUmIiIiEjEUGEiIiIiEUOFiYiIiEQMFSYiIiISMVSYiIiISMTQPCYSUgfLKnnnxRXs21ZG7pB+fH/yKBJSEsIdloiIRAgVJhIya1Z8w6wJD9FQ34DNbsPb6OWF+1/lsffn0P+U7HCHJyIiEUC3cqJMXU097y9YyZu/f5vNn28NdzgBjZ5GfnXV4zTUN2D6TLweL5hQfbCGh6f8DtM0wx2iiEhEMc06TF9NuMMIOV0xiSKri7/kV1c9Tm1VHYYBpgnnXHoW97x6B7FxsWGN7fP3vqayvKrZep/Xx5Y12yldv4vcwf3CEJmISGQxGzdjVj0EDSsBE9NxFkbyTAzHaeEOLSR0xSRKVO6v4v7L/x911fWAvygB+HjJal6479UwRuZXW1nb6vaaQyffXwUiIsczvXsxK66BhlXAkQ9yz2eYFVMwG7eENbZQUWESJd5/ZSUNbk+zWyKmz2TxM8vx+XxhisxvyLmDgj4RNC7BSd7puSGOSEQk8pi1L4F5GPAes9YHeDBr/hymqELLMoVJdXU1d9xxB/369SM2NpaCggLmzJlDY2Njm20Nw2jx65577glB5KFRsfsANnvLp7Ouuh53XUOII2oqo186P/zPi+DY2uTI6+vvuwpXQlxY4hIRiSgNn+MvRI7nhYZ/hzqasLBEHxOfz8ell17KBx98gMPhID8/n02bNjF79mxKSkr4y1/+0q7jDB06FKfTGVjOzo6ekSAFQ/v7O5S2ICsvk7h4Z4vbQun/m3sTvXIyeP2JJezfeYC+hVlc899XcNFN3wt3aCIikcGWBthpesUEwHZkW/SzRGHyxhtv8MEHHwDw+uuvc8kll/Dkk09y66238uKLL3L77bczbNiwNo+zaNEi+vfv383Rhsd5V3yHPgW92LutHJ+3abU95Z4rg95GCSWbzcZVP/shV/3sh5imGRExiYhEEiN+EqZ7WQtbfBjxV4U8nnCwRGHy1ltvAeByuZgwYQIAkyZN4tZbbwXg7bffbldhctZZZ1FbW0teXh7XX389d9xxR5MrKMdzu9243e7AclWVf1SJx+PB4/E02ffo8vHrQ8aAR9+9l9//9M+sfvcrMCElPZFrZk5kzHWjTiiusOfWjaI1t2jNC5SbVSm39jGNczGdM6DmOfy9LQzAC66JGPYJGCH8/nX1OWvvcQzTAhNIjB8/nmXLltGvXz927NgBgNfrJSbGX1dNnz6defPmBW1vGAY9evQgJyeH7du3c+jQIQCuvvpqXn01+IiV2bNnM2fOnGbr58+fT3x8/AlkJCIicnKpra1l8uTJVFZWkpycHHS/sF4xufvuu3n00Udb3WfdunUtru9IPfXxxx/zne98B8MwqK2t5dJLL+W9997jr3/9K4899ljQviYzZ87kzjvvDCxXVVWRnZ3NuHHjmn1TPR4PxcXFjB07FofD0e7YrEC5WU+05gXKzaqUm/V0dV5H7zq0JayFyfDhw7nxxhtb3SclJSVQOOzfvx+fz4fNZqOsrCywT05OTqvHOPvsswOv4+PjueKKK3jvvfcA2LFjR9DCxOl0tnirx+FwBD1JrW2zOuVmPdGaFyg3q1Ju1tNVebX3GGEtTK666iquuqrtzjzjx4/n2Wefpb6+nqVLl3LJJZewcOHCJtvB37l15syZALz77rv07duXFStWUFZWxhVXXIHdbqe+vp4333wz0DY3V/NniIiIRApLdH69/PLLGTVqFP/85z+ZOHEiBQUFbNy4EYDJkycHOr5WVlayYcMG4NtONlu2bOGmm24iISGB/Px8du7cycGDBwG46aab6Nu3bxgyEhERkZZYYoI1u93OkiVLuPXWW8nIyKCkpIScnBzuu+8+nn/++Vbbjho1ihkzZpCTk8PWrVvx+XwMHz6cefPm8Yc//CE0CYiIiEi7WOKKCUBycjK/+93v+N3vfhd0n6lTpzJ16tQm6woLC3n66ae7OToRERHpCpa4YiIiIiInBxUmIiIiEjFUmIiIiEjEsEwfk0hwdFK3liaJ8Xg81NbWUlVVFXXj2JWb9URrXqDcrEq5WU9X53X0d2dbE6SqMOmA6upqILqeSiwiIhJK1dXVpKSkBN1uiWflRAqfz8fu3btJSkpq9mTco9PV79ixo9VnAFiRcrOeaM0LlJtVKTfr6eq8TNOkurqaPn36YLMF70miKyYdYLPZ6NevX6v7JCcnR9V/zGMpN+uJ1rxAuVmVcrOersyrtSslR6nzq4iIiEQMFSYiIiISMVSYdBGn08n999/f4tOIrU65WU+05gXKzaqUm/WEKy91fhUREZGIoSsmIiIiEjFUmIiIiEjEUGEiIiIiEUOFiYiIiEQMFSadVF1dzR133EG/fv2IjY2loKCAOXPm0NjY2GZbwzBa/LrnnntCEHlTCxYsYNiwYbhcLtLS0rjyyispKSlps92TTz7JkCFDcDqdZGZm8uMf/5h9+/aFIOL260xuU6dObfHctDWxXiitWLGCCRMmkJGREYhv3rx57Wrb2fMdCp3Na/bs2UF/ptrz8xgKjz/+OKNHj6Z37944nU5yc3O58cYb2bJlS5ttI/mcQedzs8J5mzt3LmeccQapqak4nU769evHVVddxZo1a9psG+mfkZ3NLSSfkaZ0mNfrNS+44AITMB0Ohzlo0CDTZrOZgHn99de32R4wAXPo0KHm2WefHfiaN29eCKL/1rPPPhuIJS8vz0xOTjYBMzMz09yzZ0/Qdvfcc0+g3YABA0yXy2UCZlFRkVlTUxPCDILrbG433nijCZh9+/Ztcm4uvfTSEEbfut/+9rdmTEyMOXDgwECOTz/9dJvtOvs9CZXO5nX//febgNmzZ88m5+zss882GxsbQxB523Jzc03DMMyioiIzLy8vkF9WVpZZWVkZtF2knzPT7HxuVjhvl19+udm7d2/zzDPPNAcPHhz4nE9LSzMPHz4ctJ0VPiM7m1soPiNVmHTCwoULA//pFi9ebJqmaT7xxBOBdatXr261/dH9tm7dGoJoW+Z2u82ePXuagDlp0iTTNE1z165dZlJSkgmYP/3pT1tst3fvXtPhcJiAedddd5mmaZpffvmlaRiGCZiPP/54yHIIprO5mea3P3T3339/iKLtuP3795u1tbXm1q1b2/0L/ES+J6HSmbxM89tfcDfeeGP3B9lJDzzwgLl9+/bA8u233x7I8fXXX2+xjRXOmWl2LjfTtMZ5q6ura7J8bMHx73//u8U2VviMNM3O5WaaofmM1K2cTnjrrbcAcLlcTJgwAYBJkyYFtr/99tvtOs5ZZ51FfHw8p5xyCo888ghut7vrgw3i008/Zf/+/cC3sffp04eRI0cCwXN455138Hg8TdqdfvrpFBYWttoulDqb27Hmzp2L0+kkOzuba665JqIunaenp+NyuTrUpiu+J92tM3kda+HChbhcLnr37s0ll1zC559/3oXRnZhf/vKX5OTkBJbPP//8wOtgk1dZ4ZxB53I7ViSft7i4OBYtWsTIkSMZMmQIDz30EAAZGRkMHDiwxTZW+IyEzuV2rO78jFRh0gk7duwA/B+kR5+Q2KtXr8D20tLSNo/Ro0cP+vXrh9Pp5JtvvmHmzJnccMMN3RNwC47mAJCZmRl4fTSPYDl0tl0onWiMsbGx9O7dm379+rFz505effVVRowYwa5du7on4BCwwnk7EXa7naysLPr378/evXtZsmQJ55xzTkT9kjvK6/Xyhz/8AYD8/HzGjBnT4n5WPGftze0oK5y3ffv28cknn7Bu3Tp8Ph95eXm8//77JCUltbi/lc5bR3M7qrs/I1WYHOPuu+8O2hnr6Nf69etbbGt2YALdjz/+mIqKCr744gt27drF97//fQD++te/NvlPHQ4dyaMr2oVSe2L82c9+RkVFBevWraOkpCTQ+fLgwYM899xz3R1iyFnhvLVl8uTJlJWVsWnTJtatWxf4i9TtdvP73/8+zNE1VVNTwxVXXMGyZcvIyspi8eLFHZ7uO1LPWUdzs8p5mzFjBj6fj+3bt/OjH/2IrVu38qMf/Yjq6uoOHScSz1tncgvFZ6QKk2MMHz6cG2+8sdWvlJQUsrOzAdi/fz8+nw+AsrKywHGOvazZkrPPPhvDMACIj4/niiuuCGwLVWFyNAdoGvvR18Fy6Gy7UDqRGE899VQSExMDy1OmTAm8jqS/dDrKCuetswYOHEhaWlpg+aKLLiI9PR2IrHO2d+9eLrjgAhYvXszAgQNZuXIlQ4YMCbq/lc5ZR3MD65w38I+kzMnJYdasWQCsXbuWV155pcV9rXTeoGO5QWg+I1WYHOOqq67i+eefb/Wrd+/ejB8/HoD6+nqWLl0K+O+THnV0+6JFiygqKqKoqChwiWvFihW89tpreL3ewDHefPPNQNvc3NyQ5DpixIjAh8DR2Hfv3s3HH3/cJIej8T/11FMAjBkzhpiYmCbt1qxZw+bNm5u0C6fO5gZw//33U15eHlhesGBB4HX//v27O/QuM2bMGIqKipg5cybQ/u9JpDs+L4BHH320yQdicXExFRUVQOScs7Vr1zJy5EhWr17N+eefz6pVq8jPz2+yj1XPWWdyg8g/bxUVFbz44os0NDQE1h39vAf/FSKw5mdkZ3ODEH1Gdlu32ijW2Nhojho1KjBcuKioKDDUavLkyYH9nnvuuWYjcI6uS0hIME877TSzR48egX1uuummkObxzDPPtDgUsWfPnuauXbtM0/x2BNGxPbBnzpwZWD9w4MDAULgBAwa0OswslDqbG2DabDazsLDQLCgoaDL0cd++fWHKpqmFCxeaBQUFZm5ubiC+jIwMs6CgIPD/7+i2Y0c8tOd7Ek6dzevocNWcnBxz8ODBgdEPCQkJ5tq1a8OUTVPHDoE+fpqAP/7xj6ZpWvOcmWbnc4v083Z0dJjL5TJPPfVUMzs7O5BnUlKSuW3bNtM0rfkZeSK5heIzUoVJJ1VWVpq33nqr2adPH9PhcJj9+/c377vvPrOhoSGwT0uFyaZNm8wZM2aYgwcPNhMTE82UlBRz+PDh5rx580yPxxPyPF566SVz6NChptPpNFNSUsyJEyeaGzduDGxv6T+mz+cz586daxYVFZkOh8Ps2bOneeONN0bMvApHdSa3Bx980Bw1apSZmZlpxsXFmYWFheaMGTPMnTt3hiGDlh37/+r4rwsuuMA0zZZ/EZhm29+TcOpsXs8884x54YUXmr179zadTqfZv39/c8qUKeb69evDk0gLji22jv86+v/PiufMNDufW6Sft4MHD5rXXHONmZ+fb7pcLjMmJsbMzs42r7vuOvObb74J7GfFz8gTyS0Un5HGkTcXERERCTv1MREREZGIocJEREREIoYKExEREYkYKkxEREQkYqgwERERkYihwkREREQihgoTERERiRgqTERERCRiqDARERGRiKHCRESixooVK5gwYQIZGRkYhoFhGIHHsouINagwEZGo8dlnn1FcXExaWlq4QxGRTlJhIiKW8Oabb2IYBjabjX/84x8A/P3vfw+se/fdd7n++uupqqpi2bJl4Q1WRDpNhYmIWMJll13GzTffjGma3HLLLezZs4fp06cDcPvttzNmzBjS09NxuVxhjlREToQKExGxjN/+9rcUFBSwefNmzjzzTHbv3s1pp53Gww8/HO7QRKSLqDAREctITEzkpZdewmazsW/fPhwOBy+//DJOpzPcoYlIF1FhIiKWUlpais/nA8Dj8bB9+/YwRyQiXUmFiYhYxq5du5gxYwYAQ4cOBeDmm2+mvLw8jFGJSFdSYSIilmCaJlOnTuXgwYOce+65rFq1itNPP519+/Yxbdo0AF5//XUKCwsZPXp0oN19991HYWEhU6ZMCVPkItIRMeEOQESkPZ544gneeecdXC4Xzz33HHFxcbzwwgt85zvf4Y033uDPf/4zNpuNkpKSJu3Ky8spLy+nX79+YYpcRDrCME3TDHcQIiIiIqBbOSIiIhJBVJiIiIhIxFBhIiIiIhFDhYmIiIhEDBUmIiIiEjFUmIiIiEjEUGEiIiIiEUOFiYiIiEQMFSYiIiISMVSYiIiISMRQYSIiIiIR4/8H/jyg5GVthVsAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ],
      "source": [
        "plt.figure(figsize=(6,4))\n",
        "\n",
        "all= [int(sum([value(k*instance.U[i,k]) for k in instance.k])) for i in instance.i]\n",
        "for i in instance.i:\n",
        "  print(\"TAG300\", i, int(sum([value(k*instance.U[i,k]) for k in instance.k] ) ) )\n",
        "\n",
        "plt.scatter(df[\"x\"],df[\"y\"], c=all, s=  20)\n",
        "\n",
        "for k in instance.k:\n",
        "  node = int(sum([value(i*instance.mu[i,k]) for i in instance.i]))\n",
        "  xk,yk = df.loc[node-1,'x'], df.loc[node-1,'y']\n",
        "  plt.scatter(xk,yk, s=  300 , alpha =0.4, zorder=-2)\n",
        "  print(k, xk,yk)\n",
        "plt.xlabel(\"x1\", fontweight='bold')\n",
        "plt.ylabel(\"x2\", fontweight='bold')\n",
        "plt.xticks(fontweight='bold')\n",
        "plt.yticks(fontweight='bold')\n",
        "plt.grid()\n",
        "plt.show()"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.8.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}