けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

AOJ 2903 Board (AUPC 2018 day1-H)

2 変数劣モジュラ関数の和を最小カットで表すという、とても面白い問題!!!

問題へのリンク

問題概要

 H ×  W の二次元ボードが与えられる。以下のような "#" のところを最小個数の長方形で敷き詰めよ。例えば

4 10
##########
....#.....
....#.....
..........

は下のようになる。

4 10
1111111111
....2.....
....2.....
..........

制約

  •  1 \le R, C \le 25

問題の解釈

長方形を敷き詰めるのは、下図のように各 "#" を「縦」「横」に分類する作業だとみなすことができる。

----------
....|.....
....|.....
..........

そして上図で長方形が 2 個というのは以下のように解釈できる。

  • 向き付けする前時点で "#" マスは 12 マスある
  • "#" と "#" が隣接している箇所は 11 箇所あって、そのうちの 10 箇所が向きが整合している
  • 12 から「整合している箇所の個数」を引くと、2 になる

こう解釈すると、

  • "#" を「縦」か「横」に割り当てる
  • 隣接する 2 個が、「同じ属性」なら「利得」が入る

という感じになる。

2 変数劣モジュラ関数の和

0-1 変数が何個かあって、ある種の制約がたくさんくっついている問題を最小カットにするフレームワークがある。

  • p = 0 かつ q = 1 ならば、C のペナルティ
  • p = 1 かつ q = 0 ならば、C のペナルティ
  • p = 0 かつ q = 0 ならば、G の利得
  • p = 1 かつ q = 1 ならば、G の利得

といった制約は上手く扱うことができる。より一般には、 p, q の値に応じたコストを  f(p, q) と表したとき、

 f(0, 0) + f(1, 1) \le f(0, 1) + f(1, 0)

を満たしているとき、 f は 2 変数劣モジュラ関数であるという。2 変数劣モジュラ関数の和は、あるグラフのカットとして表せる。詳しくは次の記事を参照。

theory-and-me.hatenablog.com

より具体的に見て行こう。

p = 1 かつ q = 0 ならば、C のペナルティ

この場合の関数  f(p, q) の値は下図のように書ける。

 0 + 0 \le 0 + C

であるから、確かに 2 変数劣モジュラ関数である。

具体的には、下図のようなグラフのカットに対応する。

ここで、頂点  p, q について、

  • 頂点  s 側に属するとき: p, q の値が 1 である
  • 頂点  t 側に属するとき: p, q の値が 0 である

というように定めている。

p = 1 かつ q = 1 ならば、G の利得

この場合の関数  f(p, q) の値は下図のように書ける。

 0 - G \le 0 + 0

であるから、確かに 2 変数劣モジュラ関数である。

具体的にカットで表す方法を考える。 f(p, q) 全体に  G を足して、下図のような劣モジュラ関数をカットで表せればよい。

これは

  •  p = 0 のとき、 q の値によらず  G のコストがかかる
  •  p = 1 かつ  q = 0 のとき、 G のコストがかかる

というように解釈できる。よって、下図のようなカットに対応する。

コード

以上の考え方に基づいて、グラフの最小カット問題に帰着した上で最小カットを求めれば良い。

#include <bits/stdc++.h>
using namespace std;


