Browse Source

Added (currently untested) affine scaling LP solver

karmarkar
Thomas Johnson 3 years ago
parent
commit
89eddd7cf6
  1. 483
      src/affine_scaling.rs
  2. 6
      src/main.rs
  3. 9
      src/variable.rs

483
src/affine_scaling.rs

@ -0,0 +1,483 @@
// vim: ts=4 sw=4 et
use core::fmt::{Display, Formatter};
use core::mem;
use core::ops::{Add, AddAssign, Mul, MulAssign, Neg};
use std::collections::{HashMap, HashSet};
use nalgebra::{DMatrix, DVector};
pub type Fl = f64;
// In case it wasn't clear, variable names and algorithm details are ripped straight from the
// Wikipedia page. Some day I'll go back and understand the algorithm in greater depth, I promise.
#[allow(non_upper_case_globals)]
pub static DEFAULT_β: Fl = 0.5;
#[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 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 Display for Term {
fn fmt(&self, fmt: &mut Formatter) -> core::fmt::Result {
write!(fmt, "{}{}", self.coeff, self.var.get_name())
}
}
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 Display for Expr {
fn fmt(&self, fmt: &mut Formatter) -> core::fmt::Result {
for term in &self.terms {
write!(fmt, "{} + ", term)?;
}
write!(fmt, "{}", self.offset)
}
}
#[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
}
}
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);
}
}
#[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
}
}
#[derive(Debug, Clone)]
#[allow(non_snake_case)]
pub struct StandardFormLP {
vars: Vec<Variable>,
c: DVector<Fl>,
A: DMatrix<Fl>,
b: DVector<Fl>,
}
impl StandardFormLP {
pub fn from_ineqs_and_objective(
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 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(nvars, ncstr);
let mut b = DVector::zeros(ncstr);
for (i, var) in vars.iter().enumerate() {
c[i] = objective.get_coeff(var);
for (j, constraint) in constraints.iter().enumerate() {
A[(i, j)] = constraint.get_coeff(var);
}
}
for (j, constr) in constraints.into_iter().enumerate() {
b[j] = constr.get_rhs();
}
StandardFormLP { vars, 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(mut self) -> StandardFormLP {
let x0 = DVector::from_element(self.vars.len(), 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);
c[lastcol] = -1.0; // Try to minimize u
let u = Variable::from_name("u".into());
self.vars.push(u);
Self { A, c, ..self }
}
}
pub enum IterationResult {
Distance(Fl),
Indeterminate(),
Unbounded(),
}
#[allow(non_snake_case)]
pub struct AffineScaling {
problem: StandardFormLP,
ε: Fl,
β: Fl,
pub x: DVector<Fl>, // current guess
// Here be local variables that we don't want to reallocate
AT: DMatrix<Fl>,
ADk2: DMatrix<Fl>,
lm: DMatrix<Fl>,
w: DVector<Fl>,
r: DVector<Fl>,
}
impl AffineScaling {
pub fn from_feasible_point_εβ(
problem: StandardFormLP,
point: DVector<Fl>,
ε: Fl,
β: Fl,
) -> Option<Self> {
if point.len() != problem.vars.len() {
return None;
}
let shape = problem.A.shape();
#[allow(non_snake_case)]
let AT = problem.A.transpose();
Some(AffineScaling {
problem,
x: point,
ε,
β,
AT,
ADk2: DMatrix::zeros(shape.0, shape.1),
lm: DMatrix::zeros(shape.0, shape.0),
w: DVector::zeros(shape.0),
r: DVector::zeros(shape.1),
})
}
pub fn from_feasible_point(problem: StandardFormLP, point: DVector<Fl>) -> Option<Self> {
Self::from_feasible_point_εβ(problem, point, DEFAULT_ε, DEFAULT_β)
}
// Returns an upper bound on the error
pub fn iterate(&mut self) -> IterationResult {
use IterationResult::*;
self.ADk2.copy_from(&self.problem.A);
for (i, mut col) in self.ADk2.column_iter_mut().enumerate() {
col *= self.x[i] * self.x[i];
}
self.ADk2.mul_to(&self.problem.c, &mut self.w);
self.ADk2.mul_to(&self.AT, &mut self.lm);
self.lm.clone().lu().solve_mut(&mut self.w);
self.problem.A.tr_mul_to(&self.w, &mut self.r);
self.r -= &self.problem.c;
let min = self.r[self.r.imin()];
self.r.component_mul_assign(&self.x);
let rv = {
if min >= 0.0 {
Distance(self.r.sum())
} else if self.r[self.r.imax()] < 0.0 {
Unbounded()
} else {
Indeterminate()
}
};
self.r.component_mul_assign(&self.x);
self.r *= self.β / self.r.norm();
self.x -= &self.r;
rv
}
pub fn iterate_until_precise(&mut self) -> Option<Fl> {
use IterationResult::*;
loop {
match self.iterate() {
Distance(dist) if dist < self.ε => break Some(dist),
Distance(_) => (),
Unbounded() => break None,
Indeterminate() => (),
}
}
}
}
#[cfg(test)]
mod test {
#[test]
fn inequality_test() {}
}

6
src/main.rs

@ -6,12 +6,12 @@
#[macro_use]
/// Useful for disjunctive clauses.
pub mod clause;
/// Variables -- both boolean, and otherwise.
pub mod variable;
/// Affine scaling, linear systems and inequalities, etc.
pub mod affine_scaling;
/// I made this before I knew I didn't need it. :/
pub mod hash_wrapper;
mod messy_minsat;
pub mod messy_minsat;
/// Variables -- both boolean, and otherwise.
pub mod variable;
fn main() {}

9
src/variable.rs

@ -1,5 +1,6 @@
// vim: ts=4 sw=4 et
use crate::affine_scaling;
use crate::clause::Literal;
use colored::{control::ShouldColorize, Colorize};
use core::fmt::{Display, Formatter};
@ -148,6 +149,14 @@ impl<T: core::fmt::Debug + Clone> core::fmt::Debug for Variable<T> {
}
}
impl core::ops::Mul<affine_scaling::Fl> for Variable<affine_scaling::Fl> {
type Output = affine_scaling::Term;
fn mul(self, coeff: affine_scaling::Fl) -> Self::Output {
affine_scaling::Term::from_var(self, coeff)
}
}
//impl core::fmt::Debug for Variable<bool> {
// fn fmt(&self, fmt: &mut Formatter) -> core::fmt::Result {
// <Self as core::fmt::Display>::fmt(self, fmt)

Loading…
Cancel
Save