Marathon Match 87 "PopulationMapping"

2013年6月5日に初めてMarathonに参加して以来の2回目のMarathonです.
前回は,何をして良いのか分からずにただランダムに出力するということをしていました.今回は適当にしなかったので備忘録として書きます.
(まだ,結果が出てないので復習等は書いていないです)

結果

66人中34位でした.レートは 686 -> 886 でした.

問題

高さH,幅Wの格子状の領域,その領域内の総人口,人口の上限が与えられる.各々の格子は島か海のどちらか.矩形型のクエリを投げるとその領域内の島に住んでいる総人口が与えられる.
求めることは次のスコアを最大にする島を選択する.ただし,選んだ島の総人口は与えられた人口の上限以下でないといけない.

  スコア=(選んだ島の数)* 0.996^(クエリの数)

・H,W <= 500
・与えられた領域内には5〜15の街がありその周りに人口が集まる

観察

テストケースが10個与えられているので,クエリの数を無視してすべての格子内の人口が分かるとした前提で最適解を求めて観察してみます.
緑色で表示されている場所が最終的に選んだ島です.
観察すると人口の多い部分は選ばれていないので,人口の多い島を少ないクエリで求めるのをどうするかを考えることにしました.

f:id:tatanaideyo:20150828041828p:plain:w150 f:id:tatanaideyo:20150828041848p:plain:w150 f:id:tatanaideyo:20150828041905p:plain:w150
f:id:tatanaideyo:20150828042134p:plain:w150 f:id:tatanaideyo:20150828042149p:plain:w150 f:id:tatanaideyo:20150828042253p:plain:w150
f:id:tatanaideyo:20150828042217p:plain:w150 f:id:tatanaideyo:20150828042235p:plain:w150 f:id:tatanaideyo:20150828042203p:plain:w150
f:id:tatanaideyo:20150828042308p:plain:w150

方針

与えられた領域を矩形領域に分割して,最後にそれぞれの領域内の島の数を価値,人口を重さとしてナップサック問題を解き得られた解の領域内の島を選択するという方針で解きました.
領域の分割方法は次のように考えました.

与えられた領域に対して分割するかどうかの領域の集合をS1,
これ以上分割しない領域の集合をS2とする

1. S1に与えれた領域を入れる
2. S1の中で最も人口の多い領域rを選ぶ
3. rを分割するかどうかを決める

  3-1. 分割数が上限を超えている,rの人口が総人口のある割合以下,
     rの島の数が小さくないなら分割しないでS2にrを入れる
  3-2. それ以外はrを島の数で二等分して,
     片方の領域の人口を求めるてそれぞれをS1に入れる

4. S1の要素数が0ならば終了.それ以外なら2に戻る

3-2ではrの人口数が分かっているので,二等分した片方の人口が分かれば両方の人口が分かります.したがって,最終的なクエリの数は分割した領域の数に等しくなります.
ナップサック問題は総人口がint型としか分からず値が大きくなる可能性があるのと,分割数が高々40ぐらいだろうと思ったので半分全列挙で解きました(あり本を参考).

上の方針で考えれることとして,手順2の選ぶ領域の優先順位と,3の領域を分割するかどうかだと思います.選択する領域の優先順位を,

 ・島の数
 ・島の数/全体の島の数 + 人口/総人口
 ・島の数 * 人口

の大きい順に選ぶということをやったのですが一番良かったのが最も多い人口の領域を選ぶことでした.
3の領域を分割するかどうかは良くわからないまま終わりました.

テストケースの結果

青色の矩形はクエリをした領域です.
f:id:tatanaideyo:20150828055517p:plain:w150 f:id:tatanaideyo:20150828055525p:plain:w150 f:id:tatanaideyo:20150828055628p:plain:w150
f:id:tatanaideyo:20150828055632p:plain:w150 f:id:tatanaideyo:20150828055636p:plain:w150 f:id:tatanaideyo:20150828055653p:plain:w150
f:id:tatanaideyo:20150828055644p:plain:w150 f:id:tatanaideyo:20150828055648p:plain:w150 f:id:tatanaideyo:20150828055640p:plain:w150
f:id:tatanaideyo:20150828055657p:plain:w150

反省

twitterのタイムライン上の会話を見ると上位の人は人口の分布を推定するということをやっていました.この部分は考えようと思ったのですがよく分からずに放置していました.放置はいけないです.
観察の図を見ると選んでいない領域が菱型になっておりこれをどう扱うか考えていたのですがよく分からずに放置しました.
良さそうな方針を思いついても実装が重そうだとなかなか取り掛かれなかったり,実装してもあまり良くなかった場合に落ち込むので実装力は重要だと感じました.

残念なことにいくつかのテストケースで-1となっていました.調べた結果分かったことは領域内の人口を100倍にしたものをlong long 型で持っているのですが
 「100 * Population::queryRegion(0, 0, w - 1, h - 1)」
の計算等でint型でオーバーフローしてました.悔しすぎる orz

ソースコード

// ver-4 2015/8/19-
// Score: 921866.32, Rank: 13/42, Last Submission Time: 08.19.2015 01:14:08

