git.sesse.net Git - nms/blob - tsp/tsp.cpp

   1 #include <stdio.h>
   2 #include <limits.h>
   3 #include <vector>
   4 #include <set>
   5 #include <algorithm>
   6
   7 #define MIN_ROW 1
   8 #define MAX_ROW 75
   9 #define MIN_SWITCH 1
  10 #define MAX_SWITCH 6
  11 #define HEAP_MST 0
  12 #define USE_KOPT 1
  13 #define KOPT_SLICES 3
  14 #define KOPT_ITERATIONS 10000
  15
  16 static const unsigned num_cache_elem = (MAX_ROW * MAX_SWITCH * 2) * (MAX_ROW * MAX_SWITCH * 2);
  17 static unsigned short dist_cache[(MAX_ROW * MAX_SWITCH * 2) * (MAX_ROW * MAX_SWITCH * 2)],
  18         opt_dist_cache[MAX_ROW * MAX_SWITCH * MAX_ROW * MAX_SWITCH],
  19         pess_dist_cache[MAX_ROW * MAX_SWITCH * MAX_ROW * MAX_SWITCH];
  20
  21 inline unsigned short &cache(
  22         unsigned row_from, unsigned switch_from, unsigned side_from,
  23         unsigned row_to, unsigned switch_to, unsigned side_to)
  24 {
  25         return dist_cache[(row_from * MAX_SWITCH * 2 + switch_from * 2 + side_from) * (MAX_ROW * MAX_SWITCH * 2) +
  26                 row_to * MAX_SWITCH * 2 + switch_to * 2 + side_to];
  27 }
  28
  29 inline unsigned short &opt_cache(
  30         unsigned row_from, unsigned switch_from,
  31         unsigned row_to, unsigned switch_to)
  32 {
  33         return opt_dist_cache[(row_from * MAX_SWITCH + switch_from) * (MAX_ROW * MAX_SWITCH) +
  34                 row_to * MAX_SWITCH + switch_to];
  35 }
  36
  37 inline unsigned short &pess_cache(
  38         unsigned row_from, unsigned switch_from,
  39         unsigned row_to, unsigned switch_to)
  40 {
  41         return pess_dist_cache[(row_from * MAX_SWITCH + switch_from) * (MAX_ROW * MAX_SWITCH) +
  42                 row_to * MAX_SWITCH + switch_to];
  43 }
  44
  45 struct order {
  46         unsigned row, num;
  47         int side;
  48         int cost;
  49
  50         bool operator< (const order &other) const
  51         {
  52                 return (cost < other.cost);
  53         }
  54 };
  55
  56 static unsigned best_so_far = UINT_MAX;
  57 order *best_tour;
  58
  59 int distance_switch(unsigned from, unsigned to)
  60 {
  61         /* on the same side of the middle? 9.6m per switch. */
  62         if ((from > 3) == (to > 3)) {
  63                 return abs(from - to) * 96;
  64         }
  65
  66         /* have to cross the border? 25.8m from sw3->sw4 => 16.2m extra gap. */
  67         /* that's _got_ to be wrong. say it's 3m. */
  68         return abs(from - to) * 96 + 30;
  69 }
  70
  71 int distance_middle(unsigned sw, unsigned middle)
  72 {
  73         /* symmetry: 4-5-6 are just mirrored 3-2-1. */
  74         if (middle == 2) {
  75                 if (sw > 3)
  76                         sw = 7 - sw;
  77
  78                 /* estimate 25.8m/2 = 12.9m from sw3 to the middle */
  79                 return 129 + (3 - sw) * 96;
  80         }
  81
  82         /* more symmetry -- getting from 1-6 to the top is like getting from 6-1 to the bottom. */
  83         if (middle == 3) {
  84                 middle = 1;
  85                 sw = 7 - sw;
  86         }
  87
  88         /* guesstimate 4.8m extra to get to the bottom */
  89         if (sw > 3)
  90                 return 48 + 162 + (sw - 1) * 96;
  91         else
  92                 return 48 + (sw - 1) * 96;
  93 }
  94
  95 int distance_row(unsigned from, unsigned to)
  96 {
  97         /* 4.1m per double row, plus gaps */
  98         unsigned base_cost = 41 * abs(from - to);
  99
 100         if ((from <= 9) != (to <= 9))
 101                 base_cost += 25;
 102         if ((from <= 17) != (to <= 17))
 103                 base_cost += 25;
 104         if ((from <= 25) != (to <= 25))
 105                 base_cost += 25;
 106         if ((from <= 34) != (to <= 34))
 107                 base_cost += 25;
 108
 109         /* don't calculate gaps here just yet, just estimate 4.1m per double row */
 110         return base_cost;
 111 }
 112
 113 int pessimistic_distance(int row_from, int switch_from, int row_to, int switch_to)
 114 {
 115         /* we'll need to go to one of the three middles */
 116         int best2 = distance_middle(switch_from, 2) + distance_middle(switch_to, 2);
 117         int distrow = distance_row(row_from, row_to);
 118         if ((switch_from > 3) != (switch_to > 3))
 119                 return best2 + distrow;
 120         if (switch_from > 3) {
 121                 int best3 = distance_middle(switch_from, 3) + distance_middle(switch_to, 3);
 122                 return std::min(best2, best3) + distrow;
 123         } else {
 124                 int best1 = distance_middle(switch_from, 1) + distance_middle(switch_to, 1);
 125                 return std::min(best2, best1) + distrow;
 126         }
 127 }
 128
 129 int distance(int row_from, int switch_from, int side_from, int row_to, int switch_to, int side_to)
 130 {
 131         /* can we just walk directly? */
 132         if (row_from == row_to && side_from == side_to) {
 133                 return distance_switch(switch_from, switch_to);
 134         }
 135
 136         /* can we just switch sides? */
 137         if (row_from + 1 == row_to && side_from == 1 && side_to == 0) {
 138                 return distance_switch(switch_from, switch_to) + distance_row(row_from, row_to) - 31;
 139         }
 140         if (row_from == row_to + 1 && side_from == 0 && side_to == 1) {
 141                 return distance_switch(switch_from, switch_to) + distance_row(row_from, row_to) - 31;
 142         }
 143
 144         return pessimistic_distance(row_from, switch_from, row_to, switch_to);
 145 }
 146
 147 int optimistic_distance(int row_from, int switch_from, int row_to, int switch_to)
 148 {
 149         if (row_from == row_to)
 150                 return distance_switch(switch_from, switch_to);
 151         else if (abs(row_from - row_to) == 1)
 152                 return distance_switch(switch_from, switch_to) + distance_row(row_from, row_to) - 31;
 153         else
 154                 return pessimistic_distance(row_from, switch_from, row_to, switch_to);
 155 }
 156
 157 #if HEAP_MST
 158 // this is, surprisingly enough, _slower_ than the naive variant below, so it's not enabled
 159 struct prim_queue_val {
 160         std::pair<unsigned, unsigned> dst;
 161         int cost;
 162
 163         bool operator< (const prim_queue_val &other) const
 164         {
 165                 return (cost > other.cost);
 166         }
 167 };
 168
 169 // standard O(V^2 log v) prim
 170 int prim_mst(std::set<std::pair<unsigned, unsigned> > &in)
 171 {
 172         std::set<std::pair<unsigned, unsigned> > set2;
 173         std::priority_queue<prim_queue_val> queue;
 174
 175         // pick the first one
 176         std::set<std::pair<unsigned, unsigned> >::iterator start = in.begin();
 177
 178         unsigned row = start->first;
 179         unsigned num = start->second;
 180
 181         set2.insert(*start);
 182
 183         // find all the edges out from it
 184         for (std::set<std::pair<unsigned, unsigned> >::iterator j = in.begin(); j != in.end(); ++j) {
 185                 if (set2.count(*j))
 186                         continue;
 187
 188                 unsigned d = opt_cache(row, num, j->first, j->second);
 189                 prim_queue_val val = { *j, d };
 190                 queue.push(val);
 191         }
 192
 193         unsigned total_cost = 0;
 194         while (set2.size() != in.size()) {
 195 invalid:
 196                 prim_queue_val val = queue.top();
 197                 queue.pop();
 198
 199                 // check if dst is already moved
 200                 if (set2.count(val.dst))
 201                         goto invalid;
 202
 203                 unsigned row = val.dst.first;
 204                 unsigned num = val.dst.second;
 205                 set2.insert(val.dst);
 206
 207                 total_cost += val.cost;
 208
 209                 // find all the edges from this new node
 210                 for (std::set<std::pair<unsigned, unsigned> >::iterator j = in.begin(); j != in.end(); ++j) {
 211                         if (set2.count(*j))
 212                                 continue;
 213
 214                         unsigned d = opt_cache(row, num, j->first, j->second);
 215                         prim_queue_val val = { *j, d };
 216                         queue.push(val);
 217                 }
 218         }
 219
 220         return total_cost;
 221 }
 222 #else
 223 // extremely primitive O(V^3) prim
 224 int prim_mst(std::set<std::pair<unsigned, unsigned> > &set1)
 225 {
 226         std::set<std::pair<unsigned, unsigned> > set2;
 227
 228         // pick the first one
 229         std::set<std::pair<unsigned, unsigned> >::iterator start = set1.begin();
 230         set2.insert(*start);
 231         set1.erase(start);
 232
 233         unsigned total_cost = 0;
 234         while (set1.size() > 0) {
 235                 unsigned best_this_cost = UINT_MAX;
 236                 std::set<std::pair<unsigned, unsigned> >::iterator best_set1;
 237
 238                 for (std::set<std::pair<unsigned, unsigned> >::iterator i = set1.begin(); i != set1.end(); ++i) {
 239                         for (std::set<std::pair<unsigned, unsigned> >::iterator j = set2.begin(); j != set2.end(); ++j) {
 240                                 unsigned d = opt_cache(i->first, i->second, j->first, j->second);
 241                                 if (d < best_this_cost) {
 242                                         best_this_cost = d;
 243                                         best_set1 = i;
 244                                 }
 245                         }
 246                 }
 247
 248                 set2.insert(*best_set1);
 249                 set1.erase(best_set1);
 250                 total_cost += best_this_cost;
 251         }
 252
 253         return total_cost;
 254 }
 255 #endif
 256
 257 void print_tour(std::vector<std::pair<unsigned, unsigned> > &points)
 258 {
 259         std::set<std::pair<unsigned, unsigned> > points_left;
 260         unsigned total_cost = 0;
 261         for (unsigned i = 0; i < points.size(); ++i) {
 262                 points_left.insert(points[i]);
 263         }
 264
 265         for (unsigned i = 0; i < points.size(); ++i) {
 266                 printf("%2u: ", i);
 267                 if (best_tour[i].side == 0)
 268                         printf("%2u-%u (left side)  ", best_tour[i].row, best_tour[i].num);
 269                 else
 270                         printf("%2u-%u (right side) ", best_tour[i].row, best_tour[i].num);
 271                 if (i == 0) {
 272                         printf("           ");
 273                 } else {
 274                         printf("cost=%4u  ", best_tour[i].cost);
 275                         total_cost += best_tour[i].cost;
 276                 }
 277
 278                 // let's see how good the MST heuristics are
 279                 if (i != points.size() - 1) {
 280                         std::set<std::pair<unsigned, unsigned> > mst_tree = points_left;
 281                         printf("mst_bound=%5u, ", prim_mst(mst_tree));
 282
 283                         unsigned rest_cost = 0;
 284                         for (unsigned j = i + 1; j < points.size(); ++j) {
 285                                 rest_cost += best_tour[j].cost;
 286                         }
 287
 288                         printf("rest_cost=%5u", rest_cost);
 289                 }
 290
 291                 printf("\n");
 292
 293                 std::set<std::pair<unsigned, unsigned> >::iterator j = points_left.find(std::make_pair(best_tour[i].row, best_tour[i].num));
 294                 points_left.erase(j);
 295         }
 296         printf("Done (total_cost=%u)\n", total_cost);
 297 }
 298
 299 bool any_opt;
 300 std::pair<unsigned, unsigned> best_reorganization[KOPT_SLICES];
 301
 302 void optimize_specific_tour(order *tour, std::set<std::pair<unsigned, unsigned> > &fragments,
 303         std::vector<std::pair<unsigned, unsigned> > &history, unsigned cost_so_far)
 304 {
 305         if (fragments.size() == 1) {
 306                 if (cost_so_far < best_so_far) {
 307                         any_opt = true;
 308                         memcpy(best_reorganization, &history[0], sizeof(std::pair<unsigned, unsigned>)*history.size());
 309                         best_so_far = cost_so_far;
 310                 }
 311                 return;
 312         }
 313
 314         // match the first point with all the others
 315         std::set<std::pair<unsigned, unsigned> >::iterator start = fragments.begin();
 316         std::pair<unsigned, unsigned> start_val = *start;
 317
 318         std::set<std::pair<unsigned, unsigned> >::iterator i = start;
 319         ++i;
 320
 321         fragments.erase(start);
 322
 323         while (i != fragments.end()) {
 324                 std::pair<unsigned, unsigned> i_val = *i;
 325                 std::set<std::pair<unsigned, unsigned> >::iterator i_copy = i;
 326                 ++i;
 327
 328                 unsigned p1 = start_val.second, p2 = i_val.first;
 329                 history[fragments.size() - 1] = std::make_pair(p1, p2);
 330
 331                 unsigned this_cost = cache(
 332                         tour[p1].row, tour[p1].num, tour[p1].side,
 333                         tour[p2].row, tour[p2].num, tour[p2].side
 334                 );
 335
 336                 std::set<std::pair<unsigned, unsigned> > frag_copy = fragments;
 337                 frag_copy.erase(frag_copy.find(i_val));
 338                 frag_copy.insert(std::make_pair(start_val.first, i_val.second));
 339                 optimize_specific_tour(tour, frag_copy, history, cost_so_far + this_cost);
 340         }
 341
 342         fragments.insert(start_val);
 343 }
 344
 345 void piece_fragment(std::vector<std::pair<unsigned, unsigned> > &points, order *ntour, std::set<std::pair<unsigned, unsigned> > &fragments, unsigned s, unsigned j)
 346 {
 347         // find this fragment, and copy
 348         for (std::set<std::pair<unsigned, unsigned> >::iterator i = fragments.begin(); i != fragments.end(); ++i) {
 349                 if (i->first != s)
 350                         continue;
 351                 memcpy(best_tour + j, ntour + s, sizeof(order) * (i->second - i->first + 1));
 352
 353                 j += (i->second - i->first + 1);
 354
 355                 // find the connection point, if any
 356                 for (unsigned k = 0; k < KOPT_SLICES; ++k) {
 357                         if (best_reorganization[k].first == i->second) {
 358                                 piece_fragment(points, ntour, fragments, best_reorganization[k].second, j);
 359                                 best_tour[j].cost =
 360                                         cache(best_tour[j-1].row, best_tour[j-1].num, best_tour[j-1].side,
 361                                               best_tour[j].row, best_tour[j].num, best_tour[j].side);
 362                                 return;
 363                         }
 364                 }
 365         }
 366 }
 367
 368 void reorganize_tour(std::vector<std::pair<unsigned, unsigned> > &points, std::set<std::pair<unsigned, unsigned> > &fragments)
 369 {
 370         order *ntour = new order[points.size()];
 371         memcpy(ntour, best_tour, sizeof(order) * points.size());
 372
 373         piece_fragment(points, ntour, fragments, 0, 0);
 374
 375         delete[] ntour;
 376 }
 377
 378 void optimize_tour(std::vector<std::pair<unsigned, unsigned> > &points)
 379 {
 380         if (points.size() < KOPT_SLICES * 2) {
 381                 fprintf(stderr, "Warning: Too high KOPT_SLICES set (%u), can't optimize\n", KOPT_SLICES);
 382         }
 383
 384         order *tour = new order[points.size()];
 385         any_opt = false;
 386
 387         for (unsigned i = 0; i < KOPT_ITERATIONS; ++i) {
 388                 int cost = best_so_far;
 389                 memcpy(tour, best_tour, sizeof(order) * points.size());
 390
 391                 for (unsigned j = 0; j < KOPT_SLICES; ++j) {
 392                         // we break the link between RR and RR+1.
 393                         unsigned rand_remove;
 394                         do {
 395                                 rand_remove = (unsigned)((points.size() - 2) * (rand() / (RAND_MAX + 1.0)));
 396                         } while (tour[rand_remove + 1].cost == -1 || (rand_remove != 0 && tour[rand_remove].cost == -1));
 397
 398                         cost -= tour[rand_remove + 1].cost;
 399                         tour[rand_remove + 1].cost = -1;
 400                 }
 401
 402                 /* find the starts and the ends of each point */
 403                 int last_p = 0;
 404                 std::set<std::pair<unsigned, unsigned> > fragments;
 405                 for (unsigned j = 1; j < points.size(); ++j) {
 406                         if (tour[j].cost != -1)
 407                                 continue;
 408                         fragments.insert(std::make_pair(last_p, j - 1));
 409                         last_p = j;
 410                 }
 411                 fragments.insert(std::make_pair(last_p, points.size() - 1));
 412
 413                 /* brute force */
 414                 std::vector<std::pair<unsigned, unsigned> > history(KOPT_SLICES);
 415                 optimize_specific_tour(tour, fragments, history, cost);
 416
 417                 if (any_opt) {
 418                         printf("\nk-opt reduced cost to %u:\n", best_so_far);
 419
 420                         reorganize_tour(points, fragments);
 421                         print_tour(points);
 422                         optimize_tour(points);
 423
 424                         return;
 425                 }
 426         }
 427
 428         delete[] tour;
 429
 430         printf("\nk-opt made no further improvements.\n");
 431 }
 432
 433 unsigned do_tsp(std::vector<std::pair<unsigned, unsigned> > &points, std::set<std::pair<unsigned, unsigned> > &points_left, order *ord, order *temp, unsigned ind, unsigned cost_so_far)
 434 {
 435         if (cost_so_far >= best_so_far)
 436                 return UINT_MAX;
 437         if (ind == points.size()) {
 438                 memcpy(best_tour, ord, sizeof(order) * points.size());
 439                 printf("\nNew best tour found! cost=%u\n", cost_so_far);
 440                 print_tour(points);
 441                 best_so_far = cost_so_far;
 442 #if USE_KOPT
 443                 optimize_tour(points);
 444 #endif
 445                 return 0;
 446         }
 447
 448         unsigned last_row = ord[ind-1].row;
 449         unsigned last_switch = ord[ind-1].num;
 450         unsigned last_side = ord[ind-1].side;
 451
 452         /*
 453          * The minimum spanning tree gives us a reliable lower bound for what we can
 454          * achieve for the rest of the graph, so check it before doing anything else.
 455          */
 456         std::set<std::pair<unsigned, unsigned> > mst_set = points_left;
 457         mst_set.insert(std::make_pair(last_row, last_switch));
 458
 459         unsigned min_rest_cost = prim_mst(mst_set);
 460         if (cost_so_far + min_rest_cost >= best_so_far) {
 461                 return UINT_MAX;
 462         }
 463
 464         /*
 465          * Simple heuristic: always search for the closest points from this one first; that
 466          * will give us a sizable cutoff.
 467          */
 468         unsigned toi = 0;
 469
 470         for (std::set<std::pair<unsigned, unsigned> >::iterator i = points_left.begin(); i != points_left.end(); ++i) {
 471                 /* try both sides */
 472                 temp[toi].row = i->first;
 473                 temp[toi].num = i->second;
 474                 temp[toi].side = 0;
 475                 temp[toi].cost = cache(last_row, last_switch, last_side, i->first, i->second, 0);
 476                 ++toi;
 477
 478                 temp[toi].row = i->first;
 479                 temp[toi].num = i->second;
 480                 temp[toi].side = 1;
 481                 temp[toi].cost = cache(last_row, last_switch, last_side, i->first, i->second, 1);
 482                 ++toi;
 483         }
 484         std::sort(temp, temp + toi);
 485
 486         unsigned best_this_cost = UINT_MAX;
 487         for (unsigned i = 0; i < toi; ++i)
 488         {
 489                 ord[ind] = temp[i];
 490
 491                 std::set<std::pair<unsigned, unsigned> >::iterator j = points_left.find(std::make_pair(temp[i].row, temp[i].num));
 492                 points_left.erase(j);
 493                 unsigned cost_rest = do_tsp(points, points_left, ord, temp + points.size() * 2, ind + 1, cost_so_far + temp[i].cost);
 494                 points_left.insert(std::make_pair(temp[i].row, temp[i].num));
 495
 496                 best_this_cost = std::min(best_this_cost, cost_rest);
 497         }
 498
 499         return best_this_cost;
 500 }
 501
 502 int main()
 503 {
 504         std::vector<std::pair<unsigned, unsigned> > points;
 505         std::set<std::pair<unsigned, unsigned> > points_left;
 506
 507         for ( ;; ) {
 508                 unsigned row, sw;
 509                 if (scanf("%u-%u", &row, &sw) != 2)
 510                         break;
 511
 512                 if (row < MIN_ROW || row > MAX_ROW || sw < MIN_SWITCH || sw > MAX_SWITCH) {
 513                         fprintf(stderr, "%u-%u is out of bounds!\n", row, sw);
 514                         exit(1);
 515                 }
 516
 517                 points.push_back(std::make_pair(row, sw));
 518                 if (points.size() != 1)
 519                         points_left.insert(std::make_pair(row, sw));
 520         }
 521
 522         // precalculate all distances
 523         for (unsigned i = 0; i < points.size(); ++i) {
 524                 for (unsigned j = 0; j < points.size(); ++j) {
 525                         cache(points[i].first, points[i].second, 0, points[j].first, points[j].second, 0) =
 526                                 distance(points[i].first, points[i].second, 0, points[j].first, points[j].second, 0);
 527
 528                         cache(points[i].first, points[i].second, 0, points[j].first, points[j].second, 1) =
 529                                 distance(points[i].first, points[i].second, 0, points[j].first, points[j].second, 1);
 530
 531                         cache(points[i].first, points[i].second, 1, points[j].first, points[j].second, 0) =
 532                                 distance(points[i].first, points[i].second, 1, points[j].first, points[j].second, 0);
 533
 534                         cache(points[i].first, points[i].second, 1, points[j].first, points[j].second, 1) =
 535                                 distance(points[i].first, points[i].second, 1, points[j].first, points[j].second, 1);
 536
 537                         opt_cache(points[i].first, points[i].second, points[j].first, points[j].second) =
 538                                 optimistic_distance(points[i].first, points[i].second, points[j].first, points[j].second);
 539
 540                         pess_cache(points[i].first, points[i].second, points[j].first, points[j].second) =
 541                                 pessimistic_distance(points[i].first, points[i].second, points[j].first, points[j].second);
 542                 }
 543         }
 544
 545         order *ord = new order[points.size()];
 546         best_tour = new order[points.size()];
 547         order *temp = new order[points.size() * points.size() * 4];
 548
 549         /* always start at the first one, left side (hack) */
 550         ord[0].row = points[0].first;
 551         ord[0].num = points[0].second;
 552         ord[0].side = 0;
 553
 554         do_tsp(points, points_left, ord, temp, 1, 0);
 555         printf("All done.\n");
 556 }
 557
 558