燃やす埋めるは今度こそちゃんとマスターする！！！

問題概要

$H$ × $W$ の各格子点に電球がある。
$i + j$ が偶数となるような $(i, j)$ について、 $(i, j), (i, j+1), (i+1, j), (i+1, j+1)$ を頂点とする正方形の中心に装置があって、各装置を押すことで四隅の電球を点灯することができる。

ただし、装置を起動するにはそれぞれ $W$ のコストがかかる。また電球 $(i, j)$ が光ることによって得られる利得は $B_{ij}$ で与えられる。

いくつかの電源をオンにして「利得の合計 - 電気代の合計」を最大化しなさい。

制約

$2 \le H, W \le 50$

考えたこと

見るからに「燃やす埋める」だけど詰めるのにすごく時間がかかった。電球と装置との間の「報酬」「利得」を考えている限りはダメで、一見

電球がオンなのに装置がオフなのはダメ (ペナルティ INF)
装置がオンなのに電球がオフなのはダメ (ペナルティ INF)

という形でできそうなのだが、電球に対して一般に装置が 2 つあるので、例えば装置 p, q が電球 b を光らせることができるとき

電球 b がオンならば、装置 p または装置 q がオンである
装置 p または装置 q がオンならば、電球 b がオンである

といった関係性を表現しないといけなくて、とても手に負えそうにない。

電球を考えなくても良いように

そこで、装置だけで話が閉じるようにすることを考える。基本的には

装置 p を使うことで W のコストを払い、四隅の電球の分の利得 sumB を得る

という風に考えることで、装置だけの世界で考察することができる。ただし、斜め方向に隣り合う装置については

装置 p と装置 q とが電球 b を共有しているとき、両方を使う場合にはペナルティ B (電球 b の分) が発生する

という風に考えれば良い。こうすることで「燃やす埋める」の世界観に持ち込めそうに思える。

二部グラフを生かして片側を反転

しかしこのままだと「変数 p と q が同じ側にいるときにペナルティが発生する」という形の制約式になってしまう。これは燃やす埋めるで扱うことができない。燃やす埋めるで扱えるのは、

p = 0 かつ q = 1 のときにペナルティ
p = 1 かつ q = 0 のときにペナルティ
p = 0 かつ q = 0 のときに利得
p = 1 かつ q = 1 のときに利得

といった制約になる。しかし、今回、装置を頂点とした

装置 p と装置 q をともに使うとペナルティが発生するような (p, q) の間に辺をはる

としたグラフもまた、グリッドグラフであり、すなわち二部グラフであることに注意する。

そうすると具体的には

i % 2 == 0 の位置にある装置 p については、オンのとき 1、オフのとき 0
i % 2 == 1 の位置にある装置 q については、オンのとき 0、オフのとき 1

という風に i の偶奇に応じてオンオフ時の 0-1 を入れ替えてあげることで、

p = 1 かつ q = 0 のときは b のペナルティ

という式として扱うことができる。

詰め

こうして燃やす埋めるの世界観に完全に持ち込むことができた。あとは整理するだけである。

例えばサンプルの

4 4 10
100 100 100 100
100 100 100 100
1 100 100 1
1 1 1 1

の例では以下のようなグラフを作ることができる。

これに対して最小カットは以下のようになる。

図に見るように、装置 1, 2, 3 についてはオンにして、装置 4, 5 についてはオフにするのが最適だとわかる。

ここで、例えば装置 1 と装置 3 はパリティが異なるので、装置 1 がオン (変数値は 1) で装置 3 がオン (変数値は 0) だと、ペナルティ 100 が生じていることがわかる。装置 2 と装置 3 についても同様である。

最後に、求めた最小カットを、「各装置についての四隅の電球の和の総和」から引けば答えになる。

#include <iostream>
#include <vector>
#include <map>
#include <queue>
using namespace std;

// edge class (for network-flow)
template<class FLOWTYPE> struct Edge {
    int rev, from, to;
    FLOWTYPE cap, icap;
    Edge(int r, int f, int t, FLOWTYPE c) : rev(r), from(f), to(t), cap(c), icap(c) {}
};

