「与えられたグラフを強連結成分分解すると DAG になるので、DAG 上で DP する」というのが想定解だが、フローでも解けると話題になった問題！

問題へのリンク

問題概要

頂点数 $N$ の有向グラフがある。最初、すべての頂点は白色である。以下の操作を 2 回行う。

頂点 $v$ を 1 つ選ぶ
頂点 $v$ を始点とするウォーク (同じ頂点は何度通ってもよい) を 1 つ選ぶ
そのウォークに含まれる頂点をすべて黒色に塗る。

2 回の操作後の黒色の頂点の個数の最大値を求めよ。

制約

$1 \le N \le 300$

考えたこと

強連結成分に含まれる 2 つの頂点 $u, v$ については、黒色に塗るかどうかについて同期していると考えてよい。頂点 $u$ を黒色に塗るならば、同じ強連結成分内の頂点 $v$ も黒色に塗って損はないためだ。

よって、まずグラフを強連結成分分解しよう。一つの強連結成分を一つの頂点とみなして DAG を新たに構築する。DAG の各頂点 $v$ の重み $w_{v}$ を、その強連結成分に含まれる頂点数とすると、次の問題になる。

DAG が与えられる。DAG の各頂点には重みがついている。この DAG 上で以下の操作を 2 回行う

頂点 $v$ を 1 つ選ぶ
頂点 $v$ を始点とするパス (同じ頂点を二度通ることはあり得ない) を 1 つ選ぶ
そのパスに含まれる頂点をすべて黒色に塗る。

2 回の操作後の黒色の頂点の重みの総和の最大値を求めよ。

この後は、DP でやるか、フローでやるかの方針がある。

解法 1：DP

基本的には次の配列 dp を更新すれば良さそうに思える。

dp[u][v] ← 頂点 $u$ を始点とするパスと、頂点 $v$ を始点とするパスに含まれる頂点の重みの総和の最大値

更新式は一見次のようにすればよさそうだが、注意したい点がある。

頂点 $u$ を始点とする各辺 $e = (u, w)$ に対して、

$u = v$ のとき：chmax(dp[u][u], dp[u][w])
$u \neq v$ のとき：chmax(dp[u][v], dp[w][v] + weight[u])

一つは、dp[w][v] を表す最適解において、頂点 $v$ を始点とするパスが頂点 $u$ を含んでいるかもしれないという点だ。この場合、dp[w][v] + weight[u] という値は、weight[u] をダブル加算した値となってしまう。

対策 1：トポロジカルソート順の活用

これを回避するためには、トポロジカルソート順を活用しよう。ここでは、DAG の頂点番号をトポロジカルソート順にソートされているものとする (ACL の scc を使うと自然とそうなる)。

このとき、 $u \lt v$ である場合に限り、各辺 $e = (u, w)$ に対して更新式 chmax(dp[u][v], dp[w][v] + weight[u]) を適用することにすれば問題ない。今度こそ、dp[w][v] を表す最適解において、頂点 $v$ を始点とするパスが頂点 $u$ を含むことはあり得ないからだ。

しかし、そのように DP の更新を制限すると、あらゆる解を網羅できるのかが不安になって来るだろう。たとえば、次の場合に困る。

頂点 $u$ から出ている辺はない (シンクである)
頂点 $v$ から出ている辺はある
しかし、トポロジカルソート順的には $u \lt v$ である

この場合、上記の DP 更新式における「トポロジカルソート順的に若い側からのみ頂点を進めて得られる頂点組を用いて更新する」という考え方を適用すると、うまく更新できない。

対策 2：スーパーシンクを用いる

しかし、頂点 $u, v$ から到達可能な頂点集合の中に共通の頂点 $t$ があれば事情は異なる。共通の到達可能頂点 $t$ があれば、 $u \lt v$ であっても

dp[t][t] についての解から遡って
dp[t][v] についての解を得て、そこから遡って
dp[u][v] についての解を得る

という更新が可能になるからだ。そこで、スーパーシンク $t$ を用意して、すべての頂点から $t$ への辺を張ることにしよう。これで解決だ。

計算量

以上の解法の計算量は $O(N^{3})$ となる。

解法 2：フロー

2 本のパスを考える問題は、多くの場合、容量 2 のフローを考える問題へと帰着可能である。

