CODE THANKS FESTIVAL 2018 G - Sum of Cards

普通に難しくて何時間もかかってしまった。
CODE THANKS FESTIVAL 2018 G - Sum of Cards

問題

Nを正の整数とする。表に1からNまでの整数が、裏に1からNまでの整数が書かれたカードがN枚テーブルの上に存在する。現在の配置について、表に書かれている数字は1からNまでそれぞれちょうど1個存在し、裏に書かれている数字についても同じことが言える。
N枚のカードのうち任意の枚数を裏返し、見える整数の合計を最大化したい。ただし、見える整数がK種類以上存在する必要がある。条件を満たす配置における整数の合計の最大値を求めよ。

制約

1 <= K <= N <= 5000
i枚目のカードの表にはA[i]が、裏にはB[i]が書かれている。
- A[1], ..., A[N]は1からNまでの並び替えで、B[1], ..., B[N]は1からNまでの並び替え

解法

2回以上登場する数字の種類の数をuとすると、見える整数がK種類以上であることとu <= N - Kは同値。
ここで、uをN-K以下にする代わりに、実数xを決めたとき各整数を2回以上選択した場合のペナルティをxにした場合の最大化問題を考える(ラグランジュ緩和)。これは最小カット問題に帰着できる:

頂点数2 + N
1 -> 2 + i: 重み -B[i]
2 + i -> 2: 重み -A[i]
A[u] = B[v]のとき 2 + u -> 2 + v: 重みx

このネットワークにおいて1と2をカットするとき、1側の頂点がA[i]を選ぶカード、2側の頂点がB[i]を選ぶカードに対応する。
このネットワークの最小カットはDinic法などで計算できる。これをxを変えて実行し二分探索して、重みxの辺をカットする回数がN-K以下になる最大のxを見つける。
xはどのあたりの範囲で探索すればよいだろうか? まず範囲について考えよう。x=Nなら重みxの辺は全くカットされず、x=0なら全く制限がないのと同じなので、0からNまでの範囲を探索すればOKである。探索精度はどうか? これは難しいが、実験的には1/(N+1)くらいの精度になるまで探索すれば十分。
計算量はO(N^3 log N)だが、Dinic法は一般に最悪計算量よりも高速に動作するので、定数倍に気をつければこれで通る。

登場する典型

マトロイド交差
- おそらくマトロイド交差のアルゴリズムを直接適用できることは少なくて、解法のあたりをつけるくらいしかできない
ラグランジュ緩和
- 競技プログラミングにおける動的計画法更新最適化まとめ（CHT, MongeDP, AlianDP, インラインDP, きたまさ法, 行列累乗） - はまやんはまやんはまやんのAlienDPに近い考え方。Convex Hull Trickとも近いかもしれない。
  - ある操作がK回までしかできないとき、代わりにその操作のコストをxと置いてxで二分探索し、その操作をちょうどK回行うときのxを調べる
  - 精度に気をつける (今回はたまたま適当にやってうまく行ったけど、いつもうまくいく保証は多分ない)

実装上の注意点

TLEに気をつける
xを実数にすると計算誤差が怖いので、(N+1)倍して整数演算にする

提出: #4130939 (Rust)

/**
 * Dinic's algorithm for maximum flow problem.
 * Verified by: yukicoder No.177 (http://yukicoder.me/submissions/148371)
 * Min-cut (the second element of max_flow's returned values) is not verified.
 */

#[derive(Clone)]
struct Edge<T> {
    to: usize,
    cap: T,
    rev: usize, // rev is the position of the reverse edge in graph[to]
}

struct Dinic<T> {
    graph: Vec<Vec<Edge<T>>>,
    iter: Vec<usize>,
    zero: T,
}

