nagissさんの作品(ソースコード)

//------------------------------------------------------------------------------
/// @file
/// @author   ハル研究所プログラミングコンテスト実行委員会
///
/// @copyright  Copyright (c) 2019 HAL Laboratory, Inc.
/// @attention  このファイルの利用は、同梱のREADMEにある
///             利用条件に従ってください。
//------------------------------------------------------------------------------

#include <iostream>
#include <vector>
#include <queue>
#include <algorithm>
#include <numeric>
#include <cstring>
#include <cassert>
#include <bitset>
#include <map>
#include <unordered_set>
#include <unordered_map>
#include <chrono>
#include "Answer.hpp"

using namespace std;
using PositionId = unsigned char;
constexpr int INTMAX = 2147483647;
constexpr double INF = (double)INFINITY;
constexpr unsigned long long ULLMAX = 0xffffffffffffffffull;
constexpr unsigned long long XORSHIFT_SEED = 38913499813419860ull;
constexpr int BEAM_WIDTH = 25;
constexpr int NUM_LAST_EXHAUSTIVE_SEARCH_TURNS = 3;
constexpr int NUM_ANNEALING_ITERATIONS = 23000;
constexpr int NUM_TSP_EVALUATION_ITERATION = 8;
constexpr int CHOKUDAI_WIDTH = 1;
constexpr double HEIGHT_PENALTY = 0.5;
// 他のパラメータ:
// - 焼き鈍し温度スケジュール


//------------------------------------------------------------------------------
namespace hpc {

int dist[125][125];        // エサ・カメ間の距離
Action actions[1000][25];  // 各ターンのカメの行動
Point positions[125];      // 0 ~ 99: エサの位置  100 ~ 124 カメの初期位置
int height[100];           // エサの高さ
int left_food[100];        // 巡回路の次のエサ
int right_food[100];       // 巡回路の前のエサ

vector<int> hamiltonian_cycle;  // 高いエサのハミルトン閉路
int n_turtles, n_all_foods;     // カメの数とエサの数

inline int dist_(const Point& a, const Point& b) {
	return abs(a.x - b.x) + abs(a.y - b.y);
}

struct Rmq {
	// range minimum query (平方分割)
	static const int n_backets = 9, backet_size = 10, n = 90;
	int backet_min[n_backets], data[n];
	PositionId backet_min_idxs[n_backets], idxs[n];
	Rmq() {
		fill(data, data + n, INTMAX);
		fill(backet_min, backet_min + n_backets, INTMAX);
	}
	inline void chmin(const int& i, const int& x, const PositionId& pos_id) {
		if (data[i] > x) {  // 小さくなるとき
			data[i] = x;  idxs[i] = pos_id;
			if (backet_min[i / backet_size] > x) {
				backet_min[i / backet_size] = x;  backet_min_idxs[i / backet_size] = pos_id;
			}
		}
	}
	inline PositionId get_min(const int& i) {
		const int q = i / backet_size;
		int mi = INTMAX;
		PositionId res = -1;
		for (int j = 0; j < q; j++) {
			if (mi > backet_min[j]) mi = backet_min[j], res = backet_min_idxs[j];
		}
		for (int j = q * backet_size; j < i; j++) {
			if (mi > data[j]) mi = data[j], res = idxs[j];
		}
		return res;
	}
	inline void clear() {
		fill(data, data + n, INTMAX);
		fill(backet_min, backet_min + n_backets, INTMAX);
	}
};

struct UnionFind {
	// 経路圧縮・サイズ付き Union-find 木
	vector<int> p;  // 値が正なら親ノードを指す、負ならサイズを表す
	UnionFind(const int& n) : p(n, -1) {};

