{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mixed-Integer Linear Programming\n",
    "\n",
    "## Ice Cream\n",
    "\n",
    "This example was originally contributed by Joshua Lerman.\n",
    "\n",
    "An ice cream stand sells cones and popsicles. It wants to maximize its profit, but is subject to a budget.\n",
    "\n",
    "We can write this problem as a linear program:\n",
    "\n",
    "> **max** cone $\\cdot$ cone_margin + popsicle $\\cdot$ popsicle margin\n",
    "\n",
    "> *subject to*\n",
    "\n",
    "> cone $\\cdot$ cone_cost + popsicle $\\cdot$ popsicle_cost $\\le$ budget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cone_selling_price = 7.\n",
    "cone_production_cost = 3.\n",
    "popsicle_selling_price = 2.\n",
    "popsicle_production_cost = 1.\n",
    "starting_budget = 100."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This problem can be written as a cobra.Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'cone': 33.333333333333336, 'popsicle': 0.0}"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from cobra import Model, Metabolite, Reaction\n",
    "\n",
    "cone = Reaction(\"cone\")\n",
    "popsicle = Reaction(\"popsicle\")\n",
    "\n",
    "# constrainted to a budget\n",
    "budget = Metabolite(\"budget\")\n",
    "budget._constraint_sense = \"L\"\n",
    "budget._bound = starting_budget\n",
    "cone.add_metabolites({budget: cone_production_cost})\n",
    "popsicle.add_metabolites({budget: popsicle_production_cost})\n",
    "\n",
    "# objective coefficient is the profit to be made from each unit\n",
    "cone.objective_coefficient = \\\n",
    "    cone_selling_price - cone_production_cost\n",
    "popsicle.objective_coefficient = \\\n",
    "    popsicle_selling_price - popsicle_production_cost\n",
    "\n",
    "m = Model(\"lerman_ice_cream_co\")\n",
    "m.add_reactions((cone, popsicle))\n",
    "\n",
    "m.optimize().x_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In reality, cones and popsicles can only be sold in integer amounts. We can use the variable kind attribute of a cobra.Reaction to enforce this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'cone': 33.0, 'popsicle': 1.0}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cone.variable_kind = \"integer\"\n",
    "popsicle.variable_kind = \"integer\"\n",
    "m.optimize().x_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the model makes both popsicles and cones."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Restaurant Order\n",
    "\n",
    "To tackle the less immediately obvious problem from the following [XKCD comic](http://xkcd.com/287/):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"http://imgs.xkcd.com/comics/np_complete.png\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(url=r\"http://imgs.xkcd.com/comics/np_complete.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want a solution satisfying the following constraints:\n",
    "\n",
    "$\\left(\\begin{matrix}2.15&2.75&3.35&3.55&4.20&5.80\\end{matrix}\\right) \\cdot \\vec v = 15.05$\n",
    "\n",
    "$\\vec v_i \\ge 0$\n",
    "\n",
    "$\\vec v_i \\in \\mathbb{Z}$\n",
    "\n",
    "This problem can be written as a COBRA model as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'french_fries': 0.0,\n",
       " 'hot_wings': 2.0,\n",
       " 'mixed_fruit': 1.0,\n",
       " 'mozarella_sticks': 0.0,\n",
       " 'sampler_plate': 1.0,\n",
       " 'side_salad': 0.0}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "total_cost = Metabolite(\"constraint\")\n",
    "total_cost._bound = 15.05\n",
    "\n",
    "costs = {\"mixed_fruit\": 2.15, \"french_fries\": 2.75,\n",
    "         \"side_salad\": 3.35, \"hot_wings\": 3.55,\n",
    "         \"mozarella_sticks\": 4.20, \"sampler_plate\": 5.80}\n",
    "\n",
    "m = Model(\"appetizers\")\n",
    "\n",
    "for item, cost in costs.items():\n",
    "    r = Reaction(item)\n",
    "    r.add_metabolites({total_cost: cost})\n",
    "    r.variable_kind = \"integer\"\n",
    "    m.add_reaction(r)\n",
    "\n",
    "# To add to the problem, suppose we want to\n",
    "# eat as little mixed fruit as possible.\n",
    "m.reactions.mixed_fruit.objective_coefficient = 1\n",
    "    \n",
    "m.optimize(objective_sense=\"minimize\").x_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is another solution to this problem, which would have been obtained if we had maximized for mixed fruit instead of minimizing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'french_fries': 0.0,\n",
       " 'hot_wings': 0.0,\n",
       " 'mixed_fruit': 7.0,\n",
       " 'mozarella_sticks': 0.0,\n",
       " 'sampler_plate': 0.0,\n",
       " 'side_salad': 0.0}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.optimize(objective_sense=\"maximize\").x_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boolean Indicators\n",
    "\n",
    "To give a COBRA-related example, we can create boolean variables as integers, which can serve as indicators for a reaction being active in a model. For a reaction flux $v$ with lower bound -1000 and upper bound 1000, we can create a binary variable $b$ with the following constraints:\n",
    "\n",
    "$b \\in \\{0, 1\\}$\n",
    "\n",
    "$-1000 \\cdot b \\le v \\le 1000 \\cdot b$\n",
    "\n",
    "To introduce the above constraints into a cobra model, we can rewrite them as follows\n",
    "\n",
    "$v \\le b \\cdot 1000 \\Rightarrow v- 1000\\cdot b \\le 0$\n",
    "\n",
    "$-1000 \\cdot b \\le v \\Rightarrow v + 1000\\cdot b \\ge 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import cobra.test\n",
    "model = cobra.test.create_test_model(\"textbook\")\n",
    "\n",
    "# an indicator for pgi\n",
    "pgi = model.reactions.get_by_id(\"PGI\")\n",
    "# make a boolean variable\n",
    "pgi_indicator = Reaction(\"indicator_PGI\")\n",
    "pgi_indicator.lower_bound = 0\n",
    "pgi_indicator.upper_bound = 1\n",
    "pgi_indicator.variable_kind = \"integer\"\n",
    "# create constraint for v - 1000 b <= 0\n",
    "pgi_plus = Metabolite(\"PGI_plus\")\n",
    "pgi_plus._constraint_sense = \"L\"\n",
    "# create constraint for v + 1000 b >= 0\n",
    "pgi_minus = Metabolite(\"PGI_minus\")\n",
    "pgi_minus._constraint_sense = \"G\"\n",
    "\n",
    "pgi_indicator.add_metabolites({pgi_plus: -1000,\n",
    "                               pgi_minus: 1000})\n",
    "pgi.add_metabolites({pgi_plus: 1, pgi_minus: 1})\n",
    "model.add_reaction(pgi_indicator)\n",
    "\n",
    "\n",
    "# an indicator for zwf\n",
    "zwf = model.reactions.get_by_id(\"G6PDH2r\")\n",
    "zwf_indicator = Reaction(\"indicator_ZWF\")\n",
    "zwf_indicator.lower_bound = 0\n",
    "zwf_indicator.upper_bound = 1\n",
    "zwf_indicator.variable_kind = \"integer\"\n",
    "# create constraint for v - 1000 b <= 0\n",
    "zwf_plus = Metabolite(\"ZWF_plus\")\n",
    "zwf_plus._constraint_sense = \"L\"\n",
    "# create constraint for v + 1000 b >= 0\n",
    "zwf_minus = Metabolite(\"ZWF_minus\")\n",
    "zwf_minus._constraint_sense = \"G\"\n",
    "\n",
    "zwf_indicator.add_metabolites({zwf_plus: -1000,\n",
    "                               zwf_minus: 1000})\n",
    "zwf.add_metabolites({zwf_plus: 1, zwf_minus: 1})\n",
    "\n",
    "# add the indicator reactions to the model\n",
    "model.add_reaction(zwf_indicator)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a model with both these reactions active, the indicators will also be active"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PGI indicator = 1\n",
      "ZWF indicator = 1\n",
      "PGI flux = 4.86\n",
      "ZWF flux = 4.96\n"
     ]
    }
   ],
   "source": [
    "solution = model.optimize()\n",
    "print(\"PGI indicator = %d\" % solution.x_dict[\"indicator_PGI\"])\n",
    "print(\"ZWF indicator = %d\" % solution.x_dict[\"indicator_ZWF\"])\n",
    "print(\"PGI flux = %.2f\" % solution.x_dict[\"PGI\"])\n",
    "print(\"ZWF flux = %.2f\" % solution.x_dict[\"G6PDH2r\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because these boolean indicators are in the model, additional constraints can be applied on them. For example, we can prevent both reactions from being active at the same time by adding the following constraint:\n",
    "\n",
    "$b_\\text{pgi} + b_\\text{zwf} = 1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PGI indicator = 1\n",
      "ZWF indicator = 0\n",
      "PGI flux = 9.82\n",
      "ZWF flux = 0.00\n"
     ]
    }
   ],
   "source": [
    "or_constraint = Metabolite(\"or\")\n",
    "or_constraint._bound = 1\n",
    "zwf_indicator.add_metabolites({or_constraint: 1})\n",
    "pgi_indicator.add_metabolites({or_constraint: 1})\n",
    "\n",
    "solution = model.optimize()\n",
    "print(\"PGI indicator = %d\" % solution.x_dict[\"indicator_PGI\"])\n",
    "print(\"ZWF indicator = %d\" % solution.x_dict[\"indicator_ZWF\"])\n",
    "print(\"PGI flux = %.2f\" % solution.x_dict[\"PGI\"])\n",
    "print(\"ZWF flux = %.2f\" % solution.x_dict[\"G6PDH2r\"])"
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}