From 419de61a52f8a96f8db694e750404903dc5e3f1b Mon Sep 17 00:00:00 2001 From: nrealus Date: Sun, 27 Sep 2026 12:41:56 +0200 Subject: [PATCH 1/2] feat: solve a model in place Model::[try_]solve_in_place solves without consuming the model, so that it can be modified and solved again with a warm start. Its results are read with Model::status, objective_value and get_solution, as on SolvedModel. To reuse a model, Model::[try_]clear_solver resets the solver state, Model::[try_]clear_model removes the problem, and Model::[try_]overwrite replaces it (sharing the problem passing code of Model::try_new). Model::num_nz and get_column_bounds / get_column_cost, and Problem::get_column_bounds / get_column_cost, read the current model. --- src/lib.rs | 201 +++++++++++++++++++++++++++++++++++++++++++++- tests/in_place.rs | 108 +++++++++++++++++++++++++ 2 files changed, 307 insertions(+), 2 deletions(-) create mode 100644 tests/in_place.rs diff --git a/src/lib.rs b/src/lib.rs index 6a2da8b..cd51bc6 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -113,7 +113,7 @@ use std::ffi::{c_void, CStr, CString}; use std::num::{NonZeroU32, TryFromIntError}; use std::ops::{Bound, Index, RangeBounds}; use std::os::raw::c_int; -use std::ptr::null; +use std::ptr::{null, null_mut}; use highs_sys::*; @@ -234,6 +234,16 @@ where self.collower[col.index()] = low; self.colupper[col.index()] = high; } + + /// The bounds `(lower, upper)` of a column + pub fn get_column_bounds(&self, col: Col) -> (f64, f64) { + (self.collower[col.index()], self.colupper[col.index()]) + } + + /// The objective coefficient of a column + pub fn get_column_cost(&self, col: Col) -> f64 { + self.colcost[col.index()] + } } fn bound_value + Copy>(b: Bound<&N>) -> Option { @@ -423,6 +433,12 @@ impl Model { self.highs.num_rows().expect("num rows does not fit usize") } + /// Gets the number of non-zero coefficients in the constraint matrix + pub fn num_nz(&self) -> usize { + let n = unsafe { Highs_getNumNz(self.as_ptr()) }; + n.try_into().expect("num nz does not fit usize") + } + /// Create a Highs model to be optimized (but don't solve it yet). /// If the given problem is a [RowProblem], it will have to be converted to a [ColProblem] first, /// which takes an amount of time proportional to the size of the problem. @@ -438,7 +454,81 @@ impl Model { pub fn try_new>>(problem: P) -> Result { let mut highs = HighsPtr::default(); highs.make_quiet(); + Self::pass_problem(&mut highs, &problem.into())?; + Ok(Self { highs }) + } + + /// Replaces the model with the given problem, but keeps the HiGHS instance and its options. + /// + /// As with [`Model::new`], the objective sense is reset to minimisation + /// (see [`Model::set_sense`]), and a [RowProblem] is first converted to a [ColProblem]. + /// + /// # Panics + /// + /// If the problem is incoherent. + pub fn overwrite>>(&mut self, problem: P) { + self.try_overwrite(problem).expect("incoherent problem") + } + + /// Tries to replace the model with the given problem, but keeps the HiGHS instance and its options. + /// + /// As with [`Model::try_new`], the objective sense is reset to minimisation + /// (see [`Model::set_sense`]), and a [RowProblem] is first converted to a [ColProblem]. + /// + /// Returns an error if the problem is incoherent. + pub fn try_overwrite>>( + &mut self, + problem: P, + ) -> Result<(), HighsStatus> { let problem = problem.into(); + self.try_clear_model()?; + Self::pass_problem(&mut self.highs, &problem)?; + Ok(()) + } + + /// Removes all variables and constraints, but keeps the HiGHS instance and its options. + /// + /// # Panics + /// + /// If HIGHS returns an error status value. + pub fn clear_model(&mut self) { + self.try_clear_model() + .unwrap_or_else(|e| panic!("HiGHS error: {e:?}")) + } + + /// Tries to remove all variables and constraints, but keeps the HiGHS instance and its options. + /// + /// Returns the error status value if HIGHS returned an error status. + pub fn try_clear_model(&mut self) -> Result<(), HighsStatus> { + unsafe { highs_call!(Highs_clearModel(self.highs.mut_ptr())) }?; + Ok(()) + } + + /// Resets the solver state (basis and solution) but keeps the problem data, + /// so that the next solve starts from scratch. + /// + /// # Panics + /// + /// If HIGHS returns an error status value. + pub fn clear_solver(&mut self) { + self.try_clear_solver() + .unwrap_or_else(|e| panic!("HiGHS error: {e:?}")) + } + + /// Tries to reset the solver state (basis and solution) but keeps the problem data, + /// so that the next solve starts from scratch. + /// + /// Returns the error status value if HIGHS returned an error status. + pub fn try_clear_solver(&mut self) -> Result<(), HighsStatus> { + unsafe { highs_call!(Highs_clearSolver(self.highs.mut_ptr())) }?; + Ok(()) + } + + /// Pass the problem data to a HiGHS instance with an empty model + fn pass_problem( + highs: &mut HighsPtr, + problem: &Problem, + ) -> Result { log::debug!( "Adding a problem with {} variables and {} constraints to HiGHS", problem.num_cols(), @@ -484,7 +574,6 @@ impl Model { problem.matrix.avalue.as_ptr() )) } - .map(|_| Self { highs }) } } @@ -556,6 +645,70 @@ impl Model { .map(|_| SolvedModel { highs: self.highs }) } + /// Find the optimal value for the problem, keeping the model alive. + /// + /// Unlike [`Model::solve`], this does not consume the model: it can be modified and solved + /// again, and HiGHS warm-starts from the previous basis. + /// + /// Returns the resulting model status. + /// + /// # Panics + /// + /// If HIGHS returns an error status value. + pub fn solve_in_place(&mut self) -> HighsModelStatus { + self.try_solve_in_place() + .unwrap_or_else(|e| panic!("HiGHS error: {e:?}")) + } + + /// Find the optimal value for the problem, keeping the model alive. + /// + /// Unlike [`Model::try_solve`], this does not consume the model: it can be modified and + /// solved again, and HiGHS warm-starts from the previous basis. + /// + /// Returns the resulting model status, or the error status value if HIGHS returned an error status. + pub fn try_solve_in_place(&mut self) -> Result { + unsafe { highs_call!(Highs_run(self.highs.mut_ptr())) }?; + Ok(self.status()) + } + + /// The model status after the last solve ([`HighsModelStatus::NotSet`] before any solve). + pub fn status(&self) -> HighsModelStatus { + let model_status = unsafe { Highs_getModelStatus(self.highs.unsafe_mut_ptr()) }; + HighsModelStatus::try_from(model_status).unwrap() + } + + /// The objective value after the last solve. + /// + /// If an error occurs (e.g. the model is infeasible) then the returned value may be zero. + pub fn objective_value(&self) -> f64 { + unsafe { Highs_getObjectiveValue(self.as_ptr()) } + } + + /// The solution found by the last solve. + pub fn get_solution(&self) -> Solution { + let cols = self.num_cols(); + let rows = self.num_rows(); + let mut colvalue: Vec = vec![0.; cols]; + let mut coldual: Vec = vec![0.; cols]; + let mut rowvalue: Vec = vec![0.; rows]; + let mut rowdual: Vec = vec![0.; rows]; + unsafe { + Highs_getSolution( + self.highs.unsafe_mut_ptr(), + colvalue.as_mut_ptr(), + coldual.as_mut_ptr(), + rowvalue.as_mut_ptr(), + rowdual.as_mut_ptr(), + ); + } + Solution { + colvalue, + coldual, + rowvalue, + rowdual, + } + } + /// Adds a new constraint to the highs model. /// /// Returns the added row index. @@ -832,6 +985,50 @@ impl Model { } } + /// The current bounds `(lower, upper)` of a column. + /// + /// # Panics + /// + /// If the column does not exist. + pub fn get_column_bounds(&self, col: Col) -> (f64, f64) { + let (_, lower, upper) = self.column_data(col); + (lower, upper) + } + + /// The current objective coefficient of a column. + /// + /// # Panics + /// + /// If the column does not exist. + pub fn get_column_cost(&self, col: Col) -> f64 { + self.column_data(col).0 + } + + /// `(cost, lower, upper)` of a column + fn column_data(&self, col: Col) -> (f64, f64, f64) { + let index = c(col.index()); + let mut num_col: HighsInt = 0; + let mut num_nz: HighsInt = 0; + let (mut cost, mut lower, mut upper) = (0., 0., 0.); + unsafe { + highs_call!(Highs_getColsByRange( + self.as_ptr(), + index, + index, + &mut num_col, + &mut cost, + &mut lower, + &mut upper, + &mut num_nz, + null_mut(), + null_mut(), + null_mut() + )) + } + .unwrap_or_else(|e| panic!("HiGHS error: {e:?}")); + (cost, lower, upper) + } + /// Hot-starts at the initial guess. See HIGHS documentation for further details. /// /// # Panics diff --git a/tests/in_place.rs b/tests/in_place.rs new file mode 100644 index 0000000..a99124b --- /dev/null +++ b/tests/in_place.rs @@ -0,0 +1,108 @@ +use highs::{ColProblem, HighsModelStatus, RowProblem, Sense}; + +#[test] +fn solve_in_place_then_resolve() { + let mut pb = RowProblem::new(); + let x = pb.add_column(1., 0..); + let y = pb.add_column(2., 0..); + pb.add_row(..=6., [(x, 3.), (y, 1.)]); + let mut model = pb.optimise(Sense::Maximise); + assert_eq!(model.status(), HighsModelStatus::NotSet); + + assert_eq!(model.solve_in_place(), HighsModelStatus::Optimal); + assert_eq!(model.get_solution().columns(), &[0., 6.]); + assert_eq!(model.objective_value(), 12.); + + model.change_column_bounds(x, 1..); + assert_eq!(model.solve_in_place(), HighsModelStatus::Optimal); + assert_eq!(model.get_solution().columns(), &[1., 3.]); + + model.change_column_cost(y, 1.); + model.solve_in_place(); + assert_eq!(model.objective_value(), 4.); +} + +#[test] +fn add_row_and_col_then_solve_in_place() { + let mut model = ColProblem::default().optimise(Sense::Minimise); + let col = model.add_col(1., 1.., vec![]); + model.add_row(..1., vec![(col, 1.)]); + assert_eq!(model.num_nz(), 1); + assert_eq!(model.solve_in_place(), HighsModelStatus::Optimal); + assert_eq!(model.get_solution().columns(), &[1.]); + + model.change_column_bounds(col, 2..); + assert_eq!(model.solve_in_place(), HighsModelStatus::Infeasible); +} + +#[test] +fn column_getters() { + let mut model = ColProblem::default().optimise(Sense::Minimise); + let col = model.add_col(3., 1..5, vec![]); + assert_eq!(model.get_column_bounds(col), (1., 5.)); + assert_eq!(model.get_column_cost(col), 3.); + model.change_column_bounds(col, 0..); + model.change_column_cost(col, -1.); + assert_eq!(model.get_column_bounds(col), (0., f64::INFINITY)); + assert_eq!(model.get_column_cost(col), -1.); +} + +#[test] +#[should_panic] +fn column_getter_out_of_range() { + let mut other = ColProblem::default().optimise(Sense::Minimise); + let col = other.add_col(1., 0..1, vec![]); + let model = ColProblem::default().optimise(Sense::Minimise); + model.get_column_cost(col); +} + +#[test] +fn problem_column_getters() { + let mut pb = RowProblem::new(); + let col = pb.add_column(3., 1..); + assert_eq!(pb.get_column_bounds(col), (1., f64::INFINITY)); + assert_eq!(pb.get_column_cost(col), 3.); + pb.change_column_bounds(col, 0..=4); + pb.change_column_cost(col, 2.); + assert_eq!(pb.get_column_bounds(col), (0., 4.)); + assert_eq!(pb.get_column_cost(col), 2.); +} + +#[test] +fn clear_solver_then_resolve() { + let mut pb = RowProblem::new(); + pb.add_column(1., 0..10); + let mut model = pb.optimise(Sense::Maximise); + assert_eq!(model.solve_in_place(), HighsModelStatus::Optimal); + model.clear_solver(); + assert_eq!(model.status(), HighsModelStatus::NotSet); + assert_eq!(model.solve_in_place(), HighsModelStatus::Optimal); + assert_eq!(model.get_solution().columns(), &[10.]); +} + +#[test] +fn clear_model() { + let mut pb = RowProblem::new(); + let x = pb.add_column(1., 0..10); + pb.add_row(..=5, [(x, 1.)]); + let mut model = pb.optimise(Sense::Maximise); + model.clear_model(); + assert_eq!((model.num_cols(), model.num_rows()), (0, 0)); +} + +#[test] +fn overwrite() { + let mut pb = RowProblem::new(); + pb.add_column(1., 0..10); + let mut model = pb.optimise(Sense::Maximise); + model.solve_in_place(); + assert_eq!(model.get_solution().columns(), &[10.]); + + let mut other = ColProblem::new(); + let row = other.add_row(..=2.5); + other.add_integer_column(1., 0.., [(row, 1.)]); + model.overwrite(other); + model.set_sense(Sense::Maximise); + assert_eq!(model.solve_in_place(), HighsModelStatus::Optimal); + assert_eq!(model.get_solution().columns(), &[2.]); +} From 8730582b39f774b34a7e0c5893a1613195356f1c Mon Sep 17 00:00:00 2001 From: nrealus Date: Sun, 27 Sep 2026 12:48:45 +0200 Subject: [PATCH 2/2] feat: compute infeasible subsystems (IIS) Add HighsIisBoundStatus, the Iis type, and [try_]get_iis on SolvedModel and Model (via Highs_getIis), as well as Model::[try_]solve_or_iis, which solves in place and returns either the solution or an IIS. An IIS is only returned when HiGHS could check that it is infeasible: otherwise try_get_iis returns Err(HighsStatus::Warning). The per-column and per-row IIS statuses are not requested, as HiGHS copies them without having set them in that case. Requires the IIS bound status constants from highs-sys. --- Cargo.toml | 2 +- src/iis.rs | 216 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 6 +- src/status.rs | 33 ++++++++ tests/iis.rs | 122 ++++++++++++++++++++++++++++ 5 files changed, 377 insertions(+), 2 deletions(-) create mode 100644 src/iis.rs create mode 100644 tests/iis.rs diff --git a/Cargo.toml b/Cargo.toml index 628b7fa..a32e46b 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -9,5 +9,5 @@ repository = "https://github.com/rust-or/highs" keywords = ["linear-programming", "optimization", "math", "solver"] [dependencies] -highs-sys = "1.14.3" +highs-sys = "1.15.0" log = "0.4.27" diff --git a/src/iis.rs b/src/iis.rs new file mode 100644 index 0000000..3623ffe --- /dev/null +++ b/src/iis.rs @@ -0,0 +1,216 @@ +use std::convert::{TryFrom, TryInto}; +use std::ptr::null_mut; + +use highs_sys::*; + +use crate::{ + try_handle_status, Col, HighsIisBoundStatus, HighsModelStatus, HighsPtr, HighsStatus, Model, + Row, Solution, SolvedModel, +}; + +/// An infeasible subsystem (IIS) of a model: +/// bounds of some of its columns (variables) and rows (constraints) that together are infeasible. +/// Depending on the [`iis_strategy`](https://github.com/ERGO-Code/HiGHS/blob/master/docs/src/guide/advanced.md#irreducible-infeasibility-system-iis-detectionid-highs-iis) option used, +/// the subsystem could also be reduced to a minimal (irreducible) one. +/// +/// Returned by [`SolvedModel::get_iis`] and [`Model::get_iis`](crate::Model::get_iis). +/// HiGHS checks that the subsystem is infeasible before returning it: when it cannot, +/// the `try_get_iis` methods return `Err(HighsStatus::Warning)`. +/// +/// # Search strategy +/// +/// How hard HiGHS searches depends on the `iis_strategy` option +/// (see [`Model::set_option`](crate::Model::set_option) and [`iis_strategy`](https://github.com/ERGO-Code/HiGHS/blob/master/docs/src/guide/advanced.md#irreducible-infeasibility-system-iis-detectionid-highs-iis)). +/// - `0` (the default) only finds trivial subsystems: inconsistent bounds, or a single row that +/// cannot be satisfied within the bounds of its columns. Otherwise, the IIS is empty. +/// - `2` searches the whole model, and may return a subsystem that is not irreducible (minimal). +/// - `6` (`2 + 4`) also reduces the subsystem to an irreducible one, at a higher cost. +/// +/// HiGHS also flags some columns and rows as "maybe in conflict": +/// the members of a subsystem not proven irreducible, and every column and row when it finds no subsystem. +/// This flag is not exposed, and ignoring it is not unsound: the listed bounds are infeasible on +/// their own, and an empty IIS says nothing about which columns and rows are involved. +/// +/// # Unsolved models +/// +/// If the model was not solved, or its solve stopped early (e.g. on a time limit), HiGHS first +/// solves it again, with the `iis_time_limit` option (no limit by default) as time limit. +/// This changes the status and solution of the model. +/// +/// # Integer, semi-continuous and semi-integer columns +/// +/// HiGHS checks the subsystem on the continuous relaxation of the model, ignoring integrality. +/// An infeasibility that only comes from integrality is thus never returned: +/// - with the default `iis_strategy` (`0`), the IIS is empty; +/// - with a search of the whole model (e.g. `iis_strategy` `2` or `6`), the search keeps +/// integrality, but the check fails, and the `try_get_iis` methods return `Err(HighsStatus::Warning)`. +/// +/// Semi-continuous and semi-integer columns are not relaxed: whatever the `iis_strategy`, HiGHS +/// keeps their `[lower, upper]` bounds, ignoring that they may also be 0. +/// An IIS involving such columns may thus not be infeasible in the original model. +#[derive(Clone, Debug)] +pub struct Iis { + columns: Vec<(Col, HighsIisBoundStatus)>, + rows: Vec<(Row, HighsIisBoundStatus)>, +} + +impl Iis { + /// Whether the IIS has no columns and no rows: HiGHS found no infeasible subsystem + pub fn is_empty(&self) -> bool { + self.columns.is_empty() && self.rows.is_empty() + } + + /// The columns (variables) in the IIS, with the bounds that take part in it + pub fn columns(&self) -> &[(Col, HighsIisBoundStatus)] { + &self.columns + } + + /// The rows (constraints) in the IIS, with the bounds that take part in it + pub fn rows(&self) -> &[(Row, HighsIisBoundStatus)] { + &self.rows + } + + /// Whether the given column is part of the IIS + pub fn contains_column(&self, col: Col) -> bool { + self.columns.iter().any(|&(c, _)| c == col) + } + + /// Whether the given row is part of the IIS + pub fn contains_row(&self, row: Row) -> bool { + self.rows.iter().any(|&(r, _)| r == row) + } +} + +impl HighsPtr { + /// Compute an IIS of the current model with `Highs_getIis` + pub(crate) fn try_get_iis(&mut self) -> Result { + let cols = self.num_cols()?; + let rows = self.num_rows()?; + let mut iis_num_col: HighsInt = 0; + let mut iis_num_row: HighsInt = 0; + let mut col_index: Vec = vec![0; cols]; + let mut row_index: Vec = vec![0; rows]; + let mut col_bound: Vec = vec![0; cols]; + let mut row_bound: Vec = vec![0; rows]; + + // The per-column and per-row statuses are not requested: HiGHS (1.15) leaves them + // unset when it cannot check the subsystem, but still copies them. + let status = unsafe { + highs_call!(Highs_getIis( + self.mut_ptr(), + &mut iis_num_col, + &mut iis_num_row, + col_index.as_mut_ptr(), + row_index.as_mut_ptr(), + col_bound.as_mut_ptr(), + row_bound.as_mut_ptr(), + null_mut(), + null_mut() + )) + }?; + // Notably returned when HiGHS cannot check that the subsystem is infeasible + if status == HighsStatus::Warning { + return Err(status); + } + + let iis_num_col: usize = iis_num_col.try_into()?; + let iis_num_row: usize = iis_num_row.try_into()?; + col_index.truncate(iis_num_col); + col_bound.truncate(iis_num_col); + row_index.truncate(iis_num_row); + row_bound.truncate(iis_num_row); + + let columns = col_index + .into_iter() + .zip(col_bound) + .map(|(i, b)| Ok((Col(i.try_into()?), bound_status(b)))) + .collect::>()?; + let rows = row_index + .into_iter() + .zip(row_bound) + .map(|(i, b)| (Row(i), bound_status(b))) + .collect(); + Ok(Iis { columns, rows }) + } +} + +fn bound_status(raw: HighsInt) -> HighsIisBoundStatus { + HighsIisBoundStatus::try_from(raw) + .unwrap_or_else(|e| panic!("HiGHS returned an unrecognized IIS bound status: {e:?}")) +} + +impl SolvedModel { + /// Compute an infeasible subsystem of the model: see [`Iis`]. + /// + /// Meant to be called when [`SolvedModel::status`] is + /// [`Infeasible`](crate::HighsModelStatus::Infeasible). + /// + /// # Panics + /// + /// If HIGHS returns an error status value, or cannot check that the subsystem is infeasible. + pub fn get_iis(&mut self) -> Iis { + self.try_get_iis() + .unwrap_or_else(|e| panic!("HiGHS error: {e:?}")) + } + + /// Tries to compute an infeasible subsystem of the model: see [`Iis`]. + /// + /// Returns `Err(HighsStatus::Warning)` if HiGHS emits a warning, notably when it cannot check + /// that the subsystem is infeasible, or the error status value if HIGHS returned an error status. + pub fn try_get_iis(&mut self) -> Result { + self.highs.try_get_iis() + } +} + +impl Model { + /// Compute an infeasible subsystem of the model: see [`Iis`]. + /// + /// Meant to be called when [`Model::status`] is [`HighsModelStatus::Infeasible`]. + /// + /// # Panics + /// + /// If HIGHS returns an error status value, or cannot check that the subsystem is infeasible. + pub fn get_iis(&mut self) -> Iis { + self.try_get_iis() + .unwrap_or_else(|e| panic!("HiGHS error: {e:?}")) + } + + /// Tries to compute an infeasible subsystem of the model: see [`Iis`]. + /// + /// Returns `Err(HighsStatus::Warning)` if HiGHS emits a warning, notably when it cannot check + /// that the subsystem is infeasible, or the error status value if HIGHS returned an error status. + pub fn try_get_iis(&mut self) -> Result { + self.highs.try_get_iis() + } + + /// Solve the model in place, and return an infeasible subsystem if it is infeasible + /// (see [`Iis`], which may be empty), or its solution otherwise. + /// + /// The solution may not be optimal, e.g. if the solve reached a time limit: + /// see [`Model::status`]. + /// + /// # Panics + /// + /// If HIGHS returns an error status value, or cannot check that the subsystem is infeasible. + pub fn solve_or_iis(&mut self) -> Result { + self.try_solve_or_iis() + .unwrap_or_else(|e| panic!("HiGHS error: {e:?}")) + } + + /// Solve the model in place, and return an infeasible subsystem if it is infeasible + /// (see [`Iis`], which may be empty), or its solution otherwise. + /// + /// The solution may not be optimal, e.g. if the solve reached a time limit: + /// see [`Model::status`]. + /// + /// Returns `Err(HighsStatus::Warning)` if HiGHS emits a warning while computing the + /// subsystem, notably when it cannot check that it is infeasible, + /// or the error status value if HIGHS returned an error status. + pub fn try_solve_or_iis(&mut self) -> Result, HighsStatus> { + if self.try_solve_in_place()? == HighsModelStatus::Infeasible { + Ok(Err(self.try_get_iis()?)) + } else { + Ok(Ok(self.get_solution())) + } + } +} diff --git a/src/lib.rs b/src/lib.rs index cd51bc6..bdaf47a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -117,10 +117,11 @@ use std::ptr::{null, null_mut}; use highs_sys::*; +pub use iis::Iis; pub use matrix_col::{ColMatrix, Row}; pub use matrix_row::{Col, RowMatrix}; pub use options::{HighsOptionValue, TrySetOptionError}; -pub use status::{HighsModelStatus, HighsSolutionStatus, HighsStatus}; +pub use status::{HighsIisBoundStatus, HighsModelStatus, HighsSolutionStatus, HighsStatus}; /// A problem where variables are declared first, and constraints are then added dynamically. /// See [`Problem`](Problem#impl-1). @@ -266,6 +267,9 @@ macro_rules! highs_call { } } +// Declared after `highs_call!` so that the macro is in scope in this module. +mod iis; + /// A model to solve #[derive(Debug)] pub struct Model { diff --git a/src/status.rs b/src/status.rs index 4e52739..2d9c925 100644 --- a/src/status.rs +++ b/src/status.rs @@ -152,3 +152,36 @@ impl TryFrom for HighsSolutionStatus { } } } + +/// The status of a bound of a column or row that is part of an IIS. +#[derive(Clone, Copy, Debug, PartialOrd, PartialEq, Ord, Eq)] +pub enum HighsIisBoundStatus { + /// The bound was dropped from the IIS + Dropped = IIS_BOUND_STATUS_DROPPED as isize, + /// No bound + Null = IIS_BOUND_STATUS_NULL as isize, + /// The column or row is free + Free = IIS_BOUND_STATUS_FREE as isize, + /// The lower bound is part of the IIS + Lower = IIS_BOUND_STATUS_LOWER as isize, + /// The upper bound is part of the IIS + Upper = IIS_BOUND_STATUS_UPPER as isize, + /// Both bounds are part of the IIS + Boxed = IIS_BOUND_STATUS_BOXED as isize, +} + +impl TryFrom for HighsIisBoundStatus { + type Error = InvalidStatus; + + fn try_from(value: c_int) -> Result { + match value { + IIS_BOUND_STATUS_DROPPED => Ok(Self::Dropped), + IIS_BOUND_STATUS_NULL => Ok(Self::Null), + IIS_BOUND_STATUS_FREE => Ok(Self::Free), + IIS_BOUND_STATUS_LOWER => Ok(Self::Lower), + IIS_BOUND_STATUS_UPPER => Ok(Self::Upper), + IIS_BOUND_STATUS_BOXED => Ok(Self::Boxed), + n => Err(InvalidStatus(n)), + } + } +} diff --git a/tests/iis.rs b/tests/iis.rs new file mode 100644 index 0000000..79ee4bc --- /dev/null +++ b/tests/iis.rs @@ -0,0 +1,122 @@ +use highs::{Col, Row, SolvedModel}; +use highs::{ColProblem, HighsIisBoundStatus, HighsModelStatus, HighsStatus, Model, Sense}; + +use HighsIisBoundStatus::{Free, Lower, Upper}; + +fn model_with_iis_strategy(iis_strategy: Option) -> Model { + let mut model = ColProblem::new().optimise(Sense::Minimise); + if let Some(strategy) = iis_strategy { + model.set_option("iis_strategy", strategy); + } + model +} + +fn solve_infeasible(model: Model) -> SolvedModel { + let solved = model.solve(); + assert_eq!(solved.status(), HighsModelStatus::Infeasible); + solved +} + +#[test] +fn solved_model_get_iis() { + let mut model = model_with_iis_strategy(None); + let x = model.add_col(0., 1.., []); // x >= 1 + let y = model.add_col(0., 1.., []); // y >= 1 + let z = model.add_col(0., 1.., []); // z >= 1 + let r0 = model.add_row(..=1.1, [(x, 1.), (z, 1.)]); // x + z <= 1.1: infeasible with x, z >= 1 + let r1 = model.add_row(5.., [(y, 1.)]); // y >= 5: independent + + let iis = solve_infeasible(model).get_iis(); + assert_eq!(iis.columns(), &[(x, Lower), (z, Lower)]); + assert_eq!(iis.rows(), &[(r0, Upper)]); + assert!(!iis.is_empty()); + assert!(iis.contains_column(x)); + assert!(!iis.contains_column(y)); + assert!(iis.contains_column(z)); + assert!(iis.contains_row(r0)); + assert!(!iis.contains_row(r1)); +} + +/// x + y >= 3, x <= 1, y <= 1: infeasible, but not because of a single row +fn non_trivially_infeasible(iis_strategy: Option) -> (SolvedModel, [Col; 2], [Row; 3]) { + let mut model = model_with_iis_strategy(iis_strategy); + let x = model.add_col(0., 0..10, []); + let y = model.add_col(0., 0..10, []); + let r0 = model.add_row(3.., [(x, 1.), (y, 1.)]); + let r1 = model.add_row(..=1, [(x, 1.)]); + let r2 = model.add_row(..=1, [(y, 1.)]); + (solve_infeasible(model), [x, y], [r0, r1, r2]) +} + +#[test] +fn default_iis_strategy_only_finds_trivial_subsystems() { + let (mut solved, _, _) = non_trivially_infeasible(None); + assert!(solved.get_iis().is_empty()); +} + +#[test] +fn irreducible_iis_strategy() { + let (mut solved, [x, y], [r0, r1, r2]) = non_trivially_infeasible(Some(2 + 4)); + let iis = solved.get_iis(); + assert_eq!(iis.columns(), &[(x, Free), (y, Free)]); + assert_eq!(iis.rows(), &[(r0, Lower), (r1, Upper), (r2, Upper)]); +} + +/// 2y = 1 with y integer: the relaxation (y = 0.5) is feasible +fn infeasible_because_of_integrality(iis_strategy: Option) -> SolvedModel { + let mut model = model_with_iis_strategy(iis_strategy); + let y = model.add_integer_column(0., 0..10, []); + model.add_row(1.0..=1.0, [(y, 2.)]); + solve_infeasible(model) +} + +#[test] +fn iis_of_mip_infeasible_because_of_integrality_is_empty() { + assert!(infeasible_because_of_integrality(None).get_iis().is_empty()); +} + +#[test] +fn iis_of_mip_infeasible_because_of_integrality_cannot_be_checked() { + let mut solved = infeasible_because_of_integrality(Some(2)); + assert_eq!(solved.try_get_iis().err(), Some(HighsStatus::Warning)); +} + +#[test] +fn iis_of_mip_with_infeasible_relaxation() { + // x + y <= 1 with x >= 1 and y >= 1, y integer. + let mut model = model_with_iis_strategy(None); + let x = model.add_col(0., 1.., []); + let y = model.add_integer_column(0., 1..5, []); + let row = model.add_row(..=1.0, [(x, 1.), (y, 1.)]); + let iis = solve_infeasible(model).get_iis(); + assert_eq!(iis.columns(), &[(x, Lower), (y, Lower)]); + assert_eq!(iis.rows(), &[(row, Upper)]); +} + +#[test] +fn iis_of_feasible_semi_continuous_model_is_empty() { + // x in {0} U [5, 10] with x <= 1: feasible with x = 0. + let mut model = model_with_iis_strategy(None); + let x = model.add_semi_continuous_column(0., 5..10, []); + model.add_row(..=1.0, [(x, 1.)]); + let mut solved = model.solve(); + assert_eq!(solved.status(), HighsModelStatus::Optimal); + assert!(solved.get_iis().is_empty()); +} + +#[test] +fn model_solve_or_iis() { + let mut model = model_with_iis_strategy(None); + let col = model.add_col(1., 1.., []); + let r0 = model.add_row(..=2., [(col, 1.)]); + let solution = model.solve_or_iis().expect("feasible"); + assert_eq!(solution.columns(), &[1.]); + + let r1 = model.add_row(..0.5, [(col, 1.)]); // col >= 1 but row <= 0.5 + let iis = model.solve_or_iis().expect_err("infeasible"); + assert!(iis.contains_column(col)); + assert!(iis.contains_row(r1)); + assert!(!iis.contains_row(r0)); + assert_eq!(model.status(), HighsModelStatus::Infeasible); + assert!(!model.get_iis().is_empty()); +}