Sorting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Jo's $O(n^3)$ sorting formulation\n",
    "(I've renamed the `elements` predicate to `element`.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting jo-n3.lp\n"
     ]
    }
   ],
   "source": [
    "%%file jo-n3.lp\n",
    "\n",
    "#const elements=15.\n",
    "pos(1..elements).\n",
    "element(X) :- pos(Y), X = Y*2. %Can be anything, really, but went for a quick and dirty way of getting the list\n",
    "{sort(X,Y) : pos(Y) } = 1 :- element(X).\n",
    ":- sort(X1,Y), sort(X2,Y), X1 != X2.\n",
    ":- sort(X1,Y), sort(X2,Y+1), X2 < X1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When `elements` is 10, we get a quick solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time         : 0.091s (Solving: 0.09s 1st Model: 0.09s Unsat: 0.00s)\r\n",
      "CPU Time     : 0.090s\r\n"
     ]
    }
   ],
   "source": [
    "!clingo jo-n3.lp -c elements=10 -q | grep Time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But things bog down when we consider just 15 elements again. Interestingly the bottleneck is the solver rather than the grounder. Wow!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "clingo version 5.1.0\n",
      "Reading from jo-n3.lp\n",
      "Solving...\n",
      "*** Info : (clingo): INTERRUPTED by signal!\n",
      "UNKNOWN\n",
      "\n",
      "TIME LIMIT   : 1\n",
      "Models       : 0+\n",
      "Calls        : 1\n",
      "Time         : 10.001s (Solving: 9.99s 1st Model: 0.00s Unsat: 0.00s)\n",
      "CPU Time     : 9.990s\n"
     ]
    }
   ],
   "source": [
    "!clingo jo-n3.lp -c elements=15 -q --time-limit=10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adam's $O(n^4)$ fix\n",
    "I wasn't expecting that the solving time would blow up in the encoding above. I think it has to do with how the last constraint as formulated: adjacent positions must be in order. Suppose a position is left blank because we are part way through the search and we haven't assigned every element an position yet. We might have perfectly sorted items before and after the gap, but when we finally plug the gap, we'll find out that some of the items before the gap were supposed to be after the gap all along. When sorting k elements with this formulation, half-way through, there are $k/2$ blanks. Between the blanks, we can form up to about $k/2$ sorted sublists of the k variables and put them between the blanks without triggering the integrity constraint. There are $k^{k/2}$ ways of doing this, and precisely one of those is the correct way. If the solver were making random choices, these aren't good chances of stumbling across the right order.\n",
    "\n",
    "This suggests one quick fix. Let's make the last constraint even more inefficient for grounding with the hope that it improves solving time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting adam-n4.lp\n"
     ]
    }
   ],
   "source": [
    "%%file adam-n4.lp\n",
    "\n",
    "#const elements=15.\n",
    "pos(1..elements).\n",
    "element(X) :- pos(Y), X = Y*2. %Can be anything, really, but went for a quick and dirty way of getting the list\n",
    "{sort(X,Y) : pos(Y) } = 1 :- element(X).\n",
    ":- sort(X1,Y), sort(X2,Y), X1 != X2.\n",
    ":- sort(X1,Y1), sort(X2,Y2), X2 < X1, Y1 >= Y2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time         : 0.031s (Solving: 0.00s 1st Model: 0.00s Unsat: 0.00s)\r\n",
      "CPU Time     : 0.030s\r\n"
     ]
    }
   ],
   "source": [
    "!clingo adam-n4.lp -q | grep Time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can set `elements` to 50 and still get solutions with almost no search time. I won't set it higher because of that quartic scaling in grounding time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time         : 4.679s (Solving: 0.23s 1st Model: 0.23s Unsat: 0.00s)\r\n",
      "CPU Time     : 4.670s\r\n"
     ]
    }
   ],
   "source": [
    "!clingo adam-n4.lp -c elements=50 -q | grep Time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adam's $O(n^2)$ fix\n",
    "The first integrity constraint in the original formulation that checks that two different elements don't get assigned the same position generates $O(n^3)$ ground rules. If we want to be competitive with traditional sorting techniques, we should try avoid any rules that generate more than $O(n \\lg(n)))$ ground instances (or at worst $O(n^2)$).\n",
    "\n",
    "Because we're producing a predicate like `sort(X,Y)` as output, we're signed up for $O(n^2)$ scaling from the start. To do better than this, we'd have to work with a binary encoding.\n",
    "\n",
    "In the same way Jo says every X has exactly one Y, we can say every Y has exactly one X. This yields $O(n)$ instances (with a $O(n)$-sized aggregate in the head of each rule. Now `sort(X,Y)` witnesses some permutation of the input values.\n",
    "\n",
    "Only one of those permutations is a valid sorting. We want the one that puts the lowest-value inputs the lowest-possible indexes. We'll formulate this as a prioritized weighted minimization. Effectively, this will implement [selection sort](https://en.wikipedia.org/wiki/Selection_sort)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting adam-n2-min.lp\n"
     ]
    }
   ],
   "source": [
    "%%file adam-n2-min.lp\n",
    "\n",
    "#const elements=15.\n",
    "pos(1..elements).\n",
    "element(X) :- pos(Y), X = Y*2. %Can be anything, really, but went for a quick and dirty way of getting the list\n",
    "{sort(X,Y) : pos(Y) } = 1 :- element(X).\n",
    "{sort(X,Y) : element(X) } = 1 :- pos(Y).\n",
    "#minimize {X@elements-Y: sort(X,Y)}."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this encoding, we can bump up `elements` to 400. As expected for a problem that doesn't inherently require difficult search, all of the time is spent on grounding. I've set the optimization strategy to \"unsatisfiable-core based optimization\". In the default \"branch and bound\" strategy, the solver finds one solution, and then iteratively adds constraints until there are no better solutions. In the unsat-core strategy, the solver proves there is no zero cost solution, and then iteratively relaxes the problem until there is a solution (then known optimal). For this sorting task, it becomes satisfiable after one step at each priority (so only `elements` backtracks are required in total). After $O(n^2)$ grounding work, solving only takes $O(n)$ time here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "clingo version 5.1.0\n",
      "Reading from adam-n2-min.lp\n",
      "Solving...\n",
      "OPTIMUM FOUND\n",
      "\n",
      "Models       : 379\n",
      "  Optimum    : yes\n",
      "Optimization : 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 72 74 76 78 80 82 84 86 88 90 92 94 96 98 100 102 104 106 108 110 112 114 116 118 120 122 124 126 128 130 132 134 136 138 140 142 144 146 148 150 152 154 156 158 160 162 164 166 168 170 172 174 176 178 180 182 184 186 188 190 192 194 196 198 200 202 204 206 208 210 212 214 216 218 220 222 224 226 228 230 232 234 236 238 240 242 244 246 248 250 252 254 256 258 260 262 264 266 268 270 272 274 276 278 280 282 284 286 288 290 292 294 296 298 300 302 304 306 308 310 312 314 316 318 320 322 324 326 328 330 332 334 336 338 340 342 344 346 348 350 352 354 356 358 360 362 364 366 368 370 372 374 376 378 380 382 384 386 388 390 392 394 396 398 400 402 404 406 408 410 412 414 416 418 420 422 424 426 428 430 432 434 436 438 440 442 444 446 448 450 452 454 456 458 460 462 464 466 468 470 472 474 476 478 480 482 484 486 488 490 492 494 496 498 500 502 504 506 508 510 512 514 516 518 520 522 524 526 528 530 532 534 536 538 540 542 544 546 548 550 552 554 556 558 560 562 564 566 568 570 572 574 576 578 580 582 584 586 588 590 592 594 596 598 600 602 604 606 608 610 612 614 616 618 620 622 624 626 628 630 632 634 636 638 640 642 644 646 648 650 652 654 656 658 660 662 664 666 668 670 672 674 676 678 680 682 684 686 688 690 692 694 696 698 700 702 704 706 708 710 712 714 716 718 720 722 724 726 728 730 732 734 736 738 740 742 744 746 748 750 752 754 756 758 760 762 764 766 768 770 772 774 776 778 780 782 784 786 788 790 792 794 796 798 800\n",
      "Calls        : 1\n",
      "Time         : 6.289s (Solving: 3.97s 1st Model: 0.04s Unsat: 0.00s)\n",
      "CPU Time     : 6.290s\n"
     ]
    }
   ],
   "source": [
    "!clingo adam-n2-min.lp -c elements=400 --opt-strategy=usc -q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Does setting it up as optimization feel cheap?** Here's another cheap solution: I could instead set it up as heuristically informed search. We can safely stop after one solution because we know the elements were selected in order. This search encounters 0 conflicts, and is thus very fast. (Really, we're just abusing a sorting implementation inside of the heuristic system to solve the entire problem.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting adam-n2-heu.lp\n"
     ]
    }
   ],
   "source": [
    "%%file adam-n2-heu.lp\n",
    "\n",
    "#const elements=15.\n",
    "pos(1..elements).\n",
    "element(X) :- pos(Y), X = Y*2. %Can be anything, really, but went for a quick and dirty way of getting the list\n",
    "{sort(X,Y) : pos(Y) } = 1 :- element(X).\n",
    "{sort(X,Y) : element(X) } = 1 :- pos(Y). \n",
    "#heuristic sort(X,Y):element(X),pos(Y). [X@Y,level]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time         : 2.122s (Solving: 0.04s 1st Model: 0.04s Unsat: 0.00s)\r\n",
      "CPU Time     : 2.110s\r\n"
     ]
    }
   ],
   "source": [
    "!clingo adam-n2-heu.lp -c elements=400 -q | grep Time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adam's $O(n \\lg(n))$ reformulation\n",
    "I won't actually formulate this one here, but I'll sketch the idea.\n",
    "\n",
    "The big problem above was that having the `sort(X,Y)` predicate already cost us $O(n^2)$ time in grounding. We just want to say that the input elements have some permutation, and that adjacent elements in that permutation are ordered.\n",
    "\n",
    "Let's talk about this as a circuit design problem. We have $n$ input integers represented with $b$ bits each. An [omega network](https://en.wikipedia.org/wiki/Omega_network) with $\\lg(n)$ stages can form any permutation of those $n$ values. Each switch (of which there are $n\\lg(n)/2$ total) in the network nondeterministically decides to keep or swap the order of inputs -- this involves only $O(b)$ gates. To check that the final outputs of the network are ordered, let's put a comparator between adjacent outputs. Each comparator costs $O(b)$ gates, and there are $n-1$ of them. This circuit has $O(b n \\lg(n))$ total gates, so we can model in with an answer set program involving only as many ground instances. If $b = \\lg(n)$, were still competitive with something like mergesort with ran in $n \\lg(n)$ time only when comparing numbers was assumed to cost unit time.\n",
    "\n",
    "Right now, we only know this circuit can be modeled by an answer set program that *grounds* in $O(n \\lg(n))$ space. We the solving time could be exponential in this size. If we are allowed to sprikle in heuristics, we can deterministically tell the individual switch how to act correctly on the first try.\n",
    "\n",
    "Alternatively, we could have just built a [sorting network](https://en.wikipedia.org/wiki/Sorting_network) from the start and had the whole solution derived at ground time (at the cost of an extra factor of $\\lg(n)$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reflection\n",
    "If a problem can be cleanly stated and solved in deterministic [quasilinear time](https://en.wikipedia.org/wiki/Time_complexity#Quasilinear_time) by a well known algorithm, touching ASP is often just a distraction. ASP is interesting when the problem is messy (we need to reformulate it many many times) or the problem can't be scalably solved in deterministic time. ASP lets us try to formulate a nondeterminstic solution (one that is efficient if you can always make the right guesses). From there, we can hope the solver's heuristics (or the ones we can provide) provide enough guidance for efficient search.\n",
    "\n",
    "ASP usually targets problems where the gap between the best known deterministic and nondeterministic problems is exponential. I have a hunch that there are will also be problems where just a polynomial gap is big enough to make ASP attractive too. There might be problems were we know of a determistic algorithm of unattractive degree, say $O(n^3)$, where we can also find a better nondeterministic algorithm, where solutions can be verified in $O(n\\lg(n))$ time perhaps, for which the common heuristics do well.\n",
    "\n",
    "[Primality testing](https://en.wikipedia.org/wiki/AKS_primality_test) is one such area. I have no intuition for solver performance here."
   ]
  }
 ],
 "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
}