blue-green-points-biggest_circle.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyomo.environ import *\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import random \n",
    "import matplotlib.pyplot as plt\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<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",
       "      <th>c</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.562539</td>\n",
       "      <td>0.456148</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.305930</td>\n",
       "      <td>0.990781</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.480551</td>\n",
       "      <td>0.974520</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.519639</td>\n",
       "      <td>0.173766</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.528628</td>\n",
       "      <td>0.769045</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.243803</td>\n",
       "      <td>0.832328</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.379342</td>\n",
       "      <td>0.371002</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.360388</td>\n",
       "      <td>0.157756</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.793967</td>\n",
       "      <td>0.582600</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.915850</td>\n",
       "      <td>0.602913</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0.241242</td>\n",
       "      <td>0.904694</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>0.769900</td>\n",
       "      <td>0.905204</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>0.029838</td>\n",
       "      <td>0.693053</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.407071</td>\n",
       "      <td>0.066217</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>0.513599</td>\n",
       "      <td>0.421776</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>0.764609</td>\n",
       "      <td>0.639988</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>0.197498</td>\n",
       "      <td>0.254418</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>0.340360</td>\n",
       "      <td>0.190976</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>0.198594</td>\n",
       "      <td>0.404394</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>0.285602</td>\n",
       "      <td>0.811163</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>0.649476</td>\n",
       "      <td>0.991961</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>0.036390</td>\n",
       "      <td>0.274849</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>0.591636</td>\n",
       "      <td>0.235543</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>0.450146</td>\n",
       "      <td>0.672424</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>0.088437</td>\n",
       "      <td>0.430823</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           x         y  c\n",
       "0   0.562539  0.456148  1\n",
       "1   0.305930  0.990781  1\n",
       "2   0.480551  0.974520  1\n",
       "3   0.519639  0.173766  1\n",
       "4   0.528628  0.769045  1\n",
       "5   0.243803  0.832328  1\n",
       "6   0.379342  0.371002  0\n",
       "7   0.360388  0.157756  0\n",
       "8   0.793967  0.582600  0\n",
       "9   0.915850  0.602913  0\n",
       "10  0.241242  0.904694  0\n",
       "11  0.769900  0.905204  0\n",
       "12  0.029838  0.693053  1\n",
       "13  0.407071  0.066217  1\n",
       "14  0.513599  0.421776  0\n",
       "15  0.764609  0.639988  0\n",
       "16  0.197498  0.254418  1\n",
       "17  0.340360  0.190976  1\n",
       "18  0.198594  0.404394  0\n",
       "19  0.285602  0.811163  0\n",
       "20  0.649476  0.991961  1\n",
       "21  0.036390  0.274849  0\n",
       "22  0.591636  0.235543  1\n",
       "23  0.450146  0.672424  1\n",
       "24  0.088437  0.430823  0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N= 25\n",
    "df=pd.DataFrame()\n",
    "df['x']= [random.random() for i in range(N)]\n",
    "df['y']= [random.random() for i in range(N)]\n",
    "df['c']= [random.randint(0,1) for i in range(N)]\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = AbstractModel()\n",
    "model.i = RangeSet(N)\n",
    "model.x = Var(model.i, within=Binary)\n",
    "model.R = Var(bounds=(0,0.5), within=Reals)\n",
    "model.Count = Var(bounds=(1,N), within=Reals)\n",
    "\n",
    "x0, y0 = 0.5, 0.5 \n",
    "\n",
    "def D_rule(model, i):\n",
    "    return np.sqrt( (df.loc[i-1,'x']-x0)**2 + (df.loc[i-1,'y']-y0)**2 )\n",
    "model.D = Param(model.i, rule=D_rule)\n",
    "\n",
    "def C1_rule(model, i):\n",
    "    return model.D[i] <= model.R + 2*(1-model.x[i])-0.01\n",
    "model.C1 = Constraint(model.i, rule=C1_rule)\n",
    "\n",
    "def C1A_rule(model, i):\n",
    "    return model.D[i] >= 0.01+ model.R - 2*(model.x[i])\n",
    "model.C1A = Constraint(model.i, rule=C1A_rule)\n",
    "\n",
    "\n",
    "\n",
    "def C2_rule(model):\n",
    "    return sum(model.x[i]*(1-df.loc[i-1,'c']) for i in model.i ) == model.Count\n",
    "model.C2 = Constraint( rule=C2_rule)\n",
    "\n",
    "def C3_rule(model):\n",
    "    return sum(model.x[i]*(df.loc[i-1,'c'])  for i in model.i) == model.Count\n",
    "model.C3 = Constraint( rule=C3_rule)\n",
    "\n",
    "def obj_rule(model):\n",
    "    return model.R\n",
    "model.OF = Objective(rule=obj_rule,sense=maximize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time is  0.0494\n",
      "feasible\n",
      "0.464918075924774 9.0\n"
     ]
    }
   ],
   "source": [
    "opt = SolverFactory('glpk')\n",
    "instance = model.create_instance()\n",
    "\n",
    "start = time.time()\n",
    "results = opt.solve(instance) # solves and updates instance\n",
    "print('Time is ', round(time.time() - start, 4) ) \n",
    "\n",
    "if (results.solver.status == SolverStatus.ok) and (results.solver.termination_condition == TerminationCondition.optimal):\n",
    "    print('feasible')\n",
    "elif (results.solver.termination_condition == TerminationCondition.infeasible):\n",
    "    print('infeasible')\n",
    "else:\n",
    "    print ('Solver Status:',  results.solver.status)\n",
    "print(value(instance.OF) , value(instance.Count))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcEAAAGtCAYAAABujc/JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABFCUlEQVR4nO3de3Bb2Z0f+O8BefkESYkPkSIJSqTej252t2TSbvmhtqsm3V2puGbGu2sn8VRmx2mqJs3KVmWL49rsJpWazdaGqa2dRPFMu+24piZRokw8romz0YwrmYzsHtlNquUm9eqWmk09QEqUKJEiCVIkQOLsHz9e3QsQFAESuPcC9/upQvFeACSOjoD7wznnd85RWmsQERH5UcDtAhAREbmFQZCIiHyLQZCIiHyLQZCIiHyLQZCIiHyr2O0CbFV9fb3evXt3ysfm5+dRWVnpbIHyGOsrM6yv9LGuMsP6ykyq+rp06dIjrXXDRr+b90Fw9+7d+OCDD1I+dv78eZw8edLZAuUx1ldmWF/pY11lhvWVmVT1pZS6k87vsjuUiIh8i0GQiIh8i0GQiIh8i0GQiIh8i0GQiIh8i0GQiIh8i0GQiIh8i0GQiIh8i0GQiIh8i0GQiIh8y7EgqJT6gVLqoVLq6jqPK6XUv1RKjSilLiulXnGqbERE5E9OtgT/EMDrz3n8DQD7Vm9vAfgDB8pEREQ+5lgQ1Fr/DMDUc57yVQB/pMX7ALYppXY6UzoiIvIjL+0i0QIgbDsfW73vfvITlVJvQVqLaGxsxPnz51P+wUgksu5jtJZX6yu2EsPE/ATmo/OoLKlEU2UTjCLD7WJ5tr68KBt15dX3QS7wvZWZrdSXl4KgSnGfTvVErfW7AN4FgOPHj+v1thzhdiSZ8WJ9hWfC6HynE5FoBLF4DEbAQLAkiOFTwwjVhFwtmxfry6u2Wldefh/kAt9bmdlKfXkpO3QMgP3d3ArgnktlIY/ov9D/7MIHALF4DJFoBP0X+l0uGTmJ7wPKFS8FwR8D+I3VLNHPApjRWq/pCiV/GRgfeHbhM8XiMQyOD7pUInID3weUK05Okfj3AH4B4IBSakwp9VtKqVNKqVOrTzkHYBTACIDvAfhtp8pG3tXd0g0jkDjuYwQMdLV0uVQicgPfB5Qrjo0Jaq2/scHjGsDfc6g4lCf6TvThzJUza8aC+k70uV00chDfB5QrXuoOJVojVBPC8Klh9BzrQVdzF3qO9RRsMsRGwjNh9J7rRdf3utB7rhfhmfDGv1Qg+D6gXPFSdihRSqGaEE6/edrtYrgqOTtyaGIIZ66c8VUg4PuAcoEtQaI8wOxIotxgECTKA8yOJMoNBkGiPJAqO7IExfjcjmPA06dyW1wElpaAWAxYXgZWVoB4HNAp15wgInBMkMhdWkvgWlxcezMD2fIy/g98CXfv/iGi0Tj0ygrKUITyolL8wwPHgf/239J7LaXkVlwMlJRYN8NIPLffSkvlcaICxSBIlEtaAwsLwNwcMD9vtdgWF+V4aSmtltoOlOL7X/kX+NFHP8LNxzexv24/fu3Qr6GhsiGzsmgNRKNyS5dhAJWVQEVF4s/KSgmSKtWKh0T5gUGQKBvswc5+i0SkSzILGiob0HO8x7ojEACKiuRmMrs/7betvn4sBjx5IrdkgUBiUKypAbZtk2MGR8oDDIJEmdIamJlB2b17wIcfWsFuZWVzf6+kBCgrs27l5fKzuFgCnP2neVxUJAEokzKbt+VlqzUYjUqQM8cS7fdHo9Jifd6/Kx6Xf3skknh/UZEVEGtqULSwIK/NwEgewyBItBGzJTQ1JbfpaWBlBZW3bwNVVRv/flmZPC8YtAKc+bOsLLNgtlnmeCAgAaq0NP3fXVyUVu78vNzM44WF9btVV1as+gKw7coVub+mRm7btwN1dfLvJ3IRgyBRsqdPrQv41BQwO5ve75WWSrBLvuV7YokZrGtr1z4Wi1lBcW4OmJmRLwxLS2ufaw+Mt27JfcEg0NAA1NdLUMz3uqK8wyBIFI/LhfnBA7nNz2/8O+XlWKqvB44eBaqrJdiVlOS+rF5jGFbrzm5xUYLhzAwwM4P4xx+n/n2zK/XWLWmpbtsmAbG+XoKuE61k8jUGQfKnpSXg4UMJepOTMk62HqUk0NXWWreyMkTOnwfa2x0rcl4pKwOamuQGYHphAfjsZ62WotkitCftaC1dzdPTwCefSACsq5OW4s6dkoBDlGUMguQfMzNWay9VpqOpqMgKdtu3y62YH5UtM7tVGxvl3OweffRIvojMzCQ+Px6X+ycngevXpbXZ3CwBsbLS+fJTQeInmwrb7CwwNgaMj0sX3XoqKuTi3NgorQ92w+VeUZG08hoagEOHZHzx0SMrKCZ3S692reKjj6Rl3twsNwZE2gIGQSo8i4sS9MbG1k9qUUpaejt2SOBLJ8uTcsswpJW3c6ecP30qwXBiQn7au05nZ+X28ccSEHfulIAYDLpTdspbDIJUGJaX5WI5NiYtiVSrsBiGFfR27GAmoteVlwNtbXKLxaQb+/59GctNFRBv3JAu0127gJYWdmFTWvguofyltQS8sTG5OKaa1F1UJMkZra3S7cbJ2vnJMOT/sLVVvvA8eADcu7c2IM7MAJcvA9euSSDcvXtt5iqRDYMg5Z9YDLh7V9Lqnz5N/Zz6erlg7tzJFkGhKS6WANfSIgHx4UMJiA8eWAFxZUXeI3fvsnVIz8V3BOWPSAQYHZWWX6pWX1WVBL6WFulKo8JXXGwlyMRiQDgsgW9uznpOcutw1y6Zj0gEBkHyOq0lKWJ0VH4mKymxusnY7eVvhgF0dMhtagq4c0daiOu1DvfulZ4CdpH7GoMgedPysrT4bt1auzgzIBmBHR3yzZ7TGSiZOc/z6FF5H925s7Z1eOmSTK/Yu1e+RPF95EsMguQtS0vAp5/Kt/VYbO3jTU0S/OrqnC8b5R/DkFV92ttlJZo7d2T6jNk6nJ8Hhocls3TPHslE5bihr/B/m7whGpXgd+vW2vG+4mK5OLW3c+ks2jxz9Z/Dh6V7/fZt64vW4qKMGd68aQVNP64F60MMguSuWEwuSKOja9fvrKyUi1EoxG/nlD0lJcDBg9INevu2vPfMXS9iMQmEn34qCTQdHUyyKnC8stCWhMNAfz8wMAB0dwN9fRKzNrS8bAW/5G7P6mrgwAGZ1M6kBcqV4mIJhB0d8kYeGZFtoQDpjTBbi+3twL59m15cYdOfEXIEgyBtWjgMdHZK3kosBgwNAWfOyBDLuh/y5WW5sIyMrA1+VVUS/JqaGPzIOYGAtPra2iSbdGTEWm4vHrfGqPfvl8n3GSTQbOozQo5iOhRtWn+/9eEG5GckIvevEY/LN+u/+AtZANkeACsrgVdeAb70Jaask3uUkmzjL31Jmmz2uYSxmIwZ/uVfSmJNmjL6jJAr2BKkTRsYWNuYi8WAwcGkJ05OAlevrp3qUFEh365bWxn4yFt27JDbvXvypc3sJl1YAH75S/lCd/jwhlnKaX9GyDUMgrRp3d3SvWP/kBsG0NW1erKwIPvA3b+f+Ivl5Vbw49ws8rLmZumev31bEmbMN/uTJ8DPfy7j1ocPr7t7xYafEXIdr0C0aX198tk38wUMQ877/sGKXDD+8i8TA2BxsVwwvvxlGX9hAKR8EAhI8sxXviKJNPb37YMHwPnz0tORnN2M53xG+pwpOm2MLUHatFBIBvj7+6V7p6sL+PbfmUDLp9es7iNTa6tsnFpW5k5hibbKMOQ9vHu37GM4Nib3ay3zW+/fB154QVqOq1J9Rpgd6i0MgrQloRBw+jRkvO/qVeBe0vqeNTWydFVtrSvlI8q68nLg5ZeldXj9umznBciE+4sXJQi+8MKzL3zPPiPkSQyCtDXxOPDJJ3Kzb2Rrfmtua2PSCxWmmhrgc5+T5JmrV60J9xMTEhgPHZKpF3z/exqDIG3e7Czw4YfWnCpAPvC7dsl8Py47RX7Q3CwbNl+/LvMJARkfvHJFukw7O2UOLHkSgyBlTmuZUHzzZuKu3rW10g1UXe1e2YjcYBgS7FpbZRBwfl7un54GfvYzWZx7/34mg3kQgyBlJhKRnO/paeu+QEC6ftrb2fVD/lZXB5w8KcMDIyPyJdEcMrh3TxaF4Ia+nsIgSOkxM+A++iix9bdtmyQJrDNPish3AgEZDmhulh3tp6bk/vl54MIFWby7o4NfGD2CQZA2trAgrb/Hj637AgHp3tm7lx9molSqqoBXX5VxwuvXZZwwHpfjyUn58lha6nYpfY9BkJ4vHJYBfvsef9XV8gHm2B/R85mJYvX1stzakydy/+Qk8NOfAi+9JMuzkWs4SkuprazIAP/QkBUAlZItZb7wBQZAokxUVgInTkjPiWlpSRYXvXYtcYiBHMWWIK21sAB88AEwM2PdFwxK64+D+kSbYyaQ1dfL1CJzXuHoqAw1vPIKx9ZdwJYgJXrwQFK67QGwtRX44hcZAImyoaFBMkgbG637ZmbkcxcOu1Ysv2JLkITWqLh7F5ibs+4LBIAjR2StRCLKnpISWUj01i1JlInHZdhhaEjGDdk96hi2BEm6Zd5/H+XmgsCArI944gQDIFEutbfLGLu9G/T2bVQnbzxNOcMg6HfmihbmIsCAdNew+5PIGdXV8nlrbn52lzEzA7z33tqNqCnrGAT9LByWjUEXF6379u+XnUC57ieRc4qKgGPHZCK9aX5eAuHDh+6VywcYBP3qxg0ZfzDHHgwDs4cOyUoXnPxO5I59+4Djx6HNNUaXl2UaxaefuluuAsYg6DfxuAS/mzet+1a7Y2Lbt7tWLCJatXMnZl54QcblTdevy7QKJsxkHYOgnywvy/bW9jTshgZJgKmocK9cRJRgpbJSEmbsm1GPjcnwhTm/kLKCQdAvFhdl8d5J287voZCkaRdzpgyR55SWyqa9bW3WfdPT8jleWHCvXAWGQdAP5uaAv/qrxM1vDxyQdQuzsL9ZeCaM3nO96PpeF3rP9SI8wwm/RFkRCMg+hUeOWGP15m4UzBzNCjYBCt3jx8DFi9acI6WAF19M/Ha5BeGZMDrf6UQkGkEsHsPQxBDOXDmD4VPDCNWEsvIaRL7X0SFDFpcuybig2bPz2c8CNTVuly6vsSVYyO7dA95/3wqAxcUy/SFLARAA+i/0PwuAABCLxxCJRtB/oT9rr0FEAJqa5PNbVCTn0aiMEdq3OKOMMQgWqvFx2brFzCYrK5O9zRoasvoyA+MDzwKgKRaPYXB8MKuvQ0SQxbdffRUwDDlfXpYvupxLuGkMgoVofFzSqbWW82AQ+Pznc9Jt0t3SDSNgJNxnBAx0tXRl/bWICLKS04kT8sUWkC+6g4PS80MZYxAsNMkB0Nzd2j7nKIv6TvQhWBJ8FgiNgIFgSRB9J/py8npEBOtzbU5t0lrGC+/ccbdceYhBsJCkCoCf+5ykWudIqCaE4VPD6DnWg67mLvQc62FSDJETzI16q6qs+y5flv0JKW3MDi0ULgRAU6gmhNNvns756xBREnOsf2BAtmACZKf64uKsJsAVMrYEC4GLATAbwmGgt1fm7ff2cl9RooyUlMjn3b66zOXLHCNME1uC+a4AAmBnp8z7jcVkWdMzZ4DhYVnQhojSUFws3yJ/8QvZpV5ruS4UFwM7drhdOk9jSzCf3b+fOgkmTwIgAPT3WwEQkJ+RiNxPRBkwDJk8b27QG48DH3wATE25Wy6PYxDMV1NTMg8wOQDm2T6AAwNrN9COxSTjm4gyZHaNmtngKyvyIZuZcbdcHuZoEFRKva6UuqGUGlFKfTvF4zVKqf+slBpWSl1TSv2mk+XLG5GIRAlzInwwmJcBEJAFMIzEaYYwDOnZIaJNKCtLHBIxJ9RzrdGUHAuCSqkiAN8B8AaAwwC+oZQ6nPS0vwfguta6E8BJAP+PUir/ruy5tLSUuBRaaWle7wTf1ycx3AyEhiHnfZxmSLR5lZXSNWp+sKJRGS/k7hNrONkS7AIworUe1VpHAZwF8NWk52gAVUopBSAIYArAsoNl9DZzl+mnT+W8qEiaTHm8F2AoJEkwPT3yT+npYVIMUVZUVyeuNbq4KIEwGnW3XB6jtDmmlOsXUuprAF7XWn9r9fybALq11m/bnlMF4McADgKoAvA/aa3/S4q/9RaAtwCgsbHx2NmzZ1O+ZiQSQdAcJM538TiqPv4YJeZcIKUwe/BgVneDL6j6cgDrK32sq8xks76MmRlUXb8OtXqtj1VXY/bw4axso+YVqerrtddeu6S1Pr7R7zo5RUKluC85Av81AEMAvgxgD4D/qpR6T2s9m/BLWr8L4F0AOH78uD558mTKFzx//jzWeyzvmM0js4nU2Zn1ybAFVV8OYH2lj3WVmazXV3e3bKlmqq+XLdUKxFbqy8mvAmMA7J1crQCSZ3P+JoAfaTEC4BakVehvN28Cd+9a5/v2cTUIIkpfUxNw6JB1fucOcOuWe+XxECeD4EUA+5RS7avJLl+HdH3a3QXwFQBQSjUCOADA3wvhjY0BN25Y562twEF+LyCiDO3dC7S0WOfXrgGTk+6VxyMcC4Ja62UAbwP4CYCPAPyx1vqaUuqUUurU6tN+F8CrSqkrAP4CwO9orR85VUbPmZmRblBTQ4N0gxIRbUZnp2zFBFg7T8zPu1oktzm6bJrW+hyAc0n3vWM7vgfgV5wsk2fFYrLagzkXsKoKOH68oAazichhRUXAZz4DvPeeZIuaK1N8/vNrJ+z6BK+oXmSu+2fO6SkulgBYzKVeiWiLysokEJpfqCORxNWnfIZB0ItGRoAHD6zzl16y1gMkItqqbduAl1+2zh8+BK5fd604bmIQ9JpHjxITYfbsAXbudK88VHDCM2H0nutF1/e60HuuF+EZ7l3lS83NkmluGh0FJibcK49L2L/mJYuLMlBtdkvU1jITlLIqPBNG5zudiEQjiMVjGJoYwpkrZzB8ahihGi7T4zsHDgBzc1bwGxoCTp6ULlOfYEvQK8xtT8wljUpLgWPHmAhDWdV/of9ZAASAWDyGSDSC/gvcu8qXlJLhFnPXiVjMd+ODvMJ6xfXrwPS0HCslAdBH38bIGQPjA88CoCkWj2FwnHtX+ZZhAK+8ItcdAHj8WPISfIJB0AsmJhJXbzh0CKirc688VLC6W7phBBJT4Y2Aga4W7l3la7W1ieODN25YX8oLHIOg26LRxAnxO3dKMgxRDvSd6EOwJPgsEBoBA8GSIPpOcO8q39u/X4IhIN2hv/zl2h2vCxCDoNuuXLHGAcvKuCIM5VSoJoThU8PoOdaDruYu9BzrYVIMCaWkW9ScNL+wINenAsfsUDfduyc3U2enb1dtIOeEakI4/eZpt4tBXlReLrtLXLok5+PjslxjAW/wyZagW5aWgMuXrfO2NmDHDvfKQ0QEyPxB+y41V64U9PqiDIJuuXzZ6m8vLweOHHG3PEREpqNHrVWqVlYS8xYKDIOgG8bGEldmeOklrgtKRN5RVLR22kS4MFcWYhB02uIicPWqdb57t+zyTETkJTU1QEeHdX7tmpXEV0AYBJ02PGx1g1ZUJO72TETkJQcOJK4mc+2au+XJAQZBJ42NyWrtppdfZjcoEXlXUZFki5rGxgpuN3oGQacsLyduVdLRYU1MJSLyqh07JGPUdOWKJMsUCAZBp9y8KdMiAJkUf+CAu+UhIkrX0aPWHOb5eeCTT9wtTxYxCDphfj5xbdDDh9kNSkT5o7Q0MX9hZASYnXWvPFnEIOiEq1dlqyRAukBbWtwtDxFRptraEtcWvXy5ILZcYhDMtQcPEpNhjh51ryxERJullCTJmHucTk/Lsmp5jkEwl+LxxJTiXbtk7g0RUT6qqkrc5ebjj61erjzFIJhLo6PWmnuGARw86G55iIi2au9eoKREjp8+BW7fdrU4W8UgmCuLi4kZVAcOWG8cIqJ8VVwsew+abt7M630HGQRz5eOPZW4gIF0Iu3a5Wx4iomzZtUtWvAIkAI6MuFueLWAQzIVIJHGx2aNHrcFkIqJ8FwgkDu+MjkrvVx7ilTkXbt60jnfs4ALZRFR4mputRL94HLhxw93ybBKDYLZFIolpw1wZhogKkVKy8IcpHAbm5twrzyYxCGZbcitw2zbXikJElFP19UBDgxxrLbkQeYZBMJvYCiQiv7G3BicmgKkp98qyCb4PguGZMHrP9aLre13oPdeL8MwWdk9mK5CI/Ka6Gmhttc7zLFPU16s4h2fC6HynE5FoBLF4DEMTQzhz5QyGTw0jVBPK7I+xFUhEfrVvn+w1CMhSkXNzMjUsD/i6Jdh/of9ZAASAWDyGSDSC/gv9mf8xtgKJyK+CQaCpyTofHXWvLBnydRAcGB94FgBNsXgMg+ODmf0htgKJyO/sa4qOjeXNvEFfB8Hulm4YASPhPiNgoKulK7M/ZF8eja1AIvKj2lprq6V4PHEPVQ/zdRDsO9GHYEnwWSA0AgaCJUH0nehL/48sLQH37lnn9jX1iIj8xN4avH3bWjrSw3wdBEM1IQyfGkbPsR50NXeh51hP5kkxd+9aW4ls3y43IiI/amyU8UFAAuCdO+6WJw2+zg4FJBCefvP05n5Z68RtRNrbs1Im8o/wTBj9F/oxMD6A7pZu9J3oyzwzmcgrlJLW4PCwnI+OynXRw2sn+z4IbsnEhDX4W1oK7Nzpbnkor2R1ig6RV7S2ysoxS0tyfRwfB0LefT97NzznA3srsK3N0992yHuyOkWHyCsCAaCjwzr3+HQJXrU3a24OePRIjpUCdu92tTiUf7I2RYfIa3btAoqK5Hh2FnjyxNXiPA+D4GbZW4FNTUBZmWtFofyUtSk6RF5jGLLVkim8heUoc4xBcDNiscT/VCbE0CZkZYoOkVfZxwHHx4GVFffK8hwMgpsxNmb9h1ZVAXV17paH8lJWpugQeVVdHVBZKcexmCQSehCzQzfD3hXKsUDagi1N0SHyulDI2mPw7l2gpcXd8qTAlmCmZmdlrVBABn7tW4gQEZElFJLEQUASCRcW3C1PCgyCmbIvlN3UBBSzMU1ElFJZmbXzPODJBBkGwUzZg6AHm/ZERJ7S1mYdh8Oy0paHMAhmYnoaePpUjg0j8RsOERGt1dgIlJTI8dOn1vxqj2AQzIS9FbhzJ1eIIaJ1hcNAby/Q1SU/PdgT6IxAILHXzGMVwQGtdGmduGWSfSIoEZFNOAx0dkoOXSwGDA0BZ87IutIeXkYzd0Iha3/Bhw9l5x2PNCK8UYp88PixLAgLyGLZ9fXuloeIPKu/3wqAgPyMROR+X6qpAcrL5TgWk+upRzAIpsveFdrcbKX9EhElGRiwAqApFgMG/bwsbFOTdfzggXvlSMIgmI54HLh/3zpnVygRPUd3t+TO2RmGjA/6lj0Iemj1GAbBdDx6ZH2tKy/n7vFE9Fx9fbLBuhkIDUPO+/y8LGxtrVUhT58CMzPulmcVg2A6Hj60jnfuZFcoET1XKCRJMD090vrr6fFxUowpEAB27LDOPdIaZHZoOiYnrWP7fyIR0TpCIeA0l4VN1NRk5VdMTAAHDrhbHrAluLGnT621QgMB7hhBRLRZO3ZYUyNmZ63FR1zEILgReyuwrs4zc1uIiPJOcXFiQ8IDXaK8om+EXaFERNnjsSxRBsHn0ToxCHKtUCKirbEHwakpmYLmIgbB55mZsaZGlJXJLvJERLR5ZWXWjvPxOPDkiavFYRB8HnsrkMukEREl2PQi4bW11vHUVE7Kli5Hg6BS6nWl1A2l1IhS6tvrPOekUmpIKXVNKfVTJ8u3BscDiYhSMhcJ/+53gYsX5WdnZ5qB0I9BUClVBOA7AN4AcBjAN5RSh5Oesw3A7wP4G1rrIwD+B6fKt8bysuwfaGJLkIjomS0tEm7PEJ2acnWjXSdbgl0ARrTWo1rrKICzAL6a9Jy/CeBHWuu7AKC1fgi3TE9bA7bV1bJzBBERAdjiIuGVldY1NRYD5uayXr50Ke1QBFZKfQ3A61rrb62efxNAt9b6bdtzfg+AAeAIgCoA/0Jr/Ucp/tZbAN4CgMbGxmNnz55N+ZqRSATBYHBT5S0fH0fFnTsAgMXGRszv2bOpv5NPtlJffsT6Sh/rKjP5UF/hsIwY2UOIUpJEn87ycFU3bqBkdUul+Y4OLNqzRjOUqr5ee+21S1rr4xv9rpPLpqVacDM5AhcDOAbgKwDKAfxCKfW+1vpmwi9p/S6AdwHg+PHj+uTJkylf8Pz581jvsQ198IG0AAHp6G5r29zfySNbqi8fYn2lj3WVmXyor+SNg81FwtNeI7WtDbh2TY5bWoBXXtl0WbZSX052h44BsFdNK4B7KZ7z51rrea31IwA/A9DpUPkS2dN2t21zpQhERF615UXC7eOCLm6y62RL8CKAfUqpdgDjAL4OGQO0+08A/pVSqhhACYBuAP+vg2UU0ai1pl0gIF9viIgowZYWCa+ulmXUlpeBxUVgYQGoqMhq+dLhWEtQa70M4G0APwHwEYA/1lpfU0qdUkqdWn3ORwD+HMBlAIMAvq+1vupUGZ+x73NVU8P1Qokop5Ln2yUnnBQkpRL3ZrVn4zvI0a2UtNbnAJxLuu+dpPN/DuCfO1muNexdoTU1rhWDiApf8tja0BCwZw9w8KAP9h/cts2aj+1ShiibOKlwPJCIHJJqvl08nuZ8u3xnX4qSQdBD2BIkIoekmm+ndZrz7fIdg6AHLS3JIC0AFBVx0WwiyqnubpleYKeUjA8WvGBQ/rEAMD8PrKw4XgQGwWTJSTEq1fRGIqLs6OuTWGAGQsOQXLy+PnfL5YhAIDEjNBJxvAiOJsbkBft/AluBRJRj5ny7/n7pAu3qAo4c8UFSjKmqSlqBgHSJOjwExSCYbGHBOnZhzgoR+U/yfLvz510rivOqqqwd5l0YF2R3aDJ7EDQ3fiQiotxwOTmGQTCZ2SwH2BIkIso1BkEP0ZrdoURETrJniC4sOJ4hyiBot7Rk7SFYUrI2b5mIiLIrOUPU3hBx4uUdfTWvYyuQiMh5ZWXWsTlP2yEMgnYcDyQicp49CC4tOfrSDIJ2zAwlInJeaal1zJagi9gSJCJyHrtDPcLcSBdgECQicgqDoEfYl3IvKXGvHEREfsIxQY+wB0FOjyAicgZbgh4RjVrHDIJERM5gYowHxOPWRHmlgGKuLU5E5IjiYuuaG4+v3WU4hxgETWwFEhG5x6XWIIOgiUkxRETusTc+lpcde1kGQROTYoiI3BOwhSNzaMqJl33eg0qp/6iUqnOqMK5iECQick9RkXXs4E4SG7UEfx3ANaXUrzlRGFcxCBIRucceBL3SEgQQB7ADwH9USp1RSm13oEzuYBAkInKPvTvUQy3BVwFcA6AAfB3SKvxVpVSb/ZbzUjrBXun2byRERJR7LrUEnzsZTms9qJR6GcD/DuDbAJoA/DD5aRv9nbxj7nJMRETO8GhLEFrrZQBnAHwMCXgqxY2IiGjzvNgSVEopAP8AwD8BUAYJeP8FwAe5LxoREfmG1taxg/MEN+rGfB/AcUjwmwbw97XW/zbnpSIiIn8ZH7eOb9wA9u935GU3CoKfWf15DsDf1Vrfz3F5iIiIHLPRmOAMgP9Za/3XGQCJiChnXEpI3CgIvqC1/kMnCkJERD4WClnHBw869rLPDYJa6zGnCkJERD5mzwh1cK42F9AmIiL32ecGBpwLTQyCqdhTdYmIKPdcWrWLQdDk0grmRESExO5QtgRdYN9I176YNhER5R5bgi6z7xwRjbpXDiIiP2JijMvsQZAtQSIiZzExxmUMgkRE7mFL0GUMgkRE7rFfd4ud252PQdDEIEhE5J7FReu4rMyxl2UQNAUCVhNca0e38iDyuvBMGL3netH1vS70nutFeCbsdpGokMRi1phgUZGjLcHC2hF+qwzD+o+IxRz9jyDyqvBMGJ3vdCISjSAWj2FoYghnrpzB8KlhhGpCG/8Boo0sLVnHDrYCAbYEE3GuINEa/Rf6nwVAAIjFY4hEI+i/0O9yyahguNQVCjAIJrKPC9q/mRD52MD4wLMAaIrFYxgcH3SpRFRw7EGwtNTRl2YQtKuosI4XFtwrB5GHdLd0wwgYCfcZAQNdLV0ulYgKDluCHsEgSLRG34k+BEuCzwKhETAQLAmi70SfyyWjgsExQY+orLSO5+fdKweRh4RqQhg+NYyeYz3oau5Cz7EeJsVQdrnYEmT6ox1bgkQphWpCOP3mabeLQYWKY4IeYQ+CbAkSETmDY4IeUVpqTZhfXt7SbhLhMNDbC3R1yc8w5xYTEa21smL1vCmV2BhxALtDk1VWArOzcrywkDh3ME3hMNDZCUQiMt1waAg4cwYYHgZCHEYhIrJEItZxRYWjO0gAbAmulYVxwf5+KwAC8jMSkfuJiMhmbs46rqpy/OUZBJNlYVxwYGDtgjOxGDDIucVERIkYBD3GPk1iky3B7u7ExWcAOe/i3GJyGMemyfNcDoIcE0xmD4Lm2GCG+vpkDNDsEjUMIBiU+4mcwrFpygv2IFhd7fjLsyWYbNs263h2NnG34zSFQnKh6emRb+A9PbzwkPM4Nk2et7ycmBlqb4Q4hC3BZIYh44ILCxIA5+aAmpqM/0woBJzm3GJyEcemyfPsmaGVlY5nhgJsCaZmbw0+eeJWKYi2hGPT5HkujwcCDIKpMQhSAejrk7FoMxBybJo8h0HQo+zdnzMz7pXDAfbswXCY2YOFhGPT5Hn2RoYLSTEAxwRTswdBMznGhb7qXEvOHpyclHNeKAsHx6bJs+JxYHraOt++3ZViOHplV0q9rpS6oZQaUUp9+znP+4xSakUp9TUny/eMYVhZSlpveqqE1yVnD2rN7EEicsiTJ1b2fWWl4wtnmxwLgkqpIgDfAfAGgMMAvqGUOrzO8/4ZgJ84VbaUfDAuyOxBInLN1JR1XFvrWjGcbAl2ARjRWo9qraMAzgL4aorn9QL4EwAPHSzbWvYgWKDjgsweJCLXeCQIKq21My8kXZuva62/tXr+TQDdWuu3bc9pAfDvAHwZwL8G8P9prX+Y4m+9BeAtAGhsbDx29uzZlK8ZiUQQDAY3Vd7imRnUXLsGAFgpK8OTV17Z1N/xslgMuHZNeiS0BlpbI7h/P4gjR9YGR1prK+8vv2FdZabg60tr1F68CLW8DACYfvllxMvLN/3nUtXXa6+9dklrfXyj33UyMUaluC85Av8egN/RWq8olerpq7+k9bsA3gWA48eP65MnT6Z83vnz57HeYxuKx4HiYqvPursb2MJ/klcdPChjgIODwNtvn8ev//pJJsWkaUvvL59hXWWm4OtrdtaaKF9aCvzKr2zpz22lvpwMgmMA7JfXVgD3kp5zHMDZ1QBYD+BNpdSy1vpPHSmhXSAA1NVJyiQgP9vaHC9GrtmzB8+fZ1YoETnAI12hgLNjghcB7FNKtSulSgB8HcCP7U/QWrdrrXdrrXcD+CGA33YlAJp27LCOzWBIRERb46Eg6FhLUGu9rJR6G5L1WQTgB1rra0qpU6uPv+NUWdJWX28dT07KwNlzummJiCgNjx9bx34JggCgtT4H4FzSfSmDn9b67zhRpueqrpb+6qUlySKZmUnMGiUioszMzwOLi3JcXLypDQqyqfCWQcm2hgbrmF2iRERbMzFhHdfVud67xiC4EY4L+g53YyfKIXsQbGpyrxyruHboRuzjgtPTsglkMautUHE3dqIcikYT1wttbHSvLKvYEtxIaam1unk8njigSwWHu7ET5dCDB5JgCEhCTGmpu+UBg2B67OOCDx64Vw7KOa6nSpRD9q5QD7QCAQbB9Nj/s+7ft77JUMHheqpEObKykphX4YHxQIBBMD21tdY2H9Eo8OiRu+WhnOFu7EQ58uiRBEJAPlQeWRuVQTAdSgHNzdb5+Lh7ZaGc4m7sRDnisaxQE9Mc09XcDIyOyvHERMHuNk/cjZ0o67ROzKfwyHggwJZg+rZvByoq5DgWAx66u90hEVHemJ6WlbcAyQjdvt3d8tgwCGaipcU6vpe8AQYREaVkX3Giqcn1VWLsGAQzYR8XnJiwBnmJiCi15eXERoPHBtgZBDNRXQ1UVcnxykriQC8REa11/74EQkAyQj3UFQowCGbO3hpklygR0fPdvWsde3BjcgbBTNnHBR8+lHmDRES01vy8tYGuUkBrq7vlSYFBMFOVldaegvF44rccIiKy2BNiGhs9sVZoMgbBzdi92zq+fZvLqBERJdM6MQh6LCHGxCC4Gc3N1rpaT59yziARUbLJSWsH+dLSxL1ZPYRBcDOKioBdu6zzW7fcKwsRkRfZh4paWz27wpY3S5UP7EFwclIGgImISFqA9mXSPNoVCjAIbl5FReIisLdvu1YUIiJPuXVLEgcB2YXHnF/tQQyCW2FPkLl715oQSkTkV8vLwJ071vmePe6VJQ0MgltRXy9TJgD5j+cWS0Tkd3fvyiYDgFwfPbRjRCoMgluhFNDebp0zQYaI/CweBz791Drfs8dTi2WnwiC4Va2tki0KAHNzkiRDRORH9+4lTovwcEKMiUFwqwwj8T/65k33ykJE5CZ7K7C93bPTIuy8X8J8sHev9Z89NcXWIBH5z8OHwOysHBcVJSYOehiDYDaUlyeujs7WIBH5jb0VuGuXtaqWxzEIZgtbg0TkV0+eAI8eybFSQEeHq8XJBINgtrA1SER+Zb/eNTfL9TBPMAhmE1uDROQ3jx8nLpG2d697ZdkEBsFsYmuQiPzmo4+s49ZWoLravbJsAoNgtrE1SER+cf8+MD0tx4EAcPCgu+XZBAbBbGNrkIj8IB5PbAW2t+fVWKCJQTAXkluD9++7Wx4iomy7e9faQs4wgH373C3PJjEI5kJ5eeJE0WvXgJUV14pDRJRVy8uJvVx79+bNvMBkDIKrwmGgtxfo6pKf4fAW/+D+/UBJiRw/fQqMjGy5jEREnjA6CiwtyXFZWeJGAnmm2O0CeEE4DHR2ApGI7AAyNAScOQMMD29h/VfDAA4dkj8CSBBsa8vLPnMiomeWlhJXhzl40NpEIA+xJQigv98KgID8jETk/i0JhYCaGjmOx6VblIgon924YW0gXlUl0yLyGIMggIEBKwCaYjFgcHCLf1gp4OhR6/z+fWtpISKifDM9nbhr/KFDnt8vcCMMggC6u9eO6RqGjA9uWW1t4jelq1cBrbPwh4mIHBSPW8M7gOwY7/Fd49PBIAigrw8IBq1AaBhy3teXpRc4dChx493bt7P0h4mIHPLpp3L9AuR69sIL7pYnSxgEIUN3w8NAT4+0/np6tpgUk6ysTLJFTTduANFolv44EVGOzc8nTok4eLBgkvyYHboqFAJOn87hC3R0WJNLYzHpFn3llRy+YHrCM2H0X+jHwPgAeht7EZ4JI1STrehPRAXhyhXpDgUk2S+Pp0QkY0vQKYEAcOSIdT4+DkxMuFceSADsfKcT3730XVy8dxGT85PofKcT4ZmtTpIkooIxNmatgayUzCfL82QYOwZBJzU2JibJXL7sardo/4V+RKIRxOKSGquhEYlG0H9hq3NDiKggRKOJU7va261pXwWCQdBpR4/KGCEgk06vXHGtKAPjA88CoCkWj2FwfKtzQ4ioIFy/bn1RLy8HDhxwtzw5wCDoNMOQ7gTTvXtyc0F3SzeMQOLcECNgoKslG3NDiCivPXqUuH7kCy8AxYWXRsIg6IYdOxK3W7p82VqHz0F9J/oQLAk+C4QKCsGSIPpOZGtuCBHlpWgU+PBD63znzoKYE5gKg6BbjhyxUoxjMQmEDgvVhDB8ahg9x3rQ1dyFhsoGDJ8aZnYokd8NDQGLi3JcUpK48lWBYRB0S3Ex8NJL1vnEhGRhOSxUE8LpN09j4O8OIFQdYgCkLQvPhNF7rhdd3+tC77leZhvnm9u3gQcPrPOXX7byGApQ4XXw5pP6etl30FxB5upVua+A33BU2MxpN2bW8dDEEM5cOYP/cPw/uF00SsfsbGI2aEeHDN8UMLYE3XboEFBRIcexGHDpkjUplSjPJE+7icVjiEQjmJh3d04spWFlBfjlL63rT3W1XJ8KHIOg24qLpbvBnHw6NQV8/LG7ZSLapPWm3cxH510qEaXt+vXEtUGPHZNFPgpc4f8L80FtbeL8m08/lW2XiPLMetNuKksqXSoRpWViInFh/6NHZRcBH2AQ9Iq9exNTkIeGZGdfojySPO3GCBgIlgTRVNnkcsloXYuLcr0x7dyZOIWrwDEIeoVS0i1qjg8uLwMffCD99ER5InnaTc+xHgyfGoZRZGz8y+S8eFzyEMxdxcvLExfz8AFmh3qJYQCf+Qzw3nvy5pybkz2dPLDbBFG6zGk3dp/iU5dKQ8915YrkIQDyRfyVV9buMF7g2BL0mupq4MUXrfPxcW7CS0TZNzoq27uZDh2S/ASfYRD0olAosU/+2jVgetq98hBRYZmclGxQUygE7NnjXnlcxCDoVS+8YG1ZEo/L+KC5jBER0WZFIjIOqLWcb9+e2PvkMwyCXhUIAMePW/3zi4vAwIAkzBARbUYsBgwOWokwZWWSh+CD+YDr8e+/PB9UVMiEVXMi/eystAi5ogwRZUpraQHOry5cUFQEdHUBpaXulstljgZBpdTrSqkbSqkRpdS3Uzz+t5RSl1dvP1dK+StXN5WGhsSuislJV3acIKI8d+2aXD9ML79ccLvEb4ZjQVApVQTgOwDeAHAYwDeUUoeTnnYLwJe01i8C+F0A7zpVPk9rawP277fOw2Hgxg33ykNE+eX2beDWLet8/36ZFE+OtgS7AIxorUe11lEAZwF81f4ErfXPtdZmGuT7AFodLJ+3HTggGVymmzcT05uJiFIZH5f5gKadOxO/VPuc0maGUK5fSKmvAXhda/2t1fNvAujWWr+9zvP/VwAHzecnPfYWgLcAoLGx8djZs2dTvmYkEkGwkNa/i8dR/fHHMJ48kXOlMHvwIGLbt2flzxdSfcVWYpiYn8B8dB6VJZVoqmzK+qolhVRfuca6yky26suYmkL1jRvPMkGXg0HMHDki44EFJFV9vfbaa5e01sc3+l0nV4xRKe5LGYGVUq8B+C0An0/1uNb6Xax2lR4/flyfPHky5QueP38e6z2Wt774ReDCBUmSAeTNnKW+/UKpr+Q97cz1K4dPDWd10+BCqS8nsK4yk5X6evwYeP99a1f4qirg1Vdlp/gCs5X6crI7dAyA/QrUCuBe8pOUUi8C+D6Ar2qtHztUtvxRXAx0d8saf4CsLTowYGV80bp72vVf6He5ZEQOefJEpkKYmeQVFcBnP1uQAXCrnAyCFwHsU0q1K6VKAHwdwI/tT1BKtQH4EYBvaq1vOli2/FJWJoHQnEO4tAT8/OcMhKvW29NucHzQpRIROWhuTlqA5pzisjLgc5+Tn7SGY0FQa70M4G0APwHwEYA/1lpfU0qdUkqdWn3aPwJQB+D3lVJDSqkPnCpf3qmqkkmuZt/+4iID4ar19rTraulyqUREDpmfB37xC2syfEmJtADN3WloDUfnCWqtz2mt92ut92it/+nqfe9ord9ZPf6W1nq71vql1duGg5q+Vlcnk10ZCBOst6dd34k+l0tGlEOLi9ICXFqS8+JiCYBVVe6Wy+O4Yky+q69nIEyy3p522UyKyTfhmTB6z/Wi63td6D3Xi/BM2O0iUTYtLkoLcGFBzgMBuS5wMvyGuJ9gITAD4eCgJMqYgfDVV4HKSrdL54pUe9r5VXK27NDEEM5cOeP7LwYFY2EhMQAqJUMldXXulitPsCVYKNgipHUwW7aAzc3JlCl7ADx2DNixw91y5REGwULCQEgpMFu2QD15IgHQ3GItEJAWIJdDywiDYKFJFQgvXABmZtwtF7mG2bIF6NGjxCxQMwmmsdHdcuUhBsFCVF8v8wjNQGjOI7SvIE++wWzZAvPgQeLeoiUlMv7PMcBNYRAsVHV1iRPql5flg8NFtxEOA7290mDu7ZXzQsZs2QIyPg5cvGitBFNWJgGQWaCbxuzQQlZXB5w4IcHv6VNZRHd4WI4PHHC7dK4Ih4HOTiASkZ6koSHgzBmpllABxwRmyxaA27cTd4OorORE+CxgS7DQVVUBn/88UF1t3Xfzplz9fbhDfX+/FQAB+RmJyP1EnqS1bIhrD4DmYtgMgFvGIOgHZWXSImxosO4Lh2VeoTmu4BMDA1YANMViUhVEnhOLyZt2dNS6b/t2+TxzLdCsYBD0i+JiGQSz9/lNTiamWPuAfZjUZBhSNUSeEokA772XmNC2c6d0gSa/iWnTGAT9JBAAXnopcTxwdlY+aOZGvQWurw8IBq1riGHIeR8TJclDjCdPgL/6q8Q5vvv3y0T4YqZyZBODoB/t3y/ZIWp1n+PVuYSlExPulssBoZAkwfT0SOuvp6fwk2Ioz4yOovqjj6x++6IiCX4HDlifWcoafqXwq7Y22Zj30iX5sMXjCI6OAh9+CLz4ojXHsACFQsBpJkqS18Tjkvxy964kwwAy7seFsHOKLUE/a2gAvvCFxMzRsbG13TBElFtLS7ICjH0e7/btwBe/yACYYwyCfldZKVMo7P2Bs7PAz34G+KB7lMh1Dx8CP/0pMDX17K6lhgaZAlFa6mLB/IHdoSRdny+9hMiNG5I8E4/L1ImLF4G9e4GDBzkWQZRt8Tjw0UeJ0x8A4PBhRKqq5LNIOcdapmeWGhulVWifgDsyIt005m7VRLR18/My7GAPgKWlMv1hzx73yuVDDIKUqKZGxiHs+5E9fgycP8/uUaJsCIel+9O+s0tjI3DyZOKCFuQIBkFay5w9fvCgdV80Kt2jH364dskVItpYLAb88peyZOHKitwXCABHj8rnraTE1eL5FccEKTWlgH37gG3b5ENrriozNiZ7mXV2cvdqonRNT0sANHeAB2SVhmPHErOzyXFsCdLzNTRIN01rq3Xf4qKsZzg87Lu1R4kysrIiyS8XLiQGwLY2GXZgAHQdW4K0McMAXn5Z1i28fNlKkrl712oV1te7W0Yir5mclM+LPfgZhnxedu50r1yUgEGQ0tfUBNTWygf7/n25b2FBskfb24FDhwp6pRmitESjsvXR2Fji/XV18mWyvNydclFKDIKUmZIS4Phx2eH6yhUrSebWLeDBAxnkb2x0t4xEbhkbkwAYjVr3GQZw+LAsSMH5tp7DIEib09Ii32wvX5bgB0ircHBQEmaOHpXVaIj8YH5evhTatz0CgOZm+Sxw5RfPYhCkzTMX9w2H5duv2Sp8+FDmFXZ0SIYpt36hQhWPy4T3Gzfk2FReDrzwAntF8gCvTrR1oZB82D/+GLhzR+6Lx2W1mbEx6QpqaXG3jETZNjEhmZ+RiHWfUjI+fuAAv/zlCf4vUXaUlMgWTLt2SbfQ9LTcv7go86Pu3JFuIaaEU76bngauX09Y8BqAvLc7O2VuLeUNBkHKrpoa4MQJSZy5ft2aTvH4sexMsXu3bOrL1TEo38zPS2/HvXuJ9xcXy3u6o4OJL3mIQZCyTymZXN/YCNy8KZmjWsvt1i0ZQ2xvl4WCDcPt0hI9XzQKfPIJcPt24rifUvxSVwAYBCl3DAM4ckS6SK9etTLnlpflonLrlgTC9nYGQ/KelRV5j46MrF0vd+dOmRfLDOi8xyBIuRcMyhYxExPSnTQ3J/cvL0tW3eioFQyZTEBuW1mR3oqREeDp08THamsl0Wv7dnfKRlnHKw45p6lJukjv35fgZ2bVxWISHEdHZRPf3bu58gw5LxaTLs9bt9bun1lZKcGvqcmVolHuMAiSs5SSCcQ7d0ryzM2bknAAyNjL9evAp59KMGxrY8uQcm9pSb6A3b69dkH4khKZ7tDWxp3eCxSvMOQOM3mmpUW6nj75xFpoeGlJJt/fvCkXn927E3e7J8qGhQX5wnX3bmLCCyALQezZI+PZ7JUoaAyC5C6lJNC1tkowvHnT2rswFpOL1OiodKN2dMhSbURbMTcn433j45KxbBcMSi9ESwtbfj7BIEjeEAjIt+5QSL6Zj45a3aRaS1LNxIRMSO7o4EWKMhOPy1j0nTsyZzXZtm0S/JqaONfPZxgEyVsCAen+3LVLplSMjiYuSjw7KzvdX79uPa+szKXCkudFIhL4xsYSd3Yw1ddL8GtocL5s5AkMguRNSsluFDt2SPfVrVtyIVtZkcejUek6HRmR55iT8x1sHcZWYug914uB8QF0t3Sj70QfQjUhx16f1rFRq08pafHt3cslzohBkPJAVZWsS3rokFzYbt+25m/F41ZXqWFI5mlrq8znyqHwTBjXJq/hu598F7F4DEMTQzhz5QyGTw0zELrFbPWFw2sntwOys4PZ5c7eA1rFIEj5wzDk2/uePfJN/9atxEWMYzG5CN65I9mkra1yy8GqHv0X+tGu2xGLy8U2Fo8hEo2g/0I/Tr95OuuvR+uYn5f3wv37wJMnax83W327dknXJ8f7KAmDIOUfc65hc7N8+x8bk0w/c4oFIMc3b8pt+3YJhs3NWVvjcWB8ALuDuxPui8VjGBwfzMrfp+eYn5dFrO/dkzHiVCoqJOuYrT7aAIMg5bdgEDh4UCY0T09LV9j9+4ndYdPTcrt6VQLijh0yfriFbZ26W7qhZhJbFUbAQFdL16b/Jj1HJCJB7/799QNfICD/r2z1UQYYBKkwKCXjgLW1sqP3xIS0EB8+tOaCaS3dp1NTskxbeblcNBsbZf5hBpOi+0704U/+7E9gBAzE4jEYAQPBkiD6TvTl6B/oM/G4dG9OTkrgM9ebTRYISGZnc7P8P3IhdsoQgyAVnkDA6i6NRqWr9N69tZugPn0qSTa3b0sArK+3guIGXWihmhCONBxBT3UPBscH0dXSxezQrdBaAt3kJPDokWR1mpnAyQIBac2bgY9L69EW8N1Dha2kRHanaG+XgPjwIfDggVxs7V2mKyty/4MHcl5ZabUsa2ul2zWJUWQwCWYrFhYk4JmBL9U8PlNRkRX4duxg4KOs4TuJ/KOkxMoYjcdlnNAMfOaOFqb5ebmFw3JuGIlBkfPLMrOyIi3xmRm5PX6cmMiUSkWFtM7N+aJcw5NygEGQ/CkQkHHAujrZImd+XoLhw4dygU5eUDkWS2wpBgKoGR2V8aiaGpnLGAzyQg1IwDOD3cwM8OQJ6gYHNw56JSUS9Boa5CcXTScHMAgSAdL92dEhNzMpY2pKWotTU2u76uJxFM/NybJudhUVEhDtt0INjsvLEtjm5+Xn7KwEvUhk7cLUyeeA1EldnRX4qqqY0UmOYxAkShYIWN2epkjEyiydmrIW9062sCA3s8UIyIXdDI5lZXIrL7eOy8q8O8a1tGQFueSfyRvPbiQYlFbztm3WjYugk8s8+skj8phgUG5tbXIejWIuFgP27ZOsxrk5CQypWjxaW2OM6ykuXhsci4vlVlS08U/zdey3eDz1fcvL0rKNRqWb1zxOvsVia7uF06GUFfBWg95URQXw2muZ/y2iHGMQJNqMkhJEa2tlor4pHpcWoxkUNwqOdsvL8rvJCTpeFghI0K6stG5m4Evq/tWF2B1MBYFBkChbAgFZhSZ5JZqVFWkFRiKyYbB5e/rUOt5Mi8sJxcUS3CoqrEBnHpeVcQyP8h6DIFGuFRWlDo520WhiYFxaktbhykp6P5VKfQsE1t5XXCyZmIYBlJbKz5KStTfDKMyEHiIbBkEiLzADzxbWMyWizDE1i4iIfItBkIiIfItBkIiIfItBkIiIfItBkIiIfItBkIiIfItBkIiIfMvRIKiUel0pdUMpNaKU+naKx5VS6l+uPn5ZKfWKk+UjIiJ/cSwIKqWKAHwHwBsADgP4hlLqcNLT3gCwb/X2FoA/cKp8RETkP062BLsAjGitR7XWUQBnAXw16TlfBfBHWrwPYJtSaqeDZSQiIh9xctm0FgBh2/kYgO40ntMC4L79SUqptyAtRTQ2NuL8+fMpXzASiaz7GK3F+soM6yt9rKvMsL4ys5X6cjIIplpuPnl/mXSeA631uwDeBYDjx4/rkydPpnzB8+fPY73HaC3WV2ZYX+ljXWWG9ZWZrdSXk92hYwBCtvNWAPc28RwiIqKscDIIXgSwTynVrpQqAfB1AD9Oes6PAfzGapboZwHMaK3vJ/8hIiKibHCsO1RrvayUehvATwAUAfiB1vqaUurU6uPvADgH4E0AIwAWAPymU+UjIiL/UVqvGXLLK0qpSQB31nm4HsAjB4uT71hfmWF9pY91lRnWV2ZS1dcurXXDRr+Y90HweZRSH2itj7tdjnzB+soM6yt9rKvMsL4ys5X64rJpRETkWwyCRETkW4UeBN91uwB5hvWVGdZX+lhXmWF9ZWbT9VXQY4JERETPU+gtQSIionUxCBIRkW/lfRDkHoWZSaO+/tZqPV1WSv1cKdXpRjm9YqP6sj3vM0qpFaXU15wsn9ekU19KqZNKqSGl1DWl1E+dLqOXpPF5rFFK/Wel1PBqffl2ARGl1A+UUg+VUlfXeXxz13qtdd7eICvPfAqgA0AJgGEAh5Oe8yaAP4Mszv1ZAANul9vj9fUqgO2rx2+wvp5fX7bn/XfIikdfc7vcXq4vANsAXAfQtnq+w+1ye7y+/jcA/2z1uAHAFIASt8vuUn19EcArAK6u8/imrvX53hLkHoWZ2bC+tNY/11pPr56+D1nE3K/SeX8BQC+APwHw0MnCeVA69fU3AfxIa30XALTWfq6zdOpLA6hSSikAQUgQXHa2mN6gtf4Z5N+/nk1d6/M9CK63/2Cmz/GLTOvityDfrPxqw/pSSrUA+FUA7zhYLq9K5/21H8B2pdR5pdQlpdRvOFY670mnvv4VgEOQ3XSuAPj7Wuu4M8XLO5u61ju5n2AuZG2PQp9Iuy6UUq9BguDnc1oib0unvn4PwO9orVfky7qvpVNfxQCOAfgKgHIAv1BKva+1vpnrwnlQOvX11wAMAfgygD0A/qtS6j2t9WyOy5aPNnWtz/cgyD0KM5NWXSilXgTwfQBvaK0fO1Q2L0qnvo4DOLsaAOsBvKmUWtZa/6kjJfSWdD+Pj7TW8wDmlVI/A9AJwI9BMJ36+k0A/7eWQa8RpdQtAAcBDDpTxLyyqWt9vneHco/CzGxYX0qpNgA/AvBNn347t9uwvrTW7Vrr3Vrr3QB+COC3fRoAgfQ+j/8JwBeUUsVKqQoA3QA+cricXpFOfd2FtJqhlGoEcADAqKOlzB+butbndUtQc4/CjKRZX/8IQB2A319t3Sxrn65mn2Z90ap06ktr/ZFS6s8BXAYQB/B9rXXKlPdCl+b763cB/KFS6gqku+93tNa+3GJJKfXvAZwEUK+UGgPwjwEYwNau9Vw2jYiIfCvfu0OJiIg2jUGQiIh8i0GQiIh8i0GQiIh8i0GQiIh8i0GQiIh8i0GQKA8ppf6xUkorpeJKqS/b7q9XSk2sPvaeUoqfcaLn4DxBojyklCoGcAGyE0EYwAta6xml1A8B/DqAOQCdWutbLhaTyPMYBInylFJqH2Rx5QoAZyArj/zR6sO/pbX+gUtFI8obDIJEeWx1ia0/WD1dBFAG4E+11r/qXqmI8geDIFGeU0r9GYDXV08nIbuT+3J9SaJMcdCcKI8ppcqQuH1MNYBGl4pDlHcYBIny2/8F4AiAZQDjAEoB/NvVrXmIaAMMgkR5Sin1GoD/ZfX0nwL4HwGsAHgJwD9xp1RE+YVjgkR5SClVA+AKpCt0EMCJ1f3p/k8A/xCyV98XtdYXXCwmkecxCBLlIaXUvwHwtyGbh76stb65er8B4BcAjkF2IO/UWkdcKyiRxzEIEhGRb3FMkIiIfItBkIiIfItBkIiIfItBkIiIfItBkIiIfItBkIiIfItBkIiIfItBkIiIfOv/BzSpMevx+BTrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A=['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf']\n",
    "plt.figure(figsize=(7,7))\n",
    "x0, y0 = 0.5, 0.5 \n",
    "\n",
    "r= value(instance.R)\n",
    "theta = np.linspace(0,2*np.pi,100)\n",
    "X= x0+r*np.cos(theta)\n",
    "Y= y0+r*np.sin(theta)\n",
    "\n",
    "plt.plot(X,Y,c='r', lw=3, alpha = 0.3)\n",
    "\n",
    "for i in instance.i:\n",
    "    x0 = df.loc[i-1,'x']\n",
    "    y0 = df.loc[i-1,'y']\n",
    "    if df.loc[i-1,'c'] == 0 :\n",
    "        KOLOR= 'b'\n",
    "    else:\n",
    "        KOLOR= 'g'\n",
    "    plt.scatter(x0,y0, s=30 , c=KOLOR )\n",
    "\n",
    "plt.xlabel('  X  ' , fontweight = 'bold', fontsize= 14)\n",
    "plt.ylabel('  Y  ', fontweight = 'bold', fontsize= 14)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# ==========================================================\n",
      "# = Solver Results                                         =\n",
      "# ==========================================================\n",
      "# ----------------------------------------------------------\n",
      "#   Problem Information\n",
      "# ----------------------------------------------------------\n",
      "Problem: \n",
      "- Name: unknown\n",
      "  Lower bound: 0.464918075924774\n",
      "  Upper bound: 0.464918075924774\n",
      "  Number of objectives: 1\n",
      "  Number of constraints: 53\n",
      "  Number of variables: 28\n",
      "  Number of nonzeros: 128\n",
      "  Sense: maximize\n",
      "# ----------------------------------------------------------\n",
      "#   Solver Information\n",
      "# ----------------------------------------------------------\n",
      "Solver: \n",
      "- Status: ok\n",
      "  Termination condition: optimal\n",
      "  Statistics: \n",
      "    Branch and bound: \n",
      "      Number of bounded subproblems: 19\n",
      "      Number of created subproblems: 19\n",
      "  Error rc: 0\n",
      "  Time: 0.022581815719604492\n",
      "# ----------------------------------------------------------\n",
      "#   Solution Information\n",
      "# ----------------------------------------------------------\n",
      "Solution: \n",
      "- number of solutions: 0\n",
      "  number of solutions displayed: 0\n"
     ]
    }
   ],
   "source": [
    "results.write()"
   ]
  }
 ],
 "metadata": {
  "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.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}