Solving MAXSAT and saying a few words about it.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 

1200 lines
36 KiB

// vim: ts=4 sw=4 et
use crate::clause;
use core::fmt::{Display, Formatter};
use core::mem;
use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Shl};
use std::collections::{HashMap, HashSet};
use nalgebra::{DMatrix, DVector, Vector};
pub type Fl = f64;
#[allow(non_upper_case_globals)]
pub static DEFAULT_α: Fl = 0.25;
#[allow(non_upper_case_globals)]
pub static DEFAULT_ε: Fl = 1e-6;
pub type Variable = crate::variable::Variable<Fl>;
#[derive(Clone, Debug)]
pub struct Term {
coeff: Fl,
var: Variable,
}
impl Term {
pub fn get_var(&self) -> &Variable {
&self.var
}
pub fn evaluate(&self) -> Option<Fl> {
self.var.get_assignment().map(|x| x * self.coeff)
}
pub fn from_var(var: Variable, coeff: Fl) -> Term {
Term { var, coeff }
}
pub fn get_coeff(&self) -> Fl {
self.coeff
}
pub fn set_coeff(&mut self, coeff: Fl) {
self.coeff = coeff;
}
pub fn scale_coeff(&mut self, scale: Fl) {
self.coeff *= scale;
}
pub fn add_coeff(&mut self, addend: Fl) {
self.coeff += addend;
}
}
impl Add for Term {
type Output = Expr;
fn add(self, other: Term) -> Self::Output {
let mut expr = Expr::new_zero();
expr.add_term(self);
expr.add_term(other);
expr
}
}
impl Add<&Variable> for Term {
type Output = Expr;
fn add(self, other: &Variable) -> Self::Output {
let mut expr = Expr::new_zero();
expr.add_term(self);
expr.add_term(Self::from(other));
expr
}
}
impl Neg for Term {
type Output = Self;
fn neg(self) -> Self {
Self {
coeff: -self.coeff,
..self
}
}
}
impl Mul<Fl> for Term {
type Output = Term;
fn mul(mut self, factor: Fl) -> Self::Output {
self *= factor;
self
}
}
impl MulAssign<Fl> for Term {
fn mul_assign(&mut self, factor: Fl) {
self.scale_coeff(factor);
}
}
impl Shl<Expr> for Term {
type Output = LeqInequality;
fn shl(self, rhs: Expr) -> LeqInequality {
Expr::from(self) << rhs
}
}
impl Shl<Fl> for Term {
type Output = LeqInequality;
fn shl(self, rhs: Fl) -> LeqInequality {
Expr::from(self) << rhs
}
}
impl Display for Term {
fn fmt(&self, fmt: &mut Formatter) -> core::fmt::Result {
write!(fmt, "{}{}", self.coeff, self.var.get_name())
}
}
impl From<&Variable> for Term {
fn from(var: &Variable) -> Self {
Self {
coeff: 1.0,
var: var.clone(),
}
}
}
pub fn combine_like_terms<I: IntoIterator<Item = Term>>(iter: I) -> impl Iterator<Item = Term> {
let mut map = HashMap::new();
for term in iter {
map.entry(term.var.clone())
.and_modify(|c| *c += term.coeff)
.or_insert(term.coeff);
}
map.into_iter().map(|(var, coeff)| Term { coeff, var })
}
#[derive(Clone, Debug)]
pub struct Expr {
terms: Vec<Term>,
offset: Fl,
}
impl Expr {
pub fn new_zero() -> Expr {
Expr {
terms: Vec::new(),
offset: 0.0,
}
}
pub fn num_terms(&self) -> usize {
return self.terms.len();
}
pub fn iter_terms(&self) -> impl Iterator<Item = &Term> {
self.terms.iter()
}
pub fn drain_terms(&mut self) -> impl Iterator<Item = Term> {
let mut newlist = Vec::new();
mem::swap(&mut self.terms, &mut newlist);
newlist.into_iter()
}
pub fn add_term(&mut self, term: Term) {
self.terms.push(term);
}
pub fn get_offset(&self) -> Fl {
self.offset
}
pub fn set_offset(&mut self, offset: Fl) {
self.offset = offset;
}
pub fn add_offset(&mut self, addend: Fl) {
self.offset += addend;
}
pub fn scale(&mut self, factor: Fl) {
for term in &mut self.terms {
*term *= factor;
}
self.offset *= factor;
}
pub fn get_variables(&self) -> HashSet<Variable> {
self.terms.iter().map(|t| t.get_var()).cloned().collect()
}
pub fn add_vars_to(&self, set: &mut HashSet<Variable>) {
self.terms.iter().for_each(|t| {
set.insert(t.get_var().clone());
});
}
}
impl Add<Fl> for Expr {
type Output = Self;
fn add(mut self, addend: Fl) -> Self {
self += addend;
self
}
}
impl AddAssign<Fl> for Expr {
fn add_assign(&mut self, addend: Fl) {
self.add_offset(addend);
}
}
impl Add<Term> for Expr {
type Output = Self;
fn add(mut self, addend: Term) -> Self {
self += addend;
self
}
}
impl AddAssign<Term> for Expr {
fn add_assign(&mut self, addend: Term) {
self.add_term(addend);
}
}
impl Add<&Variable> for Expr {
type Output = Expr;
fn add(mut self, addend: &Variable) -> Self::Output {
self += addend;
self
}
}
impl AddAssign<&Variable> for Expr {
fn add_assign(&mut self, addend: &Variable) {
self.add_term(Term::from(addend));
}
}
impl Shl<Expr> for Expr {
type Output = LeqInequality;
fn shl(self, rhs: Expr) -> LeqInequality {
LeqInequality { lhs: self, rhs }
}
}
impl Shl<Fl> for Expr {
type Output = LeqInequality;
fn shl(self, rhs: Fl) -> LeqInequality {
LeqInequality {
lhs: self,
rhs: Expr::from(rhs),
}
}
}
impl Display for Expr {
fn fmt(&self, fmt: &mut Formatter) -> core::fmt::Result {
for term in &self.terms {
write!(fmt, "{} + ", term)?;
}
write!(fmt, "{}", self.offset)
}
}
impl From<Term> for Expr {
fn from(t: Term) -> Self {
Expr {
terms: vec![t],
offset: 0.0,
}
}
}
impl From<Fl> for Expr {
fn from(offset: Fl) -> Self {
Expr {
terms: Vec::new(),
offset,
}
}
}
impl From<&Variable> for Expr {
fn from(var: &Variable) -> Self {
Self::from(Term::from(var))
}
}
#[derive(Debug)]
pub struct NewVariableCounter(usize);
impl NewVariableCounter {
pub fn new() -> Self {
NewVariableCounter(0)
}
pub fn get_new_var(&mut self) -> Variable {
self.0 += 1;
Variable::from_name(format!("_s{}", self.0))
}
}
#[derive(Clone, Debug)]
pub struct LeqInequality {
lhs: Expr,
rhs: Expr,
}
impl LeqInequality {
/// Algebraically moves all terms to the left side, and the offset to the right side.
pub fn leftify(&mut self) {
for term in self.rhs.drain_terms() {
self.lhs += -term;
}
self.rhs += self.lhs.get_offset();
self.lhs.set_offset(0.0);
}
/// Indicates whether or not all variables are on the left and the left offset is 0.0.
pub fn is_left(&self) -> bool {
self.rhs.num_terms() == 0 && self.lhs.get_offset() == 0.0
}
/// Replaces the inequality with an equality by introducing a slack variable.
pub fn slacken(mut self, counter: &mut NewVariableCounter) -> Equality {
if !self.is_left() {
self.leftify()
}
self.lhs
.add_term(Term::from_var(counter.get_new_var(), 1.0));
let mut eq = Equality::new(self.rhs.get_offset());
for term in self.lhs.terms.into_iter() {
eq.insert_variable(term.get_var().clone(), term.get_coeff());
}
eq
}
pub fn add_vars_to(&self, set: &mut HashSet<Variable>) {
self.lhs.add_vars_to(set);
self.rhs.add_vars_to(set);
}
}
impl Display for LeqInequality {
fn fmt(&self, fmt: &mut Formatter) -> core::fmt::Result {
write!(fmt, "{} ≤ {}", self.lhs, self.rhs)
}
}
#[derive(Clone, Debug)]
pub struct LinearCombination {
dir: HashMap<Variable, Fl>,
}
impl LinearCombination {
pub fn new() -> Self {
Self {
dir: HashMap::new(),
}
}
pub fn get_coeff(&self, v: &Variable) -> Fl {
self.dir.get(v).map(|x| *x).unwrap_or(0.0)
}
pub fn add_vars_to(&self, set: &mut HashSet<Variable>) {
for var in self.dir.keys() {
set.insert(var.clone());
}
}
pub fn insert_variable(&mut self, var: Variable, coeff: Fl) {
self.dir.insert(var, coeff);
}
}
impl Display for LinearCombination {
fn fmt(&self, fmt: &mut Formatter) -> core::fmt::Result {
let mut first = true;
for (var, val) in self.dir.iter() {
if first {
if *val >= 0.0 {
write!(fmt, "{}{}", val, var.get_name())?;
} else {
write!(fmt, "-{}{}", -val, var.get_name())?;
}
first = false;
} else {
if *val >= 0.0 {
write!(fmt, " + {}{}", val, var.get_name())?;
} else {
write!(fmt, " - {}{}", -val, var.get_name())?;
}
}
}
Ok(())
}
}
#[derive(Clone, Debug)]
pub struct Equality {
lhs: LinearCombination,
rhs: Fl,
}
impl Equality {
pub fn new(offset: Fl) -> Self {
Equality {
lhs: LinearCombination::new(),
rhs: offset,
}
}
pub fn insert_variable(&mut self, var: Variable, coeff: Fl) {
self.lhs.insert_variable(var, coeff);
}
pub fn get_coeff(&self, v: &Variable) -> Fl {
self.lhs.get_coeff(v)
}
pub fn add_vars_to(&self, set: &mut HashSet<Variable>) {
self.lhs.add_vars_to(set);
}
pub fn get_rhs(&self) -> Fl {
self.rhs
}
}
impl Display for Equality {
fn fmt(&self, fmt: &mut Formatter) -> core::fmt::Result {
write!(fmt, "{} = {}", self.lhs, self.rhs)
}
}
#[derive(Debug, Clone)]
#[allow(non_snake_case)]
pub struct StandardFormLP {
pub c: DVector<Fl>,
pub A: DMatrix<Fl>,
pub b: DVector<Fl>,
}
impl StandardFormLP {
pub fn from_ineqs_and_objective(
ineqs: impl IntoIterator<Item = LeqInequality>,
objective: LinearCombination,
) -> (Vec<Variable>, StandardFormLP) {
let mut counter = NewVariableCounter::new();
let constraints = ineqs
.into_iter()
.map(|c| c.slacken(&mut counter))
.collect::<Vec<_>>();
let mut varset = HashSet::new();
objective.add_vars_to(&mut varset);
constraints.iter().for_each(|c| c.add_vars_to(&mut varset));
let vars = varset.into_iter().collect::<Vec<_>>(); // Order needs to be deterministic
let ncstr = constraints.len();
let nvars = vars.len();
let mut c = DVector::zeros(nvars);
#[allow(non_snake_case)]
let mut A = DMatrix::zeros(ncstr, nvars);
let mut b = DVector::zeros(ncstr);
for (i, constraint) in constraints.iter().enumerate() {
b[i] = constraint.get_rhs();
for (j, var) in vars.iter().enumerate() {
A[(i, j)] = constraint.get_coeff(var);
}
}
for (j, var) in vars.iter().enumerate() {
c[j] = objective.get_coeff(var);
}
(vars, StandardFormLP { c, A, b })
}
pub fn from_ineqs_and_objective_with_varset(
varset: impl IntoIterator<Item = Variable>,
ineqs: impl IntoIterator<Item = LeqInequality>,
objective: LinearCombination,
) -> StandardFormLP {
let mut counter = NewVariableCounter::new();
let constraints = ineqs
.into_iter()
.map(|c| c.slacken(&mut counter))
.collect::<Vec<_>>();
let vars = varset.into_iter().collect::<Vec<_>>();
let ncstr = constraints.len();
let nvars = vars.len();
let mut c = DVector::zeros(nvars);
#[allow(non_snake_case)]
let mut A = DMatrix::zeros(ncstr, nvars);
let mut b = DVector::zeros(ncstr);
for (i, constraint) in constraints.iter().enumerate() {
b[i] = constraint.get_rhs();
for (j, var) in vars.iter().enumerate() {
A[(i, j)] = constraint.get_coeff(var);
}
}
for (j, var) in vars.iter().enumerate() {
c[j] = objective.get_coeff(var);
}
StandardFormLP { c, A, b }
}
/// Converts the constrained maximization problem to one where we instead try to find a point
/// in the feasible region of the problem. When the last variable in the problem is zero, the
/// remaining ones are guaranteed to be a representation of a point in the feasible region,
/// with the variables in the same order they appeared in the original problem.
///
/// The all-ones vector is in the feasible region of the problem returned by this function.
pub fn to_feasible_region_problem(self) -> StandardFormLP {
let x0 = DVector::from_element(self.A.ncols(), 1.0);
let v = self.b.clone() - &self.A * x0;
let lastcol = self.A.ncols(); // or self.vars.len()
#[allow(non_snake_case)]
let mut A = self.A.insert_column(lastcol, 0.0);
A.set_column(lastcol, &v);
let mut c = DVector::zeros(lastcol + 1);
c[lastcol] = -1.0; // Try to minimize u
Self { A, c, ..self }
}
pub fn set_frp_from(&mut self, from: &Self) {
let nvars = from.c.len();
self.b.copy_from(&from.b);
self.c.rows_mut(0, nvars).fill(0.0);
self.c[nvars] = -1.0;
self.A.columns_mut(0, nvars).copy_from(&from.A);
self.A.column_mut(nvars).copy_from(&from.b);
for (i, row) in from.A.row_iter().enumerate() {
self.A[(nvars, i)] -= row.sum();
}
}
}
impl Display for StandardFormLP {
fn fmt(&self, fmt: &mut Formatter) -> core::fmt::Result {
write!(fmt, "variables: ")?;
writeln!(fmt, "")?;
writeln!(fmt, "A:")?;
write!(fmt, "{}", &self.A)?;
writeln!(fmt, "b:")?;
write!(fmt, "{}", &self.b)?;
writeln!(fmt, "c:")?;
write!(fmt, "{}", &self.c)?;
Ok(())
}
}
#[derive(Debug, Clone)]
#[allow(non_snake_case)]
pub struct HomogenousLP {
pub c: DVector<Fl>,
pub A: DMatrix<Fl>,
}
impl HomogenousLP {
pub fn project_point_inplace<Ul, Ur>(
point: &Vector<Fl, nalgebra::Dynamic, Ul>,
out: &mut Vector<Fl, nalgebra::Dynamic, Ur>,
) where
Ul: nalgebra::storage::Storage<Fl, nalgebra::Dynamic>,
Ur: nalgebra::storage::StorageMut<Fl, nalgebra::Dynamic>,
{
let n = point.len();
let a_0 = ((n + 1) as Fl).recip();
let factor = (a_0 * (point.sum() / a_0 + 1.0)).recip();
let mut sum = 0.0;
for (i, x_i) in point.iter().enumerate() {
out[i] = factor * *x_i;
sum += out[i];
}
out[n] = 1.0 - sum;
}
pub fn project_point(point: &DVector<Fl>) -> DVector<Fl> {
let mut rv = DVector::zeros(point.len() + 1);
Self::project_point_inplace(point, &mut rv);
rv
}
pub fn project_point_inv_inplace<Ul, Ur>(
point: &Vector<Fl, nalgebra::Dynamic, Ul>,
out: &mut Vector<Fl, nalgebra::Dynamic, Ur>,
) where
Ul: nalgebra::storage::Storage<Fl, nalgebra::Dynamic>,
Ur: nalgebra::storage::StorageMut<Fl, nalgebra::Dynamic>,
{
let n = point.len() - 1;
let a_0 = (point.len() as Fl).recip();
let scale = a_0 / point[n];
out.copy_from(&point.rows(0, n));
*out *= scale;
}
pub fn project_point_inv(point: &DVector<Fl>) -> DVector<Fl> {
let mut rv = DVector::zeros(point.len() - 1);
Self::project_point_inv_inplace(point, &mut rv);
rv
}
#[allow(non_snake_case)]
pub fn project_system(A: &DMatrix<Fl>, b: &DVector<Fl>, out: &mut DMatrix<Fl>) {
let n = A.ncols();
out.columns_mut(0, A.ncols()).copy_from(&(A));
out.column_mut(A.ncols()).copy_from(&(-((n + 1) as Fl) * b));
}
pub fn project_from_sflp_inplace(&mut self, sflp: &StandardFormLP) {
let len = sflp.c.len();
Self::project_system(&sflp.A, &sflp.b, &mut self.A);
self.c.rows_mut(0, len).copy_from(&sflp.c);
}
pub fn project_b_from(&mut self, b: &DVector<Fl>) {
let n1 = self.A.ncols();
self.A.column_mut(n1 - 1).copy_from(&(n1 as f64 * b));
}
pub fn project_from_sflp(sflp: StandardFormLP) -> Self {
#[allow(non_snake_case)]
let mut A = DMatrix::zeros(sflp.A.nrows(), sflp.A.ncols() + 1);
Self::project_system(&sflp.A, &sflp.b, &mut A);
let len = sflp.c.len();
let c = sflp.c.insert_row(len, 0.0);
HomogenousLP { A, c }
}
pub fn copy_from(&mut self, other: &Self) {
self.A.copy_from(&other.A);
self.c.copy_from(&other.c);
}
}
impl Display for HomogenousLP {
fn fmt(&self, fmt: &mut Formatter) -> core::fmt::Result {
write!(fmt, "variables: ")?;
writeln!(fmt, "")?;
writeln!(fmt, "A:")?;
write!(fmt, "{}", &self.A)?;
writeln!(fmt, "c:")?;
write!(fmt, "{}", &self.c)?;
Ok(())
}
}
#[derive(Debug, Clone)]
#[allow(non_snake_case)]
pub struct Karmarkar {
pub problem: HomogenousLP,
pub guess: DVector<Fl>,
ε: Fl,
nvars: usize,
nconstr: usize,
αr: Fl,
: Fl,
a0: Fl,
Dc: DVector<Fl>,
B: DMatrix<Fl>,
BT: DMatrix<Fl>,
BBT: DMatrix<Fl>,
BTBBT: DMatrix<Fl>,
BTBBTB: DMatrix<Fl>,
cp: DVector<Fl>,
}
impl Karmarkar {
pub fn from_feasible_value_αε(
problem: HomogenousLP,
value: DVector<Fl>,
α: Fl,
ε: Fl,
) -> Option<Karmarkar> {
if value.len() != problem.c.len() {
return None;
}
let n = value.len();
let r = ((n * (n - 1)) as Fl).sqrt().recip();
let a0 = ((n + 1) as Fl).recip();
#[allow(non_snake_case)]
let Ashape = problem.A.shape();
Some(Karmarkar {
problem: problem,
guess: value,
ε,
nvars: Ashape.1,
nconstr: Ashape.0,
αr: α * r,
: Ashape.1 as Fl / α,
a0,
Dc: DVector::zeros(Ashape.1),
B: DMatrix::zeros(Ashape.0 + 1, Ashape.1),
BT: DMatrix::zeros(Ashape.1, Ashape.0 + 1),
BBT: DMatrix::zeros(Ashape.0 + 1, Ashape.0 + 1),
BTBBT: DMatrix::zeros(Ashape.1, Ashape.0 + 1),
BTBBTB: DMatrix::zeros(Ashape.1, Ashape.1),
cp: DVector::zeros(Ashape.1),
})
}
pub fn from_feasible_value(problem: HomogenousLP, value: DVector<Fl>) -> Option<Karmarkar> {
Self::from_feasible_value_αε(problem, value, DEFAULT_α, DEFAULT_ε)
}
pub fn iterate(&mut self) -> IterationResult {
use IterationResult::*;
let problem = &self.problem;
#[cfg(debug_karmarkar)]
{
print!("infeasiblity: {}", &problem.A * &self.guess);
}
// First transform c, giving Dc
self.Dc.copy_from(&problem.c);
self.Dc.component_mul_assign(&self.guess);
// Then project Dc into the null space of B (A augmented with the simplex constraint)
{
#[allow(non_snake_case)]
let mut Bslice = self.B.rows_mut(0, self.nconstr);
Bslice.copy_from(&problem.A);
for (i, mut col) in Bslice.column_iter_mut().enumerate() {
col.copy_from(&(&col * self.guess[i]));
}
self.B.row_mut(self.nconstr).fill(1.0);
}
self.B.transpose_to(&mut self.BT);
self.B.mul_to(&self.BT, &mut self.BBT);
if !self.BBT.try_inverse_mut() {
return Failed();
}
self.BT.mul_to(&self.BBT, &mut self.BTBBT);
self.BTBBT.mul_to(&self.B, &mut self.BTBBTB);
for i in 0..self.nvars {
self.BTBBTB[(i, i)] -= 1.0;
}
self.BTBBTB.mul_to(&self.Dc, &mut self.cp);
// Now normalize the result and step in that direction
self.cp.normalize_mut();
self.cp *= -self.αr;
self.cp.add_scalar_mut(self.a0);
// And finally, transform back into the problem space.
self.cp.component_mul_assign(&self.guess);
let scale = 1.0 / self.cp.sum();
self.cp *= scale;
if self.cp.iter().all(|x| !x.is_nan()) {
let est_dist = self.guess.metric_distance(&self.cp) * self.;
self.guess.copy_from(&self.cp);
Distance(est_dist)
} else {
Failed()
}
}
/// Iterates until either the iteration fails (usually indicating the edge of the simplex, and
/// therefore a solution) or the desired precision is reached. If `unproject` is true, try to
/// get accurate to within ε in the unprojected space, assuming the problem this instance is
/// solving is the result of projecting a lower-dimensional one. Otherwise, just try to be
/// accurate within ε in this space.
pub fn iterate_until_failed_or_precise(&mut self, unproject: bool) {
use IterationResult::*;
loop {
match self.iterate() {
Failed() => break,
Distance(d) => {
let tolerance = if unproject {
self.guess[self.guess.len() - 1] * self.ε
} else {
self.ε
};
if d < tolerance {
break;
}
}
_ => {}
}
}
}
}
#[derive(Debug, Clone, Copy)]
pub enum IterationResult {
Distance(Fl),
Failed(),
Indeterminate(),
Unbounded(),
}
/// Defines a DSL which allows for the declaration of many variables at once. See the unit tests
/// for examples.
#[macro_export]
macro_rules! linear_vars {
[] => {};
[$varname:ident$(,)?] => {
let $varname = $crate::karmarkar::Variable::from_name(stringify!($varname).into());
};
[$varname:ident, $($varnames:ident),+] => {
let $varname = $crate::karmarkar::Variable::from_name(stringify!($varname).into());
linear_vars![$($varnames),+]
};
}
// Solves linear programs where some of the variables are constrained to either 0 or 1.
pub struct BMIPSolver {
freevars: Vec<Variable>, // Varaibles which aren't constrained except by the system.
bvars: Vec<Variable>, // Variables which are constrained to 0 or 1.
best_solution: Option<(Fl, DVector<Fl>)>, // value of objective function, solution vector
problems: Vec<StandardFormLP>,
projected_problems: Vec<Karmarkar>,
feasibility_problems: Vec<StandardFormLP>,
projected_feasibility_problems: Vec<Karmarkar>,
objective_offsets: Vec<Fl>,
feasible_points: Vec<DVector<Fl>>, // From here down, they're locals to avoid reallocation
}
impl BMIPSolver {
fn remove_var(from: StandardFormLP) -> StandardFormLP {
let nvars = from.A.ncols();
StandardFormLP {
A: from.A.remove_column(nvars - 1),
b: from.b,
c: from.c.remove_row(nvars - 1),
}
}
pub fn nvars(&self) -> usize {
self.bvars.len() + self.freevars.len()
}
pub fn from_ineqs_and_objective(
ineqs: impl IntoIterator<Item = LeqInequality>,
constrained_vars: impl IntoIterator<Item = Variable>,
objective: LinearCombination,
) -> BMIPSolver {
let constr_varset = constrained_vars.into_iter().collect::<HashSet<_>>();
let ineqs = ineqs.into_iter().collect::<Vec<_>>();
let mut varset = HashSet::new();
objective.add_vars_to(&mut varset);
ineqs.iter().for_each(|c| c.add_vars_to(&mut varset));
let freevars = varset
.difference(&constr_varset)
.cloned()
.collect::<Vec<_>>();
let bvars = constr_varset.into_iter().collect::<Vec<_>>();
let mut problems = Vec::new();
let mut projected_problems = Vec::new();
let mut feasibility_problems = Vec::new();
let mut projected_feasibility_problems = Vec::new();
let mut feasible_points = Vec::new();
let mut objective_offsets = Vec::new();
let problem = StandardFormLP::from_ineqs_and_objective_with_varset(
freevars.iter().cloned().chain(bvars.iter().cloned()),
ineqs,
objective,
);
let proj_problem = Karmarkar::from_feasible_value(
HomogenousLP::project_from_sflp(problem.clone()),
DVector::zeros(problem.A.ncols() + 1),
)
.unwrap();
let feas_problem = problem.clone().to_feasible_region_problem();
let proj_feas_problem = Karmarkar::from_feasible_value(
HomogenousLP::project_from_sflp(feas_problem.clone()),
DVector::zeros(problem.A.ncols() + 2),
)
.unwrap();
problems.push(problem);
projected_problems.push(proj_problem);
feasibility_problems.push(feas_problem);
projected_feasibility_problems.push(proj_feas_problem);
feasible_points.push(DVector::zeros(0)); // The first one should never end up being used.
objective_offsets.push(0.0);
for i in 0..bvars.len() {
let problem = Self::remove_var(problems[i].clone());
let proj_problem = Karmarkar::from_feasible_value(
HomogenousLP::project_from_sflp(problem.clone()),
DVector::zeros(problem.A.ncols() + 1),
)
.unwrap();
let feas_problem = problem.clone().to_feasible_region_problem();
let proj_feas_problem = Karmarkar::from_feasible_value(
HomogenousLP::project_from_sflp(feas_problem.clone()),
DVector::zeros(problem.A.ncols() + 2),
)
.unwrap();
let feasible_point = DVector::zeros(problem.A.ncols() + 1);
problems.push(problem);
projected_problems.push(proj_problem);
feasibility_problems.push(feas_problem);
projected_feasibility_problems.push(proj_feas_problem);
feasible_points.push(feasible_point);
}
BMIPSolver {
freevars,
bvars,
best_solution: None,
problems,
feasibility_problems,
projected_problems,
projected_feasibility_problems,
objective_offsets,
feasible_points,
}
}
fn solve_at_depth(&mut self, nassigned: usize) -> Option<Fl> {
let nvars = self.nvars();
self.feasible_points[nassigned].fill(1.0);
HomogenousLP::project_point_inv_inplace(
&self.feasible_points[nassigned],
&mut self.projected_feasibility_problems[nassigned].guess,
);
self.projected_feasibility_problems[nassigned].iterate_until_failed_or_precise(true);
HomogenousLP::project_point_inv_inplace(
&self.projected_feasibility_problems[nassigned].guess,
&mut self.feasible_points[nassigned],
);
if self.feasible_points[nassigned][nvars] > DEFAULT_ε {
// The problem as given isn't feasible
return None;
}
HomogenousLP::project_point_inplace(
&self.projected_feasibility_problems[nassigned]
.guess
.rows(0, self.nvars()),
&mut self.projected_problems[nassigned].guess,
);
self.projected_problems[nassigned].iterate_until_failed_or_precise(true);
HomogenousLP::project_point_inv_inplace(
&self.projected_problems[nassigned].guess,
&mut self.feasible_points[nassigned].rows_mut(0, nvars),
);
Some(
self.problems[nassigned]
.c
.dot(&self.feasible_points[nassigned].rows_mut(0, nvars))
+ *self.objective_offsets.last().unwrap(),
)
}
// `true` if we can prune
fn can_prune(&mut self, nassigned: usize) -> bool {
match self.best_solution {
None => false,
Some((best, _)) => match self.solve_at_depth(nassigned) {
None => true,
Some(value) => value < best,
},
}
}
fn record_maximum(&mut self, value: Fl) {
let best_solution = match &mut self.best_solution {
x @ None => {
*x = Some((value, DVector::zeros(self.bvars.len())));
x.as_mut().unwrap()
}
Some(ref mut x) => x,
};
for (i, var) in self.bvars.iter().enumerate() {
best_solution.1[i] = var.get_assignment().unwrap();
}
best_solution.0 = value;
}
fn accumulate_maximum(&mut self, nassigned: usize) {
if let Some(value) = self.solve_at_depth(nassigned) {
match &mut self.best_solution {
None => {
self.record_maximum(value);
}
Some((best, _)) => {
if value > *best {
self.record_maximum(value);
}
}
}
}
}
fn assign(&mut self, nassigned: usize, v: Fl) {
let nvars = self.nvars();
self.bvars.iter().nth_back(nassigned).unwrap().assign(v);
let newb = &self.problems[nassigned].b - v * &self.problems[nassigned].A.column(nvars - 1);
self.problems[nassigned + 1].b.copy_from(&newb);
self.projected_problems[nassigned + 1]
.problem
.project_b_from(&newb);
self.feasibility_problems[nassigned + 1].b.copy_from(&newb);
self.projected_feasibility_problems[nassigned + 1]
.problem
.project_b_from(&newb);
self.objective_offsets
.push(self.problems[nassigned].c[nvars - 1] * v);
}
fn unassign(&mut self, nassigned: usize) {
self.bvars.iter().nth_back(nassigned).unwrap().unassign();
self.objective_offsets.pop();
}
fn search_at_depth(&mut self, nassigned: usize) {
if nassigned == self.bvars.len() {
self.accumulate_maximum(nassigned);
} else if !self.can_prune(nassigned) {
self.assign(nassigned, 0.0);
self.search_at_depth(nassigned + 1);
self.unassign(nassigned);
self.assign(nassigned, 1.0);
self.search_at_depth(nassigned + 1);
self.unassign(nassigned);
}
}
pub fn assign_from_best_solution(&mut self) -> Option<Fl> {
match &self.best_solution {
None => None,
Some(sol) => {
for (i, var) in self.bvars.iter().enumerate() {
var.assign(sol.1[i]);
}
Some(sol.0)
}
}
}
pub fn solve(&mut self) -> Option<Fl> {
self.search_at_depth(0);
self.assign_from_best_solution()
}
}
fn clause_to_inequality(
clause: &clause::Clause,
varmap: &HashMap<clause::Variable, Variable>,
) -> LeqInequality {
let mut expr = Expr::from(0.0);
for literal in clause.iter() {
if literal.is_neg() {
expr += 1.0;
expr += -varmap.get(&literal.get_variable()).unwrap();
} else {
expr += varmap.get(&literal.get_variable()).unwrap();
}
}
expr << 1.0
}
/// Accepts a ClauseList and a list of weights. Entries correspond between the lists. If the weight
/// is `None`, the clause is hard; otherwise, it's soft with the given weight. If the weight list
/// is too short, the clauses with no weight assumed to be hard.
pub fn clause_list_to_constraints<T: Copy>(
clause_list: &mut clause::ClauseList,
weights: &Vec<Option<T>>,
) -> (
Vec<LeqInequality>,
LinearCombination,
HashMap<clause::Variable, Variable>,
)
where
Fl: From<T>,
{
let mut weightmap = HashMap::new();
let mut counter = 0;
for (i, clause) in clause_list.iter_mut().enumerate() {
if let Some(Some(weight)) = weights.get(i) {
counter += 1;
let b = clause::Variable::from_name(format!("_b{}", counter));
weightmap.insert(b.clone(), Fl::from(*weight));
clause.push(b.get_pos_literal());
}
}
let mut varmap = HashMap::new();
let mut objective = LinearCombination::new();
let vars = clause_list.get_vars();
for var in &vars {
let lvar = Variable::from_name(var.get_name().clone());
weightmap
.get(&var)
.map(|x| objective.insert_variable(lvar.clone(), *x));
varmap.insert(var.clone(), lvar);
}
let mut constraints = Vec::new();
for clause in clause_list.iter() {
constraints.push(clause_to_inequality(&clause, &varmap));
}
(constraints, objective, varmap)
}
pub fn ilp_maxsat<T: Copy>(mut clause_list: clause::ClauseList, weights: &Vec<Option<T>>)
where
Fl: From<T>,
{
let (ineqs, objective, varmap) = clause_list_to_constraints(&mut clause_list, weights);
let varset = varmap.values().cloned();
let mut solver = BMIPSolver::from_ineqs_and_objective(ineqs, varset, objective);
solver.solve();
for (bv, lv) in varmap {
lv.get_assignment().map(|x| x > 0.5).map(|x| bv.assign(x));
}
}
#[cfg(test)]
mod test {
use super::*;
// Designed to be a --nocapture test for human analysis. Pretty much just for debugging.
#[test]
fn inequality_test() {
linear_vars![x, y];
let mut constraints = Vec::<LeqInequality>::new();
constraints.push(&x * 50.0 + &y * 24.0 << 2400.0);
constraints.push(&x * 30.0 + &y * 33.0 << 2100.0);
constraints.push(-&x << -45.0);
constraints.push(-&y << -5.0);
let mut objective = LinearCombination::new();
objective.insert_variable(x.clone(), 1.0);
objective.insert_variable(y.clone(), 1.0);
println!("Maximize {}", objective);
println!("subject to");
for constraint in constraints.iter() {
println!("{}", constraint);
}
let (vars, sf) = StandardFormLP::from_ineqs_and_objective(constraints, objective);
println!("");
println!("{}", sf);
let frp = sf.clone().to_feasible_region_problem();
println!("");
println!("{}", frp);
let frp_feasible_point = DVector::from_element(frp.A.ncols(), 1.0);
println!("FRP feasible point: {}", frp_feasible_point);
let proj_frp = HomogenousLP::project_from_sflp(frp.clone());
println!("");
println!("{}", proj_frp);
let proj_frp_feasible_point = HomogenousLP::project_point(&frp_feasible_point);
println!(
"FRP feasible point (projected): {}",
proj_frp_feasible_point
);
println!(
"unprojected FRP feasible point: {}",
HomogenousLP::project_point_inv(&proj_frp_feasible_point)
);
println!("");
let mut km =
Karmarkar::from_feasible_value(proj_frp.clone(), proj_frp_feasible_point).unwrap();
km.iterate_until_failed_or_precise(true);
let a0 = (km.guess.len() as Fl).recip();
print!("final x': {}", km.guess);
println!("sum of elements: {}", km.guess.sum());
println!(
"infeasiblilty in projected FRP: {}",
(&proj_frp.A * &km.guess)
);
println!(
"Ax': {}",
proj_frp.A.columns(0, proj_frp.A.ncols() - 1) * km.guess.rows(0, km.guess.len() - 1)
);
println!(
"b': {}",
proj_frp.A.column_iter().last().unwrap() * *km.guess.iter().last().unwrap()
);
println!("ratio: {}", a0 / km.guess[km.guess.len() - 1]);
let feasible_point = HomogenousLP::project_point_inv(&km.guess);
print!("feasible x: {}", feasible_point);
println!(
"infeasiblilty in FRP: {}",
(&frp.A * &feasible_point - &frp.b)
);
println!("Ax: {}", &frp.A * &feasible_point);
println!("b: {}", &frp.b);
let feasible_point = feasible_point.remove_row(km.guess.len() - 2);
let proj_feasible_point = HomogenousLP::project_point(&feasible_point);
let proj = HomogenousLP::project_from_sflp(sf.clone());
let mut km = Karmarkar::from_feasible_value(proj.clone(), proj_feasible_point).unwrap();
km.iterate_until_failed_or_precise(true);
print!("final x': {}", km.guess);
let solution = HomogenousLP::project_point_inv(&km.guess);
println!("vars: {:?}", vars);
print!("final x: {}", solution);
}
}