	inline bool unite(int u, int v) {  // 既に同じ集合だった場合は false を返す
		if ((u = root(u)) == (v = root(v))) return false;
		if (p[u] > p[v]) swap(u, v);
		p[u] += p[v]; p[v] = u;
		return true;
	}
	inline bool same(const int& u, const int& v) { return root(u) == root(v); }
	inline int root(int u) { return p[u] < 0 ? u : p[u] = root(p[u]); }
	inline int size(const int& u) { return -p[root(u)]; }
};

template<class T> struct Tsp {
	// 巡回セールスマン問題を Held-Karp 下界を使った枝刈り全探索で解く
	// 参考: https://gist.github.com/spaghetti-source/c31558e07adcd2ced2d6
	struct Edge {  // 無向辺
		int from, to;
		T weight;
		inline Edge(const int& from, const int& to, const int& weight) : from(from), to(to), weight(weight) {}
	};

	struct State {
		vector<double> p;  // 緩和問題を解くときの各頂点のペナルティ  総和は常に 0
		vector<int> b;     // 各辺について  1: 使わないと仮定  -1: 使うと仮定  0: 探索対象
		vector<int> x;     // 緩和問題の解  各辺について使うなら 1 、使わないなら 0
		vector<int> d;     // 各頂点の次数
		double score;      // 緩和問題の解(下界)
		bool operator<(const State &s) const {
			return score > s.score;
		}
	};

	int n, m;  // 頂点の数・辺の数
	vector<Edge> edges;
	State best_state;
	vector<vector<int>> adj;

	inline Tsp(const int& n) : n(n), adj(n) { assert(n >= 3); }
	inline Tsp(const vector<int>& vertexes) : n(vertexes.size()), m(0), adj(n) {
		assert(n >= 3);
		edges.reserve(m);
		for (int i = 0; i < n - 1; i++) {
			for (int j = i + 1; j < n; j++) {
				add_edge(i, j, dist[vertexes[i]][vertexes[j]]);
			}
		}
	}

	inline void add_edge(const int& from, const int& to, const T& weight) {
		adj[from].emplace_back(edges.size());
		adj[to].emplace_back(edges.size());
		edges.emplace_back(from, to, weight);
		m++;
	}

	inline bool evaluate(State& s) {
		// 緩和問題を解いて下界を計算する
		// 下界が厳密解に等しい(閉路を得られた)とき true を返す
		// そうでなければ(閉路以外の 1-木なら)false を返す
		double eta;
		vector<int> id(m);  iota(id.begin(), id.end(), 0);  // 辺
		vector<double> modified_weights(m);  // 辺の優先順(辺について仮定があればその優先順、仮定が同じなら補正された辺の重み順)
		for (int iter = 0; iter < NUM_TSP_EVALUATION_ITERATION; iter++) {
			auto residue = [&](const int& i) {  // 補正された辺の重み
				if (s.b[i] == 1) return INF;  // 使わないと仮定した辺のコストは INF
				return edges[i].weight + s.p[edges[i].from] + s.p[edges[i].to];
			};

			// 使う順に辺をソート
			for (int i = 0; i < m; i++)
				modified_weights[i] = s.b[i] * (1 << 30) + s.p[edges[i].from] + s.p[edges[i].to] + edges[i].weight;
			sort(id.begin(), id.end(), [&](const int& i, const int& j) {
				return modified_weights[i] < modified_weights[j];
			});

			UnionFind uf(n);
			s.x.assign(m, 0);
			s.d.assign(n, 0);
			s.score = 0;
			for (const int& i : id) {  // すべての辺についてループ  // 途中で抜けても良いかも
				const Edge &e = edges[i];

				// 頂点 0 が関係する場合と関係しない場合
				if (((e.from == 0 || e.to == 0) && s.d[0] < 2) ||
					((e.from != 0 && e.to != 0) && uf.unite(e.from, e.to))) {
					s.d[e.from]++;  s.d[e.to]++;
					s.score += residue(i);
					s.x[i] = 1;
				}
			}
			eta = 1.25 * pow(0.64, iter);
			bool found = true;  // 構成した 1-木は閉路か?
			for (int i = 1; i < n; i++) {
				if (s.d[i] != 2) found = false;
				s.p[i] += eta * (s.d[i] - 2);
			}
			if (found) return true;  // 次数がすべて 2 なら閉路である(最適解が見つかっている)
		}  // end for iter;
		return false;  // 閉路にならなかった
	}  // end evaluate

	inline T solve() {
		// 枝刈り探索をする
		// 閉路を構成する辺の重みの和を返す
		if (adj[0].size() < 2) return INTMAX;
		State s = { vector<double>(n), vector<int>(m) };
		T upper_bound = INTMAX;  // 上界(暫定解)

		if (evaluate(s)) {
			best_state = s;
			return round(s.score);
		}

		priority_queue<State> que;  que.push(s);
		while (!que.empty()) {
			s = que.top();  que.pop();  // 最大の下界
			if (s.score >= upper_bound) continue;  // 上界より下界が大きい -> その枝は最適解でない

			// 構成した 1-木で次数が 2 を超えている頂点を取り出す
			int u = 1;
			for (; u < n; u++) if (s.d[u] > 2) break;

			// 使うと仮定したものに加えて、あといくつ u に繋がる辺を選べばいいか?
			int freedom = 2;
			for (const int& e : adj[u]) if (s.b[e] == -1) freedom--;

			const vector<int>& x = s.x;
			for (const int& e : adj[u]) {  // u に繋がるすべての辺 e に対して
				if (x[e] == 1 && s.b[e] == 0) {  // e が構成した 1-木に含まれていて、使うとも使わないとも仮定していない場合
					if (freedom-- > 0) {
						s.b[e] = 1;  // 使わないと仮定する
						if (evaluate(s)) {  // 閉路が構成された場合
							if (upper_bound > s.score) best_state = s, upper_bound = s.score;  // 暫定解の更新
						} else if (s.score < upper_bound) {
							que.push(s);  // 今までの暫定解より得られた下界が小さければまだわからない
						}
						s.b[e] = -1;  // 使うと仮定する(重複して同じ状態を探索することを防ぐため)
					} else {
						s.b[e] = 1;  // 使わないと仮定する(この辺を使う状態は、上で push された状態から探索される)
					}
				}
			}
			if (evaluate(s)) {
				if (upper_bound > s.score) best_state = s, upper_bound = s.score;
			} else if (s.score < upper_bound) {
				que.push(s);
			}
		}
		return round(upper_bound);
	}  // end solve

	inline vector<int> get_hamiltonian_cycle(const vector<int>& vertexes) {
		// 解から閉路を取り出す
		vector<Array<int, 2>> E(n);  // ハミルトン閉路の隣接する頂点
		for (int i = 0; i < m; i++) {  // 高いエサ間のすべての辺について
			if (best_state.x[i]) {  // ハミルトン閉路の一部なら
				const int& from = edges[i].from;
				const int& to = edges[i].to;
				E[from].add(to);
				E[to].add(from);
			}
		}
		vector<int> hamiltonian_cycle(n);
		int p = 0, pp = 0;
		hamiltonian_cycle[0] = vertexes[0];
		for (int i = 1; i < n; i++) {
			int v = E[p][0] != pp ? E[p][0] : E[p][1];
			hamiltonian_cycle[i] = vertexes[v];
			pp = p;  p = v;
		}
		return hamiltonian_cycle;
	}  // end get_hamiltonian_cycle
};  // end struct Tsp

template<class T> struct Kruskal {
	// クラスカル法
	struct Edge {  // 無向辺
		int from, to;
		T weight;
		//int used;
		Edge() {}
		Edge(const int& from, const int& to, const T& weight) : from(from), to(to), weight(weight)/*, used(0)*/ {}
		bool operator<(const Edge& e) const {
			return weight < e.weight;
		}
	};

	int n;
	Array<Edge, 400> edges;
	UnionFind uf;

	Kruskal() {}
	Kruskal(const int& n) : n(n), uf(n) {}

	inline void add_edge(const int& u, const int& v, const T& weight) {
		edges.add({ u, v, weight });
	}

	inline T solve() {
		sort(edges.begin(), edges.end());  // コストの小さい順
		T res = 0;
		for (auto& e : edges) if (uf.unite(e.from, e.to)) {
			res += e.weight;
			//e.used = 1;
		}
		return res;
	}
};

struct ManhattanMst {
	// マンハッタン距離での最小全域木
	Array<PositionId, 100> vertexes;
	Rmq rmq;
	Kruskal<int> kruskal;
	ManhattanMst(const int& n, Array<PositionId, 100>& vertexes) : vertexes(vertexes), kruskal(n) {}
	void construct_edges() {  // 各頂点から 8 方向の最近傍の頂点を見つける
		for (int b = 0; b < 4; b++) {
			rmq.clear();
			sort(vertexes.begin(), vertexes.end(), [&](const int& x, const int& y) {
				return positions[x].y != positions[y].y ? positions[x].y > positions[y].y : positions[x].x < positions[y].x;
			});  // y 座標の降順、 y 座標が同じなら x 座標の昇順
			for (PositionId& v : vertexes) {
				const int &x = positions[v].x, &y = positions[v].y;
				PositionId u = rmq.get_min(x + y + 1);  // 最も近い頂点
				if (u != (PositionId)(-1)) kruskal.add_edge(v, u, dist[v][u]);  // 同じ辺が重複することもあるが問題ないはず
				rmq.chmin(x + y, y - x, v);
			}
			if (b == 0 || b == 2) for (PositionId& v : vertexes) swap(positions[v].x, positions[v].y);
			else if (b == 1) for (PositionId& v : vertexes) positions[v].x = 29 - positions[v].x;
			else if (b == 3) for (PositionId& v : vertexes) positions[v].y = 29 - positions[v].y;
		}
	}
	inline int solve() {
		construct_edges();
		return kruskal.solve();
	}
};

inline void set_global_variables(const Stage& aStage, const bool& tsp = false) {
	const auto& turtlePositions = aStage.turtlePositions();
	const auto& foods = aStage.foods();
	n_turtles = turtlePositions.count();
	n_all_foods = foods.count();
	for (int idx_food = 0; idx_food < n_all_foods; idx_food++) {
		const Food& food = foods[idx_food];
		height[idx_food] = food.height();
		positions[idx_food] = food.pos();
	}
	for (int idx_turtle = 0; idx_turtle < n_turtles; idx_turtle++) {
		positions[idx_turtle + 100] = turtlePositions[idx_turtle];
	}
	for (int i = 0; i < 125; i++) {
		for (int j = 0; j < 125; j++) {
			dist[i][j] = dist_(positions[i], positions[j]);
		}
	}
	if (tsp) {
		// 高さがカメの数の半分より大きいエサを高いエサとして、
		// 高いエサをすべて食べて回る最短の巡回路を計算する
		vector<int> high_foods;  // 高い位置のエサ
		for (int i = 0; i < n_all_foods; i++)
			if (height[i] > n_turtles / 2)
				high_foods.push_back(i);  // 高い位置のエサを取り出す
		const int n_high_foods = high_foods.size();
		auto compare_by_dist_from_centor = [&](const int& a, const int& b) {
			return dist_(positions[a], Point(30, 15)) < dist_(positions[b], Point(30, 15));
		};
		if (n_high_foods >= 3) {
			Tsp<int> tsp(high_foods);
			tsp.solve();
			hamiltonian_cycle = tsp.get_hamiltonian_cycle(high_foods);
			// 中央に最も近いエサを vector の一番左に持ってくる  // 最終的にこの操作は必要なくなった
			auto center_food = min_element(hamiltonian_cycle.begin(), hamiltonian_cycle.end(), compare_by_dist_from_centor);
			rotate(hamiltonian_cycle.begin(), center_food, hamiltonian_cycle.end());
		} else {
			if (n_high_foods == 2 && !compare_by_dist_from_centor(high_foods[0], high_foods[1])) swap(high_foods[0], high_foods[1]);
			hamiltonian_cycle = high_foods;
		}
		fill(left_food, left_food + n_all_foods, 100);  fill(right_food, right_food + n_all_foods, 100);
		for (int i = 0; i < n_high_foods; i++) {
			right_food[hamiltonian_cycle[i]] = hamiltonian_cycle[(i + 1) % n_high_foods];
			left_food[hamiltonian_cycle[i]] = hamiltonian_cycle[(i - 1 + n_high_foods) % n_high_foods];
		}
	}
}


int turtle_values[25];  // エサを食べるカメを貪欲に選ぶときの評価値
int turtle_idxs[25];    // 評価値で nth_element するための一時的な配列

struct State {  // ビームサーチの状態
	bitset<101> food_is_eaten;               // エサを次に選んでいいかどうか  添字 100 は多分結局必要無くなった
	bitset<101> food_is_available;           // 食べてないエサについて、食べていいかどうか
	Array<PositionId, 25> turtle_positions;  // カメの位置
	Array<int, 25> turtle_turns;             // カメが最後にエサを食べたときのターン
	Array<char, 100> permutations;           // エサを食べた順
	bool closed;                             // 閉路の好きな頂点を使っていいかどうか
	int score;                               // スコア
	inline State() : closed(false) {
		for (int i = 0; i < n_turtles; i++) {
			turtle_positions.add(100 + i);
			turtle_turns.add(0);
		}
		for (int i = 0; i < n_all_foods; i++) food_is_available[i] = true;
	}
	inline State(const Array<char, 100>& permutations)
		: permutations(permutations) {
		for (int i = 0; i < n_turtles; i++) {
			turtle_positions.add(100 + i);
			turtle_turns.add(0);
		}
	}
	inline void close_cycle() {
		// 高いエサを封鎖する
		if (closed) return;
		for (const int& idx_food : hamiltonian_cycle) food_is_available[idx_food] = false;
		closed = true;
	}
	bool operator<(const State& right) const {
		return score < right.score;
	}
	bool operator>(const State& right) const {
		return score > right.score;
	}
};

inline int calc_penalty(const State& state) {
	// 残ったエサから盤面を評価する

	int l = 0, r = 0, u = 0, d = 0, remain_heights = 0;
	for (int i = 0; i < n_all_foods; i++) if (!state.food_is_eaten[i]) {
		l = min(l, positions[i].x - 30);
		r = max(r, positions[i].x - 30);
		u = min(u, positions[i].y - 15);
		d = max(d, positions[i].y - 15);
		remain_heights += height[i];
	}
	return (r - l + d - u) + (int)(remain_heights * HEIGHT_PENALTY);

	/*
	// 最小全域木の大きさ
	{
		Array<PositionId, 100> vertexes;
		for (int i = 0; i < n_all_foods; i++) if (!state.food_is_eaten[i]) {
			vertexes.add(i);
		}
		ManhattanMst mst(n_all_foods, vertexes);
		return mst.solve();
	}
	*/
}

inline void write_state(const State& best_state) {
	// state を actions に書き込む
	const Array<char, 100>& best_permutation = best_state.permutations;
	State state_reproduction;  // カメの行動を再現するための State
	for (int idx_food : best_permutation) {
		const int& h = height[idx_food];
		// カメを h 匹貪欲に選ぶ
		for (int idx_turtle = 0; idx_turtle < n_turtles; idx_turtle++) {
			turtle_values[idx_turtle]  // 一時的な配列
				= state_reproduction.turtle_turns[idx_turtle]
				+ dist[state_reproduction.turtle_positions[idx_turtle]][idx_food];
		}
		iota(turtle_idxs, turtle_idxs + n_turtles, 0);
		nth_element(turtle_idxs, turtle_idxs + (h - 1), turtle_idxs + n_turtles, [&](const int& x, const int& y) { return turtle_values[x] < turtle_values[y]; });
		// 最大ターン数を計算
		int max_turn = 0;
		for (int j = 0; j < h; j++) {
			const int& idx_turtle = turtle_idxs[j];
			const int turn = state_reproduction.turtle_turns[idx_turtle] + dist[state_reproduction.turtle_positions[idx_turtle]][idx_food];
			if (max_turn < turn) max_turn = turn;
		}
		// actions に書き込む
		for (int j = 0; j < h; j++) {
			const int& idx_turtle = turtle_idxs[j];
			Point pos = positions[state_reproduction.turtle_positions[idx_turtle]];
			for (; pos.x < positions[idx_food].x; actions[state_reproduction.turtle_turns[idx_turtle]++][idx_turtle] = Action_MoveRight, pos.x++);
			for (; pos.x > positions[idx_food].x; actions[state_reproduction.turtle_turns[idx_turtle]++][idx_turtle] = Action_MoveLeft, pos.x--);
			for (; pos.y < positions[idx_food].y; actions[state_reproduction.turtle_turns[idx_turtle]++][idx_turtle] = Action_MoveDown, pos.y++);
			for (; pos.y > positions[idx_food].y; actions[state_reproduction.turtle_turns[idx_turtle]++][idx_turtle] = Action_MoveUp, pos.y--);
			for (; state_reproduction.turtle_turns[idx_turtle] < max_turn; actions[state_reproduction.turtle_turns[idx_turtle]++][idx_turtle] = Action_Wait);
			state_reproduction.turtle_positions[idx_turtle] = idx_food;
		}
	}  // end for food
	assert(best_state.score == *max_element(state_reproduction.turtle_turns.begin(), state_reproduction.turtle_turns.end()));  // iota を使わないとこれに引っかかる
	return;
}

State beam_search(const int& beam_width) {
	iota(turtle_idxs, turtle_idxs + n_turtles, 0);
	vector<State> beam = { State() };
	for (int i = 0; i < n_all_foods; i++) {
		vector<State> new_beam;
		new_beam.reserve((n_all_foods - i) * beam.size());  // これもっと小さくても良いかも
		unordered_map<bitset<101>, int> mp;  // 同じ状態を保存しないようにして多様性を確保するための連想配列  値は new_beam のインデックス

		for (State& state : beam) {  // beam_width
			for (int idx_food = 0; idx_food < n_all_foods; idx_food++) {  // エサをひとつ選ぶ  // 100
				if (state.food_is_eaten[idx_food]) continue;
				const int& h = height[idx_food];
				// カメを h 匹貪欲に選ぶ
				for (int idx_turtle = 0; idx_turtle < n_turtles; idx_turtle++) {  // カメの評価値を設定する  // 25
					turtle_values[idx_turtle]  // 一時的な配列
						= state.turtle_turns[idx_turtle]
						+ dist[state.turtle_positions[idx_turtle]][idx_food];
				}
				iota(turtle_idxs, turtle_idxs + n_turtles, 0);
				nth_element(turtle_idxs, turtle_idxs + (h - 1), turtle_idxs + n_turtles, [&](const int& x, const int& y) { return turtle_values[x] < turtle_values[y]; });
				// 最大ターン数を計算
				int max_turn = 0;
				for (int j = 0; j < h; j++) {  // 新しい food を食べるカメ
					const int& idx_turtle = turtle_idxs[j];
					const int turn = state.turtle_turns[idx_turtle] + dist[state.turtle_positions[idx_turtle]][idx_food];
					if (max_turn < turn) max_turn = turn;  // この時点では新しい food を食べるカメ内での max であることに注意
				}
				// new_state を処理
				State new_state(state);
				new_state.permutations.add(idx_food);
				new_state.food_is_eaten.set(idx_food);
				for (int j = 0; j < h; j++) {  // 新しい food を食べるカメ
					const int& idx_turtle = turtle_idxs[j];
					new_state.turtle_positions[idx_turtle] = idx_food;
					new_state.turtle_turns[idx_turtle] = max_turn;
				}
				for (int j = h; j < n_turtles; j++) {  // 新しい food を食べないカメ
					const int& idx_turtle = turtle_idxs[j];
					if (max_turn < new_state.turtle_turns[idx_turtle]) max_turn = new_state.turtle_turns[idx_turtle];  // 全てのカメでの max をとる
				}
				int penalty = calc_penalty(new_state);
				new_state.score = max_turn + penalty;  // 小さいほど良い

				if (n_all_foods - i > NUM_LAST_EXHAUSTIVE_SEARCH_TURNS) {  // 通常
					if (mp.count(new_state.food_is_eaten) == 0) {
						mp[new_state.food_is_eaten] = new_beam.size();
						new_beam.push_back(new_state);
					} else if (new_beam[mp[new_state.food_is_eaten]].score > new_state.score) {
						new_beam[mp[new_state.food_is_eaten]] = new_state;
					}
				} else {  // 最後らへん
					new_beam.push_back(new_state);
				}
			}  // end for food
		}  // end for state
		// score が上位 beam_width 個の状態だけ残す
		int additional_beam_width = n_all_foods - i <= 0 ? BEAM_WIDTH * 8 : 0;
		if (new_beam.size() > beam_width + additional_beam_width) {
			nth_element(new_beam.begin(), new_beam.begin() + (beam_width + additional_beam_width - 1), new_beam.end());
			new_beam.resize(beam_width + additional_beam_width);
		}
		beam = new_beam;
	}  // end for idx_food

	State& best_state = *min_element(beam.begin(), beam.end());
	return best_state;
}

State chokudai_search(const int& n_iter, const int& chokudai_width = 1) {
	vector<priority_queue<State, vector<State>, greater<State>>> ques(n_all_foods + 1);
	for (int i = 0; i < n_all_foods + 1; i++) {
		vector<State> vec;  vec.reserve(n_all_foods * n_iter);
		ques[i] = priority_queue<State, vector<State>, greater<State>>(greater<State>(), move(vec));
	}

	vector<unordered_set<bitset<101>>> closed_(n_all_foods);

	ques[0].emplace();
	for (int it = 0; it < n_iter; it++) {
		for (int depth = 0; depth < n_all_foods; depth++) {
			priority_queue<State, vector<State>, greater<State>>& que = ques[depth];
			priority_queue<State, vector<State>, greater<State>>& next_que = ques[depth + 1];
			unordered_set<bitset<101>>& closed = closed_[depth];

			// 候補から最もいい感じのを選ぶ
			vector<State> chokudai;  chokudai.reserve(chokudai_width);
			for (int i = 0; i < chokudai_width && que.size(); i++) {
				State st = que.top();  que.pop();
				if (n_all_foods - depth > NUM_LAST_EXHAUSTIVE_SEARCH_TURNS) {  // 通常
					if (closed.count(st.food_is_eaten)) i--;
					else chokudai.push_back(st), closed.emplace(st.food_is_eaten);
				} else {  // 最後らへん
					chokudai.push_back(st);
				}
			}
			assert((int)chokudai.size() <= chokudai_width);

			for (const State& state : chokudai) {
				for (int idx_food = 0; idx_food < n_all_foods; idx_food++) {  // エサをひとつ選ぶ  // 100
					if (state.food_is_eaten[idx_food]) continue;
					const int& h = height[idx_food];
					// カメを h 匹貪欲に選ぶ
					for (int idx_turtle = 0; idx_turtle < n_turtles; idx_turtle++) {  // カメの評価値を設定する  // 25
						turtle_values[idx_turtle]  // 一時的な配列
							= state.turtle_turns[idx_turtle]
							+ dist[state.turtle_positions[idx_turtle]][idx_food];
					}
					iota(turtle_idxs, turtle_idxs + n_turtles, 0);
					nth_element(turtle_idxs, turtle_idxs + (h - 1), turtle_idxs + n_turtles, [&](const int& x, const int& y) { return turtle_values[x] < turtle_values[y]; });
					// 最大ターン数を計算
					int max_turn = 0;
					for (int j = 0; j < h; j++) {  // 新しい food を食べるカメ
						const int& idx_turtle = turtle_idxs[j];
						const int turn = state.turtle_turns[idx_turtle] + dist[state.turtle_positions[idx_turtle]][idx_food];
						if (max_turn < turn) max_turn = turn;  // この時点では新しい food を食べるカメ内での max であることに注意
					}
					// new_state を処理
					State new_state(state);
					new_state.permutations.add(idx_food);
					new_state.food_is_eaten.set(idx_food);
					for (int j = 0; j < h; j++) {  // 新しい food を食べるカメ
						const int& idx_turtle = turtle_idxs[j];
						new_state.turtle_positions[idx_turtle] = idx_food;
						new_state.turtle_turns[idx_turtle] = max_turn;
					}
					for (int j = h; j < n_turtles; j++) {  // 新しい food を食べないカメ
						const int& idx_turtle = turtle_idxs[j];
						if (max_turn < new_state.turtle_turns[idx_turtle]) max_turn = new_state.turtle_turns[idx_turtle];  // 全てのカメでの max をとる
					}
					int penalty = calc_penalty(new_state);
					new_state.score = max_turn + penalty;  // 小さいほど良い

					next_que.push(new_state);

				}  // end for food
			}  // end for state in chokudai
		}  // end for depth
	}
	const State& best_state = ques[n_all_foods].top();
	return best_state;
}

State beam_search_with_hamiltonian_cycle(const int& beam_width) {
	iota(turtle_idxs, turtle_idxs + n_turtles, 0);
	State init_state;

	vector<State> beam = { init_state };
	for (int i = 0; i < n_all_foods; i++) {
		vector<State> new_beam;
		new_beam.reserve((n_all_foods - i) * beam.size());  // これもっと小さくても良いかも
		unordered_map<bitset<101>, int> mp;  // 同じ状態を保存しないようにして多様性を確保するための連想配列  値は new_beam のインデックス

		for (State& state : beam) {  // beam_width
			for (int idx_food = 0; idx_food < n_all_foods; idx_food++) {  // エサをひとつ選ぶ  // 100
				if (state.food_is_eaten[idx_food]) continue;
				if (!state.food_is_available[idx_food]) continue;
				const int& h = height[idx_food];
				// カメを h 匹貪欲に選ぶ
				for (int idx_turtle = 0; idx_turtle < n_turtles; idx_turtle++) {  // カメの評価値を設定する  // 25
					turtle_values[idx_turtle]  // 一時的な配列
						= state.turtle_turns[idx_turtle]
						+ dist[state.turtle_positions[idx_turtle]][idx_food];
				}
				iota(turtle_idxs, turtle_idxs + n_turtles, 0);
				nth_element(turtle_idxs, turtle_idxs + (h - 1), turtle_idxs + n_turtles, [&](const int& x, const int& y) { return turtle_values[x] < turtle_values[y]; });
				// 最大ターン数を計算
				int max_turn = 0;
				for (int j = 0; j < h; j++) {  // 新しい food を食べるカメ
					const int& idx_turtle = turtle_idxs[j];
					const int turn = state.turtle_turns[idx_turtle] + dist[state.turtle_positions[idx_turtle]][idx_food];
					if (max_turn < turn) max_turn = turn;  // この時点では新しい food を食べるカメ内での max であることに注意
				}
				// new_state を処理
				State new_state(state);
				new_state.permutations.add(idx_food);
				new_state.food_is_eaten.set(idx_food);
				for (int j = 0; j < h; j++) {  // 新しい food を食べるカメ
					const int& idx_turtle = turtle_idxs[j];
					new_state.turtle_positions[idx_turtle] = idx_food;
					new_state.turtle_turns[idx_turtle] = max_turn;
				}
				for (int j = h; j < n_turtles; j++) {  // 新しい food を食べないカメ
					const int& idx_turtle = turtle_idxs[j];
					if (max_turn < new_state.turtle_turns[idx_turtle]) max_turn = new_state.turtle_turns[idx_turtle];  // 全てのカメでの max をとる
				}
				int penalty = calc_penalty(new_state);
				//if (idx_food == state.key_food) new_state.food_is_eaten[next_food[idx_food]] = false, new_state.key_food = next_food[idx_food];
				if (right_food[idx_food] != 100) {  // 高いエサの場合
					new_state.close_cycle();
					new_state.food_is_available[right_food[idx_food]] = true;
					new_state.food_is_available[left_food[idx_food]] = true;
				}
				new_state.score = max_turn + penalty;  // 小さいほど良い

				if (n_all_foods - i > NUM_LAST_EXHAUSTIVE_SEARCH_TURNS) {  // 通常
					if (mp.count(new_state.food_is_eaten) == 0) {
						mp[new_state.food_is_eaten] = new_beam.size();
						new_beam.push_back(new_state);
					} else if (new_beam[mp[new_state.food_is_eaten]].score > new_state.score) {
						new_beam[mp[new_state.food_is_eaten]] = new_state;
					}
				} else {  // 最後らへん
					new_beam.push_back(new_state);
				}
			}  // end for food
		}  // end for state
		// score が上位 beam_width 個の状態だけ残す
		int additional_beam_width = n_all_foods - i <= 0 ? BEAM_WIDTH * 8 : 0;
		if (new_beam.size() > beam_width + additional_beam_width) {
			nth_element(new_beam.begin(), new_beam.begin() + (beam_width + additional_beam_width - 1), new_beam.end());
			new_beam.resize(beam_width + additional_beam_width);
		}
		beam = new_beam;
	}  // end for idx_food

	State& best_state = *min_element(beam.begin(), beam.end());
	return best_state;
}

State chokudai_search_with_hamiltonian_cycle(const int& n_iter, const int& chokudai_width = 1) {
	// 高いエサはだいたい巡回路の順に食べる制約付きで chokudai サーチをする

	vector<priority_queue<State, vector<State>, greater<State>>> ques(n_all_foods + 1);
	for (int i = 0; i < n_all_foods + 1; i++) {
		// capacity を確保した priority_queue をつくる
		vector<State> vec;  vec.reserve(n_all_foods * n_iter);
		ques[i] = priority_queue<State, vector<State>, greater<State>>(greater<State>(), move(vec));
	}
	vector<unordered_set<bitset<101>>> closed_(n_all_foods);
	ques[0].emplace();
	for (int it = 0; it < n_iter; it++) {
		for (int depth = 0; depth < n_all_foods; depth++) {
			priority_queue<State, vector<State>, greater<State>>& que = ques[depth];
			priority_queue<State, vector<State>, greater<State>>& next_que = ques[depth + 1];
			unordered_set<bitset<101>>& closed = closed_[depth];

			// 候補から最もいい感じのを選ぶ
			vector<State> chokudai;  chokudai.reserve(chokudai_width);
			for (int i = 0; i < chokudai_width && que.size(); i++) {
				State st = que.top();  que.pop();
				if (n_all_foods - depth > NUM_LAST_EXHAUSTIVE_SEARCH_TURNS) {  // 通常
					if (closed.count(st.food_is_eaten)) i--;
					else chokudai.push_back(st), closed.emplace(st.food_is_eaten);
				} else {  // 最後らへん
					chokudai.push_back(st);
				}
			}
			assert((int)chokudai.size() <= chokudai_width);

			for (const State& state : chokudai) {
				for (int idx_food = 0; idx_food < n_all_foods; idx_food++) {  // エサをひとつ選ぶ  // 100
					if (state.food_is_eaten[idx_food]) continue;
					if (!state.food_is_available[idx_food]) continue;
					const int& h = height[idx_food];
					// カメを h 匹貪欲に選ぶ
					for (int idx_turtle = 0; idx_turtle < n_turtles; idx_turtle++) {  // カメの評価値を設定する  // 25
						turtle_values[idx_turtle]
							= state.turtle_turns[idx_turtle]
							+ dist[state.turtle_positions[idx_turtle]][idx_food];
					}
					iota(turtle_idxs, turtle_idxs + n_turtles, 0);  // これが無いと nth_element の結果が変わってエサ順列からカメの割り当てが再現できなくなる
					nth_element(turtle_idxs, turtle_idxs + (h - 1), turtle_idxs + n_turtles, [&](const int& x, const int& y) { return turtle_values[x] < turtle_values[y]; });
					// 最大ターン数を計算
					int max_turn = 0;
					for (int j = 0; j < h; j++) {  // 新しい food を食べるカメ
						const int& idx_turtle = turtle_idxs[j];
						const int turn = state.turtle_turns[idx_turtle] + dist[state.turtle_positions[idx_turtle]][idx_food];
						if (max_turn < turn) max_turn = turn;  // この時点では新しい food を食べるカメ内での max であることに注意
					}
					// new_state を処理
					State new_state(state);
					new_state.permutations.add(idx_food);
					new_state.food_is_eaten.set(idx_food);
					for (int j = 0; j < h; j++) {  // 新しい food を食べるカメ
						const int& idx_turtle = turtle_idxs[j];
						new_state.turtle_positions[idx_turtle] = idx_food;
						new_state.turtle_turns[idx_turtle] = max_turn;
					}
					for (int j = h; j < n_turtles; j++) {  // 新しい food を食べないカメ
						const int& idx_turtle = turtle_idxs[j];
						if (max_turn < new_state.turtle_turns[idx_turtle]) max_turn = new_state.turtle_turns[idx_turtle];  // 全てのカメでの max をとる
					}
					int penalty = calc_penalty(new_state);
					new_state.score = max_turn + penalty;  // 小さいほど良い
					if (right_food[idx_food] != 100) {  // 高いエサの場合
						// 初めて高いエサを食べた場合、他の高いエサを封鎖する
						new_state.close_cycle();
						// 巡回路の両隣の高いエサを開放する
						new_state.food_is_available[right_food[idx_food]] = true;
						new_state.food_is_available[left_food[idx_food]] = true;
					}

					next_que.push(new_state);

				}  // end for food
			}  // end for state in chokudai
		}  // end for depth
	}
	const State& best_state = ques[n_all_foods].top();
	return best_state;
}

int& evaluate(State& state, const int& upper_bound) {
	// state から permutations を取り出して評価を格納して返す

	const Array<char, 100>& permutation = state.permutations;
	iota(state.turtle_positions.begin(), state.turtle_positions.end(), 100);
	fill(state.turtle_turns.begin(), state.turtle_turns.end(), 0);
	for (const int& idx_food : permutation) {
		const int& h = height[idx_food];
		// カメを h 匹貪欲に選ぶ
		for (int idx_turtle = 0; idx_turtle < n_turtles; idx_turtle++) {
			turtle_values[idx_turtle]  // 一時的な配列
				= state.turtle_turns[idx_turtle]
				+ dist[state.turtle_positions[idx_turtle]][idx_food];
		}
		iota(turtle_idxs, turtle_idxs + n_turtles, 0);
		nth_element(turtle_idxs, turtle_idxs + (h - 1), turtle_idxs + n_turtles, [&](const int& x, const int& y) { return turtle_values[x] < turtle_values[y]; });
		// 最大ターン数を計算
		int max_turn = 0;
		for (int j = 0; j < h; j++) {
			const int& idx_turtle = turtle_idxs[j];
			const int turn = state.turtle_turns[idx_turtle] + dist[state.turtle_positions[idx_turtle]][idx_food];
			if (max_turn < turn) max_turn = turn;
		}

		// if (max_turn > upper_bound) { state.score = 1000; return state.score; }  // 焼き鈍すときには外す

		// 状態を更新
		for (int j = 0; j < h; j++) {
			const int& idx_turtle = turtle_idxs[j];
			state.turtle_turns[idx_turtle] = max_turn;
			state.turtle_positions[idx_turtle] = idx_food;
		}
	}  // end for food
	state.score = *max_element(state.turtle_turns.begin(), state.turtle_turns.end());
	return state.score;
}

unsigned long long xorshift() {
	// 64bit xorshift
	static unsigned long long y = XORSHIFT_SEED;
	y ^= y << 7;
	y ^= y >> 9;
	return y;
}

inline void rotate_left(Array<char, 100>& permutations, const int& l, const int& r) {  // 閉区間
	for (int i = l; i < r; i++) {
		swap(permutations[i], permutations[i + 1]);
	}
}

inline void rotate_right(Array<char, 100>& permutations, const int& l, const int& r) {  // 閉区間
	for (int i = r - 1; i >= l; i--) {
		swap(permutations[i], permutations[i + 1]);
	}
}

inline void reverse_(Array<char, 100>& permutations, const int& l, const int& r) {  // 閉区間
	reverse(&permutations[l], &permutations[r] + 1);
}

// 統計情報(結果には影響しない)
int stat_changed[100] = {};
int stat_diff[100] = {};

State& anneal(const int& n_iterations, State& state) {
	// エサを食べる順番を焼き鈍し
	// エサをひとつ選んでエサ順列の他の部分に挿入する

	const int type = 1;
	evaluate(state, 1000);
	//double start_temp = 0.000001, end_temp = 0.000001, temp;  // 山登り
	double start_temp = 1.75, mid_temp = 0.875, end_temp = 0.0001, temp;
	int no_change_cnt = 0;

	if (type == 1) {
		// 順列の変更箇所の候補をエサ同士の距離が近いものに絞る
		Array<PositionId, 100> pos_ids;
		for (PositionId i = 0; i < n_all_foods; i++) pos_ids.add(i);
		ManhattanMst mst(n_all_foods, pos_ids);  // エサ対を絞り込むのにマンハッタン距離 MST の前処理を流用する
		mst.construct_edges();
		const int n_edges = mst.kruskal.edges.count();
		vector<int> permutations_inv(n_all_foods);  // エサ i は何番目か?
		for (int i = 0; i < n_all_foods; i++) permutations_inv[state.permutations[i]] = i;
		for (int iter = 0; iter < n_iterations; iter++) {
			//temp = start_temp + (end_temp - start_temp) * (double)iter / (double)n_iterations;
			temp = iter > n_iterations / 2
				? (mid_temp * (double)(n_iterations - iter) / (double)(n_iterations / 2)
					+ end_temp * (double)(iter - n_iterations / 2) / (double)(n_iterations / 2))
				: (start_temp * (double)(n_iterations / 2 - iter) / (double)(n_iterations / 2)
					+ mid_temp * (double)iter / (double)(n_iterations / 2));
			
			int idx_edge = xorshift() % n_edges;
			int idx1 = permutations_inv[mst.kruskal.edges[idx_edge].from];
			int idx2 = permutations_inv[mst.kruskal.edges[idx_edge].to];
			if (idx1 > idx2) swap(idx1, idx2);
			int act = xorshift() % 2;
			if (act == 0) {  // rotate_left
				idx2 -= xorshift() % 2;
				idx2 += (int)(xorshift() % 3) - 1;  // 挿入先にランダム性を足す
			} else {  // rotate_right
				idx1 += xorshift() % 2;
				idx1 += (int)(xorshift() % 3) - 1;  // 挿入先にランダム性を足す
			}
			idx1 = min(n_all_foods - 1, max(0, idx1));
			idx2 = min(n_all_foods - 1, max(0, idx2));
			if (idx1 == idx2) continue;

			if (act == 0) rotate_left(state.permutations, idx1, idx2);
			else if (act == 1) rotate_right(state.permutations, idx1, idx2);

			int prev_score = state.score;
			if (evaluate(state, prev_score) <= prev_score) {  // 良くなった場合
				if (n_all_foods == 100) {
					stat_changed[idx1] += prev_score - state.score;
					stat_changed[idx2] += prev_score - state.score;
					stat_diff[idx2 - idx1] += prev_score - state.score;
				}
				for (int i = 0; i < n_all_foods; i++) permutations_inv[state.permutations[i]] = i;
				continue;
			} else if (xorshift() > exp((prev_score - state.score) / temp) * (double)ULLMAX) {  // 悪くなった場合確率的に状態を戻す
				if (act == 0) rotate_right(state.permutations, idx1, idx2);
				else if (act == 1) rotate_left(state.permutations, idx1, idx2);
				state.score = prev_score;
			}
		}
	} else {
		// 全てのエサ対を候補にする
		for (int iter = 0; iter < n_iterations; iter++) {
			temp = start_temp + (end_temp - start_temp) * (double)iter / (double)n_iterations;

			int idx1 = xorshift() % n_all_foods, idx2 = xorshift() % n_all_foods;
			if (idx1 == idx2) continue;
			if (idx1 > idx2) swap(idx1, idx2);
			int act = xorshift() % 2;

			if (act == 0) rotate_left(state.permutations, idx1, idx2);
			else if (act == 1) rotate_right(state.permutations, idx1, idx2);
			//else if (act == 2) swap(state.permutations[idx1], state.permutations[idx2]);
			//else reverse_(state.permutations, idx1, idx2);

			int prev_score = state.score;
			if (evaluate(state, prev_score) < prev_score) {
				//no_change_cnt = 0;  
				if (n_all_foods == 100) {
					stat_changed[idx1] += prev_score - state.score;
					stat_changed[idx2] += prev_score - state.score;
					stat_diff[idx2 - idx1] += prev_score - state.score;
				}
				continue;
			}

			double p = exp((prev_score - state.score) / temp);  // 悪くなる場合にも状態遷移する確率
			if (xorshift() > p * (double)ULLMAX) {  // 状態を戻す
				if (act == 0) rotate_right(state.permutations, idx1, idx2);
				else if (act == 1) rotate_left(state.permutations, idx1, idx2);
				//else if (act == 2) swap(state.permutations[idx1], state.permutations[idx2]);
				//else reverse_(state.permutations, idx1, idx2);
				state.score = prev_score;
				//if (++no_change_cnt >= n_all_foods*n_all_foods) break;
			}
		}
	}
	return state;
}

State anneal(const int& n_iterations) {
	State state;
	for (int i = 0; i < n_all_foods; i++) state.permutations.add(i);
	return anneal(n_iterations, state);
}


//------------------------------------------------------------------------------
/// コンストラクタ。
/// @detail 最初のステージ開始前に実行したい処理があればここに書きます。
Answer::Answer()
{
}

//------------------------------------------------------------------------------
/// デストラクタ。
/// @detail 最後のステージ終了後に実行したい処理があればここに書きます。
Answer::~Answer()
{
	//for (int i = 0; i < 100; i++) {
	//  cout << i << " changed:" << stat_changed[i] << " diff:" << stat_diff[i] << endl;
	//}
}

//chrono::system_clock::time_point t0;


//------------------------------------------------------------------------------
/// 各ステージ開始時に呼び出される処理。
/// @detail 各ステージに対する初期化処理が必要ならここに書きます。
/// @param aStage 現在のステージ。
void Answer::initialize(const Stage& aStage)
{
	//set_global_variables(aStage);
	set_global_variables(aStage, true);


	// 全てのエサでビームサーチをする
	//State st = beam_search(BEAM_WIDTH);


	// 焼きなまし
	//State st = anneal(NUM_ANNEALING_ITERATIONS);


	// 組み合わせる
	//State st = beam_search(BEAM_WIDTH);
	//anneal(NUM_ANNEALING_ITERATIONS, st);

	/*
	// TSP を解く
	Array<int, 100> permutations;
	vector<int> high_foods;  // 高い位置のエサ
	for (int i = 0; i < n_all_foods; i++) {
		if (height[i] > n_turtles / 2) high_foods.push_back(i);  // 高い位置のエサを取り出す
		else permutations.add(i);
	}
	const int n_high_foods = high_foods.size();
	if (n_high_foods >= 3) {
		Tsp<int> tsp(high_foods);
		tsp.solve();
		hamiltonian_cycle = tsp.get_hamiltonian_cycle(high_foods);
		// 中央に近いやつを前に持ってくる
		auto center_food = min_element(hamiltonian_cycle.begin(), hamiltonian_cycle.end(), [](const int& a, const int& b) {
			return dist_(positions[a], Point(30, 15)) < dist_(positions[b], Point(30, 15));
		});
		rotate(hamiltonian_cycle.begin(), center_food, hamiltonian_cycle.end());

		for (const int& idx_food : hamiltonian_cycle) permutations.add(idx_food);
	} else {
		for (const int& idx_food : high_foods) permutations.add(idx_food);
	}
	State st(n_turtles, permutations);
	evaluate(st);  // score が付いてないと assert に引っかかるので

	anneal(NUM_ANNEALING_ITERATIONS, st);
	*/


	// 制約付きビームサーチ
	//State st = beam_search_with_hamiltonian_cycle(BEAM_WIDTH);
	//anneal(NUM_ANNEALING_ITERATIONS, st);


	// chokudai サーチ
	//State st = chokudai_search(BEAM_WIDTH, CHOKUDAI_WIDTH);
	//anneal(NUM_ANNEALING_ITERATIONS, st);


	// 制約付き chokudai サーチ
	State st = chokudai_search_with_hamiltonian_cycle(BEAM_WIDTH, CHOKUDAI_WIDTH);
	anneal(NUM_ANNEALING_ITERATIONS, st);


	write_state(st);

}

//------------------------------------------------------------------------------
/// 各ターンのカメの行動を指定する。
/// @detail 各カメの行動を指定し、aActionsの要素に追加してください。
///         aActions[i]がi番目のカメの行動になります。
///         aActionsの要素数がaStage.turtlePositions()の要素数と異なる場合、アサートに失敗します。
/// @param[in] aStage 現在のステージ。
/// @param[out] aActions 各カメの行動を指定する配列。
void Answer::setActions(const Stage& aStage, Actions& aActions)
{
	const int& turn = aStage.turn();
	for (int i = 0; i < n_turtles; i++) {
		aActions.add(actions[turn][i]);
	}

	/*
	// 解答コードのサンプルです。

	// まだ食べられていないエサを1つ選んで目的地に設定する。
	Point targetFoodPosition;
	for (Food food : aStage.foods()) {
		if (!food.isEaten()) {
			targetFoodPosition = food.pos();
			break;
		}
	}

	// すべてのカメを目的地に向かって移動させる。
	for (Point turtlePosition : aStage.turtlePositions()) {
		if (turtlePosition.x < targetFoodPosition.x) {
			aActions.add(Action_MoveRight);
		} else if (turtlePosition.x > targetFoodPosition.x) {
			aActions.add(Action_MoveLeft);
		} else if (turtlePosition.y < targetFoodPosition.y) {
			aActions.add(Action_MoveDown);
		} else if (turtlePosition.y > targetFoodPosition.y) {
			aActions.add(Action_MoveUp);
		} else {
			// 移動する必要がなければ、その場で待機する。
			aActions.add(Action_Wait);
		}
	}
	*/
}

//------------------------------------------------------------------------------
/// 各ステージ終了時に呼び出される処理。
/// @detail 各ステージに対する終了処理が必要ならここに書きます。
/// @param aStage 現在のステージ。
void Answer::finalize(const Stage& aStage)
{
}

} // namespace
// EOF