// graph class (for network-flow)
template<class FLOWTYPE> struct Graph {
    vector<vector<Edge<FLOWTYPE> > > list;
    
    Graph(int n = 0) : list(n) { }
    void init(int n = 0) { list.clear(); list.resize(n); }
    void reset() { for (int i = 0; i < (int)list.size(); ++i) for (int j = 0; j < list[i].size(); ++j) list[i][j].cap = list[i][j].icap; }
    inline vector<Edge<FLOWTYPE> >& operator [] (int i) { return list[i]; }
    inline const size_t size() const { return list.size(); }
    
    inline Edge<FLOWTYPE> &redge(Edge<FLOWTYPE> e) {
        if (e.from != e.to) return list[e.to][e.rev];
        else return list[e.to][e.rev + 1];
    }
    
    void addedge(int from, int to, FLOWTYPE cap) {
        list[from].push_back(Edge<FLOWTYPE>((int)list[to].size(), from, to, cap));
        list[to].push_back(Edge<FLOWTYPE>((int)list[from].size() - 1, to, from, 0));
    }
    
    void add_undirected_edge(int from, int to, FLOWTYPE cap) {
        list[from].push_back(Edge<FLOWTYPE>((int)list[to].size(), from, to, cap));
        list[to].push_back(Edge<FLOWTYPE>((int)list[from].size() - 1, to, from, cap));
    }
};

template<class FLOWTYPE> struct Dinic {
    const FLOWTYPE INF = 1LL<<59; // to be set
    vector<int> level, iter;

    Dinic() { }
    void dibfs(Graph<FLOWTYPE> &G, int s) {
        level.assign((int)G.size(), -1);
        level[s] = 0;
        queue<int> que;
        que.push(s);
        while (!que.empty()) {
            int v = que.front();
            que.pop();
            for (int i = 0; i < G[v].size(); ++i) {
                Edge<FLOWTYPE> &e = G[v][i];
                if (level[e.to] < 0 && e.cap > 0) {
                    level[e.to] = level[v] + 1;
                    que.push(e.to);
                }
            }
        }
    }
    
    FLOWTYPE didfs(Graph<FLOWTYPE> &G, int v, int t, FLOWTYPE f) {
        if (v == t) return f;
        for (int &i = iter[v]; i < G[v].size(); ++i) {
            Edge<FLOWTYPE> &e = G[v][i], &re = G.redge(e);
            if (level[v] < level[e.to] && e.cap > 0) {
                FLOWTYPE d = didfs(G, e.to, t, min(f, e.cap));
                if (d > 0) {
                    e.cap -= d;
                    re.cap += d;
                    return d;
                }
            }
        }
        return 0;
    }
    
    FLOWTYPE solve(Graph<FLOWTYPE> &G, int s, int t) {
        level.assign((int)G.size(), -1); iter.assign((int)G.size(), 0);
        FLOWTYPE res = 0;
        while (true) {
            dibfs(G, s);
            if (level[t] < 0) return res;
            for (int i = 0; i < (int)iter.size(); ++i) iter[i] = 0;
            FLOWTYPE flow = 0;
            while ((flow = didfs(G, s, t, INF)) > 0) {
                res += flow;
            }
        }
    }
};

using pint = pair<int,int>;

long long N, M, W;
long long B[110][110];

map<pint,int> ma; // (i, j) -> 装置 ID

// (i, j) にある装置の四隅の電球の利得 
long long calc(int i, int j) {
    long long res = 0;
    for (int di = 0; di < 2; ++di) {
        for (int dj = 0; dj < 2; ++dj) {
            res += B[i+di][j+dj];
        }
    }
    return res;
}

int main() {
    cin >> N >> M >> W;
    ma.clear();
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < M; ++j) {
            cin >> B[i][j];
            if ( i+1 < N && j+1 < M && (i + j) % 2 == 0 ) 
                ma[pint(i, j)] = 0;            
        }
    }
    int iter = 0;
    for (auto &it : ma) it.second = iter++;
  
    // グラフを作る
    Graph<long long> G(ma.size() + 2);
    int s = ma.size(), t = s + 1;      
    for (int i = 0; i < N-1; ++i) {
        for (int j = 0; j < M-1; ++j) {
            if (!ma.count(pint(i, j))) continue;
            int id = ma[pint(i, j)];
            
            if (i % 2 == 0) { // 正規の装置
                G.addedge(s, id, calc(i, j));
                G.addedge(id, t, W);
            }
            else { // 反転した装置 
                G.addedge(s, id, W);
                G.addedge(id, t, calc(i, j));
            }
            
            // 正規から反転へとペナルティ辺をはる
            if (i % 2 == 0) {
                for (int di = -1; di <= 1; di += 2) {
                    for (int dj = -1; dj <= 1; dj += 2) {
                        int ni = i + di, nj = j + dj;
                        if (!ma.count(pint(ni, nj))) continue;
                        int id2 = ma[pint(ni, nj)];
                        int bi = i + (di == -1 ? 0 : 1);
                        int bj = j + (dj == -1 ? 0 : 1);
                        G.addedge(id, id2, B[bi][bj]);
                    }
                }
            }
        }
    }

    Dinic<long long> din;
    long long res = din.solve(G, s, t);

    long long all = 0;
    for (auto it : ma) all += calc(it.first.first, it.first.second);
    cout << all - res << endl;
}