// 2-variable submodular optimization
template<class T> struct TwoVariableSubmodularOpt {
    // edge class
    struct Edge {
        // core members
        int rev, from, to;
        T cap, icap, flow;
        
        // constructor
        Edge(int r, int f, int t, T c)
        : rev(r), from(f), to(t), cap(c), icap(c), flow(0) {}
        void reset() { cap = icap, flow = 0; }
        
        // debug
        friend ostream& operator << (ostream& s, const Edge& E) {
            return s << E.from << "->" << E.to << '(' << E.flow << '/' << E.icap << ')';
        }
    };
    
    // inner data
    int N, s, t;
    T offset;
    vector<vector<Edge>> list;
    vector<pair<int,int>> pos;
    
    // constructor
    TwoVariableSubmodularOpt() : N(0), s(0), t(0), offset(0) {}
    TwoVariableSubmodularOpt(int n) : N(n), s(n), t(n + 1), offset(0), list(n + 2) {}
    void init(int n = 0) {
        N = n, s = n, t = n + 1;
        offset = 0;
        list.assign(n + 2, Edge());
        pos.clear();
    }
    friend ostream& operator << (ostream& s, const TwoVariableSubmodularOpt &G) {
        const auto &edges = G.get_edges();
        for (const auto &e : edges) s << e << endl;
        return s;
    }

    // add 1-Variable submodular functioin
    void add_single_cost(int xi, T false_cost, T true_cost) {
        assert(0 <= xi && xi < N);
        if (false_cost >= true_cost) {
            offset += true_cost;
            add_edge(s, xi, false_cost - true_cost);
        } else {
            offset += false_cost;
            add_edge(xi, t, true_cost - false_cost);
        }
    }
    
    // add constraint (xi = T, xj = F is penalty C)
    void add_penalty(int xi, int xj, T cost) {
        assert(0 <= xi && xi < N);
        assert(0 <= xj && xj < N);
        assert(cost >= 0);
        add_edge(xi, xj, cost);
    }
    
    // add constraint (xi = T, xj = T is only cost 0, and the others are cost C)
    void add_both_true_benefit(int xi, int xj, T cost) {
        assert(0 <= xi && xi < N);
        assert(0 <= xj && xj < N);
        assert(cost >= 0);
        add_edge(xj, 0, cost);
        add_edge(xi, xj, cost);
    }
    
    // add constraint (xi = F, xj = F is only cost 0, and the others are cost C)
    void add_both_false_benefit(int xi, int xj, T cost) {
        assert(0 <= xi && xi < N);
        assert(0 <= xj && xj < N);
        assert(cost >= 0);
        add_edge(xi, cost, 0);
        add_edge(xi, xj, cost);
    }
    
    // add general 2-variable submodular function
    // (xi, xj) = (F, F): A, (F, T): B
    // (xi, xj) = (T, F): C, (T, T): D
    void add_submodular_function(int xi, int xj, T A, T B, T C, T D) {
        assert(0 <= xi && xi < N);
        assert(0 <= xj && xj < N);
        assert(B + C >= A + D);  // submodular constraint
        offset += A;
        add_single_cost(xi, 0, D - B);
        add_single_cost(xj, 0, B - A);
        add_penalty(xi, xj, B + C - A - D);
    }
    
    // solve
    T solve() {
        return dinic(s, t) + offset;
    }
    vector<bool> reconstruct() {
        vector<bool> res(N);
        queue<int> que;
        que.push(s);
        while (!que.empty()) {
            int v = que.front();
            que.pop();
            if (s < N) res[s] = true;
            for (const auto &e : list[v]) {
                if (e.cap && !res[e.to]) {
                    res[e.to] = true;
                    que.push(e.to);
                }
            }
        }
        return res;
    }
    
private:
    // add edge
    Edge &get_rev_edge(const Edge &e) {
        if (e.from != e.to) return list[e.to][e.rev];
        else return list[e.to][e.rev + 1];
    }
    Edge &get_edge(int i) {
        return list[pos[i].first][pos[i].second];
    }
    const Edge &get_edge(int i) const {
        return list[pos[i].first][pos[i].second];
    }
    vector<Edge> get_edges() const {
        vector<Edge> edges;
        for (int i = 0; i < (int)pos.size(); ++i) {
            edges.push_back(get_edge(i));
        }
        return edges;
    }
    void add_edge(int from, int to, T cap) {
        if (!cap) return;
        pos.emplace_back(from, (int)list[from].size());
        list[from].push_back(Edge((int)list[to].size(), from, to, cap));
        list[to].push_back(Edge((int)list[from].size() - 1, to, from, 0));
    }
    
    // Dinic's algorithm
    T dinic(int s, int t, T limit_flow) {
        T current_flow = 0;
        vector<int> level((int)list.size(), -1), iter((int)list.size(), 0);
        
        // Dinic BFS
        auto bfs = [&]() -> void {
            level.assign((int)list.size(), -1);
            level[s] = 0;
            queue<int> que;
            que.push(s);
            while (!que.empty()) {
                int v = que.front();
                que.pop();
                for (const Edge &e : list[v]) {
                    if (level[e.to] < 0 && e.cap > 0) {
                        level[e.to] = level[v] + 1;
                        if (e.to == t) return;
                        que.push(e.to);
                    }
                }
            }
        };
        
        // Dinic DFS
        auto dfs = [&](auto self, int v, T up_flow) {
            if (v == t) return up_flow;
            T res_flow = 0;
            for (int &i = iter[v]; i < (int)list[v].size(); ++i) {
                Edge &e = list[v][i], &re = get_rev_edge(e);
                if (level[v] >= level[e.to] || e.cap == 0) continue;
                T flow = self(self, e.to, min(up_flow - res_flow, e.cap));
                if (flow <= 0) continue;
                res_flow += flow;
                e.cap -= flow, e.flow += flow;
                re.cap += flow, re.flow -= flow;
                if (res_flow == up_flow) break;
            }
            return res_flow;
        };
        
        // flow
        while (current_flow < limit_flow) {
            bfs();
            if (level[t] < 0) break;
            iter.assign((int)iter.size(), 0);
            while (current_flow < limit_flow) {
                T flow = dfs(dfs, s, limit_flow - current_flow);
                if (!flow) break;
                current_flow += flow;
            }
        }
        return current_flow;
    };
    T dinic(int s, int t) {
        return dinic(s, t, numeric_limits<T>::max());
    }
};


void AOJ_2903() {
    int n, m;
    cin >> n >> m;
    vector<string> fi(n);
    for (int i = 0; i < n; ++i) cin >> fi[i];
    
    auto get_id = [&](int i, int j) -> int { return i * m + j; };
    
    // 0: 横, 1: 縦
    TwoVariableSubmodularOpt<int> tvs(n * m);
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
            if (fi[i][j] == '.') continue;
            tvs.add_single_cost(get_id(i, j), 1, 1);
            if (i+1 < n && fi[i+1][j] == '#') {
                // (1, 1) だけ 1 の利得 (-1 のコスト)
                tvs.add_submodular_function(get_id(i, j), get_id(i+1, j), 0, 0, 0, -1);
            }
            if (j+1 < m && fi[i][j+1] == '#') {
                // (0, 0) だけ 1 の利得 (-1 のコスト)
                tvs.add_submodular_function(get_id(i, j), get_id(i, j+1), -1, 0, 0, 0);
            }
        }
    }
    cout << tvs.solve() << endl;
}


int main() {
    AOJ_2903();
}

 

古い説明

かつては、「p = 1 かつ q = 1 ならば、G の利得」などを表すのには、スーパーノードが必要だと思い込んでいた。その方法も、それはそれで 2 変数の場合を一般化して、

「3 個以上の変数がすべて 1 ならば、G の利得」

などを表すのに使えるので、一応残しておく。

p = 1 かつ q = 1 ならば、G の利得

少し難しいが、スーパーノード用意して下図のようにする。

  • p = 1 かつ q = 1 のときのみ、0 の辺を選べる
  • p = 0 や q = 0 だと、C の辺しか選べない

という機構にしている。

仕組みとしては、例えば q = 0 を選んでおきながら 0 の辺を選ぼうとすると ∞ のコストを食らってしまう。