{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Small disks and geometric facility location: *a mixed integer conic quadratic tutorial*\n",
    "In this tutorial we are going to:\n",
    "* show examples of compact models constructed with *Fusion*,\n",
    "* write a conic model with binary variables,\n",
    "* apply the *Big-M* formulation,\n",
    "* construct a base model and use its various extensions to solve a number of different problems.\n",
    "\n",
    "### Introduction\n",
    "\n",
    "We are going to work with variations of the following geometric setup. Suppose we have $n$ of points $p_1,\\ldots,p_n\\in\\mathbb{R}^d$ and we want to find a configuration of $k$ balls with centers $x_1,\\ldots,x_k\\in\\mathbb{R}^d$ and radii $r_1,\\ldots,r_k$, respectively, which covers some/all of the points $p_i$ and satisfies additional constraints depending on the problem. For visualization purposes we will focus on the plane and disks, that is the case $d=2$, but the code and mathematical formulation remain valid in all dimensions.\n",
    "\n",
    "For example, the next snippet displays a set of points in $\\mathbb{R}^2$. Here is a sample question: how many of those points can we cover with three disks of radius at most $0.1$ each?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAEjdJREFUeJzt3X2MXFd5x/HvE4dQlwZMsWnD2oldatK6JMJ0myKhAoXQOKDGJoTiVFSJFGoFkVIJGmEUFCGjKmks8VLVf+CmUVMkMCFCwaUGq+RFFFRTb+SQ4CAT44Zm7YgsL4ZKGBLTp3/MOEw2s97dmTszd+75fiQr9945mvPkev3bM+feeyYyE0lSWc4YdQGSpOEz/CWpQIa/JBXI8JekAhn+klQgw1+SCmT4S1KBDH9JKpDhL0kFOnPUBcxl+fLluXr16lGXIUlj5f777/9+Zq6Yr11tw3/16tVMTU2NugxJGisR8d2FtKtk2iciNkTEoYg4HBFbT9PuiojIiJisol9JUm/6Dv+IWALsAC4F1gFXRsS6Lu3OBt4DfL3fPiVJ/ali5H8RcDgzj2Tmk8AuYGOXdh8GbgF+VkGfkqQ+VBH+E8BjHfvT7WNPi4j1wKrM/EIF/UmS+lRF+EeXY09/SUBEnAF8FHjfvG8UsSUipiJiamZmpoLSJEndVBH+08Cqjv2VwLGO/bOBlwP3RcSjwKuA3d0u+mbmzsyczMzJFSvmvVNJktSjKsJ/P7A2ItZExFnAZmD3qRcz88eZuTwzV2fmamAfcFlmeh+nJI1I3+GfmSeB64C9wLeAOzLzYERsi4jL+n1/SVL1KnnIKzP3AHtmHbtxjravq6JPSVLvXNtHkgpk+EtSgQx/SSqQ4S9JBTL8JalAhr8kFcjwl6QCGf6SVCDDX5IKZPhLUoEMf0kqkOEvSQUy/CWpQIa/JBXI8JekAhn+klQgw1+SCmT4S1KBDH9JKpDhL0kFMvwlqUCGvyQVyPCXpAIZ/pJUIMNfkgpk+EtSgQx/SSqQ4S9JBTL8JalAhr8kFaiS8I+IDRFxKCIOR8TWLq9fGxEPRcQDEfHViFhXRb+SpN70Hf4RsQTYAVwKrAOu7BLun8rMCzLzFcAtwEf67VeS1LsqRv4XAYcz80hmPgnsAjZ2NsjMn3TsPg/ICvqVJPXozAreYwJ4rGN/GvjD2Y0i4t3Ae4GzgNdX0K8kqUdVjPyjy7Fnjewzc0dmvhR4P/DBrm8UsSUipiJiamZmpoLSJEndVBH+08Cqjv2VwLHTtN8FbOr2QmbuzMzJzJxcsWJFBaVJkrqpIvz3A2sjYk1EnAVsBnZ3NoiItR27bwYeqaBfSVKP+p7zz8yTEXEdsBdYAtyWmQcjYhswlZm7gesi4mLgKeBHwFX99itJ6l0VF3zJzD3AnlnHbuzY/usq+pEkVcMnfCWpQIa/JBXI8JekAhn+klQgw1+SCmT4S1KBDH9JKpDhL0kFMvwlqUCGvyQVyPCXpAJVsrZPSe46cJTtew9x7PgJXrJsKddfcj6b1k+MuixJWpQiwr+qwL7rwFE+8LmHOPHULwA4evwEH/jcQwBj/wvAX2pSWRo/7XMqsI8eP0Hyy8C+68DRRb/X9r2Hng7+U0489Qu27z1UUbWjUeU5kjQeGh/+VQb2seMnFnV8XDT1l5qkuTU+/KsM7JcsW7qo4+Oiqb/UJM2t8eFfZWBff8n5LH3OkmccW/qcJVx/yfk91VYXTf2lJmlujQ//KgN70/oJbrr8AiaWLSWAiWVLuenyC8b+wmhTf6lJmlvj7/Y5FcxV3cmyaf3E2If9bFWfI0n1F5k56hq6mpyczKmpqVGXsSjeLilp1CLi/sycnK9d40f+w9LkZwAkNU/j5/yHxdslJY0Tw78i3i4paZwY/hXxdklJ48Twr4i3S0oaJ17wrYi3S0oaJ4Z/hZr4DICkZjL8pXn4/IaayPBvKAOrGj6/oaaq5IJvRGyIiEMRcTgitnZ5/b0R8XBEPBgRd0fEeVX0q+5cn786Pr+hpuo7/CNiCbADuBRYB1wZEetmNTsATGbmhcCdwC399qu5GVjV8fkNNVUVI/+LgMOZeSQznwR2ARs7G2TmvZn50/buPmBlBf1qDgZWdXx+Q01VRfhPAI917E+3j83lGuCLFfSrORhY1fH5DTVVFeEfXY51XSo0It4BTALb53h9S0RMRcTUzMxMBaWVycCqTlO/w0Gq4m6faWBVx/5K4NjsRhFxMXAD8NrM/Hm3N8rMncBOaC3pXEFtRfKBs2r5/IaaqIrw3w+sjYg1wFFgM/DnnQ0iYj3wCWBDZj5RQZ+ah4El6XT6nvbJzJPAdcBe4FvAHZl5MCK2RcRl7WbbgV8DPhsRD0TE7n77lST1rpKHvDJzD7Bn1rEbO7YvrqIfSVI1XNVTkgrk8g6SBsIlRurN8JdUOddEqj/DXxoD4zaKPt0SI3WuuySGv1Rz4ziKdomR+vOCr1Rz47hQn0uM1J/hL9XcOI6iXWKk/gx/qebGcRTtmkj155y/VHPXX3L+M+b8YTxG0S4xUm+Gv1RzLtSnQTD8pTHgKFpVc85fkgpk+EtSgQx/SSqQc/6SRmbclq1oEsNftWIYlGMcl61oEqd9VBunwuDo8RMkvwyDuw4cHXVpGoBxXLaiSRz5q2dVj9JdCbIs47hsRZM48ldPBjFKNwzKMo7LVjRJMeF/14GjvPrme1iz9d949c33OJXQp0F8ZDcMyuLib6NVRPg7l1y9QYzSDYOyuPjbaDV6zv/UnPTRLoHkXHJ/XrJsadfz2s8ofdzWsPHOpP65bMXoNDb8Z99G1o1zyb0b1EqT4xIG3qaocdfY8O82Jz2bc8m9G7dRetW8M6lcTfnE19jwn29U71xy/8ZllD4I3plUpiZ94mvsBd/Tjeq9sHR63hk1P+9MKlOTHkxrbPjPdefIx97+Cr629fUG/xy8M2phvDOpTE36xNfY8Pc2st40aWQzSP58lalJn/gaO+cPZc9J96pJI5tB8+erPOP6fcrdNHbkr940aWQjVa1Jn/gqGflHxAbg48AS4NbMvHnW668BPgZcCGzOzDur6FfVa9LIRhqEpnzi63vkHxFLgB3ApcA64MqIWDer2f8AVwOf6rc/DVaTRjaS5lbFyP8i4HBmHgGIiF3ARuDhUw0y89H2a/9XQX8asKaMbCTNrYo5/wngsY796fYxSVJNVRH+0eVY9vRGEVsiYioipmZmZvosS5I0lyrCfxpY1bG/EjjWyxtl5s7MnMzMyRUrVlRQmiSpmyrCfz+wNiLWRMRZwGZgdwXvK0kakL7DPzNPAtcBe4FvAXdk5sGI2BYRlwFExB9ExDTwNuATEXGw334lSb2r5D7/zNwD7Jl17MaO7f20poMkSTXgE76SVKBGr+0zSE35QgdJZTL8e9CkL3SQVCanfXrgsseSxp3h3wOXPZY07gz/HrjssaRxZ/j3wK/wkzTuirvgW8VdOqfae7ePpHFVVPhXeZeOyx5LGmdFTft4l44ktRQV/t6lI0ktRYW/d+lIUktR4e9dOpLUUtQFX+/SkaSWosIfvEtHkqCwaR9JUovhL0kFMvwlqUCGvyQVyPCXpAIZ/pJUIMNfkgpU3H3+klQnVSwz3wvDX5JGpMpl5hfL8FffRjVykcbd6ZaZN/xVa6McuUjjbpTLzHvBV33xC3Kk3o1ymXnDX33xC3Kk3o1ymXnDX33xC3Kk3m1aP8FNl1/AxLKlBDCxbCk3XX6Bd/uo/q6/5PxnzPmDX5AjLcaolpmvZOQfERsi4lBEHI6IrV1ef25EfKb9+tcjYnUV/Wr0RjlykdS7vkf+EbEE2AG8EZgG9kfE7sx8uKPZNcCPMvO3I2Iz8HfA2/vtW/XgF+RI46eKkf9FwOHMPJKZTwK7gI2z2mwEbm9v3wm8ISKigr4lST2oIvwngMc69qfbx7q2ycyTwI+BF81+o4jYEhFTETE1MzNTQWmSpG6qCP9uI/jsoQ2ZuTMzJzNzcsWKFRWUJknqporwnwZWdeyvBI7N1SYizgReAPywgr4lST2oIvz3A2sjYk1EnAVsBnbParMbuKq9fQVwT2Y+a+QvSRqOvu/2ycyTEXEdsBdYAtyWmQcjYhswlZm7gX8CPhkRh2mN+Df3268kqXeVPOSVmXuAPbOO3dix/TPgbVX0JUnqn8s7SFKBDH9JKpDhL0kFMvwlqUCGvyQVyPCXpAIZ/pJUIMNfkgpk+EtSgQx/SSqQ4S9JBTL8JalAhr8kFcjwl6QCGf6SVCDDX5IKZPhLUoEMf0kqkOEvSQUy/CWpQIa/JBXI8JekAhn+klQgw1+SCmT4S1KBDH9JKpDhL0kFMvwlqUCGvyQVqK/wj4hfj4h/j4hH2v994RztvhQRxyPiC/30J0mqRr8j/63A3Zm5Fri7vd/NduAv+uxLklSRfsN/I3B7e/t2YFO3Rpl5N/C/ffYlSapIv+H/G5n5OED7vy/uvyRJ0qCdOV+DiPgy8JtdXrqh6mIiYguwBeDcc8+t+u0lSW3zhn9mXjzXaxHxvYg4JzMfj4hzgCf6KSYzdwI7ASYnJ7Of95Ikza3faZ/dwFXt7auAz/f5fpKkIeg3/G8G3hgRjwBvbO8TEZMRceupRhHxH8BngTdExHREXNJnv5KkPsw77XM6mfkD4A1djk8B7+zY/6N++pEkVcsnfCWpQIa/JBXI8JekAhn+klQgw1+SCmT4S1KBDH9JKpDhL0kFMvwlqUCGvyQVyPCXpAIZ/pJUIMNfkgpk+EtSgSKznl+YFREzwHcreKvlwPcreJ+q1bUusLZe1LUuqG9tda0Lxru28zJzxXxvUtvwr0pETGXm5KjrmK2udYG19aKudUF9a6trXVBGbU77SFKBDH9JKlAJ4b9z1AXMoa51gbX1oq51QX1rq2tdUEBtjZ/zlyQ9Wwkjf0nSLI0I/4jYEBGHIuJwRGzt8vp7I+LhiHgwIu6OiPNqVNu1EfFQRDwQEV+NiHV1qa2j3RURkRExlLsfFnDOro6ImfY5eyAi3jmMuhZSW7vNn7V/3g5GxKfqUFdEfLTjfH07Io4Po64F1nZuRNwbEQfa/0bfVKPazmtnxoMRcV9ErBxSXbdFxBMR8c05Xo+I+Pt23Q9GxCsX3UlmjvUfYAnwHeC3gLOAbwDrZrX5Y+BX29vvAj5To9qe37F9GfClutTWbnc28BVgHzBZh7qAq4F/qOnP2lrgAPDC9v6L61DXrPZ/BdxWo3O2E3hXe3sd8GiNavsscFV7+/XAJ4dU22uAVwLfnOP1NwFfBAJ4FfD1xfbRhJH/RcDhzDySmU8Cu4CNnQ0y897M/Gl7dx8wlN/eC6ztJx27zwOGdRFm3traPgzcAvysZnWNwkJq+0tgR2b+CCAzn6hJXZ2uBD49hLpgYbUl8Pz29guAYzWqbR1wd3v73i6vD0RmfgX44WmabAT+JVv2Acsi4pzF9NGE8J8AHuvYn24fm8s1tH5jDsOCaouId0fEd2iF7HvqUltErAdWZeYXhlTTgupqe2v74+6dEbFqOKUtqLaXAS+LiK9FxL6I2FCTuoDWNAawBrhnCHXBwmr7EPCOiJgG9tD6ZDIMC6ntG8Bb29tvAc6OiBcNobb5LDb3nqUJ4R9djnUdPUfEO4BJYPtAK+rossuxZ9WWmTsy86XA+4EPDryqltPWFhFnAB8F3jekep7uusux2efsX4HVmXkh8GXg9oFX1bKQ2s6kNfXzOloj7FsjYlkN6jplM3BnZv5igPV0WkhtVwL/nJkraU1nfLL98zdoC6ntb4DXRsQB4LXAUeDkoAtbgMX8nXfVhPCfBjpHfivp8rExIi4GbgAuy8yf16m2DruATQOt6Jfmq+1s4OXAfRHxKK15xd1DuOg77znLzB90/B3+I/D7A65pwbW123w+M5/KzP8GDtH6ZTDquk7ZzPCmfGBhtV0D3AGQmf8J/Aqt9WtGXltmHsvMyzNzPa38IDN/PITa5rPYbHm2YVy8GPCFkTOBI7Q+yp66aPN7s9qsp3VhZ20Na1vbsf2nwFRdapvV/j6Gc8F3IefsnI7ttwD76nLOgA3A7e3t5bQ+mr9o1HW1250PPEr7+Z4anbMvAle3t3+3HWIDr3GBtS0Hzmhv/y2wbYjnbjVzX/B9M8+84Ptfi37/Yf2PDPgkvQn4djvgb2gf20ZrlA+tqYHvAQ+0/+yuUW0fBw6267r3dAE87NpmtR1K+C/wnN3UPmffaJ+z36nLOWv/Y/wI8DDwELC5DnW19z8E3Dysc7WIc7YO+Fr77/MB4E9qVNsVwCPtNrcCzx1SXZ8GHgeeojXKvwa4Fri24+dsR7vuh3r5t+kTvpJUoCbM+UuSFsnwl6QCGf6SVCDDX5IKZPhLUoEMf0kqkOEvSQUy/CWpQP8PwfgQq4m9DvAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe4f474da90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from mosek.fusion import *\n",
    "import numpy as np\n",
    "\n",
    "def data(n):\n",
    "    np.random.seed(1236)\n",
    "    return [[x,y] for x,y in np.broadcast(np.random.uniform(0.0, 1.0, n), np.random.uniform(0.0, 0.3, n))]\n",
    "def gaussianData(n):\n",
    "    np.random.seed(1236)\n",
    "    return [[x,y] for x,y in np.broadcast(np.random.normal(0.0, 1.0, n), np.random.normal(0.0, 1.0, n))]\n",
    "\n",
    "def display(P, R, X):\n",
    "\timport matplotlib.pyplot as plt\n",
    "\tplt.scatter(*zip(*P))\n",
    "\tplt.axis(\"equal\")\n",
    "\tfor i in range(len(R)):\n",
    "\t\tplt.gcf().gca().add_artist( plt.Circle((X[2*i],X[2*i+1]), R[i], fc=\"r\", color=\"r\", alpha=0.5) )\n",
    "\t\tplt.gcf().gca().plot(X[2*i],X[2*i+1], \"or\")\n",
    "\tplt.show()\n",
    "\n",
    "display(data(20), [], [])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A base model formulation in *Fusion*\n",
    "\n",
    "Our input is a collection of points $p_1\\ldots,p_n\\in\\mathbb{R}^d$ and an integer $k\\geq 1$. What we want to compute is the centers $x_1,\\ldots,x_k\\in\\mathbb{R}^d$ and the radii $r_1,\\ldots,r_k\\geq 0$ of $k$ balls $B_d(x_j;r_j)$, subject to various constraints on how the points and balls interact. The basic notion we need to model is that of which points are covered by which balls. To this end we introduce standard binary variables\n",
    "\n",
    "\\begin{equation}\n",
    "s_{ij} = \\left\\lbrace\n",
    "   \\begin{array}{cl}\n",
    "   1 & \\text{if the ball } B_d(x_j;r_j) \\text{ covers the point } p_i, \\\\\n",
    "   0 & \\text{otherwise}.\n",
    "   \\end{array}\\right.\n",
    "\\end{equation}\n",
    "\n",
    "For this to work we must ensure the implication\n",
    "\n",
    "\\begin{equation}\n",
    "s_{ij}=1 \\implies r_j\\geq \\lVert p_i-x_j \\rVert_2.\n",
    "\\end{equation}\n",
    "\n",
    "This can be achieved by introducing a constraint\n",
    "\n",
    "\\begin{equation}\n",
    "r_j\\geq \\lVert p_i-x_j \\rVert_2 - M(1-s_{ij})\n",
    "\\end{equation}\n",
    "\n",
    "where $M$, also known as *Big-M*, is large enough to guarantee that when $s_{ij}=0$ the corresponding constraint is trivially satisfied, for example the right-hand side is surely negative. The spread of the data points $p_i$ determines a safe value of $M$ for us:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bigM(P):\n",
    "    return 2*np.shape(P)[1]*max(np.amax(P,0)-np.amin(P,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we get a family of conic constraints\n",
    "\n",
    "\\begin{equation}\n",
    "(r_j + M(1-s_{ij}),\\  p_i-x_j) \\in \\mathcal{Q}^{d+1}, \\quad i=1,\\ldots,n,\\ j=1,\\ldots,k\n",
    "\\end{equation}\n",
    "\n",
    "and its direct implementation in *Fusion*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def basicModel(k, P, bigM):\n",
    "    n, dim = len(P), len(P[0])                              #in P each row is one p_i\n",
    "    M = Model(\"balls\")\n",
    "    R = M.variable(k, Domain.greaterThan(0.0))              #r_1...r_k\n",
    "    X = M.variable([k,dim], Domain.unbounded())             #each row is one x_j\n",
    "    S = M.variable([n,k], Domain.binary())\n",
    "\n",
    "    RRep = Var.repeat(R, n)\n",
    "    Penalty = Expr.flatten(Expr.mul(bigM, Expr.sub(1.0, S)))    #M*(1-S)\n",
    "    CoordDiff = Expr.sub(np.repeat(P,k,0), Var.repeat(X, n))    #each p_i meet each x_j once in Expr.sub\n",
    "\n",
    "    M.constraint(Expr.hstack(Expr.add(RRep,Penalty), CoordDiff), Domain.inQCone())\n",
    "\n",
    "    return M, R, X, S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now take this as a base and solve a number of different optimization problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 1. Minimal enclosing ball\n",
    "\n",
    "**Objective: find the smallest ball containing all of the points.** \n",
    "\n",
    "This is a classical problem, see the [bounding sphere](https://en.wikipedia.org/wiki/Bounding_sphere) and [smallest-circle problem](https://en.wikipedia.org/wiki/Smallest-circle_problem) Wikipedia articles. There is an algorithm with running time $O(n)$ in any fixed dimension, and also MOSEK can solve it efficiently for a large number of points.\n",
    "\n",
    "Since we want to include all points, the Big-M constant and the $s_{ij}$ indicators are superfluous. All we want to do is to set $k=1$ and minimize the radius of the single ball in the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnXl0XGd5/7/vnRlrJHtk2ZYsyZYdCxNZtuNsKCHgEIwxsSFh+ZUUCC2E4J4QoKf4HLYEyo/ttNBS0lAohYChtD9OKEsISxI7MUFQJ2RxYrLIlp04yiJLliVvkiONNDP3/f3x6kpXd+6+35nnc46OpNHMve/MaJ77vM/yfRjnHARBEETlIEW9AIIgCMJfyLATBEFUGGTYCYIgKgwy7ARBEBUGGXaCIIgKgww7QRBEhUGGnSAIosIgw04QBFFhkGEnCIKoMNJRnLSxsZGvWrUqilMTBEEklscee2yEc95kdb9IDPuqVauwb9++KE5NEASRWBhjL9i5H4ViCIIgKoxIPHaCiARZBl5+GRgbE19nzwKjo8DwMDAyIn4uFoFSSXwViwDnQColviQJSKeBefOAxYuBxkZgyRIglxNfCxaI79kswFjUz5aoYsiwE5XH1BRw7BgwMAA8+yzw3HPAyZPCcAPCQHM++5VOAzU1QCYjDLL2SzHwypcsA319wOSk+FlSbXyVC0FDA7B0KdDRAZxzDrBsmbgQSLRJJoKHDDuRbNRG/JlngEOHxM+AMLLz5gHz5wO1tcDCheF40qWSWFd/P9DbK9YBiAtHezuwdi0ZeyJQyLATyUKWgRdfBJ5+GnjkEWE8gVkjnssBbW3RGstUSlxIamuBRYtmby8WxUXomWfmGvuODuCyy4DOThHaIQiPkGEn4s/kpAipPP448PDDIk7OmAh3RG3EnZBOizU3NMzeViwCzz8vLlScA8uXAxs3AuvXAytWJOe5EbGCDDsRT06fFmGVhx8GnnpKhDfSaeHRNjZGvTpTegdH8cCRExjLF5DLZrBx9RJ0ttbr31l5TkuWCMN+9izw858DP/uZSMZedhlw0UXA6tUiD0AQNiDDTsSHQgHo6QF27RJGHQDq6oDmZmEAE0Dv4Cj2HDyOoiwDAMbyBew5eBwAjI27AmOzFTYAkM8D3d3Anj0iZLNxI7BpE7ByJVXdEKYk49NCVDYnTwIPPgjs3i3KEHO5xBqvB46cmDHqCkVZxgNHTlgbdi3ZrAjNACJks3cv8PvfA6tWAW95C3DBBeTFE7qQYSeiQZZFEnHPHuCxx8RtTU2JTx6O5QuObrdNOi2MPOfAqVPAt78tkrNvehNw+eWitJIgpvFs2BljWQB/BFAzfbyfc84/7/W4RPK5c/9RfG33IQycnsCyhlp8cusavKOjAXj0UeDuu4Hjx4VXmqQEqAW5bEbXiOeyGX9OwJiotFm0SCSV77oL+PWvgQ0bgK1bRSllKuXPuYjEwrhSduX2AIwxAPM552cZYxkAewF8jHP+kNFjurq6OGnFVDZ37j+Km+94ChOFEgCgpjCJKwZ68JmpQ2ivz4jOTSWWXEFoY+wAkJYkbFm71Hkoxi6yDJw4IRKvra3Ae98LnHdexVwsiVkYY49xzrus7ufZY+fiynB2+tfM9Je3qwWReL62+xAmCiWkS0VcfPQgrjq0F7VTefzvkia0b3hlKGtwVJ3iE8rxQz2vJIkwVlOTqCb6+teBV74SePe7gXPPTWSugvCGLzF2xlgKwGMAXgng3znnD/txXCK5DJ56GecdO4K3HfwjFk2cwYm6BpyqrQeK4ZzfU3WKRzpb6wM/hyENDaLDdmgI+Id/AC68ELjmGlETT1QNvuzVOOclzvmFANoAXMoYO097H8bYDYyxfYyxfcPDw36clogjnAMHDuBzj/4Prn/sN2BcxtGFzchnRPWGb7FmC8yqUyoexkSt/6pVomz0c58Dvvc9kdMgqgJfq2I456cZY90AtgF4WvO32wDcBogYu5/nJWLCCy8AP/kJcOAArmieh3sKLSiqcjhpScLG1eFUvQRWnZIkGANaWkQM/pFHgD/9SVTRXH11ReY3iFn8qIppAlCYNuq1ALYA+CfPKyOSw9SUqM741a9EQ9GqVVjNGLbUhx/jVgi8OiVJSJIolSwWgfvuEz0D27eLOniKv1ckfnjsrQB+NB1nlwD8lHP+Wx+OSySB558X2/z+flG2mJk1nFHGmjeuXlJWnQIAhaKM3sHR6GLgUZJOi8avsTHglluA170OeM97yHuvQPyoinkSwEU+rIVIEmovPZcTcrQxQjHc3YeHkZ8uuQSAfLEUWhI1tuRyQsr4oYeAJ58k770CoUJXwjnPPw988YvAnXeKLX5Mu0U7W+uRSZX/i1dNEtUMSRKVMvPmCe99507hyRMVAUkKEPZRe+kLFsTOS9eDkqgWqL33J54g771CII+dsMfQEPDlL8966TGXzlUwSpZWZRLVCEkS+RHFe//Rj8RFnEgs5LET1vT2At/4hvg5AV66Gr0kaphll1ZE0R1riOK9d3cDR48CH/3o3KEgRGIgj50whnPg/vuBr3xFKAk2N0e9Isd0ttZjy9qlMx56LpsJVrfFAUp3rBIWUrpjewdHo1uUJIl5rC+8AHzhCyKfQiQO8tgJfaamgNtvF7K6bW2J1v2OtMXfBF+12/2EMTFo++RJEX770IeASy+Nbj2EY8hjJ8o5cwb4l38R3np7e6KNepyJfWJ38WJR8fTNbwJ33CHGExKJgDx2Yi4vvgjceqsofVu1iqojPGAVP09Ed2xdnQjN3HmnaEL7m78RtxGxhgw7Mcvjj4vJPHV1syPZqhAzg2w32WlHXTLuid0Z0mmxc3viCaEYuWOHkAgmYgsZdkLw4IPAd78rEqRV5JFpDXV7Yx0ODIzpGmQAtqWA7cTPI9FudwtjoqFJkQO+6SYhMEbEEjLsBPCHPwDf/77w0rPZqFcTGnpe9ZP9Z8rup+5UtZvsNIufq7Vq4prYNaS5GRgenjXuVbyzizOUPK1mOAd+9zth1NvaqsqoA/petRFj+YKjZKdZnDzykkavKGGYf/xHkZMhYgd57NXMffcB//3fs5ohARBGA47bczipPlEMtV0jbqQuCQgvv/vwcKieuu/vw5IlohzyK18Bbr5ZqEYSsYE89mrl978PxagH3YDj5Rx2q0+UhObG1UuQ1gyINkp2Ko1RRuQLpdC89sDeh8WLhUzzP/2T6FQlYgMZ9mpk717gBz+Y1QcJiDDG09k9R+/gKHbu7cOtew5j594+9A6OGhrq89sW6naqOu1i7WytN714hKUwGej7sHix6Fb96leBwUHvxyN8gUIx1cZDDwG33SaSXg4aj9xs5cNowLFzDqPSwy1rl2LL2qWOnpfTZOfG1Uuwq+eYo7X7TeDvw5IlwMiI8Nw/8xlgqfFOhQgHMuzVxOHDwHe+A7S2OkqU2qnJ1iPoBhyzUIL6HGYe6/bL2wONdXe21qP70DDyxfKuzbAakUJphGpsFKWQt9wC/P3fC1lnIjIoFFMtDA8LhcZFi4SglwPcbuWdxKTdYHZ+9Tmibt3ftKYp0NfBiqDfhxmUUsjvfY/kByKGDHs1MDEB/Nu/iQ/bwoWOH+7WMAatrGh2fm3rvh5hecxRK0yGev62NmD/fuAXvxDltEQkUCim0pFl4Ic/FFULLkvSvGzlg2zAsbuuOLTuh92IpJcT2X55CFr6jIn/s9/8RlRcveY1wZ+TKIMMe6Xz298Cf/oT8IpXuD5EHAyjHnrrkhhDoSjj1j2Hy5KhdjVeEtHib4LbnIhvpNNC9vd73xPhGQ//e4Q7PBt2xtgKAP8FoAWADOA2zvk3vB6X8IHHHgN+/nPhQXlQaYxS08TM0GrXlc2kMFWUZxKVWoNmtd7IDaJPxELnvbZWTGS69VYxsGPx4nDOSwDwx2MvAvg45/xxxlgOwGOMsfs45wd8ODbhlhdfBP7jP4RQU8Z7LDkKTRM7hla9rp17+5AvzE3aOTFofhnEqL3+qJPFMyxaJGrbv/Ut4NOfJl3/EPFs2DnngwAGp38eY4wdBLAcABn2qBgdFZ5SbW2kSo1eDZxVNY722F4NmtPH6z0/wL4CpBFeX7dY6by3tAB9faLLeft20vcPCV9j7IyxVQAuAvCwzt9uAHADAKwkXYng4Bz48Y/FFKQVKyJbhp63vavnGLoPDWPTmqYyQ6VnzMwMrbrpRzGe2UyqzGMH7Bs0BgaO8koOhnJjZLSbSEtM92K0u2cIwFzjHtSFIVY5EUXu949/BC64ALjkkvDXUIX4ZtgZYwsA/ALADs55WecI5/w2ALcBQFdXF9VBBcX+/SJZumpVpMswUk7MF0tlhsrISDJAx8zqU5RlpKUU0pLk2qDpGXXl9p17++YYYKPdRNFALJKDz3nees/53gNDkHVKBJ2Gg2Kn8y5JIon6wx8CHR2uSm4JZ/hi2BljGQij/mPO+R1+HJNwweio0IBZulR8mALEasqQWfhDa6iMjKRT8sUSsukU0inhuTs1aEYhDGA2HKNcdNysT/289Z6znlHXnt8usdN5nz8fOHVK7CY//GEKyQSMH1UxDMBOAAc557d4XxLhCiUEk8+L9u4AMUtqApjzsxFqQ+VnUi9fLCEtSdi2vsWxYTOT2lVTlGXDsE02k0KxxA2Pob5AOMEoHBQbr9wOy5cLraJLLqGQTMD44dZtBPA+AJsZY3+e/nqLD8clnKCEYJYtC/xUZklNu8Mr1HFvoxh4WmJlrfB28KJcmJZmDWg2kzK8HwfXXVu+UEJa0jPDAnX3pxO0F5EwJJF9h7HZkMyZ8klVhH94Nuyc872cc8Y5P59zfuH0191+LI6wyegosHNnKCEYwLx6xI4nqo17b1y9BJLO1lzmQOtCdyVyTj1ixVCqxbqKJY5sWt+4M7AZz11LvlgCY6zsOamft55+ixnaC0EYksiBMH8+MDUldpckORAY1HmadJQQzORkKJPj7Sgq6hlVJXShFzLobK1H9+HhsooWmXP0n8qbnq+gakjSW4vZ81DCGEZhlaIsI50qT8gCsx60UcJV5uKikElL9pqr0ilMlWTdOLteAjg2tepuoJBM4JBhTzohV8HYUVTUK7WzEp3SK1MEjA0nAGy/vL0s3q+cz6wSRvsYs3PkCyVsW99ieRHQfWyxhBs3rTb8uzbBeX/vkO4w7XXLcmWvXaxq1Z2iDslQlUwgkGFPMvk88J//GVoIBrCvqOgkqWe2CzAypIoBs1Pad3/vEJ7qHwUHBwNDWgKKsrMwgCKgdeuew7YfozaydhKdfSPjusfpGxkve3x7Yx0ODIyV6+SU9HVyYodSJXPHHcD110e9moqDDHuS6e4W8fUQa9bteIpOS+3MdgEb2urLDJjaI9ca7fbGOlMvmIOj4LBScc/BoZnac7tIjM2s0aiKaODMBPpGxmeMtVl4Rfv4AwNjWLcsN/N4JZSj7HwSoXPT2ioal668UoRnCN8gPfakMjYG/PKXomU7RIIY2mC2C9jc2WyoJa4YbXW8+8n+M/jW/c/MGOGn+r1XiRRlPuMxu8Eo0flk/xlb5Y9Kolb7+L6RcWy/vB07tnQgk5bK4vOxT6amUmLmLmm3+w557Enl3nuBQsHRiDs/CKKr0WwX0Ds4OiexWlC1dhoZ7aI82+VpNx5uhdF4OyNkzmeakbwkNPUStwp2egH0btcL66h3DqGGcFpahArpc88Bq43zEYQzyLAnkRMngHvuEVvZCPC7q9FI26S9sa6szT5fLAnNmcPDpkZb8VbtJDsV6QIzCYN8seQocQrMGlWzMIsZipHVqxhS/q7+2U4yVS8spA5VhR7CYUzMR/3pT4GbbqKOVJ+gUEwSuesu8d0HOd44oB3dpoQenuofNWyzN6qiUTOWL2BDm7VxSk13qm5oM6/OcOr9KzuOQsm5/AAwm7Cd0hGgUcfwAfshMjsNZKGHcBobgYMHgZ6e8M5Z4ZDHnjQGBoDf/17MlgwRP9vX9WrI9So9vIZRctkMNnc2A4BuGaFCUZbRfXgYxZK/cd58oaQr7JVNp9DRsqAsKaxGucg9cOSE7sVtXkoq6wVQ7m/2HnmVMA4ExoR2++23A1/6koi9E54gw540fvlL4amH+M/v52QhoxpybUjAD5TSv2w6BcYYuEmCzs4OwIxsJgVwzInDG3nqmbSEzZ3NWLawVjd2r/a0jQysXrzfTojMblgo9Hr4RYuEbvvjj1PTkg9QKCZJ9PUBjz4aeiWMn+3rdrVk/EAx1vliydSoeyUtSdjU0YRM2t7HSTGsna31uHHTamxb36Jb9aP8rodbw2tHyiAy7fbGRuG1T02Ff+4Kgzz2JHHXXaIKJqRmJAU/29fj2PKeliSkJeao6kWNYojVwz/M0BplM0/b76EZeiGbSKti1ORywPPPA088QV67R8iwJ4UTJ8Q2NeTYOuBv+7rbChG7FSlOj280ucjJ4xUjaGeNTo1yEOWlsdNqV7Nwoaj46uqiChkPkGFPCn/6k/gesrcO+Os12tU8VyMxZjqEQmvslGlHdlAqTxSc1qurXwM7Fx4rzRw9Ym2I/aahQdS0v/QSQCM0XUOGPQkUCsDu3YGrNxpVvvjpNWqPZQczo84w1zj3Do7OaWJyihMNmWw6Nec1sNotGKu02yNxgzXcwJgoDujuBt7//qhXk1hYkEklI7q6uvi+fftCP29i2b8f+MY3AtWE0VNJVDAzIl6MjRNBLTOUaUlG6ohWZDMpdDQvmNGc0UO7a9BTrDR7DdXHuXKdKMHUG2RtdpsWO6qZiaRQAIaGxP/8/PlRryZWMMYe45x3Wd2PPPYkcM89IrEUIGbVKkbljV7KIHsHRx13cpqtHTCvVTcjXyhZPvbKdc1ztNPBgF09x/DAkRNzko/qmat6yJxjz8EhQKX/ogyyVv5udJsWp0Ou3RDJLiGTAUolUQG2aVOw56pQyLDHnaNHgWeeCTzeaBUW0TMi3YeGDcsgrWR69xw87puOy1i+gN09Q74cSw8lQaq3K9DW36tnrhpVyYhwz9znrme8zUJQ6vMDKNPUyaZT2LSmyZMRtqtKGYixX7wYuPtu4IorIskrJR0y7HHnj38UzUgBVwjYqSZR/713cNQwyWh1HKta9lw2g4JKgtYOfl0ktKiTxL2Do7Z2BUVZxm6bpY9eUWQL9DR1FI/frdE1U6VUsNqlufb4ldLHZ54B1qxxtf5qhi6FcWZ8XMgHNDcHfio7FS7q8kar5qSde/sM9cvNDP+OLR3Yfnk7NnU0uRpk7TdFWcaunmP4zh+OoPvQsO3HmV1m/LpEKxcdI9kBRWHSLXaT20bNap4HbmezwJ49ttdLzOLLJ4cx9gPG2HHG2NN+HI+YpqdHdOGFIPbV2Vov2uIN0JY32vHu9T7Edmamzp5z1gRmMylsW98yp0szTPKFkusGJi0b2haWDbrWG34t6dymoO5QNXsvvDSEOXmd9c7juWN56VIh6Ts2ZnsdhMAvl+g/AWzz6ViEwsMPA3V1oZ3OyEvOplNl1Rd2PvR6H2I7M1MVT09tSBWBrs7Wemy/vD0Zsz0NePromTKJAwbgvOX1c6QFrlzXjCvXNc+5bdv6lpldjZXsgNXfrLAjP2B2Hs8dy5IkBnAc9qd6qprwJcbOOf8jY2yVH8cippmcFK3VIYRhFJzUq9ttNNJ+iO3MTDXz9JT7uGl0CgqJCftjN8qvVyovc47Dx87qDr+2iklvXL1EV0WSAZ5moBrJD5iNKlTjS8dyNgs88gjwqlfZfwxBydPYcuSIKPlKh/sW2e1y1H7orYZOq3+3+rCbeXpqI7Vl7VJHjU5BIXMRNmGwV8liRL5YQu/gqGVFkZ7Oi8z5nEEhaYlB5vA8A1Xv/2HZwto577t6Z6a+ry8dy4sXiz6OqSkxRo+wRWjZKcbYDYyxfYyxfcPD9pNQVcv+/bHXpVbCIju2dGDr+mZbgx7sDISw6tAcyxewq+cYBs5MlEkCRIXMOealpDlhEzeYhar0kpHquakcmCm1rJ2XDmwGamdr/cz7qJZd1uZUtANUtMqVtshkgGJRyAwQtgnNHeSc3wbgNkB0noZ13kQiyyK+viQC6VQDrMrW9Dx4PU/OTrjHbunik/1nTMsPzUbdBUG+WJoTSnHTCauuS9e+Rk6mH/mpyKmHnXAZ4JPOjSSJsGRnp7fjVBEUiokj/f3A2bNiGxoiRsZbr1FlV88x7Oo5NmcC0sbVS8q233ohAKMPu3J+vwjbe8imUzMCZMrr4aYbViu1oLyGdvMJyvn9UuQ0OoeT2z2xZIkQwXvXu0jx0Sa+GHbG2O0ANgFoZIz1A/g853ynH8euSp4Ov2rUTB7AzFPUbsXTErPdjaq+kGTTKUyVZE8x6iiRGMNUSZ6p5FFej2zGWF7ACUVZdiRd7LeOu945grxwzKG2FhgeFl3YEchWJxG/qmKu9eM4xDQPPihGhYWI2dbaSaOKkbCiNvEJzNU/96tGPCpSEisbhVeUZaQl//Ikdoy6YryD0HFXE/SFQ5eeHjLsNqFQTNw4cUJ4JiFrUZttrd0OxzA63q6eY9OefTK9cz2M5psGecFSVCmNdFuC1HEP+sJRRkODcHi2bg3m+BUGGfa48dJLIo4YcizRbGsdRM14JRl1MxgYatJSIAY+Xyjh8LGz/mkUOCTUASC5nPhs5POitp0whQx73Hj++UgSRGZbazfDMaqNbDqFoszLLn4cHFMl2XIKlFvUFwx1XgSA72qPkaJ8JgYHgfZ4lLjGGTLscaO3N3DtdT2sttZq70yd9PQTiTHMS0u+JBvDRGIMm9aI6Va7e4bKYuEy58ikJPBScCqUCopomRY/1B4jh3My7DYhwx4nZBno6wt8BJ4RTrtO/QzPqC8kfk1WCgOJCY0XqwudUQw+TBS1x8Qa9poaoRvz2tdGvZLYQ4Y9ToyMiC67kGUE3GCnWQYQBruhNo2XTk2Y3q9QkmcmEiUJmaNMOyXOJDqUlssBhw5FvYpEEH8LUk0MDIjtZgKwYyB2bOmY+dnKC1drmiSNpBh1IKA687Coq6MEqk3IsMeJEBOnXmdZWpVA2hX/IkTZYrFUnnh1ilUDk8RYsHXmQUMJVNtEP6KGmCWkxKnnyTYw1+o2Ev+KM1E1qjPGsKmjaY5Yllv0hNgUsukUrlzX7Oji3Ts4ip17+3DrnsOmE7FChXOxsyVMIY89LoSYOLUj4GRH9GvgzASe6h+d4yUaef/K/d1op4SBnwGwtCSBMXsJU855mUCaG/GwjMRm3letfg8g3nMlh2Fnd2YmMRFp8rWmRsxB3bgxujUkADLsceHECaBQCCVxatZlqohYaW/Xfqh7B0dxYGBsjlHX1r1r2dzZHFvD7ifpFENH8wLXSdXNnWK4inLRZGBIp8olCxQYgBKffV85+Jxdk3oIx1i+YKvs0a56Y+jkcsKwE6aQYY8LZ84IedKAUHvgZrFYI6Ov/VDH9oMfA/KFEg4MjGHdstxMu79TNnc2zxh4oNyDVlDm1Gpr/5X3olAsF1aTOUf3oWHT9ylU9UYnZLPA0JAIyZDSoyFk2OPC2FhgFTFao+C2SUb9oXb7wberUJh0irKMJ/vPYNv6FtMa90xKspXI1jaQZTMpgJcbdDVm74WVxEGo6o1OSKXEznZykipjTCDDHhfGxkScPQDs1pxboZ5sZFbloh7vpjVabYuyljXtlYTeLFI1hZJsO5atNJAZee9avFQiGekDtTeGN1zdEEkSnxcy7IZQVUxcGBkJbBSeX9tntadtVuWiNBn1Do7i3gNDc6pvjp7OY/H86pldaUcfxiikZYTdC7UIuxljVuXS2VqPdcvKK7Se7D+D73Qfib5C5uzZaM8fc8hjjwsjIyLjHwB+1ZCrt+GdrfW6miTA7IWk+9Cwbnz35MtTntdS6WjfL7f6PGaXFat8SN/IuO7t+WIp+gqZsbFozpsQyGOPCyMjgU1hN6s5V0hLEs5vW2gaQ9Vuw83u2zs4mvjhGVGifm21fQd+YXU8s7/7NRjbFbJMht0C8tjjwqlTgRl2PeXG9sY6wwENgH4t9YGBMSxbWDtzv42rlxh67UnTfIkT2gYvv3IkWqwSoVY7vcgqZCQJOHkymnMnBDLscYBz4PRpoKUlsFPYUW7UlkRq0ZYzDpwxToIqlRtJk+CNA1vWLp3zXgVhQO2MsbMasKKUWobOvHliBiphCIVi4sDEhFB1DCh5agftdt9OnbtZs1Eum8GmjiZImlpjBtFGT+jDwLCr59icFn4jz5qBIZt2/j+Ty2bKLh56dLbWY8vapcZ3iKpqdd480dBHGEKGPQ6cPRtoc5IdnMjwAuYVFQBmQjtXrmueeUwum8HW9S2oSdG/nRHKBVXpEP1O9xFDj52DY9KFzvv2y9ttJz3N7hdZDqWmhkIxFvgSimGMbQPwDQApAN/nnH/Vj+NWDcVi5F10drb76u27VQzdbKCyUVy+mtFr3JI5tzSe3EVTm7rPwA6xa1aSJGCKKqvM8Ow6McZSAP4dwJsBrANwLWNsndfjVhWl6OPQZtt95e/q7bvZhSAjmV+kIu9ejCFhduM6TWzrVVXZidEHBmOx+MzEGT889ksBPMs5fw4AGGM/AfB2AAd8OHZ1IMuRD9gwGmatGPP7e4ewu2cIu3qOgYEhk5IMRakKMjdtkzerpokLa44/j40vPoHc5DjGaurwwMoLcGjpqqiXZRuzihanyVirebihQ4bdEj8M+3IAL6l+7wfwah+OWz3EYAKP2YdXW/rIwVEoGV+IspmUaZu8WXNTHFhz/HlsOfIIMrIwHvWT49hy5BEASIRxV7xpo4YmNzsmu/NwQ4EMuyV+GHa9fXfZp54xdgOAGwBg5cqVPpy2gpDlyGPsgPGH96l+++3jaUkCuH6bfPdhUaIW9xr3jS8+MWPUFTJyCRtffCIRhl2pZCkUyx0GOyEUr9O1AocMuyV+lCf0A1ih+r0NQNmIE875bZzzLs55V1MIwyQSBWORh2LMMIv/blvfMqfqZcvapYYJv3yhNEc7Jq7kJvVb6Y1ujyN7Dh4vex+ymZRlmaMf07VCIeIqsrjjh8f+KIBzGWPtAI4CeA+A9/pw3Ooh5v+kRlK7DEzXyzfTNLEjihU1YzV1qNcx4mM1MVA2tEEV4GMxAAAfy0lEQVT34WHd0tVMSrL0vBOhs8957D8zUePZsHPOi4yxvwWwG6Lc8Qec8x7PK6smJCmQUIxfW+oNbfW6zUgb2vSP1d5Y53hSUkpiKMnxMPoPrLxgTowdAApSCg+svCDCVdnHqNvXzk4ptgM21HAeyqSxJOPLq8M5vxvA3X4cqyoJoOPU68xK7UVhxaJa9J/Kz4xq29BWP2fCj5rDx5xLqsbFqAOzCdIkV8XokctmLC/2ZjXrsYm9cx5pl3YSoMteHJAk32PsXrbUehcFRftlU0eT5eMrQdXx0NJVsTbk2XQKmbRk25NOSxLaG+vK3td7Dwyh+9Aw8sXSjDicdlar0WMjk+6VZQrFWECvThyorfXdsHvZUhvJC+QLpXgm0iqYtE6zl8QYNq1pcjTNSJm/qn1f1d2tY/nCzKxWJSHOwFCUZTzVP+p4IEhgFItiqDVhCHnscWDBAmHYfRzQ66UN3I4Ot9ZLczsIgjBGYgxFnRCVIiNgNAhDj6f6R211txZlGX0j42UNa06HnwfK5CSwJKKu14RAHnscSKeB+fOFJ+ITXtrArYy/9sN8f6/oSCWj7i8pg2s8h3nlkf5j7O8Ix/IFw8oaLZHIQ0xNAY2N4Z83QZDHHhcWLRLyvRnnHxSzpJabZJeVDrd2uo/TChjCGiste+U9DepiakdHPzK9GDLslpBhjwuLFwN9fY4fZlX94iaxpTxGSaqp0ZvuY4WZrgyhj5VhVS7UbqQZlAtCNpPCVFF21Fug9DREWhUjScDCheGfN0GQYY8LTU3AoUOOHxZUQ4lyUbAqcbPjMaYYQ4mxSJqTUoyhlICmKKco78PAmQlXO6Zt61tm3l+7Fwe1KFykSBIlTy0gwx4XmppEUsghQTeUmHn9dqtj8sVS2SSlsKhEo75i0ezc2c2dzVi2sNZRzF27qzN6rLqkMlaaMZyTYbeADHtccLm1jGoIQu/gKHb3DNm6L4M7b50xhjQTMsDELKcnvCfZ1bs6I8nmTWusexYigXNRSUYYQlUxcSGXc9V0EcUQBCWub6fSIi1JpvfTrl0N55yMug7qC7lWtMvpcXbuFXmdLWuXlom5xdaoA2TYLSCPPS643FpGMQTBaj6qNsFmpgu+cfUS7O4ZCnWCUNLJZmbb6a3ei7QkIS0xw25gJSyzZe1SbL+83fe1+k6hIIw6SQqYQoY9LjQ2zk5SchiPVsfBlWTnrp5jgRl5M+/QKMGmt9VXdhVk1J0xVZRn5paavRfK+w+Uv/5qYqfeaMbLLwMrVljfr8ohwx4X6upEN93EhPjZBV6Fv+xiFNdnYLpG3WhXAWBmfdUAg84EGht/0yJzjl09x0yrWXLZTJkHbpZgTUxz2dgY0NkZ9SpiDxn2OLFmDbB/v2vDHpaWttV8VD20xv2BIydQKMm2uhsrBcYYUgxlMgFpSfL1dVB2Q3qlqn6Oy4sExoBzzol6FbGHkqdxYs0a4bG7JCxvrLO13nGyTW8yj53uxkpC5hy189K6U6f8MqzK8QDoTkJqb6wLPdnuK5wDy5ZFvYrYQx57nFi2zJMImNfSRyd62067Wu1qj1ghGTQ6+e31BsVYvqCbE/Hr4qu8Zzv39unu3vpGxrFl7VLfk+2haLUXCkBNDQmA2YAMe5xYtsyTyqNRiMSONxZUfL53cBTdh4d9887VRl1bfeOmvT5sGBhu3XPYcNyg+n5ukspK2M1s9+ZWasKIsHI7OHsWeMUrYjH4Pe6QYY8TNhOoRt6Rl9JHv+Lz6rXZ0SLJplOYLMqOjZgS01fWHiejbmaUldutnm9NRnKs4wLMht3CbFwLbU4qJU5tQ4Y9blgkUO2Ifrn5MPkRn9euzY6XvmlNEwDg3gNDjoxYUZbRfWgYRZnHKgSTliSsW5bzrHipfe3SEgNjzFJMTTHcXnZvTgmt0oYSp7ah5GncsEigmnlHXjDz5Hbu7TPVhekdHMXOvX3Y1XPMkZHNplMYODOB3T3HXEkO5IulWBn1XDaDdctyrma+WlGUuaVRVxtuNwlutxj97/i+O6DEqW3IY48by5ebxhCD8o7MNNjH8oWZumlteEfrpdslLUloys2rGC135XVxuvPwCwYx/g4QF+IwE6Oh7A4oceoITx47Y+wvGWM9jDGZMdbl16KqmhUrxLCNgr6hDso70np4RiihH8WDt2ppN2LL2qXoP5V3tVZAGA51a33UKDXiURh1QDQ3PX10FPceGCorcfQyo1avTFV7zFB2BydOABdfTIlTm3j12J8G8BcAvuvDWggAmDcPuPBC4Mkngebmsj/b9Y7clJ8p8flb9xw2vZ86MeZmp8DA0Nla7ynhWZRlpKX4GPY4JG/1LipFWcbuniHXEhN2E6N+V9qUMTkJdJHvaBdPHjvn/CDn3Pl0CMKcSy8F8vrerB3vyI6XZYaTgddudgpKRQiDN+/LSNgqyXh9TfRQXm83HnwsJAhkWXjqHR3hnTPhhJY8ZYzdwBjbxxjbNzw8HNZpk4nyD2wQ4uhsrcf2y9uxY0sHtl/eXuYpeU2w6kkBa1FXX1jdV4/ewVFsaIuv6FQ2Hc1ugINjx5YOnN8WzOg3p4n20BKjZpw+LYoK5s8P75wJx/ITyRjbwxh7Wufr7U5OxDm/jXPexTnvampqcr/iaiCXA849FzjjLrFo5mVZVbgA1vF2q+oLOzxw5AQ2dzZjxaJaW/cPm6h2A0oD01P9zuPiEmO2JlU58baj0PsvY2wMeM1rwjtfBWAZY+ecbwljIYSG174W+NGPgEWLHD/UbHq93a5Ao7Z3vTit9r529NXH8gX0Do5i8IzzcYCVjN0GJi0MDFeuEzkZ5b0yapRy4m1Hofc/ByVvsHZtOOerEKjcMa4o/8gu5AXMShcB512BdhNjTiYr5bIZ1xU1RDkcfI7x3ba+BYC5Dr5dAk+MmvHyy0BLi5hXQNjGa7nj/2GM9QN4DYC7GGO7/VkWgaVLxYDr8XHHD7VTuhhE8suuoVaMS2I0wFV4SW3mshns2NLhOvyUliRkJOseB/WuLDEj74w4eRLYuJHKHB3iyWPnnP8SwC99WguhhjHxD/3rX7tKGileltKsoiWI5JcdQ53NpLCpQwxJ9lPVMAzSEoPMxSxW54+d9ZSd1O9rhc66Dw+jIFvH/5VdmV5yPXGcd17UK0gcJCkQZy64ACiVZuOMLggz+WXnYlEszT6Xhlp9v2LFotqZY5l5qGFyfttC1M5LO2pAUkoXtZ6yk/i52qh3ttY7UslM0kVTl/FxoL6eRuG5gGLsceacc4CVK0V1TEODq0OEkfxyoimuLrd76ZS+Js7pieLMWLfewdFYNP8cPnbWUaWM2UQpp5K86tCKWWJcS2KmIhkxPAy85z00uNoFZNjjDGPAW94CfPe7rg07oJ/88mswghutGGU0ntnf1ceOA/liCdlMytRj1oZNjF7PDW31jjVylAuiXmJcKXFU7yYSNRVJj1JJ/P+/+tVRrySRkGGPOxdeKMSPpqaE3IAP+DkYwShhqniLRvF9M69TeWzsqmZMnGyrma9qNnc2uxI/U4ZkAPqDwSMrSQyC48eByy7z5NBUM2TY4042C2zZAuzaBbS1+XJIPwcjmNXLZzOpslF2iidpFrppb6wzPbYb7IzOY2BIp4w1z81CMYpRt7sTchJSUT8GMC4/TLQhV8O50IbZvDnqlSQWSp4mgcsvB4pFT0lUNX7qf5jFcZWwhdKer04imoUJ+kbGLY9tB23y0ux4EmOoyUiWmudG51GMul2NHqdSDBJjyQ6tOGFsTOiur14d9UoSCxn2JNDSIkq+TngbpgGIMIyR0JQbQ2ploGTOkUlLZbo2Zt6lYhjd6tAobGirn3Neo+Olpytv3M5lVRKhTjR69KQYzm9baPjezEtJleORW3HqFHDVVVS77gEKxSSFrVuBW27x1IFn1hnqNtmmjfnqYXS71VxOvXhye2Md+kbGbe0uFM9fu1b1cO1sOgUw90ZdvV6nz18vpGIUe88XS74P0NDiV0LdE4WCmEdw8cXhnrfCIMOeFNauFTW94+Omg67NMEpGMjBPHYlWzVAMDL2Do2XHN5I+mJgq4tY9h2eMi1L6qMboXGqM/q6upfcq9qW+IPoxQNpK50f57jbZbYSfCXVPDA2JnFJtPMXhkgIZ9qSQTovt6e23ux7oa2QwOLgvH14jQ83BdY1EZ2s9egbOlNWzF+W5+uHaxynnshpDp2dQ/ay00Xq1RkNQ2hvrbHvbVjo/Cl4HaGi980JJ9i2h7hpZFmWOV1wRzvkqGDLsSeJ1rxMSAy69drseZe/gaFm4YtOaJlsTmADoqjvqGYnewVHDJiWzx9nBKLRkN0FsVUWzY0v50Aej0NGBgTHbnrDeMcwuyHaOqUXPOzci1O7VwUEho0EDqz1DydMkUVsLXHONqPF1gR15gd5BMTdTHXPOF0u498CQrck7na31hl2VWiNhd+CDnnExmy9qJnZlFBbJplNlYllGwzachFYOD511PPREO0jFllSDgwEaTnYtoXWvFovCW3/b28I5X4VDHnvS2LgR+M1vRElYLufooXbkBYwMpsy5bc/Z7s7AS2u82WO1MXl12CGbToFhbq+RxFjZjuT+3iHd+LtZ2aFyUVReP788YbvhGbvHdLJrCa3EcnAQeNObhKop4Rky7EkjkxH6Gd/6lmPDDlhra/thjIwMkTLBSbmY2GnSMTIuTsJK6rXYSZb2Do4aVqfMSxuXHXYfGrYtEuZl2IXXARpGr102nUImLYVfFTM1BUiSkM8gfIEMexK5+GKheHfqlKsJS2aYGVu7hsOsBFIdDza6AKQlhqJsrrmycfUS7O45Nse8senb1XQfGrb0dLW7EbOQhllZpN0KG6/DLvT0eZwc0yjJayePEgiDg8A73wksDGbOazVChj2JpFLAtdcC//zPQkvDx0YOo2oTp52PZiWQaq1wwJ3GycCZiTKflU/frjaAdo2tejdhR8dGi1X+QblgqrVdnNSla6tY1i3LzdTzO/WuIx93p2ZiQhQCkHyAr5BhTyrr1gGdncBLL/kalzRq4nHrzVnVZLsdu2Y07Pmp/lFs7pyd/ekEZTdhpuJodHEzO1c2k5oT93daM653/wMDY770HkTO0BBw3XWuezMIfciwJxXGgHe9C/jSl0T9r4fWey1+fuj9aNrRw6jyRn27m1K9oiwjLaV0yx3Pb1to+LqYnWtTR9Oc352KsLkRbYtFF6kVY2PA4sWiIIDwFSp3TDKrVwtp04GBqFdiSFATnIw0VdS3u7145IulMh2XbetbZnYCemQz+qWR2UyqzKA6lR5wersTMbLI4FwM0rj2Wt/kqIlZyLAnHeWD8fLLUa9EFz2xKz8GKrctyurevqFt9rhuRcScXhB6B0cxVSxP0EqMlXnrbjCrvdfDiRhZZAwMAJdcAnR1Rb2SisRTKIYx9jUAbwUwBeAIgOs556f9WBhhk4YG4PrrgW9+E2hvj40iXpChgN7BUQyemSy7fcWi2jletVknrBntjXWOYuBGtf9+KTIaJbSnSvKMBo+d8YSxmYH68stCIuOv/zo2/6+VhleP/T4A53HOzwdwGMDN3pdEOKarC7j0UuDo0ahXAiD4UIBR5+TpiWLZbZ2t9di6vtmR5+60W9TIYBpV5BiFkQD96prO1nrMS5evXynT1L7eRsRiBirnImH6gQ/4XqpLzOLJY+ec36v69SEA13hbDuEKxoT3c+CA8Ibmz490OX5OaNLDLOb8ne4jYj6pSo5XrzywvbHOWCLXoCJGfV61h+y0Ychs5qnRzsBsTXYkAmIzA1UJwVx6adQrqWj8rIr5IID/8fF4hBNiFJLxc0KTHmZNVIqXrPaWjcoD3cwdVTxqdajGqb69Ei7SO7/RBdCOnK8RsamKoRBMaFgadsbYHgAtOn/6LOf8V9P3+SyAIoAfmxznBgA3AMDKlStdLZawQAnJ/PnPvs1HdUNQJY4KdrVT1OgZTDdzR3f3DKEmra/8qHjuWkOql28wG2itt6aNq5dgV88x3fub7Rj0tOwjQQnBfPSjFIIJAUvDzjnfYvZ3xth1AK4G8EbOjYUyOOe3AbgNALq6uvwZ3knMJSYhGaOWdb9CAYrBNDJ0RmhDKXrzTSXGMC8lGcbHObjp37RyvmbNSE4ugJ2t9YbPl4OX1d3HJvSiQCGYUPGUPGWMbQPwaQBv45yPW92fCAElJDM0JGRQIyCoEkftOZzuAJT7K8bWKG7d0bLA1ZrsDvZQdg9Oa/yNnq92WHcQr7cnxsaEeB2FYELDa4z9WwBqANzHxBv2EOf8Rs+rIrzR1QVs2wbs3g2sWhXJhymMlnWz8IQWtcE0SzbKnKNvZBznty10FIN3OthjLF9wrNlithOKjUSAlslJYGQEuOkmCsGEiNeqmFf6tRDCRxS5gYEBoLc30ni7U+zUv6vvoyhBaklLDOmUNFMVoz6OnVmpmzubsWxhrWENfDaTQiZlLXFrZ2B3IsW77CDLQH+/KG1cuzbq1VQVpBVTqaTTwIc+JLRkRkaAxsaoV2SJHXEs7X2KMhdx8bS+Edcjm06Zqj6qjS4AfYnbDnuiaH7nG2LrmWvhHHjxReCNbyTlxgggw17J5HLAjh3AF78InD0LLLAfO7bbOepnh6md+ne9+8icI5OScOPrV9t6XlM6SVMFrdH16iUnzsv2i4EBoKMDeO97Ka4eAWTYK53ly0WJ2de/DtTUiCSWBXZlZZ3Kz5qdz6wdXn271xp5q1mpekbXq5ecGC/bL06eFE7ERz5CAl8RQSJg1cAFF4hxei++KOKeFtgVkfJDbMpOO7y6GsSsMsQOVrNSq8oAB8H4uCi13bFDVGgRkUCGvVrYtg244goxmMMCu16xHx2mVu3w2tCIVxlgrxcGwoRiUYy5u/FG4Jxzol5NVUOGvVqQJOD97wde8QpLsTC7xs8PI2nlqWvrsb3WyAelD1/1lErACy+I2aWXXBL1aqoeirFXEzU1Yov8ta+J5NayZbp3s1vJYSSk1d5of8yZWTmgUTu8l5h1nJOZiZh6pEepBDz/PHDVVcDb3x71agiQYa8+cjng4x8Xg7AHB4HW1rK72DV+fSP6zcZGt+sRtPyAHnFMZvqViA4dWRae+pVXit4JqoCJBWTYq5GFC4FPfAL46leBY8eAlnKNNzvGz48Ye5w96DAJWuo4EGRZeOqbNomyRh/n7hLeIMNerSxaBHzqU8JzNzDuVvil4hhHDzpsgpY69h3FqF9xBXDddUBKf0wfEQ10ia1mliwRGh6LF7uavkSJSP9IVLVOqQT09Ymu0uuvJ6MeQ8iwVzuLFgGf/rTw2B0a9zBUHKuFxFwklUTpm98MvO99ZNRjCoViCBFz/9SngH/9V+GJrVxpOwlGYRR/SESuYXJSiHq99a3AX/4lJUpjDBl2QpDLCeP+wx8CDz4ojLsN+QHCP2J9kRwdFVIB27cDr389GfWYQ4admCWbFYqQK1YAP/0p0Nwc+WBsM5JU952ktZYxNCQqXj7zGSHsRcQeMuzEXCQJuPpqIR727W8D+Tx6pzKxM0pJqvtO0lrnwLmQoFi+HPi7vwOamqJeEWETMuyEPhddBHz+83jm5i/jqaefwVjdIoCxUI2SmZebpLpvN2uN3MMvFoVo3KtfDXzwg0BtbXjnJjxDVTGEMW1t+Mjqq/Hswla0nTkOaWa4hTMFRzdoVR+VC0rv4OjM73rEse7b6VqtnnvgjI/P6r58+MNk1BMIGXbClGfzEnZe8g7sXXUh2s4MIVuYBBC8AbWSBE5S3bfTtfohh+ya4WGRJP3Yx4TuC3WTJhJ61whTljXUophK4871b8B/X/QW5CbHsXTsBHI1wUbxrLzcxNR9w/laI9mNFAqi1HXxYuALXxAD0YnEQoadMOWTW9egNpMCGMP+tnX459dfhxeaV+FNteNiyx4QVl5ukpqjnK419N3I8LBQ+7zmGuD//t9EDT8n9PHkdjHGvgzg7QBkAMcBfIBzPuDHwoh48I6LlgMAvrb7EAZOT2BB61Ksve5zWDnxAvBf/wWcOiXkf32ua7aj+hjrum8NTtYapOKlOinbkGF4w4IpnHNhp+g+XrnS8/GJeMC4wfxHWw9mrJ5zPjr9898BWMc5v9HqcV1dXXzfvn2uz0vEhJMnhXF//HEhSVBnX4fdDpFXhkRIEM9dXXa5aPwM6gqTuH/d5dj28Q/g7Zeu8mfhRKAwxh7jnFvGyTx57IpRn2Y+APdXCSJ5LF4skmwPPhiI954kj9xvgnjuDxw5ARSm0DZ2Akfrm/Ddy67BQP1SPHn/c2TYKwzPGTDG2D8AeD+AMwDe4HlFRLJgDNi4EVi7Vhj3/fuB+nph9KntPD6USpg/PIgFnOOezo34Q/urUEyJj//A6YmIF0f4jWUohjG2B4CeWPdnOee/Ut3vZgBZzvnnDY5zA4AbAGDlypWveuGFF1wvmogpnAO9vcBPfiIqLJYsEQJjRHRwLvT2Jyfxj2eb8PNlF+Fk3dz3ZHlDLR64aXNECyScYDcU4ynGrjnhOQDu4pyfZ3VfirFXOLIM/PnPwO23i4qLpqZYa85UJJwDIyPA2bPAxRcD73wn7hxmuPmOpzBRKM3crTaTwlf+YsNMkpyIN6HE2Blj53LOn5n+9W0Aer0cj6gQJEkYkw0bgIceAn72M2HgW1qE0BgRLKdOia81a4B3vxtYvRpgDO+YrmJUKpyWNdTik1vXkFGvQLxWxfwCwBqIcscXANzIObec1kAee5WRzwPd3cCddwpN75YWoKbGl0NXc+XMHDifldZdvhy49lrgvPMoz1FhhFUV804vjyeqhGwW2LZNJFl/9zvg3nuBiQmRZG1ocG18Equa6CelEnD8ODA1JS6YH/mI6BqlyUZVDak7EuGRywHveIcYq/bEE8A994gkayYjtN/Tzv4dk6Tw6Dvj4yK8xRhw2WXA5s0zIReCIMNOhE9NDXDppcAllwhp2O5u4H//V3ifS5YACxbYOkySFB59QUmIjo+L3c573iNkdRsaol4ZETPIsBPRwRhwzjnAddcJidhHHwXuvltIxmYywsibxOJz2YyuEY+jwqNrOBeVLadOiWqj884Dtm4VfQMOdzhE9UD/GUQ8WLAAeMMbxDzNw4eFkX/4YTGWDRANT/Pnzwk1BKmpEimyLAz52bPi92XLxIXvVa8ScXSCsIAMOxEvJAno7BRff/VXIlTz9NPAAw+InwERq29omDNNKfFVMVNToqJlclJcvNauBV7zGvE6NDZGvToiYZBhJ+KLJAGrVomvq64S8eWDB0Vt/KFDAOfoBNC5vkF4/JmEhGA4FyWgZ8+KeDljwLx5IufQ1QWce67vgmpEdUGGnUgGjIkO1qYm4IorgJdfBo4cEfH43l7x89SUuB/nwqu3aewDrYXXGnFJEqGWRYtEA1dHh5DLXbWKYuaEb9B/EpFM5s8Hzj9ffL31rcJYnjgBDA7OGvvnnpsNbShx+HnzxFdNDTBvHnpHJrCnd9h9LbwsiwvK1JQ419SUmEYkSbPnVRvxtjYRM7dZ+UMQbiDDTlQGkjTr0auN/cmTooFnbAw4c0bUfo+MiNtPnsSzjz+DpVMlyGBQ0rIcwDOPnUTnBa3TN/DZpK3ynXNxfEkS5YaLF4sqnqYmERPP5YQAWksLGXEidMiwE5WLJAkja5J8/Ninfo26Qh7zJydQW5yExGUwzpHiMt76sdfNGm/lK5USXwsWCONdV0dNQUTsIMNOVDUtixfg6OkUxmrmqk8ub6gVnj9BJBAaZk1UNTPDulXUZlL45NY1Ea2IILxDHjtR1WiHdZOULVEJkGEnqp53XLScDDlRUVAohiAIosIgw04QBFFhkGEnCIKoMMiwEwRBVBhk2AmCICoMT8OsXZ+UsWGI4ddJpBHASNSLCJFqer70XCuXSnm+53DOm6zuFIlhTzKMsX12poRXCtX0fOm5Vi7V9nwpFEMQBFFhkGEnCIKoMMiwO+e2qBcQMtX0fOm5Vi5V9Xwpxk4QBFFhkMdOEARRYZBh9wBj7BOMMc4Yq9gx8oyxrzHGehljTzLGfskYa4h6TUHAGNvGGDvEGHuWMXZT1OsJCsbYCsbY7xljBxljPYyxj0W9pqBhjKUYY/sZY7+Nei1hQYbdJYyxFQDeBODFqNcSMPcBOI9zfj6AwwBujng9vsMYSwH4dwBvBrAOwLWMsXXRriowigA+zjlfC+AyAB+t4Oeq8DEAB6NeRJiQYXfPvwL4FMSIzIqFc34v57w4/etDANqiXE9AXArgWc75c5zzKQA/AfD2iNcUCJzzQc7549M/j0EYvIrVLGaMtQG4CsD3o15LmJBhdwFj7G0AjnLOn4h6LSHzQQD3RL2IAFgO4CXV7/2oYGOnwBhbBeAiAA9Hu5JAuRXCAZOjXkiY0KANAxhjewC06PzpswA+A+DKcFcUHGbPlXP+q+n7fBZiG//jMNcWEnrTqCt6J8YYWwDgFwB2cM5Ho15PEDDGrgZwnHP+GGNsU9TrCRMy7AZwzrfo3c4Y2wCgHcATTEynbwPwOGPsUs75sRCX6BtGz1WBMXYdgKsBvJFXZn1sP4AVqt/bAAxEtJbAYYxlIIz6jznnd0S9ngDZCOBtjLG3AMgCqGeM/T/O+V9HvK7AoTp2jzDGngfQxTmvBIGhMhhj2wDcAuD1nPPhqNcTBIyxNERi+I0AjgJ4FMB7Oec9kS4sAJjwRn4E4CTnfEfU6wmLaY/9E5zzq6NeSxhQjJ2w4lsAcgDuY4z9mTH2nagX5DfTyeG/BbAbIpn400o06tNsBPA+AJun388/T3u0RAVBHjtBEESFQR47QRBEhUGGnSAIosIgw04QBFFhkGEnCIKoMMiwEwRBVBhk2AmCICoMMuwEQRAVBhl2giCICuP/A4N/WtzCoIz2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe4f80008d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 596 ms, sys: 96 ms, total: 692 ms\n",
      "Wall time: 507 ms\n"
     ]
    }
   ],
   "source": [
    "def minimalEnclosing(P):\n",
    "    M, R, X, S = basicModel(1, P, 0.0)\n",
    "    M.constraint(S, Domain.equalsTo(1.0))\n",
    "    M.objective(ObjectiveSense.Minimize, R.index(0))\n",
    "    M.writeTask('p1.task.gz')\n",    
    "    M.solve()\n",
    "    return R.level(), X.level()\n",
    "\n",
    "%time P=gaussianData(1000); r, x = minimalEnclosing(P); display(P, r, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 2. Maximum coverage\n",
    "\n",
    "**Objective: find the largest subset which can be covered by a ball of fixed radius.** \n",
    "\n",
    "This is also a well-known and multi-facet problem, see [this paper](http://www.cs.princeton.edu/~chazelle/pubs/CirclePlacement.pdf) by Chazelle and Evanston. It can be interpreted as selecting a location for a single fixed-power transmitter so that the signal reaches as many locations as possible. In the combinatorial setting it is equivalent to finding maximal cliques in unit-disk graphs. The best known algorithms for this problem achieve running time $O(n^2\\log{n})$. \n",
    "\n",
    "In this setup we still operate with $k=1$ and we bound the radius from above. Each $s_i=s_{i1}$ determines if the $i$-th point is covered, with the optimization objective equal to $\\sum s_{i}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XucXHV5+PHPs/dNsrlvks3mCgRiIGDMlouhiogSbxApykUKtrRUf1JbabH4Qy2lLwWlP8VfS6tRqLf2h0oVo6QGCQhykyxNABMSCLnAJkuyLLlsrnt7fn98Z5LJZGZnZs85c86Z87xfr31ld+ZkzvfszD7ne57vc75fUVWMMcYkS1XYDTDGGFN+FvyNMSaBLPgbY0wCWfA3xpgEsuBvjDEJZMHfGGMSyIK/McYkkAV/Y4xJIF+Cv4gsFpENIrJRRG7K8fzHRaRLRNakvv7Mj/0aY4wZnhqvLyAi1cBdwHuADmCViCxT1XVZm/5IVa8v9nUnTpyos2bN8to8Y4xJlGefffYNVW0utJ3n4A+cCWxU1U0AInIvcDGQHfxLMmvWLNrb231onjHGJIeIbC1mOz/SPq3Aaxk/d6Qey/ZHIvK8iNwnItN92K8xxphh8iP4S47HsmeL+wUwS1VPBx4CvpfzhUSuE5F2EWnv6uryoWnGGGNy8SP4dwCZPflpwPbMDVS1W1UPp378NrAw1wup6lJVbVPVtubmgikrY4wxw+RH8F8FzBGR2SJSB1wOLMvcQERaMn68CHjRh/0aY4wZJs8DvqraLyLXAyuAauAeVV0rIrcC7aq6DPi0iFwE9ANvAh/3ul9jjDHDJ1FdzKWtrU2t2scYY0ojIs+qaluh7ewOX2OMSSAL/sYYk0AW/I0xJoEs+BtjTAJZ8DfGmASy4G+MMQlkwd8YYxLIgr8xxiSQBX9jjEkgC/7GGJNAFvyNMSaBLPgbY0wCWfA3xpgEsuBvjDEJZMHfGGMSyIK/McYkkAV/Y4xJIAv+xhiTQBb8jTEmgSz4G2NMAlnwN8aYBPIl+IvIYhHZICIbReSmIba7VERURAquLG+MMSY4noO/iFQDdwHvA+YBV4jIvBzbNQGfBn7ndZ/GGGO88aPnfyawUVU3qWovcC9wcY7t/hH4KnDIh30aY4zxwI/g3wq8lvFzR+qxI0RkATBdVX/pw/6MMcZ45EfwlxyP6ZEnRaqArwN/U/CFRK4TkXYRae/q6vKhacYYY3LxI/h3ANMzfp4GbM/4uQk4DfiNiGwBzgaW5Rr0VdWlqtqmqm3Nzc0+NM0YY0wufgT/VcAcEZktInXA5cCy9JOqukdVJ6rqLFWdBTwNXKSq7T7s2xhjzDB4Dv6q2g9cD6wAXgR+rKprReRWEbnI6+sbY4zxX40fL6Kqy4HlWY99Mc+25/mxT2OMMcNnd/gaY0wCWfA3xpgEsuBvjDEJZMHfGGMSyIK/McYkkC/VPsaUjSocPgw9Pe5r3z7375tvwu7dMDAA/f3uX4CaGqiuhtpaGDfOfY0aBU1N7mvUKKivD/eYjAmBBX8TXb298PrrsH07vPwyvPSS+7m/H6pSF62q7quqygV4kaNfmc8PDkJf39FtRdxjg4PQ0AAtLTB3Lpx4ovt+8mR34jCmQtmn20TH7t0uwG/Y4P7dtu1okK6rcz31KVNcT94vqu4qYc8eeOghWLHCPV5VBdOnuxPCySfDnDnuKsGYCmHB34RH1fXq162DJ56ArVvd4+lAP23a0R5+UERcD3/0aPeVlnlC+NWv3HYnnwznnANveQs0Nx+9ujAmhiz4m/IaHIRXXoE1a+DJJ12AVYWxY2HGjOgE1OrqY08Ig4PuRPXd77r2TpkCb387vPWt7gohKu02pkgW/E157N8P7e2wfDns3OmC64QJLuDHQVUVjB/vvlTd8dx/P/z0pzBrFrz//XDGGTZ4bGLDgr8Jjip0dMCjj7qvvj4X8GfODLtl3oi4/P+oUe4Yd+2Cf/s3aGyECy6Ac8+FSZPCbqUxQ7Lgb/w3OAjPPw+/+IVL8dTWumBYWxt2y/wncrSE9PBheOABd9zz58MHP+gGii0lZCLIgr/xj6qr0rn3Xti0CcaMcb38pAS/+nqX/x8chI0b4UtfcqmgSy+NT3rLJIYFf+OPrVvhvvtcj7+pyeXBkxL0s1VVuSud9MnwC1+ARYtgyRJLB5nIsOBvvNmxww18PvkkjBiR7KCfTcRVBQ0OwqpV8PTTbkzg/e931U3GhMiCvxme/n53Q9R//Zer3Jk5M/ia/LiqqoLWVvc7W7nSDX7/8R+7UlH7nZmQWPA3pevogO98BzZvdkGtri7sFsVDTY0bEzh4EJYuhWeegWuucRVQxpSZdTtM8fr7XTXLF7/oavVnz7bAPxyNje53t349/O//DY8/7lJDxpSR9fxNcTo64O67XRWP9fa9E4GpU49eBaxaBVdfbVcBpmys52+Gpup6pn//925w13r7/kpfBbz4Itx8M/z+92G3yCSEBX+TX3+/q9lfutRNZGZlisFIXwWMHAl33AEPPuhOusYEyJfgLyKLRWSDiGwUkZtyPP8JEXlBRNaIyOMiMs+P/ZoA9fTAnXe6GS1nznRz3ptgNTW5lNoPfwj33OPWMzAmIJ5z/iJSDdwFvAfoAFaJyDJVXZex2X+q6jdT218EfA1Y7HXfJr/7V2/jjhUb2L77IFPHNnLjhaewZEFrcf95+3b4xjegu9vq9sutrs79zn/7W+jshE99yk0dYYzP/BjwPRPYqKqbAETkXuBi4EjwV9W9GduPBOyaNkD3r97G5376Agf73FKG23Yf5HM/fQGg8AnghRfgn//ZzcMzbVrQTTVZ1nfu5YlXuuk52MusdU9z2rotnPTlz7sTgjE+8iPt0wq8lvFzR+qxY4jIp0TkFeCrwKdzvZCIXCci7SLS3tXV5UPTkumOFRuOBP60g30D3LFiw9D/8Zln4J/+yc1h39wcYAtNLus79/LQizvpOdQHImypG82TL3ex5Yab3TKWxvjIj+CfKydwXM9eVe9S1ROBvwM+n+uFVHWpqrapaluzBZ9h2777YEmPA24lrX/5FzcdgS1XGIonXummP6ve/836kTy18zB85SvuvgBjfOJH8O8Apmf8PA3YPsT29wJLfNhvYO5fvY1Ftz/M7JseYNHtD3P/6m1hN6kkU8c2lvQ4jz0G3/qWqzhpzLONCVzPob6cj++gzs2Q+tWvwtq1ZW6VqVR+5PxXAXNEZDawDbgcuDJzAxGZo6rp69YPAKFcwxYzCOopXx4RN154yjHHANBYW82NF55y/MZPPummamhtjdUqVEdy44f6aGqoZdGJE5jbMrrwf4ywpobanCeApoZaVwkkAl/7Gnz2s3BKjvfSJ56KBSKuko+tVJ57/qraD1wPrABeBH6sqmtF5NZUZQ/A9SKyVkTWADcA13jdb6nSQX3b7oMoR4N6dq9+2PnyCFmyoJXbLplP69hGBGgd28htl8w//kO+ahV885uuxx+zwH8kN47rMT/04k7Wd+4t8D+jbdGJE6jJmuitpqqKRSem7vodNcrNBnrHHW69gAAU+3cSR5V8bMMhGtGbSdra2rS9vd2311t0+8Nsy5Hzbh3byBM3nX/k59k3PZCzFEmAzbd/wLf2hG7DBrjtNpg8OXapnrsf35y3h3ztubNDaJF/irqi2bMHDhyAW26BlhZf91/s30kcVfKxZRKRZ1W1rdB2iZnbp9hB0KljG3N+QPLmy+No5053A9f48bEL/JA/N57v8Vyimjaa2zI6Zzuy23veeDjxzjvdQjE+DtAPq1ggJir52IYjMdM7FDsIeuOFp9BYW33MY3nz5XF04IC7gQtcSWcMNTXkXgs43+PZ4pY2ytXe/359gC0btrqpN/r7fdtXycUCMVLJxzYciQn+xQb1ovPlcTQw4Gbm7Ox06Z6YKpgbLyBXSWX/4CBPvNLtWxv9lK+9K3vqYM0at6COT+nbSu78VPKxDUdi0j7p4F3MSP+SBa2VEeyz/eIX7kauE04IuyWepNMipaRtMtMm+ZSSNiqnvGmuw/0w80S3xsL06W5lMI9K+TuJG7+OrVIqhhIz4Jt47e0u3TNzpltRKkHSaZPs3nO2qA4YFxzgPngQurrg85+P/Yk96rJLwcFdPUQpO1DsgG9i0j6J1t0N3/62qwxJWOCH3GmTbJlpo/Wde7n78c3c+dBL3P345tDHAgqmuRob3X0Ad93lTgQmMJVQCp5mwb/SDQ7C97/vcsIjRoTdmlAUSuc0NdRywVsmMbdldCQHg+e2jOaCt0w6MqCd2d4jxo2DN9+En/0spFYmQyVVDCWvG5g0Tz0Fq1e71aISaqg7Z7PTPEMNBodZCpqvBPQY06a59RcWLgz0DuAkq6RScOv5V7Lubtfrb2mJ3Zz8fqZeSqkO8uMegtBUV7t7N5YutfRPQCqpYsh6/pUqM92TcSNXVG9uypQ9QJtOvQDDamsp1UFDzq8TB2PHwtatLv1z5ZWFtzclqaRqKAv+lSpHusfvoBqUIFIvRaVNcFcJ2ZVBpdxDEAmW/glUpZSCW/CvRPv3u3Vgp0w5Jt0T1Xx2tjBTL8O5hyAqMq/qpuohzvjKP3PKt7/h0kEVKk4191FrqwX/SrRyJRw+fNxdvHHJZ4edein2KiFKsq/qtksDte3r6fzhcs675kMhty4YcZp+PYptteBfYR54bB0D//hNXqsZScOrm4/ptXoJqkGNFeR63YpIvZRZrqu6roYmepd+F6640C0MX2GGqrmPWvCPYlut2qeC3L96G4987XscONxPX3XtcTXqw50TJ6ja93yvCxSuazfHyHVSP1DXSNWe3fD44yG0KHhxqrmPYlut519BvnPfU1y96X/YMWr8kccyc/rDzWcHNVYw1Otee+5sC/YlEATNsRJF98hx8JOfwFlnwciRIbQsOHGquY9iW63nX0HmP/sbBqSagapjB/gye4VzW0Zz7bmz+esLTi46wAY1VhCXMYg4yBX4AQ7V1sGhQ24cqMLEqeY+im214F8pOjo4t+tldo4ae9xTXgdKvc6fX+7XTaIhf5ctLW5G156eMrcqWOnp18eNOHrs9TXRDGmZU8UDVIscyfmHtYxkNH9TpnS/+Q2nzphAdfWxmTw/Bkq9zp9f7tdNoly/S4DZE0e4wd7+fjeddwjuX72NRbc/zOybHmDR7Q/7HuwO9R1NHe4+2BfZdXmXLGg9cgUwkJpNOcx1hC34lyjoD/Kw7N8Pjz7KrHmzAxkoLWpisQi9bhLNbRnNvKlNxz2+bnuPG5ifMMHN+z8wkON/ByfoRdPjNstmlNprA74liGKtLuDm6u/rg9pa5rbUBhI8sweL06te+XECsGDvj81vHDjusSMD8+fOhi1bYMMGmDevbG0KusQxilU0Q4lSe33p+YvIYhHZICIbReSmHM/fICLrROR5EVkpIjP92G+5RemsfcTgICxf7np2AYriVMd+iNrc/V4MNYB+9+Ob2bR/EB58sOjX8+MqN+hgF7d1ecc05h6bCaO9noO/iFQDdwHvA+YBV4hIdtdiNdCmqqcD9wFf9brfMETprH3Exo2wY4dbzCNAcVv3thiVcELLPHkJ+WdudYu+97P14Sfdql8F+JWuCTo4R7GKJp/7V29jf2//cY/XVkko7fWj538msFFVN6lqL3AvcHHmBqr6iKqmr0mfBqb5sN+yi2Qv46GHoKEh8N1UYllm3E9o2SevfOWeaX0KT2/Zxd/d8M2CQdyvq9wgg3N6rpyDfQNUp+awah3bGKklFTPdsWIDfQPHv0ejGmpCaa8fwb8VeC3j547UY/lcC/y3D/stu8j1Mnp6XL6/uTnwXVViWWbcT2j5lqcc6gqge8QY5q39HTfft2bIE4BfV7mZJY6Cf8E588oEYED1yN9iMa8dRuFGvt/d7gPhfN78GPDN9UnL2QURkauANuCdeZ6/DrgOYMaMGT40zZtcs/Dddsn86MzMt2GDm6/fx1kb883hU4nz7YQ9gZxX+U5SiuY9tt6aOiYe2MP47teHHHT1847UIKZA9jKQHFbhRtTu8vWj598BTM/4eRqwPXsjEbkAuBm4SFUP53ohVV2qqm2q2tZcht7sUPLlPAGeuOl8Nt/+AZ646fySPyy+9jieeeaYhVq8GioHXollmXG/z2Coq7F8df8Ag8ApXVuG7MVH7io3i5crE78LN4r9m47a79SPnv8qYI6IzAa2AZcDxywhJCILgG8Bi1V1pw/7DFwQJWq+9jh6e2HNGpg0aVhtyaXQHD6VVpYZ57n7YeiFZ7KPLdPehlG8bduLvNh2Xt7XjvqKVV560X4WbpTyNx2136nn4K+q/SJyPbACqAbuUdW1InIr0K6qy4A7gFHAT8QNzLyqqhd53XeQgqjs8fWEsmmTu2uzxr9bNeKeAx+OUk9oUVoGs9DJK31s2XP9H6htYMa+bm4+c+KQrx/lFatuvPCUY4IuFN+L9jP9UurfdJR+p75EDlVdDizPeuyLGd9f4Md+yimI/JyvJ5Q1a3xfoSnuOfC0INceiNoymMWcvI47STTWsXDiOGbW7C5HE32VOQ43dkQt9TVV7DnYV1Iv2suJI1sky7+LZHf45lHqB6SYJdp8O6EMDro1esePL7xtCSphUNfvAJ15Isk1bXIUl8HM5biTxJtvus/QO3PWXkRSdopl14E+Gmur+fplby2pN+1n+iVqg7ilsOCfRykfkGLzfr71OHbscGWe48aVelhDSgeH37zUxaFUG2uq8pcNRsH6zr3HtDeX4Qbo7BNJvjr6INJigaeXxo6Fl16Cgwd9LRrwolAHys+0qV/pFz+vIsrNgv8Qiv2AFPuh9K3HsX27K/EMSH/GjSiH+gdCT23ks75zLw+u28FgEb+L4QTofHX02fxOi5UlvZSuBHr9dZg925/X9KCYDlQUUyxRG8QthQV/H5TyoSx0QikmfcQrr/g60JspqFW7gvDEK91FBX4YXoAu5oQRRFqsbO+BKmzbFongX0wHKqopligN4pbCpnT2gV/TPhQ9n8qGDYHN5ROnip9i2zTcAJ3vhCEZz/t5r0N6np6yvQf19S71EwHFdKCiVicfdxb8feDXh7Kom0/6+2Hr1sDWY43TNA5DtSk9xYGXAL3oxAlUyfFjHiLC4lOn+LrOcPYNdrn4/h40NcHLL/v7msNUTAcqqKkiksrSPj7wK+9XVPpo507fp3TIFKeKn0UnTsiZ868S4b3zJg87MGcOtuYyqOp7CqbQ+EIg78GIEfDqq26N3zJMDjiUYgdO45piiSIL/j7x40NZVE6zs9OVegYkTne95qpOaqip5rxTmj0F/uyTXy5+p2AK9fgDeQ9E3FdnZ+h5/zgPnHpV1DhfACz4R0hRvZ+tWwPr9afFaRoHv9saVoXPUDfYXXtugIFZNRLBH5LZqw9zdUDL+UdIUTnNzs7QL9ErWVgVPqFNMldVBd3xWL+gEoW5OqD1/COmYO+nuxvq6srXoAqQ64YpyJ3aytcDT9/dG1QKplC6LbCbvurq3DiSCUWY9y5Y8I+bXbus51+CXDdMPbhuB8CRgeLMm6jyDXiXY/rqfCmsQG/6qq93Uz2YwAyV0w/z3gUL/nEyOAh79sDoeOTjh8tLLzf7//b1Dx6Xw891Y1j6Jqp0fj1KA96B3vRVV2dpnzz8GIgtlNMPc3oIC/5xsn+/+zdH7Xml8NLLzfV/S5HePmoD3oHe9FVX5xZ0V63oz1Wp/BqILXTncphVThb846Snp+Q/0CjNP18ML73cYit18onijWzg71TbOT8Pvb1ucaD6ej+aWxH8mkSumJx+WFVOVu0TJ/v2lRT8h1qWMaq89HKL7QlXiRx3525Ub2QD/yqB8n0etu466DoW5gi/BmL9mvolCBb8QzKstXx7e0uazXOoXnRUeZleIt82DbXVx6w9/N55k3nvvMmxWY/Yr/WT830ent+21322zBF+Be0oz0dUsWmfsO6aK8Zw8on3r97GT7/7FB9Y9Rr7tvYXlb6J0yRtaV6ml0j/3xNf38SiV5+j6fABeupH0PP2d9B67lnHbR/VYJ+LH+MQ+d73/b39gd41Hkd+DcRG+c7ligz+Yd41V4xS84np4zlhr7vkLHYQNI7LMnqZXmJuy2iaXn6RSa88Q+2g+/2OPnyA0Y+vhHEjYP78QNsedfk+DyPqamAg/2I4QYhy5wz8DdpRvXO5IoO/rwulB6DUfGL6eEQVUqtJFTMIGqdJ2jJ56eW2/s9TMJgVyPr6YOXKxAf/fJ+H01vHlLXnH/XOWVpUg7ZfKjLnH8UVfzKVmk9Mt1tFODqbfOH0jV+54ljZs6e0xxMk3+dh1sSRR1f2KoMwpzQwR/nS8xeRxcA3gGrgO6p6e9bz7wDuBE4HLlfV+/zYbz5RXfEnrdR8Yvp4BkUgYx3Z7PRNvrJOL8E+bqWijBmTM9D3jWoiusmu8sn5eXhtT1mDf9Q7Z0nh+R0XkWrgLuB9wDzgChGZl7XZq8DHgf/0ur9iRHmEHUpflCJ9PINShaZ6/tnpmyDKOuNYKsq7381g1hKXfVXVPNRyGnc+9BJ3P7452u0vs/Wde/nFc9s5+/aHi6868yjs8sdhVdpVID96/mcCG1V1E4CI3AtcDKxLb6CqW1LPlSWxGOUR9rTMfGJ68OszP1qTs63p73/8g04g9/zuQUwBEKf1fI+YP59H1+9k4Sv/c6Ta54kZZ7Bh0iwgoMXQYyp9cp/YO0B/VTWvlyn3HuaUBnEZbygHP4J/K/Baxs8dwPF1dWUWl8GaYj+MSxa0smTiBbD3CZg27bjXCaKsM46logDPjZvOc23T8z4f+RNYmaRP7lU6yP5a1+suR2GE352zUiqHol4MUk5+BP9ct5wWfydS5guJXAdcBzBjxgwvbYqkXB/Skj6MTU15qzKCKOuMY6ko5G93pqifwMqh51AfooOoVHGw9ujUDuXIvfvVOSu1J2/jDUf5McrTAWR2s6YB24fzQqq6VFXbVLWtubnZh6ZFR/pDum33QZSjH9JcA9OQ58M4apS7wzfHXb5BLAYS2gIjeazv3Mvdj28umLvP1e5sUT+BlUNTQy11A/3sbRh5zLQhUSmMKEaplUNhjzdEiR/BfxUwR0Rmi0gdcDmwzIfXrSj5PqTVeebqyflhrKmBkSOhv/+4p4Io64xSqWgpg8/Z7c4Wh3sdymHRiRNo1AF2NzQdeSxKhRHFKLUnH/VikHLynPZR1X4RuR5YgSv1vEdV14rIrUC7qi4TkT8AfgaMAz4kIv+gqqd63Xec5PswDqjSWFtd/ODXuHFw8CDUHh/Y0mV86fLMX619nSde6fZUnjlUqWh2GejsiSPY/MaBQMpCSx18zmx37MpVy2Ruy2jqdo3mvt7xCESyMKKQUsu641AMUi6+1Pmr6nJgedZjX8z4fhUuHZRY+T6krRm5/6I+jBMmwKZNefcT6KpPBfbzfMfR+nq/9+tl8Dlq8/NHyQlj6vjs4nP57B99IOymDMtwKofiUgwStIqc3iGKhvqQlvRhbG6GF1/M+3S5yjOLmTvfz/0WGny23v0w9fW5DkVMWU9++Cz4l4lvH9KZM+Hw4bxPl6s8s9jX82u/Q81TlO9qZ/ueg4GloSpGVRVMnhx2KzyxnvzwWPAvI18+pK2tQ96KX67yzGLKKf3c71Czfd79+Obc89QHmIaqCOmqsalTw22HCYUF/7hpaTla7pmjUqhcM3nm2k82v/ebL3df7NWF3dyVpbfX3Tsy2n4fSVSRs3pWtIYGmDIFDhzI+XS5yjNz7ef0aWNCKQst5erCbu7K0NMDJ51kC7cnlPX84+jkk+F3v3M1/zmUq7olKlU0xVyFpOU6USR2sHj/fpg7N+xWmJBY8PdJWVcmmjMHHnssmNeOoVzjAbMnjmDd9p6C6a9ylcZGSfpkN2rndn4+roOPTtxmA6YJZMHfB2WfKbDAoG8S5boKmTqmsWCPPpYzl3pw5GQ3MEAT8IKOZHVCZ7VMOgv+Pij7TIGtrVBd7aZ5qLG3MJ9i0lJxnbl0uNInuxF9h3hjxBh66kZAQme19EPU1yIeSsVHjnK8OWWfKbCuDhYsYNPDT/HILklerjqD13x9XGcuHa70sY4+tI9fzznnyGBvEme19CruawNUdO4g30yafq/cE8ZMgY+MnsWalzp9X2Wr2Jkzo8CPlcaiNnNp0NIntSpgQ/PMI48ncVZLr+K+FnFFB/9yvTlhzBT45Y0D9A+C6NF8dTpXPVxxW7ZxqHx9saI0c2k5LDpxAo2DAxysqWf7aDdtelJntfQq7msDVHTap1xvThjzi2w8KGwZP5WJ+3ext2HUkce95KrjNvjpV74+KiWr5TC3ZTT1XTv4Yc1paFX1kYkF45CmiJpSZxSNmooO/uV8c4aauiGIcYepYxt5dupcLv39ymOCv5dcdRQGP0vJ4SctX++X2WPr+cINV/GF004LuymxFuZaxH6o6LRPFBZuCGrc4cYLT2HblFluvczUHC1ec9X5guZQwdTPMYJS005Jy9f7YmDAlQmfdFLYLYm9JQtaue2S+bSObURw07Pfdsn82FxFVXTPPwrTvQZVBrpkQSvoH9L9wq8YvXsXOnac52qfUucF8vsGqeEs2JL+f14rnvy8yzfSdwzv3AlnneWmCYmZKJZVxnlG0YoO/hD+mxPkuMOSt02DL38K/vVfYdYsz69XajD1e4xgOGknP/L1fp7EIn3HsKqbDvz888NtxzDEvawyiio++JdTrp5J4OMOZ5wBjY3uj7q+3vPLlRJM/R4jCCuH7+dJLEqD5tlXIO+cXMdJs1rhhBPK2g4/lP1GygSo6Jx/OeXL7b9rbnOw4w719fCe97jL+TIbzhjBUMLK4ft5EovCoDnkHj95/vdbePSkP4jlLJ5xL6uMIgv+PsnXM3lkfVfwg0KLFrmBvCJmtfST38E6rJp7P09ifp8Qhyv7CqRmoJ9DUs0tHd6vDsMQxo2Ulc7SPj4ZqmcS+LjD5Mkwf75b2L25Obj9ZPFzwDXzNcudHvFzAZxyLaZTSPaVxsT9u3l89gK27C9vB6EYxQzkxr2sMoos+Psk9Bs+LrwQ7rgDJk4s62V9Jdwg5edJLIgT4nBkjp+IDlKtgzwz7dTI9ZQLkcSwAAAN80lEQVSLHciNQuVepfEl+IvIYuAbQDXwHVW9Pev5euD7wEKgG7hMVbf4se+oCL1nMm8ezJ4N3d0wwercSy239PMkFoUTYuYVSPP+3ayeego94ydxW8R6yqUM5IZduVdpPOf8RaQauAt4HzAPuEJE5mVtdi2wS1VPAr4OfMXrfqMm9Bs+qqrgsstg796jC3MnVNzmKApCevxkTF0VtQMDvND2rkjegGQDueHxo+d/JrBRVTcBiMi9wMXAuoxtLgZuSX1/H/AvIiKqlRWlQu+ZzJ0Lp53mcv+TJ4fXjpBFqdwyTHNbRjO3fw/86TV85KpLw25OTqGnSxPMj2qfVuC1jJ87Uo/l3EZV+4E9wHG5CRG5TkTaRaS9q6vLh6YljAh85CNucfcyV/5ESVTKLUPX2+s+Ex/4QNgtySsKU7AklR/BP9foYnaPvphtUNWlqtqmqm3NZaxaqSizZsHZZ8Prr4fdktBEpdwydJ2dLvCPGxd2S/IKPV2aYH6kfTqA6Rk/TwO259mmQ0RqgDHAmz7s2+Ty4Q/DM88kdpnHqJRbhurQoaM3AEZc6OnShPKj578KmCMis0WkDrgcWJa1zTLgmtT3lwIPV1q+P1KmTHGln9v8XbEsLpK2QMtxVF2v/7LLYNSowtubRPLcLVTVfhG5HliBK/W8R1XXisitQLuqLgPuBn4gIhtxPf7Lve7XFHDxxdDeDrt2RfqyPyhRKLcMzeuvu8H/d7wj7JaYCPMlJ6Cqy4HlWY99MeP7Q8BH/NiXKVJjI1x3HXzpSzB6NFRXF/4/Jv4OHXJTffzpn9p7boZkc/tUspNPhsWLoaMj7JaYckinez72sUSX+priWPCvdEuWuDt+d+0KuyUmaJbuMSWw4F/p0umf3btd9Y+pTJbuMSWy4J8EJ5/s6r1fey3xUz9UpIEBV9l11VWW7jFFs+CfFJdcAqefbvn/SqMKr77q6vnf+c6wW2NixIJ/UtTUwF/8hZvy2abOqBzbt8Nb3gJXXBHLFbpMeCz4J8moUfBXfwV9fbBvX9itMV51d0NTE3zyk1CbsKkrjGcW/JNm6lT4y790a/729obdGjNcBw7AwYPwmc/AmDFht8bEkAX/JJo/H6680uX/BwYKb2+ipbfXlXV+8pMwfXrh7Y3JwYJ/Ul14Ibz3vbB1a6Knf46dvj5XtXXllbBwYditMTFmwT+pRNwg4bveBVu22AkgDvr7XWXPRz/qTt7GeGDBP8mqq+Hqq+EP/9CuAKKuv9+9Rx/+MHzwg1bZYzyz4J901dXwJ3/iTgB2BRBNfX0u8F9yiQv+FviND5K30oc5Xk2NOwHU1MDKlW41MJsiIBp6e12O/7LL3F3aFviNTyz4GyedAmpshF/+ElpboaEh7FYlW0+PuyHvmmvg3e+2wG98ZcHfHFVV5QYTW1vh7rtd/fjYsWG3Kpm6ulye/7OfhdNOC7s1pgJZzt8cSwTOPRduvtkFnwQvBB8KVXf/RVMT3HKLBX4TGAv+JreTToJ/+Ac3S6RVApVHfz9s3gynngpf+AK0tITdIlPBLPib/CZMgM99Ds4+2wWlQ4fCblHl2rfP1fB/8IPw6U/DyJFht8hUOMv5m6E1NMCf/7lbE+A//sMNDE+ZYoOPfhkcdDNz1tfDX/81LFhgv1tTFhb8TWFVVe5O4FNPhX//d1i71k0QZ9VA3uzb5ybYO+ccN12DTdBmyshT2kdExovIr0Xk5dS/4/Js9ysR2S0iv/SyPxOySZPgb//W3RPwxhtusXBbGax0g4NuUPfgQTfF9ic+YYHflJ3XnP9NwEpVnQOsTP2cyx3AH3vcl4mC6mp3FfDlL8Ps2W4sYP/+sFsVH3v2uDupFy50v8OFCy3NY0LhNfhfDHwv9f33gCW5NlLVlUCPx32ZKJk0CW680S0Of/iwC2g2IJzfvn3uRFlXBzfcYL19EzqvOf/JqtoJoKqdIjLJy4uJyHXAdQAzZszw2DQTuOpqd09AWxs8+ij87GduOoKWFltZKu3QIXevxJgxbhnNM89002gYE7KCn0IReQiYkuOpm/1ujKouBZYCtLW1WTI5Lhoa3BTDb387/PrXsHy5GwtoaUluoOvtdWMiDQ3wsY/BO97hKnqMiYiCf5mqekG+50Rkh4i0pHr9LcBOX1tn4qWpyc08ed55bn6gRx91g5sTJyajbl0V9u6FXbtceueii+A973FrJxsTMV67ZcuAa4DbU//+3HOLTPyNH+8miVuyBJ55xl0JbN3qJo2bONGVjlaSgQFXsnn4sCuBvfRSeNvb3PEaE1GiHkr1RGQC8GNgBvAq8BFVfVNE2oBPqOqfpbb7LTAXGAV0A9eq6oqhXrutrU3b29uH3TYTIQMDsH49PPggPPecC/4TJsCIEfGtdFF1s27u2uWO4Zxz4Pzz4YQT4ntMpiKIyLOq2lZwOy/BP0gW/CtUVxc89RT89rfuXgGA0aPdgGjUrwgGBmD3ble5AzBtmlsE56yzrHLHREaxwT+ho3EmNM3NLhf+oQ/Bjh2wbp07GWzc6J5vbHSBtK4u/B60qqvW2b3bDeBWVbm7nM8+G045xV29GBNTFvxNOETcHEFTprh0SU8PvPQSrFoFL77oyiOrqtyA8ciRbtC0vj64E0I60Pf0uDtvRdxj48e7YL9wIcyZY3l8UzEs+JtoaGpyAXbhwqNVM52dbgnD9evdlcGOHUdTQ4ODLkDX1bmv+np3b4HI0S9wr6Xqtu/rc4Oyvb3uC459vfHj4YwzXK++tdUN3iahSskkkgV/Ez0iLvUzZgzMnevKJdMnhHTOPT3Y2tUF3d3uq6fHzYk/OOjy8yLuRrT015gxMH26qzhqbnarlI0a5U4848dboDeJYsHfxEPmCcEY41nEyyuMMcYEwYK/McYkkAV/Y4xJIAv+xhiTQBb8jTEmgSz4G2NMAlnwN8aYBLLgb4wxCWTB3xhjEsiCvzHGJJAFf2OMSSAL/sYYk0AW/I0xJoEs+BtjTAJZ8DfGmASy4G+MMQnkKfiLyHgR+bWIvJz6d1yObd4qIk+JyFoReV5ELvOyT2OMMd557fnfBKxU1TnAytTP2Q4AV6vqqcBi4E4RGetxv8YYYzzwGvwvBr6X+v57wJLsDVT1JVV9OfX9dmAn0Oxxv8YYYzzwGvwnq2onQOrfSUNtLCJnAnXAKx73a4wxxoOCC7iLyEPAlBxP3VzKjkSkBfgBcI2qDubZ5jrgOoAZM2aU8vLGGGNKUDD4q+oF+Z4TkR0i0qKqnangvjPPdqOBB4DPq+rTQ+xrKbAUoK2tTQu1zRhjzPB4TfssA65JfX8N8PPsDUSkDvgZ8H1V/YnH/RljjPGBqA6/gy0iE4AfAzOAV4GPqOqbItIGfEJV/0xErgL+HVib8V8/rqprCrx2F7B12I07aiLwhg+vEydJO+akHS8k75iTdrww/GOeqaoFi2o8Bf84EJF2VW0Lux3llLRjTtrxQvKOOWnHC8Efs93ha4wxCWTB3xhjEigJwX9p2A0IQdKOOWnHC8k75qQdLwR8zBWf8zfGGHO8JPT8jTHGZKmY4C8ii0Vkg4hsFJHjJpgTkXoR+VHq+d+JyKzyt9I/RRzvDSKyLjWT6koRmRlGO/1U6JgztrtURDRVchxbxRyviHw09T6vFZH/LHcb/VbE53qGiDwiIqtTn+33h9FOv4jIPSKyU0R+n+d5EZH/m/p9PC8ib/Nt56oa+y+gGjdf0Am4uYOeA+ZlbfO/gG+mvr8c+FHY7Q74eN8FjEh9/8k4H2+xx5zargl4DHgaaAu73QG/x3OA1cC41M+Twm53GY55KfDJ1PfzgC1ht9vjMb8DeBvw+zzPvx/4b0CAs4Hf+bXvSun5nwlsVNVNqtoL3IubcTRT5gyk9wHvFhEpYxv9VPB4VfURVT2Q+vFpYFqZ2+i3Yt5jgH8EvgocKmfjAlDM8f45cJeq7gJQ1ZzTq8RIMceswOjU92OA7WVsn+9U9THgzSE2uRg3O4KqmxpnbGoqHc8qJfi3Aq9l/NyReiznNqraD+wBJpSldf4r5ngzXYvrPcRZwWMWkQXAdFX9ZTkbFpBi3uOTgZNF5AkReVpEFpetdcEo5phvAa4SkQ5gOfCX5WlaaEr9Wy9awYndYiJXDz67jKmYbeKi6GNJTa/RBrwz0BYFb8hjFpEq4OvAx8vVoIAV8x7X4FI/5+Gu7H4rIqep6u6A2xaUYo75CuC7qvp/ROQc4AepY845U3AFCCxuVUrPvwOYnvHzNI6/HDyyjYjU4C4Zh7rcirJijhcRuQA39fZFqnq4TG0LSqFjbgJOA34jIltw+dFlMR70LfYz/XNV7VPVzcAG3Mkgroo55mtx84mhqk8BDbg5cCpVUX/rw1EpwX8VMEdEZqdmEb0cN+NopswZSC8FHtbUiEoMFTzeVArkW7jAH/dcMBQ4ZlXdo6oTVXWWqs7CjXNcpKrt4TTXs2I+0/fjBvYRkYm4NNCmsrbSX8Uc86vAuwFE5C244N9V1laW1zLg6lTVz9nAHk0toOVVRaR9VLVfRK4HVuAqBu5R1bUicivQrqrLgLtxl4gbcT3+y8NrsTdFHu8dwCjgJ6lx7VdV9aLQGu1RkcdcMYo83hXAe0VkHTAA3Kiq3eG12psij/lvgG+LyGdw6Y+Px7gTh4j8P1zabmJqHOPvgVoAVf0mblzj/cBG3Hrof+LbvmP8ezPGGDNMlZL2McYYUwIL/sYYk0AW/I0xJoEs+BtjTAJZ8DfGmASy4G+MMQlkwd8YYxLIgr8xxiTQ/wfg3sG9wWwAJAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe4cb179c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 2s, sys: 19.3 s, total: 1min 22s\n",
      "Wall time: 5.55 s\n"
     ]
    }
   ],
   "source": [
    "def maxCoverage(P, rMax):\n",
    "    M, R, X, S = basicModel(1, P, bigM(P))\n",
    "    M.constraint(R, Domain.lessThan(rMax))\n",
    "    M.objective(ObjectiveSense.Maximize, Expr.sum(S))\n",
    "    M.writeTask('p2.task.gz')\n",        
    "    M.solve()\n",
    "    return R.level(), X.level()\n",
    "\n",
    "%time P=data(100); r, x = maxCoverage(P, 0.25); display(P, r, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 3. Smallest diameter k-circle cover\n",
    "\n",
    "**Objective: find $k$ circles which together cover all the points and have the smallest sum of diameters.**\n",
    "\n",
    "This is a different geometric facility location problem, where we want to serve the whole population with $k$ transmitters while minimizing their total range. See the paper [Minimum-Cost Coverage of Point Sets by Disks](http://jeffe.cs.illinois.edu/pubs/pdf/carrots.pdf) for more details. It can also be seen as a form of clustering witk $k$ clusters.\n",
    "\n",
    "The objective here is of course to minimize $\\sum r_j$. A point $p_i$ is covered if $\\max_j s_{ij}=1$, a condition we can linearize as: \n",
    "\n",
    "$$\\sum_j s_{ij}\\geq 1,\\ \\forall i=1,\\ldots,n.$$\n",
    "\n",
    "For $k=1$ we recover as a special case the minimal enclosing ball problem. It is quite frequent that certain locations are served by their own transmitter with range $r_j\\approx 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt0VfWZ//H3k3sgIeESSLgJIoJRrGiwtv6WWnWUWkXbaqu1s+oqHcfOME6ntTPOWG2ly2mnOnVNR1en9jLjdE1rqZ2xtIOlFbVjsVpQVERFwaogEAIkECCQ2/f3x/ccOYaEnCT7nL3P3p/XWlk552Sb/Xg4efZ3P9+bOecQEZFkKQo7ABERyT8lfxGRBFLyFxFJICV/EZEEUvIXEUkgJX8RkQRS8hcRSSAlfxGRBFLyFxFJoJKwAxjIhAkT3IwZM8IOIzydndDSAl1dUFoKZmFHJFHU3Q29vTB2LIwZo8+J8Mwzz+xyztUNdlxkk/+MGTNYu3Zt2GHkn3Pw29/C/fdDVRWMHx92RBJ1XV2wdSvMnQt/8RdQUxN2RBIiM3szm+MCKfuY2UIz22hmm8zs5mMcd6WZOTNrCuK8sdPVBf/1X/CDH0B9vRK/ZKe0FGbMgD/+EW6/Hd56K+yIpACMOPmbWTFwL/BBoBG4xswa+zmuGrgReHqk54ylffvgm9+EX/8ajjsOysvDjkgKiRlMnuzLQEuXQhLvmmVIgmj5nwlscs697pzrBB4ALu/nuK8C3wAOBXDOeNm61f/BbtoEM2dCcXHYEUmhGjfOf/3Lv8DPf+77A0T6EUTynwJsyXi+NfXaO8xsPjDNOffLAM4XL88952/VDx2CKVPUYScjN2qUv3v82c/g29+Gjo6wI5IICqLDt79s9c4mAWZWBNwNXDfoLzK7HrgeYPr06QGEFmG9vbBiBSxbBpMmwejRYUckcVJS4u8in3kGduyAG2+EukEHgEiCBNHy3wpMy3g+FdiW8bwaOAV43MzeAM4ClvfX6eucu8851+Sca6qL8wf10CH4znfgJz+BadOU+CU3zGD6dNi1C77yFXj11bAjkggJIvmvAWab2UwzKwOuBpanf+ic2+ucm+Ccm+GcmwE8BSxyziWzR2rvXvj612HNGjj+eD9SQySXJk3yAwj+8R/hySfDjkYiYsRlH+dct5ktAVYCxcAPnHMbzGwpsNY5t/zYvyFBWlvhrrv85K24l7UkWsaMgbIy+Ld/8xMIzzsv7IgkZIFM8nLOrQBW9HnttgGOPS+Icxac3bvhG9/wLf/Jk8OORpKoosIPKvj+9/2Q0Asu0ACDBNPaPvnQ0gJf+5ofy19fH3Y0kmTl5b6f6T//Ex5+2M8ol0RS8s+1nTt94j9wwNdeRcJWVuYvAD/+Mfzyl7oAJJSSfy7t2gX/9E9+dI8Sv0RJaanvd1q2DH71K10AEkjJP1f27PGJ/+BBmDgx7GhEjpa+APzoR/Doo2FHI3mm5J8LbW1w552+xq8Wv0RZaakvAf3Hf8Djj4cdjeSRkn/Q9u+Hf/5nP7qnoSHsaEQGV1YGU6f61WQ1DyAxlPyD1N3tx1Fv26bhnFJYysv9Z/a734XXXgs7GskDJf+gOAcPPgjr1/tWlEihqaiA2lq/IuiuXWFHIzmm5B+U1av9Qm3Tp2vijBSumhp/B/uv/+pHqUlsKfkHYfNmP2ty8mStxS+Fr77e7wb27/+u/QBiTMl/pHbvhrvv9i2mioqwoxEJxvTpvvN3xYrBj5WCpOQ/EocOwbe+5W+Ta2vDjkYkOOnloJct8xsOSewo+Q9Xby/cf7+/PdZ6PRJHpaV+nsq99/qtRiVWlPyH6+GH4Xe/8xNkROJq9Gg/DPTuu/2kRYkNJf/heOEFvwvX9OlQpLdQYm7CBL8U+be/7UucEgvKXEO1d6/fgnHiRO3CJckxZQps2ACPPBJ2JBIQJf+hcM4vgnX4MFRVhR2NSP6Y+QvAsmV+BrsUPCX/oVi3Dn7/ey3dIMlUXu6/vvc9lX9iQMk/W3v3+olcEyeqzi/JNXGin9So8k/BUxbLhso9Ip7KP7Gh5J+NZ59VuUckraxM5Z8YUPIfzN69fp1zlXtEjlD5p+Apmx2Lyj0i/VP5p+Ap+R+Lyj0iAysr84sZqvxTkJT8B3LggF/SVuUekYHV1fnyz29/G3YkMkTKagNZtcpfAFTuERmYmV/Y8Gc/838vUjCU/PvT1gbLl2sDdpFsVFZCR4dvMEnBUPLvz4oVvrO3rCzsSEQKQ0ODbzC1tYUdiWRJyb+v5mY/fE2tfpHslZX5BpN2/ioYSv59PfSQ34e3pCTsSEQKS0ODbzg1N4cdiWQhkORvZgvNbKOZbTKzm/v5+Q1mtt7MnjOz35lZYxDnDdxbb/l9S9XqFxm6khLfcPr5z8OORLIw4uRvZsXAvcAHgUbgmn6S+4+cc/Occ6cB3wC+OdLzBs45ePBB33mloZ0iw9PQAKtX+4aURFoQWe5MYJNz7nXnXCfwAHB55gHOucz930YDLoDzBuvVV+H55/2epSIyPEVFvgH14IO+QSWRFUTynwJsyXi+NfXau5jZX5rZZnzL/8YAzhuc3l6/LWN1tR+3LCLDN2mSb0i9+mrYkcgxBJH8+8uWR13ynXP3OudmAX8HfKnfX2R2vZmtNbO1LS0tAYSWpRdf9LMUx4/P3zlF4srMN6SWLVPrP8KCSP5bgWkZz6cCx1rp6QHgiv5+4Jy7zznX5JxrqqurCyC0LD38MIwZo1a/SFDGj/cNqjffDDsSGUAQyX8NMNvMZppZGXA1sDzzADObnfH0Q8BrAZw3GNu2wSuvwLhxYUciEh9mUFoKjz0WdiQygBEnf+dcN7AEWAm8DCxzzm0ws6Vmtih12BIz22BmzwGfBz410vMG5okn/PA0tfpFgjVpkh/5094ediTSj0BmMjnnVgAr+rx2W8bjvw7iPIHr6IBHH/Urd4pIsEpKoKcH1qyB888POxrpI9kD2tetg85OreEjkivjx/slH3p7w45E+khu8k+vQ1JbG3YkIvFVVQW7dsHGjWFHIn0kN/m//jq8/bYf5SMiuVNZCb/5TdhRSB/JTf6PPurLPeroFcmtujpfYt21K+xIJEMyk//evfDUU+roFcmHoiLfyHryybAjkQzJTP5PP+1r/sXFYUcikgwTJ8LKlX6AhURC8pK/c/DrX2spB5F8qqiAgwf9hEqJhOQl/x07YPduGD067EhEkqWsDJ55JuwoJCV5yf+ll8KOQCSZxo+HP/wBurvDjkRIYvJ/8kmoqQk7CpHkKSuDw4e12FtEJCv5793rx/cr+YuEwwxeeCHsKISA1vYpGOnNJUYwtv+V7ftYvXk37Ye6qK4o5exZ45nboIliIlkZN87ffV9xhebYhCwRyf+hdW9z58qNXLBqGe9p38PJZWOHlbBf2b6PR17eSXdqnZL2Q1088vJOgIK/AOiiJnkxapTf37e5Gerrw44m0WJf9nlo3dv8/X+vp2XXPk7a+Ue2FVXwyMs7eWX7vsH/4z5Wb979TuJP6+7tZfXm3UGFG4r0Ra39UBdw5KI2nPdI5JjSrX0NvAhd7JP/nSs30tHVw/S27RS7XnqKioedsNPJMdvXC0VcL2oSUWPGwO9/H3YUiRf75L+trQOAk5tfp8eO/O8OJ2FXV5QO6fVCEdeLmkRUTQ1s2qRNXkIW++Q/ubYSc73M3/YKrZXV77w+nIR99qzxlBS9+y0rKSri7FmFPVs4rhc1iaj031B6AIaEIvbJ/4sXz2H6ob2M6jpEZ4nftGW4CXtuwxguPGniO0mxuqKUC0+aWPAdo3G9qEmEVVTAs8+GHUWixX60zxXzp1C7YSz7V/tF3EY6kmVuw5iCT/Z9pf9/NNpH8mbMGG3wErLYJ3+A80raYf40mDIlp+cp5OGScbyoSYRVVPghnwcOaJ2tkMS+7AP4lQSrqwc/biSn0HBJkeyZ+a/t28OOJLHin/y7umDrVr+XaA5puKTIEDnn/zYlFPFP/s3N0Nt7ZIRBjmi4pMgQVVZqxE+I4p/8t2/3LYwc03BJkSGqrlbyD1H8k//mzVCS+35tDZcUGaKKCr+x0sGDYUeSSPFP/hs35ryzF+I7B0AkZ9Kdvtu2hR1JIsV7qGdXlx9OluMhnmkaLikyRL29vtP3hBPCjiRx4p38m5t9vT/Hnb0Sb4U8fyPyRo3ydf/zzgs7ksSJd/LfsSMvnb1RpIQVjDjv4RAJVVXw2mthR5FIgTSJzWyhmW00s01mdnM/P/+8mb1kZi+Y2SozOy6I8w5qdzLH2GvCWXA0fyPH0p2+fd5jyb0RJ38zKwbuBT4INALXmFljn8PWAU3OuVOBB4FvjPS8WWlp8ZtGJ4wSVnA0fyPH0iVZjfjJuyBa/mcCm5xzrzvnOoEHgMszD3DOPeacS//rPgVMDeC8g9u1C8rL83KqKFHCCo7mb+SBmdb2D0EQyX8KsCXj+dbUawNZDDwcwHkHt3t3Ilv+SljB0fyNPFHyz7sgkr/181q/vaxm9kmgCbhzgJ9fb2ZrzWxtS0vLyCNrbU1k8lfCCo7mb+SBc7B/f9hRJE4Qo322AtMynk8Fjpq1YWYXArcA5zrnDvf3i5xz9wH3ATQ1NY1smE53t18udnzyEp7W5w+W5m/kWG8v7NNghHwLIvmvAWab2UzgbeBq4BOZB5jZfOA7wELn3M4Azjm4/fuPzCBMICUsKRilpb5/TvJqxGUf51w3sARYCbwMLHPObTCzpWa2KHXYnUAV8FMze87Mlo/0vINqb09s4hcpKGVlSv4hCGSSl3NuBbCiz2u3ZTy+MIjzDIk6kEQKQ3m5kn8I4rvuwf79mjgiUgjKyvzgDMmr+C7v0NYWdgQiiZb1EiNlZX5CpnMq1eZRfFv+bW2+I0lE8m5IS4wUFfnEf+hQnqNMtvi2/Lu71YqQ2Ci0hfqOtcRIv3E7pzJtnin5i0RcIa4sOuQlRsyU/PNMyV8Kx/r1sGoV7N0LNTVwwQUwb17YUeXckFvREVBdUdpvoh9wiREl/7yLb82/p0fJP07Wr4df/MInfvDff/EL/3rMFeJCfcNaYqSnJ8dRSSYlfykMq1b5bTkzdXX512OuEBfqG9aaSGr551V8yz4SL+kWf7avx8jZs8a/q+YPhbFQ35CXGEnornthiW/yLy4OOwIJUk1N/4m+pib/seRZYhbq099sXsU7+aslER8XXOBr/Jmln9JS/3oCxH6hPueO7OoleaHkL4UhPaongaN9EkPJP6/im/xLSpT842bePCX7OFPyz6v4vtvFxRo9IFIoVPbJu/i+21VVGjcsEnGvbN/H9594nZ+u3cI5d6/moXVvhx1SYsS37DNunFr+BajQ1rCR4UsvW0FXJ8VllbzV3snf/7eftHfF/CkhRxd/8W75a5JXQRnSSpBS8NLLVpT1dNFWUQVAR1cPd67cGHJkyRDfln91tZJ/jgXdSi/ENWxk+NIX+bKebraNOfLvu62tI6yQEiW+Lf/q6rAjiLVctNILcQ0bGb700g+lPV20Vh5J/pNrK8MKKVHi2/KvqnpXzV+15GDlopU+5JUgpaCll60o7elh9yj/maksLeaLF88JObJkiG/Lv7zcbw/X3a1acg7kopU+rJUgpWClF38bVVbM/vLRTKmt5GsfmafO3jyJb8vfjFcPF/OHJ15jV8/Ra4aoljwyuWilF9oaNrqbHLm5DWOYe9oULv3CQmhsDDucRIlt8n9o3du8tuUw4w90QGokQV+qJQ9frlaaLJQ1bApxd61IUx9d3sW27HPnyo3sLKuitKd7wGNUSx6+Ya3XHiPH6vOQIXKuoJL/mjvuYcfYSfRaETvGTmLNHfeEHdKwxLblv62tg9ZRYygbIPmrljxyhdJKzwWNTAqIc/6rqv+786hZc8c9nHL7TVR2HQagvm0nNbffxBpgwS1Lwg1uiGKb/CfXVrKnsoYid/QsX9Vnj0217MFpZFJAurp84i8pjFQ07a6vvpP40yq7DjPtrq9CgSX/2JZ9vnjxHA7UjKM3Y6JXSVERC0+uZ/H/m6lkNgCNjMqORiYFZP9+mDkz7CiyNrGtZUivR1lsk/8V86ew5NpzqCotAucSV5MeLtWys5P0Po/AtLfDnMIZ17+ztm5Ir0dZYdxrDdNl7z8BPjDPLxU7alTY4RQE1bKzl+Q+j8CYwXHHhR1F1rbcdCs1GTV/gI7ScrbcdCv1IcY1HLFt+b9jzhx/aylZGahmrVq25MzkyWFHkLUFtyzhxS/fxY7aifRi7KidyItfvqvgOnshoORvZgvNbKOZbTKzm/v5+Tlm9qyZdZvZlUGcM2tz5kCHForKlmrZkjddXVBZCWPHhh3JkCy4ZQn1rc0UuV7qW5sLMvFDAMnfzIqBe4EPAo3ANWbWd6reW8B1wI9Ger4hmzxZq3sOgWrZkjft7TBrlv4+QxJEzf9MYJNz7nUAM3sAuBx4KX2Ac+6N1M/yv7vK5MlHxhLrQ5YV1bIlL/bvL6jO3rgJouwzBdiS8Xxr6rVoGDUKxo9X6Uckigqoszdugkj+/TWn3bB+kdn1ZrbWzNa2tAQ4bladviLRVECdvXETRPLfCkzLeD4V2DacX+Scu8851+Sca6qrC3DcrDp9RaKlQDt74ySI5L8GmG1mM82sDLgaWB7A7w2OOn1FokWdvaEbcfJ3znUDS4CVwMvAMufcBjNbamaLAMxsgZltBa4CvmNmG0Z63iGZNg2Ki6F74BU+RSSP2tvhPe8JO4pEC2SGr3NuBbCiz2u3ZTxegy8HhaO8HObNg40bIchykogMnUt1CWrzllDFf4Zv2llnwcGDYUchIh0dfgRefaEtiBAvsV7b511OPNF/D2i8v5Y9FhmmPXvgkktU7w9Zclr+tbV+6dh9I1+aWMsei4xAby+cemrYUSRecpI/wPvfD3v3jvjXaNljkWHq7ISyMpgxI+xIEi9Zyb+x8Uhn0who2WORYdqzBxYsKJidu+IsWcm/ocFPKhlhx6+WPRYZpsOHoakp7CiEpCV/M3jf+2D3yMozWvZYZBh6e/3f4OzZYUciJGm0T8pvSus5/NzbbNrYMexROunjNdpHZAhaW+Gkk7SrXkQkKvk/tO5tbv3DXv6uxyjp6ab9EDzy8k6AYV0AlOxFhmD/fn/nLZGQqLLPnSs30t5rrJ12MuMP+lE/GqUjkgc9PX4v7VNOCTsSSUlU8t/W5lf2fHraKZT09rwz8kejdERybOdOeO97/XwbiYREJf/JtZUAbK+ewJu19dQc8mv8a5SOSA4550f5nH9+2JFIhkQl/y9ePIfK0mIw4/Hjm6jq7NAoHZFca2/3y6rPmhV2JJIhUR2+V8z3u0veuXIjr/bMoGj0KC6aVcOJ6rgVyZ3WVli8WGv5REyikj/4C0D6IsBDwPJo7TsjEitdXVBaCqefHnYk0keiyj5HOftsX4/ss06PiASkuRk+8AGN7Y+gZCf/ujq/m9CuXWFHIhI/zvnd8845J+xIpB/JTv4AF12kTV5EcmHPHr+PxpQpYUci/VDynzMHxo2DAwfCjkQkXtrb4YMfDDsKGYCSf3ExfOhD0NISdiQi8XHwIFRXa0ZvhCn5g9/kpapK5R+RoDQ3w5VX+pE+EklK/gCVlXDVVf4DKyIjs28fTJjgG1USWUr+ae9/v//ABrDHr0hiOef3y7jmGrX6s7TmjnvYMXYSvVbEjrGTWHPHPXk5r5J/Wmmp/8Du3h3IVo8iidTaCtOna1JXltbccQ+n3H4T9W07KcJR37aTU26/KS8XgMTN8D2m00/3H9zWVj8CSLLyyvZ92thGfKOprQ1uuMEv3yyDmnbXV6nsOvyu1yq7DjPtrq/CLUtyem79C2UqKoKrr/YfYLX+s/LK9n088vLOd5bFbj/UxSMv7+SV7SqfJU5LCzQ2+t26JCsT2/ofZTjQ60FS8u+rsdF/aehnVlZv3k13n+UxtEFOAvX2+rkyH/uYFnAbgp21dUN6PUhK/n2Z+Q/wgQNa8ycLA22Eow1yEmbHDjjzTDj++LAjKShbbrqVjtLyd73WUVrOlptuzfm5lfz7c/zxsGCB/0DLMQ20EY42yEmQnh7o7IQPfzjsSArOgluW8OKX72JH7UR6MXbUTuTFL9/FghzX+0HJf2Af+Yhfjra7O+xIIu3sWeMp6dO5pw1yEmbbNjjvPL9hiwzZgluWUN/aTJHrpb61OS+JHwJK/ma20Mw2mtkmM7u5n5+Xm9lPUj9/2sxmBHHenJo8GS69FLZuDTuSSJvbMIYLT5r4Tku/uqKUC0+aqNE+SXHgAJSXq9VfgEY81NPMioF7gT8BtgJrzGy5c+6ljMMWA63OuRPM7Grgn4CPj/TcOXfppbB2rV+dUEM/BzS3YYySfRL19vpZ8TfeCDU1YUcjQxREy/9MYJNz7nXnXCfwAHB5n2MuB+5PPX4QuMCsAIYElJfDn/2Zn/Wr8o/Iu23bBmedBWecEXYkMgxBJP8pwJaM51tTr/V7jHOuG9gLHFUUNrPrzWytma1ticpQy+OPh8suU/lHJNOBA1BWBtdeq6GdBSqI5N/fv3zfGVLZHINz7j7nXJNzrqmuLvfjXLN26aXQ0ODLPyJJly73fPrTKvcUsCCS/1ZgWsbzqcC2gY4xsxKgBiicTKryj8gRKvfEQhDJfw0w28xmmlkZcDWwvM8xy4FPpR5fCTzqXIGtn6Dyj4jKPTEy4uSfquEvAVYCLwPLnHMbzGypmS1KHfZ9YLyZbQI+Dxw1HLQgqPwjSaZyT6wEsqqnc24FsKLPa7dlPD4EXBXEuUKVLv8sXQpjxkCJFkWVBFG5J1Y0w3eojj8eFi2Ct97Syp+SHG1tfsc7lXtiQ8l/OBYtgtNOU/1fkuHQIT/Y4XOfU7knRpT8h6OkBK6/3m/7GJX5CCK50NPjyz2f/jTMmhV2NBIgJf/hqqryLaHubti/P+xoRILnnC9vXnIJnH122NFIwJT8R6KhAZYs8a3/zs6woxEJ1ttvw6mnwpVXqs4fQ0r+IzVvHnziE77+r81fJC5aWvxihn/+5xrVFlNK/kG46CI491yNAJJ42L/f72Xxuc/58qbEkpJ/EMzgk5/0HWLNzWFHIzJ8nZ2wc6cvZ2pzllhT8g9Kebn/gykvh9bWsKMRGbreXtiyBa65xtf6JdaU/IM0dqy/VT5wANrbw45GJHvOwRtv+PLlwoVhRyN5oOQftJkz4Qtf8K1/DQGVQpBO/O99L1x3nUb2JISSfy6cdBL8zd/A7t1w8GDY0YgMLD2Wf/58v26VRvYkhpJ/rsyb5/c2bW6Gjo6woxE5Wjrxn3wyfPazfqlmSQwl/1yaP993Au/YoQuAREs68c+Zc2SggiSKkn+uLVjg/7iam1UCkmhIJ/7GRvjrv4aKirAjkhAo+efDggV+FFBLizqBJVzOwZtv+rLkX/2VX6ZZEknJP19OOw0+/3m/C5iGgUoY0on/jDP83aha/Imm5J9P8+bB3/6tXxtdW0FKPnV3wx//6Idz3nCDOndFyT/v5s6FL30Jior8OulaC0hyraPD1/gXLfL7UJSWhh2RRICSfxiOOw6+8hX//c03tRqo5E5rq+9r+uxn4aMfheLisCOSiFDyD0ttLXzxi3DeeX52pfYDkCA554cYA9x6K7zvfZq5K++i6XxhKiuDT30Kpk6FH/7QbwtZXR12VFLoent9mWfWLN+xO3Zs2BFJBCn5h80MLrzQL5/7rW/5+uzEiWFHJYWqs9NvLHTOOfCnf6rJWzIglX2iorHR9wPU1PhlddURLEO1f7/fevHaa2HxYiV+OSYl/yipr/cjgU491Q/L6+oKOyIpBM75DVja2+Gmm+Dii1Xfl0Ep+UfN6NG+TnvVVX4oaEtL2BFJlHV2+gEDEyf6O8d588KOSAqEav5RVFwMl10G73kPfPe7/i5g6lSNz5YjnPMNg0OH4GMf8619LccsQ6BPS5RNnw633QYPPwz/8z8wahTU1YUdlYSts9PX9mfOhM98xjcMRIZIyT/qSkv9zMzTTvN3AW+8AVOm6C4gidTalwCNqOZvZuPM7Ddm9lrqe78Dis3sV2bWZma/HMn5Ei19F/CRj6gvIIkya/tLl8KHPqTELyMy0g7fm4FVzrnZwKrU8/7cCfzpCM8l6buA22+HceN8MtCIoHhLj+TZscO39r/0JZV5JBAjTf6XA/enHt8PXNHfQc65VYDWMQ5K5l3Ajh1+Uk9PT9hRSdD27fMX+EmT1NqXwI30kzTJObcdwDm33cw0NTVf0ncB73sfLF8OTzzhJ/XU12uMd6E7eNDv/FZX5zdcOf10vwqsSIAGTf5m9ghQ38+Pbgk6GDO7HrgeYPr06UH/+niqq/OzOS++GH72M3j2Waiq8usE6SJQWA4fhu3b/b/f4sX+wq6OfcmRQZO/c+7CgX5mZs1m1pBq9TcAO0cSjHPuPuA+gKamJq1vMBRTp8KNN8LmzfCTn8Crr/qVQ7WoV/R1dfmkX1oKH/+4X+lV2ytKjo207LMc+BTw9dT3n484Ihk+MzjhBPiHf4AXX4Qf/9hPENNqodHU0+OTvnO+nn/RRfp3krwZafL/OrDMzBYDbwFXAZhZE3CDc+4zqedPAHOBKjPbCix2zq0c4bllIGZ+mn9jIzzzDDzwgL8IjBsHY8aoHBS2ri5f0+/pgXPP9bO5x48POypJGHMRXT2yqanJrV27Nuww4qGz018EVqzwK4aWlfnx4ho5kj/O+YXX9uzx7/8HPuDLOw0NYUcmMWNmzzjnmgY7Tn/9SVBW5jsPzzrLDx187DFYvdpv+jFhgl9MTnKjp8eP0z982Cf6z3zGj94ZNSrsyCThlPyTxMyvBzNzJlx5JaxZA//7v34f4cpKfyHQkMJgHDzoZ2Gb+Yvu+ef7nbVUcpOIUPJPqjFj4IILfOlh40b49a/h+ed98p8wQaNNhqOnx5d1Dh707+/HP+4Tf21t2JGJHEXJP+mKi33HcGOjb6muXg2PP+4fO+cTlzqJB9bZCbt3++9FRXDKKX5bzpNOUp+KRJo+nXJEXR1ccQVcfrlfMnjDBnjySb8ZOPi+gbFj/QUjqZzzLfs9e3yfSUWFb92fcQbMnq07JikYSv5yNDM/aWzqVD9zuLXVl4bF20RdAAAGJUlEQVSefhrWr/dJr6TED09Mwj6xvb2wd69fawf8RfKyy/xw2hkzkn0xlIKl5C+DGzvWt27POsuvJb9pk19G4umnoaPDH1NU5CcoVVUVdjJ0zv8/trf7/7d0ueuEE/xdUWOjHyarMpgUOCV/GZqKCl/XPuUUuPZa3zewbZsfQvryy/57eoXRqF8QMhP9wYM+Xuf8Hc1pp8GJJ/q7n4YGDYeV2FHyl+ErLvariNbX+7Hr4BP/sS4IzvkkW1bmS0ZlZf4rF52jzvmO2M5OP84+/Tg9nFWJXhJMyV+CdawLws6dR2a5trTArl1HHvdNys75WruZ/yoqOvI4/fPMLzhyTObvqKnxZasJE4581dT4UUxK9JJgSv6Se5kXhP6kW+jt7f5r/37/vbMTurv9xSPze1GRv1MoLvZfJSX+a/ToI2Wm6mr/XJPWRPql5C/hM/MloPJy3zIXkZxTs0hEJIGU/EVEEiiySzqbWQvwZgC/agKwK4DfE7SoxgWKbTiiGhdEN7aoxgWFHdtxzrm6wX5JZJN/UMxsbTZrW+dbVOMCxTYcUY0LohtbVOOCZMSmso+ISAIp+YuIJFASkv99YQcwgKjGBYptOKIaF0Q3tqjGBQmILfY1fxEROVoSWv4iItJHLJK/mS00s41mtsnMbu7n5583s5fM7AUzW2Vmx0UothvMbL2ZPWdmvzOzxqjElnHclWbmzCwvox+yeM+uM7OW1Hv2nJl9Jh9xZRNb6piPpT5vG8zsR1GIy8zuzni/XjWztnzElWVs083sMTNbl/obvSRCsR2XyhkvmNnjZjY1T3H9wMx2mtmLA/zczOxbqbhfMLPTh3wS51xBfwHFwGbgeKAMeB5o7HPMB4BRqcefBX4SodjGZDxeBPwqKrGljqsG/g94CmiKQlzAdcA9Ef2szQbWAWNTzydGIa4+x/8V8IMIvWf3AZ9NPW4E3ohQbD8FPpV6fD7wwzzFdg5wOvDiAD+/BHgYMOAs4OmhniMOLf8zgU3Oudedc53AA8DlmQc45x5zzh1MPX0KyMvVO8vY9mU8HQ3kqxNm0NhSvgp8AzgUsbjCkE1sfwbc65xrBXDO7YxIXJmuAX6ch7ggu9gcMCb1uAbYFqHYGoFVqceP9fPznHDO/R+w5xiHXA78p/OeAmrNrGEo54hD8p8CbMl4vjX12kAW46+Y+ZBVbGb2l2a2GZ9kb4xKbGY2H5jmnPtlnmLKKq6Uj6Zudx80s2n5CS2r2E4ETjSz1Wb2lJktjEhcgC9jADOBR/MQF2QX21eAT5rZVmAF/s4kH7KJ7Xngo6nHHwaqzWx8HmIbzFDz3lHikPz720+v39azmX0SaALuzGlEGafs57WjYnPO3eucmwX8HfClnEflHTM2MysC7ga+kKd43jl1P6/1fc9+Acxwzp0KPALcn/OovGxiK8GXfs7Dt7C/Z2a1EYgr7WrgQedcTw7jyZRNbNcA/+Gcm4ovZ/ww9fnLtWxiuwk418zWAecCbwPduQ4sC0P5N+9XHJL/ViCz5TeVfm4bzexC4BZgkXPucJRiy/AAcEVOIzpisNiqgVOAx83sDXxdcXkeOn0Hfc+cc7sz/g2/C5yR45iyji11zM+dc13OuT8CG/EXg7DjSrua/JV8ILvYFgPLAJxzvwcq8OvXhB6bc26bc+4jzrn5+PyBc25vHmIbzFBzy9Hy0XmR446REuB1/K1sutPm5D7HzMd37MyOYGyzMx5fBqyNSmx9jn+c/HT4ZvOeNWQ8/jDwVFTeM2AhcH/q8QT8rfn4sONKHTcHeIPU/J4IvWcPA9elHp+USmI5jzHL2CYARanHdwBL8/jezWDgDt8P8e4O3z8M+ffn638kx2/SJcCrqQR/S+q1pfhWPvjSQDPwXOpreYRi+xdgQyqux46VgPMdW59j85L8s3zPvpZ6z55PvWdzo/Kepf4Yvwm8BKwHro5CXKnnXwG+nq/3agjvWSOwOvXv+RxwUYRiuxJ4LXXM94DyPMX1Y2A70IVv5S8GbgBuyPic3ZuKe/1w/jY1w1dEJIHiUPMXEZEhUvIXEUkgJX8RkQRS8hcRSSAlfxGRBFLyFxFJICV/EZEEUvIXEUmg/w8lSD+xCJqDyQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe4f461f2e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 18s, sys: 30.8 s, total: 1min 49s\n",
      "Wall time: 4.23 s\n"
     ]
    }
   ],
   "source": [
    "def minDiamKCircleCover(P, k):\n",
    "    M, R, X, S = basicModel(k, P, bigM(P))\n",
    "    M.constraint(Expr.sum(S,1), Domain.greaterThan(1.0))       #sum inside each row, i.e. along dimension 1\n",
    "    M.objective(ObjectiveSense.Minimize, Expr.sum(R))\n",
    "    M.writeTask('p3.task.gz')\n",        
    "    M.solve()\n",
    "    return R.level(), X.level()\n",
    "\n",
    "%time P=data(20); r, x = minDiamKCircleCover(P, 3); display(P, r, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 4. Smallest area k-circle cover\n",
    "\n",
    "**Objective: find $k$ circles which together cover all the points and have the smallest sum of areas.**\n",
    "\n",
    "Here we have another version of the previous problem where the cost of a network is measured (roughly, up to intersections of balls) by the total area it covers, so the $\\ell_2$ norm replaces the $\\ell_1$ norm of the radii vector. As pointed out in [here](http://jeffe.cs.illinois.edu/pubs/pdf/carrots.pdf) the most realistic models of wireless transmission actually optimize for the $\\ell_p$ norm with some $p>2$. \n",
    "\n",
    "This problem is similar to the previous one, but with a different objective: we minimize $\\sqrt{\\sum_j r_j^2}$. This is expressed using a conic constraint. Optimizing the $\\ell_2$ norm results in a more balanced collection of balls, as we can see below for the same dataset as in Problem 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt4lOWZ+PHvnRMJJCQQCIQAhiIKKFZs0FpPVbB4WMVabdHVarVaraddW39rt7vbrr2qbvldttvV7pZVW9ftSq36s2pRroK6HqoWEAU5yckqIUCAQAIJOc3z++OekRBynnfmfd9578915WJmMszcmcP9Pu9zuB9xzmGMMSZasvwOwBhjTPpZ8jfGmAiy5G+MMRFkyd8YYyLIkr8xxkSQJX9jjIkgS/7GGBNBlvyNMSaCLPkbY0wE5fgdQHdGjBjhKisr/Q7DGGNCZfny5buccyN7u19gk39lZSXLli3zOwxjjAkVEflLX+7nSbePiJwnIutFZKOI3N3D/S4TESciVV48rzHGmIFJOvmLSDbwEHA+MBW4QkSmdnG/IuB24J1kn9MYY0xyvGj5nwxsdM5tds61AAuAOV3c70fAT4CDHjynMcaYJHiR/CuATzpc3xq/7VMiMh0Y55x7wYPnM8YYkyQvkr90cdunmwSISBbwU+A7vT6QyI0iskxEltXW1noQmjHGmK54kfy3AuM6XB8LbOtwvQg4HnhVRD4CPg8819Wgr3NuvnOuyjlXNXJkrzOVjDHGDJAXyX8pMElEJohIHjAXeC7xS+fcPufcCOdcpXOuEngbuNg5Z/M4jTHGJ0knf+dcG3ArsAhYCzzpnFstIveIyMXJPr4xxhjvebLIyzm3EFjY6bZ/6ua+X/TiOY0xxgyc1fYxxpgIsuRvjDERZMnfGGMiyJK/McZEkCV/Y4yJIEv+xhgTQZb8jTEmgiz5G2NMBFnyN8aYCLLkb4wxEWTJ3xhjIsiSvzHGRJAlf2OMiSBL/sYYE0GW/I0xJoIs+RtjTARZ8jfGmAiy5G+MMRFkyd8YYyLIkr8xxkSQJX9jjIkgS/7GGBNBlvyNMSaCLPkbY0wEWfI3xpgIsuRvjDERZMnfGGMiyJK/McZEkCV/Y4yJIEv+xhgTQZ4kfxE5T0TWi8hGEbm7i9/fJCKrROQ9EXlDRKZ68bzGGGMGJunkLyLZwEPA+cBU4Ioukvv/OOemOedOBH4CPJDs8xpjjBk4L1r+JwMbnXObnXMtwAJgTsc7OOfqO1wdAjgPntcYY8wA5XjwGBXAJx2ubwVO6XwnEbkFuBPIA87x4HmNMcYMkBfJX7q47YiWvXPuIeAhEbkS+AfgmiMeSORG4EaA8ePHexBawLW0QEMD7N+v/zY0wN690NgI7e3Q1gbOQW4uZGdDQQEMGwZFRfpTWKj/5uX5/ZcYY0LGi+S/FRjX4fpYYFsP918A/HtXv3DOzQfmA1RVVWVO15BzsG8fbNsGW7fC2rWwaRPU10NWFkj8+BmL6eXsbP03cbtz+tPerv9mZR26PRbTA8DRR8PkyTBuHIwZA8XFh/6/yXzNzdp4OHBAPyexmN6elQWDBulnZMgQ/WwZgzfJfykwSUQmANXAXODKjncQkUnOuQ3xqxcCG8h09fXw4YewdCmsWaOtexFN2EOG6Jdx2LDkE7RzegaxcSO8997hzzF1Kpx8MhxzDAwd6s3fZfzV0KCNiOpqWLcOPvkE6urg4EFN9B0bBh0lGhBDh0JpqTYWjj5aGwqjRunZpYmUpJO/c65NRG4FFgHZwKPOudUicg+wzDn3HHCriMwCWoE6uujyCT3nYMcOTfRvvaXJGLSrprhYv3CpIKItu0GDDn+O5mb44AM9+ABMnAhf+IIeEEaNsrOCsGhrg7/8BVau1M9Vbe2hA/zgwfozYsShs8WeJBoKe/fCq6/CH/946Axz8mT4/OdhyhR9PJPxxHVuIQREVVWVW7Zsmd9h9O7gQXj3XXjxRe3SAW1dFRcfaoX5LRbTM5F9+/R6RQWcfz6cdJIenEywOAdbtsBrr8E77xxq1Q8frsne6wN3e7seEPbv1+euqICzzoJTTtHPsQkVEVnunKvq9X6W/AeopgZefx2WLNFW9rBhmvSD3qJ2Tg8EdXV6tnD22XDmmXr6b/zV1AQrVsDChdqQyMvTVng6B/Sd04NAXZ1eP/VUOOcc+Mxngv/ZNoAl/9RwDjZsgGef1e6d7GztQgnrbJuWFti5U1t+kyfDJZfo+IB9ydPr4EHthvn97/UAEJSGRHu7fj6am2HCBJg71z4fIWDJ32uffAJPPaWDqoWF2iLLlC+Bc7Brl7b4PvtZuOwyiMJUW7+1tsLbb8OTT+pA7ujRkJ/vd1RHcg727NEzxhNO0M/HUUf5HZXphiV/r+zcqS2yN9/UL2YmD5YmBq2bmuC002DOHP17jfc2boT//E99vUeO1NlZQeecfh8aG7W78PLLwxF3xFjyT1ZzM7zwAvzhD9q9U14enAHcVIvFdEyjrQ0uuAAuuiiYLdIwOngQnntOP1fFxTqIGzaxmE41LSyEG26A447L3AZRCFnyT8amTTB/vrbKxo6FHC+WQ4RQW5sOPI4apV/yo4/2O6Jw27wZ/uM/dLpmRUX4P1f19dpdeM458LWv2cyxgLDkPxDNzfD88/oT1lZZKuzZo9NEL7wQLr7YzgL6yzntNnzkEV3cl0mfq1hMx8MqKuD226GszO+IIs+Sf39t2QK//CVs3x7t1n53EmcBZWVw00069c/0rq0Nfvc7nb5ZUZG5B84dO7Tr54474Nhj/Y4m0vqa/CPSid2DRKvsRz/SgazKSkv8XcnJ0dfm4EF9rV5//cgSAuZwjY3wr/8KL72kr12mJn7QrsH8fLjvPl2cZgIv2lmurQ2efloHdjO5VealxCrT+fP1TODyy+1g2ZX9++FnP9N+/srKaAyIDh2qCwcffli7UGfNisbfHVLR/dbu368J7L339Mtp1Q77Lj9fX7OXXtJZH9/6lvZlG3XgADzwgNbkGTcuWglw0CDtNn38cR0PmD3b74hMN6KZ/Ldvh5/+VGcqTJjQry/nupp63ty0m4aDrRTl53LaxFIml0ewYmZ2th4A1q2De+6BO+/U6bBR19SkLf6PP9bEH0V5eXoA+O//1stnn+13RKYLkUj+z66oZt6i9Wzb28QJsp8Htixi4vD8fn8519XUs3jtTtritdIbDrayeO1OgNAfAAZ0UBPRL3ltLfz4x3D33Xo9qmIx+NWvdAFX1FdIJw4Av/61roafNs3viEwnGT/g++yKar73zCqq9zZRXr+TSxc9zlsbd7Guvf/9+29u2v1p4k9oi8V4c9Nur8L1ReKg1nCwFTh0UFtXU9/L/4wbOVIPBPfeq10dUfXCC1p2efz4aHX1dGfQIP1sPPigLho0gZLxyX/eovU0tbYzun4XN739FDERduUXDihhJ5JjX28PC08OaqWluiHI/ffrvO+oWb5caz9Z4j9cYaFOCPjZz3SczQRGxif/bXubGLl/Dze9o4l/X4EOTA4kYRfld73bUXe3h4VnB7Vhw/QA8C//ortNRcWOHbpy13bE6trIkTq+9qtf2fTgAMn45D8pv51vvfM04mLsLTjUhz2QhH3axFJyOtX3ycnK4rSJKdqlK008PaglVq/Om6cbhGS69nZ49FEdALciZ90bO1Z3lfvzn/2OxMRldvJvaeGnu/7E0LaD1A0+tCPRQBP25PKhzJpS9mlSLMrPZdaUstAP9np+UBsxQqc7PvSQ7hmQyV59FdauteqnvRHRktW//vWhjWKMrzJ3to9z8MQTHNdQQ/aMKTRs3uPJ9MzJ5UNDn+w7S/w9nk5hLS/XjW9+8xu49trM7AffsQOeeEIXCGbi3+e1wYO1TtTjj8Ntt9lr5rPMTf6vvKJbLFZWMjkri8ljUr8XaZjXAHh+UBPRwc9XXtEptbNmeffYQeAcLFigZb4HDfI7mvCoqIBly2D1ajj+eL+jibTM7PZZtw4ee0w/aGmqwZ/0dMlMlJV1aLXn2rV+R+OtTZvg3XdtYVt/iejEgCee0PES45vMS/4NDdrXXFqa1hZZpq4BSFpiE/Jf/ELrv2cC5+C3v9VpjNZ10X/DhmldqOXL/Y4k0jIr+Se+lI2NWmQqjTJ1DYAnior0PVmwIDOm+n3wAaxfrwc1MzAjRujnIdMnBARYZiX/lSu1nGxFRdqfOlPXAHimogLeeAPef9/vSJLjnO7pXFJirf5kFBXp4G/YPw8hljnJv6FBS8mOHOnLXruZugbAM1lZuhHMww/rexVWW7dq7Z5hw/yOJPyKi+HFFzPjbDCEMiP5d+zu8am0cKauAfBUYaFWvXziifB+4f/3f3UVr7X6k1dSovsdRLEcSABkxlTPjRt1Z6mjjvI1jExcA+C5igrdOe2ss8K33d+BA5r8bZ9ab4jogfTVV+HrX/c7msgJf8s/Md+6sNCX7h7TT1lZenb2299qCeQQeHZFNafd/zJzvj2fP7z7Cet2NfkdUuYoK9OxIBv4Tbvwt/xXrdKWf2Wl35EESqAXnJWW6jz5lSvhxBP9jqZHiZLgTa3tzNy2nt0umy0ZsodDIOTm6naqmzfD5Ml+RxMpnjSVReQ8EVkvIhtF5O4ufn+niKwRkZUiskREvOmfaW/XVv+wYdYH20HgF5yJaAG4BQv0ix9giZLgeW2tTNm5hX35hbZ+w2vZ2bqdqkmrpJO/iGQDDwHnA1OBK0Rkaqe7rQCqnHMnAE8BP0n2eQGtElhdrQNH5lOhWHBWXKwbfAS8yuO2vdrFM35vDdkuRnuW7vVs6zc8NHy4boITkm7ATOFFy/9kYKNzbrNzrgVYAMzpeAfn3CvOucb41beB5Pf6a2uDJ5/UqZ3mMKFZcDZypL6HrQGLq4MxJQUATKr9mFiHs0tbv+GhggKd/rt9u9+RRIoXyb8C6DhXa2v8tu5cD7yY9LOuWaOlYQsLk36oTBOaBWeFhVrzf80avyPp1l2zj6UgN5tJuz9mf95gwNZvpIRz0doAKAC8SP5ddbZ3OYlbRK4CqoB53fz+RhFZJiLLamtre37Wl16yzTO6EaoFZ4WFutAnoC6ZXsH9F09hUnMdjbn5tn4jVXJydNDXpI0Xs322AuM6XB8LHHEIF5FZwPeBs5xzzV09kHNuPjAfoKqqqvtVQNu3a2vR53n9QZWS+vypUlqqdXK2bYMxY/yOpktzxuTAtNFaotqkRmGhfg5M2niR/JcCk0RkAlANzAWu7HgHEZkO/BI4zzm3M+lnfP11nSFgM3y6FZoFZyL6Xr72Gsyd63c0Xaup8TuCzFdYCH/5i87gy872O5pISLrbxznXBtwKLALWAk8651aLyD0icnH8bvOAQuB3IvKeiDw34Cc8eBAWL7Zt8zJJWZlu+tIU0MVTe/aEtxxFWGRn62yf/fv9jiQyPFnk5ZxbCCzsdNs/dbjs3TZOa9fqasC8PM8e0vgsL0/f0zVr4HOf8zuaI9XW2uctHUR01k9x6nfdM2Es77B0KeTn+x2F8Vp+vr63QbRrlyX/dAlzxdeQCVd5h9ZW3f3HNtHIPMOHw4oV+h7nBmxK6p494dqnd9Uq3b963z5tRc+cCdOmpT2MfpcYsW6ftApXy3/LFu0eCFpyMMnLzdXEH8Tpfg0NOhUxDFatguef18QP+u/zz+vtaTSgEiPOBXfcJwOF5BMd9/77Vrkzk2VlaY2XoJV6bm/3/XPX51b0kiVHrphubdXb09j676nESI+tf9vUPW3Ck0md0/ofpQFcqGS8UVqq73HQZtb4XHOmX63oRIu/r7enyIBLjFjyT5vwJP89e7QUwODBfkdiUqWgAOrrdYA1SLKyfD0g9atQX3czZdI8g2ZAJUacC0/3WgYIT/K3uh/RIBK8RVU5Ob4m/361omfOPHJMLDdXb0+jAZUYEfG9ey1KwvNKf/yx3xGYdHBOV3oGydChvlYe7Vcreto0uOiiQy394mK9nubZPgPa0zory87s0yg851jr1lkFzygoKtKFfBdd5Hckh5SW6kKvoiJfnv60iaUsXrvzsK6fHlvR06b5MrWzs36XGBHx7TWOonAkf+d0q8bhw/2OxKRaYaFO94zFgtMFMGKEr3vMhqpQX7Is+adNOJL/nj3Q3Gzz+6MgN1cT7e7dwdmox+fkDyEq1JcM5+zsPo0C0rTqxe7dVsEzSkT0PQ+K0tLgnIVkqrY2PfBb8k+bcHyiGxqCN/fbpI5zwVrmH9B9BjLK/v0wYYIdZNMoHK90Q4PvC21MGsViaV+U1KORI3W6Z1ub35FkroYGmDzZ7ygiJRx9/rW11t8fJXl5wVrolZWlrdIdO6CkxO9oMkqibEXRzm0sHLKTS46q5pLpPW0BbrwSjuRvJXUjY11NPR+s2sn7W17l5S0juWv2scFIBlOn6owzS/6eSZStaGtvp8g5VrkhvPOMFqALxHue4cLR7bN7d7hK6poBSSSDujYY1lRP9d4mvvfMKp5dUe13aJr8revRU4myFQWtzdQNHsruwcU0tbYzb5Ht5ZsO4Wj5NzREr9snIDXZe9Lveu29SCSD7OwchrRoad9EMvC9JThhgjZAgrjfQEglylOUHNzP4okzPp3Rt22vlXVOh3C0/NvaojXVMyA12XsyoHrtvUg8lkPIcoda2YFIBjk5usVkkKaghlyi9IPgWF824dPbx5QU+BVSpIQj+cdiSSf/dTX1PPLGFn62+EMeeWNLUkkq5XqqyR4Q/ao02UeJZOAEsjtM7Q1MMpgxQxcbGk+cNrGUglg7zdl5bC0uA6AgN5u7ZgdsP4cMFZ7kn4RUtFJTKiA12Xsy4HrtPUhUguzY8g9UMpgyRctO2wHAE5PLhzK7LIu1U2fgsrKpKCngvkun+d/FFxHh6PPPzh7QIq+OfdKd9WlXIb8UF3ed6NNck70nRfm5Xb6uPdZr70Xivfjz+u20tWZRUVIQnNk+oH3+s2bBH/4A48Z5PuYRObEYlcPyuff+m7l31Ci/o4mcjE3+n04j6+GsIZlWakrNnKl9/B27fnyoyd6Tflea7KPJ5UOZPCwPWlu57u5zkg3Te6efDs8/z7rqvSxev+vTvz9xNgnYAaCvdu+G446DTEz8LS36/W1r0/yVk6ONhwCNXYYj+efl9buwVld90p0l00pNqcSsngDP9klppclYLLgzasrKYNo01jz5Km05h9ehCfTZZBDt3w+zZ/sdRb89u6KaeYvWs21vE2NKCrhr1tFcMrwNqqth/XrYsAF27tREL6INV+d0X4ijj9aVzGPH6gyy/Hzf/o5wJP9hw+CTT2DIkD7/l95a9V60UlMqIDXZe5KySpMtLVBe7v3jeuXii4k9+gcoHnJESy6wZ5NBU1enCXDqVL8j6ZdnV1TzvWdW0dTaTmFzIxNWrGDf//s3Pho9iMrhBZrMi4pg/PjDPxvO6ed6/Xp49139XeJs/swzYfTotP8t4Uj+I0bApk39+i/d9Uknfmf9s93zvS+7pUXf86CaOJFPJkyh7JNN1BYOO+xXgT2bDBLndD/uG27QLpEQmbdoPbkN+7jww7eYsXU14qCuoIglBwZx/fSjuv+PItrtM2jQoc92SwssWgQLF2r31+WXQ2VlWv4OCFPy7+cMi+76pHvdSi7iOo+V+NKX3dwcnFr+XRFhyi3X0PR//p6sWIxYvBJl4M8mg6K2Vrs+jjvO70j6xzlGrlnBzauWkNPexs7C4bRnxQ9eAznjy8uDceP0YLhlC/zwhzBnDlx4YVrK2YRjqufw4f0e8B3QHqImJfP3+629PfC7tp133gyO+vJsKlt1urB9vvqovR0OHICvfjVQg5+9qquDn/+cb33wEvvzCtg+dMShxE+SZ3wiOpZUUQHPPgv//M9p2cc6HC3/oqIBfVAisfuRx1Ixf7/fsrJCsZ3fSX/7TU7auUk3IOnHeFSkVVfDOefAxIl+R9J3NTUwbx7U1zOh6ng2ras9bO2RZ2d8ubk6CLx7N9xzD9xxB5xwQvKP241wtPyHDg1XKyHEumvBpLUvW0Tf86ArKYHrrtNSz7bZUO/27dP39bLL/I6k77Ztg3vv1a7IsWOZPKY49T0KpaV65vvAA7B8uXeP24knyV9EzhOR9SKyUUTu7uL3Z4rIuyLSJiL9f+fLyvRIa1+wlEussu0orX3Zzul7HZa531VVcPLJ2qI13Wtv1xbtjTeG5yypthZ+8hP9THYYg5pcPpTrT5/A38w6hutPn5Ca3oUhQzTvPfggrF7t/ePjQfIXkWzgIeB8YCpwhYh0nr/1MXAt8D8DepKCAn0hGhuTiNT0he9jJU1N+kUbPDg9z5csEbjqKv2MBqj8RqA4p1O1Z80Kz9TOtjZ46CE4eNC/mWeDB+sZwL/9W0oKCnrR8j8Z2Oic2+ycawEWAHM63sE595FzbiUw8CI9xx6rpZ1NyqWlZdOdhgaYNCl9z+eFkhLtn927V5OFOVxNjb6nc+f6HUnfLVqkM3D8PgMtKtIz4f/6L897PrxI/hXAJx2ub43f5q1Jk6ygVhQcPKgH+rA5+mjt/6+u1i4Oo/bs0RbsLbeEZze+rVvh6ad19k0QxhrLy2HFCvjTnzx9WC+Sf1evzoAOUSJyo4gsE5FltbW1h/8yKG+ESS0Rfa/D6PTT4fzzdZqejU/plM4DB+Bv/iY821/GYvDII7pSNygHKxFdAfz44552LXqR/LcC4zpcHwtsG8gDOefmO+eqnHNVIzsv8kks97et9DJX4r0NcmmHnojo/PVTT4WPPor2AaCxUQdMb789ratWk7ZxI2zeHLxFhoMHa8+Hh61/L5L/UmCSiEwQkTxgLvCcB497uIIC7Q7Yu9fzhzYBsW+fdu+FZbC3Kzk58M1v6q5fUT0DaGzUwma33gqf/azf0fTP4sWaa4LYyzByJLz0kg5GeyDp5O+cawNuBRYBa4EnnXOrReQeEbkYQERmiMhW4HLglyIysLlLp55qg76ZrL5e3+Owy82Fm2/Wnb8++ihaZ6v792viv/12nQYbJnV1sHRp8Fr9CYMH63dkzRpPHs6TFb7OuYXAwk63/VOHy0vR7qDkTJmSeMBgHpnNwCVayIn3OOxyc+Fb39LVv4sXawXLQYP8jiq1amu1VXrXXeGr2wOa+CHYxeYGD9bPkwcrf8OxwjdhxAgd+DhwwO9IjNcOHNAWV1mZ35F4JycHvv51+MY3dLpjpp61JubxFxZqcbIwJn6AVauCX1Zk2DBYt86Ts8lw1PZJEIHTToNnntEPmo98L3ucaerq4OKLM++MTkRr2ZSXw89/rge5UaMy5+9sadGpkSeeqKt3ff5eDphzOtgb8IKC5OTo2VVtbdJrEMLV8gfd4MTnUg+h2xA+6JzTufEpLGLluylTtFjXhAm6eCgTFoPt3Anbt8PXvqZ9/GFN/KDrEZqbg7uDXGc1NUk/RLha/qA75Bx1lM4M8WnucE9lj631PwD19donHqYpgQMxcqT2h7/2GvzmN1q9dPTo8J0FtLToYrbKSt2QZWzyw3m+q6nx/n1YtSo1W7GK6ESCE09M6mHC1/IXgQsu8HXKZyDKHmeSujp9T8OWBAciKwu++EX48Y91WuuWLVq3JQxTQtvatG9/505dz/CP/5gZiR90eqqX78GqVfD884cWZe3bp9dXrUr+sfPy9EwlSeFr+YPOHS4o0BaID6vwutsi0rbwG4CWFp0Fc9JJfkeSXmVl8J3v6ODdggV6ECgt1RZi0MRi2jJub9fibBdcEJ4Vu33V2upt8l+yRB+z83MsWZJ86z8ry5NSN+Fr+YMuvZ45U+uo+8D3sseZZMcOHRDNz/c7kvQT0bGAH/xASyAMGqQLw7ZvD8bagKYm+PhjHdCdMQPuuw+uvDLzEj94f9bZXRkGr8ozeBBvOFv+oHVUXnhBvyRZfT+GeTFLJ3F/m+2TpMRA75ln+h2Jv7Ky9MznxBNhwwZtHSbmnA8bNuCd7AakrU27FBob9XkvvVQX3pVmeMMmx+NUWFzcdaL34swuFvNkzUh4k//o0bqCcOVKGDOmT//Fy83JbYtID2zfrkkvrLV8vJaVpSVMjj1WE/DSpfDGG9r6dk5n0xQXezsjxTlt4dfV6YE4J0e7VU87Tefrh2X2S7IKC709wM6cqX38Hbt+cnP19mS1tHiyx0B4kz9oq2T5cv3Q9mFVns3SCZD2du23vPRSvyMJpuHDYfZs/dmzR8cG3nlH/00kFBFd8ZmXpy3B3NzuE1h7uyaNlhadZtrUpAcb57RVf/bZMH267q0blGqW6VRe7m2ff6JfPxWzfZzTWY9JCnfyHzMGzjpLW0d9mHVgs3QCpKYGzjgjc2aLpNLw4fCFL+hPLKazg2pqdHxgyxY9OOzZo3V1OkocCGIxbdGXlOhAc1mZJvmxYzXphWVbxVQqKdHXwctJJNOmeZPsu9LH3o6ehDv5A1x0Ebz+uraGejlFtVk6AdHaqglpzpze72sOl5Wl6wVGjjxyUVx7u64gbm/X11dE75+XF9xKlUEhohvybNoU7PGNtjbNcx7EGM7ZPh2VlurUsz6seLNZOgFRUwPnneff3qiZKjsbhg7VQeLSUj1jKCnRriFL/L074YQjz56CZs8eOP74fk1y6U74kz/Al76kUwV7Kfjm++bkRmeRDBqkyd+YIKmq0qQa5G04Gxu9GTQmE7p9QKekXXstPPig1k7poZVjs3R85JzO8Pn2t7WFakyQDB2q01r//GdP+tQ9t3+/ns15tMd1ZrT8QRehzJgB2wa0g6RJh23btHV1yil+R2JM1845Rwd9g1huY/duuPBCz/YbyJzkLwJXX62zGqzef/A0NuqH9uqrrf/ZBNeECdqy9ql6QLf279exGw8bTpmT/EEHt77xDX3jgnjkjqpEd8+11+pgpDFBJQLXXaezajyon+OJWEyL6V13nadlszMr+YN2/Zx8stYjMcFQXa3vi3X3mDAYNQquuEI/t0FoRG7bpmMRHhc/zLzkL6Kt/7IyPVoaf9XW6pTO666z7h4THl/8ohbd82DTlKTs26czGa+80vPvT+Ylf9CVenfcoadLmbpvahg0NOjp8x132CpSEy7Z2XDzzdpN6Vf/f0OD/vzt36ak1HdmJn/Qwm+33Qa7dumxyyiiAAANQElEQVTovUmvlhZt9d92mxVuM+FUUgLf/a4OtG7fnt7nrq/XDavuvFNLcaRA5iZ/0KqEV1+tuw8FoT56VMRi+pr/9V/rakRjwmrkSPj7v9czgI8/TkkeWVdTzyNvbOFniz/kkdc3s3n1Zp2xePfdMHWq58+XkNnJH3Q13KxZuuelHQBSLxbT13rmTF15bUzYlZbC97+ve4hs2aKtco8kysw3HGwlt62Vkh3V/H6X8OJlN+k2nymUGSt8eyICV12lxcRee63XFcCm/z7dIKephWOadlN+wUym23x+k0kKC3UiSVUVPPywngWMHp10BdA3N+0m1tZG2YG9ZLsYzxx3Nm8ddQLly+s4/1yPYu9G5rf8QQdvrr1Wp0vZGYCnPm25NLUwdt9O3iqdyNXueJ5dmeY+UmNSTUSLv917L5x/vq64/egj7ZsfyJTQxkaKdlQzav9uPhg9kf975tW8OWE6saxstu1t8jz8zjK/5Z+QkwM33KD/vv46VFZ6Uhkv6t7ctJv29jbG7dvJsoqp/O6Ec2lrh3mL1nPJ9Aq/wzPGe4WFcPnlWk7+3Xdh4UI9EwBtaBYV6X2ysg7fU6GxUWfvtLTo7UVFvF01k8XFn6E+//DFW2NKClL+Z0Qn+YMm/uuu039ffll3w4nKNnUp0nigibENu3hn3DSenjaT9iytO5KOlosxvsrP1w12Tj1VZwNt26ZjAuvW6UY7ra2a5J3TA0F5uZ45TJqkheMqKjhz1Q5efGYVtB6qJFqQm81ds70p3taTaCV/0CPzNdfoKP6TT+pqPpuDPjCNjRzdXMdTk0/n5YkzcHLoTCodLRdjAkFEE3t5OXzuc4duj8V0nUt29uFnAR0kzo7nLVrPtr1NjCkp4K7Zx6blrDl6yR/0jfirv4KKCvjFL3RP0yDv3hNEu3dDUxP537mTP60TnA8tF2MCLbGLWi8umV7hSxepJ53eInKeiKwXkY0icncXvx8kIr+N//4dEan04nmTNn06/OAH+gYFpY5H0Dmnr1VuLvzgB5z91+dz36XTqCgpQICKkgLuu3Sa9fcbE3Dikkx4IpINfAicC2wFlgJXOOfWdLjPt4ETnHM3ichc4MvOua/19LhVVVVu2bJlScXWZ/X18O//DqtX69nAoEHped6waW7WxD91qi59T8GSc2NMckRkuXOuqrf7edHyPxnY6Jzb7JxrARYAnXfmngM8Fr/8FDBTJECTwIcO1WXcV1+txeCsJPThnNPXZOdOXbX73e9a4jcm5Lzo868APulwfSvQuXbvp/dxzrWJyD6gFNjV8U4iciNwI8D48eM9CK0fsrPh3HO1HMEjj8CHH9pZABxq7U+aBNdfH8zt7Ywx/eZF8u+qBd+52dyX++Ccmw/MB+32ST60ASgvh+99T6eCLligB4XRo6O3JiAW0+lr7e3a2p85U6fIGmMyghff5q3AuA7XxwKdN9JN3GeriOQAxcAeD547NRJnAdOmwdNP64bOQ4bo9NAA9ValhHNajfPAAZ22dvnlVpXTmAzkRfJfCkwSkQlANTAXuLLTfZ4DrgHeAi4DXnbJjjSnw+jRcMstupT7ySdhzRot8zpsWGYeBOrqdKn65Mnw1a+mrJSsMcZ/SSf/eB/+rcAiIBt41Dm3WkTuAZY5554DHgEeF5GNaIt/brLPm1af+Qz83d/B2rXwxBNaz6OwUHeoCvtBwDmds9/QAGPHwo036myesP9dxpgeJT3VM1XSOtWzP9rb9Qxg0SL44APtIiorC9/AcEuLzt5pa9N9D2bP1n+zs/2OzBiThL5O9bQRvP7KztaxgGnTdED0zTfhj3/UVcJFRdolFNTB4VhMu3YaGvRgNXs2nHGGdm8ZYyLFkn8yRo+Gr3wFLrwQ3n9fq4WuWaNdKXl5WjIiyXrfSWtpgT17dMqmiG5KfcYZcOKJWpjKGBNJlvy9kJ8Pp5yiP42NsGEDLFumP83NejAYPFjHCfLzU9ef7pw+X0ODxgHawp8xQzehmDRJ4zDGRJ4lf68NHgyf/az+XHONlnb9+GMt87phg/azZ2VpF0xuribnvLxDP70dGJzTUrHNzdqqT/wkSseWlOhCtWOP1ZLVlZU2P98YcwTLCqmUk6PTJSdOhLPP1tv279e639XVOmZQW6vdMnV1el2k+wOAc/qTGFsoLdUZR+Xluhq5vFx/Z4wxvbDkn26FhXDMMfrTWXu7HhxaW/XMoD1eJjlRDzw3V/+/zcgxxiTJkn+QZGdbwTRjTFoEdE6iMcaYVLLkb4wxEWTJ3xhjIsiSvzHGRJAlf2OMiSBL/sYYE0GW/I0xJoIs+RtjTARZ8jfGmAiy5G+MMRFkyd8YYyLIkr8xxkSQJX9jjIkgS/7GGBNBlvyNMSaCLPkbY0wEWfI3xpgIsuRvjDERZMnfGGMiyJK/McZEkCV/Y4yJIEv+xhgTQUklfxEZLiJ/FJEN8X+HdXO/l0Rkr4i8kMzzGWOM8UayLf+7gSXOuUnAkvj1rswDrk7yuYwxxngk2eQ/B3gsfvkx4JKu7uScWwI0JPlcxhhjPJJs8h/lnKsBiP9blnxIxhhjUi2ntzuIyGJgdBe/+r7XwYjIjcCNAOPHj/f64Y0xxsT1mvydc7O6+52I7BCRcudcjYiUAzuTCcY5Nx+YD1BVVeWSeSxjjDHdS7bb5zngmvjla4DfJ/l4xhhj0iDZ5H8/cK6IbADOjV9HRKpE5OHEnUTkdeB3wEwR2Sois5N8XmOMMUnotdunJ8653cDMLm5fBnyzw/UzknkeY4wx3rIVvsYYE0GW/I0xJoIs+RtjTARZ8jfGmAiy5G+MMRFkyd8YYyLIkr8xxkSQJX9jjIkgS/7GGBNBlvyNMSaCLPkbY0wEWfI3xpgIsuRvjDERZMnfGGMiSJwL5oZZIlIL/MWDhxoB7PLgcbwW1LjAYhuIoMYFwY0tqHFBuGM7yjk3srcHCWzy94qILHPOVfkdR2dBjQsstoEIalwQ3NiCGhdEIzbr9jHGmAiy5G+MMREUheQ/3+8AuhHUuMBiG4igxgXBjS2ocUEEYsv4Pn9jjDFHikLL3xhjTCcZkfxF5DwRWS8iG0Xk7i5+f6eIrBGRlSKyRESOClBsN4nIKhF5T0TeEJGpQYmtw/0uExEnImmZ/dCH1+xaEamNv2bvicg30xFXX2KL3+er8c/bahH5nyDEJSI/7fB6fSgie9MRVx9jGy8ir4jIivh39IIAxXZUPGesFJFXRWRsmuJ6VER2isgH3fxeROTn8bhXishJ/X4S51yof4BsYBPwGSAPeB+Y2uk+ZwOD45dvBn4boNiGdrh8MfBSUGKL368IeA14G6gKQlzAtcCDAf2sTQJWAMPi18uCEFen+98GPBqg12w+cHP88lTgowDF9jvgmvjlc4DH0xTbmcBJwAfd/P4C4EVAgM8D7/T3OTKh5X8ysNE5t9k51wIsAOZ0vINz7hXnXGP86ttAWo7efYytvsPVIUC6BmF6jS3uR8BPgIMBi8sPfYntBuAh51wdgHNuZ0Di6ugK4Ik0xAV9i80BQ+OXi4FtAYptKrAkfvmVLn6fEs6514A9PdxlDvBfTr0NlIhIeX+eIxOSfwXwSYfrW+O3ded69IiZDn2KTURuEZFNaJK9PSixich0YJxz7oU0xdSnuOK+Ej/dfUpExqUntD7FdgxwjIi8KSJvi8h5AYkL0G4MYALwchrigr7F9kPgKhHZCixEz0zSoS+xvQ98JX75y0CRiJSmIbbe9DfvHSETkr90cVuXrWcRuQqoAualNKIOT9nFbUfE5px7yDk3Efg74B9SHpXqMTYRyQJ+CnwnTfF8+tRd3Nb5NXseqHTOnQAsBh5LeVSqL7HloF0/X0Rb2A+LSEkA4kqYCzzlnGtPYTwd9SW2K4BfO+fGot0Zj8c/f6nWl9i+C5wlIiuAs4BqoC3VgfVBf97zLmVC8t8KdGz5jaWL00YRmQV8H7jYOdccpNg6WABcktKIDukttiLgeOBVEfkI7Vd8Lg2Dvr2+Zs653R3ew/8EPpfimPocW/w+v3fOtTrntgDr0YOB33ElzCV9XT7Qt9iuB54EcM69BeSj9Wt8j805t805d6lzbjqaP3DO7UtDbL3pb245UjoGL1I8MJIDbEZPZRODNsd1us90dGBnUgBjm9Th8kXAsqDE1un+r5KeAd++vGblHS5/GXg7KK8ZcB7wWPzyCPTUvNTvuOL3Oxb4iPj6ngC9Zi8C18YvT4knsZTH2MfYRgBZ8cs/Bu5J42tXSfcDvhdy+IDvn/v9+On6Q1L8Il0AfBhP8N+P33YP2soH7RrYAbwX/3kuQLH9K7A6HtcrPSXgdMfW6b5pSf59fM3ui79m78dfs8lBec3iX8YHgDXAKmBuEOKKX/8hcH+6Xqt+vGZTgTfj7+d7wJcCFNtlwIb4fR4GBqUprieAGqAVbeVfD9wE3NThc/ZQPO5VA/lu2gpfY4yJoEzo8zfGGNNPlvyNMSaCLPkbY0wEWfI3xpgIsuRvjDERZMnfGGMiyJK/McZEkCV/Y4yJoP8PKIUBLvKuKv0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe4d30af198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 47.7 s, sys: 20 s, total: 1min 7s\n",
      "Wall time: 2.71 s\n"
     ]
    }
   ],
   "source": [
    "def minAreaKCircleCover(P, k):\n",
    "    M, R, X, S = basicModel(k, P, bigM(P))\n",
    "    t = M.variable(1, Domain.greaterThan(0.0))             \n",
    "    M.constraint(Expr.vstack(t, R), Domain.inQCone())        #t>=sqrt(sum r_i^2)\n",
    "    M.constraint(Expr.sum(S,1), Domain.greaterThan(1.0))\n",
    "    M.objective(ObjectiveSense.Minimize, t)\n",
    "    M.writeTask('p4.task.gz')\n",        
    "    M.solve()\n",
    "    return R.level(), X.level()\n",
    "\n",
    "%time P=data(20); r, x = minAreaKCircleCover(P, 3); display(P, r, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 5. Maximum k-coverage\n",
    "\n",
    "**Objective: find the largest subset which can be covered by $k$ balls of a fixed radius.**\n",
    "\n",
    "We can combine the previous techniques to express more complicated problems, like covering the largest possible number of points with $k$ fixed-power transmitters. Here we want to maximize $\\sum_i\\max_j{s_{ij}}$, which is precisely the number of covered points, and we want to have an upper bound on $r_j$. The objective is equivalent to:\n",
    "\n",
    "$$\\max\\sum_i t_i \\quad \\textrm{subject to}\\quad t_i\\in\\{0,1\\},\\ \\ t_i \\leq \\sum_js_{ij}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt8VPWZ+PHPkxtJIBAgQSBcQhG5KK1osK5U6x3UqnhrcdfWW9e1rdV9bddV17a7a7erW3druyuvrvzUtdvWel+XtlhaEWvFS4miIggKiEJACAkhAZKQZL6/P54ZMoQJmWQu55w5z/v1you5nMx5Msw853u+3+95vuKcwxhjTLjkeR2AMcaY7LPkb4wxIWTJ3xhjQsiSvzHGhJAlf2OMCSFL/sYYE0KW/I0xJoQs+RtjTAhZ8jfGmBAq8DqA3lRUVLjq6mqvwzDGmEB54403djnnKvvazrfJv7q6mtraWq/DMMaYQBGRj5LZLi3dPiIyT0TWi8gGEbn9CNtdLiJORGrSsV9jjDEDk3LyF5F8YCFwHjADuFJEZiTYrgy4GXg91X0aY4xJTTpa/icBG5xzm5xzB4DHgIsTbPc94AdAWxr2aYwxJgXpSP5VwJa4+1ujjx0kIrOA8c65X6dhf8YYY1KUjuQvCR47uEiAiOQB9wHf6vOFRG4QkVoRqa2vr09DaMYYYxJJR/LfCoyPuz8O2BZ3vww4DnhRRDYDJwOLEw36OucWOedqnHM1lZV9zlQyxhgzQOlI/iuBKSIySUSKgAXA4tiTzrk9zrkK51y1c64aeA24yDln8ziNMcYjKSd/51wncBOwFHgPeMI5t0ZE7hKRi1J9fWOMMemXlou8nHNLgCU9HvtuL9ueno59GmOMGTir7WOMMSFkyd8YY0LIkr8xxoSQJX9jjAkhS/7GGBNClvyNMSaELPkbY0wIWfI3xpgQsuRvjDEhZMnfGGNCyJK/McaEkCV/Y4wJIUv+xhgTQpb8jTEmhCz5G2NMCFnyN8aYELLkb4wxIWTJ3xhjQsiSvzHGhJAlf2OMCSFL/sYYE0KW/I0xJoQs+RtjTAhZ8jfGmBCy5G+MMSFkyd8YY0LIkr8xxoSQJX9jjAkhS/7GGBNClvyNMSaE0pL8RWSeiKwXkQ0icnuC528UkdUi8paIvCwiM9KxX2OMMQOTcvIXkXxgIXAeMAO4MkFyf9Q5N9M5dzzwA+CHqe7XGGPMwKWj5X8SsME5t8k5dwB4DLg4fgPnXHPc3cGAS8N+jTHGDFBBGl6jCtgSd38r8NmeG4nIN4C/AYqAM9OwX2OMMQOUjpa/JHjssJa9c26hc24ycBvw7YQvJHKDiNSKSG19fX0aQjPGGJNIOpL/VmB83P1xwLYjbP8YMD/RE865Rc65GudcTWVlZRpCM8YYk0g6kv9KYIqITBKRImABsDh+AxGZEnf3AuCDNOzXGGPMAKXc5++c6xSRm4ClQD7wsHNujYjcBdQ65xYDN4nI2UAHsBu4OtX9GmOMGbh0DPjinFsCLOnx2Hfjbt+Sjv0YY4xJD7vC1xhjQsiSvzHGhJAlf2OMCSFL/sYYE0KW/I0xJoQs+RtjTAhZ8jfGmBCy5G+MMSFkyd8YY0LIkr8xxoRQWso7BF4kAvv2wd690NWl9/PyID8fhgyBwYP1vjHG5IhwJX/nYM8e2LoVPvwQ1qyBTz7Rx0S6f5zr/jf2M2wYjB7NqyWjWfhhF+8whLKjKrh13jTmz6ry+i8zxph+yf3k39nJC08vZ+VjS6jcsolRtDNzXDkTR5TC0KFQUqKJXRKtSRPlHBw4wMZ3PqBu7XIujES4yEHLoFLefHkiQ6++iDMvPR0Kgvt2PruqjnuXrmdbUytjy0u4de5UO6gZk8OCm6360tQEr77Kxv95iqa1HzM2r5CWQaVsyC9lc2M+Zx9VzrTyocm9lggMGsSLOztpGTJSH3OOoq4OpmzbwP5//Td45ddw3nnw2c9CeXnm/q4MeHZVHXc8s5rWji4A6ppaueOZ1QB2ADAmR+VW8ncONmyAZcvgT38C53h12352lR26KlhnJMKKjQ1MG5Nk8o9qaevoviPCgYIiGgqKaAAoKoLHHoPHH4eTT4Yzz4TJk498RuET9y5dfzDxx7R2dHHv0vWW/I3JUbmR/J2DtWvhF7+Abdtg0CCoqoL8fHZ98H7CXzkkkSeprLgw4e+VFRdCaSlMnKgDxm+8Aa+8AmPHwlVXwfTpvj4IbGtq7dfjxpjgC37yb26GJ56Al16CESOguvqQp4+YsPtpzuSRPP/eTjojkYOPFeTlMWfyyO6N8vNhzJjuweV77oHTToMvfQnKyvq9z2wYW15CXYJEP7a8xINoTOjFvjt79kBLS/dPR0f3TLziYh2zKyvTn+HDtQFmkhbc5O8crFwJP/0ptLVp0k8wHTOphJ2kWDfRio0NtLR1UFZcyJzJIxN3H4lo3//QofDaa/DWW3DNNXDiib47C7h17tRD+vwBSgrzuXXuVA+jMqHR1qZn7B9/rGfw69drso99T2Lf3dj327nug0D8NqNHw4wZMHUqjBsHRx2ljTGTkDjnvI4hoZqaGldbW5v4yZYWePhh7V456iidh38E67Y3J5ewM2nvXti5E2bPhmuv1esHfMRm+5is6urSJP/738Pbb+tjzmnrfehQ7brtD+dg/37tCThwQA8KhYVw+ul65l0Vns+yiLzhnKvpc7vAJf+6OvjRj2D3bv0P9VErus+DjHN6jUFFBdxyi44JGBMmzc16xv6b30Bjoyb7iorMXETZ0QE7dkBnJxxzjM7GO+44PSjksNxM/u+8A/ffrzNrKiq8CawX67Y3J+xeOnv6qMPPMurr9YP5zW/CzJlZjtSYLHMONm+G5cthxQrtoqmo6POMPa37b2zUHoMhQ2DePPizP4OR/e/6DYLcS/6vvAIPPACVlb7rMgF46OUPex1Yvv5zkw7/hb179SBw4436QTQmFzU2ws9/Dm++qY22UaO8vRiytVW/dyJw4YVw/vkaVw5JNvkHY8D35Zdh0SLtJiku9jqahHqbOtrrlNIhQ3Qw6ic/0ZbJKadkMDpjsiwS0e/tL36htydO9EcXbUkJTJigZ97PPquTMb76VTj6aK8jyzr/J/+VK32f+GGAU0pLSvTveuABbX3U9HmwNsb/PvkE/vu/4b33dNpziQ+nDBcW6gzBxkb43vfg3HPhkktCNV3U36UqN23SlvGYMb5O/KBTSgt6DFolNaW0uFhnLP3kJ9ovakxQOafds9/+Nnz0EUya5M/EH2/ECD0TWLYMvvMd2LLF64iyxr/Jv7NTZ/XEiq/53LQxQzl7+qiDLf2y4sLEg72JlJZqN1BsFpMxQdPZCU8+Cf/1XzqQOnq0P7p5kpGfrweAAwfgn/5Jp5CHgH8HfKuqXO0XvqCt/rDYtk0/hLfdlvPT0UwOaW3VpP/WW/r5DfKFVfv3a7fVFVfAF74QnANYnOAP+La3a+shTMaMgQ8+gOef1znJxhfsArgj2LcPfvxj/dxWVwcyWR6itBTGj9eSMW1tcPnlwf+beuHf5F9YmLNveq9E9MK1p56C449P6azHElZ6WLnrI9i/X7sqN23SFn+ufF8LC3V20uLFej9HDwBp6fMXkXkisl5ENojI7Qme/xsRWSsi74jIMhGZmMSLpiO04Ckq0g/fI4901zTpp1jCqmtqxdGdsJ5dVZfWUMPgSOWuQ62jAxYuhI0btY5Orn1fCwq6DwBLlngdTUaknPxFJB9YCJwHzACuFJEZPTZbBdQ45z4NPAX8INX95rSjjtJpci+/PKBft4SVPlbuOgHndN2Kd9/VLpJcS/wxBQV6RvP44zqekWPS0fI/CdjgnNvknDsAPAZcHL+Bc265c25/9O5rwLg07Dd3iWiXz6OPai2UfrKElT69lbUOdbnrP/wBli7Nra6e3hQWamNs4UKty5VD0pH8q4D4ybFbo4/15nrguTTsN7eVlOig95/+1O9ftYSVPrfOnUpJ4aGzV0Jd7nrzZi2jPm5csGf19MfgwVpl9Mc/1plNOSIdyT/RoT/h/FERuQqoAe7t5fkbRKRWRGrr9+9PtEm4VFRof2NXV9/bxrGElT7zZ1Vx96UzqSovQYCq8hLuvnRmOAd7OzrgwQe7k2GYVFTArl3wq195HUnapGO2z1ZgfNz9ccC2nhuJyNnAncDnnXPtiV7IObcIWARQM3asPy9AyKbBg/VKyXXr4Nhjk/61WGKy2T7pMX9Wlb13oLX3t2zRK3fDqKpKS1HPnp0T70E6kv9KYIqITALqgAXAn8dvICKzgAeAec65nWnYZ3iUlMDvftev5A+WsEyabd8OTz8dqkVRDlNQoEtGPvQQ/MM/BP5CzJSTv3OuU0RuApYC+cDDzrk1InIXUOucW4x28wwBnhQdIPrYOXdRqvsOhYoKXemovl7LWec656CpSa92bmrSchc7d+q/nZ3aBZafr90OI0dqieBhw7RGS1VV9mrEh0kkolOPCwtzrvxxv1VUwIcfaoPsggu8jiYlabnIyzm3BFjS47Hvxt0+Ox37CaW8PE12L7+sVQdzTWenrt0av35rc7POInFO//6iIv3puYbrunU6KA66fSSiMzNmzIBp03RQcuzY3J+RkmnvvKPvdXW115H4Q1UV/O//wqmnau2xgPLvFb6m26hR2t960UW5M8Ni92549VX47W91YZvY+q1lZbrw/UASdmwd11degRdf1MeqqnTBjhNO8H1lWN967jlNcnYQVUVFegb6+utwzjleRzNglvyDYNAgrTPyySfB7nONRGDDBq1dtHKlJpPKSu2ySQcR7faJdf04p2cRixbpe3jWWdpaC1OxwFRt26ZnYxP7vii/pz7XtA6yigo9KJ55ZmAbZJb8g8I5Xbw+iMm/rU1bSc89pwew4uLszBMX0fGAYcO0XO/SpTp1dvp0Xcd15szMLByeS156SQc6+9nq77mmdUtbB8+/p3M9cuIAMHiwXvPQz5l4fmLJPyiKirQFdtJJXkfSP+vW6dzwXbu0he9Vv3FRkZYicE6nK/7wh/qlvfpqHSfwOU8K9bW2wgsvaLdjP63Y2HAw8cd0RiKs2NiQG8kftJtyADPx/MKaPUExdKgOiAbF3r06Q+Tuu3VQt7raH4NjInrKXl2t1SjvvFPHUzo7vY6sV54V6nvzTT1jGsAMn36vaR1EFRU6GL4zmLPXLfkHRWmpdpn4/cpn5zRp3HEH/PGPWv+lvNzrqA4noutFVFbCz3+u67j6dAk/Twr1OacXNA0fPqBf723t6iOuaR00eXn6s2KF15EMiCX/oBDRn/p6ryPpXWcn/PKXWuN90CDtZvH7YNigQXoWsGuXXrjzyiua+HzEk0J927bphV0DPFsb8JrWQVNZCcuX++4zkwxL/kEiAi0tXkeRWEsL3HefDqpOnKhrEgdFbNZRZaUuR/jEE77qBvKkUF+sguUAp3emtKZ1kBQXaxdnY6PXkfSbDfgGSSTiz+RfXw/33qtfgCAv5VdSogeuJUu01Xvjjb64NuDWuVMPWU0MslCob/36lK/mnTZmaO4l+0REdCbeyGCd1Vjyj7d6NSxbBnv26PTAs87S6YB+IaIlD/xkxw4d1G1v1+mbQZefrwewt9+G//xPuOkmPSh4yJNCfe+9548B+qDYuBE+/Wmvo+gXS/4xq1drudaO6GyEPXu6y7f65QBQUKCnmH5RX6+Jv6MjENMlkyaiA9Vr18L998PNN3tewjirhfr279fJBRMmZGd/QTd0KKxZE7jyK9bnH7NsWXfij+no0Mf9QsQ/fdEtLfBv/6Yt/lwsOBc7AKxZo1UcB7ieciDV1XVPMDB9GzJEL/jyy3czSZb8Y/bs6d/jXnDOH1/Izk544AFoaMitFn9PIjoG8Oqr8Otfex1N9uzeHcjZK57Jz9fGwQCWXPWSJf+YYcP697gXnNOuH69jePJJ7SYLYqmJ/oqdATz1FLzxhtfRZIffxpWCwM8z8XphyT/mrLMOX5yhsFAf94vOTu8H4d5+W2v0hGHx7pjCQr0g7IEH9HqAXFdfn7W6/eu2N/PQyx/yo+ff56GXP2Td9mC1ng8RsORvA74xsUFdP8/2AW+vlt23Dx5+WPv4M3Txlm8rQZaWaov4Zz+Dv/7r3D7wNTRkJfn3q/ib32fiOeevyRhJsOQfb+ZMf32gesrL03r3Xnn6af2AZ2gWiO8rQY4ZA6tWaTnqoBXY648DB7JS7TTp4m9BmIkXiQRuwNeSf5A4513yf/99rcMfV9c93a1031eCFNEB7kce0ZXCvO6Cy5Surqyc2SRd/O1IM/H8kvwhcDPCrM8/KCIR/UIOoLxuytrbtSzz8OEHu3tirfTYFzXWSk+lzzYQlSAHD9b34/HHvY4kc/LzszLbJ+nib0GYiQeBWxsiWNGmIPADS7HuFi8W0K6t1St54yo8HqmVPlCBqQRZVaUF4LZt8zqSzCgszEorNunib0GYiSfi/yKGPYQi+WeilZp1zc1w3HHZ32+stG+PpRYz0UoPTCXIvDz9or/0kteRZMaIEYd3s2RA0sXfgjATTyRYxQzJ8T7/+D7pnnzVl5yMSAQmT87+fjdt0hZujzVcy4oLE76vqbTSY/8Xvpzt09OoUXz45K+47oNhbN4Xyd7qWtlQWaldW1mQVPG3IMzEs+TvHz1njiTiq77kZHhROO2FF7SuTY8BwDmTRx72/qajlR6USpDrGtpYu2Enw0ve5cNxMw6urgUE/wDgx8V3/D4TD7ydiTcAOZv8E/VJ9+S7vuTe7N2r5WKzXTJ2zx4tbZDgSt5AtdIzYMXGBlxRKWdsrOXNqukgcnB1rZxI/rl8HUOKDpvlNmk405wL3OyvnE3+fbXqfdmX3JuGBvjKV7L/hXz9de3z72UgKyit9ExoaeuAolKqmncyoekTPh4+Bsjw6lrZUlWl/+9+qSXlI4muRXn1nY9xn6lmuheTMVKQswO+R2rVB2pVoc5OTb6zZ2dtl8+uquPUf/k9D/39Qn6+cX+wBsazpKy4EEQ4kF/En3389sHHM7q6VraUlemgb1ub15H4TqIeheK2fTzaEqz+fsjh5N/bzJF5x47m+s9NCkbiB51ieeqpWRtMenZVHXc8s5qurVsZcqCVXZH84M2MyoLY56uxdCjHb3ufvEhX5lfXyqbp0wNXpTIbEvUo5EcirBoUkF6EODmb/HNiDVHndMrd6adnbZf3Ll1Pa0cXo5t3kee0hZPq/P1cFPt8lZYWk+8izMxv5e5LZwa/vz9m+nRozYEurDTrrUfBBXAVu5zt84cc6JNuatIlBbO4olKsz/rohi20F3R/0AM3MyoLDn6+Ps7ni5dMglxJ/ADjxwfuitVs6DnLraCrk66iIv7ykuDVerL/XT9raoLzz8/qoFusz3pKwxZaikoPPh6YmVFeKCyEdeu8jiK9qqp01k/AKlVmWs8ehQmde5ly2XnMP3G8x5H1X1qSv4jME5H1IrJBRG5P8PxpIvKmiHSKyOXp2GfOa2qCigr4zGeyuttb506lItLGsNYW2gt09kKgZkZ5YehQXfA8l+TnwwUX6Ewzc4hpY4Zy/ecm8ddnTeGiY0dR85X5Xoc0ICknfxHJBxYC5wEzgCtFZEaPzT4GrgEeTXV/oRCJQGMjfPWrWV84fP6sKu45uYKS6GyWQI6VZFtpqS54nmt95CedpF0/AStVnDWNjTBlSmBXtEtHn/9JwAbn3CYAEXkMuBhYG9vAObc5+lywap56pa5OB3mnT/dk92dXFcNxow8r6WB6IaJJsrkZSnJgqmdMWRnMmdPrhX6h19IC114b2Gsh0tHtUwVsibu/NfpYagJWGztt9u7VBHLFFd7FYGu49l8A13BNyhln6IwzW9D9UG1tWt7b7yUnjiAdLf9Eh70BfVJE5AbgBoAJgwfrASBMMw4iEdi5E265xds6Ibt2eVM6+kiCsIxfLib/6mo9A9yzx581f7yycyfMn394tdEASUdm3QrED3WPAwZU6Nw5t8g5V+Ocq6msqICPPw5Pi8M5/XtPPx1OOMHbWHbv9lfyjy3jF1u8I7aM3+rV3sYVLxLRNY5zjYjOONu92+tI/CM2BjJnjrdxpCgdyX8lMEVEJolIEbAAWJzyq44cCZ/6lA6khcH27Tp4dNVV3vchdnR4H0O8Iy3j5xfO6fKHueiEE7SirM38UXV1cM452S+0mGYpJ3/nXCdwE7AUeA94wjm3RkTuEpGLAERktohsBa4AHhCRNX2+sAjcdJPOdmlsTDVMf2ts1BkjX/+6P1rcsSUj/SIIy/jFCqHlosJCuP56HdAO+8yfPXt0Rbv5wZzeGS8tHerOuSXOuWOcc5Odc9+PPvZd59zi6O2VzrlxzrnBzrmRzrljk3rhESPgb/9WW3m5Ogi5e7d+ob71rUOWSfRUltZwTVpQlvHL5fGpT31Ku3/q6ryOxDuRiJ79XH99Tszq8v+ndcIE+Lu/0znUudbv2NgIBw7o3zfeR1cIDhrkry6MoCzjF+DBv6RcdJE2UPx0xpVNdXVw2mkwo+dlTMHk/+QPMGkS3HGHtpB37vQ6GiANC8Lv2KEtiTvu0BkVflJRoQclv5g5Ey68sLulP2yY3vfTbJ+8vMAt49dvJSV64WFDQ1bW+PWVPXv07//iF/3VJZqC4BR2q66G734XfvQj2LJFB6A8+k9ItKDD8+/pQanPK2Gd0/jHjNEpnaNGZTrc/quo8N+X2+/L+IkEbhm/AZkxAy67DJ5+WhtlOZIIj6i9XXsd/v7vA7da15EEo+Ufc9RR8J3vaL2bDz/0rHWaaEGHpMoeHzigcc+aBXfe6c/ED/oBz+X+60yIRMKR/EHPuj77WW3E5LquLti6Fa65Bo45xuto0ip43/DSUvjGN2DBAp0GumNH1gcneytv3GvZY+e6Y73ySp3VU1qaeFs/GDHC6wiCpatLD5Z+GoDOpLw8uO46LfmwfbvX0WROJAIffQRz58LnP+91NGkXvOQPOhvl/PPh+9/XgdLNm7O65Fxv5Y0TPt7aqq39iRPhn/8Zzjuv1zVxfWPsWP03rCU2+qulRbtACoLTi5qykhKdoVZenpvX4sQS/5w52mDLwe6tYCb/mDFj4LbbtLhSY6OenmVhHnJvS0QeUva4o0PjaWrSQbLbbtN4g2DQIJ1lZbXck9PSAscmN3s5p5SX6+e6vBy2Deiifn+KRLRBecopeobj98baAAW/qZKfryURZs6E3/wGXnpJDwAjRybsg123vZkVGxtoaeugrLiQOZNH9rtccWz7hK/T0qIHooICnRb2hS8Esxtlxgz43e9yaoArY5zTefBhNGIE3H473HeftpQnTAh2K/nAAR3LOOssvdo+h8/mxPnpYp44NTU1rra2tv+/uHcv1NbCkiU6LbS4WAdW8/IOm6UD2mJPuV59rCBbeztUVmqX1OzZWvUvqN58E+6/P6tLSAbWRx/Bv/97MA/y6dLaCo88Aq+8op+ZIF7z0NKiRQ2//GVN/gE9iInIG865mr62y73D2pAheiZw2mnwwQfw/PPwxhsQibB2XSP5FNGZX3jwPzY2S6dfyd85TfTNzbB/v77W7Nn6gZkyJTdmyowf312yIKBfgqxobdWzI79cne2VkhL4q7/SKdhPPaVdQUF5T5zTbqv8fL3gMiRdeLmX/GPy8mDqVP1paIC1a6n958eZ3LCFkfu6r1BsKyyio71Av8RFRfp7IvqBiET0NPDAAU32sUFl53Rmx2c+o90jxx6be62+igqYPFnPaHLtb0un+nq45BI7QIJ+dy68UL9zDz6o/eZVVf4+C9i3T2fhnXCCtvgDXqytP3I3+ccbORJOPZU/nNHBo7v3M7R9H0e1NDBuzw6O2ttItWvVRP/JJ90LV8Qu1y8v166cigr9IFdXa+tm6NDc/sLHSvn+x39Y8u9NJKKflVNO8ToSfznmGPje93QMbvFi7f6srPQ6qkNFIlquobgYvvlNqKnJ7e9zAuFI/lG3zp3KHc+splmG0Fw8hA8qJ1JSmM/dl86EWXGLj1lXh5o5U7+4bW36JTGHqq/XFmOIWotJGzQILr1Uk+pDD+l056FDtSHh5XcrViKmvR1OPVXLNYR0UkOokv/8aIK/d+l6tjW1Mra8hFvnTj34+EGW+FVhoV7g8uyzNvCbSGur1nU3vZswQa/Kf/ddnYTx/vs6g2bUqOyWL9+3Twdz8/L0TO3MM/UsPsTf9VAlf9ADwGHJ3vTulFPgmWfCt6RmX/bu1a7AHLvkPyMKCuD44/Wnrk6nYy9frmNpw4dnrgs1EtGE39qq+/mLv4CTTgptS7+n0CV/008VFdq18d57WlvJqIYG+MpX7IDYX1VVesXsxRfrdOIlS3T5UhGdbVNWpjP2BvK+dnbqdM3YWsrO6aSMc8/VQegcvVhroCz5m76de65Ol7WxENXerl1is2d7HUlwlZbC5z6n5RPq6/Vq+A0bYM2aQwvGRSKatIuK9LMXm4nX1dVd2DEvTx8rKNCp1jNmaLmN2MQMk5Alf9O3qVN14O6dd7TlFmaxOeFf/nLu1+/PBhHt/x81Ss8wQQ+u27drDf29e/WK+dgaArGDQXGxnpWWl+vZwogR+hrWuk+aJX/TNxHtL333Xe0/zYEl7Aasvl6vfzjjDK8jyV2DBvlvgaMcZB2WJjkjRmitk+3b/bW+bzZ1dOjBL4eLfZnwsORvkjdnjvan7tjhdSTeqKvTq3nHjfM6EmNSZsnfJC8vT1c0ikR03nSY7Nyp4x3z5nkdiTFpYcnf9M/o0fC1r2nr32/r/GZKc7N2dd18c3YvTDKh8OyqOubc8wKTbv8Nc+55gWdX1WVlv5b8Tf/NmgVXXKHzs3O9/7+9XWea3HKLXedg0u7ZVXXc8cxq6ppacUBdUyt3PLM6KwcAm+1jBuaCC3TK4yuvsG7QCFZsakxpgRxfiq3Gdt11MG2a19GYHHTv0vW0dnQd8lhrRxf3Ll2f8UoE1vI3AxNdxPvt0VN47/XVtLTqBTctbR08/95O1m1v9jjAFHV06JnNggW6PoQxGbCtqbVfj6eTJX8zcIWF3Fw2m1Wjjmbcnh2I0xXSYgvkBFZ7e3fiP/98u6rZZMzY8sTXzPT2eDpZ8jeeFGdbAAAJ2UlEQVQp+XhvJ48efx6vTfg045t2UNilg8AtbQEdDG5u1msZrrvOEr/JuFvnTqWk8NBrRkoK87l17tSM79v6/E1KxpaXUNfUylMzz2Z7WQWXrFnO7pIy8oYN8zq0/otdv3D77dbHb7Ii6TLzGWDJ36QktkBOa0cXKybNYseQEVz39hJOKTsQnDLQBw7o4PX48bqq06hRXkdkQsSrMvNpSf4iMg/4MZAPPOicu6fH84OA/wFOBBqALznnNqdj38ZbPVsurVOmMezq06le8wd49VUtvuXXyorOda/qdNllunCNzeM3ISEuxXnaIpIPvA+cA2wFVgJXOufWxm3zdeDTzrkbRWQBcIlz7ktHet2amhpXW1ubUmzGQ85pFdCHH9b66lVV/qqH09amrf0pU7R/P+zVSk3OEJE3nHM1fW2Xjpb/ScAG59ym6I4fAy4G1sZtczHwj9HbTwH3i4i4VI88xr9EdCGNf/kXXQby+ee7y/d6uR5wS4tetFVUBNdeC6ed5q+DkjFZko7kXwXErb7AVuCzvW3jnOsUkT3ASGBX/EYicgNwA8AEWzM2NwwerOWg586FFSvgd7+DTz7RrqDhw7Mzm6arS0sxt7frVbrXXqvrEwwenPl9G+NT6Uj+ib69PVv0yWyDc24RsAi02yf10IxvVFTo0n3nnafdQb/9ra7cVFgIlZXp72t3Tssv74q2L2bPhrPOgqOPDsYgtDEZlo7kvxUYH3d/HLCtl222ikgBMAxoTMO+TdAUFWmr+8QTtXTCH/6gZwTt7fp8Xp6eFfR3HdfY+q1792rid07XILj8cjj5ZL1tjDkoHcl/JTBFRCYBdcAC4M97bLMYuBp4FbgceMH6+0NORKdWXnWVdgvt2nX4Oq7xawZHIoe/RuzgEFu/9eij4dhjbf1WY5KQcvKP9uHfBCxFp3o+7JxbIyJ3AbXOucXAQ8DPRGQD2uJfkOp+TQ4R0a6fysruiqEdHXq17d692qLft0/77mPXDhQW6tlBWVn3vzZwa0zSUp7qmSk21dMYY/ov2ameNvJljDEhZMnfGGNCyJK/McaEkCV/Y4wJIUv+xhgTQpb8jTEmhCz5G2NMCFnyN8aYELLkb4wxIWTJ3xhjQsiSvzHGhJAlf2OMCSFL/sYYE0KW/I0xJoQs+RtjTAhZ8jfGmBCy5G+MMSFkyd8YY0LIkr8xxoSQJX9jjAkhS/7GGBNClvyNMSaELPkbY0wIWfI3xpgQsuRvjDEhZMnfGGNCyJK/McaEkCV/Y4wJIUv+xhgTQiklfxEZISK/F5EPov8O72W734pIk4j8OpX9GWOMSY9UW/63A8ucc1OAZdH7idwLfDnFfRljjEmTVJP/xcBPo7d/CsxPtJFzbhnQkuK+jDHGpEmqyf8o59x2gOi/o1IPyRhjTKYV9LWBiDwPjE7w1J3pDkZEbgBuAJgwYUK6X94YY0xUn8nfOXd2b8+JyA4RGeOc2y4iY4CdqQTjnFsELAKoqalxqbyWMcaY3qXa7bMYuDp6+2rg/1J8PWOMMVmQavK/BzhHRD4AzoneR0RqROTB2EYi8kfgSeAsEdkqInNT3K8xxpgU9NntcyTOuQbgrASP1wJfjbt/air7McYYk152ha8xxoSQJX9jjAkhS/7GGBNClvyNMSaELPkbY0wIWfI3xpgQsuRvjDEhZMnfGGNCyJK/McaEkCV/Y4wJIUv+xhgTQpb8jTEmhCz5G2NMCFnyN8aYEBLn/LlglojUAx+l4aUqgF1peJ1082tcYLENhF/jAv/G5te4INixTXTOVfb1Ir5N/ukiIrXOuRqv4+jJr3GBxTYQfo0L/BubX+OCcMRm3T7GGBNClvyNMSaEwpD8F3kdQC/8GhdYbAPh17jAv7H5NS4IQWw53+dvjDHmcGFo+RtjjOkhJ5K/iMwTkfUiskFEbk/w/N+IyFoReUdElonIRB/FdqOIrBaRt0TkZRGZ4ZfY4ra7XESciGRl9kMS79k1IlIffc/eEpGvZiOuZGKLbvPF6OdtjYg86oe4ROS+uPfrfRFpykZcScY2QUSWi8iq6Hf0fB/FNjGaM94RkRdFZFyW4npYRHaKyLu9PC8i8h/RuN8RkRP6vRPnXKB/gHxgI/ApoAh4G5jRY5szgNLo7a8Bj/sotqFxty8CfuuX2KLblQEvAa8BNX6IC7gGuN+nn7UpwCpgePT+KD/E1WP7bwIP++g9WwR8LXp7BrDZR7E9CVwdvX0m8LMsxXYacALwbi/Pnw88BwhwMvB6f/eRCy3/k4ANzrlNzrkDwGPAxfEbOOeWO+f2R+++BmTl6J1kbM1xdwcD2RqE6TO2qO8BPwDafBaXF5KJ7S+Bhc653QDOuZ0+iSvelcAvsxAXJBebA4ZGbw8DtvkothnAsujt5Qmezwjn3EtA4xE2uRj4H6deA8pFZEx/9pELyb8K2BJ3f2v0sd5cjx4xsyGp2ETkGyKyEU2yN/slNhGZBYx3zv06SzElFVfUZdHT3adEZHx2QksqtmOAY0RkhYi8JiLzfBIXoN0YwCTghSzEBcnF9o/AVSKyFViCnplkQzKxvQ1cFr19CVAmIiOzEFtf+pv3DpMLyV8SPJaw9SwiVwE1wL0ZjShulwkeOyw259xC59xk4Dbg2xmPSh0xNhHJA+4DvpWleA7uOsFjPd+zXwHVzrlPA88DP814VCqZ2ArQrp/T0Rb2gyJS7oO4YhYATznnujIYT7xkYrsSeMQ5Nw7tzvhZ9POXacnE9rfA50VkFfB5oA7ozHRgSejP/3lCuZD8twLxLb9xJDhtFJGzgTuBi5xz7X6KLc5jwPyMRtStr9jKgOOAF0VkM9qvuDgLg759vmfOuYa4/8P/B5yY4ZiSji26zf855zqccx8C69GDgddxxSwge10+kFxs1wNPADjnXgWK0fo1nsfmnNvmnLvUOTcLzR845/ZkIba+9De3HC4bgxcZHhgpADahp7KxQZtje2wzCx3YmeLD2KbE3b4QqPVLbD22f5HsDPgm856Nibt9CfCaX94zYB7w0+jtCvTUfKTXcUW3mwpsJnp9j4/es+eAa6K3p0eTWMZjTDK2CiAvevv7wF1ZfO+q6X3A9wIOHfD9U79fP1t/SIbfpPOB96MJ/s7oY3ehrXzQroEdwFvRn8U+iu3HwJpoXMuPlICzHVuPbbOS/JN8z+6OvmdvR9+zaX55z6Jfxh8Ca4HVwAI/xBW9/4/APdl6r/rxns0AVkT/P98CzvVRbJcDH0S3eRAYlKW4fglsBzrQVv71wI3AjXGfs4XRuFcP5LtpV/gaY0wI5UKfvzHGmH6y5G+MMSFkyd8YY0LIkr8xxoSQJX9jjAkhS/7GGBNClvyNMSaELPkbY0wI/X+2MHePB7wgogAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe4cb179a20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3min 40s, sys: 1min 43s, total: 5min 23s\n",
      "Wall time: 10.7 s\n"
     ]
    }
   ],
   "source": [
    "def maxKCoverage(P, k, rMax):\n",
    "    M, R, X, S = basicModel(k, P, bigM(P))\n",
    "    M.constraint(R, Domain.equalsTo(rMax))\n",
    "    t = M.variable(\"t\", len(P), Domain.binary())\n",
    "    M.constraint(Expr.sub(t, Expr.sum(S,1)), Domain.lessThan(0.0))     #contains one inequality for each column\n",
    "    M.objective(ObjectiveSense.Maximize, Expr.sum(t))\n",
    "    M.setSolverParam('mioConicOuterApproximation', 'on')\n",
    "    M.writeTask('p5.task.gz')\n",    
    "    M.solve()\n",
    "    return R.level(), X.level()\n",
    "\n",
    "%time P=data(20); r, x = maxKCoverage(P, 3, 0.1); display(P, r, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comments\n",
    "\n",
    "For simplicity of the presentation we ignored various issues and hints and we leave them as an exercise:\n",
    "* Checking that a solution has actually been found.\n",
    "* More generally: handling exceptions.\n",
    "* Symmetry breaking: in problems 3 and 4 we can achieve about $k$-fold speedup by setting $s_{11}=1$, i.e. declaring that the first point is covered by circle number $1$.\n",
    "* For large instances one should perhaps terminate the solver when a solution of sufficiently high quality has been located.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\"><img alt=\"Creative Commons License\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by/4.0/80x15.png\" /></a><br />This work is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\">Creative Commons Attribution 4.0 International License</a>. The **MOSEK** logo and name are trademarks of <a href=\"http://mosek.com\">Mosek ApS</a>. The code is provided as-is. Compatibility with future release of **MOSEK** or the `Fusion API` are not guaranteed. For more information contact our [support](mailto:support@mosek.com). \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}