今回は各頂点に重みがある。最小費用流問題へと帰着させるため、基本的には、

各頂点 $v$ の重みを -1 倍する (その値を $c_{v}$ とする)
さらに、各頂点 $v$ を $v_{in}$ と $v_{out}$ に分けて、その間にコスト $c_{v}$ の辺を張る (常套手段)

というような方針を考えたい。しかし、問題がある。それは、2 本のウィークが同じ頂点を通るとき、コスト $c_{v}$ をダブルカウントしてしまうことだ。そこで、頂点 $v_{in}$ と $v_{out}$ の間に 2 本の辺を張ることにして、一方は容量を 1、コストを $c_{v}$ として、他方は容量を 1、コストを 0 とする。

この方針をもとに、まとめると、次のようなフローネットワークを作れば良い。スーパーソース $s$ と、スーパーシンク $t$ も用意しておく。

$s$ から $v_{in}$ へと、容量 2、コスト 0 の辺を張る
$v_{out}$ から $t$ へと、容量 2、コスト 0 の辺を張る
各辺 $e = (u, v)$ に対して、 $u_{out}$ から $v_{in}$ へと、容量 2、コスト 0 の辺を張る
$v_{in}$ から $v_{out}$ へと、容量 1、コスト $c_{v}$ の辺を張る
$v_{in}$ から $v_{out}$ へと、容量 1、コスト 0 の辺を張る

こうして得られたフローネットワーク上で、容量が 2 のフローを流して、最小コストを求めよう。計算量は $O(N^{3})$ となる。

コード (解法 1 - DP)

強連結成分分解については、ACL を活用した。

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

int main() {
    // 入力
    int N;
    cin >> N;
    vector<vector<int>> input(N, vector<int>(N));
    for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) cin >> input[i][j];
    
    // SCC を適用するためのグラフを構築する
    atcoder::scc_graph SG(N);
    for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) {
        if (input[i][j]) SG.add_edge(i, j);
    }
    
    // SCC
    vector<vector<int>> scc = SG.scc();
    
    // SCC 後の DAG を構築 (超頂点を含む)
    int t = (int)scc.size();
    vector<int> weight((int)scc.size() + 1, 1);  // 各頂点の重み
    vector<vector<int>> dag((int)scc.size() + 1);  // DAG の構造
    for (int i = 0; i < (int)scc.size(); ++i) {
        weight[i] = (int)scc[i].size();
        dag[i].push_back(t);
        for (int j = i+1; j < (int)scc.size(); ++j) {
            bool ok = false;
            for (auto a : scc[i]) for (auto b : scc[j]) if (input[a][b]) ok = true;
            if (ok) dag[i].push_back(j);
        }
    }

    // DP (dp[u][v] := u <= v のとき、u, v の 2 頂点を始点としたときの黒色頂点重みの総和の最大値)
    int V = (int)dag.size();
    vector<vector<int>> dp(V, vector<int>(V, -1));
    auto rec = [&](auto self, int u, int v) -> int {
        if (u > v) swap(u, v);
        if (dp[u][v] != -1) return dp[u][v];
        
        int res = weight[u] + (u != v ? weight[v] : 0);
        if (u == v) {
            for (auto u2 : dag[u]) {
                res = max(res, self(self, u, u2));
            }
        } else {
            for (auto u2 : dag[u]) {
                res = max(res, self(self, min(u2, v), max(u2, v)) + weight[u]);
            }
        }
        return dp[u][v] = res;
    };
    
    int res = 0;
    for (int i = 0; i < V; ++i) for (int j = i; j < V; ++j) {
        res = max(res, rec(rec, i, j));
    }
    cout << res - 1 << endl;
}

コード (解法 2 - フロー)

ACL の最小費用流は負コストの辺に対応していなさそうだった (知らないだけかもしれない) ので、自前の最小費用流ライブラリを活用した。

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

// edge class (for network-flow)
template<class FLOWTYPE, class COSTTYPE> struct FlowCostEdge {
    // core members
    int rev, from, to;
    FLOWTYPE cap, icap, flow;
    COSTTYPE cost;
    
    // constructor
    FlowCostEdge(int rev, int from, int to, FLOWTYPE cap, COSTTYPE cost)
    : rev(rev), from(from), to(to), cap(cap), icap(cap), flow(0), cost(cost) {}
    void reset() { cap = icap, flow = 0; }
    
    // debug
    friend ostream& operator << (ostream& s, const FlowCostEdge& e) {
        return s << e.from << " -> " << e.to
        << " (" << e.flow << "/" << e.icap << ", " << e.cost << ")";
    }
};

