Browse Source

Made some small optimizations

karmarkar
Thomas Johnson 3 years ago
parent
commit
3155ea50e9
  1. 127
      src/karmarkar.rs

127
src/karmarkar.rs

@ -448,7 +448,6 @@ impl Display for Equality {
#[derive(Debug, Clone)]
#[allow(non_snake_case)]
pub struct StandardFormLP {
pub vars: Vec<Variable>,
pub c: DVector<Fl>,
pub A: DMatrix<Fl>,
pub b: DVector<Fl>,
@ -458,7 +457,7 @@ impl StandardFormLP {
pub fn from_ineqs_and_objective(
ineqs: impl IntoIterator<Item = LeqInequality>,
objective: LinearCombination,
) -> StandardFormLP {
) -> (Vec<Variable>, StandardFormLP) {
let mut counter = NewVariableCounter::new();
let constraints = ineqs
.into_iter()
@ -483,7 +482,36 @@ impl StandardFormLP {
for (j, var) in vars.iter().enumerate() {
c[j] = objective.get_coeff(var);
}
StandardFormLP { vars, c, A, b }
(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
@ -492,8 +520,8 @@ impl StandardFormLP {
/// 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);
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)]
@ -501,18 +529,25 @@ impl StandardFormLP {
A.set_column(lastcol, &v);
let mut c = DVector::zeros(lastcol + 1);
c[lastcol] = -1.0; // Try to minimize u
let u = Variable::from_name("u".into());
self.vars.push(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: ")?;
for var in &self.vars {
write!(fmt, "{} ", var.get_name())?;
}
writeln!(fmt, "")?;
writeln!(fmt, "A:")?;
write!(fmt, "{}", &self.A)?;
@ -527,7 +562,6 @@ impl Display for StandardFormLP {
#[derive(Debug, Clone)]
#[allow(non_snake_case)]
pub struct HomogenousLP {
pub vars: Vec<Variable>,
pub c: DVector<Fl>,
pub A: DMatrix<Fl>,
}
@ -572,26 +606,25 @@ impl HomogenousLP {
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_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,
}
HomogenousLP { 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)?;
@ -603,8 +636,8 @@ impl Display for HomogenousLP {
#[derive(Debug, Clone)]
#[allow(non_snake_case)]
pub struct Karmarkar {
pub problem: HomogenousLP,
pub struct Karmarkar<'a> {
pub problem: Option<&'a HomogenousLP>,
pub guess: DVector<Fl>,
ε: Fl,
nvars: usize,
@ -621,9 +654,9 @@ pub struct Karmarkar {
cp: DVector<Fl>,
}
impl Karmarkar {
impl<'a> Karmarkar<'a> {
pub fn from_feasible_value_αε(
problem: HomogenousLP,
problem: &'a HomogenousLP,
value: DVector<Fl>,
α: Fl,
ε: Fl,
@ -637,7 +670,7 @@ impl Karmarkar {
#[allow(non_snake_case)]
let Ashape = problem.A.shape();
Some(Karmarkar {
problem,
problem: Some(problem),
guess: value,
ε,
nvars: Ashape.1,
@ -655,26 +688,33 @@ impl Karmarkar {
})
}
pub fn from_feasible_value(problem: HomogenousLP, value: DVector<Fl>) -> Option<Karmarkar> {
pub fn from_feasible_value(
problem: &'a HomogenousLP,
value: DVector<Fl>,
) -> Option<Karmarkar<'a>> {
Self::from_feasible_value_αε(problem, value, DEFAULT_α, DEFAULT_ε)
}
pub fn iterate(&mut self) -> IterationResult {
use IterationResult::*;
let problem = match self.problem {
Some(x) => x,
None => return Failed(),
};
#[cfg(debug_karmarkar)]
{
print!("infeasiblity: {}", &self.problem.A * &self.guess);
print!("infeasiblity: {}", &problem.A * &self.guess);
}
// First transform c, giving Dc
self.Dc.copy_from(&self.problem.c);
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(&self.problem.A);
Bslice.copy_from(&problem.A);
for (i, mut col) in Bslice.column_iter_mut().enumerate() {
col.copy_from(&(&col * self.guess[i]));
}
@ -778,13 +818,13 @@ mod test {
println!("{}", constraint);
}
let sf = StandardFormLP::from_ineqs_and_objective(constraints, objective);
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.vars.len(), 1.0);
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!("");
@ -800,16 +840,7 @@ mod test {
);
println!("");
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 mut km = Karmarkar::from_feasible_value(&proj_frp, 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);
@ -838,23 +869,11 @@ mod test {
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();
//let mut olddist = 0.0;
// for _ in 0..100 {
// let oldx = km.guess.clone();
// if let IterationResult::Failed() = km.iterate() {
// break;
// }
// let dist = (&km.guess - oldx).norm();
// print!("x: {}", km.guess);
// println!("distance: {}", dist);
// println!("distance ratio: {}", dist / olddist);
// olddist = dist;
// }
let mut km = Karmarkar::from_feasible_value(&proj, 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: {:?}", sf.vars);
println!("vars: {:?}", vars);
print!("final x: {}", solution);
}
}
Loading…
Cancel
Save