X-Git-Url: https://git.sesse.net/?a=blobdiff_plain;f=tsp%2Ftsp.cpp;h=a5f36205624fb75533b5ab78a6f94938c779c01e;hb=065361daac62b4cdb33cf19caa534e29ac1dfde0;hp=760646ce687203687c633b927c6fa77892a300e5;hpb=8009a5bab6f840576b44a6d70e833c3c2be2a88b;p=nms diff --git a/tsp/tsp.cpp b/tsp/tsp.cpp index 760646c..a5f3620 100644 --- a/tsp/tsp.cpp +++ b/tsp/tsp.cpp @@ -8,9 +8,15 @@ #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]; +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, @@ -28,6 +34,14 @@ inline unsigned short &opt_cache( 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; @@ -80,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); @@ -113,14 +126,100 @@ 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); } +#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 > set2; + std::priority_queue queue; + + // 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) { @@ -153,16 +252,18 @@ int prim_mst(std::set > &set1) 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) { + printf("%2u: ", i); if (best_tour[i].side == 0) printf("%2u-%u (left side) ", best_tour[i].row, best_tour[i].num); else @@ -171,6 +272,7 @@ void print_tour(std::vector > &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 @@ -191,6 +293,141 @@ void print_tour(std::vector > &points) 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); +} + +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) @@ -202,20 +439,33 @@ unsigned do_tsp(std::vector > &points, 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 = ord[ind-1].row; - unsigned last_switch = ord[ind-1].num; - unsigned last_side = ord[ind-1].side; - - std::set > mst_set = points_left; - mst_set.insert(std::make_pair(last_row, last_switch)); for (std::set >::iterator i = points_left.begin(); i != points_left.end(); ++i) { /* try both sides */ @@ -231,12 +481,6 @@ unsigned do_tsp(std::vector > &points, std::setfirst, 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; @@ -292,9 +536,12 @@ int main() 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];