// graph class (for network-flow)
template<class FLOWTYPE, class COSTTYPE> struct FlowCostGraph {
    // core members
    vector<vector<FlowCostEdge<FLOWTYPE, COSTTYPE>>> list;
    vector<pair<int,int>> pos;  // pos[i] := {vertex, order of list[vertex]} of i-th edge
    
    // constructor
    FlowCostGraph(int n = 0) : list(n) { }
    void init(int n = 0) {
        list.assign(n, FlowCostEdge<FLOWTYPE, COSTTYPE>());
        pos.clear();
    }
    
    // getter
    vector<FlowCostEdge<FLOWTYPE, COSTTYPE>> &operator [] (int i) {
        return list[i];
    }
    const vector<FlowCostEdge<FLOWTYPE, COSTTYPE>> &operator [] (int i) const {
        return list[i];
    }
    size_t size() const {
        return list.size();
    }
    FlowCostEdge<FLOWTYPE, COSTTYPE> &get_rev_edge
    (const FlowCostEdge<FLOWTYPE, COSTTYPE> &e) {
        if (e.from != e.to) return list[e.to][e.rev];
        else return list[e.to][e.rev + 1];
    }
    FlowCostEdge<FLOWTYPE, COSTTYPE> &get_edge(int i) {
        return list[pos[i].first][pos[i].second];
    }
    const FlowCostEdge<FLOWTYPE, COSTTYPE> &get_edge(int i) const {
        return list[pos[i].first][pos[i].second];
    }
    vector<FlowCostEdge<FLOWTYPE, COSTTYPE>> get_edges() const {
        vector<FlowCostEdge<FLOWTYPE, COSTTYPE>> edges;
        for (int i = 0; i < (int)pos.size(); ++i) {
            edges.push_back(get_edge(i));
        }
        return edges;
    }
    
    // change edges
    void reset() {
        for (int i = 0; i < (int)list.size(); ++i) {
            for (FlowCostEdge<FLOWTYPE, COSTTYPE> &e : list[i]) e.reset();
        }
    }
    
    // add_edge
    void add_edge(int from, int to, FLOWTYPE cap, COSTTYPE cost) {
        pos.emplace_back(from, (int)list[from].size());
        list[from].push_back(FlowCostEdge<FLOWTYPE, COSTTYPE>
                             ((int)list[to].size(), from, to, cap, cost));
        list[to].push_back(FlowCostEdge<FLOWTYPE, COSTTYPE>
                           ((int)list[from].size() - 1, to, from, 0, -cost));
    }

    // debug
    friend ostream& operator << (ostream& s, const FlowCostGraph &G) {
        const auto &edges = G.get_edges();
        for (const auto &e : edges) s << e << endl;
        return s;
    }
};

