{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The simplex tableau"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Michel Bierlaire](https://people.epfl.ch/michel.bierlaire), EPFL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is Algorithm 16.4 in <a href=\"http://optimizationprinciplesalgorithms.com/\">Bierlaire (2015) Optimization: principles and algorithms, EPFL Press.</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, the routine that performs one iteration of the simplex using the tableau."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simplexTableau(tableau):\n",
    "    \"\"\"\n",
    "    :param tableau: the first simplex tableau\n",
    "    :type tableau: pandas dataframe\n",
    "    \n",
    "    :return: p, q, opt, bounded  where \n",
    "               - p is the column of the variable that must enter the basis, or None,\n",
    "               - q is the row of the variable that must leave the basis, or None,\n",
    "               - opt is True if the tableau is optimal (in this case, p and q are None)\n",
    "               - bounded is False if basic direction is unbounded (in this case, p and q are None)\n",
    "    :rtype: int, int, bool, bool\n",
    "    \"\"\"\n",
    "    mtab, ntab = tableau.shape\n",
    "    m = mtab - 1\n",
    "    n = ntab - 1\n",
    "\n",
    "    reducedCost = tableau[-1, : -1]\n",
    "    # Identify the negative reduced costs\n",
    "    negativeReducedCost = reducedCost < 0\n",
    "    if not negativeReducedCost.any():\n",
    "        # The tableau is optimal\n",
    "        return None, None, True, True\n",
    "\n",
    "    # In Python, True is larger than False. The next statement returns the \n",
    "    # index of a True entry in the array, that is the index of a negative reduced cost.\n",
    "    # It is the index of the variable that will enter the basis.\n",
    "    p = np.argmax(negativeReducedCost)\n",
    "\n",
    "    # Calculate the maximum step that can be done along the basic direction d[p]\n",
    "    xb = tableau[:-1, -1]\n",
    "    minusd = tableau[:-1, p]\n",
    "    steps = np.array([xb[k] / minusd[k] if minusd[k] > 0 else np.inf for k in range(m)])\n",
    "    q = np.argmin(steps)            \n",
    "    step = steps[q]\n",
    "\n",
    "    if step == np.inf:\n",
    "        # The tableau is unbounded\n",
    "        return None, None, False, False\n",
    "    else:\n",
    "        return p, q, False, True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, the routine to pivot the tableau."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pivoting(tableau, p, q):\n",
    "    \"\"\"\n",
    "    :param tableau: valid simplex tableau\n",
    "    :type tableau: numpy.array 2D\n",
    "    \n",
    "    :param p: columns of the pivot\n",
    "    :type p: int\n",
    "    \n",
    "    :param q: row of the pivot\n",
    "    :type q: int\n",
    "    \n",
    "    :return: pivoted tableau\n",
    "    :rtype: numpy.array 2D\n",
    "    \"\"\"\n",
    "    m, n = tableau.shape\n",
    "    if q >= m:\n",
    "        raise Exception(f'The row of the pivot ({q}) must be between 0 and {m - 1})')\n",
    "    if p >= n:\n",
    "        raise Exception(f'The column of the pivot ({p}) must be between 0 and {n - 1})')\n",
    "    thepivot = tableau[q][p]\n",
    "    if np.abs(thepivot) < np.finfo(float).eps:\n",
    "        raise Exception(f'The pivot is too close to zero: {thepivot}')\n",
    "    thepivotrow = tableau[q, :]\n",
    "    newtableau = np.empty(tableau.shape)\n",
    "    newtableau[q, :] = tableau[q, :] / thepivot\n",
    "    for i in set(range(m)) - {q}:\n",
    "        newtableau[i, :] = tableau[i, :] - tableau[i][p] * thepivotrow / thepivot\n",
    "    return newtableau"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we put everything together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simplexAlgorithmTableau(tableau):\n",
    "    \"\"\"\n",
    "    :param tableau: valid simplex tableau\n",
    "    :type tableau: numpy.array 2D\n",
    "\n",
    "    :return: tableau, optimal, unbounded, where tableau is the tableau from the last iteration,\n",
    "                                          optimal is True if the last tableau is optimal,\n",
    "                                          unbounded is True if the problem is unbounded.\n",
    "    :rtype: numpy.array 2D, bool, bool\n",
    "    \"\"\"\n",
    "    while True:\n",
    "        colPivot, rowPivot, optimal, bounded = simplexTableau(tableau)\n",
    "        if optimal:\n",
    "            return tableau, True, False\n",
    "        if not bounded:\n",
    "            return tableau, False, True\n",
    "        tableau = pivoting(tableau, colPivot, rowPivot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function is designed to print the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def printResults(res):\n",
    "    lastTableau, optimal, unbounded = res\n",
    "    if optimal:\n",
    "        print('Optimal tableau')\n",
    "        print(pd.DataFrame(lastTableau).to_string(index = False, header = False))\n",
    "    elif unbounded:\n",
    "        print('Unbounded problem')\n",
    "    else:\n",
    "        print('Incorrect output')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following initial tableau:     \\\\[\n",
    "    \\begin{array}{|cccccc|c|}\n",
    "      \\hline\n",
    "      1& 2& 2& 1& 0& 0& 20 \\\\\n",
    "      2& 1& 2& 0& 1& 0& 20 \\\\\n",
    "      2& 2& 1& 0& 0& 1& 20 \\\\\n",
    "      \\hline\n",
    "      -10& -12& -12& 0& 0& 0& 0 \\\\\n",
    "      \\hline\n",
    "    \\end{array}\n",
    "    \\\\]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = np.array([[1, 2, 2, 1, 0, 0, 20], \n",
    "              [2, 1, 2, 0, 1, 0, 20], \n",
    "              [2, 2, 1, 0, 0, 1, 20], \n",
    "              [-10, -12, -12, 0, 0, 0, 0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal tableau\n",
      " 0.0  0.0  1.0  0.4  0.4 -0.6    4.0\n",
      " 1.0  0.0  0.0 -0.6  0.4  0.4    4.0\n",
      " 0.0  1.0  0.0  0.4 -0.6  0.4    4.0\n",
      " 0.0  0.0  0.0  3.6  1.6  1.6  136.0\n"
     ]
    }
   ],
   "source": [
    "res = simplexAlgorithmTableau(T) \n",
    "printResults(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following initial tableau:     \\\\[\n",
    "    \\begin{array}{|rrrrrrrr|r|}\n",
    "      \\hline\n",
    "      1& 2& 3& 0& 1& 0& 0& 0& 3 \\\\\n",
    "      -1& 2& 6& 0& 0& 1& 0& 0& 2 \\\\\n",
    "      0& 4& 9& 0& 0& 0& 1 & 0& 5 \\\\\n",
    "      0& 0& 3& 1& 0& 0& 0& 1& 1 \\\\\n",
    "      \\hline\n",
    "      0& -8& -21& -1& 0& 0& 0& 0& -11 \\\\\n",
    "      \\hline\n",
    "    \\end{array}\n",
    "    \\\\]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = np.array([[1, 2, 3, 0, 1, 0, 0, 0, 3], \n",
    "              [-1, 2, 6, 0, 0, 1, 0, 0, 2], \n",
    "              [0, 4, 9, 0, 0, 0, 1 , 0, 5], \n",
    "              [0, 0, 3, 1, 0, 0, 0, 1, 1], \n",
    "              [0, -8, -21, -1, 0, 0, 0, 0, -11]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal tableau\n",
      " 1.0  0.0  0.0  0.500000  0.50 -0.50  0.0  0.500000  1.000000\n",
      " 0.0  1.0  0.0 -0.750000  0.25  0.25  0.0 -0.750000  0.500000\n",
      " 0.0  0.0  0.0  0.000000 -1.00 -1.00  1.0  0.000000  0.000000\n",
      " 0.0  0.0  1.0  0.333333  0.00  0.00  0.0  0.333333  0.333333\n",
      " 0.0  0.0  0.0  0.000000  2.00  2.00  0.0  1.000000  0.000000\n"
     ]
    }
   ],
   "source": [
    "res = simplexAlgorithmTableau(T) \n",
    "printResults(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Consider the linear optimization problem\n",
    "    \\\\[\n",
    "    \\min_{x \\in \\mathbb{R}^2} 3 x_1 - 2 x_2\n",
    "    \\\\]\n",
    "    subject to\n",
    "    \\\\[ \\begin{array}{rl}\n",
    "      x_1 &= 1, \\\\\n",
    "      x_1, x_2 & \\geq 0.\n",
    "    \\end{array} \\\\]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem is clearly unbounded, as $x_1$ is fixed by the constraints, and increasing $x_2$ decreases the objective function. We will verify that it is detected by the algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is a problem in standard form with \\\\[A = (1 \\; 0), \\; b = 1, \\; c = \\left(\\begin{array}{r} 3 \\\\ -2 \\end{array}\\right). \\\\]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We select $x_1$ to be in the basis, so that $B=1$ and $c_B=3$. The tableau is \n",
    "\\\\[\n",
    "\\begin{array}{| c| c |}\n",
    "\\hline\n",
    "B^{-1}A & B^{-1}b\\\\\n",
    "\\hline\n",
    "c^T -c_B^TB^{-1}A & -c_B^TB^{-1}b\\\\\n",
    "\\hline\n",
    "\\end{array}\n",
    "\\\\]\n",
    "that is\n",
    "\\\\[\n",
    "\\begin{array}{|cc| c |}\n",
    "\\hline\n",
    "1 & 0 & 1\\\\\n",
    "\\hline\n",
    "0 & -2 & -3\\\\\n",
    "\\hline\n",
    "\\end{array}\n",
    "\\\\]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that $x_2$ cannot be in the basis. Indeed, it would mean that $B=0$, which is not valid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now apply the algorithm on this tableau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = np.array([[1, 0, 1], \n",
    "              [0, -2, -3]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Unbounded problem\n"
     ]
    }
   ],
   "source": [
    "res = simplexAlgorithmTableau(T) \n",
    "printResults(res)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}