diff --git a/src/lib.rs b/src/lib.rs index 10a176e..3c9eea6 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -275,6 +275,83 @@ pub enum Sense { Minimise = OBJECTIVE_SENSE_MINIMIZE as isize, } +/// Storage layout of a quadratic objective Hessian passed to +/// [`Model::pass_hessian`]. +#[derive(Clone, Copy, Eq, PartialEq, Debug)] +pub enum HessianFormat { + /// Only the lower triangle of the symmetric Hessian is stored, in + /// compressed sparse column form. This is the usual way to give the + /// Hessian of `0.5 x' Q x`. + Triangular, + /// The full square Hessian is stored in compressed sparse column form. + Square, +} + +impl HessianFormat { + fn as_raw(self) -> HighsInt { + match self { + HessianFormat::Triangular => kHighsHessianFormatTriangular, + HessianFormat::Square => kHighsHessianFormatSquare, + } + } +} + +/// Reason a Hessian could not be uploaded by [`Model::try_pass_hessian`]. +#[derive(Clone, Copy, PartialEq, Eq, Debug)] +pub enum HessianError { + /// The dimension of `Q` (its number of columns) does not fit in HiGHS' + /// integer type. + DimensionTooLarge { + /// The dimension that was requested. + dim: usize, + }, + /// The number of stored nonzero coefficients does not fit in HiGHS' + /// integer type. + TooManyNonZeros { + /// The number of nonzeros that was requested. + nnz: usize, + }, + /// A row index does not fit in HiGHS' integer type. + IndexTooLarge { + /// Position, among the stored coefficients, of the offending entry. + entry: usize, + }, + /// HiGHS rejected the Hessian and returned this status. + Highs(HighsStatus), +} + +impl std::fmt::Display for HessianError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + let max = HighsInt::MAX; + match *self { + HessianError::DimensionTooLarge { dim } => write!( + f, + "the dimension of the quadratic objective matrix (Hessian Q) is too large: \ + got {dim} but HiGHS supports at most {max}" + ), + HessianError::TooManyNonZeros { nnz } => write!( + f, + "the Hessian Q has too many nonzero coefficients: \ + got {nnz} but HiGHS supports at most {max}" + ), + HessianError::IndexTooLarge { entry } => write!( + f, + "the row index of Hessian coefficient {entry} is too large \ + for HiGHS' integer type (at most {max})" + ), + HessianError::Highs(status) => write!(f, "HiGHS rejected the Hessian: {status:?}"), + } + } +} + +impl std::error::Error for HessianError {} + +impl From for HessianError { + fn from(status: HighsStatus) -> Self { + HessianError::Highs(status) + } +} + impl Model { /// Return pointer to underlying HiGHS model pub fn as_ptr(&self) -> *const c_void { @@ -677,6 +754,117 @@ impl Model { }?; Ok(()) } + + /// Upload a quadratic objective Hessian `Q`, turning the model into a QP + /// with objective `c'x + 0.5 x' Q x` (where `c` is the linear objective + /// already set on the columns). + /// + /// `Q` is provided column by column: `columns` yields one item per column + /// of `Q` and each column yields its stored `(row index, coefficient)` pairs. + /// For [`HessianFormat::Triangular`] store only the lower triangle. + /// For [`HessianFormat::Square`] store the full matrix. + /// Both levels can be anything iterable and the indices may be any integer type + /// that converts to HiGHS' integer type. + /// + /// HiGHS solves **convex** QPs only: `Q` should be positive semidefinite. + /// HiGHS does not check this: + /// on an indefinite `Q` it may return a wrong or + /// non-optimal solution. + /// HiGHS, however, does test for negative diagonal values. + /// Verify convexity yourself if `Q` is not PSD by construction. + /// + /// # Panics + /// + /// If HiGHS returns an error status value, or an index/size does not fit in + /// HiGHS' integer type. Use [`Model::try_pass_hessian`] to handle these as + /// a [`HessianError`] instead. + pub fn pass_hessian(&mut self, format: HessianFormat, columns: C) + where + C: IntoIterator, + E: IntoIterator, + I: TryInto, + { + self.try_pass_hessian(format, columns) + .unwrap_or_else(|e| panic!("pass_hessian failed: {e}")) + } + + /// Same as [`Model::pass_hessian`], but returns a [`HessianError`] instead + /// of panicking. An empty Hessian (no coefficients) is a no-op and leaves + /// the model linear. + /// + /// ``` + /// use highs::{RowProblem, Sense, HessianFormat, HighsModelStatus}; + /// // min x^2 + y^2 s.t. x + y = 1, x, y in [-10, 10] + /// let mut pb = RowProblem::new(); + /// let x = pb.add_column(0.0, -10.0..=10.0); + /// let y = pb.add_column(0.0, -10.0..=10.0); + /// pb.add_row(1.0..=1.0, [(x, 1.0), (y, 1.0)]); + /// let mut model = pb.optimise(Sense::Minimise); + /// // Q = diag(2, 2): one column per variable, each with its diagonal entry. + /// model + /// .try_pass_hessian(HessianFormat::Triangular, [[(0, 2.0)], [(1, 2.0)]]) + /// .unwrap(); + /// let solved = model.solve(); + /// assert_eq!(solved.status(), HighsModelStatus::Optimal); + /// let cols = solved.get_solution().columns().to_vec(); + /// assert!((cols[0] - 0.5).abs() < 1e-6); + /// assert!((cols[1] - 0.5).abs() < 1e-6); + /// ``` + pub fn try_pass_hessian( + &mut self, + format: HessianFormat, + columns: C, + ) -> Result<(), HessianError> + where + C: IntoIterator, + E: IntoIterator, + I: TryInto, + { + // Build the compressed-sparse-column arrays from the per-column + // iterators. + let mut start: Vec = Vec::new(); + let mut index: Vec = Vec::new(); + let mut value: Vec = Vec::new(); + for column in columns { + let offset = index.len(); + start.push( + offset + .try_into() + .map_err(|_| HessianError::TooManyNonZeros { nnz: offset })?, + ); + for (i, v) in column { + let i: HighsInt = i + .try_into() + .map_err(|_| HessianError::IndexTooLarge { entry: index.len() })?; + index.push(i); + value.push(v); + } + } + if value.is_empty() { + return Ok(()); + } + let dim: HighsInt = start + .len() + .try_into() + .map_err(|_| HessianError::DimensionTooLarge { dim: start.len() })?; + let nnz: HighsInt = value + .len() + .try_into() + .map_err(|_| HessianError::TooManyNonZeros { nnz: value.len() })?; + unsafe { + highs_call!(Highs_passHessian( + self.highs.mut_ptr(), + dim, + nnz, + format.as_raw(), + start.as_ptr(), + index.as_ptr(), + value.as_ptr() + )) + } + .map(|_| ()) + .map_err(HessianError::from) + } } impl From for Model { @@ -1097,4 +1285,58 @@ mod test { assert_eq!(status, highs_sys::STATUS_OK); assert_eq!(value, 2); } + + #[test] + fn test_pass_hessian_convex_qp() { + use crate::status::HighsModelStatus::Optimal; + // min x^2 + y^2 s.t. x + y = 1, x, y in [-10, 10]. Optimum at (0.5, 0.5). + let mut pb = RowProblem::new(); + let x = pb.add_column(0., -10.0..=10.0); + let y = pb.add_column(0., -10.0..=10.0); + pb.add_row(1.0..=1.0, [(x, 1.), (y, 1.)]); + let mut model = pb.optimise(Sense::Minimise); + model.make_quiet(); + // Q = diag(2, 2) for the 0.5 x'Qx convention: one column per variable. + model + .try_pass_hessian(HessianFormat::Triangular, [[(0, 2.0)], [(1, 2.0)]]) + .unwrap(); + let solved = model.solve(); + assert_eq!(solved.status(), Optimal); + let cols = solved.get_solution().columns().to_vec(); + assert!((cols[0] - 0.5).abs() < 1e-6, "x = {}", cols[0]); + assert!((cols[1] - 0.5).abs() < 1e-6, "y = {}", cols[1]); + assert!((solved.objective_value() - 0.5).abs() < 1e-6); + } + + #[test] + fn test_pass_hessian_index_overflow_is_error() { + let mut model = RowProblem::default().optimise(Sense::Minimise); + let err = model.try_pass_hessian( + HessianFormat::Triangular, + [vec![(0usize, 2.0)], vec![(usize::MAX, 2.0)]], + ); + assert!(matches!(err, Err(HessianError::IndexTooLarge { entry: 1 }))); + } + + #[test] + fn test_pass_hessian_accepts_lazy_iterators() { + use crate::status::HighsModelStatus::Optimal; + let mut pb = RowProblem::new(); + let x = pb.add_column(0., -10.0..=10.0); + let y = pb.add_column(0., -10.0..=10.0); + pb.add_row(1.0..=1.0, [(x, 1.), (y, 1.)]); + let mut model = pb.optimise(Sense::Minimise); + model.make_quiet(); + model + .try_pass_hessian( + HessianFormat::Triangular, + (0..2).map(|j| std::iter::once((j, 2.0))), + ) + .unwrap(); + let solved = model.solve(); + assert_eq!(solved.status(), Optimal); + let cols = solved.get_solution().columns().to_vec(); + assert!((cols[0] - 0.5).abs() < 1e-6, "x = {}", cols[0]); + assert!((cols[1] - 0.5).abs() < 1e-6, "y = {}", cols[1]); + } }