最適化 - 燃やす埋める問題
\(\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_additive::*;
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));
}
}