2018 TCO Marathon Match Round 4 - CakeSharing
千のバグの上に立つ
provisional score: 555108.81
provisional rank: 40th / 86
2018TCOMMR4
CakeSharing
制約
10 <= H, W <= 80
10 <= prob <= 90
5 <= NP <= 2 * (H + W)
格子点上の多角形の面積：ピックの定理
面積の計算を三角形の足し合わせでやるなら，線が折れている部分の点の情報をメモしておく必要がありそう
折れ頂点(+ 格子点上の頂点？)，rose，area などを格納した構造体 polygon を作る
時計回りとかで格納しないと駄目？(vis ではそうなっている)
凹多角形になっている場合などは 2 点選んでも分割できない？→問題の制約上凸多角形しかできないので大丈夫
valid な 2 点を選ぶと領域が分割される
S = area / numPieces
R = roses / numPieces
として，分割時に片方の領域が area S, rose R に近くなるよう点を選んでいく？
残りの area, rose にも気を配って帳尻合わせする必要がありそう
線上で死ぬ rose も考えると R <= roses / numPieces とするのが良いか
当然順序に依存する
順序依存性のある問題はビームサーチという偏見があるが，頑張って焼けないかも随時考える
まずは貪欲とか書いて投げるか
isOnPieceEdge が線形探索してるの無駄かなと思ったけど早すぎる高速化はやめたほうがいい
まずはナイーブに動くコードを書くことが大事
isRoseInsidePiece の挙動が謎
rose の中心を整数にするために座標を 2 倍してるのは理解できるが…
領域分割の状態を木構造で表現すると良さそう
area と numRoses からいい感じに"あといくつの領域に分ける必要があるか"を決める
分割した時に cut 上に格子点が少ないとその後の分割が難しくなる(最も良いのは vertical or horizontal( or diagonal) な分割)ので，木の root に近いほうではそれに応じたペナルティを課したほうが良さそう
実現可能そうな方法として
・領域は木構造
・カットの依存関係は DAG
で表現するというのが思いついた
子を持たないノード(＝端点以外から別の直線が生えていない直線)は除去することができる
直線を除去した際にどの領域をマージするかが分からないと困る
プロファイラによると isRoseInsidePiece にほぼ時間が掛かっている
曲がり角だけ調べればいい
vertical or horizontal なものを全部探索して最もいい分割を採用する貪欲(激遅)を投げた
submission 1:
238407.80 (42)
大量に TLE 吐いてるはずなのにまあまあ点が出てるってことは貪欲ちゃんと書ければ 500k とか出る気がする
今回得点計算えぐすぎるし 500k でもまあまあの位置につけそうなんだよね
countRosesInsidePiece を高速化したい
rose が piece に何個含まれているかを高速に求める
rose が属する piece の idx をメモしておく？
凸包の y 座標ごとの x 区間を求めて累積和？
いずれも cut した際に計算する
凸包の y 座標ごとの x 区間を求めて累積和を求めるコードを書いた
なんかバグってるけど投げる
submission 2:
238102.07(42) -> 197224.48(43)
まあバグってるし
まじめに計算をやってバグが取れたので投げる
submission 3:
197215.81(45) -> 276301.57(45)
点は上がれど順位上がらず
いくつかの処理は凸包の頂点のみ着目すれば良いような気がするのでそれを格納する構造体 Piece を作ろう
ようやく線形探索してる isOnPieceEdge がボトルネックになったので改善する
・point p がどの piece に属しているか (複数ありえる)
・ある piece に点 p が含まれるか (set とかで)
set にしたら逆に遅くなった
状態のコピーが重いので Undo を実装する必要が出てきた
元のピースを構成する点とか idx とか持たせるならやはり piece 構造体を作ってしまうのが正しい気がする
--- Piece 構造体に格納するべき情報
piece の idx (new!, 一応)
格子点 (元 pieces)
頂点 (new!)
area
rosesInside
numDivides (あといくつの領域に分割する必要があるか)
分割する前の piece の idx (new!)
分割する前の piece 本体 (new!)
木構造にするならもっと真面目にやらないとダメそうだけど今はこれくらいで
とりあえずリファクタリングの前に submit する
貪欲解生成時の評価関数をほぼ生のスコアから (stdDevArea + stdDevRoses) に変えたらスコアが伸びそうなのでそれで提出
submission 4:
264348.51(46) -> 349042.77(37)
good ですね
roses, cumuRoses などはコピーコストがでかいのでグローバルに置くべきか？
というか現在 cumuRoses さえあれば凸包内の roses を求めることができるので roses は完全にお払い箱だった
countRosesInsidePiece -> countRosesInsideConvexHull
pts 系は後で state の外に出してよさそう
undo をなんとか実装
無限にバグってしんどかった 実装力なさすぎる
次は vertex の実装
遅かった isOnPieceEdge などは vertex を結ぶ直線上に存在するか計算することで早くなりそう
あとで getConvexHullArea を vertex に変えて validation
バカなので makeCut に isOnConvexHullEdge がいらないことに気付かなかった
分割する piece を決めて頂点を選んでいるのだから piece id と point id さえ分かればカット可能じゃん
ああ
emplace 系関数への移行
vertex だけで十分なんじゃ？
prev 保持しとくのもなんかダサい
vertex のみのバージョンができた
領域内の roses 判定も map 使わないことで最大ケースでも 2 秒とかで全探索貪欲解を得られるようになった
ようやく焼きなましとか MM っぽいことができそう
現在 undo はスタックに突っ込んだ piece 構造体をまるごと戻すことで行っているが，これでは他の領域に変更が生じた場合に不具合が起きる
cut の情報のみから pieces をマージできるよう実装を変更する必要がある
R
0 ---------- 1
| |
|U |D
| |
| L |
3 ---------- 2
cut: {}
piece: {{0, 1, 2, 3, 0}}
R R
0 --- 4 ---- 1
| | |
|U D|U |D
| | |
| | |
3 --- 5 ---- 2
L L
cut: {(4, 5)}
piece: {{4, 5, 3, 0, 4}, {5, 4, 1, 2, 5}}
4 の位置がわかれば
{4, 5, 3, 0, 4} の 4, 5 以降の 3, 0
{5, 4, 1, 2, 5} の 5, 4 以降の 1, 2
をマージして
{3, 0, 1, 2, 3}
という元の piece を復元できる
cut がどの piece に属するか瞬時に求めるには
最初考えていた木構造で pieces・有向グラフで cut の依存関係を保持する方法が有効か
少なくともノード数のオーダーで所属領域は求まる
cut が除去可能かは O(1)
ちょっとムズそうなのでもう少し簡単に
辺の削除
・辺を explicit に共有する 2 領域を見つける
・マージ
マージを簡単にするために makeCut の piece 生成規則を統一する
木構造はちょっとむずかしいので放置
ピースのマージはなんとか実装した
特定の辺を除去することを考えると，その辺の子以降のすべての辺を除去する必要がある
再帰で書けそうだが…
しんどすぎる 助けてくれ
終了 3 時間前
countRosesInsideConvexHull がバグってた
直したらかなりスコアが上がりそう
dependency.removable が空になる問題
依存グラフの更新がどこかバグってる？
焼きなましまで行くにはあまりにも実装力がなさすぎた
結局多点スタート＋乱択でごまかすことに
#define _CRT_SECURE_NO_WARNINGS #include "bits/stdc++.h" #include <random> #include <unordered_map> #include <unordered_set> #define LOCAL #define OFFICIAL #define VISUALIZE #ifdef VISUALIZE #include <opencv2/core.hpp> #include <opencv2/highgui.hpp> #include <opencv2/imgproc.hpp> #endif // LOCAL using namespace std; //呪文 #define DUMPOUT cerr #define dump(...) DUMPOUT<<" ";DUMPOUT<<#__VA_ARGS__<<" :["<<__LINE__<<":"<<__FUNCTION__<<"]"<<endl;DUMPOUT<<" ";dump_func(__VA_ARGS__) typedef unsigned uint; typedef long long ll; typedef unsigned long long ull; typedef pair<ll, ll> pll; typedef pair<double, double> pdd; typedef pair<string, string> pss; template <typename _KTy, typename _Ty> ostream& operator << (ostream& ostr, const pair<_KTy, _Ty>& m) { ostr << "{" << m.first << ", " << m.second << "}"; return ostr; } template <typename _KTy, typename _Ty> ostream& operator << (ostream& ostr, const map<_KTy, _Ty>& m) { if (m.empty()) { ostr << "{ }"; return ostr; } ostr << "{" << *m.begin(); for (auto itr = ++m.begin(); itr != m.end(); itr++) { ostr << ", " << *itr; } ostr << "}"; return ostr; } template <typename _KTy, typename _Ty> ostream& operator << (ostream& ostr, const unordered_map<_KTy, _Ty>& m) { if (m.empty()) { ostr << "{ }"; return ostr; } ostr << "{" << *m.begin(); for (auto itr = ++m.begin(); itr != m.end(); itr++) { ostr << ", " << *itr; } ostr << "}"; return ostr; } template <typename _Ty> ostream& operator << (ostream& ostr, const vector<_Ty>& v) { if (v.empty()) { ostr << "{ }"; return ostr; } ostr << "{" << v.front(); for (auto itr = ++v.begin(); itr != v.end(); itr++) { ostr << ", " << *itr; } ostr << "}"; return ostr; } template <typename _Ty> ostream& operator << (ostream& ostr, const set<_Ty>& s) { if (s.empty()) { ostr << "{ }"; return ostr; } ostr << "{" << *(s.begin()); for (auto itr = ++s.begin(); itr != s.end(); itr++) { ostr << ", " << *itr; } ostr << "}"; return ostr; } template <typename _Ty> ostream& operator << (ostream& ostr, const unordered_set<_Ty>& s) { if (s.empty()) { ostr << "{ }"; return ostr; } ostr << "{" << *(s.begin()); for (auto itr = ++s.begin(); itr != s.end(); itr++) { ostr << ", " << *itr; } ostr << "}"; return ostr; } template <typename _Ty> ostream& operator << (ostream& ostr, const stack<_Ty>& s) { if (s.empty()) { ostr << "{ }"; return ostr; } stack<_Ty> t(s); ostr << "{" << t.top(); t.pop(); while (!t.empty()) { ostr << ", " << t.top(); t.pop(); } ostr << "}"; return ostr; } template <typename _Ty> ostream& operator << (ostream& ostr, const list<_Ty>& l) { if (l.empty()) { ostr << "{ }"; return ostr; } ostr << "{" << l.front(); for (auto itr = ++l.begin(); itr != l.end(); ++itr) { ostr << ", " << *itr; } ostr << "}"; return ostr; } template <typename _KTy, typename _Ty> istream& operator >> (istream& istr, pair<_KTy, _Ty>& m) { istr >> m.first >> m.second; return istr; } template <typename _Ty> istream& operator >> (istream& istr, vector<_Ty>& v) { for (size_t i = 0; i < v.size(); i++) istr >> v[i]; return istr; } namespace aux { // print tuple template<typename Ty, unsigned N, unsigned L> struct tp { static void print(ostream& os, const Ty& v) { os << get<N>(v) << ", "; tp<Ty, N + 1, L>::print(os, v); } }; template<typename Ty, unsigned N> struct tp<Ty, N, N> { static void print(ostream& os, const Ty& value) { os << get<N>(value); } }; } template<typename... Tys> ostream& operator<<(ostream& os, const tuple<Tys...>& t) { os << "{"; aux::tp<tuple<Tys...>, 0, sizeof...(Tys)-1>::print(os, t); os << "}"; return os; } template<typename A, size_t N, typename T> inline void Fill(A(&array)[N], const T &val) { std::fill((T*)array, (T*)(array + N), val); } void dump_func() { DUMPOUT << endl; } template <class Head, class... Tail> void dump_func(Head&& head, Tail&&... tail) { DUMPOUT << head; if (sizeof...(Tail) == 0) { DUMPOUT << " "; } else { DUMPOUT << ", "; } dump_func(std::move(tail)...); } #define PI 3.14159265358979323846 #define EPS 1e-10 #define FOR(i,a,b) for(int i=(a);i<(b);++i) #define REP(i,n) FOR(i,0,n) #define all(x) (x).begin(), (x).end() #define SZ(x) ((int)(x).size()) #define fake false #define LE(n,m) ((n) < (m) + EPS) #define GE(n,m) ((n) + EPS > (m)) #define EQ(n,m) (abs((n)-(m)) < EPS) #define X(p) (p.first) #define Y(p) (p.second) const bool debug_mode = false; double time_elapses; class timer { vector<timer> timers; int n = 0; public: #ifdef LOCAL double limit = 9.85; #else double limit = 9.85; #endif double t = 0; timer() {} timer(int size) : timers(size) {} bool elapses() const { return time() - t > limit; } void measure() { t = time() - t; ++n; } void measure(char id) { timers[id].measure(); } void print() { if (n % 2) measure(); for (int i = 0; i < 128; ++i) { if (timers[i].n) cerr << (char)i << ' ' << timers[i].t << 's' << endl; } cerr << " " << t << 's' << endl; } static double time() { #ifdef LOCAL return __rdtsc() / 2.6e9; #else unsigned long long a, d; __asm__ volatile ("rdtsc" : "=a" (a), "=d" (d)); return (d << 32 | a) / 2.8e9; #endif } } timer(128); struct SXor128 { unsigned x, y, z, w; SXor128() { x = 123456789; y = 362436069; z = 521288629; w = 88675123; } SXor128(int _w) { x = 123456789; y = 362436069; z = 521288629; w = _w; } void setSeed() { x = 123456789; y = 362436069; z = 521288629; w = 88675123; } void setSeed(int _w) { x = 123456789; y = 362436069; z = 521288629; w = _w; } unsigned nextUInt() { unsigned t = (x ^ (x << 11)); x = y; y = z; z = w; return (w = (w ^ (w >> 19)) ^ (t ^ (t >> 8))); } unsigned nextUInt(unsigned mod) { unsigned t = (x ^ (x << 11)); x = y; y = z; z = w; w = (w ^ (w >> 19)) ^ (t ^ (t >> 8)); return w % mod; } // [l, r] unsigned nextUInt(unsigned l, unsigned r) { unsigned t = (x ^ (x << 11)); x = y; y = z; z = w; w = (w ^ (w >> 19)) ^ (t ^ (t >> 8)); return w % (r - l + 1) + l; } double nextDouble() { return double(nextUInt()) / UINT_MAX; } } rnd; class TestCase { vector<string> roses; int H, W, prob, NP, NR; TestCase(int _H, int _W, int _prob, int _NP) : NR(0) { H = (_H < 0 ? rnd.nextUInt(10, 80) : _H); W = (_W < 0 ? rnd.nextUInt(10, 80) : _W); prob = (_prob < 0 ? rnd.nextUInt(10, 90) : _prob); NP = (_NP < 0 ? rnd.nextUInt(5, 2 * (H + W)) : _NP); roses = vector<string>(H, string(W, '.')); for (int i = 0; i < H; i++) { for (int j = 0; j < W; j++) { int t = rnd.nextUInt(100); roses[i][j] = (t < prob ? 'R' : '.'); if (roses[i][j] == 'R') { NR++; } } } } }; namespace NGlobals { int H, W, NP, NR; vector<string> roses; vector<vector<int>> cumuRoses; } typedef pair<int, int> Point; typedef pair<Point, Point> Edge; typedef struct Piece { int idx; vector<Point> vs; double area; int rosesInside; int numDivide; Piece() {} Piece(int _idx, vector<Point>& _vs, double _area, int _rosesInside, int _numDivide) : idx(_idx), vs(_vs), area(_area), rosesInside(_rosesInside), numDivide(_numDivide) {} }*PiecePtr; ostream& operator<<(ostream& os, const Piece& p) { os << "idx: " << p.idx << endl; os << "vs: " << p.vs << endl; os << "area: " << p.area << endl; os << "rosesInside: " << p.rosesInside << endl; os << "numDivide: " << p.numDivide << endl; return os; } int gcd(int m, int n) { if (m == 0 && n == 0) return 0; if (m < 0) m = -m; if (n < 0) n = -n; if (m == 0 || n == 0) return max(m, n); int k; do { k = m % n; m = n; n = k; } while (k != 0); return m; } bool isInside(const Point& p) { return Y(p) >= 0 && Y(p) <= NGlobals::H && X(p) >= 0 && X(p) <= NGlobals::W; } bool isOnEdge(const Edge& e, const Point& p) { int x0 = e.first.first, y0 = e.first.second; int x1 = e.second.first, y1 = e.second.second; int x = p.first, y = p.second; if ((x - x0) * (y1 - y0) != (y - y0) * (x1 - x0)) return false; if (x0 == x1) return min(y0, y1) < y && y < max(y0, y1); if (y0 == y1) return min(x0, x1) < x && x < max(x0, x1); return abs(x - x0) < abs(x1 - x0) && abs(y - y0) < abs(y1 - y0); } // edge(vs[i], vs[i + 1]) に存在するとき i を返す int isOnConvexHull(const vector<Point>& vs, const Point& p) { for (int i = 0; i < vs.size() - 1; i++) if (isOnEdge(Edge(vs[i], vs[i + 1]), p)) return i; return -1; } // 凸包の隣接する頂点がなす辺と一致するか // 順方向なら 1, 逆方向なら -1 int isCutOnPiece(const vector<Point>& vs, const Edge& cut) { Point p1 = cut.first, p2 = cut.second; // p1 ある？ int n = vs.size() - 1; int idx = -1; for (int i = 0; i < n; i++) if (vs[i] == p1) idx = i; if (idx == -1) return 0; // 順方向？ if (vs[idx + 1] == p2) return 1; // 逆方向？ if (vs[(idx + n - 1) % n] == p2) return -1; // 存在しない return 0; } Point isCutOnPiece2(const vector<Point>& vs, const Edge& cut) { Point p1 = cut.first, p2 = cut.second; // p1 ある？ int n = vs.size() - 1; int idx = -1; for (int i = 0; i < n; i++) if (vs[i] == p1) idx = i; if (idx == -1) return Point(-1, -1); // 順方向？ if (vs[idx + 1] == p2) return Point(idx, 1); // 逆方向？ if (vs[(idx + n - 1) % n] == p2) return Point((idx + n - 1) % n, -1); // 存在しない return Point(-1, -1); } int countRosesInsideConvexHull(const vector<Point>& vs) { using namespace NGlobals; int n = 0; vector<int> pos(H, -1); for (int i = 0; i < vs.size() - 1; i++) { int x1 = X(vs[i]), y1 = Y(vs[i]); int x2 = X(vs[i + 1]), y2 = Y(vs[i + 1]); int dx = x2 - x1, dy = y2 - y1; if (!dy) continue; int ymin = (min(y1, y2) << 1), ymax = (max(y1, y2) << 1); for (int y = ymin + 1; y <= ymax - 1; y += 2) { int x = (x1 << 1) + dx * (y - (y1 << 1)) / dy; x += (x & 1) + 1; int yy = (y >> 1); if (pos[yy] < 0) pos[yy] = (x >> 1); else n += abs(cumuRoses[yy][x >> 1] - cumuRoses[yy][pos[yy]]); } } return n; } int countRosesInsideConvexHull2(const vector<Point>& vs) { //dump(vs); using namespace NGlobals; int n = 0; vector<int> pos(H, -1); int res = 0; for (int i = 0; i < vs.size() - 1; i++) { int x1 = X(vs[i]), y1 = Y(vs[i]); int x2 = X(vs[i + 1]), y2 = Y(vs[i + 1]); int dx = x2 - x1, dy = y2 - y1; if (!dy) continue; // dy が負なら左端 正なら右端 int ymin = (min(y1, y2) << 1), ymax = (max(y1, y2) << 1); int d = (dy < 0 ? 1 : 0); if (dy < 0) { dy = -dy; dx = -dx; } for (int y = ymin + 1; y <= ymax - 1; y += 2) { int n = dx * (y - (y1 << 1)) + dy * (x1 << 1); int x = n / dy; int r = n % dy; int yy = (y >> 1), xx; if (!r && (x & 1)) { xx = (x + d) >> 1; } else { x += (1 << (x & 1)); xx = x >> 1; } if (pos[yy] == -1) pos[yy] = xx; else res += abs(cumuRoses[yy][pos[yy]] - cumuRoses[yy][xx]); } } return res; } double getConvexHullArea(const vector<Point>& vs) { double a = 0; for (int i = 0; i < vs.size() - 1; i++) { Point p1 = vs[i], p2 = vs[i + 1]; a += (X(vs[i + 1]) - X(vs[i])) * (Y(vs[i]) + Y(vs[i + 1])); } return abs(a) / 2.0; } void addPointsOnEdge(vector<Point>& pts, int loc, const Point p1, const Point p2) { int dx = X(p2) - X(p1), dy = Y(p2) - Y(p1); int g = gcd(dx, dy); if (g == 1) return; dx /= g; dy /= g; vector<Point> pts_new(pts.size() + g - 1); for (int i = 0; i < loc; i++) pts_new[i] = pts[i]; for (int i = 1; i < g; i++) pts_new[loc + i - 1] = Point(X(p1) + dx * i, Y(p1) + dy * i); for (int i = loc; i < pts.size(); i++) pts_new[loc + g - 1 + i] = pts[i]; pts = pts_new; } vector<vector<Point>> getAllPointsOnConvexHull(const vector<Point>& vs) { vector<vector<Point>> pts(vs.size() - 1); for (int i = 0; i < vs.size() - 1; i++) { //pts[i].push_back(vs[i]); addPointsOnEdge(pts[i], pts[i].size(), vs[i], vs[i + 1]); //pts[i].push_back(vs[i + 1]); } return pts; } double eval_combine(double stdDevArea, double stdDevRoses) { return stdDevArea + stdDevRoses; } double eval_a2r1(double stdDevArea, double stdDevRoses) { return stdDevArea + stdDevArea + stdDevRoses; } double eval_a1r2(double stdDevArea, double stdDevRoses) { return stdDevArea + stdDevRoses + stdDevRoses; } double eval_area(double stdDevArea, double stdDevRoses) { return stdDevArea * 10000 + stdDevRoses; } double eval_roses(double stdDevArea, double stdDevRoses) { return stdDevArea + stdDevRoses * 10000; } double eval_default(double stdDevArea, double stdDevRoses) { return (1.0 + stdDevArea) * (1.0 + stdDevRoses); } class StateGreedy { public: vector<Piece> pieces; vector<Edge> cuts; stack<Piece> pieceStack; double calcScore() const { double stdDevArea, stdDevRoses; double meanArea, meanRoses; stdDevArea = stdDevRoses = meanArea = meanRoses = 0; int n = pieces.size(); for (int i = 0; i < n; i++) { const Piece& p = pieces[i]; meanArea += p.area; meanRoses += p.rosesInside; stdDevArea += p.area * p.area; stdDevRoses += p.rosesInside * p.rosesInside; } meanArea /= n; meanRoses /= n; stdDevArea /= n; stdDevRoses /= n; stdDevArea = sqrt(stdDevArea - meanArea * meanArea); stdDevRoses = sqrt(stdDevRoses - meanRoses * meanRoses); return (1.0 + stdDevArea) * (1.0 + stdDevRoses); } double calcPseudoScore(function<double(double, double)> eval_func) const { using namespace NGlobals; double stdDevArea, stdDevRoses; double meanArea, meanRoses; stdDevArea = stdDevRoses = meanArea = meanRoses = 0; int n = pieces.size(); for (int i = 0; i < n; i++) { meanArea += pieces[i].area; meanRoses += pieces[i].rosesInside; } meanArea /= NP; meanRoses /= NP; for (int i = 0; i < n; i++) { const Piece& p = pieces[i]; stdDevArea += p.numDivide * pow(p.area / p.numDivide - meanArea, 2); stdDevRoses += p.numDivide * pow(1.0 * p.rosesInside / p.numDivide - meanRoses, 2); } stdDevArea = sqrt(stdDevArea / NP); stdDevRoses = sqrt(stdDevRoses / NP); return eval_func(stdDevArea, stdDevRoses); } StateGreedy() {} void init() { using namespace NGlobals; Piece piece; vector<Point> vs({ { 0, 0 },{ W, 0 },{ W, H },{ 0, H },{ 0, 0 } }); piece.idx = 0; piece.vs = vs; piece.area = getConvexHullArea(vs); piece.rosesInside = NR; piece.numDivide = NP; pieces.push_back(piece); for (int i = 0; i < 4; i++) cuts.push_back(Edge(vs[i], vs[i + 1])); } // p1, p2 は piece[idx] の境界上に存在することが保証される // pi は Edge(v[loci], v[loci + 1]) 上に存在することが保証される // loc1 < loc2 bool makeCut(int idx, int loc1, Point p1, int loc2, Point p2) { // 十分小さい領域に分割されている if (pieces[idx].numDivide == 1) return false; //dump(idx, loc1, p1, loc2, p2); Piece piece_old = pieces[idx]; vector<Point> vs_total = piece_old.vs; // 逆方向 p2 -> p1 -> ... vector<Point> vs_new(loc2 - loc1 + 3); vs_new[0] = vs_new.back() = p2; vs_new[1] = p1; for (int i = loc1 + 1; i <= loc2; i++) vs_new[i - loc1 + 1] = vs_total[i]; // これいらない？ double area_new = getConvexHullArea(vs_new); if (area_new < 0.25) return false; double area_old = piece_old.area - area_new; if (area_old < 0.25) return false; // cut できる pieceStack.push(pieces[idx]); int div_total = piece_old.numDivide, div_new, div_old; double ratio_old = area_old / (area_old + area_new); double ratio_new = 1.0 - ratio_old; div_new = max(int(round(div_total * ratio_new)), 1); div_old = div_total - div_new; int roses_new = countRosesInsideConvexHull(vs_new); pieces.push_back(Piece(pieces.size(), vs_new, area_new, roses_new, div_new)); //piece_old = pieces[idx]; //vs_total = piece_old.vs; int n = vs_total.size(); vector<Point> vs_old(n + loc1 - loc2 + 2); vs_old[0] = vs_old.back() = p1; vs_old[1] = p2; for (int i = loc2 + 1; i < n; i++) vs_old[i - loc2 + 1] = vs_total[i]; for (int i = 1; i <= loc1; i++) vs_old[n - loc2 + i] = vs_total[i]; pieces[idx].vs = vs_old; pieces[idx].area = area_old; pieces[idx].numDivide = div_old; pieces[idx].rosesInside = countRosesInsideConvexHull(vs_old); cuts.emplace_back(p1, p2); } void undoCut() { int idx = pieceStack.top().idx; pieces[idx] = pieceStack.top(); pieceStack.pop(); pieces.pop_back(); cuts.pop_back(); } }; ostream& operator<<(ostream& os, const StateGreedy& s) { os << "pieces: " << s.pieces << endl; os << "cuts: " << s.cuts << endl; return os; } #ifdef VISUALIZE class CakeSharingVisGreedy { public: CakeSharingVisGreedy() {} void visualize(const StateGreedy& state, int delay = 0) { using namespace NGlobals; int pixsize = min(800 / H, 1200 / W); int boardH = H * pixsize; int boardW = W * pixsize; cv::Mat_<cv::Vec3b> img_board(boardH, boardW, cv::Vec3b(255, 255, 255)); const vector<Piece>& pieces = state.pieces; for (int i = 0; i < pieces.size(); i++) { for (int j = 0; j < pieces[i].vs.size() - 1; j++) { Point s = pieces[i].vs[j]; Point g = pieces[i].vs[j + 1]; cv::line(img_board, cv::Point(X(s) * pixsize, Y(s) * pixsize), cv::Point(X(g) * pixsize, Y(g) * pixsize), cv::Scalar(0, 0, 0), 2, cv::LINE_AA); } } int offset = 10; cv::Mat_<cv::Vec3b> img(img_board.rows + 20, img_board.cols + 20, cv::Vec3b(200, 200, 200)); cv::Rect roi(10, 10, img_board.cols, img_board.rows); cv::Mat img_roi = img(roi); img_board.copyTo(img_roi); cv::imshow("vis", img); cv::waitKey(delay); } }; #endif class SolveGreedy { public: vector<StateGreedy> states; SolveGreedy() {} StateGreedy solve() { using namespace NGlobals; vector<function<double(double, double)>> eval_func({ eval_combine, eval_area, eval_roses, eval_default, eval_a1r2, eval_a2r1 }); double min_method_score = DBL_MAX; StateGreedy min_method_state; for (int f = 0; f < eval_func.size(); f++) { StateGreedy state; state.init(); //CakeSharingVis(state).visualize(); for (int n = 1; n <= NP - 1; n++) { //dump(n, timer.time() - timer.t, time_elapses); double min_score = DBL_MAX; int min_i, min_e1, min_e2; Point min_p1, min_p2; for (int i = 0; i < state.pieces.size(); i++) { //dump(i); if (state.pieces[i].numDivide == 1) continue; vector<vector<Point>> pts = getAllPointsOnConvexHull(state.pieces[i].vs); //dump(pts); for (int e1 = 0; e1 < pts.size() - 1; e1++) { for (int e2 = e1 + 1; e2 < pts.size(); e2++) { for (int v1 = 0; v1 < pts[e1].size(); v1++) { for (int v2 = 0; v2 < pts[e2].size(); v2++) { //dump(e1, e2, pts[e1][v1], pts[e2][v2]); Point p1 = pts[e1][v1]; Point p2 = pts[e2][v2]; if (X(p1) != X(p2) && Y(p1) != Y(p2)) continue; state.makeCut(i, e1, p1, e2, p2); double score = state.calcPseudoScore(eval_func[f]); if (score < min_score) { min_score = score; min_i = i; min_e1 = e1; min_e2 = e2; min_p1 = p1; min_p2 = p2; } state.undoCut(); } } } } } state.makeCut(min_i, min_e1, min_p1, min_e2, min_p2); //CakeSharingVis(state).visualize(); } double score = state.calcScore(); states.push_back(state); //dump(f, score); if (score < min_method_score) { min_method_score = score; min_method_state = state; } } //CakeSharingVis(state).visualize(); return min_method_state; } }; class CutGraph { public: // 辺の依存グラフ // 0, 1, 2, 3 は最初の長方形の辺 // 隣接リストで管理 vector<set<int>> G; // 親のリスト vector<set<int>> PG; set<int> removable; CutGraph() {} // 依存グラフの構築 void createDependencyGraph(const vector<Edge>& cuts) { G.resize(cuts.size()); PG.resize(cuts.size()); // -1 を入れておけば削除可能とはならない for (int i = 0; i < 4; i++) G[i].insert(-1); // 各辺について for (int i = 0; i < cuts.size(); i++) { // 各頂点について for (int j = 0; j < cuts.size(); j++) { if (isOnEdge(cuts[i], cuts[j].first)) { G[i].insert(j); PG[j].insert(i); } if (isOnEdge(cuts[i], cuts[j].second)) { G[i].insert(j); PG[j].insert(i); } } } // 子を持たないノードは削除可能 for (int i = 0; i < G.size(); i++) if (G[i].empty()) removable.insert(i); } void remove(int v) { // 親を見ていく for (auto p : PG[v]) { // 親から edgeIdx を削除する G[p].erase(v); if (G[p].empty()) removable.insert(p); } PG[v].clear(); removable.erase(v); } }; ostream& operator<<(ostream& os, const CutGraph& g) { os << "dependencyGraph: " << g.G << endl; os << "removable: " << g.removable << endl; return os; } #ifdef VISUALIZE class CakeSharingVisSA { public: CakeSharingVisSA() {} void visualize(const vector<Piece>& pieces, const vector<bool>& usedPiece, int delay = 0) { using namespace NGlobals; int pixsize = min(900 / H, 1800 / W); int boardH = H * pixsize; int boardW = W * pixsize; cv::Mat_<cv::Vec3b> img_board(boardH, boardW, cv::Vec3b(255, 255, 255)); for (int i = 0; i < pieces.size(); i++) { if (!usedPiece[i]) continue; for (int j = 0; j < pieces[i].vs.size() - 1; j++) { Point s = pieces[i].vs[j]; Point g = pieces[i].vs[j + 1]; cv::line(img_board, cv::Point(X(s) * pixsize, Y(s) * pixsize), cv::Point(X(g) * pixsize, Y(g) * pixsize), cv::Scalar(0, 0, 0), 3, cv::LINE_AA); } } int offset = 10; cv::Mat_<cv::Vec3b> img(img_board.rows + 20, img_board.cols + 20, cv::Vec3b(200, 200, 200)); cv::Rect roi(10, 10, img_board.cols, img_board.rows); cv::Mat img_roi = img(roi); img_board.copyTo(img_roi); cv::imshow("vis", img); cv::waitKey(delay); } }; #endif class StateSA { stack<int> vacantPieceIdx; stack<int> vacantEdgeIdx; Edge getPieceInfoShareEdge(const Edge& e) { Edge info; for (int i = 0; i < pieces.size(); i++) { if (!usedPiece[i]) continue; Point pos = isCutOnPiece2(pieces[i].vs, e); if (pos.first == -1) continue; (pos.second == 1 ? info.first : info.second) = Point(i, pos.first); } return info; } public: vector<Piece> pieces; vector<Edge> cuts; CutGraph dependency; vector<bool> usedPiece; vector<bool> usedEdge; list<Edge> cutsOrder; StateSA() {} StateSA(const StateGreedy& gstate) : pieces(gstate.pieces), cuts(gstate.cuts) { dependency.createDependencyGraph(cuts); usedPiece.resize(pieces.size(), true); usedEdge.resize(cuts.size(), true); for (int i = 4; i < cuts.size(); i++) cutsOrder.push_back(cuts[i]); } double calcScore() const { double stdDevArea, stdDevRoses; double meanArea, meanRoses; stdDevArea = stdDevRoses = meanArea = meanRoses = 0; int n = pieces.size(); for (int i = 0; i < n; i++) { const Piece& p = pieces[i]; meanArea += p.area; meanRoses += p.rosesInside; stdDevArea += p.area * p.area; stdDevRoses += p.rosesInside * p.rosesInside; } meanArea /= n; meanRoses /= n; stdDevArea /= n; stdDevRoses /= n; stdDevArea = sqrt(stdDevArea - meanArea * meanArea); stdDevRoses = sqrt(stdDevRoses - meanRoses * meanRoses); return (1.0 + stdDevArea) * (1.0 + stdDevRoses); } double calcPseudoScore(function<double(double, double)> eval_func) const { using namespace NGlobals; double stdDevArea, stdDevRoses; double meanArea, meanRoses; stdDevArea = stdDevRoses = meanArea = meanRoses = 0; int n = pieces.size(); for (int i = 0; i < n; i++) { if (!usedPiece[i]) continue; //dump(pieces[i].area); meanArea += pieces[i].area; meanRoses += pieces[i].rosesInside; } meanArea /= NP; meanRoses /= NP; //dump(meanArea, meanRoses); for (int i = 0; i < n; i++) { if (!usedPiece[i]) continue; const Piece& p = pieces[i]; stdDevArea += p.numDivide * pow(p.area / p.numDivide - meanArea, 2); stdDevRoses += p.numDivide * pow(1.0 * p.rosesInside / p.numDivide - meanRoses, 2); } stdDevArea = sqrt(stdDevArea / NP); stdDevRoses = sqrt(stdDevRoses / NP); //dump(stdDevArea, stdDevRoses); return eval_func(stdDevArea, stdDevRoses); } // スコアの差分を計算する // 返り値は edge と分割された piece のペア double calcScoreDiff(int pieceIdx, int loc1, Point p1, int loc2, Point p2) { int pIdx1 = pieceIdx; const Piece& piece_total = pieces[pIdx1]; double area_total = piece_total.area; double roses_total = piece_total.rosesInside; vector<Point> vs_total = piece_total.vs; vector<Point> vs_new(loc2 - loc1 + 3); vs_new[0] = vs_new.back() = p2; vs_new[1] = p1; for (int i = loc1 + 1; i <= loc2; i++) vs_new[i - loc1 + 1] = vs_total[i]; int n = vs_total.size(); vector<Point> vs_old(n + loc1 - loc2 + 2); vs_old[0] = vs_old.back() = p1; vs_old[1] = p2; for (int i = loc2 + 1; i < n; i++) vs_old[i - loc2 + 1] = vs_total[i]; for (int i = 1; i <= loc1; i++) vs_old[n - loc2 + i] = vs_total[i]; double area_new = getConvexHullArea(vs_new); double area_old = area_total - area_new; double roses_new = countRosesInsideConvexHull2(vs_new); double roses_old = countRosesInsideConvexHull2(vs_old); int div_total = piece_total.numDivide; double ratio_old = area_old / area_total; double ratio_new = 1.0 - ratio_old; int div_new = max(int(round(div_total * ratio_new)), 1); int div_old = div_total - div_new; using namespace NGlobals; double meanArea = 1.0 * H * W / NP, meanRoses = 0; for (int i = 0; i < pieces.size(); i++) if (usedPiece[i]) meanRoses += pieces[i].rosesInside; meanRoses /= NP; double old_score, new_score; { double varArea = 0, varRoses = 0; for (int i = 0; i < pieces.size(); i++) { varArea += pieces[i].numDivide * pow(pieces[i].area / pieces[i].numDivide - meanArea, 2); varRoses += pieces[i].numDivide * pow(1.0 * pieces[i].rosesInside / pieces[i].numDivide - meanRoses, 2); } old_score = (1.0 + sqrt(varArea / NP)) * (1.0 + sqrt(varRoses / NP)); meanRoses = (meanRoses * NP - roses_total + roses_old + roses_new) / NP; int div_total = pieces[pIdx1].numDivide; varArea += div_new * pow(area_new / div_new - meanArea, 2) + div_old * pow(area_old / div_old - meanArea, 2) - div_total * pow(area_total / div_total - meanArea, 2); varRoses += div_new * pow(1.0 * roses_new / div_new - meanRoses, 2) + div_old * pow(1.0 * roses_old / div_old - meanRoses, 2) - div_total * pow(1.0 * roses_total / div_total - meanRoses, 2); new_score = (1.0 + sqrt(varArea / NP)) * (1.0 + sqrt(varRoses / NP)); } return new_score - old_score; } // 返り値は edge と分割された piece のペア pair<int, pair<int, int>> makeCut(int pieceIdx, int loc1, Point p1, int loc2, Point p2) { int pIdx1 = pieceIdx; int pIdx2 = vacantPieceIdx.top(); vacantPieceIdx.pop(); usedPiece[pIdx2] = true; int eIdx = vacantEdgeIdx.top(); vacantEdgeIdx.pop(); // 依存グラフに追加 (あとで関数化) for (int i = 0; i < cuts.size(); i++) { if (!usedEdge[i]) continue; if (isOnEdge(cuts[i], p1) || isOnEdge(cuts[i], p2)) { dependency.G[i].insert(eIdx); dependency.removable.erase(i); dependency.PG[eIdx].insert(i); } } usedEdge[eIdx] = true; cuts[eIdx] = Edge(p1, p2); cutsOrder.push_back(cuts[eIdx]); Piece piece_old = pieces[pIdx1]; vector<Point> vs_total = piece_old.vs; vector<Point> vs_new(loc2 - loc1 + 3); vs_new[0] = vs_new.back() = p2; vs_new[1] = p1; for (int i = loc1 + 1; i <= loc2; i++) vs_new[i - loc1 + 1] = vs_total[i]; double area_new = getConvexHullArea(vs_new); //if (area_new < 0.25) // return false; double area_old = piece_old.area - area_new; //if (area_old < 0.25) // return false; int div_total = piece_old.numDivide, div_new, div_old; double ratio_old = area_old / (area_old + area_new); double ratio_new = 1.0 - ratio_old; div_new = max(int(round(div_total * ratio_new)), 1); div_old = div_total - div_new; int roses_new = countRosesInsideConvexHull2(vs_new); pieces[pIdx2] = Piece(pIdx2, vs_new, area_new, roses_new, div_new); int n = vs_total.size(); vector<Point> vs_old(n + loc1 - loc2 + 2); vs_old[0] = vs_old.back() = p1; vs_old[1] = p2; for (int i = loc2 + 1; i < n; i++) vs_old[i - loc2 + 1] = vs_total[i]; for (int i = 1; i <= loc1; i++) vs_old[n - loc2 + i] = vs_total[i]; int roses_old = countRosesInsideConvexHull2(vs_old); pieces[pIdx1] = Piece(pIdx1, vs_old, area_old, roses_old, div_old); return pair<int, pair<int, int>>(eIdx, pair<int, int>(pIdx1, pIdx2)); } // ピース pieceIdx に対しスコアを最小化するような分割を行う bool makeBestCut(int pieceIdx) { if (pieces[pieceIdx].numDivide == 1) return true; vector<vector<Point>> pts = getAllPointsOnConvexHull(pieces[pieceIdx].vs); //dump(pts); int numPrimeEdge = 0; for (int e = 0; e < pts.size(); e++) if (pts[e].empty()) numPrimeEdge++; // valid な分割が存在しない if (numPrimeEdge >= pts.size() - 1) return false; double min_scorediff = DBL_MAX; int min_e1, min_e2; Point min_p1, min_p2; for (int e1 = 0; e1 < pts.size() - 1; e1++) { for (int e2 = e1 + 1; e2 < pts.size(); e2++) { for (int v1 = 0; v1 < pts[e1].size(); v1++) { for (int v2 = 0; v2 < pts[e2].size(); v2++) { Point p1 = pts[e1][v1]; Point p2 = pts[e2][v2]; // 斜めでもなんでもあり． double scorediff = calcScoreDiff(pieceIdx, e1, p1, e2, p2); if (scorediff < min_scorediff) { //dump(scorediff, p1, p2); min_scorediff = scorediff; min_e1 = e1; min_e2 = e2; min_p1 = p1; min_p2 = p2; } } } } } auto info = makeCut(pieceIdx, min_e1, min_p1, min_e2, min_p2); // 再帰的に分割していく if (!makeBestCut(info.second.first)) { removeEdge(info.first); return false; } if (!makeBestCut(info.second.second)) { removeEdge(info.first); return false; } return true; } bool makeRandomCut(int pieceIdx) { if (pieces[pieceIdx].numDivide == 1) return false; auto pts = getAllPointsOnConvexHull(pieces[pieceIdx].vs); int loc1 = rnd.nextUInt(pts.size()), loc2; while (true) { loc2 = rnd.nextUInt(pts.size()); if (loc1 != loc2) break; } if (loc1 > loc2) swap(loc1, loc2); int i1 = rnd.nextUInt(pts[loc1].size()); int i2 = rnd.nextUInt(pts[loc2].size()); Point p1 = pts[loc1][i1]; Point p2 = pts[loc2][i2]; makeCut(pieceIdx, loc1, p1, loc2, p2); return true; } // p1 は e を順方向に含む // 返り値はマージされたピースの id int mergePieces(int edgeIdx) { Edge e = cuts[edgeIdx]; Edge info = getPieceInfoShareEdge(e); int pIdx1 = info.first.first; int pIdx2 = info.second.first; Piece& p1 = pieces[pIdx1]; Piece& p2 = pieces[pIdx2]; int i1 = info.first.second; int i2 = info.second.second; int n1 = p1.vs.size() - 1; int n2 = p2.vs.size() - 1; //dump(info); //dump(p1, p2); p1.area += p2.area; p1.numDivide += p2.numDivide; vector<Point> vs_new; for (int i = i1 + 2; i < i1 + n1; i++) vs_new.push_back(p1.vs[i % n1]); for (int i = i2 + 2; i < i2 + n2; i++) vs_new.push_back(p2.vs[i % n2]); vs_new.push_back(vs_new.front()); //dump(p1.vs, p2.vs, vs_new); p1.vs = vs_new; p1.rosesInside = countRosesInsideConvexHull2(p1.vs); // pieces[pIdx2] は空き vacantPieceIdx.push(pIdx2); usedPiece[pIdx2] = false; // edge の更新 usedEdge[edgeIdx] = false; vacantEdgeIdx.push(edgeIdx); dependency.remove(edgeIdx); cutsOrder.erase((find(cutsOrder.begin(), cutsOrder.end(), e))); return pIdx1; } // edgeIdx および edgeIdx に依存するすべての辺の削除 // 返り値は edge を取り除いて得られた piece の id int removeEdge(int edgeIdx) { //dump(edgeIdx); // 子を持たないノードに達するまで降りる while (!dependency.G[edgeIdx].empty()) { removeEdge(*dependency.G[edgeIdx].begin()); } // 子を持たないノードに到達 // edgeIdx を共有する 2 つのピースをマージする return mergePieces(edgeIdx); } }; ostream& operator<<(ostream& os, const StateSA& s) { os << "pieces: " << s.pieces << endl; os << "cuts: " << s.cuts << endl; os << "dependency: " << s.dependency << endl; os << "validPiece: " << s.usedPiece << endl; return os; } class SolveSA { public: StateSA state; SolveSA() {} vector<Edge> solve() { SolveGreedy gsol; StateGreedy gstate = gsol.solve(); state = StateSA(gstate); // 強引に多点スタートに vector<StateSA> initStates; for (int i = 0; i < gsol.states.size(); i++) initStates.push_back(StateSA(gsol.states[i])); using namespace NGlobals; StateSA initState = state; double finalScore = initState.calcPseudoScore(eval_default); StateSA finalState = initState; while (timer.time() - timer.t < 9.3) { state = initStates[rnd.nextUInt(initStates.size())]; while (true) { int n = state.dependency.removable.size(); if (!n) break; vector<int> v; for (auto i : state.dependency.removable) v.push_back(i); int edgeIdx = v[rnd.nextUInt(n)]; int pieceIdx = state.removeEdge(edgeIdx); if (!state.makeBestCut(pieceIdx)) { cerr << "error!!!" << endl; } } double score = state.calcPseudoScore(eval_default); if (score < finalScore) { //cerr << score << endl; finalScore = score; finalState = state; } } state = finalState; while (state.cuts.size() != 4) state.cuts.pop_back(); for (auto e : state.cutsOrder) state.cuts.push_back(e); return state.cuts; } }; class CakeSharing { vector<string> getAnswer(const vector<Edge>& cuts) const { vector<string> res; for (int i = 4; i < cuts.size(); i++) { const auto& cut = cuts[i]; string s; s += to_string(X(cut.first)); s.push_back(' '); s += to_string(Y(cut.first)); s.push_back(' '); s += to_string(X(cut.second)); s.push_back(' '); s += to_string(Y(cut.second)); res.push_back(s); } return res; } void init(vector<string>& _roses, int _NP) { using namespace NGlobals; roses = _roses; H = roses.size(); W = roses[0].size(); NP = _NP; cumuRoses = vector<vector<int>>(H, vector<int>(W + 2, 0)); NR = 0; for (int i = 0; i < H; i++) { for (int j = 0; j < W; j++) { if (roses[i][j] == '.') continue; NR++; cumuRoses[i][j + 1]++; } } // cumulate for (int i = 0; i < H; i++) for (int j = 0; j < W + 1; j++) cumuRoses[i][j + 1] += cumuRoses[i][j]; } public: vector<string> cut(vector<string> _roses, int _NP) { timer.measure(); init(_roses, _NP); //cerr << state << endl; SolveSA ssol; //ssol.solve(); vector<string> ret = getAnswer(ssol.solve()); return ret; } }; // -------8<------- end of solution submitted to the website -------8<------- template<class T> void getVector(vector<T>& v) { for (int i = 0; i < v.size(); ++i) cin >> v[i]; } #ifdef OFFICIAL int main() { string filename = "testcase3.txt"; //ofstream ofs(filename); //ifstream ifs(filename); CakeSharing cs; int H = 0; //ifs >> H; cin >> H; vector<string> roses(H); int NP; //ifs >> roses; //ifs >> NP; cin >> roses; cin >> NP; //ofs << H << endl; //for (int i = 0; i < H; i++) // ofs << roses[i] << endl; //ofs << NP << endl; vector<string> ret = cs.cut(roses, NP); cout << ret.size() << endl; for (int i = 0; i < (int)ret.size(); ++i) cout << ret[i] << endl; cout.flush(); //cerr << timer.time() - timer.t << endl; //cerr << time_elapses << endl; } #endif
