CG-PMP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A column generation algorithm for the $p$-median problem (PMP) [1]\n",
    "\n",
    "Given a graph $G = (V,E)$...\n",
    "\n",
    "$$\n",
    "  \\begin{align}\n",
    "    \\min & \\sum\\limits_{i \\in S}\\sum\\limits_{j \\in V} d_{ij} x_{ij} \\\\\n",
    "    \\text{subject to} & \\nonumber \\\\\n",
    "    & \\sum\\limits_{j \\in V} x_{ij} = 1, & \\forall i \\in V \\\\\n",
    "    & \\sum\\limits_{j \\in V} x_{jj} = p, \\\\\n",
    "    & x_{ij} \\leq x_{jj}, & \\forall i, j \\in V \\\\\n",
    "    & x_{ij} \\in \\{0, 1\\}, & \\forall i, j \\in V.\n",
    "  \\end{align}\n",
    "$$\n",
    "\n",
    "..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "import random\n",
    "import sys\n",
    "\n",
    "import numpy as np\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Read data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0  7 10 17 13 29 17 25 20 17 19 23 33 31 24]\n",
      " [ 7  0 15 16 12 28 22 24 15 11 12 16 31 25 17]\n",
      " [10 15  0  9  7 21 15 17 12 16 19 23 25 23 24]\n",
      " [17 16  9  0  8 14  8 10 11 15 20 24 18 16 24]\n",
      " [13 12  7  8  0 18 14 14  9 11 14 18 22 20 19]\n",
      " [29 28 21 14 18  0 16  6  7 22 17 13 10  4 12]\n",
      " [17 22 15  8 14 16  0 12 17 21 26 27 20 18 26]\n",
      " [25 24 17 10 14  6 12  0 11 21 21 17 10  8 16]\n",
      " [20 15 12 11  9  7 17 11  0  6 11  8 15  5  7]\n",
      " [17 11 16 15 11 22 21 21  6  0  7 11 26 20 12]\n",
      " [19 12 19 20 14 17 26 21 11  7  0  6 21 15  7]\n",
      " [23 16 23 24 18 13 27 17  8 11  6  0 17 11  3]\n",
      " [33 31 25 18 22 10 20 10 15 26 21 17  0 12 11]\n",
      " [31 25 23 16 20  4 18  8  5 20 15 11 12  0 10]\n",
      " [24 17 24 24 19 12 26 16  7 12  7  3 11 10  0]]\n"
     ]
    }
   ],
   "source": [
    "# Distance matrix\n",
    "distMtx = np.array([row[1:] for row in list(csv.reader(open(\"./data/scotland_yard_adj_mtx.csv\", \"r\"), delimiter=\",\"))[1:]]).astype(\"int\")\n",
    "\n",
    "print(distMtx)\n",
    "\n",
    "EPS = 0.00001\n",
    "\n",
    "# Number of vertices\n",
    "n = len(distMtx)\n",
    "\n",
    "# number of medians\n",
    "p = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAggAAAGKCAYAAABpbLktAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAQvUlEQVR4nO3dP04jaRrA4demAydgiZYQHm0HqDOLiACpD9ASiW/ABTgF8im4ABdYOfENRiJAG1jORgSrbbeQBsk4cWLXBj30wL54BryG8p/nyVyu7n5VcuOfXVUflaIoigAAeKRa9gAAwPIRCABAIhAAgEQgAACJQAAAEoEAACQCAQBIBAIAkHyY9w9Op9P49u1bbG9vR6VSWeRMAMAbKYoiRqNR/PLLL1Gtzv6eYO5A+PbtW3z69GnePw4AlOjf//53/OMf/5j5/NyBsL29/fMf2NnZmfevAQDe0f39fXz69Onn+/gscwfCw2mFnZ0dgQAAK+bvLg9wkSIAkAgEACARCABAIhAAgEQgAACJQAAAEoEAACQCAQBI5l4oaR1MpkVc3dzF7Wgce9u1OD7Yja2q3yuxCI4twGrb2EDo9gbR7vRjMBz/3Nao1+K81YyTw0aJk60+xxZg9W3kKYZubxBnl9dP3sAiIr4Px3F2eR3d3qCkyVafYwuwHjYuECbTItqdfhTPPPewrd3px2T63B78FccWYH1sXCBc3dylT7ePFRExGI7j6ubu/YZaE44twPrYuEC4Hc1+A5tnP/7k2AKsj40LhL3t2kL340+OLcD62LhAOD7YjUa9FrNuuKvEjyvujw9233OsteDYAqyPjQuErWolzlvNiIj0Rvbw+LzVdM/+HBxbgPWxcYEQEXFy2IiL06PYrz/9qnu/XouL0yP36v8fHFuA9VApimKue87u7++jXq/HcDiMnZ2dRc/1Lqz293YcW4Dl9NL3741dSTHix1fiXz5/LHuMteTYAqy2jTzFAAD8NYEAACQCAQBIBAIAkAgEACARCABAIhAAgEQgAACJQAAAEoEAACQCAQBINvp3MQD8P/xSMt7CsryuBALAHLq9QbQ7/RgMxz+3Neq1OG81/Vpz5rZMryunGABeqdsbxNnl9ZMf4hER34fjOLu8jm5vUNJkrLJle10JBIBXmEyLaHf6UTzz3MO2dqcfk+lze8DzlvF1JRAAXuHq5i59wnusiIjBcBxXN3fvNxQrbxlfVwIB4BVuR7N/iM+zH0Qs5+tKIAC8wt52baH7QcRyvq4EAsArHB/sRqNei1k3nVXix1Xnxwe77zkWK24ZX1cCAeAVtqqVOG81IyLSD/OHx+etpvUQeJVlfF0JBIBXOjlsxMXpUezXn37du1+vxcXpkXUQmMuyva4qRVHMdc/E/f191Ov1GA6HsbOzs+i5AJbesqx4x3p569fVS9+/raQIMKetaiW+fP5Y9hismWV5XTnFAAAkAgEASAQCAJAIBAAgEQgAQCIQAIBEIAAAiUAAABKBAAAkAgEASAQCAJAIBAAgEQgAQCIQAIBEIAAAiUAAABKBAAAkAgEASAQCAJAIBAAgEQgAQCIQAIBEIAAAiUAAABKBAAAkAgEASAQCAJAIBAAgEQgAQCIQAIBEIAAAiUAAABKBAAAkAgEASAQCAJAIBAAg+VD2ACyPybSIq5u7uB2NY2+7FscHu7FVrZQ9FgAlEAhERES3N4h2px+D4fjntka9FuetZpwcNkqcDIAyOMVAdHuDOLu8fhIHERHfh+M4u7yObm9Q0mQAlEUgbLjJtIh2px/FM889bGt3+jGZPrcHAOtKIGy4q5u79M3BY0VEDIbjuLq5e7+hACidQNhwt6PZcTDPfgCsB4Gw4fa2awvdD4D1IBA23PHBbjTqtZh1M2MlftzNcHyw+55jAVAygbDhtqqVOG81IyJSJDw8Pm81rYcAsGEEAnFy2IiL06PYrz89jbBfr8XF6ZF1EAA2kIWSiIgfkfC1uW8lRQAiQiDwyFa1El8+fyx7DACWgFMMAEAiEACARCAAAIlAAAASgQAAJAIBAEgEAgCQCAQAIBEIAEAiEACARCAAAIlAAAASgQAAJAIBAEgEAgCQCAQAIBEIAEAiEACARCAAAIlAAAASgQAAJAIBAEgEAgCQCAQAIBEIAEAiEACARCAAAIlAAAASgQAAJAIBAEgEAgCQCAQAIBEIAEAiEACARCAAAIlAAAASgQAAJAIBAEgEAgCQCAQAIBEIAEAiEACARCAAAMmHsgcAlttkWsTVzV3cjsaxt12L44Pd2KpWyh4LeGMCAZip2xtEu9OPwXD8c1ujXovzVjNODhslTga8NacYgGd1e4M4u7x+EgcREd+H4zi7vI5ub1DSZMB7EAhAMpkW0e70o3jmuYdt7U4/JtPn9gDWgUAAkqubu/TNwWNFRAyG47i6uXu/oVhZk2kRv/72e/zzX/+JX3/7XViuCNcgAMntaHYczLMfm8t1LKvLNwhAsrddW+h+bCbXsaw2gQAkxwe70ajXYtbNjJX48Snw+GD3PcdihbiOZfUJBCDZqlbivNWMiEiR8PD4vNW0HgIzuY5l9QkE4Fknh424OD2K/frT0wj79VpcnB45f8xfch3L6nORIjDTyWEjvjb3raTIq7mOZfUJBOAvbVUr8eXzx7LHYMU8XMfyfTh+9jqESvz4Nsp1LMvLKQYAFs51LKtPIADwJlzHstqcYgDgzbiOZXUJBADelOtYVpNTDABAIhAAgEQgAACJQAAAEoEAACQCAQBIBAIAkAgEACARCABAIhAAgEQgAACJQAAAEoEAACQCAQBIBAIAkAgEACARCABAIhAAgEQgAACJQAAAEoEAACQCAQBIBAIAkAgEACARCABAIhAAgEQgAACJQAAAEoEAACQCAQBIBAIAkAgEACARCABAIhAAgEQgAACJQAAAkg9lDwCLNJkWcXVzF7ejcext1+L4YDe2qpWyxwJYOQKBtdHtDaLd6cdgOP65rVGvxXmrGSeHjRInA1g9S3WKYTIt4tfffo9//us/8etvv8dkWpQ9Uukck5fp9gZxdnn9JA4iIr4Px3F2eR3d3qCkyQBW09J8g+DTX+aYvMxkWkS704/n0qmIiEpEtDv9+Nrcd7oB4IWW4hsEn/4yx+Tlrm7u0nF6rIiIwXAcVzd37zcUwIorPRD+7tNfxI9Pf5v01bpj8jq3o9lxMM9+ACxBIPj0lzkmr7O3XVvofgAsQSD49Jc5Jq9zfLAbjXotZl1dUIkf124cH+y+51gAK630QPDpL3NMXmerWonzVjMiIkXCw+PzVtMFigCvUHog+PSXOSavd3LYiIvTo9ivP42m/XotLk6P3PUB8Eql3+b48Onv7PI6KhFPLszb1E9/jsl8Tg4b8bW5byVFgAWoFEUx16Xw9/f3Ua/XYzgcxs7Ozv89iHv+M8cEgEV76fv30gRChHX0n+OYALBIL33/Lv0Uw2Nb1Up8+fyx7DGWimMCQBlKv0gRAFg+AgEASAQCAJAIBAAgEQgAQCIQAIBEIAAAiUAAABKBAAAkAgEASAQCAJAIBAAgEQgAQCIQAIBEIAAAiUAAABKBAAAkAgEASAQCAJAIBAAgEQgAQCIQAIBEIAAAiUAAABKBAAAkAgEASAQCAJAIBAAgEQgAQCIQAIBEIAAAiUAAABKBAAAkAgEASAQCAJAIBAAgEQgAQCIQAIBEIAAAiUAAABKBAAAkAgEASAQCAJAIBAAgEQgAQCIQAIBEIAAAiUAAABKBAAAkH8oeAID1NpkWcXVzF7ejcext1+L4YDe2qpWyx+JvCAQA3ky3N4h2px+D4fjntka9FuetZpwcNkqcjL/jFAMAb6LbG8TZ5fWTOIiI+D4cx9nldXR7g5Im4yUEAgALN5kW0e70o3jmuYdt7U4/JtPn9mAZCAQAFu7q5i59c/BYERGD4Tiubu7ebyheRSAAsHC3o9lxMM9+vD+BAMDC7W3XFrof708gALBwxwe70ajXYtbNjJX4cTfD8cHue47FKwgEABZuq1qJ81YzIiJFwsPj81bTeghLTCAA8CZODhtxcXoU+/WnpxH267W4OD2yDsKSs1ASAG/m5LARX5v7VlJcQQIBgDe1Va3El88fyx6DV3KKAQBIBAIAkAgEACARCABAIhAAgEQgAACJQAAAEoEAACQCAQBIBAIAkAgEACARCABAIhAAgEQgAACJQAAAEoEAACQCAQBIPpQ9APDDZFrE1c1d3I7Gsbddi+OD3diqVsoeC9hQAgGWQLc3iHanH4Ph+Oe2Rr0W561mnBw2SpwM2FROMUDJur1BnF1eP4mDiIjvw3GcXV5HtzcoaTJgkwkEKNFkWkS704/imecetrU7/ZhMn9sD4O0IBCjR1c1d+ubgsSIiBsNxXN3cvd9QACEQoFS3o9lxMM9+AIsiEKBEe9u1he4HsCgCAUp0fLAbjXotZt3MWIkfdzMcH+y+51gAAgHKtFWtxHmrGRGRIuHh8XmraT0E4N0JBCjZyWEjLk6PYr/+9DTCfr0WF6dH1kEASmGhJFgCJ4eN+Nrct5IisDQEAiyJrWolvnz+WPYYABHhFAMA8AyBAAAkAgEASAQCAJAIBAAgEQgAQCIQAIDEOgjPmEwLC9YAsNEEwv/o9gbR7vRjMPzz1+s26rU4bzUteQvAxnCK4ZFubxBnl9dP4iAi4vtwHGeX19HtDUqaDADel0D4w2RaRLvTj+KZ5x62tTv9mEyf2wMA1otA+MPVzV365uCxIiIGw3Fc3dy931AAUBKB8Ifb0ew4mGc/AFhlAuEPe9u1he4HAKtMIPzh+GA3GvVazLqZsRI/7mY4Pth9z7EAoBQC4Q9b1Uqct5oRESkSHh6ft5rWQwBgIwiER04OG3FxehT79aenEfbrtbg4PbIOAgAbw0JJ/+PksBFfm/tWUgRgowmEZ2xVK/Hl88eyxwCA0jjFAAAkAgEASAQCAJAIBAAgEQgAQCIQAIBEIAAAiUAAABKBAAAkAgEASAQCAJAIBAAgEQgAQCIQAIBEIAAAiUAAABKBAAAkAgEASAQCAJAIBAAgEQgAQCIQAIBEIAAAiUAAABKBAAAkAgEASD6UPQAArKLJtIirm7u4HY1jb7sWxwe7sVWtlD3WwggEAHilbm8Q7U4/BsPxz22Nei3OW804OWyUONniOMUAAK/Q7Q3i7PL6SRxERHwfjuPs8jq6vUFJky2WQACAF5pMi2h3+lE889zDtnanH5Ppc3usFoEAAC90dXOXvjl4rIiIwXAcVzd37zfUGxEIAPBCt6PZcTDPfstMIADAC+1t1xa63zITCADwQscHu9Go12LWzYyV+HE3w/HB7nuO9SYEAgC80Fa1EuetZkREioSHx+et5lqshyAQYElMpkX8+tvv8c9//Sd+/e33tbgKGtbRyWEjLk6PYr/+9DTCfr0WF6dHa7MOgoWSYAlswqIrsE5ODhvxtbm/1ispVoqimOtjyv39fdTr9RgOh7Gzs7PouWBjPCy68r//ER9+zKzTJxKgfC99/3aKAUq0SYuuAKtFIECJNmnRFWC1CAQo0SYtugKsFoEAJdqkRVeA1SIQoESbtOgKsFoEApRokxZdAVaLQICSbcqiK8BqsVASLIFNWHQFWC0CAZbEVrUSXz5/LHsMgIhwigEAeIZAAAASgQAAJAIBAEgEAgCQCAQAIBEIAEAiEACARCAAAMncKykWRREREff39wsbBgB4Ww/v2w/v47PMHQij0SgiIj59+jTvXwEAlGQ0GkW9Xp/5fKX4u4SYYTqdxrdv32J7ezsqFb9QBgBWQVEUMRqN4pdffolqdfaVBnMHAgCwvlykCAAkAgEASAQCAJAIBAAgEQgAQCIQAIBEIAAAiUAAABKBAAAkAgEASAQCAJAIBAAg+S+bYTk83UzkwgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Vertices coordinates (x,y)\n",
    "coords = np.array(list(csv.reader(open(\"./data/coords.csv\", \"r\"), delimiter=\",\"))).astype(\"float\")\n",
    "\n",
    "x1, y1 = coords.T\n",
    "plt.scatter(x1, y1)\n",
    "plt.gca().axes.get_xaxis().set_visible(False)\n",
    "plt.gca().axes.get_yaxis().set_visible(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generate an initial set of columns (clusters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 1 0 0 0 0 0 0 1 1 1 1 1 0 0]\n",
      " [0 0 1 0 1 0 1 0 1 1 1 1 1 0 0]\n",
      " [1 0 0 1 0 0 0 1 0 0 0 1 1 0 0]\n",
      " [1 0 1 0 0 0 0 0 0 0 0 0 1 0 0]\n",
      " [0 1 0 1 1 1 1 1 1 1 1 1 1 0 0]\n",
      " [1 1 1 0 0 1 1 0 0 0 1 0 0 0 0]\n",
      " [0 1 0 1 1 0 1 0 1 0 0 0 0 0 1]\n",
      " [0 0 1 1 1 0 0 1 1 0 0 0 0 0 0]\n",
      " [1 0 0 0 1 0 0 1 1 1 1 0 0 0 0]\n",
      " [0 0 1 0 1 0 0 1 0 1 1 1 0 1 0]\n",
      " [0 0 0 0 0 0 0 1 0 1 0 0 1 0 1]\n",
      " [1 1 0 1 1 1 1 1 1 1 0 1 0 0 1]\n",
      " [0 0 1 0 1 1 1 0 1 0 0 1 1 1 1]\n",
      " [1 1 1 1 1 1 0 0 1 1 1 0 1 0 1]\n",
      " [1 1 1 0 0 1 0 0 0 1 1 0 1 1 0]]\n",
      "[50, 92, 95, 68, 86, 55, 84, 57, 100, 81, 69, 83, 115, 23, 56]\n"
     ]
    }
   ],
   "source": [
    "################################################################################\n",
    "# Compute clusters' costs\n",
    "################################################################################\n",
    "def computeClusterCost(cluster):\n",
    "    # find best median\n",
    "    median = -1\n",
    "    minDist = float('inf')\n",
    "    for i in range(len(cluster)):\n",
    "        if (cluster[i] == 1):\n",
    "            dist = 0\n",
    "            for j in range(len(cluster)):\n",
    "                if (cluster[j] == 1):\n",
    "                    dist += distMtx[i][j]\n",
    "            if (dist < minDist):\n",
    "                minDist = dist\n",
    "                median = i\n",
    "    \n",
    "    return median, minDist\n",
    "    \n",
    "# Initial 2 * n random clusters\n",
    "initClusters = [[random.randint(0, 1) for i in range(n)] for j in range(n)]\n",
    "initClustersCosts = [computeClusterCost(cluster)[1] for cluster in initClusters]\n",
    "initClusters = np.array(initClusters)\n",
    "\n",
    "# Print initial clusters\n",
    "print(np.transpose(initClusters))\n",
    "print(initClustersCosts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Restricted master problem (RMP)\n",
    "\n",
    "$$\n",
    "  \\begin{align}\n",
    "    \\min & \\sum\\limits_{j = 1}^{m} c_j y_j \\\\\n",
    "    \\text{subject to} & \\nonumber \\\\\n",
    "    & \\sum\\limits_{j = 1}^{m} a_{ij} y_j \\geq 1, & \\forall i \\in V, \\\\\n",
    "    & \\sum\\limits_{j = 1}^{m} y_j = p, \\\\\n",
    "    & y_{j} \\in [0, 1], & j = 1, \\dots, m.\n",
    "  \\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Building our RMP model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model RMP\n",
    "rmpModel = gp.Model(\"rmp\")\n",
    "\n",
    "# Create decision variables y_j\n",
    "y = rmpModel.addVars(len(initClusters), obj=initClustersCosts, name=\"y\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Adding constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Covering constraints\n",
    "coverConstrs = rmpModel.addConstrs((gp.quicksum(initClusters[i][j] * y[i]\n",
    "                                                for i in range(len(initClusters))) >= 1 \n",
    "                                                for j in range(n)), name=\"c1_\")\n",
    "\n",
    "# p-medians constraint\n",
    "pConstr = rmpModel.addConstr(gp.quicksum(y[i] for i in range(len(initClusters))) == p, name=\"c2_\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subproblem (pricing)\n",
    "\n",
    "All columns $i$ $[a_{i1} \\cdots a_{in} | 1]$ satisfying the following inequality can be added to the RMP:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "  \\label{eq:rc}\n",
    "  \\min\\limits_{a_{ij} \\in \\{0, 1\\}}\\left(\\sum\\limits_{j = 1}^n (d_{ij} - \\pi_j)a_{ij}\\right) - \\pi_{n+1} < 0,\n",
    "\\end{equation}\n",
    "$$\n",
    "where coefficients $a{ij}$ are defined by\n",
    "$$\n",
    "\\begin{equation*}\n",
    "    a_{ij} = \n",
    "    \\begin{cases}\n",
    "        1, & \\text{if } d_{ij} - \\pi_j < 0 \\\\\n",
    "        0, & \\text{otherwise}.\n",
    "    \\end{cases}\n",
    "\\end{equation*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Retrieving the dual costs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################################################################\n",
    "# Retrieve duals from the RMP\n",
    "################################################################################\n",
    "def retrieveDualCosts(coverConstrs, pConstr):\n",
    "    duals = []\n",
    "    for constr in coverConstrs.select():\n",
    "        duals.append(constr.Pi)\n",
    "\n",
    "    duals.append(pConstr.Pi)\n",
    "    \n",
    "    return duals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generating columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################################################################\n",
    "# Function to generate a new column to the RMP\n",
    "################################################################################\n",
    "def generateColumn(i, duals):\n",
    "    rc = 0\n",
    "    cluster = np.zeros(n)\n",
    "    # RC of the n first constraints (cover constraints)\n",
    "    for j in range(n):\n",
    "        if (distMtx[i][j] - duals[j] < 0):\n",
    "            cluster[j] = 1\n",
    "            rc += distMtx[i][j] - duals[j]\n",
    "\n",
    "    # RC of the p-medians constraint\n",
    "    rc -= duals[n]\n",
    "    \n",
    "    # If negative RC, then compute cluster cost\n",
    "    clusterCost = float('inf')\n",
    "    if (rc < 0):\n",
    "        clusterCost = computeClusterCost(cluster)[1]\n",
    "    \n",
    "    return cluster, clusterCost\n",
    "\n",
    "################################################################################\n",
    "# Add a column to the RMP\n",
    "################################################################################\n",
    "def appendNewColumn(cluster, clusterCost):\n",
    "    \n",
    "    # Create a new variable related to the new column\n",
    "    new_y = rmpModel.addVar(obj=clusterCost, name=\"y_\" + str(rmpModel.NumVars))\n",
    "    \n",
    "    # Add the new y's cover constrs coefficient (first |V| constraints)\n",
    "    for i in range(n):\n",
    "        if (cluster[i] == 1):\n",
    "            rmpModel.chgCoeff(coverConstrs[i], new_y, 1)\n",
    "    \n",
    "    # Add the new y's p constraint (|V|th constraint)\n",
    "    rmpModel.chgCoeff(pConstr, new_y, 1)\n",
    "    \n",
    "    return\n",
    "\n",
    "################################################################################\n",
    "# Generates and adds new columns to the RMP\n",
    "################################################################################\n",
    "def addNewColumns(duals):\n",
    "    \n",
    "    nbFound = 0\n",
    "    for i in range(n):\n",
    "        cluster, clusterCost =  generateColumn(i, duals)\n",
    "        # If the column has negative RC, then add it to the RMP\n",
    "        if (clusterCost < float('inf')):\n",
    "            appendNewColumn(cluster, clusterCost)\n",
    "            nbFound += 1\n",
    "\n",
    "    print(\" added \" + str(nbFound) + \" newcolumns\", end='')\n",
    "    return nbFound"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Column generation main routine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMP initial ObjVal: 235.00\n",
      "\n",
      "CG(0): added 15 newcolumns, RMP ObjVal: 197.12\n",
      "CG(1): added 15 newcolumns, RMP ObjVal: 160.00\n",
      "CG(2): added 15 newcolumns, RMP ObjVal: 120.43\n",
      "CG(3): added 15 newcolumns, RMP ObjVal: 102.57\n",
      "CG(4): added 15 newcolumns, RMP ObjVal: 93.67\n",
      "CG(5): added 15 newcolumns, RMP ObjVal: 75.50\n",
      "CG(6): added 15 newcolumns, RMP ObjVal: 71.00\n",
      "CG(7): added 6 newcolumns, RMP ObjVal: 69.00\n",
      "CG(8): added 6 newcolumns, RMP ObjVal: 68.00\n",
      "CG(9): added 1 newcolumns, RMP ObjVal: 68.00\n",
      "CG(10): added 2 newcolumns, RMP ObjVal: 67.00\n",
      "CG(11): added 1 newcolumns, RMP ObjVal: 67.00\n",
      "CG(12): added 1 newcolumns, RMP ObjVal: 67.00\n",
      "CG(13): added 3 newcolumns, RMP ObjVal: 67.00\n",
      "CG(14): added 0 newcolumns\n",
      "\n",
      "RMP # of variables: 140\n"
     ]
    }
   ],
   "source": [
    "# Solve the model for the first time\n",
    "rmpModel.Params.OutputFlag = 0\n",
    "rmpModel.optimize()\n",
    "print(\"RMP initial ObjVal: %.2f\\n\" % rmpModel.ObjVal)\n",
    "\n",
    "duals = retrieveDualCosts(coverConstrs, pConstr)\n",
    "\n",
    "i = 0;\n",
    "maxIter = 1000\n",
    "print(f'CG({i}):', end='')\n",
    "while (i < maxIter and addNewColumns(duals)):\n",
    "    rmpModel.optimize()\n",
    "    duals = retrieveDualCosts(coverConstrs, pConstr)\n",
    "    i += 1\n",
    "    print(\", RMP ObjVal: %.2f\" % rmpModel.ObjVal)\n",
    "    print(f'CG({i}):', end='')\n",
    "    \n",
    "print(f'\\n\\nRMP # of variables: {rmpModel.NumVars}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Imposing integrality on y variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set parameter OutputFlag to value 1\n"
     ]
    }
   ],
   "source": [
    "rmpModel.Params.OutputFlag = 1\n",
    "\n",
    "# Change y vars to binary to solve the original PMP\n",
    "for var in rmpModel.getVars():\n",
    "    var.VType = GRB.BINARY\n",
    "    \n",
    "# Change the sense of cover constraints to partitioning constraints (<= to =)\n",
    "for constr in rmpModel.getConstrs():\n",
    "    constr.Sense = '='"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving the set partitioning PMP\n",
    "\n",
    "$$\n",
    "  \\begin{align}\n",
    "    \\min & \\sum\\limits_{j = 1}^{2^n} c_j y_j \\\\\n",
    "    \\text{subject to} & \\nonumber \\\\\n",
    "    & \\sum\\limits_{j = 1}^{2^n} a_{ij} y_j = 1, & \\forall i \\in V, \\\\\n",
    "    & \\sum\\limits_{j = 1}^{2^n} y_j = p, \\\\\n",
    "    & y_{j} \\in \\{0, 1\\}, & j = 1, \\dots, 2^n.\n",
    "  \\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (linux64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 16 rows, 140 columns and 644 nonzeros\n",
      "Model fingerprint: 0xbcfa6acb\n",
      "Variable types: 0 continuous, 140 integer (140 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [3e+00, 1e+02]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [1e+00, 5e+00]\n",
      "Presolve removed 0 rows and 42 columns\n",
      "Presolve time: 0.01s\n",
      "Presolved: 16 rows, 98 columns, 394 nonzeros\n",
      "Variable types: 0 continuous, 98 integer (98 binary)\n",
      "\n",
      "Root relaxation: objective 6.700000e+01, 24 iterations, 0.00 seconds (0.00 work units)\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "*    0     0               0      67.0000000   67.00000  0.00%     -    0s\n",
      "\n",
      "Explored 1 nodes (24 simplex iterations) in 0.03 seconds (0.00 work units)\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 1: 67 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 6.700000000000e+01, best bound 6.700000000000e+01, gap 0.0000%\n"
     ]
    }
   ],
   "source": [
    "# Solve it again, but now as an integer linear program (ILP)\n",
    "rmpModel.optimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABKC0lEQVR4nO3deVhUZcMG8HtmkEUE3BFZ5BWXMg0VEBVxCzMtN8TMXTN3E0HNrNQyt1whozTStBIxEbVyzYXCNFkUS19TcQFEQNFkU1mG8/3hC58o4AzOzDPD3L/rmj84nBluOOLcnOc5z5FJkiSBiIiISBC56ABERERk3FhGiIiISCiWESIiIhKKZYSIiIiEYhkhIiIioVhGiIiISCiWESIiIhKKZYSIiIiEMhEdQBXFxcW4efMmrKysIJPJRMchIiIiFUiShJycHDRu3BhyecXnPwyijNy8eROOjo6iYxAREVEVpKSkwMHBocLPG0QZsbKyAvDom7G2thachoiIiFSRnZ0NR0fH0vfxihhEGSkZmrG2tmYZISIiMjDPmmLBCaxEREQkFMsIERERCcUyQkREREKxjBAREZFQLCNEREQkFMsIERERCcUyQkREREKxjBAREZFQLCM6kpubi0GDBuHll1/GoEGDkJubKzpStaFUKhEVFYVt27YhKioKSqVSdCQiIlKDQazAaug6dOiA2NjY0o///vtvWFlZwcPDAzExMQKTGb7IyEj4+/vjxo0bpdscHBwQHBwMX19fgcmIiEhVPDOiZU8WkcfFxsaiQ4cOOk5UfURGRsLPz69MEQGA1NRU+Pn5ITIyUlAyIiJSh0ySJEl0iGfJzs6GjY0NsrKyDOreNLm5uc+8ORAApKeno1atWjpIVH0olUq0atUKqamp5X5eJpPBwcEB165dg0Kh0HE6IiICVH//5jCNFo0aNUql/Ro1aqTlJMZHkiSkpKQgOjoa3bt3Fx2HiIgqwWEaLbpy5YroCEYvLS1NdAQiInoGnhnRIhcXF/z999/P3O/111/H9u3bdZCo+vj999/Rt2/fZ+5nZ2engzRERPQ8OGdEi1SdM5KTk8M5I2pSKpVwdnZGamoqyvsnzDkjRETiqfr+zWEaLapVqxY8PDwq3cfDw4NFpAoUCgWCg4MBPCoejyv5OCgoiEWEiMgAsIxoWUxMDNq1a1fu57jOyPPx9fVFREQE7O3ty2x3cHBAREQE1xkhIjIQHKbRgRMnTsDLywvm5uZo3rw5XFxc8P333/OMiIYolUpER0cjLS0NdnZ28Pb25hkRIiI9wEt79cjFixcBAF26dMGvv/4qOE31o1AoePkuEZEB4zCNDpSUkZYtWwpOQkREpH9YRnSAZYSIiKhiLCM6wDJCRERUMZYRLSsqKkJiYiIAlhEiIqLysIxo2fXr11FYWAgLCws4OjqKjkNERKR3WEa0rGSIpnnz5pDL+eMmIiJ6Et8dtYzzRYiIiCrHMqJlLCNERESVYxnRMpYRIiKiyrGMaBnLCBERUeVYRrQoOzsb6enpAFhGiIiIKsJ702hRyVmRRo0aGeQN/ohIf/CGkKQNDx48wJw5c3D58mU0b94cK1euhIWFhc5z8MyIFnGIhog0ITIyEs7OzujRoweGDx+OHj16wNnZGZGRkaKjkQEbOHAgatasiZCQEBw6dAghISGoWbMmBg4cqPMsLCNaxDJCRM8rMjISfn5+uHHjRpntqamp8PPzYyGhKhk4cCD27NlT7uf27Nmj80LCYRotYhkhouehVCrh7+8PSZKe+pwkSZDJZPD394ePjw+HbEhlDx48qLCIlNizZw8ePHigsyEblhEtYhkhoucRHR391BmRx0mShBs3bsDGxkaHqchYzJkzB1988YVOvhaHabSkuLgYly9fBsAyQkRVk5aWJjoCGbGS9zBd4JkRLUlJScGDBw9Qo0YNODs7i45DRAbIzs5Opf327duHrl27ajkNVRcBAQEIDQ195n7NmzfXQZpHZFJ5g5F6Jjs7GzY2NsjKyjKYS2QPHTqE3r1744UXXsCFCxdExyEiA6RUKuHs7IzU1NRy543IZDI4ODjg2rVrnDNCKnvw4AFq1qz5zP3u37//3HNGVH3/5jCNlnC+CBE9L4VCgeDgYACPisfjSj4OCgpiESG1WFhYYMCAAZXuM2DAAJ2uN8IyoiUsI0SkCb6+voiIiIC9vX2Z7Q4ODoiIiICvr6+gZGTIdu/eXWEhGTBgAHbv3q3TPJwzoiUsI0SkKb6+vhgwYABXYCWN2r17t96swMoyoiUsI0SkSQqFAt27dxcdg6oZCwsLnV2+WxkO02hBXl4eUlJSALCMEBERPQvLiBaUXJtdt25d1K9fX3AaIiIi/cYyogUcoiEiIlIdy4gWsIwQERGpjmVECy5dugSAZYSIiEgVLCNawDMjREREqmMZ0TBJklhGiIiI1MAyomHp6enIycmBXC6Hi4uL6DhERER6j2VEw0rOivznP/+BmZmZ4DRERET6j2VEwzhEQ0REpB6WEQ1jGSEiIlIPy4iGsYwQERGph2VEw1hGiIiI1MMyokH5+fm4du0aAJYRIiIiVbGMaNCVK1dQXFwMKysrNGrUSHQcIiIig8AyokGPD9HIZDLBaYiIiAxDlcpISEgInJ2dYW5uDk9PT8TExFS6f1BQEFq2bAkLCws4OjoiICAADx8+rFJgfcb5IkREROpTu4xs374dgYGBWLhwIU6fPg1XV1f07t0bt27dKnf/sLAwvP/++1i4cCEuXLiAjRs3Yvv27fjggw+eO7y+YRkhIiJSn9plZM2aNZgwYQLGjRuHVq1aYf369ahZsyY2bdpU7v4nTpyAl5cXhg8fDmdnZ7z66qsYNmzYM8+mGCKWESIiIvWpVUYKCgoQHx8PHx+f/38BuRw+Pj44efJkuc/p3Lkz4uPjS8vH1atXsW/fPvTt27fCr5Ofn4/s7OwyD0PAMkJERKQ+E3V2zszMhFKphK2tbZnttra2+Oeff8p9zvDhw5GZmYkuXbpAkiQUFRVh8uTJlQ7TLFu2DJ988ok60YTLzMzE3bt3AQDNmzcXnIaIiMhwaP1qmqioKCxduhRffvklTp8+jcjISOzduxeffvpphc+ZN28esrKySh8pKSnajvncSs6KODk5oWbNmoLTEBERGQ61zozUr18fCoUCGRkZZbZnZGRUuK7G/PnzMWrUKLzzzjsAgDZt2iAvLw8TJ07Ehx9+CLn86T5kZmZmcHe85RANERFR1ah1ZsTU1BRubm44cuRI6bbi4mIcOXIEnTp1Kvc59+/ff6pwKBQKAIAkSerm1VssI0RERFWj1pkRAAgMDMSYMWPg7u6ODh06ICgoCHl5eRg3bhwAYPTo0bC3t8eyZcsAAP369cOaNWvQrl07eHp6IjExEfPnz0e/fv1KS0l1wDJCRERUNWqXkaFDh+L27dtYsGAB0tPT0bZtWxw4cKB0UmtycnKZMyEfffQRZDIZPvroI6SmpqJBgwbo168flixZornvQg+wjBAREVWNTDKAsZLs7GzY2NggKysL1tbWouM8paioCDVr1kRhYSGSkpLg5OQkOhIREZFwqr5/8940GnDt2jUUFhbCwsICDg4OouMQEREZFJYRDSgZomnRokW5VwcRERFRxfjOqQGcL0JERFR1LCMawDJCRERUdSwjGsAyQkREVHUsIxrAMkJERFR1LCPPKSsrq3R5/BYtWghOQ0REZHhYRp5TyVkROzs7vVwDhYiISN+xjDwnDtEQERE9H7WXg6eyWEaqRqlUIjo6GmlpabCzs4O3t3e1ulcRERGpjmXkObGMqC8yMhL+/v64ceNG6TYHBwcEBwfD19dXYDIiIhKBwzTPiWVEPZGRkfDz8ytTRAAgNTUVfn5+iIyMFJSMiIhEYRl5DsXFxbh8+TIAlhFVKJVK+Pv7o7x7M5ZsmzlzJpRKpa6jERGRQCwjzyE5ORkPHz6EqakpnJ2dRcfRe9HR0U+dEXmcJElISUlBdHS0DlMREZFoLCPPoWSIplmzZpx8qYK0tDSN7kdERNUDy8hz4HwR9djZ2Wl0PyIiqh5YRp4Dy4h6vL294eDgAJlMVu7nZTIZHB0d4e3treNkREQkEsvIc2AZUY9CoUBwcDAAVFhIgoKCOORFRGRkWEaeA8uI+nx9fREREQF7e/unPrdhwwauM0JEZIRYRqooLy+v9MoQlhH1+Pr64vr16zh27BjCwsLw8ssvA0DpZdJERGRcZFJ5iz7omezsbNjY2CArK0tvbkZ35swZtG/fHvXr18ft27dFxzFoe/fuxRtvvAFra2skJyfDxsZGdCQiItIAVd+/eWakijhEozl9+vTBSy+9hOzsbHz99dei4xARkY6xjFQRy4jmyOVyzJ49G8CjCaz5+fmCExERkS6xjFQRy4hmDR8+HI0bN8bNmzcRFhYmOg4REekQy0gVsYxolqmpKQICAgAAK1euRHFxseBERESkKywjVSBJEi5dugSAZUSTJk6cCGtra1y4cAH79u0THYeIiHSEZaQKbt68idzcXCgUCjRt2lR0nGrD2toaU6ZMAQCsWLFCcBoiItIVlpEqKBmiadq0KUxNTQWnqV5mzJiBGjVqIDo6GidPnhQdh4iIdIBlpApKykiLFi0EJ6l+GjdujFGjRgF4NHeEiIiqP5aRKuDkVe0qucx39+7dpXNziIio+mIZqQKWEe168cUX0b9/f0iShNWrV4uOQ0REWsYyUgUsI9o3Z84cAMCWLVuQnp4uOA0REWkTy4iaHj58iOvXrwNgGdEmLy8vdOrUCfn5+Vi3bp3oOEREpEUsI2pKTEyEJEmwtraGra2t6DjVlkwmw3vvvQcA+PLLL5GTkyM4ERERaQvLiJoeH6KRyWSC01Rv/fv3R4sWLXDv3j1s3LhRdBwiItISlhE1cb6I7sjl8tK5I2vWrEFhYaHgREREpA0sI2piGdGtkSNHwtbWFikpKdi+fbvoOEREpAUsI2piGdEtc3Nz+Pv7A3i0RLwkSYITERGRprGMqEGSJJYRASZPnoxatWrh77//xsGDB0XHISIiDWMZUUNmZibu3bsHmUyG5s2bi45jNOrUqYOJEycC4BLxRETVEcuIGkrOijg5OcHCwkJwGuMyc+ZMmJiY4OjRo4iLixMdh4iINIhlRA0cohHH0dERw4YNA8CzI0RE1Q3LiBpYRsQqucw3IiICV65cEZyGiIg0hWVEDSwjYrVp0wZ9+vRBcXEx1qxZIzoOERFpCMuIGlhGxCtZIv7bb7/F7du3BachIiJNYBlRUWFhYenQAMuION26dYO7uzsePHiAkJAQ0XGIiEgDWEZUdO3aNRQVFaFmzZqwt7cXHcdoPX4DvS+++AJ5eXmCExER0fNiGVFRyRBNixYtIJfzxyaSr68vmjZtijt37uDbb78VHYeIiJ4T31VVxPki+kOhUGD27NkAgNWrV6OoqEhwIiIieh4sIypiGdEvY8eORYMGDXD9+nXs3LlTdBwiInoOLCMqYhnRLxYWFpg+fToA3kCPiMjQsYyoiGVE/0ybNg01a9bE6dOncfToUdFxiIioilhGVHDv3j3cunULwKMJrKQf6tWrh/HjxwN4dHaEiIgME8uICkrOijRu3BhWVlaC09DjAgMDoVAocOjQISQkJIiOQ0REVcAyogIO0egvZ2dnDBkyBACwatUqwWmIiKgqWEZUwDKi30puoBceHo6kpCTBaYiISF0sIypgGdFv7du3h4+PD5RKJdauXSs6DhERqYllRAUsI/qvZIn40NBQ3L17V3AaIiJSB8vIMyiVSly+fBkAy4g+8/HxQdu2bXH//n189dVXouMQEZEaWEaeITk5Gfn5+TAzM0OTJk1Ex6EKyGSy0rkjn3/+OR48eCA4ERERqYpl5BlKhmiaNWsGhUIhOA1VZsiQIWjSpAlu3bqF7777TnQcIiJSEcvIM3C+iOGoUaMGAgMDATy6zFepVApOREREqmAZeQaWEcMyfvx41K1bF4mJidi9e7foOEREpAKWkWdgGTEslpaWmDp1KgDeQI+IyFCwjDwDy4jheffdd2FmZoaYmBhER0eLjkNqUCqViIqKwrZt2xAVFcWhNiIjUaUyEhISAmdnZ5ibm8PT0xMxMTGV7n/v3j1MmzYNdnZ2MDMzQ4sWLbBv374qBdal3NxcpKamAmAZMSQNGzbEuHHjAPAGeoYkMjISzs7O6NGjB4YPH44ePXrA2dkZkZGRoqMRkZapXUa2b9+OwMBALFy4EKdPn4arqyt69+5delfbJxUUFKBXr164fv06IiIicPHiRYSGhsLe3v65w2vbpUuXAAANGjRAnTp1BKchdcyaNQsymQx79+7FuXPnRMehZ4iMjISfnx9u3LhRZntqair8/PxYSIiqObXLyJo1azBhwgSMGzcOrVq1wvr161GzZk1s2rSp3P03bdqEu3fvYvfu3fDy8oKzszO6desGV1fX5w6vbRyiMVzNmjXD4MGDAfAGevpOqVTC39+/3Pk9JdtmzpzJIRuiakytMlJQUID4+Hj4+Pj8/wvI5fDx8cHJkyfLfc5PP/2ETp06Ydq0abC1tUXr1q2xdOnSSv9jyc/PR3Z2dpmHCCwjhq1kEbSwsLCn/uIm/REdHV3p8ZEkCSkpKZz/QyrhvCPDpFYZyczMhFKphK2tbZnttra2SE9PL/c5V69eRUREBJRKJfbt24f58+dj9erVWLx4cYVfZ9myZbCxsSl9ODo6qhNTY1hGDFuHDh3QrVs3FBYWIjg4WHQcqkBaWppG9yPjxXlHhkvrV9MUFxejYcOG+Prrr+Hm5oahQ4fiww8/xPr16yt8zrx585CVlVX6SElJ0XbMcrGMGL6SG+ht2LAB9+7dExuGymVnZ6fR/cg4cd6RYVOrjNSvXx8KhQIZGRlltmdkZKBRo0blPsfOzg4tWrQos5T6iy++iPT0dBQUFJT7HDMzM1hbW5d56JokSaUTWFlGDFefPn3QunVr5OTkYMOGDaLjUDm8vb3h4OAAmUxW4T62trbw9vbWYSoyJJx3ZPjUKiOmpqZwc3PDkSNHSrcVFxfjyJEj6NSpU7nP8fLyQmJiIoqLi0u3Xbp0CXZ2djA1Na1ibO1LTU1FXl4eTExM0LRpU9FxqIpkMhlmz54NAAgKCkJ+fr7gRPQkhUJROoxWUSG5e/cutm3bpstYZEA478jwqT1MExgYiNDQUGzZsgUXLlzAlClTkJeXV7quw+jRozFv3rzS/adMmYK7d+/C398fly5dwt69e7F06VJMmzZNc9+FFpQM0TRt2hQ1atQQnIaex7Bhw2Bvb4/09HRs3bpVdBwqh6+vLyIiIp665N/e3h7t2rVDYWEhRo0ahRkzZqCwsFBQStJXnHdk+NQuI0OHDsWqVauwYMECtG3bFgkJCThw4EDppNbk5OQyB9zR0REHDx5EbGwsXn75ZcyYMQP+/v54//33NfddaAHni1QfpqamCAgIAACsXLmyzFk60h++vr64fv06jh07hrCwMBw7dgxJSUmIi4vDggULAADr1q1Dz549K5wwT8aJ844Mn0wygJt3ZGdnw8bGBllZWTqbP+Lv74/PP/8cs2fPxsqVK3XyNUl7srOz4eTkhKysLOzZswf9+/cXHYnU9NNPP2HUqFHIzs5G48aNERERUeHwMBkXpVIJZ2dnpKamljtvRCaTwcHBAdeuXSszf5G0T9X3b96bpgI8M1K9WFtbY8qUKQC4RLyh6t+/P2JiYtCqVSvcvHkT3bp1w/r163kzRKp03lHJx0FBQSwieoxlpAIsI9XPjBkzYGpqij/++AMnTpwQHYeqoGXLlvjzzz/h5+eHwsJCTJkyBePHj8fDhw9FRyPBSuYdNW7cuMx2BwcHREREwNfXV1AyUgXLSDkePHiApKQkACwj1YmdnR1GjRoFABx6M2BWVlb48ccf8dlnn0Eul+Pbb7+Ft7c3kpOTRUcjwXx9ffHf//639ON9+/bh2rVrLCIGgGWkHImJiZAkCbVr10aDBg1ExyENKrnMd8+ePfjnn38Ep6GqkslkeO+993Dw4EHUq1cPcXFxcHNzw9GjR0VHI8EeH4rp2rUrh2YMBMtIOR4foqlsISYyPC+88AIGDBgASZKwevVq0XHoOfn4+CAuLg7t27dHZmYmevXqhVWrVnEeCZGBYRkpB+eLVG8lN9D77rvvuO5ANeDs7Izjx49jzJgxKC4uxpw5c/DWW28hNzdXdDQiUhHLSDlYRqo3Ly8vdO7cGQUFBfj8889FxyENsLCwwLfffouQkBCYmJjgxx9/RKdOnXD58mXR0YhIBSwj5WAZqf5KbqD31VdfIScnR3Aa0gSZTIapU6ciKioKjRo1wrlz5+Dh4YFffvlFdDQiegaWkSdIksQyYgT69euHli1bIisrC6GhoaLjkAZ5eXnh9OnT8PLyQlZWFvr164ePP/6YK+8S6TGWkSfcunULWVlZkMlkaNasmeg4pCVyubx07sjatWsrvIM0GSY7OzscPXq09B5Yn3zyCQYMGIB79+6JDUZE5WIZeULJWRFnZ2eYm5sLTkPaNHLkSDRq1Ag3btxAeHi46DikYaampvjiiy+wefNmmJub45dffoGHhwfOnTsnOhoRPYFl5AkcojEeZmZm8Pf3B/BoETReDlo9jRkzBn/88QeaNGmCxMREeHp64scffxQdi4gewzLyBJYR4zJ58mTUqlUL586dw4EDB0THIS1p37494uLi4OPjg/v372Po0KGYM2cOioqKREcjIrCMPIVlxLjUrl0bkyZNAsAb6FV39evXx/79+zF37lwAwKpVq9C7d2/cvn1bcDIiYhl5QkkZadGiheAkpCv+/v4wMTFBVFQUYmJiRMchLTIxMcHy5cuxY8cOWFpa4ujRo3Bzc0NcXJzoaERGjWXkMQUFBbh69SoAnhkxJo6Ojhg+fDgA3kDPWPj5+SEmJgbNmzdHSkoKunTpgk2bNomORWS0WEYec/XqVSiVSlhaWsLe3l50HNKhkhvoRUZGIjExUXAa0oVWrVohNjYW/fv3R35+PsaPH48pU6bwMm8iAVhGHvP4EA1vkGdc2rRpg759+6K4uBhr1qwRHYd0xMbGBrt27cKnn34KmUyG9evXo1u3bkhNTRUdjciosIw8hpNXjVvJEvHffvstbt26JTgN6YpcLsdHH32EX375BbVr18aff/4JNzc3REdHi45GZDRYRh7DMmLcunbtCg8PDzx8+BBffPGF6DikY3379kVcXBzatGmDjIwM9OzZE+vWreP6M0Q6wDLyGJYR4yaTyUrPjoSEhCAvL09wItI1FxcXnDx5Em+99RaKioowY8YMjBkzBvfv3xcdjahaYxl5DMsIDRo0CC4uLrh79y6vrjBSlpaWCAsLw5o1a6BQKPD999/Dy8sL165dEx2NqNpiGfmfu3fvIjMzEwDXGDFmCoWi9Mqa1atXc4VOIyWTyRAQEIDDhw+jQYMGSEhIgLu7Ow4dOiQ6GlG1xDLyPyVnRezt7VGrVi3BaUikMWPGoEGDBkhKSsKOHTtExyGBunfvjvj4eHh4eODu3bt47bXXsGzZMs4jIdIwlpH/uXTpEgAO0RBgYWGBd999F8CjJeL5xmPcHB0d8fvvv+Odd96BJEn44IMPMHjwYGRnZ4uORlRtsIz8D+eL0OOmTp2KmjVrIiEhAYcPHxYdhwQzNzdHaGgovv76a5iammLXrl3w9PTEP//8IzoaUbXAMvI/LCP0uHr16uGdd94BwCXi6f9NmDABv/32G+zt7fHPP/+gQ4cO2L17t+hYRAaPZeR/WEboSQEBAVAoFPj1119x5swZ0XFIT3Ts2BHx8fHo2rUrcnJyMGjQIHz44YdQKpWioxEZLJYRAEqlsvR+JCwjVMLZ2RlvvvkmAJ4dobJsbW1x+PBhzJw5EwCwdOlSvP7667h7967YYEQGimUEQFJSEvLz82FmZgYnJyfRcUiPzJkzBwDw448/cp0JKqNGjRpYu3YtfvjhB1hYWODgwYNwd3fH2bNnRUcjMjgsI/j/IZrmzZtDoVAITkP6pF27dujVqxeUSiXWrl0rOg7poREjRuDkyZP4z3/+g2vXrqFTp07YunWr6FhEBoVlBJwvQpUrWSJ+48aNuHPnjuA0pI9cXV0RFxeH3r1748GDBxg5ciRmzpyJwsJC0dGIDALLCFhGqHKvvPIK2rVrh/v37+PLL78UHYf0VN26dbF37158+OGHAIDg4GD4+PggIyNDcDIi/ccyApYRqpxMJiudO/L555/jwYMHghORvlIoFFi8eDF27doFKysr/P7772jfvj3+/PNP0dGI9BrLCFhG6NmGDBmCJk2aIDMzE5s3bxYdh/TcwIEDERsbixdffBE3b95E165dsWHDBq7mS1QBoy8jOTk5uHnzJgCWEaqYiYkJZs2aBeDRDfS4pgQ9S8uWLXHq1Cn4+vqisLAQkydPxoQJE/Dw4UPR0Yj0jtGXkZJ70jRs2BC1a9cWG4b02ttvv426deviypUr2LVrl+g4ZACsrKwQERGBZcuWQS6XY+PGjejatStSUlJERyPSK0ZfRjhEQ6qytLTEtGnTAPAGeqQ6mUyG999/H/v370fdunURGxsLNzc3HDt2THQ0Ir3BMsIyQmqYPn06zM3NERsbi99++010HDIgr776KuLi4tC2bVvcvn0bvXr1wpo1a1hqicAywjJCamnYsCHGjRsH4NHZESJ1/Oc//8Eff/yBUaNGQalUYtasWRg+fDjy8vJERyMSimWEZYTUFBgYCLlcjv379+Ott95C7969MX36dF7ySyqpWbMmtmzZgnXr1sHExATh4eHo2LFj6f2xiIyRTDKAc4TZ2dmwsbFBVlYWrK2tNfa6xcXFsLKywv3793Hx4kW0aNFCY69N1Zu9vX3pVViPGzBgAG8pTyqLjo7GkCFDkJGRgdq1a2Pr1q3o27ev6FgGLS8vD7Vq1QIA5ObmwtLSUnAi46bq+7fRnhlRKpWIiIjA/fv3IZfLeYM8PPqZREVFYdu2bYiKiuLlqxUYOHBguUUEAPbs2YOBAwfqNhAZLG9vb5w+fRqdOnXCvXv38MYbb+DTTz9FcXGx6GhEOmWUZSQyMhLOzs4YOnQogEdnSJo3b47IyEjBycQp+Zn06NEDw4cPR48ePeDs7GzUP5PyPHjwAHv27Kl0nz179nDIhlTWuHFjREVFYcqUKZAkCQsWLMDAgQORlZUlOhqRzhhdGYmMjISfnx9u3LhRZntqair8/PyM8s2XPxPVlSwLr6n9iADA1NQUX375JTZt2gQzMzP8/PPP8PDwwPnz50VHI9IJoyojSqUS/v7+5V5KV7Jt5syZRjU8wZ+Jei5fvqzR/YgeN27cOBw/fhyOjo64fPkyPD09sWPHDtGxiLTOqMpIdHT0U3/9P06SJKSkpCA6OlqHqcTiz0Q9zZs31+h+RE9yd3dHfHw8evbsiby8PLz55puYO3cuioqKREcj0hqjKiNpaWka3a864M9EPStXrtTofkTladCgAQ4ePIjZs2cDeLSmzWuvvYbMzEzByYi0w6jKiJ2dnUb3qw74M1GPhYUFBgwYUOk+AwYMgIWFhY4SUXVlYmKClStXYvv27bC0tMSRI0fg5uaG+Ph40dGINM6oyoi3tzccHBwgk8nK/bxMJoOjoyO8vb11nEycZ/1MgEeT65o0aaLDVPpt9+7dFRYSrjNCmvbmm2/izz//RLNmzZCcnAwvLy9s2bJFdCwijTKqMqJQKBAcHAwAT735lnwcFBQEhUKh82yiVPYzKVFQUAA3NzfeqfYxu3fvxv379zFt2jS8+uqrmDZtGu7fv88iQlrRunVrxMbG4o033kB+fj7Gjh2LadOmoaCgQHQ0Is2QDEBWVpYEQMrKytLI6+3cuVNycHCQAJQ+HB0dpZ07d2rk9Q1RRT+TkJAQycPDo3TblClTpPv374uOS2SUlEql9PHHH5f+Pnbu3FlKTU0VHUuv5Obmlv58cnNzRccxeqq+fxvtcvBKpRLR0dFIS0uDnZ0dvL29jeqMSHkq+pkUFBRg/vz5pTeGa9OmDcLDw9GqVSvBiYmM0y+//IKRI0ciKysLjRo1wo4dO9ClSxfRsfQCl4PXL6q+fxttGSH1HTp0CKNGjcKtW7dgYWGBoKAgTJgwodL5JkSkHYmJiRg0aBDOnTsHExMTBAUFYerUqUb/+8gyol94bxrSuFdffRV//fUXXn31VTx48ACTJk3C0KFDce/ePdHRiIxOs2bNcPLkSbz55psoKirC9OnTMXbsWN6KgAwSywipxdbWFvv378eKFStgYmKCHTt2oG3btjhx4oToaERGp1atWggPD8eqVasgl8vx3XffwcvLC9evXxcdjUgtLCOkNrlcjjlz5uDEiRNo2rQpkpKS0LVrVyxZsoTLxhPpmEwmw6xZs/Drr7+ifv36OHPmDNzd3fHrr7+KjkakMpYRqjIPDw+cOXMGw4cPh1KpxEcffYRevXrh5s2boqMRGZ2ePXsiPj4e7u7uuHPnDl577TV89tln5d53ikjfsIzQc7G2tsYPP/yAzZs3w9LSEseOHYOrqyv27t0rOhqR0XFyckJ0dDTGjRuH4uJivP/++xgyZAhycnJERyOqFMsIPTeZTIYxY8YgPj4e7dq1Q2ZmJt544w0EBAQgPz9fdDwio2Jubo6NGzdi/fr1qFGjBnbu3AlPT09cvHhRdDSiCrGMkMa0bNkSJ0+ehL+/P4BHq9l26tQJly5dEpyMyLjIZDJMmjQJv/32Gxo3bowLFy6gQ4cO2LNnj+hoROViGSGNMjMzQ1BQEH7++WfUq1cPZ86cQfv27bFlyxaOXRPpWKdOnRAfHw9vb29kZ2dj4MCBmD9/Pieak95hGSGteOONN3D27Fl0794deXl5GDt2LEaNGoXs7GzR0YiMSqNGjXDkyBHMmDEDALB48WL069cP//77r+BkRP+PZYS0xt7eHocPH8bixYuhUCiwdetWtG/fHrGxsaKjERmVGjVqIDg4GN9//z3Mzc2xf/9+uLu746+//hIdjQgAywhpmUKhwIcffojff/8dTk5OuHLlCjp37oxVq1ahuLhYdDwiozJy5EicOHECzs7OuHr1Kjp27Iht27aJjkXEMkK60blzZyQkJGDw4MEoKirCnDlz8PrrryMjI0N0NCKj0q5dO8TFxZXe1mH48OEIDAxEYWGh6GhkxKpURkJCQuDs7Axzc3N4enoiJiZGpeeFh4dDJpNh4MCBVfmyZODq1KmDHTt2YMOGDTA3N8eBAwfg6urKlSKJdKxevXrYt28fPvjgAwDA2rVr0atXL/5xQMKoXUa2b9+OwMBALFy4EKdPn4arqyt69+6NW7duVfq869evY/bs2fD29q5yWDJ8MpkMEydORFxcHF566SVkZGTg1Vdfxdy5c/mXGZEOKRQKLFmyBJGRkahVqxZ+++03uLm54dSpU6KjkRFSu4ysWbMGEyZMwLhx49CqVSusX78eNWvWxKZNmyp8jlKpxIgRI/DJJ5+gadOmzxWYqoeXXnoJsbGxmDx5MgBgxYoV6NKlC65evSo4GZFxGTRoEGJiYtCyZUukpqaia9euCA0NFR2LjIxaZaSgoADx8fHw8fH5/xeQy+Hj44OTJ09W+LxFixahYcOGGD9+vEpfJz8/H9nZ2WUeVP1YWFjgq6++QkREBGrXro2YmBi0a9cO4eHhoqMRGZUXX3wRMTExGDhwIAoKCjBx4kRMnDiRKyiTzqhVRjIzM6FUKmFra1tmu62tLdLT08t9zvHjx7Fx40a1mvayZctgY2NT+nB0dFQnJhmYwYMH4+zZs/Dy8kJ2djaGDRuGt99+G3l5eaKjERkNa2tr7Ny5E0uWLIFMJkNoaCi6du2KGzduiI5GRkCrV9Pk5ORg1KhRCA0NRf369VV+3rx585CVlVX6SElJ0WJK0gdOTk6IiorC/PnzIZPJ8O2338LNzQ0JCQmioxEZDblcjg8++AD79+9HnTp1EBMTg/bt2+O3334THY2qObXKSP369aFQKJ6acZ2RkYFGjRo9tf+VK1dw/fp19OvXDyYmJjAxMcF3332Hn376CSYmJrhy5Uq5X8fMzAzW1tZlHlT9mZiYYNGiRTh69Cjs7e1x8eJFeHp6Yt26dVxKnkiHevfujbi4OLi6uuL27dt45ZVXEBQUxN9D0hq1yoipqSnc3Nxw5MiR0m3FxcU4cuQIOnXq9NT+L7zwAv7++28kJCSUPvr3748ePXogISGBwy9Uru7duyMhIQH9+vVDQUEBZsyYgYEDB+LOnTuioxEZjaZNm+LEiRMYMWIElEolAgICMHLkSA6fklaoPUwTGBiI0NBQbNmyBRcuXMCUKVOQl5eHcePGAQBGjx6NefPmAXh0K+vWrVuXedSuXRtWVlZo3bo1TE1NNfvdULVRv3597NmzB59//jlMTU3x008/wdXVFVFRUaKjERmNmjVr4vvvv0dwcDAUCgXCwsLQuXPnCs9qE1WV2mVk6NChWLVqFRYsWIC2bdsiISEBBw4cKJ3UmpycjLS0NI0HJeMjk8nw7rvv4tSpU6WXHfbs2RMLFixAUVGR6HhERkEmk2HGjBk4evQoGjZsiL/++gvu7u7Yv3+/6GhUjcgkAxgEzM7Oho2NDbKysjh/xEjl5eVhxowZpevZeHl5ISwsDE5OToKTERmP1NRUDB48GKdOnYJMJsOiRYvwwQcfQC7XnzuL5OXloVatWgCA3NxcWFpaCk5k3FR9/9aff0FElbC0tMTGjRsRFhYGKysr/PHHH3B1dUVkZKToaERGw97eHr/99hsmTZoESZIwf/58+Pr6IisrS3Q0MnAsI2RQhg0bhoSEBHTo0AH37t3D4MGDMWXKFDx48EB0NCKjYGZmhvXr1+Obb76Bqakp9uzZgw4dOuC///2v6GhkwFhGyOA0bdoUx48fx9y5cwEA69evR4cOHXD+/HnByYiMx/jx43H8+HE4ODjg0qVL8PT0xM6dO0XHIgPFMkIGqUaNGli+fDkOHToEW1tbnDt3Dh4eHvj666+5FgKRjnh4eCA+Ph7du3dHbm4u/Pz88P7770OpVIqORgaGZYQMWq9evXD27Fn07t0bDx48wKRJk/Dmm2/i33//FR2NyCg0bNgQv/76K2bNmgUA+Oyzz/Daa69xXSBSC8sIGTxbW1vs27cPK1euhImJCSIiItC2bVucOHFCdDQio2BiYoJVq1Zh27ZtqFmzJg4fPgw3NzecPn1adDQyECwjVC3I5XLMnj0bJ06cgIuLC5KTk9G1a1csWbKEp4yJdOStt97Cn3/+CRcXFyQlJcHLywvfffed6FhkAFhGqFrx8PDA6dOnS5ew/uijj9CrVy+kpqaKjkZkFNq0aYPY2Fj07dsXDx8+xJgxYzB9+nQUFBSIjkZ6jGWEqh1ra2t8//332Lx5MywtLXHs2DG4urril19+ER2NyCjUqVMHP//8MxYuXAgACAkJwSuvvMLVualCLCNULclkMowZMwanT59Gu3btcOfOHfTr1w8zZ85Efn6+6HhE1Z5cLsfHH3+Mn376CdbW1jh+/Djc3Nw4l4vKxTJC1VqLFi1w8uRJzJw5EwAQHByMTp064dKlS2KDERmJfv36ITY2Fq1atUJaWhq6d++OL7/8kpfgUxksI1TtmZmZYe3atfjll19Qv359nDlzBu3bt8fmzZv5HyKRDrRo0QKnTp3CkCFDUFhYiGnTpuHtt9/myslUimWEjMbrr7+Os2fPokePHsjLy8O4ceMwcuRIZGdni45GVO3VqlUL27dvx4oVKyCXy7F582Z4e3sjKSlJdDTSAywjZFQaN26MX3/9FUuWLIFCoUBYWBjatWuH2NhY0dGIqj2ZTIY5c+bg4MGDqFevHuLj4+Hm5obDhw+LjkaCsYyQ0VEoFPjggw/w+++/o0mTJrh69So6d+6MlStXori4WHQ8omrPx8cH8fHxaN++Pe7cuYPevXtj5cqVHDY1YiwjZLQ6d+6MhIQE+Pn5oaioCO+99x769u2LjIwM0dGIqr0mTZrg+PHjGDt2LIqLi/Hee+9h6NChyM3NFR2NBGAZIaNWu3Zt/Pjjj/j6669hYWGBgwcPwtXVFYcOHRIdjajas7CwwKZNmxASEgITExPs2LEDnp6evNrNCLGMkNGTyWSYMGEC4uLi0Lp1a2RkZKB3796YO3cuV40k0jKZTIapU6fit99+g52dHf773//Cw8MDP//8s+hopEMsI0T/06pVK8TExGDq1KkAgBUrVsDb2xtXr14VnIyo+uvcuTPi4+Ph5eWF7Oxs9O/fHwsXLuQ8LiPBMkL0GAsLC4SEhCAyMhK1a9dGTEwM2rZti23btomORlTt2dnZ4ejRo5g+fToAYNGiRejXrx/+/fdfwclI21hGiMoxaNAgnD17Fl26dEFOTg6GDx+Ot99+m5PriLTM1NQU69atw+bNm2Fubo59+/bBw8MDf//9t+hopEUsI0QVcHJywrFjx7BgwQLI5XJ8++23cHd3R0JCguhoRNXemDFj8Mcff6BJkya4cuUKOnbsiPDwcNGxSEtYRogqYWJigk8++QRHjx6Fvb09Ll68CE9PT3z++edcE4FIy9q3b4/4+Hj4+Pjg/v37GDZsGGbPno2ioiLR0UjDWEaIVNCtWzecPXsW/fv3R0FBAfz9/TFgwABkZmaKjkZUrdWrVw8HDhzA3LlzAQCrV69Gr169cOvWLcHJSJNYRohUVK9ePezevRvr1q2DmZkZfv75Z7i6uiIqKkp0NKJqTaFQYPny5YiIiECtWrUQFRUFNzc33sahGmEZIVKDTCbD9OnTcerUKbzwwgu4efMmevbsifnz5/PUMZGWDR48GKdOnUKLFi1w48YNdOnSBRs3bhQdizSAZYSoClxdXREXF4fx48dDkiQsXrwY3bt35x1IibSsZD2gkiHTd955B5MnT0Z+fj4AICsrq3RfHx+fMh+T/pJJBjALLzs7GzY2NsjKyoK1tbXoOERlhIeHY9KkScjOzkbt2rXxzTffYPDgwaJjEVVrxcXFWLp0KRYsWABJktCxY0ekpaWV+weBi4sLEhMTBaQkVd+/WUaINODq1asYPnw4Tp06BQCYNGkS1q5dCwsLC8HJiKq3/fv3Y/jw4bh3716l+7GQiKHq+zeHaYg0oGnTpoiOji6d8b9hwwZ4eHjg/PnzgpMRVW99+vTB0aNHn7nflStXOGSjx1hGiDSkRo0aWL58OQ4dOoRGjRrh/PnzcHd3x4YNG7gmCZEWvfvuuyrt9/rrr2s5CVUVywiRhvXq1Qtnz57Fa6+9hocPH2Ly5MkYMmQI769BpCXJycka3Y90j2WESAsaNmyIvXv3YvXq1ahRowZ27tyJtm3b4o8//hAdjajacXJy0uh+pHssI0RaIpfLERgYiBMnTsDFxQXJycno1q0bFi9eDKVSKToeUbWxd+9eje5HuscyQqRl7u7uOH36NEaMGAGlUon58+fDx8cHqampoqMRVQs2NjZwcXGpdB8XFxfY2NjoKBGpi2WESAesra3xww8/YMuWLbC0tERUVBRcXV3x888/i45GVC0kJiZWWEjMzc1x+fJlHScidbCMEOnQ6NGjcfr0abRr1w537txB//794e/vX7p6JBFVXWJiIu7duwcvLy84OjrCzc0N5ubmePjwITZv3iw6HlWCZYRIx1q0aIGTJ08iICAAAPD555+jY8eOuHjxouBkRIbPxsYGx48fR3JyMuLi4rBo0SIAwKxZs5CRkSE4HVWEZYRIADMzM6xZswZ79+5F/fr1kZCQADc3N2zevJlrkhBpUEBAANq1a4d///0XM2fOFB2HKsAyQiRQ3759cfbsWfTs2RN5eXkYN24cRowYgezsbNHRiKoFExMThIaGQi6XIzw8HPv27RMdicrBMkIkWOPGjXHo0CEsXboUCoUC27ZtQ7t27RATEyM6GlG14ObmVjosOmXKFOTm5gpORE9iGSHSAwqFAvPmzUN0dDSaNGmCq1evwsvLCytWrEBxcbHoeEQG75NPPoGzszOSk5Mxf/580XHoCSwjRHqkU6dOSEhIwJAhQ1BUVIS5c+eiT58+SE9PFx2NyKBZWlpi/fr1AB5NGo+NjRWciB7HMkKkZ2rXro3t27fj66+/hoWFBQ4dOgRXV1ccPHhQdDQig9a7d2+MGDECxcXFeOedd1BYWCg6Ev0PywiRHpLJZJgwYQLi4uLQpk0b3Lp1C6+99hree+89FBQUiI5HZLDWrl2LevXq4a+//sLq1atFx6H/YRkh0mOtWrXCqVOnMHXqVADAypUr0aVLF1y5ckVwMiLD1KBBA6xZswbAo3kkiYmJghMRwDJCpPcsLCwQEhKCXbt2oU6dOoiNjUW7du0QFhYmOhqRQRo1ahR69eqFhw8fYtKkSVzbRw+wjBAZiIEDB+Ls2bPw9vZGTk4ORowYgXHjxvEyRSI1yWQyrF+/HhYWFjh69Ci2bNkiOpLRYxkhMiCOjo44evQoFi5cCLlcjs2bN8PNzQ1nzpwRHY3IoDRt2hSffPIJgEdLxd+6dUtwIuPGMkJkYExMTPDxxx/j6NGjsLe3x6VLl9CxY0cEBwfzdDORGgICAtC2bVvcvXuXS8ULxjJCZKC6deuGs2fPYsCAASgoKMDMmTPRr18/3L59W3Q0IoNgYmKCb775BnK5HNu2bcP+/ftFRzJaLCNEBqxevXrYtWsXvvjiC5iZmWHv3r1wdXXFsWPHREcjMghubm6lZ0UmT57MOViCsIwQGTiZTIZp06bh1KlTeOGFF5CWloZXXnkF8+fPR1FRkeh4RHpv0aJFXCpeMJYRomrC1dUVcXFxeOeddyBJEhYvXoxu3bohKSlJdDQivcal4sVjGSGqRiwtLREaGorw8HBYW1vjxIkTaNu2LXbu3Ck6GpFe41LxYrGMEFVDQ4cORUJCAjw9PXHv3j34+flh0qRJuH//PgAgKysLXbp0gZOTE7p06YKsrCzBiYnEW7t2LerWrYu//vqrdJVW0g2ZZADXAmZnZ8PGxgZZWVmwtrYWHYfIYBQWFmLBggX47LPPIEkSXnrpJWRnZyMlJeWpfV1cXLg0Nhm9LVu2YOzYsTA3N8fff/+NZs2aiY5k0FR9/+aZEaJqrEaNGli2bBkOHTqERo0a4fz58+UWEQC4cuUK/+Mlozd69Gj4+Pjg4cOHmDx5Mtfu0RGWESIj4OPjg+jo6Gfud+XKFQ7ZkFF7fKn4I0eOcKl4HeEwDVE1kpOTg6SkpDKP69evIykpCadPn1ZpUp6XlxeOHz+ug7RE+mvFihWYO3cu6tatiwsXLqBhw4aiIxkkVd+/WUaIDIQkScjMzHyqbDz++Pfff5/76zg6OiI5OVkDiYkMV1FRETw8PJCQkIBhw4bxLtlVpOr7t4kOMxFRJZRKJdLS0sqczXj8kZycXHo1TGXq1q2LJk2aPPVYtGgRzp49+8znOzk5aeLbITJoJUvFd+jQAdu2bcOoUaPQp08f0bGqLZ4ZIdKR/Px8pKSkVHhWIyUlRaUVU+3s7MotGyUPKyurcp+XlZWF2rVrP/P17927BxsbG3W/PaJqadasWVizZg2aNGmCc+fOoVatWqIjGRQO0xDpWHnzNR5/pKenP3NmvomJCRwdHSssGo6OjjAzM6tyxmbNmuHKlSsVfp6X9xKVlZeXh5deeglJSUkICAjg+iNqYhkh0iBJknDnzp0Ki8b169dVmq9hYWFRplw4OzuX+djOzg4KhUKr30tFhYRFhKh8Bw4cQJ8+fSCXy/Hnn3/Cw8NDdCSDodUyEhISgpUrVyI9PR2urq5Yt24dOnToUO6+oaGh+O6773Du3DkAj+6QuHTp0gr3Lw/LCGnb4/M1KnqoMl+jTp06lQ6h1K9fHzKZTAffUeVu3rwJe3t7AEDHjh1x4MABDs0QVWLEiBEICwuDq6srYmNjUaNGDdGRDILWysj27dsxevRorF+/Hp6enggKCsKOHTtw8eLFci99GjFiBLy8vNC5c2eYm5vjs88+w65du3D+/PnS/ww19c0QVaSgoOCp+RqPTxK9ceOGSpe9NmrU6KmzGarM19A3eXl5pWPfubm5sLS0FJyISL/dunULL774Iu7evYvly5dj7ty5oiMZBK2VEU9PT3h4eOCLL74AABQXF8PR0RHvvvsu3n///Wc+X6lUok6dOvjiiy8wevRolb4mywg9S25ubqVnNdLS0lSar+Hg4FDpfA1zc3MdfUfaxTJCpD4uFa8+rVzaW1BQgPj4eMybN690m1wuh4+PD06ePKnSa9y/fx+FhYWoW7duhfvk5+cjPz+/9OPs7Gx1YlI1I0kS7t69W+Elr0lJSbh79+4zX+fJ+RpPPho3bqz1+RpEZLhGjx6NH374AYcPH8bkyZPx66+/6sWwa3WgVhnJzMyEUqmEra1tme22trb4559/VHqNuXPnonHjxvDx8alwn2XLluGTTz5RJxoZsOLi4mfO18jLy3vm69SuXbvcklEyrKIv8zWIyDCVLBXfpk2b0qXix44dKzpWtaDTRc+WL1+O8PBwREVFVXq6e968eQgMDCz9ODs7G46OjrqISFpQ3nyNJ9fXUHW+RmVnNjiER0Ta5uLigo8//hhz587FrFmz0LdvXy4VrwFqlZH69etDoVAgIyOjzPaMjAw0atSo0ueuWrUKy5cvx+HDh/Hyyy9Xuq+ZmdlzraVAupWXl1fpWY2bN28+c76GQqEoM1/jyUmi1Wm+BhEZtoCAAISFheHs2bMICAjA1q1bRUcyeGqVEVNTU7i5ueHIkSMYOHAggEen2I8cOYLp06dX+LwVK1ZgyZIlOHjwINzd3Z8rMOnW4/M1KnrcuXPnma9jbm7+zPkaJia8OwER6b8aNWrgm2++gaenJ8LCwjBy5EguFf+c1P7fPzAwEGPGjIG7uzs6dOiAoKAg5OXlYdy4cQAeTfCxt7fHsmXLAACfffYZFixYgLCwMDg7OyM9PR0AUKtWLS6rqweKi4uRnp5e7uWumpivUfJo0KAB52sQUbXh7u4Of39/rF27FlOmTOFS8c9J7TIydOhQ3L59GwsWLEB6ejratm2LAwcOlE5qTU5OhlwuL93/q6++QkFBAfz8/Mq8zsKFC/Hxxx8/X/pqQKlUIjo6GmlpabCzs4O3t7dGr+goKCjAjRs3Kp2vUVBQ8MzXsbW1rbRscMEsIjI2ixYtQmRkJJKSkrBgwQIuFf8cuBy8QJGRkfD398eNGzdKtzk4OCA4OBi+vr4qvUZeXh6Sk5MrvOS1KvM1nnw4OTlxvkY1w3VGiDSDS8VXjvem0XORkZHw8/N7qiiUDGVERERg0KBB+Pfffyudr5GZmfnMr2Vubg4nJ6cKL3nlfA3jwzJCpDlcKr5iLCN6TKlUwtnZucwZkSeZmJjAzMxMpfkaNjY2lQ6hNGzYkPM1qAyWESLNeXyp+M8++wzvvfee6Eh6g2VEj0VFRaFHjx4q79+wYcNyL3flfA2qKpYRIs3avHkzxo0bB3Nzc5w7dw4uLi6iI+kFrSwHT5qRlpam0n6rVq3C1KlTYWFhoeVERET0PMaMGYMffvgBR44cwaRJk7hUvJrkz96FNM3Ozk6l/dzc3FhEiIgMQMlS8ebm5jhy5Ai+++470ZEMCsuIAN7e3nBwcKiwNctkMjg6OsLb21vHyYiIqKqaNWtWumRFYGAgbt26JTaQAWEZEUChUCA4OBgAniokJR8HBQXxDrJERAYmMDAQrq6uuHv3LgICAkTHMRgsI4L4+voiIiIC9vb2ZbY7ODggIiJC5XVGiIhIf5QsFS+XyxEWFob9+/eLjmQQeDWNYNpegZWoPLyahki7AgMDsXbtWjRp0gTnz5832t8xXtpLRBViGSHSrtzcXLRu3RpJSUkIDAzE6tWrRUcSQtX3bw7TEBERaVitWrXw1VdfAXg0BzAuLk5wIv3GMkJERKQFffr0wbBhw1BcXIwJEyagsLBQdCS9xTJCRESkJUFBQahbty4SEhKwdu1a0XH0FssIERGRljRs2LB0vsjChQtx5coVwYn0E8sIERGRFo0ZMwavvPIKHj58iEmTJj11t3ZiGSEiItIqLhX/bCwjREREWvbkUvG3b98WG0jPsIwQERHpAJeKrxjLCBERkQ7UqFEDoaGhkMvl2Lp1Kw4ePCg6kt5gGSEiItIRDw8PzJgxAwAwefJk5OXlCU6kH1hGiIiIdOjTTz9FkyZNcP36dSxYsEB0HL3AMkJERKRDXCr+aSwjREREOsal4stiGSEiIhLg8aXig4KCRMcRimWEiIhIgCeXir969argROKwjBAREQkyZswY9OzZEw8ePDDqpeJZRoiIiASRyWTYsGEDzM3NcfjwYXz//feiIwnBMkJERCRQs2bNsHDhQgDGu1Q8ywgREZFgs2bNgqurK+7cuWOUS8WzjBAREQlm7EvFs4wQERHpAWNeKp5lhIiISE98+umncHJywvXr10vnkRgDlhEiIiI98fhS8WvXrkV8fLzgRLrBMkJERKRH+vbti7feeqt0qfiioiLRkbSOZYSIiEjPBAUFoU6dOjhz5gzWrl0rOo7WsYwQERHpGVtbW6NaKp5lhIiISA+NHTu2dKn4iRMn4tixY9i2bRuioqKgVCpFx9MomWQAC+FnZ2fDxsYGWVlZsLa2Fh2HyODl5eWhVq1aAIDc3FxYWloKTkRE5UlMTESrVq1QWFhYZruDgwOCg4Ph6+srKJlqVH3/5pkRIiIiPfXXX389VUQAIDU1FX5+foiMjBSQSvNYRoiIiPSQUqmEv79/uZ8rGdSYOXNmtRiyYRkhIiLSQ9HR0bhx40aFn5ckCSkpKYiOjtZhKu1gGSEiItJDaWlpGt1Pn7GMEBER6SE7OzuN7qfPWEaIiIj0kLe3NxwcHCCTycr9vEwmg6OjI7y9vXWcTPNYRoiIiPSQQqFAcHAwADxVSEo+DgoKgkKh0Hk2TWMZITJCj8++//3336vFbHyi6sjX1xcRERGwt7cvs93BwQERERF6v86IqrjoGZGRiYyMxIwZM5Camlq6zVAWUCIyVkqlEtHR0UhLS4OdnR28vb0N4oyIqu/fLCNERiQyMhJ+fn548te+5JRvdfpLi4jE4wqsRFRGyQJK5f39Ud0WUCIiw8IyQmQkjGkBJSIyLCwjREbCmBZQIiLDwjJCZCSMaQElIjIsLCNERsKYFlAiIsPCMkJkJIxpASUiMiwsI0RGxFgWUCIiw8J1RoiMkKEuoEREhkXV928THWYiIj2hUCjQvXt30TGIiABwmIaIiIgEYxkhIiIioVhGiIiISCiWESIiIhKKZYSIiIiEYhkhIiIioVhGiIiISCiWESIiIhKKZYSIiIiEMogVWEtWrM/OzhachIiIiFRV8r79rDvPGEQZycnJAQA4OjoKTkJERETqysnJgY2NTYWfN4gb5RUXF+PmzZuwsrJ66tbnVL7s7Gw4OjoiJSWFNxc0EDxmhonHzfDwmOmOJEnIyclB48aNIZdXPDPEIM6MyOVyODg4iI5hkKytrfnLZmB4zAwTj5vh4THTjcrOiJTgBFYiIiISimWEiIiIhGIZqabMzMywcOFCmJmZiY5CKuIxM0w8boaHx0z/GMQEViIiIqq+eGaEiIiIhGIZISIiIqFYRoiIiEgolhEiIiISimXEgIWEhMDZ2Rnm5ubw9PRETExMhfuGhobC29sbderUQZ06deDj41Pp/qQd6hyzx4WHh0Mmk2HgwIHaDUhPUfeY3bt3D9OmTYOdnR3MzMzQokUL7Nu3T0dpqYS6xy0oKAgtW7aEhYUFHB0dERAQgIcPH+ooLUEigxQeHi6ZmppKmzZtks6fPy9NmDBBql27tpSRkVHu/sOHD5dCQkKkM2fOSBcuXJDGjh0r2djYSDdu3NBxcuOl7jErce3aNcne3l7y9vaWBgwYoJuwJEmS+scsPz9fcnd3l/r27SsdP35cunbtmhQVFSUlJCToOLlxU/e4bd26VTIzM5O2bt0qXbt2TTp48KBkZ2cnBQQE6Di58WIZMVAdOnSQpk2bVvqxUqmUGjduLC1btkyl5xcVFUlWVlbSli1btBWRnlCVY1ZUVCR17txZ+uabb6QxY8awjOiYusfsq6++kpo2bSoVFBToKiKVQ93jNm3aNKlnz55ltgUGBkpeXl5azUn/j8M0BqigoADx8fHw8fEp3SaXy+Hj44OTJ0+q9Br3799HYWEh6tatq62Y9JiqHrNFixahYcOGGD9+vC5i0mOqcsx++ukndOrUCdOmTYOtrS1at26NpUuXQqlU6iq20avKcevcuTPi4+NLh3KuXr2Kffv2oW/fvjrJTAZyozwqKzMzE0qlEra2tmW229ra4p9//lHpNebOnYvGjRuX+YUl7anKMTt+/Dg2btyIhIQEHSSkJ1XlmF29ehVHjx7FiBEjsG/fPiQmJmLq1KkoLCzEwoULdRHb6FXluA0fPhyZmZno0qULJElCUVERJk+ejA8++EAXkQmcwGqUli9fjvDwcOzatQvm5uai41A5cnJyMGrUKISGhqJ+/fqi45CKiouL0bBhQ3z99ddwc3PD0KFD8eGHH2L9+vWio1EloqKisHTpUnz55Zc4ffo0IiMjsXfvXnz66aeioxkNnhkxQPXr14dCoUBGRkaZ7RkZGWjUqFGlz121ahWWL1+Ow4cP4+WXX9ZmTHqMusfsypUruH79Ovr161e6rbi4GABgYmKCixcvwsXFRbuhjVxVfs/s7OxQo0YNKBSK0m0vvvgi0tPTUVBQAFNTU61mpqodt/nz52PUqFF45513AABt2rRBXl4eJk6ciA8//BByOf9u1zb+hA2Qqakp3NzccOTIkdJtxcXFOHLkCDp16lTh81asWIFPP/0UBw4cgLu7uy6i0v+oe8xeeOEF/P3330hISCh99O/fHz169EBCQgIcHR11Gd8oVeX3zMvLC4mJiaXFEQAuXboEOzs7FhEdqcpxu3///lOFo6RQSrx9m26InkFLVRMeHi6ZmZlJmzdvlv773/9KEydOlGrXri2lp6dLkiRJo0aNkt5///3S/ZcvXy6ZmppKERERUlpaWukjJydH1LdgdNQ9Zk/i1TS6p+4xS05OlqysrKTp06dLFy9elH755RepYcOG0uLFi0V9C0ZJ3eO2cOFCycrKStq2bZt09epV6dChQ5KLi4v05ptvivoWjA6HaQzU0KFDcfv2bSxYsADp6elo27YtDhw4UDppKzk5uUzT/+qrr1BQUAA/P78yr7Nw4UJ8/PHHuoxutNQ9ZiSeusfM0dERBw8eREBAAF5++WXY29vD398fc+fOFfUtGCV1j9tHH30EmUyGjz76CKmpqWjQoAH69euHJUuWiPoWjI5MkngOioiIiMThn2FEREQkFMsIERERCcUyQkREREKxjBAREZFQLCNEREQkFMsIERERCcUyQkREREKxjBAREZFQLCNEREQkFMsIERERCcUyQkREREKxjBAREZFQ/wdL81MwHOSMhwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Dirty code to print the ILP solution\n",
    "for y in rmpModel.getVars():\n",
    "    if y.X > EPS: \n",
    "        cluster = np.zeros(n)\n",
    "        for i in range(n):\n",
    "            if (rmpModel.getCoeff(coverConstrs[i], y)):\n",
    "                cluster[i] = 1\n",
    "        median = computeClusterCost(cluster)[0]\n",
    "        for i in range(n):\n",
    "            if (rmpModel.getCoeff(coverConstrs[i], y) == 1):\n",
    "                plt.plot([x1[median], x1[i]], [y1[median], y1[i]], color='black', marker='o')\n",
    "            \n",
    "plt.show()"
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}