最適化 - 燃やす埋める問題

\(\def\I{\mathbb I}\) Project selection problem とも.

問題

\(n\) 個の Bool 変数 \(x_0, x_1, \ldots, x_{n-1}\) があるとする. 特に実装上ではこの添字 \(0,1,\ldots,n-1\) で変数を指す. これに対して次の目標関数を考える. \[\Theta(x ; b_0, b_1, c) = \sum_i b_0^i \I(x_i) + \sum_i b_1^i (1 - \I(x_i)) + \sum_i \sum_j c_{ij} \I(x_i) (1 - \I(x_j))\]

ここで

  • \(b_0, b_1\) はそれぞれ長さ \(n\) の整数列
  • \(c\)\(n \times n\) の整数行列
  • \(\I(x)\) は Bool 値 \(x\) に対するインディケータで真のとき \(1\), 偽のとき \(0\)

つまり,

  • \(b_0^i\)\(x_i\) が真のときのコスト
  • \(b_1^i\)\(x_i\) が偽のときのコスト
  • \(c_{ij}\)\(x_i\) が真で \(x_j\) が偽のときのコスト

与えられた \(b, c\) に対して \(\Theta\) を最小化する \(x\) の割当を求める問題を燃やす埋める問題という.

解法

最小カット問題に帰着する.

ライブラリの使い方

cost! または gain! マクロで制約を付け足していく.

/// サンプル
let n = 4;  // 使う変数の個数
let mut solver = MoyasuUmeruSolver::new(n);

solver += cost!(1; if 0); // 変数 0 が true のときのコストが1
solver += cost!(2; if not 0); // false のときのコストが2

solver += cost!(inf; if 1); // 無限大のコスト

solver += gain!(1; if not 1); // 逆に獲得できる利益. 内部ではマイナスにしてコストとして扱う

solver += cost!(100; if 0 and not 1);  // 0 が true, 1 が false のときのコスト
solver += cost!(inf; if not 0 and 1);  // 0 が false, 1 が true のときは無限のコスト

solver.min_cost();  // 最小化されたコスト

solver.max_gain();  // 逆に最大化した獲得量. 最小コストのマイナスを返す
use crate::algebra::group::*;
use crate::algebra::hyper::*;
use crate::graph::directed::dinic::*;

#[derive(Debug, Clone)]
pub struct MoyasuUmeruSolver {
    size: usize,
    cost1: Vec<Vec<Hyper<i64>>>,
    cost2: Vec<(usize, usize, Hyper<i64>)>,
}
impl MoyasuUmeruSolver {
    pub fn new(size: usize) -> Self {
        Self {
            size,
            cost1: vec![vec![Hyper::zero(); 2]; size],
            cost2: vec![],
        }
    }
    /// +cost when var `i` is `b`.
    fn add_constraint_node(&mut self, i: usize, b: bool, cost: Hyper<i64>) {
        let j = if b { 0 } else { 1 };
        self.cost1[i][j] = cost;
    }
    /// +cost when var `i` is true AND var `j` is false.
    fn add_constraint_edge(&mut self, i: usize, j: usize, cost: Hyper<i64>) {
        self.cost2.push((i, j, cost));
    }
    pub fn min_cost(&self) -> Hyper<i64> {
        let mut base: Hyper<i64> = Hyper::zero();
        let s = self.size;
        let t = self.size + 1;
        let mut g = vec![vec![]; self.size + 2];
        for i in 0..self.size {
            match (self.cost1[i][0], self.cost1[i][1]) {
                (b1, b2) if b1 > b2 => {
                    base = base + b2;
                    g[s].push((i, b1 - b2));
                }
                (b1, b2) if b1 < b2 => {
                    base = base + b1;
                    g[i].push((t, b2 - b1));
                }
                (b, _) => base = base + b,
            }
        }
        for &(i, j, cost) in self.cost2.iter() {
            g[j].push((i, cost));
        }
        base + Dinic::new(s, t, &g).maxflow()
    }
    pub fn max_gain(&self) -> Hyper<i64> {
        -self.min_cost()
    }
}
pub enum MoyasuUmeruCost {
    Node(usize, bool, Hyper<i64>),
    Edge(usize, usize, Hyper<i64>),
}
impl std::ops::AddAssign<MoyasuUmeruCost> for MoyasuUmeruSolver {
    fn add_assign(&mut self, constraint: MoyasuUmeruCost) {
        match constraint {
            MoyasuUmeruCost::Node(i, b, cost) => {
                self.add_constraint_node(i, b, cost);
            }
            MoyasuUmeruCost::Edge(i, j, cost) => {
                self.add_constraint_edge(i, j, cost);
            }
        }
    }
}
#[macro_export]
macro_rules! cost {
    (inf ; $( $rest:tt )*) => {
        cost!(@ Hyper::Inf ; $($rest)*)
    };
    ($cost:expr ; $( $rest:tt )*) => {
        cost!(@ Hyper::Real($cost) ; $($rest)*)
    };
    (@ $cost:expr ; if not $i:tt) => {
        MoyasuUmeruCost::Node($i, false, $cost)
    };
    (@ $cost:expr ; if not $i:tt and $j:tt) => {
        MoyasuUmeruCost::Edge($j, $i, $cost)
    };
    (@ $cost:expr ; if $i:tt and not $j:tt) => {
        MoyasuUmeruCost::Edge($i, $j, $cost)
    };
    (@ $cost:expr ; if $i:tt) => {
        MoyasuUmeruCost::Node($i, true, $cost)
    };
}
#[macro_export]
macro_rules! gain {
    (inf ; $( $rest:tt )*) => {
        cost!(@ Hyper::NegInf ; $($rest)*)
    };
    ($cost:expr ; $( $rest:tt )*) => {
        cost!(-$cost ; $($rest)*)
    }
}

#[cfg(test)]
mod test_umeru_moyasu {
    use crate::algebra::hyper::Hyper;
    use crate::opt::umeru_moyasu::{MoyasuUmeruCost, MoyasuUmeruSolver};

    #[test]
    fn test_simple() {
        let mut solver = MoyasuUmeruSolver::new(2);
        solver += cost!(5; if 0);
        solver += cost!(3; if not 0);
        assert_eq!(solver.min_cost(), Hyper::Real(3));

        solver += gain!(100; if 1);
        assert_eq!(solver.max_gain(), Hyper::Real(97));

        solver += cost!(inf; if 1 and not 0);
        assert_eq!(solver.max_gain(), Hyper::Real(95));
    }

    /// https://en.wikipedia.org/wiki/Max-flow_min-cut_theorem#Project_selection_problem
    #[test]
    fn test_wikipedia() {
        let mut solver = MoyasuUmeruSolver::new(6);
        // reward from projects
        solver += gain!(100; if 0);
        solver += gain!(200; if 1);
        solver += gain!(150; if 2);
        // cost for machines
        solver += cost!(200; if 3);
        solver += cost!(100; if 4);
        solver += cost!(50; if 5);
        // required machine for each projects
        solver += cost!(inf; if 0 and not 3);
        solver += cost!(inf; if 0 and not 4);
        solver += cost!(inf; if 1 and not 4);
        solver += cost!(inf; if 2 and not 5);

        assert_eq!(solver.max_gain(), Hyper::Real(200));
    }
}