// Seed-0, Area = 383, Queries = 21, Score = 352.082857960968
// Seed-1, Area = 749, Queries = 26, Score = 674.876983219446
// Seed-2, Area = 1128, Queries = 15, Score = 1062.182583362645
// Seed-3, Area = 421, Queries = 33, Score = 368.8419112942407
// Seed-4, Area = 9343, Queries = 20, Score = 8623.292496170909
// Seed-5, Area = 17718, Queries = 16, Score = 16617.43972337878
// Seed-6, Area = 1745, Queries = 21, Score = 1604.1373032425306
// Seed-7, Area = 42510, Queries = 18, Score = 39551.157368636836
// Seed-8, Area = 821, Queries = 20, Score = 757.7569452377519
// Seed-9, Area = 11212, Queries = 13, Score = 10642.765389098811

#include <bits/stdc++.h>

using namespace std;

typedef long long  ll;
typedef pair<int, int> P; // first: y座標, second: x座標


// 左上と右下の座標値を持つ矩形領域の人口を返す
// (関数内で座標値を問題文通りに変換)
class Population {
public:

    static int queryRegion(int x1, int y1, int x2, int y2);
};
int Population::queryRegion(int x1, int y1, int x2, int y2)
{
    printf("?\n%d %d %d %d\n", x1, y1, x2, y2);
    fflush(stdout);
    scanf("%d", &x1);

    return x1;
}

// 左上と右下の座標値を持つ矩形
struct Rect {
    P l, r;
    ll pop; // 領域内の人口
    int land; // 領域内のlandの数
    int h, w; // h: 高さ,w: 幅
    Rect() {};
    Rect(P l, P r, ll pop, int land) : l(l), r(r), pop(pop), land(land) {
        h = r.first - l.first + 1;
        w = r.second - l.second + 1;
    };
};
bool operator<(const Rect &a, const Rect &b) {
    return (a.land < b.land) || (a.land == b.land && a.pop < b.pop);
}

class PopulationMapping {
public:

    int max_percentage, total_population, h, w; // 50 <= h,w <= 500
    vector<string> world_map;  // '.': ocean, 'X': land
    int land_area_sum; // 島の総数
    int lim_partition; // 分割数の上限
    // 分割基準:population > lim_pop * ration_pop なら分割
    ll lim_pop;
    double ratio_pop;

    int min_land; // landの数がmin_land以下ならば分割をしない
    vector<vector<int>> land; // land[i][j] := Rect((0,0),(i,j))内のland数

    // 初期化
    void Init(int maxPercentage, vector<string> &worldMap,
              int totalPopulation);

    // 領域を分割する
    vector<Rect> DivideArea();

    // 左上lと右下rの矩形の領域にあるlandの数を返す
    int CountArea(const P &l, const P &r);

    // 価値を領域内のland 数,重みを領域内のpopulationとして最大化問題を
    // ナップサック問題(半分全列挙)で解く
    // 返す値は最適解に含まれない値(true)
    vector<bool> SolveKnapsackProblem(const vector<Rect> &rect);

    // world_mapと同様な形式をreturnする
    // map[i][j] := 'X': select, '.': unselect
    vector<string> mapPopulation(int maxPercentage, vector<string> worldMap,
                                 int totalPopulation);

private:
};

//////////////////////////////////////////////////////////////////////////

void PopulationMapping::Init(int maxPercentage, vector<string> &worldMap,
                             int totalPopulation)
{
    max_percentage = maxPercentage;
    total_population = totalPopulation;
    world_map = worldMap;

    h = world_map.size();
    w = world_map[0].size();
    lim_pop = max_percentage * total_population;

    land_area_sum = 0;
    for (int i = 0; i < h; ++i)
        for (int j = 0; j < w; ++j)
            if (world_map[i][j] == 'X')
                ++land_area_sum;


    // パラメータ設定
    lim_partition = 32;
    ratio_pop = 0.23;
    min_land = land_area_sum * min(((double)land_area_sum / (h * w)), 0.035);
    
    // land[i][j] := 左上(0,0)から右下(i,j)までの領域にあるlandの数の計算
    land.resize(h + 1);
    for (int i = 0; i <= h; ++i)
        land[i].resize(w + 1, 0);
    for (int i = 1; i <= h; ++i)
        for (int j = 1; j <= w; ++j) {
            if (world_map[i - 1][j - 1] == 'X')
                ++land[i][j];
            land[i][j] += land[i][j - 1];
            land[i][j] += land[i - 1][j];
            land[i][j] -= land[i - 1][j - 1];
        }
}

int PopulationMapping::CountArea(const P &l, const P &r)
{
    return land[r.first + 1][r.second + 1] - land[r.first + 1][l.second]
        - land[l.first][r.second + 1] + land[l.first][l.second];
}

