{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook we show how to define an integer problem using the Optimizer API and visualize the solution. We also mention how in this case an infeasibility certificate for the linear relaxation can be given a clear combinatorial interpretation.\n", "\n", "# Exact planar covering" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "from mosek import *\n", "import numpy as np\n", "from itertools import product\n", "import matplotlib.pyplot as plt\n", "import matplotlib.patches as patches" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem formulation\n", "\n", "In the *exact planar covering problem* we are given an $n\\times m$ rectangle (possibly with holes) and a collection of shapes (bricks). We ask if the rectangle can be tightly covered (without overlaps) with the given shapes. For example, can a $21\\times 21$ square be divided into $1\\times 8$ and $1\\times 9$ rectangles (allowing rotations)? Variants of the problem involve limited or unlimited number of bricks, maximizing the covered area, counting the coverings, etc. We assume that the shapes are built from unit squares and only consider grid-aligned coverings. See for instance the articles on [Polyominos](https://en.wikipedia.org/wiki/Polyomino) and [Tetrominos](https://en.wikipedia.org/wiki/Tetromino).\n", "\n", "A shape is defined as a list of unit squares, or rather offsets with respect to one fixed square, for example:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# The shape of a rectangle \n", "def shape_rectangle(a, b):\n", " return list(product(range(0, a), range(0, b)))\n", "\n", "# Shapes of a subset of Tetris blocks\n", "shapes_tetris = [\n", " [(0,0), (0,1), (0,2), (-1,1)],\n", " [(0,0), (0,1), (0,2), (1,1)],\n", " [(0,0), (0,1), (-1,1), (-2,1)],\n", " [(0,0), (0,1), (1,1), (1,2)],\n", " [(0,0), (0,1), (-1,1), (-1,2)],\n", " [(0,0), (1,0), (1,1), (1,2)]\n", "]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When a shape is actually placed in the rectangle, we say it is *anchored*. Not all positions are suitable for anchoring a shape - it may not stick out of the rectangle." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Anchor a shape at a given point p,q\n", "# return the list of coordinates it occupies\n", "# or return None if it would go outside the board or cover a forbidden spot (from noncov)\n", "def anchorShape(shp, n, m, p, q, noncov=[]):\n", " pts = [(p + x, q + y) for x,y in shp]\n", " if all(0<= x and x 0:\n", " task.appendcons(t)\n", " task.putconboundslice(numcon, numcon + t, [boundkey.up] * t, [rep] * t, [rep] * t)\n", " for k in range(t):\n", " task.putarow(numcon + k, [ encodeBrick(n,m,t,p,q,k) for p,q in product(range(n), range(m)) ], [1.0] * (n*m))\n", "\n", " # Optimize and get the results back\n", " if logging:\n", " task.set_Stream(streamtype.log, logging)\n", " task.putdouparam(dparam.mio_max_time, timelimit)\n", " task.optimize()\n", "\n", " prosta = task.getprosta(soltype.itg)\n", " if prosta == prosta.prim_infeas:\n", " print(\"No covering\\nLooking for infeasibility certificate for the relaxation\")\n", " attemptCertificate(n, m, noncov, task)\n", " else:\n", " xx = np.zeros(numvar, dtype=float)\n", " task.getxx(soltype.itg, xx)\n", " sol = [(p,q,k) for p,q,k in product(range(n), range(m), range(t)) if xx[encodeBrick(n,m,t,p,q,k)] > 0.8]\n", " display(n, m, sol, T, ['blue', 'yellow', 'green', 'red', 'violet', 'orange'])\n", " if not exact:\n", " print(\"Covered area {0}, best bound found {1}, total board area {2}\".format(\n", " int(task.getprimalobj(soltype.itg)), \n", " int(task.getdouinf(dinfitem.mio_obj_bound)),\n", " n*m-len(noncov)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The code above models and solves the problem. It is easier to add exactly $nm$ linear constraints even if some fields are excluded. In this case the corresponding constraint is fixed to zero (this follows anyway from the variable bounds in such case). Plotting the result is done with the function shown next." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Plot a solution\n", "def display(n, m, sol, T, col):\n", " fig,ax = plt.subplots(1)\n", " # Plot all small squares for each brick\n", " for p,q,k in sol:\n", " for x,y in anchorShape(T[k], n, m, p, q):\n", " ax.add_patch(patches.Rectangle((x,y), 1, 1, linewidth=0, facecolor=col[k]))\n", " # Plot grid\n", " xs, ys = np.linspace(0, n, n+1), np.linspace(0, m, m+1)\n", " for x in xs: plt.plot([x, x], [ys[0], ys[-1]], color='black')\n", " for y in ys: plt.plot([xs[0], xs[-1]], [y, y], color='black') \n", " ax.axis([0, n, 0, m])\n", " ax.axis('off')\n", " ax.set_aspect('equal')\n", " plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Examples\n", "\n", "From the introduction: a covering of the $21\\times 21$ square with $1\\times 8$ and $1\\times 9$ rectangles." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQgAAAD8CAYAAACLgjpEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAB0FJREFUeJzt3DFu2+gWBeDrhwHsTlOmzBK4BG3lrextxdqBlzDppjS7pPIUJhA8Zo75UyFFiv6+hkUu7jUE50ARTvTw9vZWAP/mP1v/AMB+CQggEhBAJCCASEAAkYAAIgEBRAICiAQEEP2x9Q/wbx6eHt7qy8TQ38Nzqblh9vSjqpsYexmeS82tsfOz3l5j5xFvv769PTSs3GdA1Jeq+u/EzP+G51Jzw2z3rep5Yuw8PJeaW2PnZ729xs6j3m7hnxhAJCCASEAAkYAAIgEBRAICiAQEEAkIIHrY43dSalK6vaedR7x9303KDfV1qkvLS/zY16UhnE5PdazfrL3fXmPnUW832GdAbFi1rm9dNZVVv1yabndfq56n1lXV+fz+nJpdeu5ot9fYedTbLXwGAUQCAogEBBAJCCASEEAkIIBIQACRJuV49sepmpomj33T7dNTVddQcnkZyitTs0vPHe32GjuPePv1VZNyF/r+VJdLew1uenbG3GNfl78aTn9/f0zOLj03c+fpqWEfi9pnQByoSdm2b9hZ1Xa7da7lZ6xa/rWc+5o37tRKXXZnC59BAJGAACIBAUQCAogEBBAJCCASEECkSTmeXbhJ2bZv2FnVdrt1ruVnrFr+tZz7mjfu1EpdbqcmJYezaSu1aW6YbWmwbtxKbbXPgNCkXGbuYE3KTV/LprlhtvV3o2q717KRzyCASEAAkYAAIgEBRAICiAQEEAkIIBIQQKRqPZ5Vtb7d3NydW76WTXPDbOvvRtVmr+Xb97aqtXcQQKRqPZ5Vtb7d3NydqtbL7mzgHQQQCQggEhBAJCCASEAAkYAAIgEBRJqU41lNytvNzd2pSbnYTk1K4LdpUo5nNSlvNzd3pyblsjsbeAcBRAICiAQEEAkIIBIQQCQggEhAAJEm5XhWk/J2c3N3alIutlOTEvhtmpTjWU3K283N3alJuezOBt5BAJGAACIBAUQCAogEBBAJCCASEEAkIIBI1Xo8q2p9u7m5O1WtF9upag38NlXr8ayq9e3m5u5UtV52ZwPvIIBIQACRgAAiAQFEAgKIBAQQCQgg0qQcz2pS3m5u7k5NysV2alICv02TcjyrSXm7ubk7NSmX3dnAOwggEhBAJCCASEAAkYAAIgEBRAICiAQEEKlaj2dVrW83N3enqvViO1ur1vtsUm7psX9vwn3k74//GI5inwGxZdV6yZ1b14NVrW88N8yqWgOfgYAAIgEBRAICiAQEEAkIIBIQQKRJuebOrdt/mpQ3nhtmNSmpH6ehuPORl4k/h33bZ0DcQ5Oy9cttqxrm5szOmNOkvPHcMKtJCXwGAgKIBAQQCQggEhBAJCCASEAAkSbltTtbv7uyqmFuzqwm5TRNyqk5TUrq9FTVfZ2ee3l9f07NLj03e+dTX1038X2hVfUy/H2emp0z13+v6e8qrTrc95XuMyA0KReZ67pLPTecPg8rp2aXnruXnedz1eWvOlYrtZHPIIBIQACRgAAiAQFEAgKIBAQQCQggEhBApGp97c47qFqfTn11Dad/NgpvO3cvO382Kadv30ttvbVqLSCu3XkHAfE5b6+x83j/r+W+/y+GqvUN5452e42d5+N9AXAjn0EAkYAAIgEBRAICiAQEEAkIIBIQQKQode1ORamd3l5j5+ctSnkHAUSalNfu1KTc6e01dp41KQHGBAQQCQggEhBAJCCASEAAkYAAIk3Ka3dqUu709ho7NSkBfqFJee1OTcqd3l5j51mTEmBMQACRgAAiAQFEAgKIBAQQCQggEhBApGp97U5V653eXmOnqjXAL1Str92par3T22vsPKtaA4wJCCASEEAkIIBIQACRgAAiAQFEmpTX7tSk3OntNXZqUgL8QpPy2p2alDu9vcbOsyYlwJiAACIBAUQCAogEBBAJCCASEECkSXntTk3Knd5eY6cmJcAvNCmv3alJudPba+w8a1ICjAkIIBIQQCQggEhAAJGAACIBAUQCAohUra/dqWq909tr7Py8VetdNilPT1Xd149nXl7fn0vNzd751FfXXT6eG37/pubmzM6Z6/vT5F34yC4Douuqnp8/njmf359Lza2xc+vbl0tLHbxK1bphTtUa4P8JCCASEEAkIIBIQACRgAAiAQFEu2xS/vnnw1s3UW77WRhaZm6NnVvffi9KaVIuMvdJm5TeQQCRJuWKO7e+rUm51M6zJiXAmIAAIgEBRAICiAQEEAkIIBIQQKRJueLOrW9rUmpSprm7/k5KlnGqvrpq+D7M4Tk1u/Tcvex8qar6UdV9mzxdL1XVP07P3YtdBoQm5UJzFz3KLW5fvpQmJXB8AgKIBAQQCQggEhBAJCCASEAAkYAAIlXrFXdufrtXtN7idv9Yh6la7zIggH3wTwwgEhBAJCCASEAAkYAAIgEBRAICiAQEEAkIIBIQQCQggEhAAJGAACIBAUQCAogEBBAJCCASEEAkIIBIQACRgAAiAQFEAgKI/gEfVdTxnB/kAQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "n, m = 21, 21\n", "T = [shape_rectangle(1,8), shape_rectangle(8,1), shape_rectangle(1,9), shape_rectangle(9,1)]\n", "t = len(T)\n", "model(n, m, t, T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another rectangle and set of blocks. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAANYAAAD8CAYAAAAL1Fp+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAABg9JREFUeJzt3TFuG0cYBeBRECDq6NKljsAj6Co5ma8i3UBHcDqXYsdUTEEBQQBS/h3u48zsfl/D4gkrSvDDmnwY8eF0OjVgWb/1fgKwRooFAYoFAYoFAYoFAYoFAYoFAYoFAYoFAb/f85t9+fJw2u8vZ29v50f59fxwbK19vZy31lr78fF47Wt+tLZ7HPtnHD1/fz89XE7/667F2u9be3m5nD0/nx/l1/PX7621Py/nrbXWvn08Xvuab63tn8b+GWfIK/xXEAIUCwIUCwIUCwIUCwIe7nnQ0dvt3m6fPa++3e6OBQF2rIlyO9YYeYU7FgQoFgQoFgQoFgQoFgTYsSbK7Vj9czsWdGTHmii3Y42RV7hjQYBiQYBiQYBiQYBiQYBiQYCBeKLcQNw/NxBDRwbiiXID8Rh5hTsWBCgWBCgWBCgWBCgWBNixJsrtWP3zIT8fi/4Ox4+37S85nh/kn+cVdqyJ8iV2LPkCeYHXWBCgWBCgWBCgWBCgWBBgx5ooX2LHkt+Wn47OY0E3dqyJcjvWIHmBOxYEKBYEKBYEKBYEKBYE2LEmyu1Y/fPqjnXXYj08Ppxu+qH/3rXWrvzLax//Mtec/3FQrM55tVj3Pej4td22Ify1b629XPmC54/HFedfX+1YI+QFXmNBgGJBgGJBgGJBgGJBgLfbZ8q93d49dx4LOrJjzZTbscbIC9yxIECxIECxIECxIECxIMCONVNux+qe27GgIzvWTLkda4y8wB0LAhQLAhQLAhQLAhQLAuxYM+V2rO65HQs6smPNlNuxxsgL3LEgQLEgQLEgQLEgQLEgwI41U27H6p7bsaAjO9ZMuR1rjLzAHQsCFAsCFAsCFAsCFAsC7Fgz5Xas7nl1x7rv2+3wE7vH1vZPl7O39/Nj77zCjjVTvoEda//U2svL5fj5+fzYO6/wGgsCFAsCFAsCFAsCFAsC7Fgz5RvYsXaPre2v/ArePn5FPfP3d+exoBs71ky5Hau11j+vcMeCAMWCAMWCAMWCAMWCADvWTLkdq7U2x47lPBZDORxbe/1+JTyeH3rnFXasmfIN7FhT5AVeY0GAYkGAYkGAYkGAYkGAHWumfAM71uj5Kv+u4G53aPv968Xs33Fv5fnTxfj8NYW/i3f4hS2G/2+qHWuGszq98uo1Xr+3/jvQ7HmB11gQoFgQoFgQoFgQoFgQMNWONfpZnZ559RqHYxt6Jxo9r+5Y7lgQYMdaSV69hh1rgbzAHQsCFAsCFAsCFAsCFAsC7FgryavXsGPdlq/yPNbouv5NvF+4Bnl2rAXzrhvRPb6HvMxrLAhQLAhQLAhQLAhQLAiwYy2Yd92I7vE95M5jQU92rAVzO9ZG8gJ3LAhQLAhQLAhQLAhQLAiwYy2Y27HWn6/yPFbX804rOeu0e/z887Na65uv5fO7ptqx5J/kxWvYAhfIC7zGggDFggDFggDFggDFgoCpdiz5J3nxGrbA23LnsaAjO9Za8uI17FgL5AXuWBCgWBCgWBCgWBCgWBBgx1pLXryGHeu2fJXnsbidM233YcdaSz7Cc9hKXuA1FgQoFgQoFgQoFgQoFgTYsdaSj/AcNpA7jwUd2bHWko/wHLaSF7hjQYBiQYBiQYBiQYBiQYBiQYCBeC35CM9hA7mBGDoyEK8lH+E5bCUvcMeCAMWCAMWCAMWCAMWCADvWWvIRnsMGcn+ws4Pd4/nzpy55ez8/pvLqNQ4r+qOYI7NjLZj3/FC36jWm+GC30fMCr7EgQLEgQLEgQLEgQLEgwI61YN7zQ92q1xj9g91Gz53Hgo7sWAvmdqyN5AXuWBCgWBCgWBCgWBCgWBBgx1owt2OtP1/leaye552qObQ22Y7VcycaPa9ew461QF7gNRYEKBYEKBYEKBYEKBYETLVj9dyJRs+r17Bj3ZY7jwUd2bFWklevYcdaIC9wx4IAxYIAxYIAxYIAxYIAO9ZK8uo17Fi35dUd667Fgq3wX0EIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwIUCwI+AdD0CTkQruVQQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "n, m = 22, 27\n", "T = [shape_rectangle(8,2), shape_rectangle(5,2), shape_rectangle(1,7)]\n", "t = len(T)\n", "model(n, m, t, T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we consider an example with a small subset of Tetris blocks we defined earlier. These blocks are:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAB+CAYAAAAEA/ugAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAABRtJREFUeJzt2DFOI2kexuG3FiRWHVCETcRewFLdAF+AQ0yfo8M+R/chuIAdkJfkYFOIIKQIdgZpWrUBeEa7UrdYCX81f+3zJJVY9eLPrp+Qu3meA0Adf1v6DwDgfyPcAMUIN0Axwg1QjHADFCPcAMUIN0Axwg1QjHADFHN8iJuenXXzMBzizj82jsn0W5KPbXfz8Hptvbvk9kPSHyXDRePdJONdMn8/zep81XR3d7/L0/PTlGRsOvxi/zS13h76pG/8KCd5eaNTssR5D0nf/3nkrYyZ58fura8+SLiHIdlsDnHnH1uvk+1tkl/a7ubb67X17pLb35LhQ7L53Hg3yfpL8vuvq1x/um66e/X1Kje3N+M8z+umw0m6rtskSevtrus2Q3K5aTn6ap1kmzQ/75ezHi6TTcvZvLzjt/NTCUAxwg1QjHADFCPcAMUIN0Axwg1QjHADFCPcAMUIN0Axwg1QjHADFCPcAMUIN0Axwg1QjHADFCPcAMUIN0Axwg1QjHADFCPcAMUIN0Axwg1QjHADFCPcAMUIN0Axwg1QjHADFCPcAMV08zy/+03Pzrp5GN79tj81jsn0W5KPbXfz8Hptvbvk9kPSHyXDRePdJONdMn8/zep81XR3d7/L0/PTlGRsOvxi/zS13h76pG/8KCd5eaNTssR5D0nf/3nkrYyZ58fura8+PuSf0lr/92T4R9vN8TGZpky5W+iBPkm/wO7/pdOT0351vrpsvbu736V7fsqQNN1e4gv9n/o+GRqf9/Lv+i0OEu5hSDabQ9z5x9brl+sSu9ttxnme122Xk67rNvmYy/zSePhbMnxINp8b7yZZf0l+/3WV60/XTXevvl4lSfPd/fbx7U02jXfXr9fWu/vtbYYF1tev16V238Zv3ADFCDdAMcINUIxwAxQj3ADFCDdAMcINUIxwAxQj3ADFCDdAMcINUIxwAxQj3ADFCDdAMcINUIxwAxQj3ADFCDdAMcINUIxwAxQj3ADFCDdAMcINUIxwAxQj3ADFCDdAMcINUEw3z/O73/TsrJuH4d1v+1Pj+HJdYneaMiUZ2y4nSYacpM/HxqsPSX+UDBeNd5OMd8n8/TSr81XT3d39Lkma7+63u+enNP5q//GFbr27357SL7C+1LseM8+P3VtffXzIP6W1acq03TYP6BLf6z/0R8nwoe3meNR27789PT9NN7c3zT/n05PTvvHm4qZk2i71T0mmJNsFnuf+L/85HyTcw5BsNoe484+t18l2m3Ge53XL3a7rNknSene/PVzkcvO57e76y8u19e5+e/vPZT7n1fnq8vrTdcvZJMnV16sc395k03h3nWSb9medLPdcvewOl1nktN/Ob9wAxQg3QDHCDVCMcAMUI9wAxQg3QDHCDVCMcAMUI9wAxQg3QDHCDVCMcAMUI9wAxQg3QDHCDVCMcAMUI9wAxQg3QDHCDVCMcAMUI9wAxQg3QDHCDVCMcAMUI9wAxQg3QDHCDVCMcAMU083z/O43PTvr5mF499v+1Dgm05Qpydh2Oft32no3SYb+Q/rhou3oePc63nh3vz39a5nP+fTktF+drxrPJrv7XbrnpzR+pDImmbLIWSfLPVdD0vdZ4LTn+bF766sPEm4ADsdPJQDFCDdAMcINUIxwAxQj3ADFCDdAMcINUIxwAxQj3ADFCDdAMcINUIxwAxQj3ADFCDdAMcINUIxwAxQj3ADFCDdAMcINUIxwAxQj3ADFCDdAMcINUMy/AZsdEUPrCMDlAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Covered area 24, best bound found 24, total board area 33\n" ] } ], "source": [ "model(11, 3, len(shapes_tetris), shapes_tetris, noncov=[], exact=False, rep=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now ask MOSEK for the maximal-area covering of a sample rectangle with holes by our Tetris blocks. You may want to enable logging to track the mixed-integer optimizer's progress." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAALIAAAD8CAYAAADT2P50AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAACixJREFUeJzt3UFOJEcaxfGXAyNoFhWsRrDCFygpb0BdoA9hn8NLn8N9iL4Auah9SSzYwsbd8sZkScagBuUssmBaGsnEFzjD2c//34ZNfoqoroeVEs8RzTAMAr51//q7NwD8FQgyLBBkWCDIsECQYYEgwwJBhgWCDAsEGRb2ay7WHDaDTgIDn3c/p5z5LC200PJ0mb3E5adLbbWddl8lMzXWeMNM2pPas/yRzY10+/vQ5DxbNcg6kfR94PkPu59TznyQls1SH3/4mL3E+5/faz2sZ/lZJl/jDTPtkXTxY/7I6qf8Z3m1gAWCDAsEGRYIMiwQZFggyLBAkGGBIMMCQYYFggwLTc3/i5quRWBfJTN0LepIh1L7Xf7zm1up79XrRpvAMq0OlKJ7i1pooWUTCL8u1ext1R7lr7HZk/qn2L7Sg9TeBNaQpCNNvq+pVQ1y20oXF/nPr1ZS12kzDMMqd6Zpmgud6Hzq0pCk8Mz+u3W4NNPdKVQaam+ki/wltJKks3iZJ7Sv571RGgL+HEGGBYIMCwQZFggyLBBkWCDIsECQYYEgwwJBhoWqpaHj42Zo2/znN5td10IFXYuJS0OSwjPN3jZcmumfFCoNpQcp8E/8v67FlPt63ptLaahESkptq/Pc5zcbSX1+cWb8Irfaf7fO3lOzJw1Pi+zn3yJSAor8tr9V2osXjaY0+9KQVDDT5RdnVlJRaebxj3jRqKQ0pKvgZ1Gd0pBUbyYH78iwQJBhgSDDAkGGBYIMCwQZFggyLBBkWCDIsECQYWH2pSFp/NN2aKbPL86UlmaGp3jRqKQ0pLvgZ1Gd0pBUZ8amNKReUvd3b+L/NXvxolEVR5ICYVHgVKKv9Xfqu6tgK1EpdVfhX7Ms8y8NBQpAUrw4s5LmXZqJloYqfZbuquAEKLXnZd/m63hHhgWCDAsEGRYIMiwQZFggyLBAkGGBIMMCQYYFggwL8y8NBQpAUrw4M/fSTLg0VOmz9HcFJ0Appei3OQy3JqWhEpHiTGlp5ml3s1Gu3XVe4ZkDqYscmfVvSf8JrPHLeCxZSQFo7E7kmvYcJM/SUKA4U1yaKbieS5p45sN4j2GVK+CKC0AlM6/jHRkWCDIsEGRYIMiwQJBhgSDDAkGGBYIMCwQZFjy7FoG+QXHXoOBWI2nimc9SOowfaFN0c1ZBb+JlNDDj07VIKvv3yu0b/CL192W9CRdlN2f1agMn54z/TUoFu8sz/66Fpp1ZraTuWvPrTZTMFHYtpDodmE6t6FoAf4IgwwJBhgWCDAsEGRYIMiwQZFggyLBAkGGBIMPC/EtDKrjVKTCz2Yxdi9kVgEpmCktD0rQ3Z0lfdy3+oaWhvk/qutiHT6kPrZEOx45C9gq3Gr/IwOEuG0n9QWhbkqSFFlo2edegXepS0ja+SKWbs5LiRaNcsy8NjSEODGmltu1CpSGpUmnmROHS0LJZ6uMPH7Mef//ze+2frWd3c9ZbZ3LwjgwLBBkWCDIsEGRYIMiwQJBhgSDDAkGGBYIMCwQZFmZfGur7eNEkpT5UGpIqlWYOFC4NLbTQ8jSza/HpUs3htspnkerM3A7D/EpDJQWgKiqVZiIFIOm5BBQz3C/0eJO/xnB/qaakaDQzldtvpSfNxGbCpaFKpZnH0/wCkDSWgCSFSkOR559n9q/XlIaAOSDIsECQYYEgwwJBhgWCDAsEGRYIMiwQZFggyLBQtTTUNMdDjapJuDRUqTQzHOQXgKSxBCQpVBqKPP880zxsKQ1FjAELnDTz0n6boZJr0+4n2ssb9VLfBe/ZS4rfNdYr7W52ypW/JU4aWqmsNFTwWR5v5lkaWl+vN8MwrHJnmqa5aKXzi+xVuJ4MyEKQYYEgwwJBhgWCDAsEGRYIMiwQZFggyLDAAS2lXYsUP9RluJ9n12L7sO1V8CdqbnUKKOlnxBdRWQMmqDncav9snf/8b2P4I7YP2359vQ6FcnGwSMvT5XnuwOWnS/UPzWS9iRKz71pI087UWOMtM5F+RmlvYnm6PI/3M/Y1pyNaeEeGBYIMCwQZFggyLBBkWCDIsECQYYEgwwJBhgWCDAuzLw1JBTcuBWZqrPGWmUjRqLQAtOtaZA+M6zSqcUSLTWmo6Jaiw3neUtTfS911YOBekrZaD/lFIx0o6UTZBSB9DuznK4uDQcvTx+znLz8Nu/BPY/aloZJDTfbP1rMsDXXXkr7Pn9GH3c/cmejzu5llM+3BMc8z5UWj1/GODAsEGRYIMiwQZFggyLBAkGGBIMMCQYYFggwLBBkWZl8aKjmdpznczrI01N9LOsmfeelB5M5En9/NLDTtCUjPMyVFI5vSUIlI0ahqyehBvW6CR1Mlpfa7vIc3t7tflqBmb6v9d4ETkPak4Sl2AtLULEtDUuwmpEjJ6HlfUskNVYrfntTqPFKAKikmtUfSxY/5I6ufpMc/4t8LpSHgFQQZFggyLBBkWCDIsECQYYEgwwJBhgWCDAsEGRYsS0NS7EqvSMnoeV9SQWmoV/wasKQUKUCVFJPSntSe5Y9sbsauBaWhgBpXepXoe/VdFwullJLU5p8CtLudLiIdSrklI2ksGulLaAlJ0vah2XUnck13ypD0DZSGpOmv9CopDZUUgMYQBxbSSm3bTX9q0q/x0lB31YrryYC/GEGGBYIMCwQZFggyLBBkWCDIsECQYYEgwwJBhoXZl4ak6a/0KikNlRSAdl2LwMhGKfXTn5r0JV4a6u+SuJ4soKSck9I25RaNmt9qFoAq+CLp1+DzBZJ6teqynx/rT6lssQyzLw1VOZ2nUgFoFJsJl4YKCkBSwcxV/NN3Ki0avY53ZFggyLBAkGGBIMMCQYYFggwLBBkWCDIsEGRYmH3XosqhJpV6Ey+jgZlw16KgNyEVzNzFP/34J2qDrsVsDzUZg5+9xhiYXu1ZoGvwUrQJ6qVApUE6ii+hO0lXBXNBJf2MXLVLQxxqEv0s6rInVpJ0Vqc3Ic3peBbekWGCIMMCQYYFggwLBBkWCDIsEGRYIMiwQJBhgSDDQtXSUNM0t+JQk/whbXb9hNynJR3VKQBJdWZuh2F+paG5HmpSdAjMUfy0kSlLM/90lYNc5yagaGmo6BCYM53PrZyzkigNAd8yggwLBBkWCDIsEGRYIMiwQJBhgSDDAkGGBYIMC5VLQ8fDHE/nKTrN6EhpbuUcSkPV9L3UFbTfwqtESkDPocw/aegmuqM3OJKUG8zCffVS30V/kRUvTZWsk/ug5UlDUpe9TnEBSJXKOYESUOm+uquC0pR0fpG/zO5Wp5LvPw/vyLBAkGGBIMMCQYYFggwLBBkWCDIsEGRYIMiwQJBhwfKkobHTEe5a5K9Q8UqvSAmodF/9XUFpSkrx68ni6wzDcJzzYNUgA1Ph1QIWCDIsEGRYIMiwQJBhgSDDAkGGBYIMCwQZFggyLBBkWCDIsECQYYEgwwJBhgWCDAsEGRYIMiwQZFggyLBAkGGBIMMCQYaF/wKopSmOPqSTWgAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Covered area 172, best bound found 172, total board area 177\n" ] } ], "source": [ "# Define a stream printer to grab output from MOSEK\n", "def streamprinter(text):\n", " sys.stdout.write(text)\n", " sys.stdout.flush()\n", "n, m = 11, 17\n", "T = shapes_tetris\n", "t = len(T)\n", "noncov = [(0,0), (1,3), (9,13), (8,8), (7,7), (5,5), (4,4), (3,3), (3,1), (8,12)]\n", "model(n, m, t, T, noncov, exact = False, rep = 0, timelimit = 20.0, logging = None)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Combinatorial interpretation of infeasibility\n", "\n", "In some cases the integer problem is declared as infeasible already because its linear relaxation is infeasible. This case deserves some attention in our example. The linear relaxation of the exact covering problem has the form\n", "\n", "$$\n", "\\begin{array}{l}\n", "Ax = b, \\\\\n", "x\\geq 0,\n", "\\end{array}\n", "$$\n", "\n", "where $A$ is the adjacency matrix discussed previously and $b$ is the characteristic vector of available positions on the board.\n", "\n", "Standard duality for linear programing and Farkas lemma imply that a [certificate of primal infeasibility](http://docs.mosek.com/modeling-cookbook/linear.html#infeasibility-in-linear-optimization) is a vector $y$ satisfying\n", "\n", "$$\n", "\\begin{array}{l}\n", "A^Ty \\geq 0, \\\\\n", "b^Ty < 0.\n", "\\end{array}\n", "$$\n", "\n", "It means that an infeasibility certificate is an assignment of a real number to every position on the board so that:\n", "\n", "* every possible placement of a single brick covers positions with non-negative sum\n", "* the sum of all numbers on the board is negative\n", "\n", "It is combinatorially obvious that the existence of such an assignment implies an exact covering does not exist (unfortunately not vice-versa since the covering problem is NP-hard in general; in other words integer infeasibility does not imply continuous infeasibility).\n", "\n", "It is very easy to compute a relaxed infeasibility certificate, if one exists. All we need to do is reoptimize the task changing all integer variables to continuous ones and read off the certificate stored in the dual variables corresponding to constraints." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# Check if the linear relaxation is infeasible\n", "# And if so, print the infeasibility certificate\n", "# as a labeling of the rectangle grid.\n", "def attemptCertificate(n, m, noncov, task):\n", " # Now we make the problem continuous\n", " task.putvartypelist(range(task.getnumvar()), [variabletype.type_cont] * task.getnumvar())\n", " task.optimize()\n", " if task.getprosta(soltype.itr) == prosta.prim_infeas:\n", " # Read the dual variables containing the certificate\n", " y = np.zeros(n * m, dtype=float)\n", " task.getyslice(soltype.itr, 0, n * m, y)\n", " for p in range(n):\n", " print(' '.join(' ' if (p,q) in noncov else '{: 3.1f}'.format(y[encodePos(n, m, p, q)])for q in range(m)))\n", " print('Certificate with sum = {0}'.format(sum(y)))\n", " else:\n", " print('No certificate found')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Certificate example\n", "\n", "Let us use MOSEK to solve the following [puzzle from cut-the-knot](https://www.cut-the-knot.org/blue/Defective12x12Square.shtml). Can the $12\\times 12$ square with three corners removed be covered using $1\\times 3$ and $3\\times 1$ tiles?\n", "\n", "We solve this problem as follows:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "No covering\n", "Looking for infeasibility certificate for the relaxation\n", " -1.0 -0.0 1.0 -1.0 -0.0 1.0 -1.0 -0.0 1.0 -1.0 \n", "-1.0 1.0 -0.0 -1.0 1.0 -0.0 -1.0 1.0 -0.0 -1.0 1.0 -0.0\n", "-0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0\n", " 1.0 -1.0 -0.0 1.0 -1.0 -0.0 1.0 -1.0 -0.0 1.0 -1.0 -0.0\n", "-1.0 1.0 -0.0 -1.0 1.0 -0.0 -1.0 1.0 -0.0 -1.0 1.0 -0.0\n", "-0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0\n", " 1.0 -1.0 -0.0 1.0 -1.0 -0.0 1.0 -1.0 -0.0 1.0 -1.0 -0.0\n", "-1.0 1.0 -0.0 -1.0 1.0 -0.0 -1.0 1.0 -0.0 -1.0 1.0 -0.0\n", "-0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 -0.0 -0.0\n", " 1.0 -1.0 -0.0 1.0 -1.0 -0.0 1.0 -1.0 0.0 1.0 -1.0 -0.0\n", "-1.0 1.0 -0.0 -1.0 1.0 -0.0 -1.0 1.0 -0.0 -1.0 1.0 -0.0\n", " -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0\n", "Certificate with sum = -1.0\n" ] } ], "source": [ "model(n = 12, m = 12,\n", " t = 2, T = [shape_rectangle(1,3), shape_rectangle(3,1)],\n", " noncov = [(0, 0), (0, 11), (11, 0)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We displayed the infeasibility certificate. Every $1\\times 3$ or $3\\times 1$ brick covers fields with nonnegative sum (in this case, in fact, exactly $0$), while the sum of all numbers is $-1$. That proves a covering does not exist." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Acknowledgement\n", "\n", "\n", "Thanks to Jaroslaw Wroblewski for inspiring problems and discussions originating from his newsletter [Trapez](http://www.math.uni.wroc.pl/~jwr/trapez/)." ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "
This work is licensed under a Creative Commons Attribution 4.0 International License. The **MOSEK** logo and name are trademarks of Mosek ApS. The code is provided as-is. Compatibility with future release of **MOSEK** or the Fusion API are not guaranteed. For more information contact our [support](mailto:support@mosek.com). " ] }, { "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 }