From 469cfb62f8ebf8d9c316dac5be52fe4035776d6b Mon Sep 17 00:00:00 2001 From: Stefan Kroboth Date: Fri, 11 Feb 2022 18:01:05 +0100 Subject: [PATCH] Different approach, no need for ndarray --- argmin/Cargo.toml | 3 - argmin/examples/simplex.rs | 2 +- argmin/src/solver/simplex/mod.rs | 346 +++++++++++++++++++++---------- 3 files changed, 242 insertions(+), 109 deletions(-) diff --git a/argmin/Cargo.toml b/argmin/Cargo.toml index 5bac6a82a..e0f98d7e5 100644 --- a/argmin/Cargo.toml +++ b/argmin/Cargo.toml @@ -33,9 +33,6 @@ slog = { version = "2.4.1", optional = true } slog-term = { version = "2.8.1", optional = true } slog-async = { version = "2.7.0", optional = true } slog-json = { version = "2.5.0", optional = true } -# remove -ndarray = { version = "0.15", features = ["serde-1"] } -ndarray-linalg = { version = "0.14", features = ["netlib"] } [dev-dependencies] approx = "0.5.0" diff --git a/argmin/examples/simplex.rs b/argmin/examples/simplex.rs index 887434fd1..1a9e969a0 100644 --- a/argmin/examples/simplex.rs +++ b/argmin/examples/simplex.rs @@ -70,7 +70,7 @@ fn run() -> Result<(), Error> { let init_param = vec![]; let res = Executor::new(problem, solver, init_param) // .add_observer(ArgminSlogLogger::term(), ObserverMode::Always) - .max_iters(10) + .max_iters(3) .run()?; // Wait a second (lets the logger flush everything before printing again) diff --git a/argmin/src/solver/simplex/mod.rs b/argmin/src/solver/simplex/mod.rs index 42f5c35df..1e3c7f0aa 100644 --- a/argmin/src/solver/simplex/mod.rs +++ b/argmin/src/solver/simplex/mod.rs @@ -16,7 +16,6 @@ // IterState, OpWrapper, SerializeAlias, Solver, TerminationReason, // }; use crate::core::*; -use ndarray::prelude::*; #[cfg(feature = "serde1")] use serde::{Deserialize, Serialize}; @@ -36,17 +35,13 @@ where Simplex { tableau: None } } - pub(super) fn tableau_ref(&self) -> &Array2 { - self.tableau.as_ref().unwrap().tableau_ref() - } - pub(super) fn cost(&self) -> F { self.tableau.as_ref().unwrap().cost() } - pub(super) fn variables(&self) -> Vec { - self.tableau.as_ref().unwrap().variables() - } + // pub(super) fn variables(&self) -> Vec { + // self.tableau.as_ref().unwrap().variables() + // } pub(super) fn step(&mut self) { self.tableau.as_mut().unwrap().step() @@ -62,92 +57,131 @@ where } } -#[derive(Clone)] +#[derive(Clone, Debug)] #[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))] -struct Tableau { +struct Constraints { c: Vec, + coeff: Vec, b: Vec, - a: Vec>, - tableau: Array2, - basis_indices: Vec, + cost: F, num_constraints: usize, num_variables: usize, + constraints_counter: usize, } -impl Tableau +impl Constraints where F: ArgminFloat, { - #[allow(non_snake_case)] - pub fn new(c: Vec, b: Vec, A: Vec>) -> Tableau { - let num_constraints = b.len(); - let num_variables = c.len(); - let mut arr = Array2::zeros((num_constraints + 1, num_variables + 1 + num_constraints)); - arr.slice_mut(s![0, 1..=num_variables]) - .assign(&Array::from_iter(c.iter().cloned())); - arr.slice_mut(s![1..=num_constraints, 0]) - .assign(&Array::from_iter(b.iter().cloned())); - for (i, row) in A.iter().cloned().enumerate() { - arr.slice_mut(s![i + 1, 1..=num_variables]) - .assign(&Array::from_iter(row.into_iter())); + pub fn new(num_constraints: usize, num_variables: usize) -> Constraints { + Constraints { + c: vec![F::from_f64(0.0).unwrap(); num_variables + num_constraints], + coeff: vec![ + F::from_f64(0.0).unwrap(); + (num_constraints + num_variables) * num_constraints + ], + b: Vec::with_capacity(num_constraints), + cost: F::from_f64(0.0).unwrap(), + num_constraints, + num_variables, + constraints_counter: 0, } - // slack variables - for i in 1..=num_constraints { - arr.slice_mut(s![i, num_variables + i]) - .fill(F::from_f64(1.0).unwrap()); + } + + pub fn add_variable_coefficients(&mut self, vars: Vec) -> Result<(), Error> { + if vars.len() > self.num_variables { + return Err(ArgminError::InvalidParameter { + text: "Too many variables.".to_string(), + } + .into()); } - println!("ARR: {:?}", arr); + self.c[..self.num_variables].clone_from_slice(&vars); - Tableau { - c, - b, - a: A.to_vec(), - tableau: arr, - basis_indices: ((num_variables + 1)..=(num_variables + num_constraints)).collect(), - num_constraints, - num_variables, + Ok(()) + } + + pub fn add_constraint(&mut self, constraint: &[F], b: F) -> Result<(), Error> { + if self.constraints_counter >= self.num_constraints { + return Err(ArgminError::InvalidParameter { + text: "Too many constraints.".to_string(), + } + .into()); + } + for (i, v) in constraint.iter().cloned().enumerate() { + self.coeff + [self.constraints_counter * (self.num_constraints + self.num_variables) + i] = v; } + self.coeff[self.constraints_counter * (self.num_constraints + self.num_variables) + + self.num_variables + + self.constraints_counter] = F::from_f64(1.0).unwrap(); + self.b.push(b); + self.constraints_counter += 1; + Ok(()) } - pub fn tableau_ref(&self) -> &Array2 { - &self.tableau + pub fn get_column(&self, idx: usize) -> Vec { + let mut out = Vec::with_capacity(self.num_constraints); + for i in 0..self.num_constraints { + out.push(self.coeff[idx + i * (self.num_constraints + self.num_variables)]); + } + out } - pub fn cost(&self) -> F { - -self.tableau[[0, 0]] + pub fn get_element(&self, row: usize, column: usize) -> F { + self.coeff[row * (self.num_constraints + self.num_variables) + column] } - pub fn variables(&self) -> Vec { - self.basis_indices - .iter() - .map(|&i| { - if i >= self.num_constraints { - F::zero() - } else { - self.tableau[[i + 1, 0]] - } - }) - .collect() + pub fn normalize_row(&mut self, row: usize, factor: F) { + for i in 0..(self.num_constraints + self.num_variables) { + self.coeff[row * (self.num_constraints + self.num_variables) + i] = + self.coeff[row * (self.num_constraints + self.num_variables) + i] / factor; + } + self.b[row] = self.b[row] / factor; } - pub(super) fn step(&mut self) { - // Obtain next variable to enter basis - let (i_enter, _): (usize, F) = self.tableau.slice(s![0, 1..]).iter().enumerate().fold( - (0, F::infinity()), - |(iacc, acc), (i, &x)| if x < acc { (i, x) } else { (iacc, acc) }, - ); - println!("next basis idx = {:?}", i_enter); + pub fn transform_non_pivot_rows(&mut self, pivot_row: usize, pivot_column: usize) { + let row: Vec = ((pivot_row * (self.num_constraints + self.num_variables)) + ..((pivot_row + 1) * (self.num_constraints + self.num_variables))) + .into_iter() + .map(|i| self.coeff[i]) + .collect(); + let b_f = self.b[pivot_row]; + for i in 0..self.num_constraints { + if i == pivot_row { + continue; + } + let factor = self.get_element(i, pivot_column); + for (j, &elem) in row + .iter() + .enumerate() + .take(self.num_variables + self.num_constraints) + { + self.coeff[(i * (self.num_constraints + self.num_variables)) + j] = self.coeff + [(i * (self.num_constraints + self.num_variables)) + j] + - factor * elem; + } + self.b[i] = self.b[i] - factor * b_f; + } + let factor = self.c[pivot_column]; + for (j, &elem) in row + .iter() + .enumerate() + .take(self.num_variables + self.num_constraints) + { + self.c[j] = self.c[j] - factor * elem; + } + self.cost = self.cost - factor * b_f; + } - // Ratio test + fn get_pivot_row_idx(&self, column: usize) -> usize { let (i_leaves, _): (usize, F) = self - .tableau - .slice(s![1.., i_enter + 1]) + .get_column(column) .iter() - .zip(self.tableau.slice(s![1.., 0]).iter()) - .filter(|(&pc, _)| pc > F::from_f64(0.0).unwrap()) - .map(|(&pc, &bc): (&F, &F)| bc / pc) + .zip(self.b.iter()) .enumerate() + .filter(|(_, (&pc, _))| pc > F::from_f64(0.0).unwrap()) + .map(|(i, (&pc, &bc)): (usize, (&F, &F))| (i, bc / pc)) .fold((0, F::infinity()), |(iacc, acc), (i, x)| { if x < acc { (i, x) @@ -155,32 +189,135 @@ where (iacc, acc) } }); - println!("leaving basis idx = {:?}", i_leaves); - - let pivot_element = self.tableau[[i_leaves + 1, i_enter + 1]]; - println!("pivot element = {:?}", pivot_element); - - // normalize row of leaving variable - self.tableau - .slice_mut(s![i_leaves + 1, ..]) - .mapv_inplace(|x| x / pivot_element); - println!("{:?}", self.tableau); - - let pivot_row = self.tableau.slice(s![i_leaves + 1, ..]).to_owned(); - - // do thing with other rows - for i in 0..=self.num_constraints { - println!("{:?}, {:?}", i, self.tableau[[i, i_enter + 1]]); - if i != i_leaves + 1 { - let m = &self.tableau.slice(s![i, ..]) - - pivot_row.mapv(|x| x * self.tableau[[i, i_enter + 1]]); - self.tableau.slice_mut(s![i, ..]).assign(&m); - } + i_leaves + } + + fn get_pivot_column_idx(&self) -> usize { + let (i_enter, _): (usize, F) = + self.c + .iter() + .enumerate() + .fold((0, F::infinity()), |(iacc, acc), (i, &x)| { + if x < acc { + (i, x) + } else { + (iacc, acc) + } + }); + i_enter + } + + fn get_cost(&self) -> F { + -self.cost + } +} + +#[derive(Clone)] +#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))] +struct Tableau { + constraints: Constraints, + basis_indices: Vec, + num_constraints: usize, + num_variables: usize, + cost: F, +} + +impl Tableau +where + F: ArgminFloat, +{ + #[allow(non_snake_case)] + pub fn new(c: Vec, b: Vec, A: Vec>) -> Result, Error> { + let num_constraints = b.len(); + let num_variables = c.len(); + let mut constraints = Constraints::new(num_constraints, num_variables); + for (row, b_i) in A.iter().zip(b.into_iter()) { + constraints.add_constraint(row, b_i)?; } - println!("{:?}", self.tableau); + constraints.add_variable_coefficients(c)?; + println!("{:?}", constraints); + + Ok(Tableau { + constraints, + basis_indices: ((num_variables + 1)..=(num_variables + num_constraints)).collect(), + num_constraints, + num_variables, + cost: F::neg_infinity(), + }) + } + + pub fn cost(&self) -> F { + self.constraints.get_cost() + } + + // pub fn variables(&self) -> Vec { + // // let inds: Vec = self + // // .basis_indices + // // .into_iter() + // // .filter(|i| i < &self.num_variables) + // // .collect(); + // // let relevant_rows = self.tableau.slice(s![1.., inds]); + // // unimplemented!() + // self.basis_indices + // .iter() + // .filter(|&i| i < &self.num_variables) + // // .map(|&i| { + // // self.tableau + // // .slice(s![1.., i]) + // // .iter() + // // .enumerate() + // // .filter(|(_, &v)| v >= F::from_f64(0.0).unwrap()) + // // .map(|(j, _)| j) + // // .collect() + // // }) + // // .map(|j| self.tableau.slice(s![j + 1, 0])) + // .map(|&i| self.tableau[[i + 1, 0]]) + // // // .map(|&i| { + // // // if i >= self.num_variables { + // // // F::zero() + // // // } else { + // // // self.tableau[[i + 1, 0]] + // // // } + // // // }) + // .collect() + // } - self.basis_indices[i_leaves] = i_enter; - println!("new basis: {:?}", self.basis_indices); + fn pivot_column(&self) -> usize { + self.constraints.get_pivot_column_idx() + } + + fn pivot_row(&self, column: usize) -> usize { + self.constraints.get_pivot_row_idx(column) + } + + fn get_pivot_element(&self, row: usize, column: usize) -> F { + self.constraints.get_element(row, column) + } + + fn normalize_pivot_row(&mut self, row: usize, pivot_element: F) { + self.constraints.normalize_row(row, pivot_element); + } + + pub fn transform_non_pivot_rows(&mut self, pivot_row: usize, pivot_column: usize) { + self.constraints + .transform_non_pivot_rows(pivot_row, pivot_column) + } + + pub(super) fn step(&mut self) { + // Obtain next variable to enter basis + let pivot_column = self.pivot_column(); + println!("pivot_column = {:?}", pivot_column); + let pivot_row = self.pivot_row(pivot_column); + println!("pivot_row = {:?}", pivot_row); + let pivot_element = self.get_pivot_element(pivot_row, pivot_column); + println!("pivot_element = {:?}", pivot_element); + self.normalize_pivot_row(pivot_row, pivot_element); + println!("constraints = {:?}", self.constraints); + println!("b = {:?}", self.constraints.b); + println!("pivoting"); + self.transform_non_pivot_rows(pivot_row, pivot_column); + println!("constraints = {:?}", self.constraints); + println!("b = {:?}", self.constraints.b); } } @@ -196,10 +333,10 @@ where op: &mut OpWrapper, _state: &LinearProgramState, ) -> Result>>, Error> { - self.tableau = Some(Tableau::new(op.c()?, op.b()?, op.A()?)); + self.tableau = Some(Tableau::new(op.c()?, op.b()?, op.A()?)?); Ok(Some( ArgminIterData::new() - .param(self.variables()) + // .param(self.variables()) .cost(self.cost()), )) } @@ -211,24 +348,23 @@ where _state: &LinearProgramState, ) -> Result>, Error> { self.step(); - println!("{:?}", self.variables()); Ok(ArgminIterData::new() - .param(self.variables()) + // .param(self.variables()) .cost(self.cost())) } fn terminate(&mut self, _state: &LinearProgramState) -> TerminationReason { - if self - .tableau_ref() - .slice(s![0, 1..]) - .iter() - .cloned() - .map(|x| x.is_sign_positive()) - .fold(true, |acc, x| acc & x) - { - return TerminationReason::NoImprovementPossible; - } + // if self + // .tableau_ref() + // .slice(s![0, 1..]) + // .iter() + // .cloned() + // .map(|x| x.is_sign_positive()) + // .fold(true, |acc, x| acc & x) + // { + // return TerminationReason::NoImprovementPossible; + // } TerminationReason::NotTerminated } }