// min-cost max-flow (<= limit_flow), slope ver.
template<class FLOWTYPE, class COSTTYPE> vector<pair<FLOWTYPE, COSTTYPE>>
MinCostFlowSlope(FlowCostGraph<FLOWTYPE, COSTTYPE> &G, int S, int T, FLOWTYPE limit_flow)
{
    // result values
    FLOWTYPE cur_flow = 0;
    COSTTYPE cur_cost = 0, pre_cost = -1;
    vector<pair<FLOWTYPE, COSTTYPE>> res;
    res.emplace_back(cur_flow, cur_cost);
    
    // intermediate values
    vector<bool> seen((int)G.size(), false);
    vector<COSTTYPE> dist((int)G.size(), numeric_limits<COSTTYPE>::max());
    vector<int> prevv((int)G.size(), -1), preve((int)G.size(), -1);
    
    // dual
    auto dual_step = [&]() -> bool {
        seen.assign((int)G.size(), false);
        dist.assign((int)G.size(), numeric_limits<COSTTYPE>::max());
        seen[S] = true, dist[S] = 0;
        while (true) {
            bool update = false;
            for (int v = 0; v < (int)G.size(); ++v) {
                if (!seen[v]) continue;
                for (int i = 0; i < G[v].size(); ++i) {
                    const FlowCostEdge<FLOWTYPE, COSTTYPE> &e = G[v][i];
                    if (e.cap > 0 && (!seen[e.to] || dist[e.to] > dist[v] + e.cost)) {
                        dist[e.to] = dist[v] + e.cost;
                        prevv[e.to] = v;
                        preve[e.to] = i;
                        seen[e.to] = true;
                        update = true;
                    }
                }
            }
            if (!update) break;
        }
        return seen[T];
    };
    
    // primal
    auto primal_step = [&]() -> void {
        FLOWTYPE flow = limit_flow - cur_flow;
        COSTTYPE cost = dist[T];
        for (int v = T; v != S; v = prevv[v]) {
            flow = min(flow, G[prevv[v]][preve[v]].cap);
        }
        for (int v = T; v != S; v = prevv[v]) {
            FlowCostEdge<FLOWTYPE, COSTTYPE> &e = G[prevv[v]][preve[v]];
            FlowCostEdge<FLOWTYPE, COSTTYPE> &re = G.get_rev_edge(e);
            e.cap -= flow, e.flow += flow;
            re.cap += flow, re.flow -= flow;
        }
        cur_flow += flow;
        cur_cost += flow * cost;
        if (pre_cost == cost) res.pop_back();
        res.emplace_back(cur_flow, cur_cost);
        pre_cost = cur_cost;
    };
    
    // primal-dual
    while (cur_flow < limit_flow) {
        if (!dual_step()) break;
        primal_step();
    }
    return res;
}

// min-cost max-flow, slope ver.
template<class FLOWTYPE, class COSTTYPE> vector<pair<FLOWTYPE, COSTTYPE>>
MinCostFlowSlope(FlowCostGraph<FLOWTYPE, COSTTYPE> &G, int S, int T)
{
    return MinCostFlowSlope(G, S, T, numeric_limits<FLOWTYPE>::max());
}

// min-cost max-flow (<= limit_flow)
template<class FLOWTYPE, class COSTTYPE> pair<FLOWTYPE, COSTTYPE>
MinCostFlow(FlowCostGraph<FLOWTYPE, COSTTYPE> &G, int S, int T, FLOWTYPE limit_flow)
{
    return MinCostFlowSlope(G, S, T, limit_flow).back();
}

// min-cost max-flow (<= limit_flow)
template<class FLOWTYPE, class COSTTYPE> pair<FLOWTYPE, COSTTYPE>
MinCostFlow(FlowCostGraph<FLOWTYPE, COSTTYPE> &G, int S, int T)
{
    return MinCostFlow(G, S, T, numeric_limits<FLOWTYPE>::max());
}

int main() {
    // 入力
    int N;
    cin >> N;
    vector<vector<int>> input(N, vector<int>(N));
    for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) cin >> input[i][j];
    
    // SCC を適用するためのグラフを構築する
    atcoder::scc_graph SG(N);
    for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) {
        if (input[i][j]) SG.add_edge(i, j);
    }
    
    // SCC
    vector<vector<int>> scc = SG.scc();
    
    // SCC 後の DAG を構築
    int V = (int)scc.size();
    vector<int> weight(V);  // 各頂点の重み
    vector<vector<int>> dag(V);  // DAG の構造
    for (int i = 0; i < V; ++i) {
        weight[i] = (int)scc[i].size();
        for (int j = i+1; j < (int)scc.size(); ++j) {
            bool ok = false;
            for (auto a : scc[i]) for (auto b : scc[j]) if (input[a][b]) ok = true;
            if (ok) dag[i].push_back(j);
        }
    }

    // フローネットワークの構築
    int s = V*2, t = V*2+1;
    FlowCostGraph<int, int> G(V*2 + 2);
    for (int v = 0; v < V; ++v) {
        G.add_edge(s, v, 2, 0);
        G.add_edge(v+V, t, 2, 0);
        for (auto w : dag[v]) G.add_edge(v+V, w, 2, 0);
        G.add_edge(v, v+V, 1, -weight[v]);
        G.add_edge(v, v+V, 1, 0);
    }
    
    auto [flow, cost] = MinCostFlow(G, s, t, 2);
    cout << -cost << endl;
}