source/blog/posts/2010/08/14-memory-allocation/tsp.d

   1 /**\r
   2 A D implementation of the "tsp" Olden benchmark, the traveling\r
   3 salesman problem.\r
   4 \r
   5 R. Karp, "Probabilistic analysis of partitioning algorithms for the\r
   6 traveling-salesman problem in the plane."  Mathematics of Operations Research\r
   7 2(3):209-224, August 1977.\r
   8 \r
   9 Converted to D and optimized by leonardo maffi, V.1.0, Oct 29 2009.\r
  10 \r
  11 Removed output unless an option is passed by Leandro Lucarella, 2010-08-04.\r
  12 Downloaded from http://www.fantascienza.net/leonardo/js/\r
  13                 http://www.fantascienza.net/leonardo/js/dolden_tsp.zip\r
  14 */\r
  15 \r
  16 version (Tango) {\r
  17     import tango.stdc.stdio: printf, fprintf, stderr;\r
  18     import tango.stdc.time: CLOCKS_PER_SEC, clock;\r
  19     import tango.math.Math: sqrt, log;\r
  20 \r
  21     import Integer = tango.text.convert.Integer;\r
  22     alias Integer.parse toInt;\r
  23 } else {\r
  24     import std.c.stdio: printf, fprintf, stderr;\r
  25     import std.c.time: CLOCKS_PER_SEC, clock;\r
  26     import std.math: sqrt, log;\r
  27 \r
  28     version (D_Version2) {\r
  29         import std.conv: to;\r
  30         alias to!(int, char[]) toInt;\r
  31     } else {\r
  32         import std.conv: toInt;\r
  33     }\r
  34 }\r
  35 \r
  36 \r
  37 double myclock() {\r
  38     return clock() / cast(double)CLOCKS_PER_SEC;\r
  39 }\r
  40 \r
  41 \r
  42 /**\r
  43 Basic uniform random generator: Minimal Standard in Park and\r
  44 Miller (1988): "Random Number Generators: Good Ones Are Hard to\r
  45 Find", Comm. of the ACM, 31, 1192-1201.\r
  46 Parameters: m = 2^31-1, a=48271.\r
  47 \r
  48 Very simple portable random number generator adapted\r
  49 from Pascal code by Jesper Lund:\r
  50 http://www.gnu-pascal.de/crystal/gpc/en/mail1390.html\r
  51 */\r
  52 final class Random {\r
  53     const int m = int.max;\r
  54     const int a = 48271;\r
  55     const int q = m / a;\r
  56     const int r = m % a;\r
  57 \r
  58     public int seed;\r
  59 \r
  60     this(int the_seed) {\r
  61         this.seed = the_seed;\r
  62     }\r
  63 \r
  64     public double nextDouble() {\r
  65         int k = seed / q;\r
  66         seed = a * (seed - k * q) - r * k;\r
  67         if (seed < 1)\r
  68         seed += m;\r
  69         return cast(double)seed / m;\r
  70     }\r
  71 \r
  72     public int nextInt(int max) {\r
  73         int n = max + 1;\r
  74         int k = cast(int)(n * this.nextDouble());\r
  75         return (k == n) ? k - 1 : k;\r
  76     }\r
  77 }\r
  78 \r
  79 \r
  80 \r
  81 /**\r
  82 * A class that represents a node in a binary tree.  Each node represents\r
  83 * a city in the TSP benchmark.\r
  84 **/\r
  85 final class Tree {\r
  86     /**\r
  87     * The number of nodes (cities) in this subtree\r
  88     **/\r
  89     private int sz;\r
  90 \r
  91     /**\r
  92     * The coordinates that this node represents\r
  93     **/\r
  94     private double x, y;\r
  95 \r
  96     /**\r
  97     * Left and right child of tree\r
  98     **/\r
  99     private Tree left, right;\r
 100 \r
 101     /**\r
 102     * The next pointer in a linked list of nodes in the subtree.  The list\r
 103     * contains the order of the cities to visit.\r
 104     **/\r
 105     private Tree next;\r
 106 \r
 107     /**\r
 108     * The previous pointer in a linked list of nodes in the subtree. The list\r
 109     * contains the order of the cities to visit.\r
 110     **/\r
 111     private Tree prev;\r
 112 \r
 113     static Random rnd;\r
 114 \r
 115     // used by the random number generator\r
 116     private const double M_E2  = 7.3890560989306502274;\r
 117     private const double M_E3  = 20.08553692318766774179;\r
 118     private const double M_E6  = 403.42879349273512264299;\r
 119     private const double M_E12 = 162754.79141900392083592475;\r
 120 \r
 121     public static void initRnd(int seed) {\r
 122         rnd = new Random(seed);\r
 123     }\r
 124 \r
 125     /**\r
 126     * Construct a Tree node (a city) with the specified size\r
 127     * @param size the number of nodes in the (sub)tree\r
 128     * @param x the x coordinate of the city\r
 129     * @param y the y coordinate of the city\r
 130     * @param left the left subtree\r
 131     * @param right the right subtree\r
 132     **/\r
 133     this(int size, double x, double y, Tree l, Tree r) {\r
 134         sz = size;\r
 135         this.x = x;\r
 136         this.y = y;\r
 137         left = l;\r
 138         right = r;\r
 139         next = null;\r
 140         prev = null;\r
 141     }\r
 142 \r
 143     /**\r
 144     * Find Euclidean distance between this node and the specified node.\r
 145     * @param b the specified node\r
 146     * @return the Euclidean distance between two tree nodes.\r
 147     **/\r
 148     double distance(Tree b) {\r
 149         return sqrt((x - b.x) * (x - b.x) + (y - b.y) * (y - b.y));\r
 150     }\r
 151 \r
 152     /**\r
 153     * Create a list of tree nodes.  Requires root to be the tail of the list.\r
 154     * Only fills in next field, not prev.\r
 155     * @return the linked list of nodes\r
 156     **/\r
 157     Tree makeList() {\r
 158         Tree myleft, myright, tleft, tright;\r
 159         Tree retval = this;\r
 160 \r
 161         // head of left list\r
 162         if (left !is null)\r
 163             myleft = left.makeList();\r
 164         else\r
 165             myleft = null;\r
 166 \r
 167         // head of right list\r
 168         if (right !is null)\r
 169             myright = right.makeList();\r
 170         else\r
 171             myright = null;\r
 172 \r
 173         if (myright !is null) {\r
 174             retval = myright;\r
 175             right.next = this;\r
 176         }\r
 177 \r
 178         if (myleft !is null) {\r
 179             retval = myleft;\r
 180             if (myright !is null)\r
 181                 left.next = myright;\r
 182             else\r
 183                 left.next = this;\r
 184         }\r
 185         next = null;\r
 186 \r
 187         return retval;\r
 188     }\r
 189 \r
 190     /**\r
 191     * Reverse the linked list.  Assumes that there is a dummy "prev"\r
 192     * node at the beginning.\r
 193     **/\r
 194     void reverse() {\r
 195         Tree prev = this.prev;\r
 196         prev.next = null;\r
 197         this.prev = null;\r
 198         Tree back = this;\r
 199         Tree tmp = this;\r
 200 \r
 201         // reverse the list for the other nodes\r
 202         Tree next;\r
 203         for (Tree t = this.next; t !is null; back = t, t = next) {\r
 204             next = t.next;\r
 205             t.next = back;\r
 206             back.prev = t;\r
 207         }\r
 208 \r
 209         // reverse the list for this node\r
 210         tmp.next = prev;\r
 211         prev.prev = tmp;\r
 212     }\r
 213 \r
 214     /**\r
 215     * Use closest-point heuristic from Cormen, Leiserson, and Rivest.\r
 216     * @return a\r
 217     **/\r
 218     Tree conquer() {\r
 219         // create the list of nodes\r
 220         Tree t = makeList();\r
 221 \r
 222         // Create initial cycle\r
 223         Tree cycle = t;\r
 224         t = t.next;\r
 225         cycle.next = cycle;\r
 226         cycle.prev = cycle;\r
 227 \r
 228         // Loop over remaining points\r
 229         Tree donext;\r
 230         for ( ; t !is null; t = donext) {\r
 231             donext = t.next; /* value won't be around later */\r
 232             Tree min = cycle;\r
 233             double mindist = t.distance(cycle);\r
 234             for (Tree tmp = cycle.next; tmp != cycle; tmp = tmp.next) {\r
 235                 double test = tmp.distance(t);\r
 236                 if (test < mindist) {\r
 237                     mindist = test;\r
 238                     min = tmp;\r
 239                 }\r
 240             }\r
 241 \r
 242             Tree next = min.next;\r
 243             Tree prev = min.prev;\r
 244 \r
 245             double mintonext = min.distance(next);\r
 246             double mintoprev = min.distance(prev);\r
 247             double ttonext = t.distance(next);\r
 248             double ttoprev = t.distance(prev);\r
 249 \r
 250             if ((ttoprev - mintoprev) < (ttonext - mintonext)) {\r
 251                 // insert between min and prev\r
 252                 prev.next = t;\r
 253                 t.next = min;\r
 254                 t.prev = prev;\r
 255                 min.prev = t;\r
 256             } else {\r
 257                 next.prev = t;\r
 258                 t.next = next;\r
 259                 min.next = t;\r
 260                 t.prev = min;\r
 261             }\r
 262         }\r
 263 \r
 264         return cycle;\r
 265     }\r
 266 \r
 267     /**\r
 268     * Merge two cycles as per Karp.\r
 269     * @param a a subtree with one cycle\r
 270     * @param b a subtree with the other cycle\r
 271     **/\r
 272     Tree merge(Tree a, Tree b) {\r
 273         // Compute location for first cycle\r
 274         Tree min = a;\r
 275         double mindist = distance(a);\r
 276         Tree tmp = a;\r
 277         for (a = a.next; a != tmp; a = a.next) {\r
 278             double test = distance(a);\r
 279             if (test < mindist) {\r
 280                 mindist = test;\r
 281                 min = a;\r
 282             }\r
 283         }\r
 284 \r
 285         Tree next = min.next;\r
 286         Tree prev = min.prev;\r
 287         double mintonext = min.distance(next);\r
 288         double mintoprev = min.distance(prev);\r
 289         double ttonext   = distance(next);\r
 290         double ttoprev   = distance(prev);\r
 291 \r
 292         Tree p1, n1;\r
 293         double tton1, ttop1;\r
 294         if ((ttoprev - mintoprev) < (ttonext - mintonext)) {\r
 295             // would insert between min and prev\r
 296             p1 = prev;\r
 297             n1 = min;\r
 298             tton1 = mindist;\r
 299             ttop1 = ttoprev;\r
 300         } else {\r
 301             // would insert between min and next\r
 302             p1 = min;\r
 303             n1 = next;\r
 304             ttop1 = mindist;\r
 305             tton1 = ttonext;\r
 306         }\r
 307 \r
 308         // Compute location for second cycle\r
 309         min = b;\r
 310         mindist = distance(b);\r
 311         tmp = b;\r
 312         for (b = b.next; b != tmp; b = b.next) {\r
 313             double test = distance(b);\r
 314             if (test < mindist) {\r
 315                 mindist = test;\r
 316                 min = b;\r
 317             }\r
 318         }\r
 319 \r
 320         next = min.next;\r
 321         prev = min.prev;\r
 322         mintonext = min.distance(next);\r
 323         mintoprev = min.distance(prev);\r
 324         ttonext = this.distance(next);\r
 325         ttoprev = this.distance(prev);\r
 326 \r
 327         Tree p2, n2;\r
 328         double tton2, ttop2;\r
 329         if ((ttoprev - mintoprev) < (ttonext - mintonext)) {\r
 330             // would insert between min and prev\r
 331             p2 = prev;\r
 332             n2 = min;\r
 333             tton2 = mindist;\r
 334             ttop2 = ttoprev;\r
 335         } else {\r
 336             // would insert between min andn ext\r
 337             p2 = min;\r
 338             n2 = next;\r
 339             ttop2 = mindist;\r
 340             tton2 = ttonext;\r
 341         }\r
 342 \r
 343         // Now we have 4 choices to complete:\r
 344         // 1:t,p1 t,p2 n1,n2\r
 345         // 2:t,p1 t,n2 n1,p2\r
 346         // 3:t,n1 t,p2 p1,n2\r
 347         // 4:t,n1 t,n2 p1,p2\r
 348         double n1ton2 = n1.distance(n2);\r
 349         double n1top2 = n1.distance(p2);\r
 350         double p1ton2 = p1.distance(n2);\r
 351         double p1top2 = p1.distance(p2);\r
 352 \r
 353         mindist = ttop1 + ttop2 + n1ton2;\r
 354         int choice = 1;\r
 355 \r
 356         double test = ttop1 + tton2 + n1top2;\r
 357         if (test < mindist) {\r
 358             choice = 2;\r
 359             mindist = test;\r
 360         }\r
 361 \r
 362         test = tton1 + ttop2 + p1ton2;\r
 363         if (test < mindist) {\r
 364             choice = 3;\r
 365             mindist = test;\r
 366         }\r
 367 \r
 368         test = tton1 + tton2 + p1top2;\r
 369         if (test < mindist)\r
 370             choice = 4;\r
 371 \r
 372         switch (choice) {\r
 373             case 1:\r
 374                 // 1:p1,this this,p2 n2,n1 -- reverse 2!\r
 375                 n2.reverse();\r
 376                 p1.next = this;\r
 377                 this.prev = p1;\r
 378                 this.next = p2;\r
 379                 p2.prev = this;\r
 380                 n2.next = n1;\r
 381                 n1.prev = n2;\r
 382                 break;\r
 383 \r
 384             case 2:\r
 385                 // 2:p1,this this,n2 p2,n1 -- OK\r
 386                 p1.next = this;\r
 387                 this.prev = p1;\r
 388                 this.next = n2;\r
 389                 n2.prev = this;\r
 390                 p2.next = n1;\r
 391                 n1.prev = p2;\r
 392                 break;\r
 393 \r
 394             case 3:\r
 395                 // 3:p2,this this,n1 p1,n2 -- OK\r
 396                 p2.next = this;\r
 397                 this.prev = p2;\r
 398                 this.next = n1;\r
 399                 n1.prev = this;\r
 400                 p1.next = n2;\r
 401                 n2.prev = p1;\r
 402                 break;\r
 403 \r
 404             default: // case 4:\r
 405                 // 4:n1,this this,n2 p2,p1 -- reverse 1!\r
 406                 n1.reverse();\r
 407                 n1.next = this;\r
 408                 this.prev = n1;\r
 409                 this.next = n2;\r
 410                 n2.prev = this;\r
 411                 p2.next = p1;\r
 412                 p1.prev = p2;\r
 413                 break;\r
 414         }\r
 415 \r
 416         return this;\r
 417     } // end merge\r
 418 \r
 419     /**\r
 420     * Compute TSP for the tree t. Use conquer for problems <= sz\r
 421     * @param sz the cutoff point for using conquer vs. merge\r
 422     **/\r
 423     Tree tsp(int sz) {\r
 424         if (this.sz <= sz)\r
 425             return conquer();\r
 426         Tree leftval  = left.tsp(sz);\r
 427         Tree rightval = right.tsp(sz);\r
 428         return merge(leftval, rightval);\r
 429     }\r
 430 \r
 431     /**\r
 432     * Print the list of cities to visit from the current node.  The\r
 433     * list is the result of computing the TSP problem.\r
 434     * The list for the root node (city) should contain every other node\r
 435     * (city).\r
 436     **/\r
 437     void printVisitOrder() {\r
 438         printf("x = %.15f y = %.15f\n", x, y);\r
 439         for (Tree tmp = next; tmp != this; tmp = tmp.next)\r
 440             printf("x = %.15f y = %.15f\n", tmp.x, tmp.y);\r
 441     }\r
 442 \r
 443     /**\r
 444     * Computes the total length of the current tour.\r
 445     **/\r
 446     double tourLength() {\r
 447         double total = 0.0;\r
 448         Tree precedent = next;\r
 449         Tree current = next.next;\r
 450         if (current == this)\r
 451             return total;\r
 452 \r
 453         do {\r
 454             total += current.distance(precedent);\r
 455             precedent = current;\r
 456             current = current.next;\r
 457         } while (precedent != this);\r
 458 \r
 459         total += current.distance(this);\r
 460         return total;\r
 461     }\r
 462 \r
 463 \r
 464     // static methods ===============================================\r
 465 \r
 466     /**\r
 467     * Return an estimate of median of n values distributed in [min, max)\r
 468     * @param min the minimum value\r
 469     * @param max the maximum value\r
 470     * @param n\r
 471     * @return an estimate of median of n values distributed in [min, max)\r
 472     **/\r
 473     private static double median(double min, double max, int n) {\r
 474         // get random value in [0.0, 1.0)\r
 475         double t = rnd.nextDouble();\r
 476 \r
 477         double retval;\r
 478         if (t > 0.5)\r
 479             retval = log(1.0 - (2.0 * (M_E12 - 1) * (t - 0.5) / M_E12)) / 12.0;\r
 480         else\r
 481             retval = -log(1.0 - (2.0 * (M_E12 - 1) * t / M_E12)) / 12.0;\r
 482 \r
 483         // We now have something distributed on (-1.0, 1.0)\r
 484         retval = (retval + 1.0) * (max - min) / 2.0;\r
 485         return retval + min;\r
 486     }\r
 487 \r
 488     /**\r
 489     * Get double uniformly distributed over [min,max)\r
 490     * @return double uniformily distributed over [min,max)\r
 491     **/\r
 492     private static double uniform(double min, double max) {\r
 493         // get random value between [0.0, 1.0)\r
 494         double retval = rnd.nextDouble();\r
 495         retval = retval * (max - min);\r
 496         return retval + min;\r
 497     }\r
 498 \r
 499     /**\r
 500     * Builds a 2D tree of n nodes in specified range with dir as primary\r
 501     * axis (false for x, true for y)\r
 502     *\r
 503     * @param n the size of the subtree to create\r
 504     * @param dir the primary axis\r
 505     * @param min_x the minimum x coordinate\r
 506     * @param max_x the maximum x coordinate\r
 507     * @param min_y the minimum y coordinate\r
 508     * @param max_y the maximum y coordinate\r
 509     * @return a reference to the root of the subtree\r
 510     **/\r
 511     public static Tree buildTree(int n, bool dir, double min_x,\r
 512                                  double max_x, double min_y, double max_y) {\r
 513         if (n == 0)\r
 514             return null;\r
 515 \r
 516         Tree left, right;\r
 517         double x, y;\r
 518         if (dir) {\r
 519             dir = !dir;\r
 520             double med = median(min_x, max_x, n);\r
 521             left = buildTree(n/2, dir, min_x, med, min_y, max_y);\r
 522             right = buildTree(n/2, dir, med, max_x, min_y, max_y);\r
 523             x = med;\r
 524             y = uniform(min_y, max_y);\r
 525         } else {\r
 526             dir = !dir;\r
 527             double med = median(min_y,max_y,n);\r
 528             left = buildTree(n/2, dir, min_x, max_x, min_y, med);\r
 529             right = buildTree(n/2, dir, min_x, max_x, med, max_y);\r
 530             y = med;\r
 531             x = uniform(min_x, max_x);\r
 532         }\r
 533 \r
 534         return new Tree(n, x, y, left, right);\r
 535     }\r
 536 }\r
 537 \r
 538 \r
 539 public class TSP {\r
 540     /**\r
 541     * Number of cities in the problem.\r
 542     **/\r
 543     private static int cities;\r
 544 \r
 545     /**\r
 546     * Set to true if the result should be printed\r
 547     **/\r
 548     private static bool printResult = false;\r
 549 \r
 550     /**\r
 551     * Set to true to print informative messages\r
 552     **/\r
 553     private static bool printMsgs = false;\r
 554 \r
 555     /**\r
 556     * The main routine which creates a tree and traverses it.\r
 557     * @param args the arguments to the program\r
 558     **/\r
 559     public static void main(char[] args[]) {\r
 560         if (!parseCmdLine(args))\r
 561             return;\r
 562 \r
 563         if (printMsgs)\r
 564             printf("Building tree of size: %d\n", nextPow2(cities+1) - 1);\r
 565 \r
 566         Tree.initRnd(1);\r
 567 \r
 568         auto t0 = myclock();\r
 569         Tree t = Tree.buildTree(cities, false, 0.0, 1.0, 0.0, 1.0);\r
 570 \r
 571         auto t1 = myclock();\r
 572         t.tsp(150);\r
 573         auto t2 = myclock();\r
 574 \r
 575         if (printMsgs)\r
 576             printf("Total tour length: %.15f\n", t.tourLength());\r
 577         auto t3 = myclock();\r
 578 \r
 579         if (printResult)\r
 580             // if the user specifies, print the final result\r
 581             t.printVisitOrder();\r
 582 \r
 583         if (printMsgs) {\r
 584             printf("Tsp build time: %.2f\n", t1 - t0);\r
 585             printf("Tsp computing time: %.2f\n", t2 - t1);\r
 586             printf("Tsp total time: %.2f\n", t3 - t0);\r
 587         }\r
 588     }\r
 589 \r
 590     private static /*unsigned*/ int nextPow2(/*unsigned*/ int x) {\r
 591         if (x < 0)\r
 592             throw new Exception("x must be >= 0");\r
 593         x = x - 1;\r
 594         x = x | (x >>  1);\r
 595         x = x | (x >>  2);\r
 596         x = x | (x >>  4);\r
 597         x = x | (x >>  8);\r
 598         x = x | (x >> 16);\r
 599        return x + 1;\r
 600     }\r
 601 \r
 602     /**\r
 603     * Parse the command line options.\r
 604     * @param args the command line options.\r
 605     **/\r
 606     private static final bool parseCmdLine(char[] args[]) {\r
 607         int i = 1;\r
 608 \r
 609         while (i < args.length && args[i][0] == '-') {\r
 610             char[] arg = args[i++];\r
 611 \r
 612             if (arg == "-c") {\r
 613                 if (i < args.length)\r
 614                     cities = toInt(args[i++]);\r
 615                 else\r
 616                     throw new Exception("-c requires the size of tree");\r
 617                 if (cities < 1)\r
 618                     throw new Exception("Number of cities must be > 0");\r
 619             } else if (arg == "-p") {\r
 620                 printResult = true;\r
 621             } else if (arg == "-m") {\r
 622                 printMsgs = true;\r
 623             } else if (arg == "-h") {\r
 624                 return usage();\r
 625             }\r
 626         }\r
 627 \r
 628         if (cities == 0)\r
 629             return usage();\r
 630 \r
 631         return true;\r
 632     }\r
 633 \r
 634     /**\r
 635     * The usage routine which describes the program options.\r
 636     **/\r
 637     private static final bool usage() {\r
 638         fprintf(stderr, "usage: tsp_d -c <num> [-p] [-m] [-h]\n");\r
 639         fprintf(stderr, "  -c number of cities (rounds up to the next power of 2 minus 1)\n");\r
 640         fprintf(stderr, "  -p (print the final result)\n");\r
 641         fprintf(stderr, "  -m (print informative messages)\n");\r
 642         fprintf(stderr, "  -h (print this message)\n");\r
 643         return false;\r
 644     }\r
 645 }\r
 646 \r
 647 \r
 648 void main(char[][] args) {\r
 649     TSP.main(args);\r
 650 }\r