//------------------------------------------------------------------------------
/// @file
/// @author ハル研究所プログラミングコンテスト実行委員会
///
/// @copyright (C)HAL Laboratory, Inc.
/// @attention このファイルの利用は、同梱のREADMEにある
/// 利用条件に従ってください。
//------------------------------------------------------------------------------
// 本作品について (nari)
// 本作品は複数のヘッダファイルに分割して開発され，それらを1ファイルにまとめて
// 提出されています．1ファイルにまとめるスクリプトをあまり賢く書けなかったため，
// 複数回同じ記述が存在するなど若干読みづらい箇所がありますが，ご了承ください．
// (using によるタイプエイリアスや，標準ライブラリの include など．)
#include "Answer.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <vector>
// #include "local/config.hpp"
// ==============================
// local/config.hpp
// ==============================
#ifndef __CONFIG_HPP
#define __CONFIG_HPP
#include <vector>
namespace hpc {
// ==============================================
// パラメータ一覧
// =============================================
// for specific stage executor
#define SPECIFIC -1
// for corner points
#define CORNER_DELTA 0.001
#define CORNER_LESS false
// for augment points
#define AUG_STEP 0.1
#define AUG_THRESHOLD 0.5
// for graph cons
#define GRAPH_PRUNE 25.0f
// for edge
#define EDGE_MAX_LEN 80
// for apsp
int APSP_NEIGHBOR = 30;
#define APSP_MAX_LEN 300
// for atom optimizer
#define ATOM_MAX_LEN 300
bool ATOM_DEVIATE = true;
bool ATOM_STRICT = false;
bool ATOM_SKIP = true;
std::vector<bool> ATOM_CROSS({1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0});
float ATOM_CROSS_START = 0.0257152;
float ATOM_CROSS_STEP = 0.07683342;
float ATOM_CROSS_STOP = ATOM_CROSS_START + ATOM_CROSS_STEP * 9.5;
int ATOM_OPT_NUM = ATOM_CROSS.size();
std::vector<float> ATOM_LENS({0.9584536, 1.65483455, 2.2260893});
std::vector<int> ATOM_NS({16, 4, 15});
// ==============================================
// パラメータ一覧 ここまで
// =============================================
} // namespace hpc
#endif
// ==============================
// end of local/config.hpp
// ==============================
// #include "local/point.hpp"
// ==============================
// local/point.hpp
// ==============================
#ifndef __POINT_HPP
#define __POINT_HPP
#include <algorithm>
#include <cmath>
#include <vector>
namespace hpc {
// 点に関する構造の再定義と複数のメソッド定義
using Point = hpc::Vector2;
Point operator+(const Point &a, const Point &b) {
return Point(a.x + b.x, a.y + b.y);
}
Point operator-(const Point &a, const Point &b) {
return Point(a.x - b.x, a.y - b.y);
}
Point operator*(const Point &a, double b) { return Point(a.x * b, a.y * b); }
Point operator*(double b, const Point &a) { return Point(a.x * b, a.y * b); }
Point operator/(const Point &a, double b) { return Point(a.x / b, a.y / b); }
float abs(const Point &a) { return std::sqrt(a.x * a.x + a.y * a.y); }
float norm(const Point &a) { return a.x * a.x + a.y * a.y; }
bool in_stage(const Point &a) {
return 0.0 <= a.x && a.x < 50.0 && 0.0 <= a.y && a.y < 50.0;
}
Point clamp_into_stage(const Point &a) {
return Point(std::max(0.01f, std::min(49.9f, a.x)),
std::max(0.01f, std::min(49.9f, a.y)));
}
float get_nearest_distance(const std::vector<Point> &pos, const Point &p) {
float ret = 1e18;
for (const Point &q : pos) {
ret = std::min(ret, norm(p - q));
}
return std::sqrt(ret);
}
Point get_nearest_point(const std::vector<Point> &pos, const Point &p) {
float best_norm = 1e18;
Point best = pos[0];
for (const Point &q : pos) {
float v = norm(p - q);
if (v < best_norm) {
best_norm = v;
best = q;
}
}
return best;
}
} // namespace hpc
#endif
// ==============================
// end of local/point.hpp
// ==============================
// #include "local/calc_length.hpp"
// ==============================
// local/calc_length.hpp
// ==============================
#ifndef __CALC_LENGTH_HPP
#define __CALC_LENGTH_HPP
// #include "config.hpp"
// #include "point.hpp"
#include <algorithm>
#include <cmath>
#include <utility>
#include <vector>
namespace hpc {
using std::pair;
using std::vector;
using dist_type = pair<int, float>;
namespace calc_field {
Terrain terrain[50][50];
float tcoeff[50][50];
// 地形情報をコピーする
void initialize_length_field(const Stage &stage) {
for (int y = 0; y < 50; y++) {
for (int x = 0; x < 50; x++) {
terrain[y][x] = stage.terrain(Point(x, y));
tcoeff[y][x] = Parameter::JumpTerrianCoefficient[(int)terrain[y][x]];
}
}
}
// 直線的に移動した時の距離を計算する
pair<int, float> calc_length(int k, Point p, const Point &q,
int th = EDGE_MAX_LEN) {
float power = std::pow(1.1, k);
float vx = q.x - p.x;
float vy = q.y - p.y;
float vlen = std::sqrt(vx * vx + vy * vy);
vx /= vlen;
vy /= vlen;
float dist = vlen;
for (int len = 1; len <= th; len++) {
float coeff = tcoeff[(int)p.y][(int)p.x] * power;
dist -= coeff;
if (dist <= 0.0f) {
return std::make_pair(len, -dist);
}
p.x += vx * coeff;
p.y += vy * coeff;
}
return std::make_pair(250, 0.0f);
}
} // namespace calc_field
} // namespace hpc
#endif
// ==============================
// end of local/calc_length.hpp
// ==============================
// #include "local/graph.hpp"
// ==============================
// local/graph.hpp
// ==============================
#ifndef __GRAPH_HPP
#define __GRAPH_HPP
#include <utility>
#include <vector>
using std::pair;
using std::vector;
namespace hpc {
using graph_type = vector<vector<vector<int>>>;
graph_type new_graph(int n, int m, int initial_value = 0) {
return graph_type(n, vector<vector<int>>(m, vector<int>(m, initial_value)));
}
using dist_type = pair<int, float>;
using adjacent_graph_type = vector<vector<vector<pair<int, dist_type>>>>;
adjacent_graph_type new_adjacent_graph(int n, int m) {
return adjacent_graph_type(n, vector<vector<pair<int, dist_type>>>(m));
}
} // namespace hpc
#endif
// ==============================
// end of local/graph.hpp
// ==============================
// #include "local/graph_cons.hpp"
// ==============================
// local/graph_cons.hpp
// ==============================
#ifndef __GRAPH_CONS_HPP
#define __GRAPH_CONS_HPP
// #include "calc_length.hpp"
// #include "graph.hpp"
// #include "point.hpp"
#include <cmath>
#include <queue>
#include <set>
#include <utility>
#include <vector>
namespace hpc {
using std::pair;
using std::priority_queue;
using std::vector;
Terrain get_terrain(const Stage &stage, int x, int y) {
return calc_field::terrain[y][x];
}
bool diff_terrain(const Stage &stage, int ax, int ay, int bx, int by) {
return get_terrain(stage, ax, ay) != get_terrain(stage, bx, by);
}
// 隣接グラフを生成する
// 頂点の次数が大体 neighbor_size ぐらいになるようにする
// 隣接グラフは各巻物数を考慮してそれぞれ計算される
adjacent_graph_type create_neighbor_graph(const Stage &stage,
const vector<Point> &pos,
const int max_len,
const int neighbor_size) {
int n = stage.scrolls().count();
int m = (int)pos.size();
auto g = new_adjacent_graph(n, m);
for (int i = 0; i < m; i++) {
g[0][i].reserve(neighbor_size + 10);
float power = 1.0;
vector<pair<pair<Point, Point>, dist_type>> ps;
ps.reserve(m - 1);
for (int j = 0; j < m; j++) {
if (j != i && abs(pos[i] - pos[j]) < GRAPH_PRUNE) {
Point v = pos[j] - pos[i];
v = v / abs(v);
ps.emplace_back(std::make_pair(pos[i], v),
std::make_pair(j, abs(pos[j] - pos[i])));
}
}
Terrain oter = stage.terrain(pos[i]);
for (int len = 1; len <= max_len; len++) {
vector<pair<pair<Point, Point>, dist_type>> nxt;
nxt.reserve(ps.size());
for (auto P : ps) {
int j = P.second.first;
float rest = P.second.second;
Point p = P.first.first;
Point v = P.first.second;
float coeff = Parameter::JumpTerrianCoefficient[(int)stage.terrain(p)];
rest -= power * coeff;
if (rest <= 0.0f) {
g[0][i].emplace_back(j, std::make_pair(len, -rest));
} else if (len == 4 && stage.terrain(p) != oter) {
continue;
} else {
p = p + v * (coeff * power);
nxt.emplace_back(std::make_pair(p, v), std::make_pair(j, rest));
}
}
if ((int)g[0][i].size() >= neighbor_size) {
break;
}
ps.swap(nxt);
}
}
for (int k = 1; k < n; k++) {
for (int i = 0; i < m; i++) {
g[k][i].reserve(g[0][i].size());
for (auto P : g[0][i]) {
int to = P.first;
g[k][i].emplace_back(to, calc_field::calc_length(k, pos[i], pos[to]));
}
}
}
return g;
} // namespace hpc
} // namespace hpc
#endif
// ==============================
// end of local/graph_cons.hpp
// ==============================
// #include "local/augment_points.hpp"
// ==============================
// local/augment_points.hpp
// ==============================
#ifndef __AUGMENT_POINTS_HPP
#define __AUGMENT_POINTS_HPP
// #include "config.hpp"
// #include "point.hpp"
#include <algorithm>
#include <cmath>
#include <vector>
namespace hpc {
using std::vector;
// 各点から近くの線分に垂線を下ろした時の点を計算する
vector<Point> get_augment_points(const Stage &stage, const vector<Point> &pos) {
const static float dy[4] = {-1.0, 0.0, +1.0, 0.0};
const static float dx[4] = {0.0, +1.0, 0.0, -1.0};
vector<Point> ret;
for (const Point &p : pos) {
Terrain terrain = stage.terrain(p);
if (terrain == Terrain::Plain) {
continue;
}
for (int dir = 0; dir < 4; dir++) {
Terrain t = terrain;
Point q = p;
bool onemore = true;
while (true) {
q.x += dx[dir] * AUG_STEP;
q.y += dy[dir] * AUG_STEP;
if (!in_stage(q)) {
break;
}
Terrain t2 = stage.terrain(q);
if (t2 > t) {
q.x -= dx[dir] * AUG_STEP;
q.y -= dy[dir] * AUG_STEP;
if (get_nearest_distance(pos, q) > AUG_THRESHOLD) {
ret.push_back(q);
}
if (onemore) {
t = t2;
onemore = false;
} else {
break;
}
}
if (t2 < t) {
if (get_nearest_distance(pos, q) > AUG_THRESHOLD) {
ret.push_back(q);
}
t = t2;
if (!onemore) {
break;
}
}
}
}
}
return ret;
}
vector<Point> get_augment_points2(const Stage &stage,
const vector<Point> &pos) {
const static float dy[4] = {-1.0, 0.0, +1.0, 0.0};
const static float dx[4] = {0.0, +1.0, 0.0, -1.0};
vector<Point> ret;
for (const Point &p : pos) {
Terrain terrain = stage.terrain(p);
if (terrain == Terrain::Plain) {
continue;
}
for (int dir = 0; dir < 4; dir++) {
Point v(dx[dir] * AUG_STEP, dy[dir] * AUG_STEP);
Terrain t = terrain;
Point q = p;
bool onemore = true;
while (true) {
q = q + v;
if (!in_stage(q)) {
break;
}
Terrain t2 = stage.terrain(q);
if (t2 > t) {
q = q - v;
if (get_nearest_distance(pos, q) > AUG_THRESHOLD) {
ret.push_back(q);
ret.push_back(q + v);
}
if (onemore) {
t = t2;
onemore = false;
} else {
break;
}
}
if (t2 < t) {
if (get_nearest_distance(pos, q) > AUG_THRESHOLD) {
ret.push_back(q);
ret.push_back(q - v);
}
t = t2;
if (!onemore) {
break;
}
}
}
}
}
return ret;
}
} // namespace hpc
#endif
// ==============================
// end of local/augment_points.hpp
// ==============================
// #include "local/corners.hpp"
// ==============================
// local/corners.hpp
// ==============================
#ifndef __CORNERS_HPP
#define __CORNERS_HPP
// #include "point.hpp"
#include <vector>
namespace hpc {
using std::vector;
// フィールド上の「角」になる点を計算する
// 「角」としては隣接する2x2マスの中に90度の角があるところを候補とする
// 各「角」に対し，点の数を絞るため，より速い地形に点を置くようにする
vector<Point> get_corners(const Stage &stage, float delta = 0.05,
bool less = true) {
vector<Point> ret;
for (int x = 0; x < 49; x++) {
for (int y = 0; y < 49; y++) {
Point a(x + 1.0 - delta, y + 1.0 - delta);
Point b(x + 1.0 + delta, y + 1.0 - delta);
Point c(x + 1.0 + delta, y + 1.0 + delta);
Point d(x + 1.0 - delta, y + 1.0 + delta);
auto at = stage.terrain(a);
auto bt = stage.terrain(b);
auto ct = stage.terrain(c);
auto dt = stage.terrain(d);
bool corner = false;
corner |= at != bt && bt != ct;
corner |= bt != ct && ct != dt;
corner |= ct != dt && dt != at;
corner |= dt != at && at != bt;
if (!corner) {
continue;
}
for (int _ = 0; _ < 4; _++) {
// check
if (at == bt && bt == ct && ct != dt) {
if (at < dt) {
// abc:light, d:heavy
ret.push_back(b);
} else {
// abc:heavy, d:light
ret.push_back(d);
}
break;
} else if (at == bt && at != ct && at != dt && ct != dt) {
if (at < ct && at < dt) {
// ab:light, c,d:heavy
ret.push_back(ct < dt ? b : a);
} else if (at > ct && at > dt) {
// ab:heavy, c,d:light
ret.push_back(ct < dt ? c : d);
} else {
// ab:middle
if (ct < dt) {
if (!less) {
ret.push_back(b);
}
ret.push_back(c);
} else {
if (!less) {
ret.push_back(a);
}
ret.push_back(d);
}
}
break;
}
// swap
std::swap(a, b);
std::swap(b, c);
std::swap(c, d);
std::swap(at, bt);
std::swap(bt, ct);
std::swap(ct, dt);
}
}
}
return ret;
}
} // namespace hpc
#endif
// ==============================
// end of local/corners.hpp
// ==============================
// #include "local/apsp.hpp"
// ==============================
// local/apsp.hpp
// ==============================
#ifndef __APSP_HPP
#define __APSP_HPP
// #include "calc_length.hpp"
// #include "config.hpp"
// #include "graph.hpp"
// #include "point.hpp"
#include <algorithm>
#include <stack>
#include <utility>
#include <vector>
namespace hpc {
using std::vector;
using paths_type = vector<vector<vector<vector<int>>>>;
// 巻物(0 ~ n-1)と初期地点(n)の間の最短経路をそれぞれ求める
// アルゴリズムは Dijkstra 法を高速にした Dial 法を使用している
// また計算にはターン数の他に「移動時に余った距離」も考慮するようにしている
std::pair<graph_type, paths_type> calc_apsp(int n, const vector<Point> &pos,
const adjacent_graph_type &g) {
int m = pos.size();
graph_type g_out = new_graph(n, n + 1, 1000);
paths_type path_out(
n, vector<vector<vector<int>>>(n + 1, vector<vector<int>>(n + 1)));
// dijkstra
vector<int> dist(m);
vector<float> restlen(m, 0.0);
vector<bool> visited(m);
vector<int> bef(m);
std::stack<int, vector<int>> que[APSP_MAX_LEN];
for (int kk = 0; kk < n; kk++) {
for (int i = 0; i < n + 1; i++) {
if (i == n && kk != 0) {
for (int j = 0; j < n + 1; j++) {
g_out[kk][i][j] = APSP_MAX_LEN;
}
continue;
}
// using Dial's algorithm
// O(md + maxlen)
std::fill(dist.begin(), dist.end(), APSP_MAX_LEN - 1);
std::fill(restlen.begin(), restlen.end(), 0.0);
std::fill(visited.begin(), visited.end(), false);
std::fill(bef.begin(), bef.end(), -1);
for (int q = 0; q < APSP_MAX_LEN; q++) {
while (que[q].size()) {
que[q].pop();
}
}
dist[i] = 0;
que[0].push(i);
for (int dial = 0; dial < APSP_MAX_LEN; dial++) {
if (que[dial].size() == 0) {
continue;
}
while (que[dial].size()) {
int j = que[dial].top();
que[dial].pop();
if (visited[j]) {
continue;
}
visited[j] = true;
for (auto P : g[kk][j]) {
int to = P.first;
auto dd = P.second;
int nd = dist[j] + dd.first;
float nr = restlen[j] + dd.second;
if (nd < dist[to] || (nd == dist[to] && nr > restlen[to])) {
dist[to] = nd;
restlen[to] = nr;
bef[to] = j;
que[nd].push(to);
}
}
}
// early break
int visited_count = 0;
for (int j = 0; j < n; j++) {
visited_count += visited[j];
}
if (visited_count == n || visited_count >= n*4/5 + 1) {
break;
}
}
for (int j = 0; j < n + 1; j++) {
g_out[kk][i][j] = dist[j];
auto &v = path_out[kk][i][j];
v.clear();
if (dist[j] >= APSP_MAX_LEN - 1) {
continue;
}
int last = j;
v.push_back(last);
while (last != i) {
last = bef[last];
v.push_back(last);
}
std::reverse(v.begin(), v.end());
}
}
}
return std::make_pair(g_out, path_out);
}
} // namespace hpc
#endif
// ==============================
// end of local/apsp.hpp
// ==============================
// #include "local/atom_optimizer.hpp"
// ==============================
// local/atom_optimizer.hpp
// ==============================
#ifndef __ATOM_OPTIMIZER_HPP
#define __ATOM_OPTIMIZER_HPP
// #include "calc_length.hpp"
// #include "config.hpp"
// #include "point.hpp"
#include <cmath>
#include <stack>
#include <vector>
namespace hpc {
using std::pair;
using std::vector;
// atom optimizer: 経路の点を少し揺らして最短経路が更新されるかを試す
// corner_opt フラグが立ってる場合は十文字に，そうでない場合は円形に揺らす
// 命名は揺らし方を原子(atom)の電子配置から思いついたため
vector<Point> atom_optimize(const Stage &stage, const vector<Point> &origin,
float deviate = 0.0, bool corner_opt = false,
int round = 0) {
// construct delta vector
auto rot = [&](int j, int k, int n) {
float th = 2 * M_PI / n * k + (ATOM_DEVIATE ? deviate * j : 0.0f);
return Point(std::cos(th), std::sin(th));
};
vector<Point> delta;
delta.push_back(Point(0, 0));
if (!corner_opt) {
for (int i = 0; i < ATOM_NS[round % ATOM_NS.size()]; i++) {
delta.push_back(ATOM_LENS[round % ATOM_LENS.size()] *
rot(round + 1, i, ATOM_NS[round % ATOM_NS.size()]));
}
} else {
for (float len = ATOM_CROSS_START; len < ATOM_CROSS_STOP;
len += ATOM_CROSS_STEP) {
for (int i = 0; i < 4; i++) {
delta.push_back(len * rot(0, i, 4));
}
}
}
int n = stage.scrolls().count();
int m = origin.size();
int u = delta.size();
using calc_field::calc_length;
int scroll_count = 0;
vector<bool> gotten(n, false);
vector<Point> cands;
vector<bool> valid;
vector<int> cands_count;
vector<int> scroll_ids;
for (int pi = 0; pi < (int)origin.size(); pi++) {
const Point &p = origin[pi];
int scroll_id = -1;
for (int i = 0; i < n; i++) {
if (!gotten[i] &&
static_cast<int>(stage.scrolls()[i].pos().x) ==
static_cast<int>(p.x) &&
static_cast<int>(stage.scrolls()[i].pos().y) ==
static_cast<int>(p.y)) {
gotten[i] = true;
scroll_count++;
scroll_id = i;
}
}
Terrain oter = stage.terrain(p);
float power = std::pow(1.1, scroll_count);
for (const Point &d : delta) {
auto q = clamp_into_stage(p + power * d);
cands.push_back(q);
cands_count.push_back(scroll_count);
scroll_ids.push_back(
scroll_id != -1 &&
static_cast<int>(stage.scrolls()[scroll_id].pos().x) ==
static_cast<int>(q.x) &&
static_cast<int>(stage.scrolls()[scroll_id].pos().y) ==
static_cast<int>(q.y)
? scroll_id
: -1);
valid.push_back(oter == stage.terrain(q));
}
}
vector<vector<pair<int, pair<int, float>>>> g(cands.size());
for (int i = 0; i < (int)cands.size(); i++) {
if (!valid[i] && ATOM_STRICT) {
continue;
}
const Point &p = cands[i];
int k = cands_count[i];
int b = i / u;
for (int j = 0; j < u; j++) {
int to = u * b + j;
if (!valid[to] && ATOM_STRICT) {
continue;
}
auto len = calc_length(k, p, cands[to]);
if (len.first > 0) {
g[i].emplace_back(to, len);
}
}
for (int j = 0; j < u; j++) {
int to = u * (b + 1) + j;
if (!valid[to] && ATOM_STRICT) {
continue;
}
if (to < (int)cands.size()) {
if (k == cands_count[to] ||
(k + 1 == cands_count[to] && scroll_ids[to] != -1)) {
g[i].emplace_back(to, calc_length(k, p, cands[to]));
}
}
}
// skipper
if (ATOM_SKIP) {
for (int j = 0; j < u; j++) {
int to = u * (b + 2) + j;
if (!valid[to] && ATOM_STRICT) {
continue;
}
if (to < (int)cands.size()) {
if (k == cands_count[to] ||
(k + 1 == cands_count[to] && scroll_ids[to] != -1)) {
g[i].emplace_back(to, calc_length(k, p, cands[to]));
}
}
}
}
}
// using Dial's algorithm
// O(md + maxlen)
m = cands.size();
vector<int> dist(m, ATOM_MAX_LEN - 1);
vector<float> restlen(m, 0.0);
vector<bool> visited(m, false);
vector<int> bef(m, -1);
std::stack<int, vector<int>> que[ATOM_MAX_LEN];
dist[0] = 0;
que[0].push(0);
for (int dial = 0; dial < ATOM_MAX_LEN; dial++) {
if (que[dial].size() == 0) {
continue;
}
bool goal = false;
int goal_pos = -1;
float best_nr = 0.0;
while (que[dial].size()) {
int j = que[dial].top();
que[dial].pop();
if (visited[j]) {
continue;
}
visited[j] = true;
if (cands_count[j] == n) {
// goal
goal = true;
if (restlen[j] > best_nr) {
goal_pos = j;
best_nr = restlen[j];
}
}
for (auto nei : g[j]) {
int to = nei.first;
int nd = dist[j] + nei.second.first;
float nr = restlen[j];
// if (scroll_ids[to] != -1) {
nr += nei.second.second;
// }
if (nd < dist[to] || (nd == dist[to] && nr > restlen[to])) {
dist[to] = nd;
restlen[to] = nr;
bef[to] = j;
que[nd].push(to);
}
}
}
if (goal) {
vector<Point> ret;
int last = goal_pos;
ret.push_back(cands[last]);
while (last != 0) {
last = bef[last];
if (abs(ret.back() - cands[last]) > 1e-8) {
ret.push_back(cands[last]);
}
}
std::reverse(ret.begin(), ret.end());
return ret;
}
}
// fail
return origin;
}
} // namespace hpc
#endif
// ==============================
// end of local/atom_optimizer.hpp
// ==============================
// #include "local/tsp.hpp"
// ==============================
// local/tsp.hpp
// ==============================
#ifndef __TSP_HPP
#define __TSP_HPP
// #include "graph.hpp"
#include <algorithm>
#include <cassert>
#include <fstream>
#include <iomanip>
#include <numeric>
#include <random>
#include <sstream>
#include <vector>
namespace hpc {
int calc_score(int n, const graph_type &g, const std::vector<int> &order) {
assert(order[0] == n);
int ret = 0;
for (int k = 0; k < n; k++) {
int p = order[k];
int q = order[k + 1];
ret += g[k][p][q];
}
return ret;
}
// 最近傍法によるTSP計算
std::vector<int> nearest_neighbor(int n, const graph_type &g) {
std::vector<bool> visited(n, false);
std::vector<int> order(n + 1);
order[0] = n;
for (int k = 0; k < n; k++) {
int bef = order[k];
int to = -1;
for (int i = 0; i < n; i++) {
if (visited[i]) {
continue;
}
if (to == -1 || g[k][bef][i] < g[k][bef][to]) {
to = i;
}
}
order[k + 1] = to;
visited[to] = true;
}
return order;
}
// 3-opt 法によるTSP計算
std::vector<int> three_opt(int n, const graph_type &g, std::vector<int> order) {
int score = calc_score(n, g, order);
while (true) {
bool updated = false;
for (int i = 1; i <= n; i++) {
for (int j = i + 1; j <= n + 1; j++) {
for (int k = j; k <= n + 1; k++) {
auto tmp = order;
int tmp_score = 0;
// l
std::reverse(tmp.begin() + i, tmp.begin() + j);
tmp_score = calc_score(n, g, tmp);
if (tmp_score < score) {
score = tmp_score;
order = tmp;
updated = true;
}
// r
std::reverse(tmp.begin() + j, tmp.begin() + k);
tmp_score = calc_score(n, g, tmp);
if (tmp_score < score) {
score = tmp_score;
order = tmp;
updated = true;
}
// lr
std::reverse(tmp.begin() + i, tmp.begin() + k);
tmp_score = calc_score(n, g, tmp);
if (tmp_score < score) {
score = tmp_score;
order = tmp;
updated = true;
}
// l
std::reverse(tmp.begin() + i, tmp.begin() + j);
tmp_score = calc_score(n, g, tmp);
if (tmp_score < score) {
score = tmp_score;
order = tmp;
updated = true;
}
// r
std::reverse(tmp.begin() + j, tmp.begin() + k);
tmp_score = calc_score(n, g, tmp);
if (tmp_score < score) {
score = tmp_score;
order = tmp;
updated = true;
}
// l
std::reverse(tmp.begin() + i, tmp.begin() + j);
tmp_score = calc_score(n, g, tmp);
if (tmp_score < score) {
score = tmp_score;
order = tmp;
updated = true;
}
// lr
std::reverse(tmp.begin() + i, tmp.begin() + k);
tmp_score = calc_score(n, g, tmp);
if (tmp_score < score) {
score = tmp_score;
order = tmp;
updated = true;
}
}
}
}
if (!updated) {
break;
}
}
return order;
}
// 局所最適解からの脱出用メソッド(キック)
template <typename RNG>
std::vector<int> bridge_kick(int n, const graph_type &g, std::vector<int> order,
RNG &rng) {
for (int _ = 0; _ < 2; _++) {
int i = rng() % n + 1;
int j = rng() % (n - 1) + 1;
if (j >= i) {
j++;
} else {
std::swap(i, j);
}
std::reverse(order.begin() + i, order.begin() + j + 1);
}
return order;
}
// 最近傍法・3-opt・キックによるTSP計算
std::vector<int> nn_three_opt(int n, const graph_type &g) {
assert((int)g.size() == n && (int)g[0].size() == n + 1 &&
(int)g[0][0].size() == n + 1);
if (n <= 8) {
std::vector<int> order(n + 1);
order[0] = n;
std::iota(order.begin() + 1, order.end(), 0);
std::vector<int> best_order = order;
int best_score = 1000;
do {
int score = calc_score(n, g, order);
if (score < best_score) {
best_order = order;
best_score = score;
}
} while (std::next_permutation(order.begin() + 1, order.end()));
return best_order;
}
// 1. Nearest Neighbor
auto order = nearest_neighbor(n, g);
// 2. 3-opt
order = three_opt(n, g, order);
// 3. kick and 3-opt
auto best_order = order;
int best_score = calc_score(n, g, order);
std::mt19937 rng(252521);
for (int i = 0; i < 30; i++) {
order = bridge_kick(n, g, order, rng);
order = three_opt(n, g, order);
int score = calc_score(n, g, order);
if (score < best_score) {
best_score = score;
best_order = order;
}
}
return best_order;
}
} // namespace hpc
#endif
// ==============================
// end of local/tsp.hpp
// ==============================
namespace hpc {
//------------------------------------------------------------------------------
// constructor / destructor
int stage_count;
int stage_turn;
int total_turn;
Answer::Answer() {
stage_count = 0;
stage_turn = 0;
total_turn = 0;
}
Answer::~Answer() { fprintf(stderr, "total turn: %d\n", total_turn); }
//------------------------------------------------------------------------------
// stage begin
// 移動経路を記述する点のリストと，その順序を格納する配列
vector<Point> pos;
vector<int> order;
// std::vector<Vector2> Answer::initialize(const Stage &aStage) {
void Answer::initialize(const Stage &aStage) {
// 初期化
order.clear();
pos.clear();
int n = aStage.scrolls().count();
calc_field::initialize_length_field(aStage);
// 候補点群の計算(巻物，初期地点，角，補助点)
for (int i = 0; i < n; i++) {
pos.push_back(aStage.scrolls()[i].pos());
}
pos.push_back(aStage.rabbit().pos());
auto corners = get_corners(aStage, CORNER_DELTA, CORNER_LESS);
pos.insert(pos.end(), corners.begin(), corners.end());
auto augments = get_augment_points2(aStage, pos);
pos.insert(pos.end(), augments.begin(), augments.end());
// デバッグ用: SPECIFIC に指定されたステージ以外では切り上げ
if (stage_count != SPECIFIC && SPECIFIC >= 0) {
order.clear();
for (int i = 0; i < n; i++) {
order.push_back(i);
}
// return pos;
return;
}
// 候補点群を線分で結ぶ有向グラフの構築
int m = pos.size();
auto g = create_neighbor_graph(aStage, pos, APSP_MAX_LEN, APSP_NEIGHBOR);
// 有向グラフの逆辺の追加
vector<vector<int>> rg(m);
for (int i = 0; i < m; i++) {
for (auto P : g[0][i]) {
rg[P.first].push_back(i);
}
}
for (int k = 0; k < n; k++) {
for (int i = 0; i < m; i++) {
g[k][i].reserve(g[k][i].size() + rg[i].size());
for (int j : rg[i]) {
g[k][i].emplace_back(j, calc_field::calc_length(k, pos[i], pos[j]));
}
}
}
// 有向グラフにおける巻物間の最短経路計算
// 最短経路長の表とその際に通る点群のリストを返す
auto apsp_out = calc_apsp(n, pos, g);
auto g2 = apsp_out.first;
auto path = apsp_out.second;
// 巡回セールスマン問題による最短経路計算
// 初期地点から始まってすべての巻物を回収する最短経路
auto tsp_order = nn_three_opt(n, g2);
// 実際に通る点群を計算
order.clear();
order.push_back(n);
for (int kk = 0; kk < n; kk++) {
int from = tsp_order[kk];
int to = tsp_order[kk + 1];
auto &v = path[kk][from][to];
for (int i = 1; i < (int)v.size(); i++) {
order.push_back(v[i]);
}
}
// 経路の最適化
vector<Point> points;
for (int i : order) {
points.push_back(pos[i]);
}
for (int i = 0; i < ATOM_OPT_NUM; i++) {
points = atom_optimize(aStage, points, i, ATOM_CROSS[i], i);
}
// 整理
pos = points;
order.clear();
for (int i = 0; i < (int)points.size(); i++) {
order.push_back(i);
}
// スタック構造として使用するために順序を逆転させる
std::reverse(order.begin(), order.end());
// return pos;
return;
}
Vector2 Answer::getTargetPos(const Stage &aStage) {
stage_turn++;
total_turn++;
// エラー
if (order.size() == 0) {
fprintf(stderr, "order is null!\n");
return aStage.rabbit().pos();
}
// 指定された点に向けて移動する
auto &p = aStage.rabbit().pos();
auto &q = pos[order.back()];
// 指定された点に到達していたらスタックからポップする
if (abs(p - q) < 1e-8) {
order.pop_back();
}
return pos[order.back()];
}
void Answer::finalize(const Stage &aStage) {
if (SPECIFIC == -1 || SPECIFIC == stage_count) {
fprintf(stderr, "stage %3d: %3d\n", stage_count, stage_turn);
}
stage_count++;
stage_turn = 0;
}
} // namespace hpc
// EOF