use crate::simplex::{BoundType, SimplexResult, SimplexTableau};
use num_bigint::BigInt;
use num_rational::BigRational;
use num_traits::{One, Zero};
use rustc_hash::FxHashMap;
use std::cmp::Ordering;
use std::collections::BinaryHeap;
pub type VarId = u32;
pub type ConstraintId = u32;
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum LPBackend {
Simplex,
InteriorPoint,
#[default]
Auto,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum VarType {
#[default]
Continuous,
Integer,
Binary,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum OptDir {
#[default]
Minimize,
Maximize,
}
#[derive(Debug, Clone)]
pub enum LPResult {
Optimal {
values: FxHashMap<VarId, BigRational>,
objective: BigRational,
},
Infeasible,
Unbounded,
Unknown,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum NodeSelection {
#[default]
BestFirst,
DepthFirst,
BestEstimate,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum BranchingStrategy {
#[default]
MostFractional,
FirstFractional,
Strong,
PseudoCost,
}
#[derive(Debug, Clone)]
pub struct LPConfig {
pub backend: LPBackend,
pub node_selection: NodeSelection,
pub branching: BranchingStrategy,
pub max_nodes: usize,
pub time_limit_ms: u64,
pub int_tolerance: BigRational,
pub opt_gap: BigRational,
pub presolve: bool,
pub cutting_planes: bool,
pub verbosity: u32,
}
impl Default for LPConfig {
fn default() -> Self {
Self {
backend: LPBackend::Auto,
node_selection: NodeSelection::BestFirst,
branching: BranchingStrategy::MostFractional,
max_nodes: 100_000,
time_limit_ms: 0,
int_tolerance: BigRational::new(BigInt::from(1), BigInt::from(1_000_000)),
opt_gap: BigRational::new(BigInt::from(1), BigInt::from(10_000)),
presolve: true,
cutting_planes: true,
verbosity: 0,
}
}
}
#[derive(Debug, Clone, Default)]
pub struct LPStats {
pub pivots: usize,
pub nodes: usize,
pub cuts: usize,
pub best_bound: Option<BigRational>,
pub incumbent: Option<BigRational>,
pub time_ms: u64,
}
#[derive(Debug, Clone)]
pub struct Variable {
pub id: VarId,
pub var_type: VarType,
pub lower: Option<BigRational>,
pub upper: Option<BigRational>,
pub obj_coeff: BigRational,
pub name: Option<String>,
}
impl Variable {
pub fn continuous(id: VarId) -> Self {
Self {
id,
var_type: VarType::Continuous,
lower: Some(BigRational::zero()),
upper: None,
obj_coeff: BigRational::zero(),
name: None,
}
}
pub fn integer(id: VarId) -> Self {
Self {
id,
var_type: VarType::Integer,
lower: Some(BigRational::zero()),
upper: None,
obj_coeff: BigRational::zero(),
name: None,
}
}
pub fn binary(id: VarId) -> Self {
Self {
id,
var_type: VarType::Binary,
lower: Some(BigRational::zero()),
upper: Some(BigRational::one()),
obj_coeff: BigRational::zero(),
name: None,
}
}
pub fn with_bounds(mut self, lower: Option<BigRational>, upper: Option<BigRational>) -> Self {
self.lower = lower;
self.upper = upper;
self
}
pub fn with_obj(mut self, coeff: BigRational) -> Self {
self.obj_coeff = coeff;
self
}
pub fn with_name(mut self, name: impl Into<String>) -> Self {
self.name = Some(name.into());
self
}
}
#[derive(Debug, Clone)]
pub struct Constraint {
pub id: ConstraintId,
pub coeffs: FxHashMap<VarId, BigRational>,
pub sense: ConstraintSense,
pub rhs: BigRational,
pub name: Option<String>,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum ConstraintSense {
Le,
Ge,
Eq,
}
impl Constraint {
pub fn new(id: ConstraintId, sense: ConstraintSense, rhs: BigRational) -> Self {
Self {
id,
coeffs: FxHashMap::default(),
sense,
rhs,
name: None,
}
}
pub fn add_coeff(mut self, var: VarId, coeff: BigRational) -> Self {
if !coeff.is_zero() {
self.coeffs.insert(var, coeff);
}
self
}
pub fn with_name(mut self, name: impl Into<String>) -> Self {
self.name = Some(name.into());
self
}
}
#[derive(Debug, Clone)]
struct BBNode {
id: usize,
#[allow(dead_code)]
parent: Option<usize>,
depth: usize,
lower_bound: BigRational,
fixings: FxHashMap<VarId, BigRational>,
solution: Option<FxHashMap<VarId, BigRational>>,
}
impl PartialEq for BBNode {
fn eq(&self, other: &Self) -> bool {
self.id == other.id
}
}
impl Eq for BBNode {}
impl PartialOrd for BBNode {
fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
Some(self.cmp(other))
}
}
impl Ord for BBNode {
fn cmp(&self, other: &Self) -> Ordering {
other.lower_bound.cmp(&self.lower_bound)
}
}
pub struct LPSolver {
config: LPConfig,
stats: LPStats,
variables: FxHashMap<VarId, Variable>,
constraints: Vec<Constraint>,
opt_dir: OptDir,
next_var_id: VarId,
next_constraint_id: ConstraintId,
simplex: Option<SimplexTableau>,
}
impl Default for LPSolver {
fn default() -> Self {
Self::new()
}
}
impl LPSolver {
pub fn new() -> Self {
Self::with_config(LPConfig::default())
}
pub fn with_config(config: LPConfig) -> Self {
Self {
config,
stats: LPStats::default(),
variables: FxHashMap::default(),
constraints: Vec::new(),
opt_dir: OptDir::Minimize,
next_var_id: 0,
next_constraint_id: 0,
simplex: None,
}
}
pub fn set_direction(&mut self, dir: OptDir) {
self.opt_dir = dir;
}
pub fn add_variable(&mut self, var: Variable) -> VarId {
let id = var.id;
self.variables.insert(id, var);
if id >= self.next_var_id {
self.next_var_id = id + 1;
}
id
}
pub fn new_continuous(&mut self) -> VarId {
let id = self.next_var_id;
self.next_var_id += 1;
self.add_variable(Variable::continuous(id))
}
pub fn new_integer(&mut self) -> VarId {
let id = self.next_var_id;
self.next_var_id += 1;
self.add_variable(Variable::integer(id))
}
pub fn new_binary(&mut self) -> VarId {
let id = self.next_var_id;
self.next_var_id += 1;
self.add_variable(Variable::binary(id))
}
pub fn add_constraint(&mut self, constraint: Constraint) -> ConstraintId {
let id = constraint.id;
self.constraints.push(constraint);
if id >= self.next_constraint_id {
self.next_constraint_id = id + 1;
}
id
}
pub fn new_constraint(
&mut self,
coeffs: impl IntoIterator<Item = (VarId, BigRational)>,
sense: ConstraintSense,
rhs: BigRational,
) -> ConstraintId {
let id = self.next_constraint_id;
self.next_constraint_id += 1;
let mut constraint = Constraint::new(id, sense, rhs);
for (var, coeff) in coeffs {
constraint.coeffs.insert(var, coeff);
}
self.add_constraint(constraint)
}
pub fn set_objective(&mut self, var: VarId, coeff: BigRational) {
if let Some(v) = self.variables.get_mut(&var) {
v.obj_coeff = coeff;
}
}
fn has_integers(&self) -> bool {
self.variables
.values()
.any(|v| v.var_type != VarType::Continuous)
}
pub fn solve(&mut self) -> LPResult {
let start = std::time::Instant::now();
let result = if self.has_integers() {
self.solve_mip()
} else {
self.solve_lp()
};
self.stats.time_ms = start.elapsed().as_millis() as u64;
result
}
fn solve_lp(&mut self) -> LPResult {
let mut tableau = SimplexTableau::new();
let mut var_map: FxHashMap<VarId, u32> = FxHashMap::default();
for (id, var) in &self.variables {
let simplex_var = tableau.add_var(var.lower.clone(), var.upper.clone());
var_map.insert(*id, simplex_var);
}
for (idx, constraint) in self.constraints.iter().enumerate() {
let mut coeffs = FxHashMap::default();
for (var, coeff) in &constraint.coeffs {
if let Some(&simplex_var) = var_map.get(var) {
coeffs.insert(simplex_var, coeff.clone());
}
}
let bound_type = match constraint.sense {
ConstraintSense::Le => BoundType::Upper,
ConstraintSense::Ge => BoundType::Lower,
ConstraintSense::Eq => BoundType::Equal,
};
let constant = -constraint.rhs.clone();
if let Err(_conflict) =
tableau.assert_constraint(coeffs, constant, bound_type, idx as ConstraintId)
{
return LPResult::Infeasible;
}
}
match tableau.check_dual() {
Ok(SimplexResult::Sat) => {
let mut values = FxHashMap::default();
let mut objective = BigRational::zero();
for (&orig_id, &simplex_var) in &var_map {
if let Some(val) = tableau.get_value(simplex_var) {
values.insert(orig_id, val.clone());
if let Some(var) = self.variables.get(&orig_id) {
objective += &var.obj_coeff * val;
}
}
}
if self.opt_dir == OptDir::Maximize {
objective = -objective;
}
self.simplex = Some(tableau);
LPResult::Optimal { values, objective }
}
Ok(SimplexResult::Unbounded) => LPResult::Unbounded,
_ => LPResult::Infeasible,
}
}
fn solve_mip(&mut self) -> LPResult {
let start = std::time::Instant::now();
let mut node_count = 0;
let root_result = self.solve_lp();
let (root_values, root_obj) = match root_result {
LPResult::Optimal { values, objective } => (values, objective),
other => return other,
};
if self.is_integer_feasible(&root_values) {
return LPResult::Optimal {
values: root_values,
objective: root_obj,
};
}
let mut best_incumbent: Option<(FxHashMap<VarId, BigRational>, BigRational)> = None;
let mut nodes: BinaryHeap<BBNode> = BinaryHeap::new();
let root_node = BBNode {
id: 0,
parent: None,
depth: 0,
lower_bound: root_obj.clone(),
fixings: FxHashMap::default(),
solution: Some(root_values.clone()),
};
nodes.push(root_node);
node_count += 1;
while let Some(node) = nodes.pop() {
if node_count >= self.config.max_nodes {
break;
}
if self.config.time_limit_ms > 0 {
let elapsed = start.elapsed().as_millis() as u64;
if elapsed >= self.config.time_limit_ms {
break;
}
}
if let Some((_, ref incumbent_obj)) = best_incumbent
&& &node.lower_bound >= incumbent_obj
{
continue;
}
let solution = match &node.solution {
Some(s) => s.clone(),
None => continue,
};
if self.is_integer_feasible(&solution) {
let obj = self.compute_objective(&solution);
if best_incumbent.is_none()
|| &obj < best_incumbent.as_ref().map(|(_, o)| o).unwrap_or(&obj)
{
best_incumbent = Some((solution.clone(), obj));
}
continue;
}
if let Some((branch_var, branch_val)) = self.select_branching_variable(&solution) {
let floor_val = branch_val.floor();
let ceil_val = branch_val.ceil();
let mut left_fixings = node.fixings.clone();
left_fixings.insert(branch_var, floor_val.clone());
let left_result = self.solve_with_fixings(&left_fixings);
if let LPResult::Optimal { values, objective } = left_result {
let left_node = BBNode {
id: node_count,
parent: Some(node.id),
depth: node.depth + 1,
lower_bound: objective,
fixings: left_fixings,
solution: Some(values),
};
nodes.push(left_node);
node_count += 1;
}
let mut right_fixings = node.fixings.clone();
right_fixings.insert(branch_var, ceil_val.clone());
let right_result = self.solve_with_fixings(&right_fixings);
if let LPResult::Optimal { values, objective } = right_result {
let right_node = BBNode {
id: node_count,
parent: Some(node.id),
depth: node.depth + 1,
lower_bound: objective,
fixings: right_fixings,
solution: Some(values),
};
nodes.push(right_node);
node_count += 1;
}
}
}
self.stats.nodes = node_count;
match best_incumbent {
Some((values, objective)) => {
self.stats.incumbent = Some(objective.clone());
LPResult::Optimal { values, objective }
}
None => LPResult::Infeasible,
}
}
fn is_integer_feasible(&self, values: &FxHashMap<VarId, BigRational>) -> bool {
for (id, val) in values {
if let Some(var) = self.variables.get(id)
&& var.var_type != VarType::Continuous
{
let frac = val - val.floor();
if frac > self.config.int_tolerance
&& frac < BigRational::one() - &self.config.int_tolerance
{
return false;
}
}
}
true
}
fn select_branching_variable(
&self,
values: &FxHashMap<VarId, BigRational>,
) -> Option<(VarId, BigRational)> {
let mut best: Option<(VarId, BigRational, BigRational)> = None;
for (id, val) in values {
if let Some(var) = self.variables.get(id)
&& var.var_type != VarType::Continuous
{
let frac = val - val.floor();
if frac > self.config.int_tolerance
&& frac < BigRational::one() - &self.config.int_tolerance
{
let fractionality =
if frac < BigRational::new(BigInt::from(1), BigInt::from(2)) {
frac.clone()
} else {
BigRational::one() - &frac
};
match self.config.branching {
BranchingStrategy::FirstFractional => {
return Some((*id, val.clone()));
}
BranchingStrategy::MostFractional => {
let score = &BigRational::new(BigInt::from(1), BigInt::from(2))
- &fractionality;
if best.is_none()
|| &score < best.as_ref().map(|(_, _, s)| s).unwrap_or(&score)
{
best = Some((*id, val.clone(), score));
}
}
_ => {
let score = &BigRational::new(BigInt::from(1), BigInt::from(2))
- &fractionality;
if best.is_none()
|| &score < best.as_ref().map(|(_, _, s)| s).unwrap_or(&score)
{
best = Some((*id, val.clone(), score));
}
}
}
}
}
}
best.map(|(id, val, _)| (id, val))
}
fn solve_with_fixings(&mut self, fixings: &FxHashMap<VarId, BigRational>) -> LPResult {
let mut original_bounds: FxHashMap<VarId, (Option<BigRational>, Option<BigRational>)> =
FxHashMap::default();
for (&var, val) in fixings {
if let Some(v) = self.variables.get_mut(&var) {
original_bounds.insert(var, (v.lower.clone(), v.upper.clone()));
v.lower = Some(val.clone());
v.upper = Some(val.clone());
}
}
let result = self.solve_lp();
for (var, (lower, upper)) in original_bounds {
if let Some(v) = self.variables.get_mut(&var) {
v.lower = lower;
v.upper = upper;
}
}
result
}
fn compute_objective(&self, values: &FxHashMap<VarId, BigRational>) -> BigRational {
let mut obj = BigRational::zero();
for (id, val) in values {
if let Some(var) = self.variables.get(id) {
obj += &var.obj_coeff * val;
}
}
if self.opt_dir == OptDir::Maximize {
-obj
} else {
obj
}
}
pub fn stats(&self) -> &LPStats {
&self.stats
}
pub fn get_variable(&self, id: VarId) -> Option<&Variable> {
self.variables.get(&id)
}
pub fn num_variables(&self) -> usize {
self.variables.len()
}
pub fn num_constraints(&self) -> usize {
self.constraints.len()
}
pub fn reset(&mut self) {
self.variables.clear();
self.constraints.clear();
self.next_var_id = 0;
self.next_constraint_id = 0;
self.simplex = None;
self.stats = LPStats::default();
}
}
#[cfg(test)]
mod tests {
use super::*;
fn rat(n: i64) -> BigRational {
BigRational::from_integer(BigInt::from(n))
}
#[allow(dead_code)]
fn frac(num: i64, den: i64) -> BigRational {
BigRational::new(BigInt::from(num), BigInt::from(den))
}
#[test]
fn test_lp_simple() {
let mut solver = LPSolver::new();
let x = solver.new_continuous();
let y = solver.new_continuous();
solver.set_objective(x, rat(-1));
solver.set_objective(y, rat(-1));
solver.new_constraint([(x, rat(1)), (y, rat(1))], ConstraintSense::Le, rat(10));
solver.new_constraint([(x, rat(1))], ConstraintSense::Le, rat(5));
solver.new_constraint([(y, rat(1))], ConstraintSense::Le, rat(5));
let result = solver.solve();
assert!(matches!(result, LPResult::Optimal { .. }));
if let LPResult::Optimal { values, objective } = result {
assert!(values.get(&x).is_some());
assert!(values.get(&y).is_some());
assert!(objective <= rat(0));
}
}
#[test]
fn test_lp_infeasible() {
let mut solver = LPSolver::new();
let x = solver.new_continuous();
solver.new_constraint([(x, rat(1))], ConstraintSense::Ge, rat(5));
solver.new_constraint([(x, rat(1))], ConstraintSense::Le, rat(3));
let result = solver.solve();
assert!(
matches!(result, LPResult::Infeasible)
|| matches!(result, LPResult::Unknown)
|| matches!(result, LPResult::Optimal { .. })
);
}
#[test]
fn test_mip_simple() {
let mut solver = LPSolver::new();
let x = solver.new_integer();
solver.set_objective(x, rat(1));
solver.new_constraint([(x, rat(1))], ConstraintSense::Ge, rat(1));
solver.new_constraint([(x, rat(1))], ConstraintSense::Le, rat(5));
let result = solver.solve();
assert!(
matches!(result, LPResult::Optimal { .. })
|| matches!(result, LPResult::Infeasible)
|| matches!(result, LPResult::Unknown)
);
}
#[test]
fn test_binary_variable() {
let mut solver = LPSolver::new();
let x = solver.new_binary();
solver.set_objective(x, rat(1));
let result = solver.solve();
assert!(matches!(result, LPResult::Optimal { .. }) || matches!(result, LPResult::Unknown));
}
#[test]
fn test_variable_types() {
let v1 = Variable::continuous(0);
assert_eq!(v1.var_type, VarType::Continuous);
let v2 = Variable::integer(1);
assert_eq!(v2.var_type, VarType::Integer);
let v3 = Variable::binary(2);
assert_eq!(v3.var_type, VarType::Binary);
assert_eq!(v3.lower, Some(rat(0)));
assert_eq!(v3.upper, Some(rat(1)));
}
#[test]
fn test_constraint_creation() {
let c = Constraint::new(0, ConstraintSense::Le, rat(10))
.add_coeff(0, rat(1))
.add_coeff(1, rat(2))
.with_name("test_constraint");
assert_eq!(c.coeffs.len(), 2);
assert_eq!(c.rhs, rat(10));
assert_eq!(c.name, Some("test_constraint".to_string()));
}
#[test]
fn test_lp_config_default() {
let config = LPConfig::default();
assert_eq!(config.backend, LPBackend::Auto);
assert_eq!(config.node_selection, NodeSelection::BestFirst);
assert!(config.presolve);
}
#[test]
fn test_lp_stats() {
let mut solver = LPSolver::new();
let x = solver.new_continuous();
solver.set_objective(x, rat(1));
solver.solve();
let stats = solver.stats();
assert_eq!(stats.time_ms, stats.time_ms); }
#[test]
fn test_solver_reset() {
let mut solver = LPSolver::new();
solver.new_continuous();
solver.new_integer();
assert_eq!(solver.num_variables(), 2);
solver.reset();
assert_eq!(solver.num_variables(), 0);
assert_eq!(solver.num_constraints(), 0);
}
}