An example of modeling and solving integer linear constraints in ASP using a few different approaches

Integer Programming.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the programs below, we'll be reasoning about the relations between a bunch of integer variables with values taken from potentially large domains. There are several ways to do this in ASP. Rather than code for a specific strategy, we'll allow modeling [Integer Programs](https://en.wikipedia.org/wiki/Integer_programming). We'll implement a few different IP solvers later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Formulation\n",
    "\n",
    "Our goal is to model integer linear constraints like this:\n",
    "$$w_a V_a + w_b V_b + \\cdots +  w_z V_z = k$$\n",
    "\n",
    "For convenience in implementing the solvers, we'll require that the domain of all of the integer variables is given by some lower and upper bounds. Additionally, we would like to use relations other than equality. We'll allow disequality and inequality relations as well.\n",
    "\n",
    "## Relations\n",
    "These terms should only occur at the head of rules (or in facts which are rules with no body).\n",
    "\n",
    "`ip_var(Name)`\n",
    " - Define a variable that will be used in the constraints.\n",
    "\n",
    "`ip_domain(Lower,Upper)`\n",
    " - Define shared lower and upper bounds for all integer variables in the program.\n",
    "\n",
    "`ip_weight(Constraint,VarName,Weight)`\n",
    " - Define one weighted term in a constraint.\n",
    " \n",
    "`ip_rel(Constraint,Rel)`\n",
    " - Define the type of constraint: `Rel` is `eq`/`ne`/`lt`/`gt`/`le`/`ge`\n",
    "\n",
    "`ip_const(Constraint,Constant)`\n",
    " - Define the final constant in the constraint."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An example integer program\n",
    "Let's model this system of (in)equations:\n",
    "\n",
    "$$\\begin{align}\n",
    "\\text{(a)} && x-2y & = 3 \\\\\n",
    "\\text{(b)} && b & \\gt 2 \\\\\n",
    "\\\\\n",
    "\\text{where} && x,y & \\in [-10,+10]\n",
    "\\end{align}$$\n",
    "\n",
    "By inspection, this system has exactly one solution where $x=9$ and $y=3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting ip-example.lp\n"
     ]
    }
   ],
   "source": [
    "%%file ip-example.lp\n",
    "     \n",
    "% x,y in [-10,+10]\n",
    "ip_var(x).\n",
    "ip_var(y).\n",
    "ip_domain(-10,10).\n",
    "\n",
    "% a:[x-2y = 3]\n",
    "ip_weight(a,x,1).\n",
    "ip_weight(a,y,-2).\n",
    "ip_rel(a,eq).\n",
    "ip_const(a,3).\n",
    "\n",
    "% b:[y > 2]\n",
    "ip_weight(b,y,1).\n",
    "ip_rel(b,gt).\n",
    "ip_const(b,2)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Direct manual encoding in ASP\n",
    "Here's how we might model the system above directly using familiar constructs in ASP. It's easy to understand, but it entails some bloat during grounding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting ip-manual-1.lp\n"
     ]
    }
   ],
   "source": [
    "%%file ip-manual-1.lp\n",
    "\n",
    "% guess values for x and y in the domain\n",
    "1 { x_val(-10..10) } 1.\n",
    "1 { y_val(-10..10) } 1.\n",
    "\n",
    "% enforce constraints\n",
    ":- x_val(X), y_val(Y), not X-2*Y = 3.\n",
    ":- y_val(Y), not Y > 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "clingo version 5.1.0\r\n",
      "Reading from ip-manual-1.lp\r\n",
      "Solving...\r\n",
      "Answer: 1\r\n",
      "x_val(9) y_val(3)\r\n",
      "SATISFIABLE\r\n",
      "\r\n",
      "Models       : 1\r\n",
      "Calls        : 1\r\n",
      "Time         : 0.003s (Solving: 0.00s 1st Model: 0.00s Unsat: 0.00s)\r\n",
      "CPU Time     : 0.000s\r\n"
     ]
    }
   ],
   "source": [
    "!clingo ip-manual-1.lp 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#delayed(1).\r\n",
      "#delayed(2).\r\n",
      ":-y_val(-10),x_val(-10).\r\n",
      ":-y_val(-10),x_val(-9).\r\n",
      ":-y_val(-10),x_val(-8).\r\n",
      ":-y_val(-10),x_val(-7).\r\n",
      ":-y_val(-10),x_val(-6).\r\n",
      ":-y_val(-10),x_val(-5).\r\n",
      ":-y_val(-10),x_val(-4).\r\n",
      ":-y_val(-10),x_val(-3).\r\n"
     ]
    }
   ],
   "source": [
    "!clingo ip-manual-1.lp 0 --text | head"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This works, but we end up building a huge collections of X and Y for which the relations do *not* hold. We might try only representing the cases where they *do* hold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting ip-manual-2.lp\n"
     ]
    }
   ],
   "source": [
    "%%file ip-manual-2.lp\n",
    "\n",
    "% guess values for x and y in the domain\n",
    "1 { x_val(-10..10) } 1.\n",
    "1 { y_val(-10..10) } 1.\n",
    "\n",
    "% enforce constraints\n",
    "a_ok :- x_val(X), y_val(Y), X-2*Y = 3.\n",
    ":- not a_ok.\n",
    "    \n",
    "b_ok :- y_val(Y), Y > 2.\n",
    ":- not b_ok."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "clingo version 5.1.0\r\n",
      "Reading from ip-manual-2.lp\r\n",
      "Solving...\r\n",
      "Answer: 1\r\n",
      "a_ok y_val(3) x_val(9) b_ok\r\n",
      "SATISFIABLE\r\n",
      "\r\n",
      "Models       : 1\r\n",
      "Calls        : 1\r\n",
      "Time         : 0.002s (Solving: 0.00s 1st Model: 0.00s Unsat: 0.00s)\r\n",
      "CPU Time     : 0.000s\r\n"
     ]
    }
   ],
   "source": [
    "!clingo ip-manual-2.lp 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#delayed(1).\r\n",
      "#delayed(2).\r\n",
      "a_ok:-y_val(-6),x_val(-9).\r\n",
      "a_ok:-y_val(-5),x_val(-7).\r\n",
      "a_ok:-y_val(-4),x_val(-5).\r\n",
      "a_ok:-y_val(-3),x_val(-3).\r\n",
      "a_ok:-y_val(-2),x_val(-1).\r\n",
      "a_ok:-y_val(-1),x_val(1).\r\n",
      "a_ok:-y_val(0),x_val(3).\r\n",
      "a_ok:-y_val(1),x_val(5).\r\n"
     ]
    }
   ],
   "source": [
    "!clingo ip-manual-2.lp --text | head"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For an linear equation over $k$ variables with domain size $d$ this approach will take $O(d^k)$ time to produce $O(d^{k-1})$ ground instances. Even if $k$ is only ever $2$, we can't afford large $d$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# IP Solvers (in ASP)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## One-hot Encoding\n",
    "The `ip_assign(Var,Val)` predicate will directly map symbolic variables to integer values. For each integer variable, one `ip_assign/2` term will be present in each answer set. This encoding uses weight constraints in a fairly straightforward manner. Grounding will take $O(\\text{vars} \\cdot \\text{vals})$ time, so try to keep the size of the domain down."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting ip-onehot.lp\n"
     ]
    }
   ],
   "source": [
    "%%file ip-onehot.lp\n",
    "\n",
    "1 { ip_assign(Var,Low..High):ip_domain(Low,High) } 1 :- ip_var(Var).\n",
    "\n",
    ":-  ip_rel(Constraint,eq),\n",
    "    ip_const(Constraint,Constant),\n",
    "    not #sum { Weight*Val,Var: ip_assign(Var,Val), ip_weight(Constraint,Var,Weight) } = Constant.\n",
    "\n",
    ":-  ip_rel(Constraint,ne),\n",
    "    ip_const(Constraint,Constant),\n",
    "    not #sum { Weight*Val,Var: ip_assign(Var,Val), ip_weight(Constraint,Var,Weight) } != Constant.\n",
    "\n",
    ":-  ip_rel(Constraint,lt),\n",
    "    ip_const(Constraint,Constant),\n",
    "    not #sum { Weight*Val,Var: ip_assign(Var,Val), ip_weight(Constraint,Var,Weight) } < Constant.\n",
    "    \n",
    ":-  ip_rel(Constraint,gt),\n",
    "    ip_const(Constraint,Constant),\n",
    "    not #sum { Weight*Val,Var: ip_assign(Var,Val), ip_weight(Constraint,Var,Weight) } > Constant.\n",
    "    \n",
    ":-  ip_rel(Constraint,le),\n",
    "    ip_const(Constraint,Constant),\n",
    "    not #sum { Weight*Val,Var: ip_assign(Var,Val), ip_weight(Constraint,Var,Weight) } <= Constant.\n",
    "    \n",
    ":-  ip_rel(Constraint,ge),\n",
    "    ip_const(Constraint,Constant),\n",
    "    not #sum { Weight*Val,Var: ip_assign(Var,Val), ip_weight(Constraint,Var,Weight) } >= Constant.\n",
    "\n",
    "#show ip_assign/2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "clingo version 5.1.0\r\n",
      "Reading from ip-example.lp ...\r\n",
      "Solving...\r\n",
      "Answer: 1\r\n",
      "ip_assign(x,9) ip_assign(y,3)\r\n",
      "SATISFIABLE\r\n",
      "\r\n",
      "Models       : 1\r\n",
      "Calls        : 1\r\n",
      "Time         : 0.005s (Solving: 0.00s 1st Model: 0.00s Unsat: 0.00s)\r\n",
      "CPU Time     : 0.000s\r\n"
     ]
    }
   ],
   "source": [
    "!clingo ip-example.lp ip-onehot.lp 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bitvector Encoding\n",
    "The `ip_hot(Var,Index)` predicate notes which bits of the value associated with `Var` should be set to one. For each integer variable, several `ip_hot/2` terms could be present in each answer set. This encoding computes the weights for use in the weight constraints from the bit indexes. I don't trust Clingo with numbers beyond $2^{30}$, however. Grounding takes $O(\\text{vars} \\cdot \\log_2(\\text{vals}))$, so let those domains grow! Clingo probably doesn't propagate these pathological weight constraints well, so try to avoid using a large number of integer variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting ip-bitvector.lp\n"
     ]
    }
   ],
   "source": [
    "%%file ip-bitvector.lp\n",
    "\n",
    "#script (python)\n",
    "import math\n",
    "def bits_needed_to_represent(x):\n",
    "    return int(math.ceil(math.log(x.number + 1)/math.log(2)))\n",
    "#end.\n",
    "\n",
    "ip_index(0..@bits_needed_to_represent(High-Low)-1) :- ip_domain(Low,High).\n",
    "    \n",
    "{ ip_hot(V,I) } :- ip_var(V), ip_index(I).\n",
    "    \n",
    ":- ip_domain(Low,High), ip_var(V), not Low #sum { 2**I: ip_hot(V,I) } High.\n",
    "    \n",
    ":-  ip_rel(Constraint,eq),\n",
    "    ip_const(Constraint,Constant),\n",
    "    not #sum {W*(2**I),V: ip_hot(V,I),ip_weight(Constraint,V,W)} = Constant.\n",
    "\n",
    ":-  ip_rel(Constraint,ne),\n",
    "    ip_const(Constraint,Constant),\n",
    "    not #sum {W*(2**I),V: ip_hot(V,I),ip_weight(Constraint,V,W)} != Constant.\n",
    "\n",
    ":-  ip_rel(Constraint,lt),\n",
    "    ip_const(Constraint,Constant),\n",
    "    not #sum {W*(2**I),V: ip_hot(V,I),ip_weight(Constraint,V,W)} < Constant.\n",
    "\n",
    ":-  ip_rel(Constraint,gt),\n",
    "    ip_const(Constraint,Constant),\n",
    "    not #sum {W*(2**I),V: ip_hot(V,I),ip_weight(Constraint,V,W)} > Constant.\n",
    "\n",
    ":-  ip_rel(Constraint,le),\n",
    "    ip_const(Constraint,Constant),\n",
    "    not #sum {W*(2**I),V: ip_hot(V,I),ip_weight(Constraint,V,W)} <= Constant.\n",
    "\n",
    ":-  ip_rel(Constraint,ge),\n",
    "    ip_const(Constraint,Constant),\n",
    "    not #sum {W*(2**I),V: ip_hot(V,I),ip_weight(Constraint,V,W)} >= Constant.\n",
    "\n",
    "#show ip_hot/2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "clingo version 5.1.0\r\n",
      "Reading from ip-example.lp ...\r\n",
      "Solving...\r\n",
      "Answer: 1\r\n",
      "ip_hot(x,0) ip_hot(x,3) ip_hot(y,0) ip_hot(y,1)\r\n",
      "SATISFIABLE\r\n",
      "\r\n",
      "Models       : 1\r\n",
      "Calls        : 1\r\n",
      "Time         : 0.031s (Solving: 0.00s 1st Model: 0.00s Unsat: 0.00s)\r\n",
      "CPU Time     : 0.030s\r\n"
     ]
    }
   ],
   "source": [
    "!clingo ip-example.lp ip-bitvector.lp 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clingcon Encoding\n",
    "This encoding hands off the integer variables and consrtaints to the theory solver in `clingcon`. The rules in this encoding simply construct the relevant theory terms. Even though this will parse and run in `clingo`, you'll need `clingcon` to enforce the constraints. I've read that `clingcon` is very smart about both large domains and large collections of variables. Go wild!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting ip-clingcon.lp\n"
     ]
    }
   ],
   "source": [
    "%%file ip-clingcon.lp\n",
    "\n",
    "#include \"csp.lp\".\n",
    "&dom {Low..High} = V :- ip_var(V), ip_domain(Low,High).\n",
    "\n",
    "&sum {W*V:ip_weight(Constraint,V,W)}  = Constant :- ip_rel(Constraint,eq), ip_const(Constraint,Constant).\n",
    "&sum {W*V:ip_weight(Constraint,V,W)} != Constant :- ip_rel(Constraint,ne), ip_const(Constraint,Constant).\n",
    "&sum {W*V:ip_weight(Constraint,V,W)}  < Constant :- ip_rel(Constraint,lt), ip_const(Constraint,Constant).\n",
    "&sum {W*V:ip_weight(Constraint,V,W)}  > Constant :- ip_rel(Constraint,gt), ip_const(Constraint,Constant).\n",
    "&sum {W*V:ip_weight(Constraint,V,W)} <= Constant :- ip_rel(Constraint,le), ip_const(Constraint,Constant).\n",
    "&sum {W*V:ip_weight(Constraint,V,W)} >= Constant :- ip_rel(Constraint,ge), ip_const(Constraint,Constant).\n",
    "    \n",
    "&show {V/0:ip_var(V)}. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "clingcon version 3.2.1\r\n",
      "Reading from stdin\r\n",
      "Solving...\r\n",
      "Answer: 1\r\n",
      "ip_domain(-10,10) ip_var(x) ip_var(y) ip_const(a,3) ip_const(b,2) ip_rel(a,eq) ip_rel(b,gt) ip_weight(a,x,1) ip_weight(b,y,1) ip_weight(a,y,-2) y=3 x=9\r\n",
      "SATISFIABLE\r\n",
      "\r\n",
      "Models       : 1\r\n",
      "Calls        : 1\r\n",
      "Time         : 0.004s (Solving: 0.00s 1st Model: 0.00s Unsat: 0.00s)\r\n",
      "CPU Time     : 0.000s\r\n"
     ]
    }
   ],
   "source": [
    "# running clingcon in Docker\n",
    "!cat ip-example.lp ip-clingcon.lp | docker run --rm -i rndmcnlly/clingcon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "clingcon version 3.2.1\r\n",
      "Reading from stdin\r\n",
      "Solving...\r\n",
      "Answer: 1\r\n",
      "ip_domain(-10,10) ip_var(x) ip_var(y) ip_const(a,3) ip_const(b,2) ip_rel(a,eq) ip_rel(b,gt) ip_weight(a,x,1) ip_weight(b,y,1) ip_weight(a,y,-2) y=3 x=9\r\n",
      "SATISFIABLE\r\n",
      "\r\n",
      "Models       : 1\r\n",
      "Calls        : 1\r\n",
      "Time         : 0.006s (Solving: 0.00s 1st Model: 0.00s Unsat: 0.00s)\r\n",
      "CPU Time     : 0.010s\r\n"
     ]
    }
   ],
   "source": [
    "# running clingcon using cmd.io\n",
    "!cat ip-example.lp ip-clingcon.lp | ssh cmd rndmcnlly/clingcon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "clingcon version 3.2.1\r\n",
      "Reading from stdin\r\n",
      "Solving...\r\n",
      "Answer: 1\r\n",
      "ip_domain(-10,10) ip_var(x) ip_var(y) ip_const(a,3) ip_const(b,2) ip_rel(a,eq) ip_rel(b,gt) ip_weight(a,x,1) ip_weight(b,y,1) ip_weight(a,y,-2) y=3 x=9\r\n",
      "SATISFIABLE\r\n",
      "\r\n",
      "Models       : 1\r\n",
      "Calls        : 1\r\n",
      "Time         : 0.087s (Solving: 0.00s 1st Model: 0.00s Unsat: 0.00s)\r\n",
      "CPU Time     : 0.080s\r\n"
     ]
    }
   ],
   "source": [
    "# running clingcon on Google Cloud Functions\n",
    "!cat ip-example.lp ip-clingcon.lp | python3 clingcon_client_for_python3.py"
   ]
  }
 ],
 "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}