#![deny(missing_debug_implementations, missing_docs)]
#[macro_use]
extern crate log;
mod helpers;
mod lu;
mod mps;
mod ordering;
mod solver;
mod sparse;
use solver::Solver;
#[derive(Clone, Copy, Debug)]
pub enum OptimizationDirection {
Minimize,
Maximize,
}
#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub struct Variable(pub(crate) usize);
impl Variable {
pub fn idx(&self) -> usize {
self.0
}
}
#[derive(Clone, Debug)]
pub struct LinearExpr {
vars: Vec<usize>,
coeffs: Vec<f64>,
}
impl LinearExpr {
pub fn empty() -> Self {
Self {
vars: vec![],
coeffs: vec![],
}
}
pub fn add(&mut self, var: Variable, coeff: f64) {
self.vars.push(var.0);
self.coeffs.push(coeff);
}
}
#[doc(hidden)]
#[derive(Clone, Copy, Debug)]
pub struct LinearTerm(Variable, f64);
impl From<(Variable, f64)> for LinearTerm {
fn from(term: (Variable, f64)) -> Self {
LinearTerm(term.0, term.1)
}
}
impl<'a> From<&'a (Variable, f64)> for LinearTerm {
fn from(term: &'a (Variable, f64)) -> Self {
LinearTerm(term.0, term.1)
}
}
impl<I: IntoIterator<Item = impl Into<LinearTerm>>> From<I> for LinearExpr {
fn from(iter: I) -> Self {
let mut expr = LinearExpr::empty();
for term in iter {
let LinearTerm(var, coeff) = term.into();
expr.add(var, coeff);
}
expr
}
}
impl std::iter::FromIterator<(Variable, f64)> for LinearExpr {
fn from_iter<I: IntoIterator<Item = (Variable, f64)>>(iter: I) -> Self {
let mut expr = LinearExpr::empty();
for term in iter {
expr.add(term.0, term.1)
}
expr
}
}
impl std::iter::Extend<(Variable, f64)> for LinearExpr {
fn extend<I: IntoIterator<Item = (Variable, f64)>>(&mut self, iter: I) {
for term in iter {
self.add(term.0, term.1)
}
}
}
#[derive(Clone, Copy, Debug)]
pub enum ComparisonOp {
Eq,
Le,
Ge,
}
#[derive(Clone, Debug, PartialEq)]
pub enum Error {
Infeasible,
Unbounded,
}
impl std::fmt::Display for Error {
fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
let msg = match self {
Error::Infeasible => "problem is infeasible",
Error::Unbounded => "problem is unbounded",
};
msg.fmt(f)
}
}
impl std::error::Error for Error {}
#[derive(Clone)]
pub struct Problem {
direction: OptimizationDirection,
obj_coeffs: Vec<f64>,
var_mins: Vec<f64>,
var_maxs: Vec<f64>,
constraints: Vec<(CsVec, ComparisonOp, f64)>,
}
impl std::fmt::Debug for Problem {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
f.debug_struct("Problem")
.field("direction", &self.direction)
.field("num_vars", &self.obj_coeffs.len())
.field("num_constraints", &self.constraints.len())
.finish()
}
}
type CsVec = sprs::CsVecI<f64, usize>;
impl Problem {
pub fn new(direction: OptimizationDirection) -> Self {
Problem {
direction,
obj_coeffs: vec![],
var_mins: vec![],
var_maxs: vec![],
constraints: vec![],
}
}
pub fn add_var(&mut self, obj_coeff: f64, (min, max): (f64, f64)) -> Variable {
let var = Variable(self.obj_coeffs.len());
let obj_coeff = match self.direction {
OptimizationDirection::Minimize => obj_coeff,
OptimizationDirection::Maximize => -obj_coeff,
};
self.obj_coeffs.push(obj_coeff);
self.var_mins.push(min);
self.var_maxs.push(max);
var
}
pub fn add_constraint(&mut self, expr: impl Into<LinearExpr>, cmp_op: ComparisonOp, rhs: f64) {
let expr = expr.into();
self.constraints.push((
CsVec::new(self.obj_coeffs.len(), expr.vars, expr.coeffs),
cmp_op,
rhs,
));
}
pub fn solve(&self) -> Result<Solution, Error> {
let mut solver = Solver::try_new(
&self.obj_coeffs,
&self.var_mins,
&self.var_maxs,
&self.constraints,
)?;
solver.initial_solve()?;
Ok(Solution {
num_vars: self.obj_coeffs.len(),
direction: self.direction,
solver,
})
}
}
#[derive(Clone)]
pub struct Solution {
direction: OptimizationDirection,
num_vars: usize,
solver: solver::Solver,
}
impl std::fmt::Debug for Solution {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
f.debug_struct("Solution")
.field("direction", &self.direction)
.field("num_vars", &self.num_vars)
.field("num_constraints", &self.solver.num_constraints())
.field("objective", &self.objective())
.finish()
}
}
impl Solution {
pub fn objective(&self) -> f64 {
match self.direction {
OptimizationDirection::Minimize => self.solver.cur_obj_val,
OptimizationDirection::Maximize => -self.solver.cur_obj_val,
}
}
pub fn var_value(&self, var: Variable) -> &f64 {
assert!(var.0 < self.num_vars);
self.solver.get_value(var.0)
}
pub fn iter(&self) -> SolutionIter {
SolutionIter {
solution: self,
var_idx: 0,
}
}
pub fn add_constraint(
mut self,
expr: impl Into<LinearExpr>,
cmp_op: ComparisonOp,
rhs: f64,
) -> Result<Self, Error> {
let expr = expr.into();
self.solver.add_constraint(
CsVec::new(self.num_vars, expr.vars, expr.coeffs),
cmp_op,
rhs,
)?;
Ok(self)
}
pub fn fix_var(mut self, var: Variable, val: f64) -> Result<Self, Error> {
assert!(var.0 < self.num_vars);
self.solver.fix_var(var.0, val)?;
Ok(self)
}
pub fn unfix_var(mut self, var: Variable) -> (Self, bool) {
assert!(var.0 < self.num_vars);
let res = self.solver.unfix_var(var.0);
(self, res)
}
pub fn add_gomory_cut(mut self, var: Variable) -> Result<Self, Error> {
assert!(var.0 < self.num_vars);
self.solver.add_gomory_cut(var.0)?;
Ok(self)
}
}
impl std::ops::Index<Variable> for Solution {
type Output = f64;
fn index(&self, var: Variable) -> &Self::Output {
self.var_value(var)
}
}
#[derive(Debug, Clone)]
pub struct SolutionIter<'a> {
solution: &'a Solution,
var_idx: usize,
}
impl<'a> Iterator for SolutionIter<'a> {
type Item = (Variable, &'a f64);
fn next(&mut self) -> Option<Self::Item> {
if self.var_idx < self.solution.num_vars {
let var_idx = self.var_idx;
self.var_idx += 1;
Some((Variable(var_idx), self.solution.solver.get_value(var_idx)))
} else {
None
}
}
}
impl<'a> IntoIterator for &'a Solution {
type Item = (Variable, &'a f64);
type IntoIter = SolutionIter<'a>;
fn into_iter(self) -> Self::IntoIter {
self.iter()
}
}
pub use mps::MpsFile;
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn optimize() {
let mut problem = Problem::new(OptimizationDirection::Maximize);
let v1 = problem.add_var(3.0, (12.0, f64::INFINITY));
let v2 = problem.add_var(4.0, (5.0, f64::INFINITY));
problem.add_constraint(&[(v1, 1.0), (v2, 1.0)], ComparisonOp::Le, 20.0);
problem.add_constraint(&[(v2, -4.0), (v1, 1.0)], ComparisonOp::Ge, -20.0);
let sol = problem.solve().unwrap();
assert_eq!(sol[v1], 12.0);
assert_eq!(sol[v2], 8.0);
assert_eq!(sol.objective(), 68.0);
}
#[test]
fn empty_expr_constraints() {
let trivial = [
(LinearExpr::empty(), ComparisonOp::Eq, 0.0),
(LinearExpr::empty(), ComparisonOp::Ge, -1.0),
(LinearExpr::empty(), ComparisonOp::Le, 1.0),
];
let mut problem = Problem::new(OptimizationDirection::Minimize);
let _ = problem.add_var(1.0, (0.0, f64::INFINITY));
for (expr, op, b) in trivial.iter().cloned() {
problem.add_constraint(expr, op, b);
}
assert_eq!(problem.solve().map(|s| s.objective()), Ok(0.0));
{
let mut sol = problem.solve().unwrap();
for (expr, op, b) in trivial.iter().cloned() {
sol = sol.add_constraint(expr, op, b).unwrap();
}
assert_eq!(sol.objective(), 0.0);
}
let infeasible = [
(LinearExpr::empty(), ComparisonOp::Eq, 12.0),
(LinearExpr::empty(), ComparisonOp::Ge, 34.0),
(LinearExpr::empty(), ComparisonOp::Le, -56.0),
];
for (expr, op, b) in infeasible.iter().cloned() {
let mut cloned = problem.clone();
cloned.add_constraint(expr, op, b);
assert_eq!(cloned.solve().map(|_| "solved"), Err(Error::Infeasible));
}
for (expr, op, b) in infeasible.iter().cloned() {
let sol = problem.solve().unwrap().add_constraint(expr, op, b);
assert_eq!(sol.map(|_| "solved"), Err(Error::Infeasible));
}
let _ = problem.add_var(-1.0, (0.0, f64::INFINITY));
assert_eq!(problem.solve().map(|_| "solved"), Err(Error::Unbounded));
}
#[test]
fn free_variables() {
let mut problem = Problem::new(OptimizationDirection::Maximize);
let v1 = problem.add_var(1.0, (0.0, f64::INFINITY));
let v2 = problem.add_var(2.0, (f64::NEG_INFINITY, f64::INFINITY));
problem.add_constraint(&[(v1, 1.0), (v2, 1.0)], ComparisonOp::Le, 4.0);
problem.add_constraint(&[(v1, 1.0), (v2, 1.0)], ComparisonOp::Ge, 2.0);
problem.add_constraint(&[(v1, 1.0), (v2, -1.0)], ComparisonOp::Ge, 0.0);
let sol = problem.solve().unwrap();
assert_eq!(sol[v1], 2.0);
assert_eq!(sol[v2], 2.0);
assert_eq!(sol.objective(), 6.0);
}
#[test]
fn fix_unfix_var() {
let mut problem = Problem::new(OptimizationDirection::Maximize);
let v1 = problem.add_var(1.0, (0.0, 3.0));
let v2 = problem.add_var(2.0, (0.0, 3.0));
problem.add_constraint(&[(v1, 1.0), (v2, 1.0)], ComparisonOp::Le, 4.0);
problem.add_constraint(&[(v1, 1.0), (v2, 1.0)], ComparisonOp::Ge, 1.0);
let orig_sol = problem.solve().unwrap();
{
let mut sol = orig_sol.clone().fix_var(v1, 0.5).unwrap();
assert_eq!(sol[v1], 0.5);
assert_eq!(sol[v2], 3.0);
assert_eq!(sol.objective(), 6.5);
sol = sol.unfix_var(v1).0;
assert_eq!(sol[v1], 1.0);
assert_eq!(sol[v2], 3.0);
assert_eq!(sol.objective(), 7.0);
}
{
let mut sol = orig_sol.clone().fix_var(v2, 2.5).unwrap();
assert_eq!(sol[v1], 1.5);
assert_eq!(sol[v2], 2.5);
assert_eq!(sol.objective(), 6.5);
sol = sol.unfix_var(v2).0;
assert_eq!(sol[v1], 1.0);
assert_eq!(sol[v2], 3.0);
assert_eq!(sol.objective(), 7.0);
}
}
#[test]
fn add_constraint() {
let mut problem = Problem::new(OptimizationDirection::Minimize);
let v1 = problem.add_var(2.0, (0.0, f64::INFINITY));
let v2 = problem.add_var(1.0, (0.0, f64::INFINITY));
problem.add_constraint(&[(v1, 1.0), (v2, 1.0)], ComparisonOp::Le, 4.0);
problem.add_constraint(&[(v1, 1.0), (v2, 1.0)], ComparisonOp::Ge, 2.0);
let orig_sol = problem.solve().unwrap();
{
let sol = orig_sol
.clone()
.add_constraint(&[(v1, -1.0), (v2, 1.0)], ComparisonOp::Le, 0.0)
.unwrap();
assert_eq!(sol[v1], 1.0);
assert_eq!(sol[v2], 1.0);
assert_eq!(sol.objective(), 3.0);
}
{
let sol = orig_sol
.clone()
.fix_var(v2, 1.5)
.unwrap()
.add_constraint(&[(v1, -1.0), (v2, 1.0)], ComparisonOp::Le, 0.0)
.unwrap();
assert_eq!(sol[v1], 1.5);
assert_eq!(sol[v2], 1.5);
assert_eq!(sol.objective(), 4.5);
}
{
let sol = orig_sol
.clone()
.add_constraint(&[(v1, -1.0), (v2, 1.0)], ComparisonOp::Ge, 3.0)
.unwrap();
assert_eq!(sol[v1], 0.0);
assert_eq!(sol[v2], 3.0);
assert_eq!(sol.objective(), 3.0);
}
}
#[test]
fn gomory_cut() {
let mut problem = Problem::new(OptimizationDirection::Minimize);
let v1 = problem.add_var(0.0, (0.0, f64::INFINITY));
let v2 = problem.add_var(-1.0, (0.0, f64::INFINITY));
problem.add_constraint(&[(v1, 3.0), (v2, 2.0)], ComparisonOp::Le, 6.0);
problem.add_constraint(&[(v1, -3.0), (v2, 2.0)], ComparisonOp::Le, 0.0);
let mut sol = problem.solve().unwrap();
assert_eq!(sol[v1], 1.0);
assert_eq!(sol[v2], 1.5);
assert_eq!(sol.objective(), -1.5);
sol = sol.add_gomory_cut(v2).unwrap();
assert!(f64::abs(sol[v1] - 2.0 / 3.0) < 1e-8);
assert_eq!(sol[v2], 1.0);
assert_eq!(sol.objective(), -1.0);
sol = sol.add_gomory_cut(v1).unwrap();
assert!(f64::abs(sol[v1] - 1.0) < 1e-8);
assert_eq!(sol[v2], 1.0);
assert_eq!(sol.objective(), -1.0);
}
}