X-Git-Url: https://git.sesse.net/?a=blobdiff_plain;f=tsp%2Ftsp.cpp;h=a5f36205624fb75533b5ab78a6f94938c779c01e;hb=40d6c7e02855fa3bb3cf35805efe853051c88ea5;hp=a87585f6f321695ced5b6406bc08802df7daac09;hpb=19bc8938f3298412f9e1bc295cae941f802951a5;p=nms diff --git a/tsp/tsp.cpp b/tsp/tsp.cpp index a87585f..a5f3620 100644 --- a/tsp/tsp.cpp +++ b/tsp/tsp.cpp @@ -4,15 +4,50 @@ #include #include +#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 pt; + unsigned row, num; int side; -}; -struct try_order { - order o; int cost; - bool operator< (const try_order &other) const + bool operator< (const order &other) const { return (cost < other.cost); } @@ -59,8 +94,36 @@ 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 pessimistic_distance(int row_from, int switch_from, int row_to, int 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); + if ((switch_from > 3) != (switch_to > 3)) + return best2 + distrow; + if (switch_from > 3) { + int best3 = distance_middle(switch_from, 3) + distance_middle(switch_to, 3); + return std::min(best2, best3) + distrow; + } else { + int best1 = distance_middle(switch_from, 1) + distance_middle(switch_to, 1); + return std::min(best2, best1) + distrow; + } } int distance(int row_from, int switch_from, int side_from, int row_to, int switch_to, int side_to) @@ -69,37 +132,98 @@ int distance(int row_from, int switch_from, int side_from, int row_to, int switc 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 == -1) { - return distance_switch(switch_from, switch_to); + 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 == -1 && side_to == 1) { - 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) + distance_row(row_from, row_to) - 31; } - - /* we'll need to go to one of the three middles */ - int best1 = distance_middle(switch_from, 1) + distance_middle(switch_to, 1); - int best2 = distance_middle(switch_from, 2) + distance_middle(switch_to, 2); - int best3 = distance_middle(switch_from, 3) + distance_middle(switch_to, 3); - return std::min(std::min(best1, best2), best3) + distance_row(row_from, row_to); -} + + 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, -1, row_to, switch_to, -1); + return pessimistic_distance(row_from, switch_from, row_to, switch_to); } -// extremely primitive O(V^2) prim -int prim_mst(std::vector > &points, try_order *temp, unsigned toi) +#if HEAP_MST +// this is, surprisingly enough, _slower_ than the naive variant below, so it's not enabled +struct prim_queue_val { + std::pair 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 > &in) { - std::set > set1, set2; + std::set > set2; + std::priority_queue queue; - for (unsigned i = 0; i < toi; ++i) - set1.insert(points[temp[i].o.pt]); + // pick the first one + std::set >::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 >::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 >::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 > &set1) +{ + std::set > set2; // pick the first one std::set >::iterator start = set1.begin(); @@ -113,7 +237,7 @@ int prim_mst(std::vector > &points, try_order *tem for (std::set >::iterator i = set1.begin(); i != set1.end(); ++i) { for (std::set >::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; @@ -128,33 +252,185 @@ int prim_mst(std::vector > &points, try_order *tem return total_cost; } - +#endif void print_tour(std::vector > &points) { + std::set > 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) { - if (best_tour[i].side == -1) - printf("%2u-%u (left side) ", points[best_tour[i].pt].first, - points[best_tour[i].pt].second); + printf("%2u: ", i); + if (best_tour[i].side == 0) + printf("%2u-%u (left side) ", best_tour[i].row, best_tour[i].num); else - printf("%2u-%u (right side) ", points[best_tour[i].pt].first, - points[best_tour[i].pt].second); + printf("%2u-%u (right side) ", best_tour[i].row, best_tour[i].num); if (i == 0) { - printf("\n"); + printf(" "); } else { - unsigned cost = distance( - points[best_tour[i-1].pt].first, - points[best_tour[i-1].pt].second, - best_tour[i-1].side, - points[best_tour[i].pt].first, - points[best_tour[i].pt].second, - best_tour[i].side); - printf("cost=%4u\n", cost); + printf("cost=%4u ", best_tour[i].cost); + total_cost += best_tour[i].cost; + } + + // let's see how good the MST heuristics are + if (i != points.size() - 1) { + std::set > mst_tree = points_left; + printf("mst_bound=%5u, ", prim_mst(mst_tree)); + + unsigned rest_cost = 0; + for (unsigned j = i + 1; j < points.size(); ++j) { + rest_cost += best_tour[j].cost; + } + + printf("rest_cost=%5u", rest_cost); } + + printf("\n"); + + std::set >::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); } -unsigned do_tsp(std::vector > &points, order *ord, try_order *temp, unsigned ind, unsigned cost_so_far) +bool any_opt; +std::pair best_reorganization[KOPT_SLICES]; + +void optimize_specific_tour(order *tour, std::set > &fragments, + std::vector > &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)*history.size()); + best_so_far = cost_so_far; + } + return; + } + + // match the first point with all the others + std::set >::iterator start = fragments.begin(); + std::pair start_val = *start; + + std::set >::iterator i = start; + ++i; + + fragments.erase(start); + + while (i != fragments.end()) { + std::pair i_val = *i; + std::set >::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 > 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 > &points, order *ntour, std::set > &fragments, unsigned s, unsigned j) +{ + // find this fragment, and copy + for (std::set >::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 > &points, std::set > &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 > &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 > 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 > 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 > &points, std::set > &points_left, order *ord, order *temp, unsigned ind, unsigned cost_so_far) { if (cost_so_far >= best_so_far) return UINT_MAX; @@ -163,54 +439,60 @@ unsigned do_tsp(std::vector > &points, order *ord, 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 > 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 = points[ord[ind-1].pt].first; - unsigned last_switch = points[ord[ind-1].pt].second; - unsigned last_side = ord[ind-1].side; - for (unsigned i = 0; i < points.size(); ++i) { - /* taken earlier? */ - for (unsigned j = 0; j < ind; ++j) { - if (ord[j].pt == i) - goto taken; - } - + for (std::set >::iterator i = points_left.begin(); i != points_left.end(); ++i) { /* try both sides */ - temp[toi].o.pt = i; - temp[toi].o.side = -1; - temp[toi].cost = distance(last_row, last_switch, last_side, - points[i].first, points[i].second, -1); + temp[toi].row = i->first; + temp[toi].num = i->second; + temp[toi].side = 0; + temp[toi].cost = cache(last_row, last_switch, last_side, i->first, i->second, 0); ++toi; - temp[toi].o.pt = i; - temp[toi].o.side = +1; - temp[toi].cost = distance(last_row, last_switch, last_side, - points[i].first, points[i].second, +1); + temp[toi].row = i->first; + temp[toi].num = i->second; + temp[toi].side = 1; + temp[toi].cost = cache(last_row, last_switch, last_side, i->first, i->second, 1); ++toi; - -taken: - 1; - } - - unsigned min_rest_cost = prim_mst(points, temp, toi); - if (cost_so_far + min_rest_cost >= best_so_far) { - return UINT_MAX; } - std::sort(temp, temp + toi); unsigned best_this_cost = UINT_MAX; for (unsigned i = 0; i < toi; ++i) { - ord[ind] = temp[i].o; - unsigned cost_rest = do_tsp(points, ord, temp + points.size() * 2, ind + 1, cost_so_far + temp[i].cost); + ord[ind] = temp[i]; + + std::set >::iterator j = points_left.find(std::make_pair(temp[i].row, temp[i].num)); + points_left.erase(j); + unsigned cost_rest = do_tsp(points, points_left, ord, temp + points.size() * 2, ind + 1, cost_so_far + temp[i].cost); + points_left.insert(std::make_pair(temp[i].row, temp[i].num)); + best_this_cost = std::min(best_this_cost, cost_rest); } @@ -220,24 +502,56 @@ taken: int main() { std::vector > points; + std::set > points_left; for ( ;; ) { unsigned row, sw; 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()]; - try_order *temp = new try_order[points.size() * points.size() * 4]; + order *temp = new order[points.size() * points.size() * 4]; /* always start at the first one, left side (hack) */ - ord[0].pt = 0; - ord[0].side = -1; + ord[0].row = points[0].first; + ord[0].num = points[0].second; + ord[0].side = 0; - do_tsp(points, ord, temp, 1, 0); + do_tsp(points, points_left, ord, temp, 1, 0); printf("All done.\n"); }