]> git.sesse.net Git - nms/blobdiff - tsp/tsp.cpp
Add logging to mbd.
[nms] / tsp / tsp.cpp
index 587573ba6fa93966ab341e80ec17ec733a38752b..a5f36205624fb75533b5ab78a6f94938c779c01e 100644 (file)
@@ -4,6 +4,44 @@
 #include <set>
 #include <algorithm>
 
+#define MIN_ROW 1
+#define MAX_ROW 75
+#define MIN_SWITCH 1
+#define MAX_SWITCH 6
+#define HEAP_MST 0
+#define USE_KOPT 1
+#define KOPT_SLICES 3
+#define KOPT_ITERATIONS 10000
+
+static const unsigned num_cache_elem = (MAX_ROW * MAX_SWITCH * 2) * (MAX_ROW * MAX_SWITCH * 2);
+static unsigned short dist_cache[(MAX_ROW * MAX_SWITCH * 2) * (MAX_ROW * MAX_SWITCH * 2)],
+       opt_dist_cache[MAX_ROW * MAX_SWITCH * MAX_ROW * MAX_SWITCH],
+       pess_dist_cache[MAX_ROW * MAX_SWITCH * MAX_ROW * MAX_SWITCH];
+
+inline unsigned short &cache(
+       unsigned row_from, unsigned switch_from, unsigned side_from,
+       unsigned row_to, unsigned switch_to, unsigned side_to)
+{
+       return dist_cache[(row_from * MAX_SWITCH * 2 + switch_from * 2 + side_from) * (MAX_ROW * MAX_SWITCH * 2) +
+               row_to * MAX_SWITCH * 2 + switch_to * 2 + side_to];
+}
+
+inline unsigned short &opt_cache(
+       unsigned row_from, unsigned switch_from,
+       unsigned row_to, unsigned switch_to)
+{
+       return opt_dist_cache[(row_from * MAX_SWITCH + switch_from) * (MAX_ROW * MAX_SWITCH) +
+               row_to * MAX_SWITCH + switch_to];
+}
+
+inline unsigned short &pess_cache(
+       unsigned row_from, unsigned switch_from,
+       unsigned row_to, unsigned switch_to)
+{
+       return pess_dist_cache[(row_from * MAX_SWITCH + switch_from) * (MAX_ROW * MAX_SWITCH) +
+               row_to * MAX_SWITCH + switch_to];
+}
+
 struct order {
        unsigned row, num;
        int side;
@@ -56,25 +94,24 @@ int distance_middle(unsigned sw, unsigned middle)
 
 int distance_row(unsigned from, unsigned to)
 {
+       /* 4.1m per double row, plus gaps */
+       unsigned base_cost = 41 * abs(from - to);
+       
+       if ((from <= 9) != (to <= 9))
+               base_cost += 25;
+       if ((from <= 17) != (to <= 17))
+               base_cost += 25;
+       if ((from <= 25) != (to <= 25))
+               base_cost += 25;
+       if ((from <= 34) != (to <= 34))
+               base_cost += 25;
+
        /* don't calculate gaps here just yet, just estimate 4.1m per double row */
-       return 41 * abs(from - to);
+       return base_cost;
 }
 
-int distance(int row_from, int switch_from, int side_from, int row_to, int switch_to, int side_to)
+int pessimistic_distance(int row_from, int switch_from, int row_to, int switch_to)
 {
-       /* can we just walk directly? */
-       if (row_from == row_to && side_from == side_to) {
-               return distance_switch(switch_from, switch_to);
-       }
-
-       /* can we just switch sides? */
-       if (row_from + 1 == row_to && side_from == 1 && side_to == 0) {
-               return distance_switch(switch_from, switch_to);
-       }
-       if (row_from == row_to + 1 && side_from == 0 && side_to == 1) {
-               return distance_switch(switch_from, switch_to);
-       }
-       
        /* we'll need to go to one of the three middles */
        int best2 = distance_middle(switch_from, 2) + distance_middle(switch_to, 2);
        int distrow = distance_row(row_from, row_to);
@@ -89,15 +126,101 @@ int distance(int row_from, int switch_from, int side_from, int row_to, int switc
        }
 }
 
+int distance(int row_from, int switch_from, int side_from, int row_to, int switch_to, int side_to)
+{
+       /* can we just walk directly? */
+       if (row_from == row_to && side_from == side_to) {
+               return distance_switch(switch_from, switch_to);
+       }
+       
+       /* can we just switch sides? */
+       if (row_from + 1 == row_to && side_from == 1 && side_to == 0) {
+               return distance_switch(switch_from, switch_to) + distance_row(row_from, row_to) - 31;
+       }
+       if (row_from == row_to + 1 && side_from == 0 && side_to == 1) {
+               return distance_switch(switch_from, switch_to) + distance_row(row_from, row_to) - 31;
+       }
+
+       return pessimistic_distance(row_from, switch_from, row_to, switch_to);
+}      
+
 int optimistic_distance(int row_from, int switch_from, int row_to, int switch_to)
 {
-       if (abs(row_from - row_to) == 1)
+       if (row_from == row_to)
                return distance_switch(switch_from, switch_to);
+       else if (abs(row_from - row_to) == 1)
+               return distance_switch(switch_from, switch_to) + distance_row(row_from, row_to) - 31;
        else
-               return distance(row_from, switch_from, 0, row_to, switch_to, 0);
+               return pessimistic_distance(row_from, switch_from, row_to, switch_to);
 }
 
-// extremely primitive O(V^2) prim
+#if HEAP_MST
+// this is, surprisingly enough, _slower_ than the naive variant below, so it's not enabled
+struct prim_queue_val {
+       std::pair<unsigned, unsigned> dst;
+       int cost;
+
+       bool operator< (const prim_queue_val &other) const
+       {
+               return (cost > other.cost);
+       }
+};
+
+// standard O(V^2 log v) prim
+int prim_mst(std::set<std::pair<unsigned, unsigned> > &in)
+{
+       std::set<std::pair<unsigned, unsigned> > set2;
+       std::priority_queue<prim_queue_val> queue;
+
+       // pick the first one
+       std::set<std::pair<unsigned, unsigned> >::iterator start = in.begin();
+       
+       unsigned row = start->first;
+       unsigned num = start->second;
+
+       set2.insert(*start);
+       
+       // find all the edges out from it
+       for (std::set<std::pair<unsigned, unsigned> >::iterator j = in.begin(); j != in.end(); ++j) {
+               if (set2.count(*j))
+                       continue;
+               
+               unsigned d = opt_cache(row, num, j->first, j->second);
+               prim_queue_val val = { *j, d };
+               queue.push(val);
+       }
+
+       unsigned total_cost = 0;
+       while (set2.size() != in.size()) {
+invalid:
+               prim_queue_val val = queue.top();
+               queue.pop();
+               
+               // check if dst is already moved
+               if (set2.count(val.dst))
+                       goto invalid;
+       
+               unsigned row = val.dst.first;
+               unsigned num = val.dst.second;
+               set2.insert(val.dst);
+
+               total_cost += val.cost;
+
+               // find all the edges from this new node
+               for (std::set<std::pair<unsigned, unsigned> >::iterator j = in.begin(); j != in.end(); ++j) {
+                       if (set2.count(*j))
+                               continue;
+                       
+                       unsigned d = opt_cache(row, num, j->first, j->second);
+                       prim_queue_val val = { *j, d };
+                       queue.push(val);
+               }
+       }
+
+       return total_cost;
+}
+#else
+// extremely primitive O(V^3) prim
 int prim_mst(std::set<std::pair<unsigned, unsigned> > &set1)
 {
        std::set<std::pair<unsigned, unsigned> > set2;
@@ -114,7 +237,7 @@ int prim_mst(std::set<std::pair<unsigned, unsigned> > &set1)
                
                for (std::set<std::pair<unsigned, unsigned> >::iterator i = set1.begin(); i != set1.end(); ++i) {
                        for (std::set<std::pair<unsigned, unsigned> >::iterator j = set2.begin(); j != set2.end(); ++j) {
-                               unsigned d = optimistic_distance(i->first, i->second, j->first, j->second);
+                               unsigned d = opt_cache(i->first, i->second, j->first, j->second);
                                if (d < best_this_cost) {
                                        best_this_cost = d;
                                        best_set1 = i;
@@ -129,16 +252,18 @@ int prim_mst(std::set<std::pair<unsigned, unsigned> > &set1)
 
        return total_cost;
 }
-
+#endif
 
 void print_tour(std::vector<std::pair<unsigned, unsigned> > &points)
 {
        std::set<std::pair<unsigned, unsigned> > points_left;
+       unsigned total_cost = 0;
        for (unsigned i = 0; i < points.size(); ++i) {
                points_left.insert(points[i]);
        }
        
        for (unsigned i = 0; i < points.size(); ++i) {
+               printf("%2u: ", i);
                if (best_tour[i].side == 0)
                        printf("%2u-%u (left side)  ", best_tour[i].row, best_tour[i].num);
                else
@@ -147,6 +272,7 @@ void print_tour(std::vector<std::pair<unsigned, unsigned> > &points)
                        printf("           ");
                } else {
                        printf("cost=%4u  ", best_tour[i].cost);
+                       total_cost += best_tour[i].cost;
                }
 
                // let's see how good the MST heuristics are
@@ -167,6 +293,141 @@ void print_tour(std::vector<std::pair<unsigned, unsigned> > &points)
                std::set<std::pair<unsigned, unsigned> >::iterator j = points_left.find(std::make_pair(best_tour[i].row, best_tour[i].num));
                points_left.erase(j);
        }
+       printf("Done (total_cost=%u)\n", total_cost);
+}
+
+bool any_opt;
+std::pair<unsigned, unsigned> best_reorganization[KOPT_SLICES];
+
+void optimize_specific_tour(order *tour, std::set<std::pair<unsigned, unsigned> > &fragments,
+       std::vector<std::pair<unsigned, unsigned> > &history, unsigned cost_so_far)
+{
+       if (fragments.size() == 1) {
+               if (cost_so_far < best_so_far) {
+                       any_opt = true;
+                       memcpy(best_reorganization, &history[0], sizeof(std::pair<unsigned, unsigned>)*history.size());
+                       best_so_far = cost_so_far;
+               }
+               return;
+       }
+
+       // match the first point with all the others
+       std::set<std::pair<unsigned, unsigned> >::iterator start = fragments.begin();
+       std::pair<unsigned, unsigned> start_val = *start;
+
+       std::set<std::pair<unsigned, unsigned> >::iterator i = start;
+       ++i;
+
+       fragments.erase(start);
+
+       while (i != fragments.end()) {
+               std::pair<unsigned, unsigned> i_val = *i;
+               std::set<std::pair<unsigned, unsigned> >::iterator i_copy = i;
+               ++i;
+               
+               unsigned p1 = start_val.second, p2 = i_val.first;
+               history[fragments.size() - 1] = std::make_pair(p1, p2);
+               
+               unsigned this_cost = cache(
+                       tour[p1].row, tour[p1].num, tour[p1].side,
+                       tour[p2].row, tour[p2].num, tour[p2].side
+               );
+
+               std::set<std::pair<unsigned, unsigned> > frag_copy = fragments;
+               frag_copy.erase(frag_copy.find(i_val));
+               frag_copy.insert(std::make_pair(start_val.first, i_val.second));
+               optimize_specific_tour(tour, frag_copy, history, cost_so_far + this_cost);
+       }
+       
+       fragments.insert(start_val);
+}
+
+void piece_fragment(std::vector<std::pair<unsigned, unsigned> > &points, order *ntour, std::set<std::pair<unsigned, unsigned> > &fragments, unsigned s, unsigned j)
+{
+       // find this fragment, and copy
+       for (std::set<std::pair<unsigned, unsigned> >::iterator i = fragments.begin(); i != fragments.end(); ++i) {
+               if (i->first != s)
+                       continue;
+               memcpy(best_tour + j, ntour + s, sizeof(order) * (i->second - i->first + 1));
+               
+               j += (i->second - i->first + 1);
+
+               // find the connection point, if any
+               for (unsigned k = 0; k < KOPT_SLICES; ++k) {
+                       if (best_reorganization[k].first == i->second) {
+                               piece_fragment(points, ntour, fragments, best_reorganization[k].second, j);
+                               best_tour[j].cost = 
+                                       cache(best_tour[j-1].row, best_tour[j-1].num, best_tour[j-1].side,
+                                             best_tour[j].row, best_tour[j].num, best_tour[j].side);
+                               return;
+                       }
+               }
+       }
+}
+
+void reorganize_tour(std::vector<std::pair<unsigned, unsigned> > &points, std::set<std::pair<unsigned, unsigned> > &fragments)
+{
+       order *ntour = new order[points.size()];
+       memcpy(ntour, best_tour, sizeof(order) * points.size());
+
+       piece_fragment(points, ntour, fragments, 0, 0);
+
+       delete[] ntour;
+}
+
+void optimize_tour(std::vector<std::pair<unsigned, unsigned> > &points)
+{
+       if (points.size() < KOPT_SLICES * 2) {
+               fprintf(stderr, "Warning: Too high KOPT_SLICES set (%u), can't optimize\n", KOPT_SLICES);
+       }
+       
+       order *tour = new order[points.size()];
+       any_opt = false;
+       
+       for (unsigned i = 0; i < KOPT_ITERATIONS; ++i) {
+               int cost = best_so_far;
+               memcpy(tour, best_tour, sizeof(order) * points.size());
+               
+               for (unsigned j = 0; j < KOPT_SLICES; ++j) {
+                       // we break the link between RR and RR+1.
+                       unsigned rand_remove;
+                       do {
+                               rand_remove = (unsigned)((points.size() - 2) * (rand() / (RAND_MAX + 1.0)));
+                       } while (tour[rand_remove + 1].cost == -1 || (rand_remove != 0 && tour[rand_remove].cost == -1));
+                       
+                       cost -= tour[rand_remove + 1].cost;
+                       tour[rand_remove + 1].cost = -1;
+               }
+
+               /* find the starts and the ends of each point */
+               int last_p = 0;
+               std::set<std::pair<unsigned, unsigned> > fragments;
+               for (unsigned j = 1; j < points.size(); ++j) {
+                       if (tour[j].cost != -1)
+                               continue;
+                       fragments.insert(std::make_pair(last_p, j - 1));
+                       last_p = j;
+               }
+               fragments.insert(std::make_pair(last_p, points.size() - 1));
+               
+               /* brute force */
+               std::vector<std::pair<unsigned, unsigned> > history(KOPT_SLICES);
+               optimize_specific_tour(tour, fragments, history, cost);
+
+               if (any_opt) {
+                       printf("\nk-opt reduced cost to %u:\n", best_so_far);
+
+                       reorganize_tour(points, fragments);
+                       print_tour(points);
+                       optimize_tour(points);
+
+                       return;
+               }
+       }
+
+       delete[] tour;
+
+       printf("\nk-opt made no further improvements.\n");
 }
 
 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)
@@ -178,41 +439,48 @@ unsigned do_tsp(std::vector<std::pair<unsigned, unsigned> > &points, std::set<st
                printf("\nNew best tour found! cost=%u\n", cost_so_far);
                print_tour(points);
                best_so_far = cost_so_far;
+#if USE_KOPT           
+               optimize_tour(points);
+#endif
                return 0;
        }
+       
+       unsigned last_row = ord[ind-1].row;
+       unsigned last_switch = ord[ind-1].num;
+       unsigned last_side = ord[ind-1].side;
+
+       /* 
+        * The minimum spanning tree gives us a reliable lower bound for what we can
+        * achieve for the rest of the graph, so check it before doing anything else.
+        */
+       std::set<std::pair<unsigned, unsigned> > mst_set = points_left;
+       mst_set.insert(std::make_pair(last_row, last_switch));
+       
+       unsigned min_rest_cost = prim_mst(mst_set);
+       if (cost_so_far + min_rest_cost >= best_so_far) {
+               return UINT_MAX;
+       }
 
        /* 
         * Simple heuristic: always search for the closest points from this one first; that
         * will give us a sizable cutoff.
         */
        unsigned toi = 0;
-       unsigned last_row = ord[ind-1].row;
-       unsigned last_switch = ord[ind-1].num;
-       unsigned last_side = ord[ind-1].side;
-       
-       std::set<std::pair<unsigned, unsigned> > mst_set = points_left;
-       mst_set.insert(std::make_pair(last_row, last_switch));
        
        for (std::set<std::pair<unsigned, unsigned> >::iterator i = points_left.begin(); i != points_left.end(); ++i) {
                /* try both sides */
                temp[toi].row = i->first;
                temp[toi].num = i->second;
                temp[toi].side = 0;
-               temp[toi].cost = distance(last_row, last_switch, last_side, i->first, i->second, 0);
+               temp[toi].cost = cache(last_row, last_switch, last_side, i->first, i->second, 0);
                ++toi;
 
                temp[toi].row = i->first;
                temp[toi].num = i->second;
                temp[toi].side = 1;
-               temp[toi].cost = distance(last_row, last_switch, last_side, i->first, i->second, 1);
+               temp[toi].cost = cache(last_row, last_switch, last_side, i->first, i->second, 1);
                ++toi;
        }
-
-       unsigned min_rest_cost = prim_mst(mst_set);
-       if (cost_so_far + min_rest_cost >= best_so_far) {
-               return UINT_MAX;
-       }
-       
        std::sort(temp, temp + toi);
 
        unsigned best_this_cost = UINT_MAX;
@@ -241,11 +509,39 @@ int main()
                if (scanf("%u-%u", &row, &sw) != 2)
                        break;
 
+               if (row < MIN_ROW || row > MAX_ROW || sw < MIN_SWITCH || sw > MAX_SWITCH) {
+                       fprintf(stderr, "%u-%u is out of bounds!\n", row, sw);
+                       exit(1);
+               }
+
                points.push_back(std::make_pair(row, sw));
                if (points.size() != 1)
                        points_left.insert(std::make_pair(row, sw));
        }
 
+       // precalculate all distances
+       for (unsigned i = 0; i < points.size(); ++i) {
+               for (unsigned j = 0; j < points.size(); ++j) {
+                       cache(points[i].first, points[i].second, 0, points[j].first, points[j].second, 0) =
+                               distance(points[i].first, points[i].second, 0, points[j].first, points[j].second, 0);
+                       
+                       cache(points[i].first, points[i].second, 0, points[j].first, points[j].second, 1) =
+                               distance(points[i].first, points[i].second, 0, points[j].first, points[j].second, 1);
+                       
+                       cache(points[i].first, points[i].second, 1, points[j].first, points[j].second, 0) =
+                               distance(points[i].first, points[i].second, 1, points[j].first, points[j].second, 0);
+                       
+                       cache(points[i].first, points[i].second, 1, points[j].first, points[j].second, 1) =
+                               distance(points[i].first, points[i].second, 1, points[j].first, points[j].second, 1);
+                       
+                       opt_cache(points[i].first, points[i].second, points[j].first, points[j].second) =
+                               optimistic_distance(points[i].first, points[i].second, points[j].first, points[j].second);
+                       
+                       pess_cache(points[i].first, points[i].second, points[j].first, points[j].second) =
+                               pessimistic_distance(points[i].first, points[i].second, points[j].first, points[j].second);
+               }
+       }
+       
        order *ord = new order[points.size()];
        best_tour = new order[points.size()];
        order *temp = new order[points.size() * points.size() * 4];