Browse Source

Karmarkar's algorithm go brrrrr

karmarkar
Thomas Johnson 3 years ago
parent
commit
bcf0661921
  1. 1
      Cargo.toml
  2. 394
      src/karmarkar.rs
  3. 5
      src/main.rs
  4. 16
      src/variable.rs

1
Cargo.toml

@ -9,4 +9,3 @@ edition = "2018"
[dependencies]
colored = "1.9.3"
nalgebra = "0.21.0"
ndarray = "0.13.1"

394
src/affine_scaling.rs → src/karmarkar.rs

@ -6,14 +6,13 @@ use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Shl};
use std::collections::{HashMap, HashSet};
use nalgebra::{DMatrix, DVector};
use ndarray::{Array1, Array2};
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;
pub static DEFAULT_α: Fl = 0.25;
#[allow(non_upper_case_globals)]
pub static DEFAULT_ε: Fl = 1e-6;
@ -450,9 +449,9 @@ impl Display for Equality {
#[allow(non_snake_case)]
pub struct StandardFormLP {
pub vars: Vec<Variable>,
c: DVector<Fl>,
A: DMatrix<Fl>,
b: DVector<Fl>,
pub c: DVector<Fl>,
pub A: DMatrix<Fl>,
pub b: DVector<Fl>,
}
impl StandardFormLP {
@ -508,37 +507,6 @@ impl StandardFormLP {
}
}
pub fn ineqs_and_objective_to_simplex_tableau(
ineqs: impl IntoIterator<Item = LeqInequality>,
objective: LinearCombination,
) -> (Array2<Fl>, Array1<Fl>, Array1<Fl>) {
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 = Array1::zeros(nvars);
#[allow(non_snake_case)]
let mut A = Array2::zeros((ncstr, nvars));
let mut b = Array1::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);
}
(A, b, c)
}
impl Display for StandardFormLP {
fn fmt(&self, fmt: &mut Formatter) -> core::fmt::Result {
write!(fmt, "variables: ")?;
@ -556,149 +524,194 @@ impl Display for StandardFormLP {
}
}
#[derive(Debug, Clone, Copy)]
pub enum IterationResult {
Distance(Fl),
Indeterminate(),
Unbounded(),
#[derive(Debug, Clone)]
#[allow(non_snake_case)]
pub struct HomogenousLP {
pub vars: Vec<Variable>,
pub c: DVector<Fl>,
pub A: DMatrix<Fl>,
}
impl HomogenousLP {
pub fn project_point_inplace(point: &DVector<Fl>, out: &mut DVector<Fl>) {
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(point: &DVector<Fl>, out: &mut DVector<Fl>) {
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(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 {
vars: Vec::new(),
A,
c,
}
}
}
impl Display for HomogenousLP {
fn fmt(&self, fmt: &mut Formatter) -> core::fmt::Result {
write!(fmt, "variables: ")?;
for var in &self.vars {
write!(fmt, "{} ", var.get_name())?;
}
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 AffineScaling {
problem: StandardFormLP,
pub struct Karmarkar {
pub problem: HomogenousLP,
pub guess: DVector<Fl>,
ε: 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>,
nvars: usize,
nconstr: usize,
αr: 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,
β: Fl,
) -> Option<Self> {
if point.len() != problem.vars.len() {
) -> Option<Karmarkar> {
if value.len() != problem.c.len() {
return None;
}
let shape = problem.A.shape();
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 AT = problem.A.transpose();
Some(AffineScaling {
let Ashape = problem.A.shape();
Some(Karmarkar {
problem,
x: point,
guess: value,
ε,
β,
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),
nvars: Ashape.1,
nconstr: Ashape.0,
αr: α * r,
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_point(problem: StandardFormLP, point: DVector<Fl>) -> Option<Self> {
Self::from_feasible_point_εβ(problem, point, DEFAULT_ε, DEFAULT_β)
pub fn from_feasible_value(problem: HomogenousLP, value: DVector<Fl>) -> Option<Karmarkar> {
Self::from_feasible_value_αε(problem, value, DEFAULT_α, DEFAULT_ε)
}
// Returns an upper bound on the error
pub fn iterate(&mut self) -> IterationResult {
use IterationResult::*;
#[cfg(debug_affine_scaling)]
#[cfg(debug_karmarkar)]
{
print!("current infeasiblity:");
print!("{}", &self.problem.b - &self.problem.A * &self.x);
}
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];
// print!("infeasiblity: {}", &self.problem.A * &self.guess);
}
#[cfg(debug_affine_scaling)]
// First transform c, giving Dc
self.Dc.copy_from(&self.problem.c);
self.Dc.component_mul_assign(&self.guess);
// Then project Dc into the null space of B (A augmented with the simplex constraint)
{
print!("ADk2:");
print!("{}", &self.ADk2);
}
self.ADk2.mul_to(&self.problem.c, &mut self.w);
self.ADk2.mul_to(&self.AT, &mut self.lm);
#[cfg(debug_affine_scaling)]
let old_rhs = self.w.clone();
self.lm.clone().lu().solve_mut(&mut self.w);
#[cfg(debug_affine_scaling)]
{
print!("error in w:");
print!("{}", &self.lm * &self.w - old_rhs);
print!("w:");
print!("{}", &self.w);
#[allow(non_snake_case)]
let mut Bslice = self.B.rows_mut(0, self.nconstr);
Bslice.copy_from(&self.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.problem.A.tr_mul_to(&self.w, &mut self.r);
self.r -= &self.problem.c;
#[cfg(debug_affine_scaling)]
{
print!("r:");
print!("{}", &self.r);
self.B.transpose_to(&mut self.BT);
self.B.mul_to(&self.BT, &mut self.BBT);
if !self.BBT.try_inverse_mut() {
return Failed();
}
let min = self.r[self.r.imin()];
self.r.component_mul_assign(&self.x); // Reuse r no less than twice.
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();
#[cfg(debug_affine_scaling)]
{
print!("step:");
print!("{}", -&self.r);
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.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() => (),
}
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()) {
self.guess.copy_from(&self.cp);
Indeterminate()
} else {
Failed()
}
}
}
impl Display for AffineScaling {
fn fmt(&self, fmt: &mut Formatter) -> core::fmt::Result {
write!(fmt, "{}", &self.problem)?;
writeln!(fmt, "ε = {}", self.ε)?;
writeln!(fmt, "β = {}", self.β)?;
writeln!(fmt, "x:")?;
write!(fmt, "{}", &self.x)?;
Ok(())
}
#[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
@ -707,10 +720,10 @@ impl Display for AffineScaling {
macro_rules! linear_vars {
[] => {};
[$varname:ident$(,)?] => {
let $varname = $crate::affine_scaling::Variable::from_name(stringify!($varname).into());
let $varname = $crate::karmarkar::Variable::from_name(stringify!($varname).into());
};
[$varname:ident, $($varnames:ident),+] => {
let $varname = $crate::affine_scaling::Variable::from_name(stringify!($varname).into());
let $varname = $crate::karmarkar::Variable::from_name(stringify!($varname).into());
linear_vars![$($varnames),+]
};
}
@ -774,28 +787,77 @@ mod test {
println!("{}", constraint);
}
let tableau =
ineqs_and_objective_to_simplex_tableau(constraints.clone(), objective.clone());
println!("tableau: A = ");
println!("{}", tableau.0);
println!("b = {}", tableau.1);
println!("c = {}", tableau.2);
let sf = StandardFormLP::from_ineqs_and_objective(constraints, objective);
println!("");
println!("{}", sf);
let frp = sf.clone().to_feasible_region_problem();
println!("");
println!("{}", frp);
let len = frp.vars.len();
let mut fras =
AffineScaling::from_feasible_point(frp, DVector::from_element(len, 1.0)).unwrap();
let frp_feasible_point = DVector::from_element(frp.vars.len(), 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!("");
println!("{}", fras);
for _ in 0..0 {
println!("{:?}", fras.iterate());
println!("{}", fras.x);
let mut km =
Karmarkar::from_feasible_value(proj_frp.clone(), proj_frp_feasible_point).unwrap();
loop {
let oldx = km.guess.clone();
if let IterationResult::Failed() = km.iterate() {
break;
}
// print!("x: {}", km.guess);
// println!("distance: {}", (&km.guess - oldx).norm());
}
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();
for _ in 0..1000 {
let oldx = km.guess.clone();
if let IterationResult::Failed() = km.iterate() {
break;
}
print!("x: {}", km.guess);
println!("distance: {}", (&km.guess - oldx).norm());
}
print!("final x': {}", km.guess);
let solution = HomogenousLP::project_point_inv(&km.guess);
println!("vars: {:?}", sf.vars);
print!("final x: {}", solution);
}
}

5
src/main.rs

@ -4,7 +4,6 @@
#![feature(non_ascii_idents)]
extern crate nalgebra;
extern crate ndarray;
#[macro_use]
/// Variables -- both boolean, and otherwise.
@ -12,10 +11,10 @@ pub mod variable;
#[macro_use]
/// Useful for disjunctive clauses.
pub mod clause;
/// 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;
/// Affine scaling, linear systems and inequalities, etc.
pub mod karmarkar;
pub mod messy_minsat;
fn main() {}

16
src/variable.rs

@ -1,7 +1,7 @@
// vim: ts=4 sw=4 et
use crate::affine_scaling;
use crate::clause::Literal;
use crate::karmarkar;
use colored::{control::ShouldColorize, Colorize};
use core::fmt::{Display, Formatter};
use core::hash::{Hash, Hasher};
@ -149,19 +149,19 @@ 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;
impl core::ops::Mul<karmarkar::Fl> for &Variable<karmarkar::Fl> {
type Output = karmarkar::Term;
fn mul(self, coeff: affine_scaling::Fl) -> Self::Output {
affine_scaling::Term::from_var(self.clone(), coeff)
fn mul(self, coeff: karmarkar::Fl) -> Self::Output {
karmarkar::Term::from_var(self.clone(), coeff)
}
}
impl core::ops::Neg for &Variable<affine_scaling::Fl> {
type Output = affine_scaling::Term;
impl core::ops::Neg for &Variable<karmarkar::Fl> {
type Output = karmarkar::Term;
fn neg(self) -> Self::Output {
affine_scaling::Term::from_var(self.clone(), -1.0)
karmarkar::Term::from_var(self.clone(), -1.0)
}
}

Loading…
Cancel
Save