impl<T> Dinic<T>
    where T: Clone,
          T: Copy,
          T: Ord,
          T: std::ops::AddAssign,
          T: std::ops::SubAssign,
{
    fn bfs(&self, s: usize, level: &mut [Option<usize>]) {
        let n = level.len();
        for i in 0 .. n {
            level[i] = None;
        }
        let mut que = std::collections::VecDeque::new();
        level[s] = Some(0);
        que.push_back(s);
        while let Some(v) = que.pop_front() {
            for e in self.graph[v].iter() {
	        if e.cap > self.zero && level[e.to] == None {
	            level[e.to] = Some(level[v].unwrap() + 1);
	            que.push_back(e.to);
                }
            }
	}
    }
    /* search augment path by dfs.
     * if f == None, f is treated as infinity.
     */
    fn dfs(&mut self, v: usize, t: usize, f: Option<T>, level: &mut [Option<usize>]) -> T {
        if v == t {
            return f.unwrap();
        }
        while self.iter[v] < self.graph[v].len() {
            let i = self.iter[v];
            let e = self.graph[v][i].clone();
            if e.cap > self.zero && level[v] < level[e.to] {
                let newf = std::cmp::min(f.unwrap_or(e.cap), e.cap);
                let d = self.dfs(e.to, t, Some(newf), level);
                if d > self.zero {
                    self.graph[v][i].cap -= d;
                    self.graph[e.to][e.rev].cap += d;
                    return d;
                }
            }
            self.iter[v] += 1;
        }
        self.zero
    }
    pub fn new(n: usize, zero: T) -> Self {
        Dinic {
            graph: vec![Vec::new(); n],
            iter: vec![0; n],
            zero: zero,
        }
    }
    pub fn add_edge(&mut self, from: usize, to: usize, cap: T) {
        let added_from = Edge { to: to, cap: cap,
                            rev: self.graph[to].len() };
        let added_to = Edge { to: from, cap: self.zero,
                            rev: self.graph[from].len() };
        self.graph[from].push(added_from);
        self.graph[to].push(added_to);
    }
    pub fn max_flow(&mut self, s: usize, t: usize) -> (T, Vec<usize>) {
        let mut flow = self.zero;
        let n = self.graph.len();
        let mut level = vec![None; n];
        loop {
            self.bfs(s, &mut level);
            if level[t] == None {
                let ret = (0 .. n).filter(|&i| level[i] == None)
                    .collect();
                return (flow, ret);
            }
            self.iter.clear();
            self.iter.resize(n, 0);
            loop {
                let f = self.dfs(s, t, None, &mut level);
                if f <= self.zero { break; }
                flow += f;
            }
        }
    }
}

fn calc(n: usize, a: &[i64], b: &[i64], pairs: &[(usize, usize)], mid: i64)
        -> (i64, i64) {
    let bias = 5001; // should be > n
    let mut din = Dinic::new(2 + n, 0);
    let ni = n as i64;
    for i in 0 .. n {
        din.add_edge(0, 2 + i, (ni - b[i]) * bias);
        din.add_edge(2 + i, 1, (ni - a[i]) * bias);
        let (u, v) = pairs[i];
        din.add_edge(2 + u, 2 + v, mid);
    }
    let (_cost, t_side) = din.max_flow(0, 1);
    let mut kind = vec![0; n];
    for t in t_side {
        if t >= 2 {
            kind[t - 2] = 1;
        }
    }
    let mut used = 0;
    for i in 0 .. n {
        let (u, v) = pairs[i];
        if (kind[u], kind[v]) == (0, 1) { used += 1; }
    }
    let mut realcost = 0;
    for i in 0 .. n {
        if kind[i] == 1 {
            realcost += b[i];
        } else {
            realcost += a[i];
        }
    }
    (realcost, used)
}

fn solve() {
    let out = std::io::stdout();
    let mut out = BufWriter::new(out.lock());
    macro_rules! puts {
        ($format:expr) => (write!(out,$format).unwrap());
        ($format:expr, $($args:expr),+) => (write!(out,$format,$($args),*).unwrap())
    }
    input! {
        n: usize,
        k: usize,
        a: [i64; n],
        b: [i64; n],
    }
    let mut fail = -1;
    let mut pass: i64 = 5001 * 5001;
    let mut pairs = vec![(0, 0); n];
    for i in 0 .. n {
        pairs[a[i] as usize - 1].0 = i;
        pairs[b[i] as usize - 1].1 = i;
    }
    while pass - fail > 1 {
        let mid = (pass + fail) / 2;
        let (_cost, pair) = calc(n, &a, &b, &pairs, mid);
        if pair <= (n - k) as i64 {
            pass = mid;
        } else {
            fail = mid;
        }
    }
    let (cost, _) = calc(n, &a, &b, &pairs, pass);
    puts!("{}\n", cost);
}

まとめ

解説とはやり方が違う。この記事の解法は完全に牛刀っぽい。ラグランジュ緩和は結構面白いテクニックなので、これを読んだ人は習得するといいことがあるかもしれません。