vector<Rect> PopulationMapping::DivideArea()
{
    vector<Rect> rect;
    priority_queue<Rect> que;
    int num_area = 0;
    bool add;

    // Population::queryRegion をlong long型にキャストしないとオーバーフローになる
    que.push(Rect(P(0, 0), P(h - 1, w - 1),
                  100 * (ll)Population::queryRegion(0, 0, w - 1, h - 1),
                  CountArea(P(0, 0), P(h - 1, w - 1))));

    while (!que.empty()) {
        Rect now = que.top();
        que.pop();
        add = true;

        // 分割するかどうか
        if (num_area + que.size() < lim_partition && lim_pop * ratio_pop < now.pop
            && (now.h != 1 && now.w != 1) && (min_land < now.land)) {
            
            P tmp_r = now.l, tmp_l;
            const int lim_land = now.land / 2;

            if (now.h >= now.w) {
                // y軸方向に分割
                tmp_r.second = now.r.second;
                while (CountArea(now.l, tmp_r) < lim_land)
                    ++tmp_r.first;

                tmp_l.first = tmp_r.first + 1;
                tmp_l.second = now.l.second;
            }
            else {
                // x軸方向に分割
                tmp_r.first = now.r.first;
                while (CountArea(now.l, tmp_r) < lim_land)
                    ++tmp_r.second;
                
                tmp_l.first = now.l.first;
                tmp_l.second = tmp_r.second + 1;
            }

            // 分割
            if (tmp_r != now.r) {
                add = false;
                ll tmp_pop = 100 * Population::queryRegion(now.l.second, now.l.first,
                                                           tmp_r.second, tmp_r.first);

                que.push(Rect(now.l, tmp_r, tmp_pop, CountArea(now.l, tmp_r)));
                que.push(Rect(tmp_l, now.r, now.pop - tmp_pop, CountArea(tmp_l, now.r)));
            }
        }

        if (add) {
            rect.emplace_back(now);
            ++num_area;
        }
    }

    return rect;
}

vector<bool> PopulationMapping::SolveKnapsackProblem(const vector<Rect> &rect)
{
    const int n = rect.size();

    // 前半部分を全列挙
    tuple<ll, int, int> ps[1 << (n / 2)]; // (重さ,価値,組合せ)
    int n2 = n / 2;
    for (int s = 0; s < 1 << n2; ++s) {
        ll tmp_w = 0;
        int tmp_v = 0;
        for (int x = 0; x < n2; ++x) {
            if (s >> x & 1) {
                tmp_w += rect[x].pop;
                tmp_v += rect[x].land;
            }
        }
        ps[s] = make_tuple(tmp_w, tmp_v, s);
    }

    // 無駄な要素を取り除く
    sort(ps, ps + (1 << n2));
    int m = 1;
    for (int i = 1; i < 1 << n2; ++i) {
        if (get<1>(ps[m - 1]) < get<1>(ps[i])) {
            ps[m++] = ps[i];
        }
    }

    // 後ろ半分を全列挙して解を求める
    const ll lim = max_percentage * total_population;
    const int INF = h * w * 2;
    ll opt_v = 0;
    int sol1 = 0, sol2 = 0;
    for (int s = 0; s < 1 << (n - n2); ++s) {
        ll tmp_w = 0;
        int tmp_v = 0;
        for (int x = 0; x < n - n2; ++x) {
            if (s >> x & 1) {
                tmp_w += rect[n2 + x].pop;
                tmp_v += rect[n2 + x].land;
            }
        }

        if (tmp_w <= lim) {
            auto x = lower_bound(ps, ps + m, make_tuple(lim - tmp_w, INF, 0)) - 1;
            if (opt_v < get<1>(*x) + tmp_v) {
                opt_v = get<1>(*x) + tmp_v;
                sol1 = get<2>(*x);
                sol2 = s;
            }
        }
    }

    vector<bool> res(n, true);
    for (int x = 0; x < n2; ++x)
        if (sol1 >> x & 1)
            res[x] = false;
    for (int x = 0; x < n - n2; ++x)
        if (sol2 >> x & 1)
            res[n2 + x] = false;

    return res;
}

vector<string> PopulationMapping::mapPopulation(int maxPercentage,
                                                vector<string> worldMap,
                                                int totalPopulation)
{
    Init(maxPercentage, worldMap, totalPopulation);

    vector<Rect> rect = DivideArea();
    const int n = rect.size();

    vector<bool> not_sol = SolveKnapsackProblem(rect);

    for (int k = 0; k < n; ++k)
        if (not_sol[k]) {
            for (int i = rect[k].l.first; i <= rect[k].r.first; ++i)
                for (int j = rect[k].l.second; j <= rect[k].r.second; ++j)
                    world_map[i][j] = '.';
        }

    return world_map;
}
////////////////////////////////////////////////////////////////////////

int main()
{
    int h, max_percentage, total_population;

    // Input
    scanf("%d %d", &max_percentage, &h);
    vector<string> world_map(h);
    for (int i = 0; i < h; ++i)
        cin >> world_map[i];
    scanf("%d", &total_population);

    // Solve
    PopulationMapping pm;
    auto ret = pm.mapPopulation(max_percentage, world_map, total_population);

    // Output
    printf("%d\n", (int)ret.size());
    for (auto x : ret)
        cout << x << "\n";
    fflush(stdout);

    return 0;
}