From 820b0360688c9fe887949399e487a1aaad92fd92 Mon Sep 17 00:00:00 2001 From: Shurikal Date: Fri, 26 Sep 2025 19:45:40 +0200 Subject: [PATCH 1/8] Added first draft for quadratic solver --- Cargo.toml | 1 + src/expression.rs | 8 + src/lib.rs | 13 + src/quadratic_expression.rs | 711 ++++++++++++++++++++++++++++++++ src/quadratic_problem.rs | 96 +++++ src/solvers/clarabel.rs | 443 ++++++++++++++++++++ tests/quadratic_solver_tests.rs | 253 ++++++++++++ 7 files changed, 1525 insertions(+) create mode 100644 src/quadratic_expression.rs create mode 100644 src/quadratic_problem.rs create mode 100644 tests/quadratic_solver_tests.rs diff --git a/Cargo.toml b/Cargo.toml index cc2307f..6dd4bad 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -35,6 +35,7 @@ all_default_solvers = [ minilp = [ "microlp", ] # minilp is not maintained anymore, we use the microlp fork instead +enable_quadratic = ["clarabel"] # Enable quadratic programming support (requires clarabel solver) [dependencies] coin_cbc = { version = "0.1", optional = true, default-features = false } diff --git a/src/expression.rs b/src/expression.rs index 7496fc2..95a83ba 100644 --- a/src/expression.rs +++ b/src/expression.rs @@ -74,6 +74,14 @@ impl FormatWithVars for LinearExpression { } } +impl Clone for LinearExpression { + fn clone(&self) -> Self { + LinearExpression { + coefficients: self.coefficients.clone(), + } + } +} + /// Represents an affine expression, such as `2x + 3` or `x + y + z` pub struct Expression { pub(crate) linear: LinearExpression, diff --git a/src/lib.rs b/src/lib.rs index c2d9959..2a1bc9d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -69,6 +69,10 @@ pub use affine_expression_trait::IntoAffineExpression; pub use cardinality_constraint_solver_trait::CardinalityConstraintSolver; pub use constraint::Constraint; pub use expression::Expression; +#[cfg(feature = "enable_quadratic")] +pub use quadratic_expression::{QuadraticAffineExpression, QuadraticConstraint, QuadraticExpression, VariablePair}; +#[cfg(feature = "enable_quadratic")] +pub use quadratic_problem::{QuadraticUnsolvedProblem, IntoQuadraticExpression}; #[cfg(not(any( feature = "coin_cbc", feature = "microlp", @@ -82,6 +86,9 @@ pub use solvers::clarabel::clarabel as default_solver; #[cfg_attr(docsrs, doc(cfg(feature = "clarabel")))] #[cfg(feature = "clarabel")] pub use solvers::clarabel::clarabel; +#[cfg_attr(docsrs, doc(cfg(all(feature = "clarabel", feature = "enable_quadratic"))))] +#[cfg(all(feature = "clarabel", feature = "enable_quadratic"))] +pub use solvers::clarabel::clarabel_quadratic; #[cfg_attr(docsrs, doc(cfg(feature = "coin_cbc")))] #[cfg(feature = "coin_cbc")] pub use solvers::coin_cbc::coin_cbc; @@ -190,3 +197,9 @@ mod cardinality_constraint_solver_trait; pub mod constraint; pub mod solvers; mod variables_macro; +/// Quadratic expression types and operations +#[cfg(feature = "enable_quadratic")] +pub mod quadratic_expression; +/// Quadratic problem definition and solving +#[cfg(feature = "enable_quadratic")] +pub mod quadratic_problem; diff --git a/src/quadratic_expression.rs b/src/quadratic_expression.rs new file mode 100644 index 0000000..ca983ef --- /dev/null +++ b/src/quadratic_expression.rs @@ -0,0 +1,711 @@ +use std::fmt::{Debug, Formatter}; +use std::ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub, SubAssign}; + +use fnv::FnvHashMap as HashMap; + +use crate::affine_expression_trait::IntoAffineExpression; + +use crate::expression::{Expression, LinearExpression}; +use crate::variable::{FormatWithVars, Variable}; +use crate::Solution; + +/// Represents a quadratic term as a pair of variables and their coefficient +/// For example, 3*x*y would be represented as QuadraticTerm { var1: x, var2: y, coeff: 3.0 } +/// For x^2, both var1 and var2 would be x +#[derive(Clone, Copy, PartialEq, Eq, Hash, Debug)] +pub struct VariablePair { + /// First variable in the pair + pub var1: Variable, + /// Second variable in the pair + pub var2: Variable, +} + +impl VariablePair { + /// Create a new variable pair, ensuring consistent ordering for commutativity (x*y = y*x) + pub fn new(var1: Variable, var2: Variable) -> Self { + // Ensure consistent ordering for commutativity (x*y = y*x) + if var1.index() <= var2.index() { + VariablePair { var1, var2 } + } else { + VariablePair { var1: var2, var2: var1 } + } + } + + /// Create a variable pair representing a square term (var^2) + pub fn square(var: Variable) -> Self { + VariablePair { var1: var, var2: var } + } +} + +/// A quadratic expression without linear or constant components +pub struct QuadraticExpression { + pub(crate) quadratic_coefficients: HashMap, +} + +impl QuadraticExpression { + /// Create a new empty quadratic expression + pub fn new() -> Self { + QuadraticExpression { + quadratic_coefficients: HashMap::default(), + } + } + + /// Create a quadratic expression with preallocated capacity + pub fn with_capacity(capacity: usize) -> Self { + QuadraticExpression { + quadratic_coefficients: HashMap::with_capacity_and_hasher(capacity, Default::default()), + } + } + + /// Add a quadratic term to this expression + pub fn add_quadratic_term(&mut self, var1: Variable, var2: Variable, coefficient: f64) { + let pair = VariablePair::new(var1, var2); + *self.quadratic_coefficients.entry(pair).or_default() += coefficient; + } + + /// Get the coefficient for a specific variable pair + pub fn get_quadratic_coefficient(&self, var1: Variable, var2: Variable) -> f64 { + let pair = VariablePair::new(var1, var2); + self.quadratic_coefficients.get(&pair).copied().unwrap_or(0.0) + } + + /// Check if the expression is empty (all coefficients are zero) + pub fn is_empty(&self) -> bool { + self.quadratic_coefficients.is_empty() || + self.quadratic_coefficients.values().all(|&c| c.abs() < f64::EPSILON) + } + + /// Iterate over all non-zero quadratic terms + pub fn iter(&self) -> impl Iterator + '_ { + self.quadratic_coefficients.iter().filter_map(|(&pair, &coeff)| { + if coeff.abs() >= f64::EPSILON { + Some((pair, coeff)) + } else { + None + } + }) + } + + /// Evaluate the quadratic expression given variable values + pub fn eval_with(&self, values: &S) -> f64 { + self.quadratic_coefficients + .iter() + .map(|(pair, &coeff)| { + let val1 = values.value(pair.var1); + let val2 = values.value(pair.var2); + coeff * val1 * val2 + }) + .sum() + } +} + +impl Default for QuadraticExpression { + fn default() -> Self { + Self::new() + } +} + +impl Clone for QuadraticExpression { + fn clone(&self) -> Self { + QuadraticExpression { + quadratic_coefficients: self.quadratic_coefficients.clone(), + } + } +} + +/// A complete quadratic expression containing quadratic, linear, and constant terms +/// Represents expressions like: x^2 + 2*x*y + 3*x + 4 +pub struct QuadraticAffineExpression { + pub(crate) quadratic: QuadraticExpression, + pub(crate) linear: LinearExpression, + pub(crate) constant: f64, +} + +impl QuadraticAffineExpression { + /// Create a new empty quadratic affine expression + pub fn new() -> Self { + QuadraticAffineExpression { + quadratic: QuadraticExpression::new(), + linear: LinearExpression { + coefficients: HashMap::default(), + }, + constant: 0.0, + } + } + + /// Create a quadratic affine expression with preallocated capacity + pub fn with_capacity(linear_capacity: usize, quadratic_capacity: usize) -> Self { + QuadraticAffineExpression { + quadratic: QuadraticExpression::with_capacity(quadratic_capacity), + linear: LinearExpression { + coefficients: HashMap::with_capacity_and_hasher(linear_capacity, Default::default()), + }, + constant: 0.0, + } + } + + /// Create from an affine expression + pub fn from_affine(expr: E) -> Self { + let constant = expr.constant(); + let linear_coeffs = expr.linear_coefficients().into_iter().collect(); + QuadraticAffineExpression { + quadratic: QuadraticExpression::new(), + linear: LinearExpression { coefficients: linear_coeffs }, + constant, + } + } + + /// Create from a quadratic expression + pub fn from_quadratic(expr: QuadraticExpression) -> Self { + QuadraticAffineExpression { + quadratic: expr, + linear: LinearExpression { + coefficients: HashMap::default(), + }, + constant: 0.0, + } + } + + /// Add a quadratic term + pub fn add_quadratic_term(&mut self, var1: Variable, var2: Variable, coefficient: f64) { + self.quadratic.add_quadratic_term(var1, var2, coefficient); + } + + /// Add a linear term + pub fn add_linear_term(&mut self, var: Variable, coefficient: f64) { + *self.linear.coefficients.entry(var).or_default() += coefficient; + } + + /// Add a constant term + pub fn add_constant(&mut self, value: f64) { + self.constant += value; + } + + /// Get the linear part as an Expression + pub fn linear_part(&self) -> Expression { + Expression { + linear: LinearExpression { + coefficients: self.linear.coefficients.clone(), + }, + constant: self.constant, + } + } + + /// Check if the expression contains any quadratic terms + pub fn is_purely_affine(&self) -> bool { + self.quadratic.is_empty() + } + + /// Evaluate the complete expression + pub fn eval_with(&self, values: &S) -> f64 { + self.quadratic.eval_with(values) + + self.linear_part().eval_with(values) + } + + /// Get the coefficient of a quadratic term + pub fn get_quadratic_coefficient(&self, var1: Variable, var2: Variable) -> f64 { + self.quadratic.get_quadratic_coefficient(var1, var2) + } + + /// Get the coefficient of a linear term + pub fn get_linear_coefficient(&self, var: Variable) -> f64 { + self.linear.coefficients.get(&var).copied().unwrap_or(0.0) + } + + /// Get the constant term + pub fn get_constant(&self) -> f64 { + self.constant + } + + /// Get a reference to the quadratic part + pub fn quadratic_part(&self) -> &QuadraticExpression { + &self.quadratic + } + + /// Creates a constraint indicating that this quadratic expression + /// is lesser than or equal to the right hand side + pub fn leq(self, rhs: RHS) -> QuadraticConstraint + where + Self: Sub, + { + let diff = self - rhs; + QuadraticConstraint { + expression: diff, + is_equality: false, + } + } + + /// Creates a constraint indicating that this quadratic expression + /// is greater than or equal to the right hand side + pub fn geq(self, rhs: RHS) -> QuadraticConstraint + where + RHS: Sub, + { + let diff = rhs - self; + QuadraticConstraint { + expression: diff, + is_equality: false, + } + } + + /// Creates a constraint indicating that this quadratic expression + /// is equal to the right hand side + pub fn eq(self, rhs: RHS) -> QuadraticConstraint + where + Self: Sub, + { + let diff = self - rhs; + QuadraticConstraint { + expression: diff, + is_equality: true, + } + } +} + +impl Default for QuadraticAffineExpression { + fn default() -> Self { + Self::new() + } +} + +impl Clone for QuadraticAffineExpression { + fn clone(&self) -> Self { + QuadraticAffineExpression { + quadratic: self.quadratic.clone(), + linear: self.linear.clone(), + constant: self.constant, + } + } +} + +/// A quadratic constraint +pub struct QuadraticConstraint { + /// The quadratic expression that defines the constraint + pub expression: QuadraticAffineExpression, + /// Whether this is an equality constraint (true) or inequality constraint (false) + pub is_equality: bool, +} + +// Conversion implementations +impl From for QuadraticAffineExpression { + fn from(value: f64) -> Self { + QuadraticAffineExpression { + quadratic: QuadraticExpression::new(), + linear: LinearExpression { + coefficients: HashMap::default(), + }, + constant: value, + } + } +} + +impl From for QuadraticAffineExpression { + fn from(var: Variable) -> Self { + let mut linear_coeffs = HashMap::default(); + linear_coeffs.insert(var, 1.0); + QuadraticAffineExpression { + quadratic: QuadraticExpression::new(), + linear: LinearExpression { coefficients: linear_coeffs }, + constant: 0.0, + } + } +} + +impl From for QuadraticAffineExpression { + fn from(expr: Expression) -> Self { + QuadraticAffineExpression { + quadratic: QuadraticExpression::new(), + linear: expr.linear, + constant: expr.constant, + } + } +} + +impl From for QuadraticAffineExpression { + fn from(quad_expr: QuadraticExpression) -> Self { + QuadraticAffineExpression::from_quadratic(quad_expr) + } +} + +// Arithmetic operations for QuadraticExpression +impl> MulAssign for QuadraticExpression { + fn mul_assign(&mut self, rhs: N) { + let factor = rhs.into(); + for value in self.quadratic_coefficients.values_mut() { + *value *= factor; + } + } +} + +impl> Mul for QuadraticExpression { + type Output = QuadraticExpression; + + fn mul(mut self, rhs: N) -> Self::Output { + self.mul_assign(rhs); + self + } +} + +impl> Div for QuadraticExpression { + type Output = QuadraticExpression; + + fn div(mut self, rhs: N) -> Self::Output { + self.mul_assign(1. / rhs.into()); + self + } +} + +impl AddAssign for QuadraticExpression { + fn add_assign(&mut self, rhs: QuadraticExpression) { + for (pair, coeff) in rhs.quadratic_coefficients { + *self.quadratic_coefficients.entry(pair).or_default() += coeff; + } + } +} + +impl Add for QuadraticExpression { + type Output = QuadraticExpression; + + fn add(mut self, rhs: QuadraticExpression) -> Self::Output { + self.add_assign(rhs); + self + } +} + +impl SubAssign for QuadraticExpression { + fn sub_assign(&mut self, rhs: QuadraticExpression) { + for (pair, coeff) in rhs.quadratic_coefficients { + *self.quadratic_coefficients.entry(pair).or_default() -= coeff; + } + } +} + +impl Sub for QuadraticExpression { + type Output = QuadraticExpression; + + fn sub(mut self, rhs: QuadraticExpression) -> Self::Output { + self.sub_assign(rhs); + self + } +} + +impl Neg for QuadraticExpression { + type Output = Self; + + fn neg(mut self) -> Self::Output { + self *= -1; + self + } +} + +// Arithmetic operations for QuadraticAffineExpression +impl> MulAssign for QuadraticAffineExpression { + fn mul_assign(&mut self, rhs: N) { + let factor = rhs.into(); + self.quadratic.mul_assign(factor); + for value in self.linear.coefficients.values_mut() { + *value *= factor; + } + self.constant *= factor; + } +} + +impl> Mul for QuadraticAffineExpression { + type Output = QuadraticAffineExpression; + + fn mul(mut self, rhs: N) -> Self::Output { + self.mul_assign(rhs); + self + } +} + +impl> Div for QuadraticAffineExpression { + type Output = QuadraticAffineExpression; + + fn div(mut self, rhs: N) -> Self::Output { + self.mul_assign(1. / rhs.into()); + self + } +} + +impl AddAssign for QuadraticAffineExpression { + fn add_assign(&mut self, rhs: QuadraticAffineExpression) { + self.quadratic.add_assign(rhs.quadratic); + for (var, coeff) in rhs.linear.coefficients { + *self.linear.coefficients.entry(var).or_default() += coeff; + } + self.constant += rhs.constant; + } +} + +impl Add for QuadraticAffineExpression { + type Output = QuadraticAffineExpression; + + fn add(mut self, rhs: QuadraticAffineExpression) -> Self::Output { + self.add_assign(rhs); + self + } +} + +impl SubAssign for QuadraticAffineExpression { + fn sub_assign(&mut self, rhs: QuadraticAffineExpression) { + self.quadratic.sub_assign(rhs.quadratic); + for (var, coeff) in rhs.linear.coefficients { + *self.linear.coefficients.entry(var).or_default() -= coeff; + } + self.constant -= rhs.constant; + } +} + +impl Sub for QuadraticAffineExpression { + type Output = QuadraticAffineExpression; + + fn sub(mut self, rhs: QuadraticAffineExpression) -> Self::Output { + self.sub_assign(rhs); + self + } +} + +impl Neg for QuadraticAffineExpression { + type Output = Self; + + fn neg(mut self) -> Self::Output { + self *= -1; + self + } +} + +// Mixed type operations +impl AddAssign for QuadraticAffineExpression { + fn add_assign(&mut self, rhs: E) { + let constant = rhs.constant(); + for (var, coeff) in rhs.linear_coefficients() { + *self.linear.coefficients.entry(var).or_default() += coeff; + } + self.constant += constant; + } +} + +impl Add for QuadraticAffineExpression { + type Output = QuadraticAffineExpression; + + fn add(mut self, rhs: E) -> Self::Output { + self.add_assign(rhs); + self + } +} + +impl SubAssign for QuadraticAffineExpression { + fn sub_assign(&mut self, rhs: E) { + let constant = rhs.constant(); + for (var, coeff) in rhs.linear_coefficients() { + *self.linear.coefficients.entry(var).or_default() -= coeff; + } + self.constant -= constant; + } +} + +impl Sub for QuadraticAffineExpression { + type Output = QuadraticAffineExpression; + + fn sub(mut self, rhs: E) -> Self::Output { + self.sub_assign(rhs); + self + } +} + +// Multiplication operations that create quadratic terms +impl Mul for Variable { + type Output = QuadraticExpression; + + fn mul(self, rhs: Variable) -> Self::Output { + let mut quadratic = QuadraticExpression::new(); + quadratic.add_quadratic_term(self, rhs, 1.0); + quadratic + } +} + +// Specific implementations for Variable multiplication with other types +impl Mul for Variable { + type Output = QuadraticAffineExpression; + + fn mul(self, rhs: Expression) -> Self::Output { + let mut result = QuadraticAffineExpression::new(); + let constant = rhs.constant; + + // Linear term (constant * variable) + if constant.abs() >= f64::EPSILON { + result.add_linear_term(self, constant); + } + + // Quadratic terms (variable * other_variables) + for (var, coeff) in rhs.linear.coefficients { + result.add_quadratic_term(self, var, coeff); + } + + result + } +} + +impl Mul for Expression { + type Output = QuadraticAffineExpression; + + fn mul(self, rhs: Variable) -> Self::Output { + rhs * self + } +} + +// Expression multiplication (returns quadratic expression) +impl Mul for Expression { + type Output = QuadraticAffineExpression; + + fn mul(self, rhs: Expression) -> Self::Output { + let mut result = QuadraticAffineExpression::new(); + + let c1 = self.constant; + let c2 = rhs.constant; + + // Constant term (c1 * c2) + result.constant = c1 * c2; + + // Linear terms from self * constant of rhs + for (var, coeff) in self.linear.coefficients.iter() { + result.add_linear_term(*var, *coeff * c2); + } + + // Linear terms from constant of self * rhs + for (var, coeff) in rhs.linear.coefficients.iter() { + result.add_linear_term(*var, c1 * *coeff); + } + + // Quadratic terms (variable from self * variable from rhs) + for (var1, coeff1) in self.linear.coefficients.iter() { + for (var2, coeff2) in rhs.linear.coefficients.iter() { + result.add_quadratic_term(*var1, *var2, *coeff1 * *coeff2); + } + } + + result + } +} + +// Numeric multiplication +macro_rules! impl_mul_numeric { + ($($t:ty),*) => {$( + impl Mul for $t { + type Output = QuadraticExpression; + + fn mul(self, mut rhs: QuadraticExpression) -> Self::Output { + rhs *= self; + rhs + } + } + + impl Mul for $t { + type Output = QuadraticAffineExpression; + + fn mul(self, mut rhs: QuadraticAffineExpression) -> Self::Output { + rhs *= self; + rhs + } + } + )*} +} +impl_mul_numeric!(f64, i32, f32, u32, u16, u8, i16, i8); + +// Display and formatting support +impl FormatWithVars for QuadraticExpression { + fn format_with(&self, f: &mut Formatter<'_>, mut variable_format: FUN) -> std::fmt::Result + where + FUN: FnMut(&mut Formatter<'_>, Variable) -> std::fmt::Result, + { + let mut first = true; + for (pair, coeff) in self.iter() { + if first { + first = false; + } else { + write!(f, " + ")?; + } + + if (coeff - 1.).abs() > f64::EPSILON { + write!(f, "{} ", coeff)?; + } + + if pair.var1 == pair.var2 { + // Square term: x^2 + variable_format(f, pair.var1)?; + write!(f, "²")?; + } else { + // Cross term: x*y + variable_format(f, pair.var1)?; + write!(f, "*")?; + variable_format(f, pair.var2)?; + } + } + + if first { + write!(f, "0")?; + } + Ok(()) + } +} + +impl Debug for QuadraticExpression { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + self.format_debug(f) + } +} + +impl FormatWithVars for QuadraticAffineExpression { + fn format_with(&self, f: &mut Formatter<'_>, mut variable_format: FUN) -> std::fmt::Result + where + FUN: FnMut(&mut Formatter<'_>, Variable) -> std::fmt::Result, + { + let has_quadratic = !self.quadratic.is_empty(); + let has_linear = !self.linear.coefficients.is_empty(); + let has_constant = self.constant.abs() >= f64::EPSILON; + + if !has_quadratic && !has_linear && !has_constant { + return write!(f, "0"); + } + + let mut first = true; + + // Write quadratic terms + if has_quadratic { + self.quadratic.format_with(f, &mut variable_format)?; + first = false; + } + + // Write linear terms + for (&var, &coeff) in &self.linear.coefficients { + if coeff.abs() >= f64::EPSILON { + if !first { + write!(f, " + ")?; + } + first = false; + + if (coeff - 1.).abs() > f64::EPSILON { + write!(f, "{} ", coeff)?; + } + variable_format(f, var)?; + } + } + + // Write constant term + if has_constant { + if !first { + write!(f, " + ")?; + } + write!(f, "{}", self.constant)?; + } + + Ok(()) + } +} + +impl Debug for QuadraticAffineExpression { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + self.format_debug(f) + } +} \ No newline at end of file diff --git a/src/quadratic_problem.rs b/src/quadratic_problem.rs new file mode 100644 index 0000000..be25a41 --- /dev/null +++ b/src/quadratic_problem.rs @@ -0,0 +1,96 @@ +use crate::quadratic_expression::QuadraticAffineExpression; +use crate::variable::ProblemVariables; +use crate::solvers::ObjectiveDirection; + +/// A problem with a quadratic objective function +pub struct QuadraticUnsolvedProblem { + /// The quadratic objective function to optimize + pub objective: QuadraticAffineExpression, + /// Whether to maximize or minimize the objective + pub direction: ObjectiveDirection, + /// The variables in the problem + pub variables: ProblemVariables, +} + +/// Trait for quadratic expressions that can be used as objective functions +pub trait IntoQuadraticExpression { + /// Transform the value into a concrete QuadraticAffineExpression + fn into_quadratic_expression(self) -> QuadraticAffineExpression; +} + +impl IntoQuadraticExpression for QuadraticAffineExpression { + fn into_quadratic_expression(self) -> QuadraticAffineExpression { + self + } +} + +impl IntoQuadraticExpression for crate::QuadraticExpression { + fn into_quadratic_expression(self) -> QuadraticAffineExpression { + QuadraticAffineExpression::from_quadratic(self) + } +} + +impl IntoQuadraticExpression for T { + fn into_quadratic_expression(self) -> QuadraticAffineExpression { + QuadraticAffineExpression::from_affine(self) + } +} + +impl ProblemVariables { + /// Create a quadratic optimization problem to maximize the given quadratic objective + #[cfg(feature = "enable_quadratic")] + pub fn maximise_quadratic( + self, + objective: E, + ) -> QuadraticUnsolvedProblem { + let objective_expr = objective.into_quadratic_expression(); + // Validate that objective variables are within bounds + for (var, _) in objective_expr.linear.coefficients.iter() { + assert!(var.index() < self.len(), + "Variable in objective function is not part of this problem"); + } + for (pair, _) in objective_expr.quadratic.quadratic_coefficients.iter() { + assert!(pair.var1.index() < self.len() && pair.var2.index() < self.len(), + "Variable in quadratic objective function is not part of this problem"); + } + QuadraticUnsolvedProblem { + objective: objective_expr, + direction: ObjectiveDirection::Maximisation, + variables: self, + } + } + + /// Create a quadratic optimization problem to minimize the given quadratic objective + #[cfg(feature = "enable_quadratic")] + pub fn minimise_quadratic( + self, + objective: E, + ) -> QuadraticUnsolvedProblem { + let objective_expr = objective.into_quadratic_expression(); + // Validate that objective variables are within bounds + for (var, _) in objective_expr.linear.coefficients.iter() { + assert!(var.index() < self.len(), + "Variable in objective function is not part of this problem"); + } + for (pair, _) in objective_expr.quadratic.quadratic_coefficients.iter() { + assert!(pair.var1.index() < self.len() && pair.var2.index() < self.len(), + "Variable in quadratic objective function is not part of this problem"); + } + QuadraticUnsolvedProblem { + objective: objective_expr, + direction: ObjectiveDirection::Minimisation, + variables: self, + } + } +} + +impl QuadraticUnsolvedProblem { + /// Solve this quadratic problem using the given solver function + #[cfg(all(feature = "clarabel", feature = "enable_quadratic"))] + pub fn using(self, solver_factory: F) -> S + where + F: FnOnce(Self) -> S, + { + solver_factory(self) + } +} \ No newline at end of file diff --git a/src/solvers/clarabel.rs b/src/solvers/clarabel.rs index d47e849..da14a02 100644 --- a/src/solvers/clarabel.rs +++ b/src/solvers/clarabel.rs @@ -3,6 +3,10 @@ use crate::affine_expression_trait::IntoAffineExpression; use crate::expression::LinearExpression; use crate::variable::UnsolvedProblem; +#[cfg(feature = "enable_quadratic")] +use crate::quadratic_problem::QuadraticUnsolvedProblem; +#[cfg(feature = "enable_quadratic")] +use crate::quadratic_expression::VariablePair; use crate::{ constraint::ConstraintReference, solvers::{ObjectiveDirection, ResolutionError, Solution, SolverModel}, @@ -59,6 +63,61 @@ pub fn clarabel(to_solve: UnsolvedProblem) -> ClarabelProblem { p } +/// The [clarabel](https://clarabel.org/stable/) solver for quadratic problems, +/// to be used with [QuadraticUnsolvedProblem::using]. +#[cfg(feature = "enable_quadratic")] +pub fn clarabel_quadratic(to_solve: QuadraticUnsolvedProblem) -> ClarabelQuadraticProblem { + let QuadraticUnsolvedProblem { + objective, + direction, + variables, + } = to_solve; + let coef = if direction == ObjectiveDirection::Maximisation { + -1. + } else { + 1. + }; + + // Linear objective vector + let mut objective_vector = vec![0.; variables.len()]; + for (var, obj_coeff) in objective.linear.coefficients { + objective_vector[var.index()] = obj_coeff * coef; + } + + // Quadratic objective matrix (P matrix) + let mut quadratic_matrix_builder = CscQuadraticMatrixBuilder::new(variables.len()); + for (pair, quad_coeff) in objective.quadratic.quadratic_coefficients { + quadratic_matrix_builder.add_term(pair, quad_coeff * coef); + } + + let constraints_matrix_builder = CscMatrixBuilder::new(variables.len()); + let mut settings = DefaultSettingsBuilder::default(); + settings.verbose(false).tol_feas(1e-9); + let mut p = ClarabelQuadraticProblem { + objective: objective_vector, + quadratic_matrix_builder, + constraints_matrix_builder, + constraint_values: Vec::new(), + variables: variables.len(), + settings, + cones: Vec::new(), + }; + + // add trivial constraints embedded in the variable definitions + for (var, def) in variables.iter_variables_with_def() { + if def.is_integer { + panic!("Clarabel doesn't support integer variables") + } + if def.min != f64::NEG_INFINITY { + p.add_constraint(var >> def.min); + } + if def.max != f64::INFINITY { + p.add_constraint(var << def.max); + } + } + p +} + /// A clarabel model pub struct ClarabelProblem { constraints_matrix_builder: CscMatrixBuilder, @@ -69,6 +128,18 @@ pub struct ClarabelProblem { cones: Vec>, } +/// A clarabel quadratic model +#[cfg(feature = "enable_quadratic")] +pub struct ClarabelQuadraticProblem { + constraints_matrix_builder: CscMatrixBuilder, + quadratic_matrix_builder: CscQuadraticMatrixBuilder, + constraint_values: Vec, + objective: Vec, + variables: usize, + settings: DefaultSettingsBuilder, + cones: Vec>, +} + impl ClarabelProblem { /// Access the problem settings pub fn settings(&mut self) -> &mut DefaultSettingsBuilder { @@ -247,6 +318,157 @@ fn fast_flatten_vecs(vecs: Vec>) -> Vec { result } +#[cfg(feature = "enable_quadratic")] +impl ClarabelQuadraticProblem { + /// Access the problem settings + pub fn settings(&mut self) -> &mut DefaultSettingsBuilder { + &mut self.settings + } + + /// Convert the problem into a clarabel solver + /// panics if the problem is not valid + pub fn into_solver(self) -> DefaultSolver { + let settings = self.settings.build().expect("Invalid clarabel settings"); + let quadratic_objective = &self.quadratic_matrix_builder.build(); + let objective = &self.objective; + let constraints = &self.constraints_matrix_builder.build(); + let constraint_values = &self.constraint_values; + let cones = &self.cones; + DefaultSolver::new( + quadratic_objective, + objective, + constraints, + constraint_values, + cones, + settings, + ).expect("Invalid clarabel problem. This is likely a bug in good_lp. Problems should always have coherent dimensions.") + } +} + +#[cfg(feature = "enable_quadratic")] +impl SolverModel for ClarabelQuadraticProblem { + type Solution = ClarabelSolution; + type Error = ResolutionError; + + fn solve(self) -> Result { + let mut solver = self.into_solver(); + solver.solve(); + match solver.solution.status { + e @ (SolverStatus::PrimalInfeasible | SolverStatus::AlmostPrimalInfeasible) => { + eprintln!("Clarabel error: {:?}", e); + Err(ResolutionError::Infeasible) + } + SolverStatus::Solved + | SolverStatus::AlmostSolved + | SolverStatus::AlmostDualInfeasible + | SolverStatus::DualInfeasible => Ok(ClarabelSolution { + solution: solver.solution, + }), + SolverStatus::Unsolved => Err(ResolutionError::Other("Unsolved")), + SolverStatus::MaxIterations => Err(ResolutionError::Other("Max iterations reached")), + SolverStatus::MaxTime => Err(ResolutionError::Other("Time limit reached")), + SolverStatus::NumericalError => Err(ResolutionError::Other("Numerical error")), + SolverStatus::InsufficientProgress => Err(ResolutionError::Other("No progress")), + SolverStatus::CallbackTerminated => Err(ResolutionError::Other("Callback terminated")), + } + } + + fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { + self.constraints_matrix_builder + .add_row(constraint.expression.linear); + let index = self.constraint_values.len(); + self.constraint_values.push(-constraint.expression.constant); + // Cones indicate the type of constraint. We only support nonnegative and equality constraints. + // To avoid creating a new cone for each constraint, we merge them. + let next_cone = if constraint.is_equality { + ZeroConeT(1) + } else { + NonnegativeConeT(1) + }; + let prev_cone = self.cones.last_mut(); + match (prev_cone, next_cone) { + (Some(ZeroConeT(a)), ZeroConeT(b)) => *a += b, + (Some(NonnegativeConeT(a)), NonnegativeConeT(b)) => *a += b, + (_, next_cone) => self.cones.push(next_cone), + }; + ConstraintReference { index } + } + + fn name() -> &'static str { + "Clarabel Quadratic" + } +} + +/// Builder for symmetric quadratic matrices in CSC format +#[cfg(feature = "enable_quadratic")] +struct CscQuadraticMatrixBuilder { + /// Indicates the row index of the corresponding element in `nzval` + rowval: Vec>, + /// All non-zero values in the matrix, in column-major order + nzval: Vec>, + n_vars: usize, +} + +#[cfg(feature = "enable_quadratic")] +impl CscQuadraticMatrixBuilder { + fn new(n_vars: usize) -> Self { + Self { + rowval: vec![Vec::new(); n_vars], + nzval: vec![Vec::new(); n_vars], + n_vars, + } + } + + fn add_term(&mut self, pair: VariablePair, coeff: f64) { + let i = pair.var1.index(); + let j = pair.var2.index(); + + if i == j { + // Diagonal term: multiply by 2 because Clarabel minimizes 1/2 * x^T * P * x + // So we need P_ii = 2 * coeff to get coeff * x_i^2 in the objective + self.rowval[j].push(i); + self.nzval[j].push(2.0 * coeff); + } else { + // Off-diagonal term: the coefficient coeff represents the full cross term coeff * x_i * x_j + // In matrix form, this becomes coeff/2 * (P_ij + P_ji) * x_i * x_j + // Since P is symmetric (P_ij = P_ji), we get coeff * P_ij * x_i * x_j + // But Clarabel uses 1/2 * x^T * P * x, so we need P_ij = 2 * coeff + self.rowval[j].push(i); + self.nzval[j].push(coeff); + + self.rowval[i].push(j); + self.nzval[i].push(coeff); + } + } + + fn build(mut self) -> clarabel::algebra::CscMatrix { + // Sort each column by row index to maintain proper CSC format + for col in 0..self.n_vars { + let mut pairs: Vec<(usize, f64)> = self.rowval[col].iter() + .zip(self.nzval[col].iter()) + .map(|(&row, &val)| (row, val)) + .collect(); + pairs.sort_by_key(|&(row, _)| row); + + self.rowval[col] = pairs.iter().map(|&(row, _)| row).collect(); + self.nzval[col] = pairs.iter().map(|&(_, val)| val).collect(); + } + + let mut colptr = Vec::with_capacity(self.n_vars + 1); + colptr.push(0); + for col in &self.rowval { + colptr.push(colptr.last().unwrap() + col.len()); + } + clarabel::algebra::CscMatrix::new( + self.n_vars, + self.n_vars, + colptr, + fast_flatten_vecs(self.rowval), + fast_flatten_vecs(self.nzval), + ) + } +} + #[cfg(test)] mod tests { @@ -277,4 +499,225 @@ mod tests { assert_eq!(matrix.get_entry((1, 1)), Some(4.)); assert_eq!(matrix.get_entry((1, 2)), Some(5.)); } + + #[cfg(feature = "enable_quadratic")] + mod quadratic_tests { + use super::*; + use crate::{QuadraticAffineExpression, VariablePair}; + + #[test] + fn test_csc_quadratic_matrix_builder_diagonal() { + variables! {vars: + x; + y; + } + + let mut builder = CscQuadraticMatrixBuilder::new(2); + + // Add diagonal terms: x^2 with coefficient 2.0, y^2 with coefficient 3.0 + builder.add_term(VariablePair::new(x, x), 2.0); + builder.add_term(VariablePair::new(y, y), 3.0); + + let matrix = builder.build(); + + // Check diagonal entries (scaled by 2 for Clarabel's 1/2 * x^T * P * x formulation) + assert_eq!(matrix.get_entry((0, 0)), Some(4.0)); // 2.0 * 2 = 4.0 + assert_eq!(matrix.get_entry((1, 1)), Some(6.0)); // 3.0 * 2 = 6.0 + + // Check off-diagonal entries are zero + assert_eq!(matrix.get_entry((0, 1)), None); + assert_eq!(matrix.get_entry((1, 0)), None); + } + + #[test] + fn test_csc_quadratic_matrix_builder_off_diagonal() { + variables! {vars: + x; + y; + } + + let mut builder = CscQuadraticMatrixBuilder::new(3); + + // Add off-diagonal term: 2*x*y (coefficient 2.0) + builder.add_term(VariablePair::new(x, y), 2.0); + + let matrix = builder.build(); + + // Off-diagonal terms: coefficient is stored directly (no factor of 2 for off-diagonal) + assert_eq!(matrix.get_entry((0, 1)), Some(2.0)); // (x,y) entry + assert_eq!(matrix.get_entry((1, 0)), Some(2.0)); // (y,x) entry (symmetric) + + // Other entries should be empty + assert_eq!(matrix.get_entry((0, 0)), None); + assert_eq!(matrix.get_entry((1, 1)), None); + assert_eq!(matrix.get_entry((2, 2)), None); + } + + #[test] + fn test_csc_quadratic_matrix_builder_mixed() { + variables! {vars: + x; + y; + } + + let mut builder = CscQuadraticMatrixBuilder::new(2); + + // Add mixed terms: x^2 + 2*x*y + y^2 (like (x+y)^2) + builder.add_term(VariablePair::new(x, x), 1.0); + builder.add_term(VariablePair::new(x, y), 2.0); + builder.add_term(VariablePair::new(y, y), 1.0); + + let matrix = builder.build(); + + // Check the resulting matrix (accounting for Clarabel's 1/2 factor): + // For (x+y)^2 = x^2 + 2*x*y + y^2, we expect: + // [2.0 2.0] (diagonal terms scaled by 2, off-diagonal terms as-is) + // [2.0 2.0] + assert_eq!(matrix.get_entry((0, 0)), Some(2.0)); // x^2 (scaled by 2) + assert_eq!(matrix.get_entry((0, 1)), Some(2.0)); // x*y coefficient + assert_eq!(matrix.get_entry((1, 0)), Some(2.0)); // y*x coefficient + assert_eq!(matrix.get_entry((1, 1)), Some(2.0)); // y^2 (scaled by 2) + } + + #[test] + fn test_simple_quadratic_problem() { + variables! {vars: + x; + y; + } + + // Create a simple quadratic objective: minimize x^2 + y^2 + let mut quadratic_obj = QuadraticAffineExpression::new(); + quadratic_obj.add_quadratic_term(x, x, 1.0); // x^2 + quadratic_obj.add_quadratic_term(y, y, 1.0); // y^2 + + // Create the problem + let problem = vars + .minimise_quadratic(quadratic_obj) + .using(clarabel_quadratic); + + // The problem should be solvable (unconstrained minimum at (0,0)) + let solution = problem.solve().expect("Should solve successfully"); + + // Check that solution is near the origin + assert!((solution.value(x)).abs() < 1e-6); + assert!((solution.value(y)).abs() < 1e-6); + } + + #[test] + fn test_quadratic_problem_with_constraints() { + variables! {vars: + x; + y; + } + + // Create quadratic objective: minimize x^2 + y^2 + let mut quadratic_obj = QuadraticAffineExpression::new(); + quadratic_obj.add_quadratic_term(x, x, 1.0); // x^2 + quadratic_obj.add_quadratic_term(y, y, 1.0); // y^2 + + // Solve with constraint x + y >= 2 + let solution = vars + .minimise_quadratic(quadratic_obj) + .using(clarabel_quadratic) + .with(x + y >> 2.0) + .solve() + .expect("Should solve successfully"); + + // The optimal solution should be x = y = 1 (closest point to origin on line x+y=2) + assert!((solution.value(x) - 1.0).abs() < 1e-3); + assert!((solution.value(y) - 1.0).abs() < 1e-3); + + // Check constraint is satisfied + assert!(solution.value(x) + solution.value(y) >= 1.99); + } + + #[test] + fn test_quadratic_problem_mixed_terms() { + variables! {vars: + x; + y; + } + + // Create quadratic objective: minimize (x-1)^2 + (y-2)^2 + 2*x*y + // Expanded: x^2 - 2x + 1 + y^2 - 4y + 4 + 2xy + // = x^2 + y^2 + 2xy - 2x - 4y + 5 + let mut quadratic_obj = QuadraticAffineExpression::new(); + quadratic_obj.add_quadratic_term(x, x, 1.0); // x^2 + quadratic_obj.add_quadratic_term(y, y, 1.0); // y^2 + quadratic_obj.add_quadratic_term(x, y, 2.0); // 2*x*y + quadratic_obj.add_linear_term(x, -2.0); // -2x + quadratic_obj.add_linear_term(y, -4.0); // -4y + quadratic_obj.add_constant(5.0); // +5 + + let solution = vars + .minimise_quadratic(quadratic_obj) + .using(clarabel_quadratic) + .solve() + .expect("Should solve successfully"); + + // This is a more complex quadratic, but should still solve + // We mainly check that it doesn't crash and produces a reasonable solution + assert!(solution.value(x).is_finite()); + assert!(solution.value(y).is_finite()); + } + + // Note: More complex single-variable quadratic tests are covered in integration tests + + #[test] + fn test_quadratic_problem_with_linear_constraints() { + variables! {vars: + x; + y; + z; + } + + // Minimize x^2 + y^2 + z^2 subject to x + y + z = 3 + let mut quadratic_obj = QuadraticAffineExpression::new(); + quadratic_obj.add_quadratic_term(x, x, 1.0); + quadratic_obj.add_quadratic_term(y, y, 1.0); + quadratic_obj.add_quadratic_term(z, z, 1.0); + + let solution = vars + .minimise_quadratic(quadratic_obj) + .using(clarabel_quadratic) + .with((x + y + z).eq(3.0)) // Equality constraint + .solve() + .expect("Should solve successfully"); + + // Optimal solution should be x = y = z = 1 + assert!((solution.value(x) - 1.0).abs() < 1e-3); + assert!((solution.value(y) - 1.0).abs() < 1e-3); + assert!((solution.value(z) - 1.0).abs() < 1e-3); + } + + #[test] + fn test_quadratic_expression_evaluation() { + variables! {vars: + x; + y; + } + + // Create quadratic expression: x^2 + 2*x*y + y^2 (which is (x+y)^2) + let mut quadratic_obj = QuadraticAffineExpression::new(); + quadratic_obj.add_quadratic_term(x, x, 1.0); + quadratic_obj.add_quadratic_term(x, y, 2.0); + quadratic_obj.add_quadratic_term(y, y, 1.0); + + let solution = vars + .minimise_quadratic(quadratic_obj) + .using(clarabel_quadratic) + .with(x + y >> 2.0) // x + y >= 2 + .solve() + .expect("Should solve successfully"); + + // At the solution, (x+y)^2 should equal the sum of coefficients times the values + let x_val = solution.value(x); + let y_val = solution.value(y); + let expected_obj_value = x_val * x_val + 2.0 * x_val * y_val + y_val * y_val; + let actual_obj_value = (x_val + y_val) * (x_val + y_val); + + assert!((expected_obj_value - actual_obj_value).abs() < 1e-6); + } + } } diff --git a/tests/quadratic_solver_tests.rs b/tests/quadratic_solver_tests.rs new file mode 100644 index 0000000..c0d292d --- /dev/null +++ b/tests/quadratic_solver_tests.rs @@ -0,0 +1,253 @@ +//! Integration tests for quadratic programming functionality +//! These tests require both 'clarabel' and 'enable_quadratic' features + +#[cfg(all(feature = "clarabel", feature = "enable_quadratic"))] +mod quadratic_integration_tests { + use good_lp::{ + clarabel_quadratic, variables, QuadraticAffineExpression, ResolutionError, Solution, + SolverModel, + }; + + #[test] + fn test_constrained_quadratic_with_bounds() { + variables! { + vars: + -2.0 <= x <= 2.0; + } + + // Minimize x^2 - x + let mut objective = QuadraticAffineExpression::new(); + + objective.add_quadratic_term(x, x, 1.0); + objective.add_linear_term(x, -1.0); + + // print objective for debugging + println!("Objective: {:?}", objective); + + let solution = vars + .minimise_quadratic(objective) + .using(clarabel_quadratic) + .solve() + .expect("Bounded quadratic should solve"); + + println!("Solution: x={:.3}", solution.value(x)); + + // Unconstrained minimum is at x=0.5, which satisfies the bounds + assert!((solution.value(x) - 0.5).abs() < 1e-3); + } + + #[test] + fn test_quadratic_with_equality_constraints() { + variables! { + vars: + x; + y; + z; + } + + // Minimize x^2 + y^2 + z^2 subject to x + y + z = 6 and x - y + 2z = 4 + let mut objective = QuadraticAffineExpression::new(); + objective.add_quadratic_term(x, x, 1.0); + objective.add_quadratic_term(y, y, 1.0); + objective.add_quadratic_term(z, z, 1.0); + + let solution = vars + .minimise_quadratic(objective) + .using(clarabel_quadratic) + .with((x + y + z).eq(6.0)) // First equality constraint + .with((x - y + 2.0 * z).eq(4.0)) // Second equality constraint + .solve() + .expect("Quadratic with equality constraints should solve"); + + // Global minimum is at x=2, y=2, z=2 + assert!((solution.value(x) - 2.0).abs() < 1e-3); + assert!((solution.value(y) - 2.0).abs() < 1e-3); + assert!((solution.value(z) - 2.0).abs() < 1e-3); + + // Verify constraints are satisfied + let sum = solution.value(x) + solution.value(y) + solution.value(z); + let diff = solution.value(x) - solution.value(y) + 2.0 * solution.value(z); + + assert!( + (sum - 6.0).abs() < 1e-3, + "First constraint violated: {} ≠ 6", + sum + ); + assert!( + (diff - 4.0).abs() < 1e-3, + "Second constraint violated: {} ≠ 4", + diff + ); + + println!( + "Solution: x={:.3}, y={:.3}, z={:.3}", + solution.value(x), + solution.value(y), + solution.value(z) + ); + } + + #[test] + fn test_quadratic_with_inequality_constraints() { + variables! { + vars: + x; + y; + } + + // Minimize x^2 + y^2 subject to x + y ≥ 4 + let mut objective = QuadraticAffineExpression::new(); + objective.add_quadratic_term(x, x, 1.0); + objective.add_quadratic_term(y, y, 1.0); + + let solution = vars + .minimise_quadratic(objective) + .using(clarabel_quadratic) + .with((x + y).geq(4.0)) // inequality constraint + .solve() + .expect("Quadratic with inequality should solve"); + + // By symmetry, the minimum occurs at x = 2, y = 2 + assert!((solution.value(x) - 2.0).abs() < 1e-3); + assert!((solution.value(y) - 2.0).abs() < 1e-3); + + let sum = solution.value(x) + solution.value(y); + assert!(sum >= 4.0 - 1e-3, "Constraint violated: {} < 4", sum); + + println!( + "Solution: x={:.3}, y={:.3}", + solution.value(x), + solution.value(y) + ); + } + + #[test] + fn test_quadratic_with_mixed_constraints() { + variables! { + vars: + x; + y; + } + + // Minimize (x-1)^2 + (y-3)^2 subject to x + y = 5, x ≥ 0 + let mut objective = QuadraticAffineExpression::new(); + objective.add_quadratic_term(x, x, 1.0); + objective.add_quadratic_term(y, y, 1.0); + objective.add_linear_term(x, -2.0); + objective.add_linear_term(y, -6.0); + objective.add_constant(10.0); + + let solution = vars + .minimise_quadratic(objective) + .using(clarabel_quadratic) + .with((x + y).eq(5.0)) // equality + .with((1.0 *x).geq(0.0)) // inequality + .solve() + .expect("Quadratic with mixed constraints should solve"); + + assert!((solution.value(x) - 1.5).abs() < 1e-3); + assert!((solution.value(y) - 3.5).abs() < 1e-3); + + println!( + "Solution: x={:.3}, y={:.3}", + solution.value(x), + solution.value(y) + ); + } + + + #[test] + fn test_quadratic_maximization_problem() { + variables! { + vars: + 0.0 <= x <= 10.0; + 0.0 <= y <= 10.0; + } + + // Maximize a concave quadratic: -x^2 - y^2 + 8x + 6y + let mut objective = QuadraticAffineExpression::new(); + objective.add_quadratic_term(x, x, -1.0); // -x^2 (concave) + objective.add_quadratic_term(y, y, -1.0); // -y^2 (concave) + objective.add_linear_term(x, 8.0); // +8x + objective.add_linear_term(y, 6.0); // +6y + + let solution = vars + .maximise_quadratic(objective) + .using(clarabel_quadratic) + .solve() + .expect("Quadratic maximization should solve"); + + // Unconstrained maximum would be at x=4, y=3 (where derivatives are zero) + assert!((solution.value(x) - 4.0).abs() < 1e-3); + assert!((solution.value(y) - 3.0).abs() < 1e-3); + + // Verify bounds + assert!(solution.value(x) >= -1e-6); + assert!(solution.value(x) <= 10.0 + 1e-6); + assert!(solution.value(y) >= -1e-6); + assert!(solution.value(y) <= 10.0 + 1e-6); + } + + #[test] + fn test_infeasible_quadratic_problem() { + variables! { + vars: + x; + y; + } + + let mut objective = QuadraticAffineExpression::new(); + objective.add_quadratic_term(x, x, 1.0); + objective.add_quadratic_term(y, y, 1.0); + + // Create conflicting constraints: x + y >= 5 and x + y <= 2 + let result = vars + .minimise_quadratic(objective) + .using(clarabel_quadratic) + .with(x + y >> 5.0) // x + y >= 5 + .with(x + y << 2.0) // x + y <= 2 (conflicts with above) + .solve(); + + // Should return an infeasibility error + match result { + Err(ResolutionError::Infeasible) => { + println!("Correctly detected infeasible problem"); + } + _ => panic!("Should have detected infeasible problem"), + } + } + + #[test] + fn test_quadratic_problem_scaling() { + // Test that the solver handles different scales of quadratic problems + variables! { + vars: + x; + y; + } + + // Large-scale problem: minimize 1000000*x^2 + 1000000*y^2 + let mut objective = QuadraticAffineExpression::new(); + objective.add_quadratic_term(x, x, 1_000_000.0); + objective.add_quadratic_term(y, y, 1_000_000.0); + + let solution = vars + .minimise_quadratic(objective) + .using(clarabel_quadratic) + .with(x + y >> 1e-3) // Very small constraint + .solve() + .expect("Scaled quadratic should solve"); + + // Should still find optimal solution near x=y=5e-4 + assert!((solution.value(x) - 5e-4).abs() < 1e-6); + assert!((solution.value(y) - 5e-4).abs() < 1e-6); + } +} + +// Provide a message when quadratic features are not enabled +#[cfg(not(all(feature = "clarabel", feature = "enable_quadratic")))] +#[test] +fn quadratic_features_not_enabled() { + println!("Quadratic programming tests require both 'clarabel' and 'enable_quadratic' features"); + println!("Run with: cargo test --features 'clarabel,enable_quadratic'"); +} From 981c39e84870dcfe184cd8d1aa7487eede79d5bd Mon Sep 17 00:00:00 2001 From: Shurikal Date: Sat, 27 Sep 2025 08:18:22 +0200 Subject: [PATCH 2/8] Reduced redundant code --- src/solvers/clarabel.rs | 297 +++++++++++++++++++--------------------- 1 file changed, 144 insertions(+), 153 deletions(-) diff --git a/src/solvers/clarabel.rs b/src/solvers/clarabel.rs index da14a02..56d20c7 100644 --- a/src/solvers/clarabel.rs +++ b/src/solvers/clarabel.rs @@ -20,6 +20,94 @@ use clarabel::solver::SupportedConeT::{self, *}; use clarabel::solver::{DefaultSolution, SolverStatus}; use clarabel::solver::{DefaultSolver, IPSolver}; +impl ClarabelProblemCore { + /// Create a new ClarabelProblemCore with common setup + fn new(variables_len: usize) -> Self { + let mut settings = DefaultSettingsBuilder::default(); + settings.verbose(false).tol_feas(1e-9); + + Self { + constraints_matrix_builder: CscMatrixBuilder::new(variables_len), + constraint_values: Vec::new(), + objective: vec![0.; variables_len], + variables: variables_len, + settings, + cones: Vec::new(), + } + } + + /// Add variable bound constraints from variable definitions + fn add_variable_bounds(variables: &crate::variable::ProblemVariables, problem: &mut P) { + for (var, def) in variables.iter_variables_with_def() { + if def.is_integer { + panic!("Clarabel doesn't support integer variables") + } + if def.min != f64::NEG_INFINITY { + problem.add_constraint(var >> def.min); + } + if def.max != f64::INFINITY { + problem.add_constraint(var << def.max); + } + } + } + + /// Get direction coefficient (-1 for maximization, 1 for minimization) + fn direction_coef(direction: ObjectiveDirection) -> f64 { + if direction == ObjectiveDirection::Maximisation { + -1. + } else { + 1. + } + } + + /// Common solver status handling + fn handle_solver_status(status: SolverStatus, solution: DefaultSolution) -> Result { + match status { + e @ (SolverStatus::PrimalInfeasible | SolverStatus::AlmostPrimalInfeasible) => { + eprintln!("Clarabel error: {:?}", e); + Err(ResolutionError::Infeasible) + } + SolverStatus::Solved + | SolverStatus::AlmostSolved + | SolverStatus::AlmostDualInfeasible + | SolverStatus::DualInfeasible => Ok(ClarabelSolution { solution }), + SolverStatus::Unsolved => Err(ResolutionError::Other("Unsolved")), + SolverStatus::MaxIterations => Err(ResolutionError::Other("Max iterations reached")), + SolverStatus::MaxTime => Err(ResolutionError::Other("Time limit reached")), + SolverStatus::NumericalError => Err(ResolutionError::Other("Numerical error")), + SolverStatus::InsufficientProgress => Err(ResolutionError::Other("No progress")), + SolverStatus::CallbackTerminated => Err(ResolutionError::Other("Callback terminated")), + } + } + + /// Common constraint addition logic + fn add_constraint_impl(&mut self, constraint: Constraint) -> ConstraintReference { + self.constraints_matrix_builder + .add_row(constraint.expression.linear); + let index = self.constraint_values.len(); + self.constraint_values.push(-constraint.expression.constant); + // Cones indicate the type of constraint. We only support nonnegative and equality constraints. + // To avoid creating a new cone for each constraint, we merge them. + let next_cone = if constraint.is_equality { + ZeroConeT(1) + } else { + NonnegativeConeT(1) + }; + let prev_cone = self.cones.last_mut(); + match (prev_cone, next_cone) { + (Some(ZeroConeT(a)), ZeroConeT(b)) => *a += b, + (Some(NonnegativeConeT(a)), NonnegativeConeT(b)) => *a += b, + (_, next_cone) => self.cones.push(next_cone), + }; + ConstraintReference { index } + } +} + +/// Trait to abstract over different Clarabel problem types +trait ClarabelProblemLike { + fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference; +} + /// The [clarabel](https://oxfordcontrol.github.io/ClarabelDocs/stable/) solver, /// to be used with [UnsolvedProblem::using]. pub fn clarabel(to_solve: UnsolvedProblem) -> ClarabelProblem { @@ -28,38 +116,18 @@ pub fn clarabel(to_solve: UnsolvedProblem) -> ClarabelProblem { direction, variables, } = to_solve; - let coef = if direction == ObjectiveDirection::Maximisation { - -1. - } else { - 1. - }; - let mut objective_vector = vec![0.; variables.len()]; + let coef = ClarabelProblemCore::direction_coef(direction); + let mut core = ClarabelProblemCore::new(variables.len()); + + // Set up linear objective for (var, obj) in objective.linear_coefficients() { - objective_vector[var.index()] = obj * coef; - } - let constraints_matrix_builder = CscMatrixBuilder::new(variables.len()); - let mut settings = DefaultSettingsBuilder::default(); - settings.verbose(false).tol_feas(1e-9); - let mut p = ClarabelProblem { - objective: objective_vector, - constraints_matrix_builder, - constraint_values: Vec::new(), - variables: variables.len(), - settings, - cones: Vec::new(), - }; - // add trivial constraints embedded in the variable definitions - for (var, def) in variables.iter_variables_with_def() { - if def.is_integer { - panic!("Clarabel doesn't support integer variables") - } - if def.min != f64::NEG_INFINITY { - p.add_constraint(var >> def.min); - } - if def.max != f64::INFINITY { - p.add_constraint(var << def.max); - } + core.objective[var.index()] = obj * coef; } + + let mut p = ClarabelProblem { core }; + + // Add variable bounds + ClarabelProblemCore::add_variable_bounds(&variables, &mut p); p } @@ -72,54 +140,32 @@ pub fn clarabel_quadratic(to_solve: QuadraticUnsolvedProblem) -> ClarabelQuadrat direction, variables, } = to_solve; - let coef = if direction == ObjectiveDirection::Maximisation { - -1. - } else { - 1. - }; + let coef = ClarabelProblemCore::direction_coef(direction); + let mut core = ClarabelProblemCore::new(variables.len()); - // Linear objective vector - let mut objective_vector = vec![0.; variables.len()]; + // Set up linear objective for (var, obj_coeff) in objective.linear.coefficients { - objective_vector[var.index()] = obj_coeff * coef; + core.objective[var.index()] = obj_coeff * coef; } - // Quadratic objective matrix (P matrix) + // Set up quadratic objective matrix let mut quadratic_matrix_builder = CscQuadraticMatrixBuilder::new(variables.len()); for (pair, quad_coeff) in objective.quadratic.quadratic_coefficients { quadratic_matrix_builder.add_term(pair, quad_coeff * coef); } - let constraints_matrix_builder = CscMatrixBuilder::new(variables.len()); - let mut settings = DefaultSettingsBuilder::default(); - settings.verbose(false).tol_feas(1e-9); let mut p = ClarabelQuadraticProblem { - objective: objective_vector, + core, quadratic_matrix_builder, - constraints_matrix_builder, - constraint_values: Vec::new(), - variables: variables.len(), - settings, - cones: Vec::new(), }; - // add trivial constraints embedded in the variable definitions - for (var, def) in variables.iter_variables_with_def() { - if def.is_integer { - panic!("Clarabel doesn't support integer variables") - } - if def.min != f64::NEG_INFINITY { - p.add_constraint(var >> def.min); - } - if def.max != f64::INFINITY { - p.add_constraint(var << def.max); - } - } + // Add variable bounds + ClarabelProblemCore::add_variable_bounds(&variables, &mut p); p } -/// A clarabel model -pub struct ClarabelProblem { +/// Common fields shared between linear and quadratic Clarabel problems +struct ClarabelProblemCore { constraints_matrix_builder: CscMatrixBuilder, constraint_values: Vec, objective: Vec, @@ -128,33 +174,33 @@ pub struct ClarabelProblem { cones: Vec>, } +/// A clarabel model +pub struct ClarabelProblem { + core: ClarabelProblemCore, +} + /// A clarabel quadratic model #[cfg(feature = "enable_quadratic")] pub struct ClarabelQuadraticProblem { - constraints_matrix_builder: CscMatrixBuilder, + core: ClarabelProblemCore, quadratic_matrix_builder: CscQuadraticMatrixBuilder, - constraint_values: Vec, - objective: Vec, - variables: usize, - settings: DefaultSettingsBuilder, - cones: Vec>, } impl ClarabelProblem { /// Access the problem settings pub fn settings(&mut self) -> &mut DefaultSettingsBuilder { - &mut self.settings + &mut self.core.settings } /// Convert the problem into a clarabel solver /// panics if the problem is not valid pub fn into_solver(self) -> DefaultSolver { - let settings = self.settings.build().expect("Invalid clarabel settings"); - let quadratic_objective = &CscMatrix::zeros((self.variables, self.variables)); - let objective = &self.objective; - let constraints = &self.constraints_matrix_builder.build(); - let constraint_values = &self.constraint_values; - let cones = &self.cones; + let settings = self.core.settings.build().expect("Invalid clarabel settings"); + let quadratic_objective = &CscMatrix::zeros((self.core.variables, self.core.variables)); + let objective = &self.core.objective; + let constraints = &self.core.constraints_matrix_builder.build(); + let constraint_values = &self.core.constraint_values; + let cones = &self.core.cones; DefaultSolver::new( quadratic_objective, objective, @@ -166,6 +212,12 @@ impl ClarabelProblem { } } +impl ClarabelProblemLike for ClarabelProblem { + fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { + self.core.add_constraint_impl(constraint) + } +} + impl SolverModel for ClarabelProblem { type Solution = ClarabelSolution; type Error = ResolutionError; @@ -173,45 +225,11 @@ impl SolverModel for ClarabelProblem { fn solve(self) -> Result { let mut solver = self.into_solver(); solver.solve(); - match solver.solution.status { - e @ (SolverStatus::PrimalInfeasible | SolverStatus::AlmostPrimalInfeasible) => { - eprintln!("Clarabel error: {:?}", e); - Err(ResolutionError::Infeasible) - } - SolverStatus::Solved - | SolverStatus::AlmostSolved - | SolverStatus::AlmostDualInfeasible - | SolverStatus::DualInfeasible => Ok(ClarabelSolution { - solution: solver.solution, - }), - SolverStatus::Unsolved => Err(ResolutionError::Other("Unsolved")), - SolverStatus::MaxIterations => Err(ResolutionError::Other("Max iterations reached")), - SolverStatus::MaxTime => Err(ResolutionError::Other("Time limit reached")), - SolverStatus::NumericalError => Err(ResolutionError::Other("Numerical error")), - SolverStatus::InsufficientProgress => Err(ResolutionError::Other("No progress")), - SolverStatus::CallbackTerminated => Err(ResolutionError::Other("Callback terminated")), - } + ClarabelProblemCore::handle_solver_status(solver.solution.status, solver.solution) } fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { - self.constraints_matrix_builder - .add_row(constraint.expression.linear); - let index = self.constraint_values.len(); - self.constraint_values.push(-constraint.expression.constant); - // Cones indicate the type of constraint. We only support nonnegative and equality constraints. - // To avoid creating a new cone for each constraint, we merge them. - let next_cone = if constraint.is_equality { - ZeroConeT(1) - } else { - NonnegativeConeT(1) - }; - let prev_cone = self.cones.last_mut(); - match (prev_cone, next_cone) { - (Some(ZeroConeT(a)), ZeroConeT(b)) => *a += b, - (Some(NonnegativeConeT(a)), NonnegativeConeT(b)) => *a += b, - (_, next_cone) => self.cones.push(next_cone), - }; - ConstraintReference { index } + self.core.add_constraint_impl(constraint) } fn name() -> &'static str { @@ -322,18 +340,18 @@ fn fast_flatten_vecs(vecs: Vec>) -> Vec { impl ClarabelQuadraticProblem { /// Access the problem settings pub fn settings(&mut self) -> &mut DefaultSettingsBuilder { - &mut self.settings + &mut self.core.settings } /// Convert the problem into a clarabel solver /// panics if the problem is not valid pub fn into_solver(self) -> DefaultSolver { - let settings = self.settings.build().expect("Invalid clarabel settings"); + let settings = self.core.settings.build().expect("Invalid clarabel settings"); let quadratic_objective = &self.quadratic_matrix_builder.build(); - let objective = &self.objective; - let constraints = &self.constraints_matrix_builder.build(); - let constraint_values = &self.constraint_values; - let cones = &self.cones; + let objective = &self.core.objective; + let constraints = &self.core.constraints_matrix_builder.build(); + let constraint_values = &self.core.constraint_values; + let cones = &self.core.cones; DefaultSolver::new( quadratic_objective, objective, @@ -345,6 +363,13 @@ impl ClarabelQuadraticProblem { } } +#[cfg(feature = "enable_quadratic")] +impl ClarabelProblemLike for ClarabelQuadraticProblem { + fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { + self.core.add_constraint_impl(constraint) + } +} + #[cfg(feature = "enable_quadratic")] impl SolverModel for ClarabelQuadraticProblem { type Solution = ClarabelSolution; @@ -353,45 +378,11 @@ impl SolverModel for ClarabelQuadraticProblem { fn solve(self) -> Result { let mut solver = self.into_solver(); solver.solve(); - match solver.solution.status { - e @ (SolverStatus::PrimalInfeasible | SolverStatus::AlmostPrimalInfeasible) => { - eprintln!("Clarabel error: {:?}", e); - Err(ResolutionError::Infeasible) - } - SolverStatus::Solved - | SolverStatus::AlmostSolved - | SolverStatus::AlmostDualInfeasible - | SolverStatus::DualInfeasible => Ok(ClarabelSolution { - solution: solver.solution, - }), - SolverStatus::Unsolved => Err(ResolutionError::Other("Unsolved")), - SolverStatus::MaxIterations => Err(ResolutionError::Other("Max iterations reached")), - SolverStatus::MaxTime => Err(ResolutionError::Other("Time limit reached")), - SolverStatus::NumericalError => Err(ResolutionError::Other("Numerical error")), - SolverStatus::InsufficientProgress => Err(ResolutionError::Other("No progress")), - SolverStatus::CallbackTerminated => Err(ResolutionError::Other("Callback terminated")), - } + ClarabelProblemCore::handle_solver_status(solver.solution.status, solver.solution) } fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { - self.constraints_matrix_builder - .add_row(constraint.expression.linear); - let index = self.constraint_values.len(); - self.constraint_values.push(-constraint.expression.constant); - // Cones indicate the type of constraint. We only support nonnegative and equality constraints. - // To avoid creating a new cone for each constraint, we merge them. - let next_cone = if constraint.is_equality { - ZeroConeT(1) - } else { - NonnegativeConeT(1) - }; - let prev_cone = self.cones.last_mut(); - match (prev_cone, next_cone) { - (Some(ZeroConeT(a)), ZeroConeT(b)) => *a += b, - (Some(NonnegativeConeT(a)), NonnegativeConeT(b)) => *a += b, - (_, next_cone) => self.cones.push(next_cone), - }; - ConstraintReference { index } + self.core.add_constraint_impl(constraint) } fn name() -> &'static str { From d1187fcf14631629cf99c35541ef8b56bffe34ae Mon Sep 17 00:00:00 2001 From: Shurikal Date: Sat, 27 Sep 2025 08:22:57 +0200 Subject: [PATCH 3/8] Fixed wrong comment --- src/solvers/clarabel.rs | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/solvers/clarabel.rs b/src/solvers/clarabel.rs index 56d20c7..82db4b9 100644 --- a/src/solvers/clarabel.rs +++ b/src/solvers/clarabel.rs @@ -420,10 +420,9 @@ impl CscQuadraticMatrixBuilder { self.rowval[j].push(i); self.nzval[j].push(2.0 * coeff); } else { - // Off-diagonal term: the coefficient coeff represents the full cross term coeff * x_i * x_j - // In matrix form, this becomes coeff/2 * (P_ij + P_ji) * x_i * x_j - // Since P is symmetric (P_ij = P_ji), we get coeff * P_ij * x_i * x_j - // But Clarabel uses 1/2 * x^T * P * x, so we need P_ij = 2 * coeff + // Off-diagonal term: add the coefficient to both (i,j) and (j,i) positions + // to create a symmetric matrix. Each entry gets the original coefficient + // since Clarabel multiplies the entire quadratic form by 1/2. self.rowval[j].push(i); self.nzval[j].push(coeff); From 422111e153925bca55440cd6068360b764e497ef Mon Sep 17 00:00:00 2001 From: Shurikal Date: Sat, 27 Sep 2025 08:27:14 +0200 Subject: [PATCH 4/8] Applied cargo fmt --- src/lib.rs | 15 ++-- src/quadratic_expression.rs | 94 +++++++++++++++---------- src/quadratic_problem.rs | 28 +++++--- src/solvers/clarabel.rs | 119 ++++++++++++++++++-------------- tests/quadratic_solver_tests.rs | 3 +- 5 files changed, 153 insertions(+), 106 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 2a1bc9d..5f772f1 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -70,9 +70,11 @@ pub use cardinality_constraint_solver_trait::CardinalityConstraintSolver; pub use constraint::Constraint; pub use expression::Expression; #[cfg(feature = "enable_quadratic")] -pub use quadratic_expression::{QuadraticAffineExpression, QuadraticConstraint, QuadraticExpression, VariablePair}; +pub use quadratic_expression::{ + QuadraticAffineExpression, QuadraticConstraint, QuadraticExpression, VariablePair, +}; #[cfg(feature = "enable_quadratic")] -pub use quadratic_problem::{QuadraticUnsolvedProblem, IntoQuadraticExpression}; +pub use quadratic_problem::{IntoQuadraticExpression, QuadraticUnsolvedProblem}; #[cfg(not(any( feature = "coin_cbc", feature = "microlp", @@ -86,7 +88,10 @@ pub use solvers::clarabel::clarabel as default_solver; #[cfg_attr(docsrs, doc(cfg(feature = "clarabel")))] #[cfg(feature = "clarabel")] pub use solvers::clarabel::clarabel; -#[cfg_attr(docsrs, doc(cfg(all(feature = "clarabel", feature = "enable_quadratic"))))] +#[cfg_attr( + docsrs, + doc(cfg(all(feature = "clarabel", feature = "enable_quadratic"))) +)] #[cfg(all(feature = "clarabel", feature = "enable_quadratic"))] pub use solvers::clarabel::clarabel_quadratic; #[cfg_attr(docsrs, doc(cfg(feature = "coin_cbc")))] @@ -195,11 +200,11 @@ pub mod variable; mod affine_expression_trait; mod cardinality_constraint_solver_trait; pub mod constraint; -pub mod solvers; -mod variables_macro; /// Quadratic expression types and operations #[cfg(feature = "enable_quadratic")] pub mod quadratic_expression; /// Quadratic problem definition and solving #[cfg(feature = "enable_quadratic")] pub mod quadratic_problem; +pub mod solvers; +mod variables_macro; diff --git a/src/quadratic_expression.rs b/src/quadratic_expression.rs index ca983ef..40a6037 100644 --- a/src/quadratic_expression.rs +++ b/src/quadratic_expression.rs @@ -27,13 +27,19 @@ impl VariablePair { if var1.index() <= var2.index() { VariablePair { var1, var2 } } else { - VariablePair { var1: var2, var2: var1 } + VariablePair { + var1: var2, + var2: var1, + } } } /// Create a variable pair representing a square term (var^2) pub fn square(var: Variable) -> Self { - VariablePair { var1: var, var2: var } + VariablePair { + var1: var, + var2: var, + } } } @@ -66,24 +72,32 @@ impl QuadraticExpression { /// Get the coefficient for a specific variable pair pub fn get_quadratic_coefficient(&self, var1: Variable, var2: Variable) -> f64 { let pair = VariablePair::new(var1, var2); - self.quadratic_coefficients.get(&pair).copied().unwrap_or(0.0) + self.quadratic_coefficients + .get(&pair) + .copied() + .unwrap_or(0.0) } /// Check if the expression is empty (all coefficients are zero) pub fn is_empty(&self) -> bool { - self.quadratic_coefficients.is_empty() || - self.quadratic_coefficients.values().all(|&c| c.abs() < f64::EPSILON) + self.quadratic_coefficients.is_empty() + || self + .quadratic_coefficients + .values() + .all(|&c| c.abs() < f64::EPSILON) } /// Iterate over all non-zero quadratic terms pub fn iter(&self) -> impl Iterator + '_ { - self.quadratic_coefficients.iter().filter_map(|(&pair, &coeff)| { - if coeff.abs() >= f64::EPSILON { - Some((pair, coeff)) - } else { - None - } - }) + self.quadratic_coefficients + .iter() + .filter_map(|(&pair, &coeff)| { + if coeff.abs() >= f64::EPSILON { + Some((pair, coeff)) + } else { + None + } + }) } /// Evaluate the quadratic expression given variable values @@ -138,7 +152,10 @@ impl QuadraticAffineExpression { QuadraticAffineExpression { quadratic: QuadraticExpression::with_capacity(quadratic_capacity), linear: LinearExpression { - coefficients: HashMap::with_capacity_and_hasher(linear_capacity, Default::default()), + coefficients: HashMap::with_capacity_and_hasher( + linear_capacity, + Default::default(), + ), }, constant: 0.0, } @@ -150,7 +167,9 @@ impl QuadraticAffineExpression { let linear_coeffs = expr.linear_coefficients().into_iter().collect(); QuadraticAffineExpression { quadratic: QuadraticExpression::new(), - linear: LinearExpression { coefficients: linear_coeffs }, + linear: LinearExpression { + coefficients: linear_coeffs, + }, constant, } } @@ -198,8 +217,7 @@ impl QuadraticAffineExpression { /// Evaluate the complete expression pub fn eval_with(&self, values: &S) -> f64 { - self.quadratic.eval_with(values) + - self.linear_part().eval_with(values) + self.quadratic.eval_with(values) + self.linear_part().eval_with(values) } /// Get the coefficient of a quadratic term @@ -305,7 +323,9 @@ impl From for QuadraticAffineExpression { linear_coeffs.insert(var, 1.0); QuadraticAffineExpression { quadratic: QuadraticExpression::new(), - linear: LinearExpression { coefficients: linear_coeffs }, + linear: LinearExpression { + coefficients: linear_coeffs, + }, constant: 0.0, } } @@ -532,17 +552,17 @@ impl Mul for Variable { fn mul(self, rhs: Expression) -> Self::Output { let mut result = QuadraticAffineExpression::new(); let constant = rhs.constant; - + // Linear term (constant * variable) if constant.abs() >= f64::EPSILON { result.add_linear_term(self, constant); } - + // Quadratic terms (variable * other_variables) for (var, coeff) in rhs.linear.coefficients { result.add_quadratic_term(self, var, coeff); } - + result } } @@ -561,30 +581,30 @@ impl Mul for Expression { fn mul(self, rhs: Expression) -> Self::Output { let mut result = QuadraticAffineExpression::new(); - + let c1 = self.constant; let c2 = rhs.constant; - + // Constant term (c1 * c2) result.constant = c1 * c2; - + // Linear terms from self * constant of rhs for (var, coeff) in self.linear.coefficients.iter() { result.add_linear_term(*var, *coeff * c2); } - + // Linear terms from constant of self * rhs for (var, coeff) in rhs.linear.coefficients.iter() { result.add_linear_term(*var, c1 * *coeff); } - + // Quadratic terms (variable from self * variable from rhs) for (var1, coeff1) in self.linear.coefficients.iter() { for (var2, coeff2) in rhs.linear.coefficients.iter() { result.add_quadratic_term(*var1, *var2, *coeff1 * *coeff2); } } - + result } } @@ -626,11 +646,11 @@ impl FormatWithVars for QuadraticExpression { } else { write!(f, " + ")?; } - + if (coeff - 1.).abs() > f64::EPSILON { write!(f, "{} ", coeff)?; } - + if pair.var1 == pair.var2 { // Square term: x^2 variable_format(f, pair.var1)?; @@ -642,7 +662,7 @@ impl FormatWithVars for QuadraticExpression { variable_format(f, pair.var2)?; } } - + if first { write!(f, "0")?; } @@ -664,19 +684,19 @@ impl FormatWithVars for QuadraticAffineExpression { let has_quadratic = !self.quadratic.is_empty(); let has_linear = !self.linear.coefficients.is_empty(); let has_constant = self.constant.abs() >= f64::EPSILON; - + if !has_quadratic && !has_linear && !has_constant { return write!(f, "0"); } - + let mut first = true; - + // Write quadratic terms if has_quadratic { self.quadratic.format_with(f, &mut variable_format)?; first = false; } - + // Write linear terms for (&var, &coeff) in &self.linear.coefficients { if coeff.abs() >= f64::EPSILON { @@ -684,14 +704,14 @@ impl FormatWithVars for QuadraticAffineExpression { write!(f, " + ")?; } first = false; - + if (coeff - 1.).abs() > f64::EPSILON { write!(f, "{} ", coeff)?; } variable_format(f, var)?; } } - + // Write constant term if has_constant { if !first { @@ -699,7 +719,7 @@ impl FormatWithVars for QuadraticAffineExpression { } write!(f, "{}", self.constant)?; } - + Ok(()) } } @@ -708,4 +728,4 @@ impl Debug for QuadraticAffineExpression { fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { self.format_debug(f) } -} \ No newline at end of file +} diff --git a/src/quadratic_problem.rs b/src/quadratic_problem.rs index be25a41..f6841bf 100644 --- a/src/quadratic_problem.rs +++ b/src/quadratic_problem.rs @@ -1,6 +1,6 @@ use crate::quadratic_expression::QuadraticAffineExpression; -use crate::variable::ProblemVariables; use crate::solvers::ObjectiveDirection; +use crate::variable::ProblemVariables; /// A problem with a quadratic objective function pub struct QuadraticUnsolvedProblem { @@ -46,12 +46,16 @@ impl ProblemVariables { let objective_expr = objective.into_quadratic_expression(); // Validate that objective variables are within bounds for (var, _) in objective_expr.linear.coefficients.iter() { - assert!(var.index() < self.len(), - "Variable in objective function is not part of this problem"); + assert!( + var.index() < self.len(), + "Variable in objective function is not part of this problem" + ); } for (pair, _) in objective_expr.quadratic.quadratic_coefficients.iter() { - assert!(pair.var1.index() < self.len() && pair.var2.index() < self.len(), - "Variable in quadratic objective function is not part of this problem"); + assert!( + pair.var1.index() < self.len() && pair.var2.index() < self.len(), + "Variable in quadratic objective function is not part of this problem" + ); } QuadraticUnsolvedProblem { objective: objective_expr, @@ -69,12 +73,16 @@ impl ProblemVariables { let objective_expr = objective.into_quadratic_expression(); // Validate that objective variables are within bounds for (var, _) in objective_expr.linear.coefficients.iter() { - assert!(var.index() < self.len(), - "Variable in objective function is not part of this problem"); + assert!( + var.index() < self.len(), + "Variable in objective function is not part of this problem" + ); } for (pair, _) in objective_expr.quadratic.quadratic_coefficients.iter() { - assert!(pair.var1.index() < self.len() && pair.var2.index() < self.len(), - "Variable in quadratic objective function is not part of this problem"); + assert!( + pair.var1.index() < self.len() && pair.var2.index() < self.len(), + "Variable in quadratic objective function is not part of this problem" + ); } QuadraticUnsolvedProblem { objective: objective_expr, @@ -93,4 +101,4 @@ impl QuadraticUnsolvedProblem { { solver_factory(self) } -} \ No newline at end of file +} diff --git a/src/solvers/clarabel.rs b/src/solvers/clarabel.rs index 82db4b9..1e8e8f5 100644 --- a/src/solvers/clarabel.rs +++ b/src/solvers/clarabel.rs @@ -2,11 +2,11 @@ use crate::affine_expression_trait::IntoAffineExpression; use crate::expression::LinearExpression; -use crate::variable::UnsolvedProblem; -#[cfg(feature = "enable_quadratic")] -use crate::quadratic_problem::QuadraticUnsolvedProblem; #[cfg(feature = "enable_quadratic")] use crate::quadratic_expression::VariablePair; +#[cfg(feature = "enable_quadratic")] +use crate::quadratic_problem::QuadraticUnsolvedProblem; +use crate::variable::UnsolvedProblem; use crate::{ constraint::ConstraintReference, solvers::{ObjectiveDirection, ResolutionError, Solution, SolverModel}, @@ -25,7 +25,7 @@ impl ClarabelProblemCore { fn new(variables_len: usize) -> Self { let mut settings = DefaultSettingsBuilder::default(); settings.verbose(false).tol_feas(1e-9); - + Self { constraints_matrix_builder: CscMatrixBuilder::new(variables_len), constraint_values: Vec::new(), @@ -37,7 +37,10 @@ impl ClarabelProblemCore { } /// Add variable bound constraints from variable definitions - fn add_variable_bounds(variables: &crate::variable::ProblemVariables, problem: &mut P) { + fn add_variable_bounds( + variables: &crate::variable::ProblemVariables, + problem: &mut P, + ) { for (var, def) in variables.iter_variables_with_def() { if def.is_integer { panic!("Clarabel doesn't support integer variables") @@ -61,7 +64,10 @@ impl ClarabelProblemCore { } /// Common solver status handling - fn handle_solver_status(status: SolverStatus, solution: DefaultSolution) -> Result { + fn handle_solver_status( + status: SolverStatus, + solution: DefaultSolution, + ) -> Result { match status { e @ (SolverStatus::PrimalInfeasible | SolverStatus::AlmostPrimalInfeasible) => { eprintln!("Clarabel error: {:?}", e); @@ -118,14 +124,14 @@ pub fn clarabel(to_solve: UnsolvedProblem) -> ClarabelProblem { } = to_solve; let coef = ClarabelProblemCore::direction_coef(direction); let mut core = ClarabelProblemCore::new(variables.len()); - + // Set up linear objective for (var, obj) in objective.linear_coefficients() { core.objective[var.index()] = obj * coef; } - + let mut p = ClarabelProblem { core }; - + // Add variable bounds ClarabelProblemCore::add_variable_bounds(&variables, &mut p); p @@ -142,23 +148,23 @@ pub fn clarabel_quadratic(to_solve: QuadraticUnsolvedProblem) -> ClarabelQuadrat } = to_solve; let coef = ClarabelProblemCore::direction_coef(direction); let mut core = ClarabelProblemCore::new(variables.len()); - + // Set up linear objective for (var, obj_coeff) in objective.linear.coefficients { core.objective[var.index()] = obj_coeff * coef; } - + // Set up quadratic objective matrix let mut quadratic_matrix_builder = CscQuadraticMatrixBuilder::new(variables.len()); for (pair, quad_coeff) in objective.quadratic.quadratic_coefficients { quadratic_matrix_builder.add_term(pair, quad_coeff * coef); } - + let mut p = ClarabelQuadraticProblem { core, quadratic_matrix_builder, }; - + // Add variable bounds ClarabelProblemCore::add_variable_bounds(&variables, &mut p); p @@ -195,7 +201,11 @@ impl ClarabelProblem { /// Convert the problem into a clarabel solver /// panics if the problem is not valid pub fn into_solver(self) -> DefaultSolver { - let settings = self.core.settings.build().expect("Invalid clarabel settings"); + let settings = self + .core + .settings + .build() + .expect("Invalid clarabel settings"); let quadratic_objective = &CscMatrix::zeros((self.core.variables, self.core.variables)); let objective = &self.core.objective; let constraints = &self.core.constraints_matrix_builder.build(); @@ -346,7 +356,11 @@ impl ClarabelQuadraticProblem { /// Convert the problem into a clarabel solver /// panics if the problem is not valid pub fn into_solver(self) -> DefaultSolver { - let settings = self.core.settings.build().expect("Invalid clarabel settings"); + let settings = self + .core + .settings + .build() + .expect("Invalid clarabel settings"); let quadratic_objective = &self.quadratic_matrix_builder.build(); let objective = &self.core.objective; let constraints = &self.core.constraints_matrix_builder.build(); @@ -409,11 +423,11 @@ impl CscQuadraticMatrixBuilder { n_vars, } } - + fn add_term(&mut self, pair: VariablePair, coeff: f64) { let i = pair.var1.index(); let j = pair.var2.index(); - + if i == j { // Diagonal term: multiply by 2 because Clarabel minimizes 1/2 * x^T * P * x // So we need P_ii = 2 * coeff to get coeff * x_i^2 in the objective @@ -425,25 +439,26 @@ impl CscQuadraticMatrixBuilder { // since Clarabel multiplies the entire quadratic form by 1/2. self.rowval[j].push(i); self.nzval[j].push(coeff); - + self.rowval[i].push(j); self.nzval[i].push(coeff); } } - + fn build(mut self) -> clarabel::algebra::CscMatrix { // Sort each column by row index to maintain proper CSC format for col in 0..self.n_vars { - let mut pairs: Vec<(usize, f64)> = self.rowval[col].iter() + let mut pairs: Vec<(usize, f64)> = self.rowval[col] + .iter() .zip(self.nzval[col].iter()) .map(|(&row, &val)| (row, val)) .collect(); pairs.sort_by_key(|&(row, _)| row); - + self.rowval[col] = pairs.iter().map(|&(row, _)| row).collect(); self.nzval[col] = pairs.iter().map(|&(_, val)| val).collect(); } - + let mut colptr = Vec::with_capacity(self.n_vars + 1); colptr.push(0); for col in &self.rowval { @@ -501,19 +516,19 @@ mod tests { x; y; } - + let mut builder = CscQuadraticMatrixBuilder::new(2); - + // Add diagonal terms: x^2 with coefficient 2.0, y^2 with coefficient 3.0 builder.add_term(VariablePair::new(x, x), 2.0); builder.add_term(VariablePair::new(y, y), 3.0); - + let matrix = builder.build(); - + // Check diagonal entries (scaled by 2 for Clarabel's 1/2 * x^T * P * x formulation) assert_eq!(matrix.get_entry((0, 0)), Some(4.0)); // 2.0 * 2 = 4.0 assert_eq!(matrix.get_entry((1, 1)), Some(6.0)); // 3.0 * 2 = 6.0 - + // Check off-diagonal entries are zero assert_eq!(matrix.get_entry((0, 1)), None); assert_eq!(matrix.get_entry((1, 0)), None); @@ -525,18 +540,18 @@ mod tests { x; y; } - + let mut builder = CscQuadraticMatrixBuilder::new(3); - + // Add off-diagonal term: 2*x*y (coefficient 2.0) builder.add_term(VariablePair::new(x, y), 2.0); - + let matrix = builder.build(); - + // Off-diagonal terms: coefficient is stored directly (no factor of 2 for off-diagonal) assert_eq!(matrix.get_entry((0, 1)), Some(2.0)); // (x,y) entry assert_eq!(matrix.get_entry((1, 0)), Some(2.0)); // (y,x) entry (symmetric) - + // Other entries should be empty assert_eq!(matrix.get_entry((0, 0)), None); assert_eq!(matrix.get_entry((1, 1)), None); @@ -549,23 +564,23 @@ mod tests { x; y; } - + let mut builder = CscQuadraticMatrixBuilder::new(2); - + // Add mixed terms: x^2 + 2*x*y + y^2 (like (x+y)^2) builder.add_term(VariablePair::new(x, x), 1.0); builder.add_term(VariablePair::new(x, y), 2.0); builder.add_term(VariablePair::new(y, y), 1.0); - + let matrix = builder.build(); - + // Check the resulting matrix (accounting for Clarabel's 1/2 factor): // For (x+y)^2 = x^2 + 2*x*y + y^2, we expect: // [2.0 2.0] (diagonal terms scaled by 2, off-diagonal terms as-is) // [2.0 2.0] assert_eq!(matrix.get_entry((0, 0)), Some(2.0)); // x^2 (scaled by 2) assert_eq!(matrix.get_entry((0, 1)), Some(2.0)); // x*y coefficient - assert_eq!(matrix.get_entry((1, 0)), Some(2.0)); // y*x coefficient + assert_eq!(matrix.get_entry((1, 0)), Some(2.0)); // y*x coefficient assert_eq!(matrix.get_entry((1, 1)), Some(2.0)); // y^2 (scaled by 2) } @@ -578,8 +593,8 @@ mod tests { // Create a simple quadratic objective: minimize x^2 + y^2 let mut quadratic_obj = QuadraticAffineExpression::new(); - quadratic_obj.add_quadratic_term(x, x, 1.0); // x^2 - quadratic_obj.add_quadratic_term(y, y, 1.0); // y^2 + quadratic_obj.add_quadratic_term(x, x, 1.0); // x^2 + quadratic_obj.add_quadratic_term(y, y, 1.0); // y^2 // Create the problem let problem = vars @@ -588,7 +603,7 @@ mod tests { // The problem should be solvable (unconstrained minimum at (0,0)) let solution = problem.solve().expect("Should solve successfully"); - + // Check that solution is near the origin assert!((solution.value(x)).abs() < 1e-6); assert!((solution.value(y)).abs() < 1e-6); @@ -603,8 +618,8 @@ mod tests { // Create quadratic objective: minimize x^2 + y^2 let mut quadratic_obj = QuadraticAffineExpression::new(); - quadratic_obj.add_quadratic_term(x, x, 1.0); // x^2 - quadratic_obj.add_quadratic_term(y, y, 1.0); // y^2 + quadratic_obj.add_quadratic_term(x, x, 1.0); // x^2 + quadratic_obj.add_quadratic_term(y, y, 1.0); // y^2 // Solve with constraint x + y >= 2 let solution = vars @@ -617,7 +632,7 @@ mod tests { // The optimal solution should be x = y = 1 (closest point to origin on line x+y=2) assert!((solution.value(x) - 1.0).abs() < 1e-3); assert!((solution.value(y) - 1.0).abs() < 1e-3); - + // Check constraint is satisfied assert!(solution.value(x) + solution.value(y) >= 1.99); } @@ -633,12 +648,12 @@ mod tests { // Expanded: x^2 - 2x + 1 + y^2 - 4y + 4 + 2xy // = x^2 + y^2 + 2xy - 2x - 4y + 5 let mut quadratic_obj = QuadraticAffineExpression::new(); - quadratic_obj.add_quadratic_term(x, x, 1.0); // x^2 - quadratic_obj.add_quadratic_term(y, y, 1.0); // y^2 - quadratic_obj.add_quadratic_term(x, y, 2.0); // 2*x*y - quadratic_obj.add_linear_term(x, -2.0); // -2x - quadratic_obj.add_linear_term(y, -4.0); // -4y - quadratic_obj.add_constant(5.0); // +5 + quadratic_obj.add_quadratic_term(x, x, 1.0); // x^2 + quadratic_obj.add_quadratic_term(y, y, 1.0); // y^2 + quadratic_obj.add_quadratic_term(x, y, 2.0); // 2*x*y + quadratic_obj.add_linear_term(x, -2.0); // -2x + quadratic_obj.add_linear_term(y, -4.0); // -4y + quadratic_obj.add_constant(5.0); // +5 let solution = vars .minimise_quadratic(quadratic_obj) @@ -671,7 +686,7 @@ mod tests { let solution = vars .minimise_quadratic(quadratic_obj) .using(clarabel_quadratic) - .with((x + y + z).eq(3.0)) // Equality constraint + .with((x + y + z).eq(3.0)) // Equality constraint .solve() .expect("Should solve successfully"); @@ -697,7 +712,7 @@ mod tests { let solution = vars .minimise_quadratic(quadratic_obj) .using(clarabel_quadratic) - .with(x + y >> 2.0) // x + y >= 2 + .with(x + y >> 2.0) // x + y >= 2 .solve() .expect("Should solve successfully"); @@ -706,7 +721,7 @@ mod tests { let y_val = solution.value(y); let expected_obj_value = x_val * x_val + 2.0 * x_val * y_val + y_val * y_val; let actual_obj_value = (x_val + y_val) * (x_val + y_val); - + assert!((expected_obj_value - actual_obj_value).abs() < 1e-6); } } diff --git a/tests/quadratic_solver_tests.rs b/tests/quadratic_solver_tests.rs index c0d292d..b7703b1 100644 --- a/tests/quadratic_solver_tests.rs +++ b/tests/quadratic_solver_tests.rs @@ -141,7 +141,7 @@ mod quadratic_integration_tests { .minimise_quadratic(objective) .using(clarabel_quadratic) .with((x + y).eq(5.0)) // equality - .with((1.0 *x).geq(0.0)) // inequality + .with((1.0 * x).geq(0.0)) // inequality .solve() .expect("Quadratic with mixed constraints should solve"); @@ -155,7 +155,6 @@ mod quadratic_integration_tests { ); } - #[test] fn test_quadratic_maximization_problem() { variables! { From 75d662d5ba8589ae700f9ace3bf641ca28fb59e1 Mon Sep 17 00:00:00 2001 From: Shurikal Date: Sat, 27 Sep 2025 08:49:29 +0200 Subject: [PATCH 5/8] Tried to minimize changes --- src/solvers/clarabel.rs | 318 ++++++++++++++++++++-------------------- 1 file changed, 161 insertions(+), 157 deletions(-) diff --git a/src/solvers/clarabel.rs b/src/solvers/clarabel.rs index 1e8e8f5..8a2e59d 100644 --- a/src/solvers/clarabel.rs +++ b/src/solvers/clarabel.rs @@ -2,11 +2,11 @@ use crate::affine_expression_trait::IntoAffineExpression; use crate::expression::LinearExpression; -#[cfg(feature = "enable_quadratic")] -use crate::quadratic_expression::VariablePair; +use crate::variable::UnsolvedProblem; #[cfg(feature = "enable_quadratic")] use crate::quadratic_problem::QuadraticUnsolvedProblem; -use crate::variable::UnsolvedProblem; +#[cfg(feature = "enable_quadratic")] +use crate::quadratic_expression::VariablePair; use crate::{ constraint::ConstraintReference, solvers::{ObjectiveDirection, ResolutionError, Solution, SolverModel}, @@ -20,100 +20,6 @@ use clarabel::solver::SupportedConeT::{self, *}; use clarabel::solver::{DefaultSolution, SolverStatus}; use clarabel::solver::{DefaultSolver, IPSolver}; -impl ClarabelProblemCore { - /// Create a new ClarabelProblemCore with common setup - fn new(variables_len: usize) -> Self { - let mut settings = DefaultSettingsBuilder::default(); - settings.verbose(false).tol_feas(1e-9); - - Self { - constraints_matrix_builder: CscMatrixBuilder::new(variables_len), - constraint_values: Vec::new(), - objective: vec![0.; variables_len], - variables: variables_len, - settings, - cones: Vec::new(), - } - } - - /// Add variable bound constraints from variable definitions - fn add_variable_bounds( - variables: &crate::variable::ProblemVariables, - problem: &mut P, - ) { - for (var, def) in variables.iter_variables_with_def() { - if def.is_integer { - panic!("Clarabel doesn't support integer variables") - } - if def.min != f64::NEG_INFINITY { - problem.add_constraint(var >> def.min); - } - if def.max != f64::INFINITY { - problem.add_constraint(var << def.max); - } - } - } - - /// Get direction coefficient (-1 for maximization, 1 for minimization) - fn direction_coef(direction: ObjectiveDirection) -> f64 { - if direction == ObjectiveDirection::Maximisation { - -1. - } else { - 1. - } - } - - /// Common solver status handling - fn handle_solver_status( - status: SolverStatus, - solution: DefaultSolution, - ) -> Result { - match status { - e @ (SolverStatus::PrimalInfeasible | SolverStatus::AlmostPrimalInfeasible) => { - eprintln!("Clarabel error: {:?}", e); - Err(ResolutionError::Infeasible) - } - SolverStatus::Solved - | SolverStatus::AlmostSolved - | SolverStatus::AlmostDualInfeasible - | SolverStatus::DualInfeasible => Ok(ClarabelSolution { solution }), - SolverStatus::Unsolved => Err(ResolutionError::Other("Unsolved")), - SolverStatus::MaxIterations => Err(ResolutionError::Other("Max iterations reached")), - SolverStatus::MaxTime => Err(ResolutionError::Other("Time limit reached")), - SolverStatus::NumericalError => Err(ResolutionError::Other("Numerical error")), - SolverStatus::InsufficientProgress => Err(ResolutionError::Other("No progress")), - SolverStatus::CallbackTerminated => Err(ResolutionError::Other("Callback terminated")), - } - } - - /// Common constraint addition logic - fn add_constraint_impl(&mut self, constraint: Constraint) -> ConstraintReference { - self.constraints_matrix_builder - .add_row(constraint.expression.linear); - let index = self.constraint_values.len(); - self.constraint_values.push(-constraint.expression.constant); - // Cones indicate the type of constraint. We only support nonnegative and equality constraints. - // To avoid creating a new cone for each constraint, we merge them. - let next_cone = if constraint.is_equality { - ZeroConeT(1) - } else { - NonnegativeConeT(1) - }; - let prev_cone = self.cones.last_mut(); - match (prev_cone, next_cone) { - (Some(ZeroConeT(a)), ZeroConeT(b)) => *a += b, - (Some(NonnegativeConeT(a)), NonnegativeConeT(b)) => *a += b, - (_, next_cone) => self.cones.push(next_cone), - }; - ConstraintReference { index } - } -} - -/// Trait to abstract over different Clarabel problem types -trait ClarabelProblemLike { - fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference; -} - /// The [clarabel](https://oxfordcontrol.github.io/ClarabelDocs/stable/) solver, /// to be used with [UnsolvedProblem::using]. pub fn clarabel(to_solve: UnsolvedProblem) -> ClarabelProblem { @@ -122,18 +28,38 @@ pub fn clarabel(to_solve: UnsolvedProblem) -> ClarabelProblem { direction, variables, } = to_solve; - let coef = ClarabelProblemCore::direction_coef(direction); - let mut core = ClarabelProblemCore::new(variables.len()); - - // Set up linear objective + let coef = if direction == ObjectiveDirection::Maximisation { + -1. + } else { + 1. + }; + let mut objective_vector = vec![0.; variables.len()]; for (var, obj) in objective.linear_coefficients() { - core.objective[var.index()] = obj * coef; + objective_vector[var.index()] = obj * coef; + } + let constraints_matrix_builder = CscMatrixBuilder::new(variables.len()); + let mut settings = DefaultSettingsBuilder::default(); + settings.verbose(false).tol_feas(1e-9); + let mut p = ClarabelProblem { + objective: objective_vector, + constraints_matrix_builder, + constraint_values: Vec::new(), + variables: variables.len(), + settings, + cones: Vec::new(), + }; + // add trivial constraints embedded in the variable definitions + for (var, def) in variables.iter_variables_with_def() { + if def.is_integer { + panic!("Clarabel doesn't support integer variables") + } + if def.min != f64::NEG_INFINITY { + p.add_constraint(var >> def.min); + } + if def.max != f64::INFINITY { + p.add_constraint(var << def.max); + } } - - let mut p = ClarabelProblem { core }; - - // Add variable bounds - ClarabelProblemCore::add_variable_bounds(&variables, &mut p); p } @@ -146,32 +72,53 @@ pub fn clarabel_quadratic(to_solve: QuadraticUnsolvedProblem) -> ClarabelQuadrat direction, variables, } = to_solve; - let coef = ClarabelProblemCore::direction_coef(direction); - let mut core = ClarabelProblemCore::new(variables.len()); + let coef = if direction == ObjectiveDirection::Maximisation { + -1. + } else { + 1. + }; - // Set up linear objective + let mut objective_vector = vec![0.; variables.len()]; for (var, obj_coeff) in objective.linear.coefficients { - core.objective[var.index()] = obj_coeff * coef; + objective_vector[var.index()] = obj_coeff * coef; } - // Set up quadratic objective matrix let mut quadratic_matrix_builder = CscQuadraticMatrixBuilder::new(variables.len()); for (pair, quad_coeff) in objective.quadratic.quadratic_coefficients { quadratic_matrix_builder.add_term(pair, quad_coeff * coef); } + let constraints_matrix_builder = CscMatrixBuilder::new(variables.len()); + let mut settings = DefaultSettingsBuilder::default(); + settings.verbose(false).tol_feas(1e-9); + let mut p = ClarabelQuadraticProblem { - core, + objective: objective_vector, + constraints_matrix_builder, + constraint_values: Vec::new(), + variables: variables.len(), + settings, + cones: Vec::new(), quadratic_matrix_builder, }; - // Add variable bounds - ClarabelProblemCore::add_variable_bounds(&variables, &mut p); + // add trivial constraints embedded in the variable definitions + for (var, def) in variables.iter_variables_with_def() { + if def.is_integer { + panic!("Clarabel doesn't support integer variables") + } + if def.min != f64::NEG_INFINITY { + p.add_constraint(var >> def.min); + } + if def.max != f64::INFINITY { + p.add_constraint(var << def.max); + } + } p } -/// Common fields shared between linear and quadratic Clarabel problems -struct ClarabelProblemCore { +/// A clarabel model +pub struct ClarabelProblem { constraints_matrix_builder: CscMatrixBuilder, constraint_values: Vec, objective: Vec, @@ -180,37 +127,33 @@ struct ClarabelProblemCore { cones: Vec>, } -/// A clarabel model -pub struct ClarabelProblem { - core: ClarabelProblemCore, -} - /// A clarabel quadratic model #[cfg(feature = "enable_quadratic")] pub struct ClarabelQuadraticProblem { - core: ClarabelProblemCore, + constraints_matrix_builder: CscMatrixBuilder, + constraint_values: Vec, + objective: Vec, + variables: usize, + settings: DefaultSettingsBuilder, + cones: Vec>, quadratic_matrix_builder: CscQuadraticMatrixBuilder, } impl ClarabelProblem { /// Access the problem settings pub fn settings(&mut self) -> &mut DefaultSettingsBuilder { - &mut self.core.settings + &mut self.settings } /// Convert the problem into a clarabel solver /// panics if the problem is not valid pub fn into_solver(self) -> DefaultSolver { - let settings = self - .core - .settings - .build() - .expect("Invalid clarabel settings"); - let quadratic_objective = &CscMatrix::zeros((self.core.variables, self.core.variables)); - let objective = &self.core.objective; - let constraints = &self.core.constraints_matrix_builder.build(); - let constraint_values = &self.core.constraint_values; - let cones = &self.core.cones; + let settings = self.settings.build().expect("Invalid clarabel settings"); + let quadratic_objective = &CscMatrix::zeros((self.variables, self.variables)); + let objective = &self.objective; + let constraints = &self.constraints_matrix_builder.build(); + let constraint_values = &self.constraint_values; + let cones = &self.cones; DefaultSolver::new( quadratic_objective, objective, @@ -220,11 +163,27 @@ impl ClarabelProblem { settings, ).expect("Invalid clarabel problem. This is likely a bug in good_lp. Problems should always have coherent dimensions.") } -} -impl ClarabelProblemLike for ClarabelProblem { - fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { - self.core.add_constraint_impl(constraint) + /// Add a constraint to the problem + pub fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { + self.constraints_matrix_builder + .add_row(constraint.expression.linear); + let index = self.constraint_values.len(); + self.constraint_values.push(-constraint.expression.constant); + // Cones indicate the type of constraint. We only support nonnegative and equality constraints. + // To avoid creating a new cone for each constraint, we merge them. + let next_cone = if constraint.is_equality { + ZeroConeT(1) + } else { + NonnegativeConeT(1) + }; + let prev_cone = self.cones.last_mut(); + match (prev_cone, next_cone) { + (Some(ZeroConeT(a)), ZeroConeT(b)) => *a += b, + (Some(NonnegativeConeT(a)), NonnegativeConeT(b)) => *a += b, + (_, next_cone) => self.cones.push(next_cone), + }; + ConstraintReference { index } } } @@ -235,11 +194,28 @@ impl SolverModel for ClarabelProblem { fn solve(self) -> Result { let mut solver = self.into_solver(); solver.solve(); - ClarabelProblemCore::handle_solver_status(solver.solution.status, solver.solution) + match solver.solution.status { + e @ (SolverStatus::PrimalInfeasible | SolverStatus::AlmostPrimalInfeasible) => { + eprintln!("Clarabel error: {:?}", e); + Err(ResolutionError::Infeasible) + } + SolverStatus::Solved + | SolverStatus::AlmostSolved + | SolverStatus::AlmostDualInfeasible + | SolverStatus::DualInfeasible => Ok(ClarabelSolution { + solution: solver.solution, + }), + SolverStatus::Unsolved => Err(ResolutionError::Other("Unsolved")), + SolverStatus::MaxIterations => Err(ResolutionError::Other("Max iterations reached")), + SolverStatus::MaxTime => Err(ResolutionError::Other("Time limit reached")), + SolverStatus::NumericalError => Err(ResolutionError::Other("Numerical error")), + SolverStatus::InsufficientProgress => Err(ResolutionError::Other("No progress")), + SolverStatus::CallbackTerminated => Err(ResolutionError::Other("Callback terminated")), + } } fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { - self.core.add_constraint_impl(constraint) + self.add_constraint(constraint) } fn name() -> &'static str { @@ -350,22 +326,18 @@ fn fast_flatten_vecs(vecs: Vec>) -> Vec { impl ClarabelQuadraticProblem { /// Access the problem settings pub fn settings(&mut self) -> &mut DefaultSettingsBuilder { - &mut self.core.settings + &mut self.settings } /// Convert the problem into a clarabel solver /// panics if the problem is not valid pub fn into_solver(self) -> DefaultSolver { - let settings = self - .core - .settings - .build() - .expect("Invalid clarabel settings"); + let settings = self.settings.build().expect("Invalid clarabel settings"); let quadratic_objective = &self.quadratic_matrix_builder.build(); - let objective = &self.core.objective; - let constraints = &self.core.constraints_matrix_builder.build(); - let constraint_values = &self.core.constraint_values; - let cones = &self.core.cones; + let objective = &self.objective; + let constraints = &self.constraints_matrix_builder.build(); + let constraint_values = &self.constraint_values; + let cones = &self.cones; DefaultSolver::new( quadratic_objective, objective, @@ -375,12 +347,27 @@ impl ClarabelQuadraticProblem { settings, ).expect("Invalid clarabel problem. This is likely a bug in good_lp. Problems should always have coherent dimensions.") } -} -#[cfg(feature = "enable_quadratic")] -impl ClarabelProblemLike for ClarabelQuadraticProblem { - fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { - self.core.add_constraint_impl(constraint) + /// Add a constraint to the quadratic problem + pub fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { + self.constraints_matrix_builder + .add_row(constraint.expression.linear); + let index = self.constraint_values.len(); + self.constraint_values.push(-constraint.expression.constant); + // Cones indicate the type of constraint. We only support nonnegative and equality constraints. + // To avoid creating a new cone for each constraint, we merge them. + let next_cone = if constraint.is_equality { + ZeroConeT(1) + } else { + NonnegativeConeT(1) + }; + let prev_cone = self.cones.last_mut(); + match (prev_cone, next_cone) { + (Some(ZeroConeT(a)), ZeroConeT(b)) => *a += b, + (Some(NonnegativeConeT(a)), NonnegativeConeT(b)) => *a += b, + (_, next_cone) => self.cones.push(next_cone), + }; + ConstraintReference { index } } } @@ -392,11 +379,28 @@ impl SolverModel for ClarabelQuadraticProblem { fn solve(self) -> Result { let mut solver = self.into_solver(); solver.solve(); - ClarabelProblemCore::handle_solver_status(solver.solution.status, solver.solution) + match solver.solution.status { + e @ (SolverStatus::PrimalInfeasible | SolverStatus::AlmostPrimalInfeasible) => { + eprintln!("Clarabel error: {:?}", e); + Err(ResolutionError::Infeasible) + } + SolverStatus::Solved + | SolverStatus::AlmostSolved + | SolverStatus::AlmostDualInfeasible + | SolverStatus::DualInfeasible => Ok(ClarabelSolution { + solution: solver.solution, + }), + SolverStatus::Unsolved => Err(ResolutionError::Other("Unsolved")), + SolverStatus::MaxIterations => Err(ResolutionError::Other("Max iterations reached")), + SolverStatus::MaxTime => Err(ResolutionError::Other("Time limit reached")), + SolverStatus::NumericalError => Err(ResolutionError::Other("Numerical error")), + SolverStatus::InsufficientProgress => Err(ResolutionError::Other("No progress")), + SolverStatus::CallbackTerminated => Err(ResolutionError::Other("Callback terminated")), + } } fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { - self.core.add_constraint_impl(constraint) + self.add_constraint(constraint) } fn name() -> &'static str { From 68aae107f11f5a736cc448a550b76082c82727a3 Mon Sep 17 00:00:00 2001 From: Shurikal Date: Sat, 27 Sep 2025 09:29:22 +0200 Subject: [PATCH 6/8] Reduced clutter in quadratic_expression --- src/expression.rs | 8 + src/quadratic_expression.rs | 607 +++--------------------------------- 2 files changed, 48 insertions(+), 567 deletions(-) diff --git a/src/expression.rs b/src/expression.rs index 95a83ba..438e388 100644 --- a/src/expression.rs +++ b/src/expression.rs @@ -82,6 +82,14 @@ impl Clone for LinearExpression { } } +impl Debug for LinearExpression { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + f.debug_struct("LinearExpression") + .field("coefficients", &self.coefficients) + .finish() + } +} + /// Represents an affine expression, such as `2x + 3` or `x + y + z` pub struct Expression { pub(crate) linear: LinearExpression, diff --git a/src/quadratic_expression.rs b/src/quadratic_expression.rs index 40a6037..ba7e2f8 100644 --- a/src/quadratic_expression.rs +++ b/src/quadratic_expression.rs @@ -1,17 +1,14 @@ -use std::fmt::{Debug, Formatter}; -use std::ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub, SubAssign}; +use std::fmt::Debug; +use std::ops::Mul; use fnv::FnvHashMap as HashMap; use crate::affine_expression_trait::IntoAffineExpression; - use crate::expression::{Expression, LinearExpression}; -use crate::variable::{FormatWithVars, Variable}; +use crate::variable::Variable; use crate::Solution; -/// Represents a quadratic term as a pair of variables and their coefficient -/// For example, 3*x*y would be represented as QuadraticTerm { var1: x, var2: y, coeff: 3.0 } -/// For x^2, both var1 and var2 would be x +/// Represents a pair of variables in a quadratic term #[derive(Clone, Copy, PartialEq, Eq, Hash, Debug)] pub struct VariablePair { /// First variable in the pair @@ -23,7 +20,6 @@ pub struct VariablePair { impl VariablePair { /// Create a new variable pair, ensuring consistent ordering for commutativity (x*y = y*x) pub fn new(var1: Variable, var2: Variable) -> Self { - // Ensure consistent ordering for commutativity (x*y = y*x) if var1.index() <= var2.index() { VariablePair { var1, var2 } } else { @@ -33,14 +29,6 @@ impl VariablePair { } } } - - /// Create a variable pair representing a square term (var^2) - pub fn square(var: Variable) -> Self { - VariablePair { - var1: var, - var2: var, - } - } } /// A quadratic expression without linear or constant components @@ -56,50 +44,12 @@ impl QuadraticExpression { } } - /// Create a quadratic expression with preallocated capacity - pub fn with_capacity(capacity: usize) -> Self { - QuadraticExpression { - quadratic_coefficients: HashMap::with_capacity_and_hasher(capacity, Default::default()), - } - } - /// Add a quadratic term to this expression pub fn add_quadratic_term(&mut self, var1: Variable, var2: Variable, coefficient: f64) { let pair = VariablePair::new(var1, var2); *self.quadratic_coefficients.entry(pair).or_default() += coefficient; } - /// Get the coefficient for a specific variable pair - pub fn get_quadratic_coefficient(&self, var1: Variable, var2: Variable) -> f64 { - let pair = VariablePair::new(var1, var2); - self.quadratic_coefficients - .get(&pair) - .copied() - .unwrap_or(0.0) - } - - /// Check if the expression is empty (all coefficients are zero) - pub fn is_empty(&self) -> bool { - self.quadratic_coefficients.is_empty() - || self - .quadratic_coefficients - .values() - .all(|&c| c.abs() < f64::EPSILON) - } - - /// Iterate over all non-zero quadratic terms - pub fn iter(&self) -> impl Iterator + '_ { - self.quadratic_coefficients - .iter() - .filter_map(|(&pair, &coeff)| { - if coeff.abs() >= f64::EPSILON { - Some((pair, coeff)) - } else { - None - } - }) - } - /// Evaluate the quadratic expression given variable values pub fn eval_with(&self, values: &S) -> f64 { self.quadratic_coefficients @@ -128,7 +78,6 @@ impl Clone for QuadraticExpression { } /// A complete quadratic expression containing quadratic, linear, and constant terms -/// Represents expressions like: x^2 + 2*x*y + 3*x + 4 pub struct QuadraticAffineExpression { pub(crate) quadratic: QuadraticExpression, pub(crate) linear: LinearExpression, @@ -147,44 +96,6 @@ impl QuadraticAffineExpression { } } - /// Create a quadratic affine expression with preallocated capacity - pub fn with_capacity(linear_capacity: usize, quadratic_capacity: usize) -> Self { - QuadraticAffineExpression { - quadratic: QuadraticExpression::with_capacity(quadratic_capacity), - linear: LinearExpression { - coefficients: HashMap::with_capacity_and_hasher( - linear_capacity, - Default::default(), - ), - }, - constant: 0.0, - } - } - - /// Create from an affine expression - pub fn from_affine(expr: E) -> Self { - let constant = expr.constant(); - let linear_coeffs = expr.linear_coefficients().into_iter().collect(); - QuadraticAffineExpression { - quadratic: QuadraticExpression::new(), - linear: LinearExpression { - coefficients: linear_coeffs, - }, - constant, - } - } - - /// Create from a quadratic expression - pub fn from_quadratic(expr: QuadraticExpression) -> Self { - QuadraticAffineExpression { - quadratic: expr, - linear: LinearExpression { - coefficients: HashMap::default(), - }, - constant: 0.0, - } - } - /// Add a quadratic term pub fn add_quadratic_term(&mut self, var1: Variable, var2: Variable, coefficient: f64) { self.quadratic.add_quadratic_term(var1, var2, coefficient); @@ -200,83 +111,35 @@ impl QuadraticAffineExpression { self.constant += value; } - /// Get the linear part as an Expression - pub fn linear_part(&self) -> Expression { - Expression { + /// Create from a quadratic expression + pub fn from_quadratic(expr: QuadraticExpression) -> Self { + QuadraticAffineExpression { + quadratic: expr, linear: LinearExpression { - coefficients: self.linear.coefficients.clone(), + coefficients: HashMap::default(), }, - constant: self.constant, + constant: 0.0, } } - /// Check if the expression contains any quadratic terms - pub fn is_purely_affine(&self) -> bool { - self.quadratic.is_empty() + /// Create from an affine expression + pub fn from_affine(expr: E) -> Self { + let constant = expr.constant(); + let linear_coeffs = expr.linear_coefficients().into_iter().collect(); + QuadraticAffineExpression { + quadratic: QuadraticExpression::new(), + linear: LinearExpression { + coefficients: linear_coeffs, + }, + constant, + } } /// Evaluate the complete expression pub fn eval_with(&self, values: &S) -> f64 { - self.quadratic.eval_with(values) + self.linear_part().eval_with(values) - } - - /// Get the coefficient of a quadratic term - pub fn get_quadratic_coefficient(&self, var1: Variable, var2: Variable) -> f64 { - self.quadratic.get_quadratic_coefficient(var1, var2) - } - - /// Get the coefficient of a linear term - pub fn get_linear_coefficient(&self, var: Variable) -> f64 { - self.linear.coefficients.get(&var).copied().unwrap_or(0.0) - } - - /// Get the constant term - pub fn get_constant(&self) -> f64 { - self.constant - } - - /// Get a reference to the quadratic part - pub fn quadratic_part(&self) -> &QuadraticExpression { - &self.quadratic - } - - /// Creates a constraint indicating that this quadratic expression - /// is lesser than or equal to the right hand side - pub fn leq(self, rhs: RHS) -> QuadraticConstraint - where - Self: Sub, - { - let diff = self - rhs; - QuadraticConstraint { - expression: diff, - is_equality: false, - } - } - - /// Creates a constraint indicating that this quadratic expression - /// is greater than or equal to the right hand side - pub fn geq(self, rhs: RHS) -> QuadraticConstraint - where - RHS: Sub, - { - let diff = rhs - self; - QuadraticConstraint { - expression: diff, - is_equality: false, - } - } - - /// Creates a constraint indicating that this quadratic expression - /// is equal to the right hand side - pub fn eq(self, rhs: RHS) -> QuadraticConstraint - where - Self: Sub, - { - let diff = self - rhs; - QuadraticConstraint { - expression: diff, - is_equality: true, - } + self.quadratic.eval_with(values) + + self.linear.coefficients.iter().map(|(&var, &coeff)| coeff * values.value(var)).sum::() + + self.constant } } @@ -296,6 +159,21 @@ impl Clone for QuadraticAffineExpression { } } +impl Debug for QuadraticAffineExpression { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "QuadraticAffineExpression {{ quadratic: {:?}, linear: {:?}, constant: {} }}", + self.quadratic, self.linear, self.constant) + } +} + +impl Debug for QuadraticExpression { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_struct("QuadraticExpression") + .field("quadratic_coefficients", &self.quadratic_coefficients) + .finish() + } +} + /// A quadratic constraint pub struct QuadraticConstraint { /// The quadratic expression that defines the constraint @@ -304,33 +182,6 @@ pub struct QuadraticConstraint { pub is_equality: bool, } -// Conversion implementations -impl From for QuadraticAffineExpression { - fn from(value: f64) -> Self { - QuadraticAffineExpression { - quadratic: QuadraticExpression::new(), - linear: LinearExpression { - coefficients: HashMap::default(), - }, - constant: value, - } - } -} - -impl From for QuadraticAffineExpression { - fn from(var: Variable) -> Self { - let mut linear_coeffs = HashMap::default(); - linear_coeffs.insert(var, 1.0); - QuadraticAffineExpression { - quadratic: QuadraticExpression::new(), - linear: LinearExpression { - coefficients: linear_coeffs, - }, - constant: 0.0, - } - } -} - impl From for QuadraticAffineExpression { fn from(expr: Expression) -> Self { QuadraticAffineExpression { @@ -341,200 +192,7 @@ impl From for QuadraticAffineExpression { } } -impl From for QuadraticAffineExpression { - fn from(quad_expr: QuadraticExpression) -> Self { - QuadraticAffineExpression::from_quadratic(quad_expr) - } -} - -// Arithmetic operations for QuadraticExpression -impl> MulAssign for QuadraticExpression { - fn mul_assign(&mut self, rhs: N) { - let factor = rhs.into(); - for value in self.quadratic_coefficients.values_mut() { - *value *= factor; - } - } -} - -impl> Mul for QuadraticExpression { - type Output = QuadraticExpression; - - fn mul(mut self, rhs: N) -> Self::Output { - self.mul_assign(rhs); - self - } -} - -impl> Div for QuadraticExpression { - type Output = QuadraticExpression; - - fn div(mut self, rhs: N) -> Self::Output { - self.mul_assign(1. / rhs.into()); - self - } -} - -impl AddAssign for QuadraticExpression { - fn add_assign(&mut self, rhs: QuadraticExpression) { - for (pair, coeff) in rhs.quadratic_coefficients { - *self.quadratic_coefficients.entry(pair).or_default() += coeff; - } - } -} - -impl Add for QuadraticExpression { - type Output = QuadraticExpression; - - fn add(mut self, rhs: QuadraticExpression) -> Self::Output { - self.add_assign(rhs); - self - } -} - -impl SubAssign for QuadraticExpression { - fn sub_assign(&mut self, rhs: QuadraticExpression) { - for (pair, coeff) in rhs.quadratic_coefficients { - *self.quadratic_coefficients.entry(pair).or_default() -= coeff; - } - } -} - -impl Sub for QuadraticExpression { - type Output = QuadraticExpression; - - fn sub(mut self, rhs: QuadraticExpression) -> Self::Output { - self.sub_assign(rhs); - self - } -} - -impl Neg for QuadraticExpression { - type Output = Self; - - fn neg(mut self) -> Self::Output { - self *= -1; - self - } -} - -// Arithmetic operations for QuadraticAffineExpression -impl> MulAssign for QuadraticAffineExpression { - fn mul_assign(&mut self, rhs: N) { - let factor = rhs.into(); - self.quadratic.mul_assign(factor); - for value in self.linear.coefficients.values_mut() { - *value *= factor; - } - self.constant *= factor; - } -} - -impl> Mul for QuadraticAffineExpression { - type Output = QuadraticAffineExpression; - - fn mul(mut self, rhs: N) -> Self::Output { - self.mul_assign(rhs); - self - } -} - -impl> Div for QuadraticAffineExpression { - type Output = QuadraticAffineExpression; - - fn div(mut self, rhs: N) -> Self::Output { - self.mul_assign(1. / rhs.into()); - self - } -} - -impl AddAssign for QuadraticAffineExpression { - fn add_assign(&mut self, rhs: QuadraticAffineExpression) { - self.quadratic.add_assign(rhs.quadratic); - for (var, coeff) in rhs.linear.coefficients { - *self.linear.coefficients.entry(var).or_default() += coeff; - } - self.constant += rhs.constant; - } -} - -impl Add for QuadraticAffineExpression { - type Output = QuadraticAffineExpression; - - fn add(mut self, rhs: QuadraticAffineExpression) -> Self::Output { - self.add_assign(rhs); - self - } -} - -impl SubAssign for QuadraticAffineExpression { - fn sub_assign(&mut self, rhs: QuadraticAffineExpression) { - self.quadratic.sub_assign(rhs.quadratic); - for (var, coeff) in rhs.linear.coefficients { - *self.linear.coefficients.entry(var).or_default() -= coeff; - } - self.constant -= rhs.constant; - } -} - -impl Sub for QuadraticAffineExpression { - type Output = QuadraticAffineExpression; - - fn sub(mut self, rhs: QuadraticAffineExpression) -> Self::Output { - self.sub_assign(rhs); - self - } -} - -impl Neg for QuadraticAffineExpression { - type Output = Self; - - fn neg(mut self) -> Self::Output { - self *= -1; - self - } -} - -// Mixed type operations -impl AddAssign for QuadraticAffineExpression { - fn add_assign(&mut self, rhs: E) { - let constant = rhs.constant(); - for (var, coeff) in rhs.linear_coefficients() { - *self.linear.coefficients.entry(var).or_default() += coeff; - } - self.constant += constant; - } -} - -impl Add for QuadraticAffineExpression { - type Output = QuadraticAffineExpression; - - fn add(mut self, rhs: E) -> Self::Output { - self.add_assign(rhs); - self - } -} - -impl SubAssign for QuadraticAffineExpression { - fn sub_assign(&mut self, rhs: E) { - let constant = rhs.constant(); - for (var, coeff) in rhs.linear_coefficients() { - *self.linear.coefficients.entry(var).or_default() -= coeff; - } - self.constant -= constant; - } -} - -impl Sub for QuadraticAffineExpression { - type Output = QuadraticAffineExpression; - - fn sub(mut self, rhs: E) -> Self::Output { - self.sub_assign(rhs); - self - } -} - -// Multiplication operations that create quadratic terms +// Essential multiplication operations for creating quadratic terms impl Mul for Variable { type Output = QuadraticExpression; @@ -544,188 +202,3 @@ impl Mul for Variable { quadratic } } - -// Specific implementations for Variable multiplication with other types -impl Mul for Variable { - type Output = QuadraticAffineExpression; - - fn mul(self, rhs: Expression) -> Self::Output { - let mut result = QuadraticAffineExpression::new(); - let constant = rhs.constant; - - // Linear term (constant * variable) - if constant.abs() >= f64::EPSILON { - result.add_linear_term(self, constant); - } - - // Quadratic terms (variable * other_variables) - for (var, coeff) in rhs.linear.coefficients { - result.add_quadratic_term(self, var, coeff); - } - - result - } -} - -impl Mul for Expression { - type Output = QuadraticAffineExpression; - - fn mul(self, rhs: Variable) -> Self::Output { - rhs * self - } -} - -// Expression multiplication (returns quadratic expression) -impl Mul for Expression { - type Output = QuadraticAffineExpression; - - fn mul(self, rhs: Expression) -> Self::Output { - let mut result = QuadraticAffineExpression::new(); - - let c1 = self.constant; - let c2 = rhs.constant; - - // Constant term (c1 * c2) - result.constant = c1 * c2; - - // Linear terms from self * constant of rhs - for (var, coeff) in self.linear.coefficients.iter() { - result.add_linear_term(*var, *coeff * c2); - } - - // Linear terms from constant of self * rhs - for (var, coeff) in rhs.linear.coefficients.iter() { - result.add_linear_term(*var, c1 * *coeff); - } - - // Quadratic terms (variable from self * variable from rhs) - for (var1, coeff1) in self.linear.coefficients.iter() { - for (var2, coeff2) in rhs.linear.coefficients.iter() { - result.add_quadratic_term(*var1, *var2, *coeff1 * *coeff2); - } - } - - result - } -} - -// Numeric multiplication -macro_rules! impl_mul_numeric { - ($($t:ty),*) => {$( - impl Mul for $t { - type Output = QuadraticExpression; - - fn mul(self, mut rhs: QuadraticExpression) -> Self::Output { - rhs *= self; - rhs - } - } - - impl Mul for $t { - type Output = QuadraticAffineExpression; - - fn mul(self, mut rhs: QuadraticAffineExpression) -> Self::Output { - rhs *= self; - rhs - } - } - )*} -} -impl_mul_numeric!(f64, i32, f32, u32, u16, u8, i16, i8); - -// Display and formatting support -impl FormatWithVars for QuadraticExpression { - fn format_with(&self, f: &mut Formatter<'_>, mut variable_format: FUN) -> std::fmt::Result - where - FUN: FnMut(&mut Formatter<'_>, Variable) -> std::fmt::Result, - { - let mut first = true; - for (pair, coeff) in self.iter() { - if first { - first = false; - } else { - write!(f, " + ")?; - } - - if (coeff - 1.).abs() > f64::EPSILON { - write!(f, "{} ", coeff)?; - } - - if pair.var1 == pair.var2 { - // Square term: x^2 - variable_format(f, pair.var1)?; - write!(f, "²")?; - } else { - // Cross term: x*y - variable_format(f, pair.var1)?; - write!(f, "*")?; - variable_format(f, pair.var2)?; - } - } - - if first { - write!(f, "0")?; - } - Ok(()) - } -} - -impl Debug for QuadraticExpression { - fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { - self.format_debug(f) - } -} - -impl FormatWithVars for QuadraticAffineExpression { - fn format_with(&self, f: &mut Formatter<'_>, mut variable_format: FUN) -> std::fmt::Result - where - FUN: FnMut(&mut Formatter<'_>, Variable) -> std::fmt::Result, - { - let has_quadratic = !self.quadratic.is_empty(); - let has_linear = !self.linear.coefficients.is_empty(); - let has_constant = self.constant.abs() >= f64::EPSILON; - - if !has_quadratic && !has_linear && !has_constant { - return write!(f, "0"); - } - - let mut first = true; - - // Write quadratic terms - if has_quadratic { - self.quadratic.format_with(f, &mut variable_format)?; - first = false; - } - - // Write linear terms - for (&var, &coeff) in &self.linear.coefficients { - if coeff.abs() >= f64::EPSILON { - if !first { - write!(f, " + ")?; - } - first = false; - - if (coeff - 1.).abs() > f64::EPSILON { - write!(f, "{} ", coeff)?; - } - variable_format(f, var)?; - } - } - - // Write constant term - if has_constant { - if !first { - write!(f, " + ")?; - } - write!(f, "{}", self.constant)?; - } - - Ok(()) - } -} - -impl Debug for QuadraticAffineExpression { - fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { - self.format_debug(f) - } -} From e3782f8247ad95f6664fecd73c51d513cde81c9b Mon Sep 17 00:00:00 2001 From: Shurikal Date: Tue, 30 Sep 2025 07:46:17 +0200 Subject: [PATCH 7/8] Implement quadratic expression trait and refactor quadratic problem handling - Introduced `IntoQuadraticExpression` trait for converting various types (e.g., `Expression`, `Variable`, numeric constants) into quadratic expressions. - Removed the `QuadraticUnsolvedProblem` struct and integrated its functionality into the existing `UnsolvedProblem` struct. - Updated solvers (`clarabel`, `coin_cbc`, `lp_solve`, `microlp`) to handle quadratic objectives appropriately, with checks and panic messages for unsupported cases. - Refactored the `clarabel` solver to support both linear and quadratic objectives, consolidating the logic for building the objective matrix. - Enhanced the `ProblemVariables` struct with new methods for creating optimization problems with quadratic objectives. - Updated tests to validate the new quadratic functionality and ensure compatibility with existing linear problem-solving methods. --- src/affine_expression_trait.rs | 15 +- src/expression.rs | 581 ++++++++++++++++++++++++++++-- src/lib.rs | 14 +- src/quadratic_expression.rs | 204 ----------- src/quadratic_expression_trait.rs | 54 +++ src/quadratic_problem.rs | 104 ------ src/solvers/clarabel.rs | 427 ++++++++-------------- src/solvers/coin_cbc.rs | 5 + src/solvers/lpsolve.rs | 3 + src/solvers/microlp.rs | 3 + src/variable.rs | 60 +++ tests/quadratic_solver_tests.rs | 75 ++-- tests/variables.rs | 2 +- 13 files changed, 884 insertions(+), 663 deletions(-) delete mode 100644 src/quadratic_expression.rs create mode 100644 src/quadratic_expression_trait.rs delete mode 100644 src/quadratic_problem.rs diff --git a/src/affine_expression_trait.rs b/src/affine_expression_trait.rs index 9a20364..9a6dd35 100644 --- a/src/affine_expression_trait.rs +++ b/src/affine_expression_trait.rs @@ -2,9 +2,14 @@ //! You can implement this trait if you want to implement your own //! variant of the [Expression](crate::Expression) type, optimized for your use case. use crate::expression::LinearExpression; +#[cfg(feature = "enable_quadratic")] +use crate::expression::QuadraticExpression; use crate::{Expression, Solution, Variable}; /// An element that can be expressed as a linear combination of variables plus a constant +/// +/// This trait can only be used with expressions that contain no quadratic terms. +/// Attempting to use it with an Expression containing quadratic terms will panic. pub trait IntoAffineExpression { /// The iterator returned by [`linear_coefficients`](IntoAffineExpression::linear_coefficients). type Iter: IntoIterator; @@ -26,9 +31,13 @@ pub trait IntoAffineExpression { Self: Sized, { let constant = self.constant(); - let coefficients = self.linear_coefficients().into_iter().collect(); + let linear_coefficients = self.linear_coefficients().into_iter().collect(); Expression { - linear: LinearExpression { coefficients }, + #[cfg(feature = "enable_quadratic")] + quadratic: QuadraticExpression::new(), + linear: LinearExpression { + coefficients: linear_coefficients, + }, constant, } } @@ -92,6 +101,8 @@ macro_rules! impl_affine_for_num { fn into_expression(self) -> Expression { Expression { + #[cfg(feature = "enable_quadratic")] + quadratic: QuadraticExpression::new(), linear: LinearExpression { coefficients: std::default::Default::default() }, constant: f64::from(self), } diff --git a/src/expression.rs b/src/expression.rs index 438e388..c9915ce 100644 --- a/src/expression.rs +++ b/src/expression.rs @@ -5,9 +5,124 @@ use fnv::FnvHashMap as HashMap; use crate::affine_expression_trait::IntoAffineExpression; use crate::constraint; +#[cfg(feature = "enable_quadratic")] +use crate::IntoQuadraticExpression; use crate::variable::{FormatWithVars, Variable}; use crate::{Constraint, Solution}; +/// Represents a pair of variables in a quadratic term +#[cfg(feature = "enable_quadratic")] +#[derive(Clone, Copy, PartialEq, Eq, Hash, Debug)] +pub struct VariablePair { + /// First variable in the pair + pub var1: Variable, + /// Second variable in the pair + pub var2: Variable, +} + +#[cfg(feature = "enable_quadratic")] +impl VariablePair { + /// Create a new variable pair, ensuring consistent ordering for commutativity (x*y = y*x) + pub fn new(var1: Variable, var2: Variable) -> Self { + if var1.index() <= var2.index() { + VariablePair { var1, var2 } + } else { + VariablePair { + var1: var2, + var2: var1, + } + } + } +} + +/// A pure quadratic expression containing only quadratic terms (similar to LinearExpression) +#[cfg(feature = "enable_quadratic")] +pub struct QuadraticExpression { + pub(crate) coefficients: HashMap, +} + +#[cfg(feature = "enable_quadratic")] +impl QuadraticExpression { + /// Create a new empty quadratic expression + pub fn new() -> Self { + QuadraticExpression { + coefficients: HashMap::default(), + } + } + + /// Add a quadratic term + pub fn add_quadratic_term(&mut self, var1: Variable, var2: Variable, coefficient: f64) { + let pair = VariablePair::new(var1, var2); + *self.coefficients.entry(pair).or_default() += coefficient; + } + + /// Evaluate the quadratic terms with given variable values + pub fn eval_with(&self, values: &S) -> f64 { + self.coefficients + .iter() + .map(|(pair, &coeff)| { + let val1 = values.value(pair.var1); + let val2 = values.value(pair.var2); + coeff * val1 * val2 + }) + .sum::() + } +} + +#[cfg(feature = "enable_quadratic")] +impl Default for QuadraticExpression { + fn default() -> Self { + Self::new() + } +} + +#[cfg(feature = "enable_quadratic")] +impl Clone for QuadraticExpression { + fn clone(&self) -> Self { + QuadraticExpression { + coefficients: self.coefficients.clone(), + } + } +} + +#[cfg(feature = "enable_quadratic")] +impl Debug for QuadraticExpression { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_struct("QuadraticExpression") + .field("coefficients", &self.coefficients) + .finish() + } +} + +#[cfg(feature = "enable_quadratic")] +impl FormatWithVars for QuadraticExpression { + fn format_with(&self, f: &mut Formatter<'_>, mut variable_format: FUN) -> std::fmt::Result + where + FUN: FnMut(&mut Formatter<'_>, Variable) -> std::fmt::Result, + { + let mut first = true; + for (pair, &coeff) in &self.coefficients { + if coeff != 0f64 { + if first { + first = false; + } else { + write!(f, " + ")?; + } + if (coeff - 1.).abs() > f64::EPSILON { + write!(f, "{} ", coeff)?; + } + variable_format(f, pair.var1)?; + write!(f, "*")?; + variable_format(f, pair.var2)?; + } + } + if first { + write!(f, "0")?; + } + Ok(()) + } +} + /// An linear expression without a constant component pub struct LinearExpression { pub(crate) coefficients: HashMap, @@ -91,7 +206,12 @@ impl Debug for LinearExpression { } /// Represents an affine expression, such as `2x + 3` or `x + y + z` +/// if quadratic features are disabled. Else it can also represent +/// `x^2 + y^2 + xy + 2x + 3`. pub struct Expression { + #[cfg(feature = "enable_quadratic")] + pub(crate) quadratic: QuadraticExpression, + pub(crate) linear: LinearExpression, pub(crate) constant: f64, } @@ -101,6 +221,12 @@ impl IntoAffineExpression for Expression { #[inline] fn linear_coefficients(self) -> Self::Iter { + #[cfg(feature = "enable_quadratic")] + { + if !self.quadratic.coefficients.is_empty() { + panic!("Cannot convert quadratic expression to affine expression: expression contains quadratic terms"); + } + } self.linear.linear_coefficients() } @@ -117,6 +243,12 @@ impl<'a> IntoAffineExpression for &'a Expression { #[inline] fn linear_coefficients(self) -> Self::Iter { + #[cfg(feature = "enable_quadratic")] + { + if !self.is_affine() { + panic!("Cannot convert quadratic expression to affine expression: expression contains quadratic terms"); + } + } (&self.linear).linear_coefficients() } @@ -128,13 +260,28 @@ impl<'a> IntoAffineExpression for &'a Expression { impl PartialEq for Expression { fn eq(&self, other: &Self) -> bool { - self.constant.eq(&other.constant) && self.linear.coefficients.eq(&other.linear.coefficients) + let base_eq = self.constant == other.constant + && self.linear.coefficients == other.linear.coefficients; + + #[cfg(feature = "enable_quadratic")] + { + base_eq && self.quadratic.coefficients == other.quadratic.coefficients + } + + #[cfg(not(feature = "enable_quadratic"))] + { + base_eq + } } } impl Clone for Expression { fn clone(&self) -> Self { Expression { + #[cfg(feature = "enable_quadratic")] + quadratic: QuadraticExpression { + coefficients: self.quadratic.coefficients.clone(), + }, linear: LinearExpression { coefficients: self.linear.coefficients.clone(), }, @@ -160,6 +307,10 @@ impl Expression { /// for `capacity` coefficients. pub fn with_capacity(capacity: usize) -> Self { Expression { + #[cfg(feature = "enable_quadratic")] + quadratic: QuadraticExpression { + coefficients: HashMap::with_capacity_and_hasher(capacity, Default::default()), + }, linear: LinearExpression { coefficients: HashMap::with_capacity_and_hasher(capacity, Default::default()), }, @@ -177,6 +328,12 @@ impl Expression { source.into_expression() } + /// Create a concrete expression struct from anything that can be expressed as a quadratic expression + #[cfg(feature = "enable_quadratic")] + pub fn from_other_quadratic(source: E) -> Self { + source.into_quadratic_expression() + } + /// Creates a constraint indicating that this expression /// is lesser than or equal to the right hand side pub fn leq(self, rhs: RHS) -> Constraint @@ -214,7 +371,55 @@ impl Expression { /// See [IntoAffineExpression::eval_with] pub fn eval_with(&self, values: &S) -> f64 { - IntoAffineExpression::eval_with(self, values) + #[cfg(feature = "enable_quadratic")] + { + // Manually evaluate linear terms to avoid IntoAffineExpression panic with quadratic terms + let mut result = self.constant; + for (&var, &coeff) in &self.linear.coefficients { + result += coeff * values.value(var); + } + result += self.quadratic.eval_with(values); + result + } + #[cfg(not(feature = "enable_quadratic"))] + { + IntoAffineExpression::eval_with(self, values) + } + } + + /// Add a quadratic term to this expression + #[cfg(feature = "enable_quadratic")] + pub fn add_quadratic_term(&mut self, var1: Variable, var2: Variable, coefficient: f64) { + self.quadratic.add_quadratic_term(var1, var2, coefficient); + } + + /// Add a linear term to this expression + pub fn add_linear_term(&mut self, var: Variable, coefficient: f64) { + *self.linear.coefficients.entry(var).or_default() += coefficient; + } + + /// Add a constant term to this expression + pub fn add_constant(&mut self, value: f64) { + self.constant += value; + } + + /// Returns true if this expression contains no quadratic terms + /// (either because quadratic features are disabled or all quadratic coefficients are zero) + pub fn is_affine(&self) -> bool { + #[cfg(feature = "enable_quadratic")] + { + self.quadratic.coefficients.is_empty() + || self + .quadratic + .coefficients + .values() + .all(|&coeff| coeff.abs() < f64::EPSILON) + } + + #[cfg(not(feature = "enable_quadratic"))] + { + true + } } } @@ -239,15 +444,102 @@ pub fn add, RHS: IntoAffineExpression>(lhs: LHS, rhs: RHS) add_mul(lhs, rhs, 1.) } +// Macro for implementing arithmetic operations between Expression and numeric types +macro_rules! impl_expr_ops_for_num { + ($($num:ty),*) => {$( + impl Mul<$num> for Expression { + type Output = Expression; + fn mul(mut self, rhs: $num) -> Self::Output { + self *= rhs; + self + } + } + + impl Div<$num> for Expression { + type Output = Expression; + fn div(mut self, rhs: $num) -> Self::Output { + self *= 1.0 / f64::from(rhs); + self + } + } + )*}; +} + +// Macro for implementing reverse operations (num op Expression) +macro_rules! impl_num_ops_for_expr { + ($($num:ty),*) => {$( + impl Add for $num { + type Output = Expression; + fn add(self, rhs: Expression) -> Self::Output { + rhs + self + } + } + + impl Mul for $num { + type Output = Expression; + fn mul(self, mut rhs: Expression) -> Self::Output { + rhs *= self; + rhs + } + } + )*}; +} + +// Apply the macros to generate implementations for numeric types +impl_expr_ops_for_num!(f64, f32, u32, u16, u8, i32, i16, i8); +impl_num_ops_for_expr!(f64, f32, u32, u16, u8, i32, i16, i8); + impl FormatWithVars for Expression { - fn format_with(&self, f: &mut Formatter<'_>, variable_format: FUN) -> std::fmt::Result + fn format_with(&self, f: &mut Formatter<'_>, mut variable_format: FUN) -> std::fmt::Result where FUN: FnMut(&mut Formatter<'_>, Variable) -> std::fmt::Result, { - self.linear.format_with(f, variable_format)?; + let mut has_terms = false; + + // Format linear terms first + for (&var, &coeff) in &self.linear.coefficients { + if coeff != 0f64 { + if has_terms { + write!(f, " + ")?; + } + if (coeff - 1.).abs() > f64::EPSILON { + write!(f, "{} ", coeff)?; + } + variable_format(f, var)?; + has_terms = true; + } + } + + // Format quadratic terms + #[cfg(feature = "enable_quadratic")] + { + for (pair, &coeff) in &self.quadratic.coefficients { + if coeff != 0f64 { + if has_terms { + write!(f, " + ")?; + } + if (coeff - 1.).abs() > f64::EPSILON { + write!(f, "{} ", coeff)?; + } + variable_format(f, pair.var1)?; + write!(f, "*")?; + variable_format(f, pair.var2)?; + has_terms = true; + } + } + } + + // Format constant if self.constant.abs() >= f64::EPSILON { - write!(f, " + {}", self.constant)?; + if has_terms { + write!(f, " + {}", self.constant)?; + } else { + write!(f, "{}", self.constant)?; + } + } else if !has_terms { + write!(f, "0")?; } + Ok(()) } } @@ -280,6 +572,12 @@ impl> MulAssign for Expression { #[inline] fn mul_assign(&mut self, rhs: N) { let factor = rhs.into(); + #[cfg(feature = "enable_quadratic")] + { + for value in self.quadratic.coefficients.values_mut() { + *value *= factor; + } + } for value in self.linear.coefficients.values_mut() { *value *= factor } @@ -287,64 +585,166 @@ impl> MulAssign for Expression { } } -impl> Mul for Expression { +// Macro for implementing additional reverse operations for extended numeric types +macro_rules! impl_extended_num_ops { + ($($num:ty),*) => {$( + impl Add for $num { + type Output = Expression; + fn add(self, rhs: Variable) -> Self::Output { + Expression::from(self) + rhs + } + } + + impl Sub for $num { + type Output = Expression; + fn sub(self, rhs: Variable) -> Self::Output { + Expression::from(self) - rhs + } + } + + impl Sub for $num { + type Output = Expression; + fn sub(self, rhs: Expression) -> Self::Output { + Expression::from(self) - rhs + } + } + )*}; +} + +impl_extended_num_ops!(f64, f32, u32, u16, u8, i32, i16, i8); + +// Specific implementations for Expression operations to handle quadratic terms +impl Add for Expression { type Output = Expression; - #[inline] - fn mul(mut self, rhs: N) -> Self::Output { - self.mul_assign(rhs); + fn add(mut self, rhs: Expression) -> Self::Output { + // Add quadratic terms + #[cfg(feature = "enable_quadratic")] + { + for (var_pair, coeff) in rhs.quadratic.coefficients { + self.add_quadratic_term(var_pair.var1, var_pair.var2, coeff); + } + } + // Add linear terms + for (var, coeff) in rhs.linear.coefficients { + self.add_linear_term(var, coeff); + } + // Add constant + self.add_constant(rhs.constant); self } } -impl> Div for Expression { +impl Sub for Expression { type Output = Expression; - #[inline] - fn div(mut self, rhs: N) -> Self::Output { - self.mul_assign(1. / rhs.into()); + fn sub(mut self, rhs: Expression) -> Self::Output { + // Subtract quadratic terms + #[cfg(feature = "enable_quadratic")] + { + for (var_pair, coeff) in rhs.quadratic.coefficients { + self.add_quadratic_term(var_pair.var1, var_pair.var2, -coeff); + } + } + // Subtract linear terms + for (var, coeff) in rhs.linear.coefficients { + self.add_linear_term(var, -coeff); + } + // Subtract constant + self.add_constant(-rhs.constant); self } } -macro_rules! impl_mul { - ($($t:ty),*) =>{$( - impl Mul for $t { + +macro_rules! impl_var_num_ops { + ($($num:ty),*) => {$( + impl Add<$num> for Variable { type Output = Expression; + fn add(self, rhs: $num) -> Self::Output { + Expression::from(self) + rhs + } + } - fn mul(self, mut rhs: Expression) -> Self::Output { - rhs *= self; - rhs + impl Sub<$num> for Variable { + type Output = Expression; + fn sub(self, rhs: $num) -> Self::Output { + Expression::from(self) - rhs } } - )*} + )*}; +} + +impl_var_num_ops!(f64, f32, u32, u16, u8, i32, i16, i8); + +// Variable operations with Expression +impl Add for Variable { + type Output = Expression; + fn add(self, rhs: Expression) -> Self::Output { + Expression::from(self) + rhs + } +} + +impl Sub for Variable { + type Output = Expression; + fn sub(self, rhs: Expression) -> Self::Output { + Expression::from(self) - rhs + } +} + +impl Add for Variable { + type Output = Expression; + fn add(self, rhs: Variable) -> Self::Output { + Expression::from(self) + rhs + } +} + +impl Sub for Variable { + type Output = Expression; + fn sub(self, rhs: Variable) -> Self::Output { + Expression::from(self) - rhs + } +} + +// Operations with Variable references +impl Add<&Variable> for Variable { + type Output = Expression; + fn add(self, rhs: &Variable) -> Self::Output { + Expression::from(self) + *rhs + } +} + +impl Sub<&Variable> for Variable { + type Output = Expression; + fn sub(self, rhs: &Variable) -> Self::Output { + Expression::from(self) - *rhs + } } -impl_mul!(f64, i32); -macro_rules! impl_ops_local { - ($( $typename:ident : $([generic $generic:ident])? $other:ident),*) => {$( - impl<$($generic: IntoAffineExpression)?> Sub<$other> - for $typename { +// Operations with Option +impl Add> for Variable { + type Output = Expression; + fn add(self, rhs: Option) -> Self::Output { + Expression::from(self) + Expression::from_other_affine(rhs) + } +} + +// Macro for implementing Expression operations with affine types +macro_rules! impl_expr_affine_ops { + ($($affine_type:ty),*) => {$( + impl Sub<$affine_type> for Expression { type Output = Expression; - fn sub(self, rhs: $other) -> Self::Output { sub(self, rhs) } + fn sub(self, rhs: $affine_type) -> Self::Output { sub(self, rhs) } } - impl<$($generic: IntoAffineExpression)?> Add<$other> - for $typename { + impl Add<$affine_type> for Expression { type Output = Expression; - fn add(self, rhs: $other) -> Self::Output { add(self, rhs) } + fn add(self, rhs: $affine_type) -> Self::Output { add(self, rhs) } } - )*} + )*}; } -impl_ops_local!( - Expression: [generic RHS] RHS, - Variable: [generic RHS] RHS, - f64: Expression, - f64: Variable, - i32: Expression, - i32: Variable -); +impl_expr_affine_ops!(Variable, f64, f32, u32, u16, u8, i32, i16, i8); macro_rules! impl_conv { ( $( $typename:ident ),* ) => {$( @@ -353,7 +753,7 @@ macro_rules! impl_conv { } )*} } -impl_conv!(f64, i32, Variable); +impl_conv!(f64, f32, u32, u16, u8, i32, i16, i8, Variable); impl std::iter::Sum for Expression { fn sum>(iter: I) -> Self { @@ -366,6 +766,89 @@ impl std::iter::Sum for Expression { } } +// Essential multiplication operations for creating quadratic terms +#[cfg(feature = "enable_quadratic")] +impl Mul for Variable { + type Output = Expression; + + fn mul(self, rhs: Variable) -> Self::Output { + let mut expr = Expression::default(); + expr.add_quadratic_term(self, rhs, 1.0); + expr + } +} + +// Multiplication operations to create quadratic terms from linear expressions +#[cfg(feature = "enable_quadratic")] +impl Mul for Expression { + type Output = Expression; + + fn mul(self, rhs: Variable) -> Self::Output { + let mut result = Expression::default(); + + // Each linear term * variable becomes a quadratic term + for (var, coeff) in self.linear.coefficients { + result.add_quadratic_term(var, rhs, coeff); + } + + // Constant * variable becomes a linear term + if self.constant.abs() >= f64::EPSILON { + result.add_linear_term(rhs, self.constant); + } + + result + } +} + +#[cfg(feature = "enable_quadratic")] +impl Mul for Variable { + type Output = Expression; + + fn mul(self, rhs: Expression) -> Self::Output { + rhs * self + } +} + +// Expression * Expression multiplication for creating complex quadratic terms +// This enables syntax like: (x - 5.0) * (x + 3.0) or (x + y) * (x - y) +#[cfg(feature = "enable_quadratic")] +impl Mul for Expression { + type Output = Expression; + + fn mul(self, rhs: Expression) -> Self::Output { + let mut result = Expression::default(); + + // Expand (a₁x₁ + a₂x₂ + ... + c₁) * (b₁y₁ + b₂y₂ + ... + c₂) + + // Linear terms × Linear terms = Quadratic terms + for (&var1, &coeff1) in &self.linear.coefficients { + for (&var2, &coeff2) in &rhs.linear.coefficients { + result.add_quadratic_term(var1, var2, coeff1 * coeff2); + } + } + + // Linear terms × Constant = Linear terms + for (&var, &coeff) in &self.linear.coefficients { + if rhs.constant.abs() >= f64::EPSILON { + result.add_linear_term(var, coeff * rhs.constant); + } + } + + for (&var, &coeff) in &rhs.linear.coefficients { + if self.constant.abs() >= f64::EPSILON { + result.add_linear_term(var, self.constant * coeff); + } + } + + // Constant × Constant = Constant + if self.constant.abs() >= f64::EPSILON && rhs.constant.abs() >= f64::EPSILON { + result.add_constant(self.constant * rhs.constant); + } + + result + } +} + #[cfg(test)] mod tests { use std::collections::HashMap; @@ -387,6 +870,24 @@ mod tests { let mut values = HashMap::new(); values.insert(a, 100); values.insert(b, -1); - assert_eq!((a + 3 * (b + 3)).eval_with(&values), 106.) + assert_eq!((a + 3.0_f64 * (b + 3.0_f64 )).eval_with(&values), 106.) + } + + #[cfg(feature = "enable_quadratic")] + #[allow(clippy::float_cmp)] + #[test] + fn expression_manipulation_quadratic() { + variables! {vars: v0; v1; } + let expression_0 = (3.0_f64 - v0) - v1 * v1; + let expression_1 = (-1.0) * v0 + (-1.0) * v1 * v1 + 3.0; + let expression_0_str = format!("{:?}", expression_0); + assert_eq!(expression_0_str, "-1 v0 + -1 v1*v1 + 3"); + println!("{:?}", expression_0); + println!("{:?}", expression_1); + assert_eq!( + expression_0, expression_1 + ); + let values = HashMap::from([(v0, 100.), (v1, -1.)]); + assert_eq!(expression_0.eval_with(&values), -98.0) } } diff --git a/src/lib.rs b/src/lib.rs index 5f772f1..50b8958 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -70,11 +70,9 @@ pub use cardinality_constraint_solver_trait::CardinalityConstraintSolver; pub use constraint::Constraint; pub use expression::Expression; #[cfg(feature = "enable_quadratic")] -pub use quadratic_expression::{ - QuadraticAffineExpression, QuadraticConstraint, QuadraticExpression, VariablePair, -}; +pub use expression::VariablePair; #[cfg(feature = "enable_quadratic")] -pub use quadratic_problem::{IntoQuadraticExpression, QuadraticUnsolvedProblem}; +pub use quadratic_expression_trait::IntoQuadraticExpression; #[cfg(not(any( feature = "coin_cbc", feature = "microlp", @@ -92,8 +90,6 @@ pub use solvers::clarabel::clarabel; docsrs, doc(cfg(all(feature = "clarabel", feature = "enable_quadratic"))) )] -#[cfg(all(feature = "clarabel", feature = "enable_quadratic"))] -pub use solvers::clarabel::clarabel_quadratic; #[cfg_attr(docsrs, doc(cfg(feature = "coin_cbc")))] #[cfg(feature = "coin_cbc")] pub use solvers::coin_cbc::coin_cbc; @@ -200,11 +196,9 @@ pub mod variable; mod affine_expression_trait; mod cardinality_constraint_solver_trait; pub mod constraint; -/// Quadratic expression types and operations #[cfg(feature = "enable_quadratic")] -pub mod quadratic_expression; +mod quadratic_expression_trait; + /// Quadratic problem definition and solving -#[cfg(feature = "enable_quadratic")] -pub mod quadratic_problem; pub mod solvers; mod variables_macro; diff --git a/src/quadratic_expression.rs b/src/quadratic_expression.rs deleted file mode 100644 index ba7e2f8..0000000 --- a/src/quadratic_expression.rs +++ /dev/null @@ -1,204 +0,0 @@ -use std::fmt::Debug; -use std::ops::Mul; - -use fnv::FnvHashMap as HashMap; - -use crate::affine_expression_trait::IntoAffineExpression; -use crate::expression::{Expression, LinearExpression}; -use crate::variable::Variable; -use crate::Solution; - -/// Represents a pair of variables in a quadratic term -#[derive(Clone, Copy, PartialEq, Eq, Hash, Debug)] -pub struct VariablePair { - /// First variable in the pair - pub var1: Variable, - /// Second variable in the pair - pub var2: Variable, -} - -impl VariablePair { - /// Create a new variable pair, ensuring consistent ordering for commutativity (x*y = y*x) - pub fn new(var1: Variable, var2: Variable) -> Self { - if var1.index() <= var2.index() { - VariablePair { var1, var2 } - } else { - VariablePair { - var1: var2, - var2: var1, - } - } - } -} - -/// A quadratic expression without linear or constant components -pub struct QuadraticExpression { - pub(crate) quadratic_coefficients: HashMap, -} - -impl QuadraticExpression { - /// Create a new empty quadratic expression - pub fn new() -> Self { - QuadraticExpression { - quadratic_coefficients: HashMap::default(), - } - } - - /// Add a quadratic term to this expression - pub fn add_quadratic_term(&mut self, var1: Variable, var2: Variable, coefficient: f64) { - let pair = VariablePair::new(var1, var2); - *self.quadratic_coefficients.entry(pair).or_default() += coefficient; - } - - /// Evaluate the quadratic expression given variable values - pub fn eval_with(&self, values: &S) -> f64 { - self.quadratic_coefficients - .iter() - .map(|(pair, &coeff)| { - let val1 = values.value(pair.var1); - let val2 = values.value(pair.var2); - coeff * val1 * val2 - }) - .sum() - } -} - -impl Default for QuadraticExpression { - fn default() -> Self { - Self::new() - } -} - -impl Clone for QuadraticExpression { - fn clone(&self) -> Self { - QuadraticExpression { - quadratic_coefficients: self.quadratic_coefficients.clone(), - } - } -} - -/// A complete quadratic expression containing quadratic, linear, and constant terms -pub struct QuadraticAffineExpression { - pub(crate) quadratic: QuadraticExpression, - pub(crate) linear: LinearExpression, - pub(crate) constant: f64, -} - -impl QuadraticAffineExpression { - /// Create a new empty quadratic affine expression - pub fn new() -> Self { - QuadraticAffineExpression { - quadratic: QuadraticExpression::new(), - linear: LinearExpression { - coefficients: HashMap::default(), - }, - constant: 0.0, - } - } - - /// Add a quadratic term - pub fn add_quadratic_term(&mut self, var1: Variable, var2: Variable, coefficient: f64) { - self.quadratic.add_quadratic_term(var1, var2, coefficient); - } - - /// Add a linear term - pub fn add_linear_term(&mut self, var: Variable, coefficient: f64) { - *self.linear.coefficients.entry(var).or_default() += coefficient; - } - - /// Add a constant term - pub fn add_constant(&mut self, value: f64) { - self.constant += value; - } - - /// Create from a quadratic expression - pub fn from_quadratic(expr: QuadraticExpression) -> Self { - QuadraticAffineExpression { - quadratic: expr, - linear: LinearExpression { - coefficients: HashMap::default(), - }, - constant: 0.0, - } - } - - /// Create from an affine expression - pub fn from_affine(expr: E) -> Self { - let constant = expr.constant(); - let linear_coeffs = expr.linear_coefficients().into_iter().collect(); - QuadraticAffineExpression { - quadratic: QuadraticExpression::new(), - linear: LinearExpression { - coefficients: linear_coeffs, - }, - constant, - } - } - - /// Evaluate the complete expression - pub fn eval_with(&self, values: &S) -> f64 { - self.quadratic.eval_with(values) - + self.linear.coefficients.iter().map(|(&var, &coeff)| coeff * values.value(var)).sum::() - + self.constant - } -} - -impl Default for QuadraticAffineExpression { - fn default() -> Self { - Self::new() - } -} - -impl Clone for QuadraticAffineExpression { - fn clone(&self) -> Self { - QuadraticAffineExpression { - quadratic: self.quadratic.clone(), - linear: self.linear.clone(), - constant: self.constant, - } - } -} - -impl Debug for QuadraticAffineExpression { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - write!(f, "QuadraticAffineExpression {{ quadratic: {:?}, linear: {:?}, constant: {} }}", - self.quadratic, self.linear, self.constant) - } -} - -impl Debug for QuadraticExpression { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - f.debug_struct("QuadraticExpression") - .field("quadratic_coefficients", &self.quadratic_coefficients) - .finish() - } -} - -/// A quadratic constraint -pub struct QuadraticConstraint { - /// The quadratic expression that defines the constraint - pub expression: QuadraticAffineExpression, - /// Whether this is an equality constraint (true) or inequality constraint (false) - pub is_equality: bool, -} - -impl From for QuadraticAffineExpression { - fn from(expr: Expression) -> Self { - QuadraticAffineExpression { - quadratic: QuadraticExpression::new(), - linear: expr.linear, - constant: expr.constant, - } - } -} - -// Essential multiplication operations for creating quadratic terms -impl Mul for Variable { - type Output = QuadraticExpression; - - fn mul(self, rhs: Variable) -> Self::Output { - let mut quadratic = QuadraticExpression::new(); - quadratic.add_quadratic_term(self, rhs, 1.0); - quadratic - } -} diff --git a/src/quadratic_expression_trait.rs b/src/quadratic_expression_trait.rs new file mode 100644 index 0000000..bf67c38 --- /dev/null +++ b/src/quadratic_expression_trait.rs @@ -0,0 +1,54 @@ +//! A quadratic expression trait for expressions that can be converted to quadratic form. +//! A quadratic expression contains linear terms, quadratic terms, and constants. + +#[cfg(feature = "enable_quadratic")] +use crate::{Expression, Variable}; + +/// An element that can be expressed as a quadratic expression +/// (combination of quadratic terms, linear terms, and constants) +#[cfg(feature = "enable_quadratic")] +pub trait IntoQuadraticExpression { + /// Convert this element into a quadratic expression + fn into_quadratic_expression(self) -> Expression + where + Self: Sized; +} + +// Implement the trait for Expression (which can now contain quadratic terms) +#[cfg(feature = "enable_quadratic")] +impl IntoQuadraticExpression for Expression { + fn into_quadratic_expression(self) -> Expression { + self + } +} + +// Implement the trait for variables +#[cfg(feature = "enable_quadratic")] +impl IntoQuadraticExpression for Variable { + fn into_quadratic_expression(self) -> Expression { + Expression::from(self) + } +} + +// Implement the trait for numeric constants +#[cfg(feature = "enable_quadratic")] +impl IntoQuadraticExpression for f64 { + fn into_quadratic_expression(self) -> Expression { + Expression::from(self) + } +} + +#[cfg(feature = "enable_quadratic")] +impl IntoQuadraticExpression for i32 { + fn into_quadratic_expression(self) -> Expression { + (self as f64).into_quadratic_expression() + } +} + +// Blanket implementation for references +#[cfg(feature = "enable_quadratic")] +impl IntoQuadraticExpression for &T { + fn into_quadratic_expression(self) -> Expression { + self.clone().into_quadratic_expression() + } +} diff --git a/src/quadratic_problem.rs b/src/quadratic_problem.rs deleted file mode 100644 index f6841bf..0000000 --- a/src/quadratic_problem.rs +++ /dev/null @@ -1,104 +0,0 @@ -use crate::quadratic_expression::QuadraticAffineExpression; -use crate::solvers::ObjectiveDirection; -use crate::variable::ProblemVariables; - -/// A problem with a quadratic objective function -pub struct QuadraticUnsolvedProblem { - /// The quadratic objective function to optimize - pub objective: QuadraticAffineExpression, - /// Whether to maximize or minimize the objective - pub direction: ObjectiveDirection, - /// The variables in the problem - pub variables: ProblemVariables, -} - -/// Trait for quadratic expressions that can be used as objective functions -pub trait IntoQuadraticExpression { - /// Transform the value into a concrete QuadraticAffineExpression - fn into_quadratic_expression(self) -> QuadraticAffineExpression; -} - -impl IntoQuadraticExpression for QuadraticAffineExpression { - fn into_quadratic_expression(self) -> QuadraticAffineExpression { - self - } -} - -impl IntoQuadraticExpression for crate::QuadraticExpression { - fn into_quadratic_expression(self) -> QuadraticAffineExpression { - QuadraticAffineExpression::from_quadratic(self) - } -} - -impl IntoQuadraticExpression for T { - fn into_quadratic_expression(self) -> QuadraticAffineExpression { - QuadraticAffineExpression::from_affine(self) - } -} - -impl ProblemVariables { - /// Create a quadratic optimization problem to maximize the given quadratic objective - #[cfg(feature = "enable_quadratic")] - pub fn maximise_quadratic( - self, - objective: E, - ) -> QuadraticUnsolvedProblem { - let objective_expr = objective.into_quadratic_expression(); - // Validate that objective variables are within bounds - for (var, _) in objective_expr.linear.coefficients.iter() { - assert!( - var.index() < self.len(), - "Variable in objective function is not part of this problem" - ); - } - for (pair, _) in objective_expr.quadratic.quadratic_coefficients.iter() { - assert!( - pair.var1.index() < self.len() && pair.var2.index() < self.len(), - "Variable in quadratic objective function is not part of this problem" - ); - } - QuadraticUnsolvedProblem { - objective: objective_expr, - direction: ObjectiveDirection::Maximisation, - variables: self, - } - } - - /// Create a quadratic optimization problem to minimize the given quadratic objective - #[cfg(feature = "enable_quadratic")] - pub fn minimise_quadratic( - self, - objective: E, - ) -> QuadraticUnsolvedProblem { - let objective_expr = objective.into_quadratic_expression(); - // Validate that objective variables are within bounds - for (var, _) in objective_expr.linear.coefficients.iter() { - assert!( - var.index() < self.len(), - "Variable in objective function is not part of this problem" - ); - } - for (pair, _) in objective_expr.quadratic.quadratic_coefficients.iter() { - assert!( - pair.var1.index() < self.len() && pair.var2.index() < self.len(), - "Variable in quadratic objective function is not part of this problem" - ); - } - QuadraticUnsolvedProblem { - objective: objective_expr, - direction: ObjectiveDirection::Minimisation, - variables: self, - } - } -} - -impl QuadraticUnsolvedProblem { - /// Solve this quadratic problem using the given solver function - #[cfg(all(feature = "clarabel", feature = "enable_quadratic"))] - pub fn using(self, solver_factory: F) -> S - where - F: FnOnce(Self) -> S, - { - solver_factory(self) - } -} diff --git a/src/solvers/clarabel.rs b/src/solvers/clarabel.rs index 8a2e59d..81d2ead 100644 --- a/src/solvers/clarabel.rs +++ b/src/solvers/clarabel.rs @@ -1,12 +1,12 @@ //! A solver that uses [clarabel](https://oxfordcontrol.github.io/ClarabelDocs/stable/), a pure rust solver. use crate::affine_expression_trait::IntoAffineExpression; -use crate::expression::LinearExpression; -use crate::variable::UnsolvedProblem; #[cfg(feature = "enable_quadratic")] -use crate::quadratic_problem::QuadraticUnsolvedProblem; +use crate::expression::VariablePair; +use crate::expression::LinearExpression; #[cfg(feature = "enable_quadratic")] -use crate::quadratic_expression::VariablePair; +use crate::expression::Expression; +use crate::variable::UnsolvedProblem; use crate::{ constraint::ConstraintReference, solvers::{ObjectiveDirection, ResolutionError, Solution, SolverModel}, @@ -22,83 +22,62 @@ use clarabel::solver::{DefaultSolver, IPSolver}; /// The [clarabel](https://oxfordcontrol.github.io/ClarabelDocs/stable/) solver, /// to be used with [UnsolvedProblem::using]. +/// Automatically handles both linear and quadratic objectives. pub fn clarabel(to_solve: UnsolvedProblem) -> ClarabelProblem { - let UnsolvedProblem { - objective, - direction, - variables, - } = to_solve; - let coef = if direction == ObjectiveDirection::Maximisation { + let coef = if to_solve.direction == ObjectiveDirection::Maximisation { -1. } else { 1. }; + + let objective = to_solve.objective; + + // Now we can safely move the other fields + let variables = to_solve.variables; + let mut objective_vector = vec![0.; variables.len()]; - for (var, obj) in objective.linear_coefficients() { - objective_vector[var.index()] = obj * coef; - } - let constraints_matrix_builder = CscMatrixBuilder::new(variables.len()); - let mut settings = DefaultSettingsBuilder::default(); - settings.verbose(false).tol_feas(1e-9); - let mut p = ClarabelProblem { - objective: objective_vector, - constraints_matrix_builder, - constraint_values: Vec::new(), - variables: variables.len(), - settings, - cones: Vec::new(), - }; - // add trivial constraints embedded in the variable definitions - for (var, def) in variables.iter_variables_with_def() { - if def.is_integer { - panic!("Clarabel doesn't support integer variables") - } - if def.min != f64::NEG_INFINITY { - p.add_constraint(var >> def.min); - } - if def.max != f64::INFINITY { - p.add_constraint(var << def.max); + + // Handle quadratic case + #[cfg(feature = "enable_quadratic")] + let quadratic_matrix_builder = if !objective.is_affine() { + // Use quadratic objective coefficients + for (var, obj_coeff) in objective.linear.coefficients { + objective_vector[var.index()] = obj_coeff * coef; } - } - p -} -/// The [clarabel](https://clarabel.org/stable/) solver for quadratic problems, -/// to be used with [QuadraticUnsolvedProblem::using]. -#[cfg(feature = "enable_quadratic")] -pub fn clarabel_quadratic(to_solve: QuadraticUnsolvedProblem) -> ClarabelQuadraticProblem { - let QuadraticUnsolvedProblem { - objective, - direction, - variables, - } = to_solve; - let coef = if direction == ObjectiveDirection::Maximisation { - -1. + let mut qmb = CscMatrixBuilder::new_square(variables.len()); + for (&pair, &quad_coeff) in &objective.quadratic.coefficients { + qmb.add_quadratic_term(pair, quad_coeff * coef); + } + Some(qmb) } else { - 1. + // Use linear objective coefficients + for (var, obj) in objective.linear_coefficients() { + objective_vector[var.index()] = obj * coef; + } + None }; - let mut objective_vector = vec![0.; variables.len()]; - for (var, obj_coeff) in objective.linear.coefficients { - objective_vector[var.index()] = obj_coeff * coef; - } - - let mut quadratic_matrix_builder = CscQuadraticMatrixBuilder::new(variables.len()); - for (pair, quad_coeff) in objective.quadratic.quadratic_coefficients { - quadratic_matrix_builder.add_term(pair, quad_coeff * coef); + // Handle linear-only case (when quadratic features are disabled) + #[cfg(not(feature = "enable_quadratic"))] + { + for (var, obj) in objective.linear_coefficients() { + objective_vector[var.index()] = obj * coef; + } } let constraints_matrix_builder = CscMatrixBuilder::new(variables.len()); let mut settings = DefaultSettingsBuilder::default(); settings.verbose(false).tol_feas(1e-9); - let mut p = ClarabelQuadraticProblem { + let mut p = ClarabelProblem { objective: objective_vector, constraints_matrix_builder, constraint_values: Vec::new(), variables: variables.len(), settings, cones: Vec::new(), + #[cfg(feature = "enable_quadratic")] quadratic_matrix_builder, }; @@ -117,7 +96,7 @@ pub fn clarabel_quadratic(to_solve: QuadraticUnsolvedProblem) -> ClarabelQuadrat p } -/// A clarabel model +/// A clarabel model that can handle both linear and quadratic objectives pub struct ClarabelProblem { constraints_matrix_builder: CscMatrixBuilder, constraint_values: Vec, @@ -125,18 +104,8 @@ pub struct ClarabelProblem { variables: usize, settings: DefaultSettingsBuilder, cones: Vec>, -} - -/// A clarabel quadratic model -#[cfg(feature = "enable_quadratic")] -pub struct ClarabelQuadraticProblem { - constraints_matrix_builder: CscMatrixBuilder, - constraint_values: Vec, - objective: Vec, - variables: usize, - settings: DefaultSettingsBuilder, - cones: Vec>, - quadratic_matrix_builder: CscQuadraticMatrixBuilder, + #[cfg(feature = "enable_quadratic")] + quadratic_matrix_builder: Option, } impl ClarabelProblem { @@ -149,13 +118,23 @@ impl ClarabelProblem { /// panics if the problem is not valid pub fn into_solver(self) -> DefaultSolver { let settings = self.settings.build().expect("Invalid clarabel settings"); - let quadratic_objective = &CscMatrix::zeros((self.variables, self.variables)); + + #[cfg(feature = "enable_quadratic")] + let quadratic_objective = if let Some(qmb) = self.quadratic_matrix_builder { + qmb.build() + } else { + CscMatrix::zeros((self.variables, self.variables)) + }; + + #[cfg(not(feature = "enable_quadratic"))] + let quadratic_objective = CscMatrix::zeros((self.variables, self.variables)); + let objective = &self.objective; let constraints = &self.constraints_matrix_builder.build(); let constraint_values = &self.constraint_values; let cones = &self.cones; DefaultSolver::new( - quadratic_objective, + &quadratic_objective, objective, constraints, constraint_values, @@ -163,28 +142,6 @@ impl ClarabelProblem { settings, ).expect("Invalid clarabel problem. This is likely a bug in good_lp. Problems should always have coherent dimensions.") } - - /// Add a constraint to the problem - pub fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { - self.constraints_matrix_builder - .add_row(constraint.expression.linear); - let index = self.constraint_values.len(); - self.constraint_values.push(-constraint.expression.constant); - // Cones indicate the type of constraint. We only support nonnegative and equality constraints. - // To avoid creating a new cone for each constraint, we merge them. - let next_cone = if constraint.is_equality { - ZeroConeT(1) - } else { - NonnegativeConeT(1) - }; - let prev_cone = self.cones.last_mut(); - match (prev_cone, next_cone) { - (Some(ZeroConeT(a)), ZeroConeT(b)) => *a += b, - (Some(NonnegativeConeT(a)), NonnegativeConeT(b)) => *a += b, - (_, next_cone) => self.cones.push(next_cone), - }; - ConstraintReference { index } - } } impl SolverModel for ClarabelProblem { @@ -215,7 +172,24 @@ impl SolverModel for ClarabelProblem { } fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { - self.add_constraint(constraint) + self.constraints_matrix_builder + .add_row(constraint.expression.linear); + let index = self.constraint_values.len(); + self.constraint_values.push(-constraint.expression.constant); + // Cones indicate the type of constraint. We only support nonnegative and equality constraints. + // To avoid creating a new cone for each constraint, we merge them. + let next_cone = if constraint.is_equality { + ZeroConeT(1) + } else { + NonnegativeConeT(1) + }; + let prev_cone = self.cones.last_mut(); + match (prev_cone, next_cone) { + (Some(ZeroConeT(a)), ZeroConeT(b)) => *a += b, + (Some(NonnegativeConeT(a)), NonnegativeConeT(b)) => *a += b, + (_, next_cone) => self.cones.push(next_cone), + }; + ConstraintReference { index } } fn name() -> &'static str { @@ -270,6 +244,7 @@ struct CscMatrixBuilder { nzval: Vec>, n_rows: usize, n_cols: usize, + is_quadratic: bool, } impl CscMatrixBuilder { @@ -279,8 +254,22 @@ impl CscMatrixBuilder { nzval: vec![Vec::new(); n_cols], n_rows: 0, n_cols, + is_quadratic: false, + } + } + + /// Create a new builder for a square matrix (used for quadratic matrices) + #[cfg(feature = "enable_quadratic")] + fn new_square(n_vars: usize) -> Self { + Self { + rowval: vec![Vec::new(); n_vars], + nzval: vec![Vec::new(); n_vars], + n_rows: n_vars, + n_cols: n_vars, + is_quadratic: true, } } + fn add_row(&mut self, row: LinearExpression) { for (var, value) in row.linear_coefficients() { self.rowval[var.index()].push(self.n_rows); @@ -288,7 +277,47 @@ impl CscMatrixBuilder { } self.n_rows += 1; } - fn build(self) -> clarabel::algebra::CscMatrix { + + /// Add a quadratic term to the matrix (for quadratic objectives) + #[cfg(feature = "enable_quadratic")] + fn add_quadratic_term(&mut self, pair: VariablePair, coeff: f64) { + let i = pair.var1.index(); + let j = pair.var2.index(); + + if i == j { + // Diagonal term: multiply by 2 because Clarabel minimizes 1/2 * x^T * P * x + // So we need P_ii = 2 * coeff to get coeff * x_i^2 in the objective + self.rowval[j].push(i); + self.nzval[j].push(2.0 * coeff); + } else { + // Off-diagonal term: add the coefficient to both (i,j) and (j,i) positions + // to create a symmetric matrix. Each entry gets the original coefficient + // since Clarabel multiplies the entire quadratic form by 1/2. + self.rowval[j].push(i); + self.nzval[j].push(coeff); + + self.rowval[i].push(j); + self.nzval[i].push(coeff); + } + } + + fn build(mut self) -> clarabel::algebra::CscMatrix { + // For quadratic matrices, sort columns to maintain proper CSC format + #[cfg(feature = "enable_quadratic")] + if self.is_quadratic { + for col in 0..self.n_cols { + let mut pairs: Vec<(usize, f64)> = self.rowval[col] + .iter() + .zip(self.nzval[col].iter()) + .map(|(&row, &val)| (row, val)) + .collect(); + pairs.sort_by_key(|&(row, _)| row); + + self.rowval[col] = pairs.iter().map(|&(row, _)| row).collect(); + self.nzval[col] = pairs.iter().map(|&(_, val)| val).collect(); + } + } + let mut colptr = Vec::with_capacity(self.n_cols + 1); colptr.push(0); for col in &self.rowval { @@ -322,161 +351,7 @@ fn fast_flatten_vecs(vecs: Vec>) -> Vec { result } -#[cfg(feature = "enable_quadratic")] -impl ClarabelQuadraticProblem { - /// Access the problem settings - pub fn settings(&mut self) -> &mut DefaultSettingsBuilder { - &mut self.settings - } - - /// Convert the problem into a clarabel solver - /// panics if the problem is not valid - pub fn into_solver(self) -> DefaultSolver { - let settings = self.settings.build().expect("Invalid clarabel settings"); - let quadratic_objective = &self.quadratic_matrix_builder.build(); - let objective = &self.objective; - let constraints = &self.constraints_matrix_builder.build(); - let constraint_values = &self.constraint_values; - let cones = &self.cones; - DefaultSolver::new( - quadratic_objective, - objective, - constraints, - constraint_values, - cones, - settings, - ).expect("Invalid clarabel problem. This is likely a bug in good_lp. Problems should always have coherent dimensions.") - } - /// Add a constraint to the quadratic problem - pub fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { - self.constraints_matrix_builder - .add_row(constraint.expression.linear); - let index = self.constraint_values.len(); - self.constraint_values.push(-constraint.expression.constant); - // Cones indicate the type of constraint. We only support nonnegative and equality constraints. - // To avoid creating a new cone for each constraint, we merge them. - let next_cone = if constraint.is_equality { - ZeroConeT(1) - } else { - NonnegativeConeT(1) - }; - let prev_cone = self.cones.last_mut(); - match (prev_cone, next_cone) { - (Some(ZeroConeT(a)), ZeroConeT(b)) => *a += b, - (Some(NonnegativeConeT(a)), NonnegativeConeT(b)) => *a += b, - (_, next_cone) => self.cones.push(next_cone), - }; - ConstraintReference { index } - } -} - -#[cfg(feature = "enable_quadratic")] -impl SolverModel for ClarabelQuadraticProblem { - type Solution = ClarabelSolution; - type Error = ResolutionError; - - fn solve(self) -> Result { - let mut solver = self.into_solver(); - solver.solve(); - match solver.solution.status { - e @ (SolverStatus::PrimalInfeasible | SolverStatus::AlmostPrimalInfeasible) => { - eprintln!("Clarabel error: {:?}", e); - Err(ResolutionError::Infeasible) - } - SolverStatus::Solved - | SolverStatus::AlmostSolved - | SolverStatus::AlmostDualInfeasible - | SolverStatus::DualInfeasible => Ok(ClarabelSolution { - solution: solver.solution, - }), - SolverStatus::Unsolved => Err(ResolutionError::Other("Unsolved")), - SolverStatus::MaxIterations => Err(ResolutionError::Other("Max iterations reached")), - SolverStatus::MaxTime => Err(ResolutionError::Other("Time limit reached")), - SolverStatus::NumericalError => Err(ResolutionError::Other("Numerical error")), - SolverStatus::InsufficientProgress => Err(ResolutionError::Other("No progress")), - SolverStatus::CallbackTerminated => Err(ResolutionError::Other("Callback terminated")), - } - } - - fn add_constraint(&mut self, constraint: Constraint) -> ConstraintReference { - self.add_constraint(constraint) - } - - fn name() -> &'static str { - "Clarabel Quadratic" - } -} - -/// Builder for symmetric quadratic matrices in CSC format -#[cfg(feature = "enable_quadratic")] -struct CscQuadraticMatrixBuilder { - /// Indicates the row index of the corresponding element in `nzval` - rowval: Vec>, - /// All non-zero values in the matrix, in column-major order - nzval: Vec>, - n_vars: usize, -} - -#[cfg(feature = "enable_quadratic")] -impl CscQuadraticMatrixBuilder { - fn new(n_vars: usize) -> Self { - Self { - rowval: vec![Vec::new(); n_vars], - nzval: vec![Vec::new(); n_vars], - n_vars, - } - } - - fn add_term(&mut self, pair: VariablePair, coeff: f64) { - let i = pair.var1.index(); - let j = pair.var2.index(); - - if i == j { - // Diagonal term: multiply by 2 because Clarabel minimizes 1/2 * x^T * P * x - // So we need P_ii = 2 * coeff to get coeff * x_i^2 in the objective - self.rowval[j].push(i); - self.nzval[j].push(2.0 * coeff); - } else { - // Off-diagonal term: add the coefficient to both (i,j) and (j,i) positions - // to create a symmetric matrix. Each entry gets the original coefficient - // since Clarabel multiplies the entire quadratic form by 1/2. - self.rowval[j].push(i); - self.nzval[j].push(coeff); - - self.rowval[i].push(j); - self.nzval[i].push(coeff); - } - } - - fn build(mut self) -> clarabel::algebra::CscMatrix { - // Sort each column by row index to maintain proper CSC format - for col in 0..self.n_vars { - let mut pairs: Vec<(usize, f64)> = self.rowval[col] - .iter() - .zip(self.nzval[col].iter()) - .map(|(&row, &val)| (row, val)) - .collect(); - pairs.sort_by_key(|&(row, _)| row); - - self.rowval[col] = pairs.iter().map(|&(row, _)| row).collect(); - self.nzval[col] = pairs.iter().map(|&(_, val)| val).collect(); - } - - let mut colptr = Vec::with_capacity(self.n_vars + 1); - colptr.push(0); - for col in &self.rowval { - colptr.push(colptr.last().unwrap() + col.len()); - } - clarabel::algebra::CscMatrix::new( - self.n_vars, - self.n_vars, - colptr, - fast_flatten_vecs(self.rowval), - fast_flatten_vecs(self.nzval), - ) - } -} #[cfg(test)] mod tests { @@ -512,7 +387,7 @@ mod tests { #[cfg(feature = "enable_quadratic")] mod quadratic_tests { use super::*; - use crate::{QuadraticAffineExpression, VariablePair}; + use crate::VariablePair; #[test] fn test_csc_quadratic_matrix_builder_diagonal() { @@ -521,11 +396,11 @@ mod tests { y; } - let mut builder = CscQuadraticMatrixBuilder::new(2); + let mut builder = CscMatrixBuilder::new_square(2); // Add diagonal terms: x^2 with coefficient 2.0, y^2 with coefficient 3.0 - builder.add_term(VariablePair::new(x, x), 2.0); - builder.add_term(VariablePair::new(y, y), 3.0); + builder.add_quadratic_term(VariablePair::new(x, x), 2.0); + builder.add_quadratic_term(VariablePair::new(y, y), 3.0); let matrix = builder.build(); @@ -545,10 +420,10 @@ mod tests { y; } - let mut builder = CscQuadraticMatrixBuilder::new(3); + let mut builder = CscMatrixBuilder::new_square(3); // Add off-diagonal term: 2*x*y (coefficient 2.0) - builder.add_term(VariablePair::new(x, y), 2.0); + builder.add_quadratic_term(VariablePair::new(x, y), 2.0); let matrix = builder.build(); @@ -569,12 +444,12 @@ mod tests { y; } - let mut builder = CscQuadraticMatrixBuilder::new(2); + let mut builder = CscMatrixBuilder::new_square(2); // Add mixed terms: x^2 + 2*x*y + y^2 (like (x+y)^2) - builder.add_term(VariablePair::new(x, x), 1.0); - builder.add_term(VariablePair::new(x, y), 2.0); - builder.add_term(VariablePair::new(y, y), 1.0); + builder.add_quadratic_term(VariablePair::new(x, x), 1.0); + builder.add_quadratic_term(VariablePair::new(x, y), 2.0); + builder.add_quadratic_term(VariablePair::new(y, y), 1.0); let matrix = builder.build(); @@ -596,14 +471,12 @@ mod tests { } // Create a simple quadratic objective: minimize x^2 + y^2 - let mut quadratic_obj = QuadraticAffineExpression::new(); + let mut quadratic_obj = Expression::default(); quadratic_obj.add_quadratic_term(x, x, 1.0); // x^2 quadratic_obj.add_quadratic_term(y, y, 1.0); // y^2 // Create the problem - let problem = vars - .minimise_quadratic(quadratic_obj) - .using(clarabel_quadratic); + let problem = vars.minimise(quadratic_obj).using(clarabel); // The problem should be solvable (unconstrained minimum at (0,0)) let solution = problem.solve().expect("Should solve successfully"); @@ -621,14 +494,14 @@ mod tests { } // Create quadratic objective: minimize x^2 + y^2 - let mut quadratic_obj = QuadraticAffineExpression::new(); + let mut quadratic_obj = Expression::default(); quadratic_obj.add_quadratic_term(x, x, 1.0); // x^2 quadratic_obj.add_quadratic_term(y, y, 1.0); // y^2 // Solve with constraint x + y >= 2 let solution = vars - .minimise_quadratic(quadratic_obj) - .using(clarabel_quadratic) + .minimise(quadratic_obj) + .using(clarabel) .with(x + y >> 2.0) .solve() .expect("Should solve successfully"); @@ -651,7 +524,7 @@ mod tests { // Create quadratic objective: minimize (x-1)^2 + (y-2)^2 + 2*x*y // Expanded: x^2 - 2x + 1 + y^2 - 4y + 4 + 2xy // = x^2 + y^2 + 2xy - 2x - 4y + 5 - let mut quadratic_obj = QuadraticAffineExpression::new(); + let mut quadratic_obj = Expression::default(); quadratic_obj.add_quadratic_term(x, x, 1.0); // x^2 quadratic_obj.add_quadratic_term(y, y, 1.0); // y^2 quadratic_obj.add_quadratic_term(x, y, 2.0); // 2*x*y @@ -660,8 +533,8 @@ mod tests { quadratic_obj.add_constant(5.0); // +5 let solution = vars - .minimise_quadratic(quadratic_obj) - .using(clarabel_quadratic) + .minimise(quadratic_obj) + .using(clarabel) .solve() .expect("Should solve successfully"); @@ -682,14 +555,14 @@ mod tests { } // Minimize x^2 + y^2 + z^2 subject to x + y + z = 3 - let mut quadratic_obj = QuadraticAffineExpression::new(); + let mut quadratic_obj = Expression::default(); quadratic_obj.add_quadratic_term(x, x, 1.0); quadratic_obj.add_quadratic_term(y, y, 1.0); quadratic_obj.add_quadratic_term(z, z, 1.0); let solution = vars - .minimise_quadratic(quadratic_obj) - .using(clarabel_quadratic) + .minimise(quadratic_obj) + .using(clarabel) .with((x + y + z).eq(3.0)) // Equality constraint .solve() .expect("Should solve successfully"); @@ -708,14 +581,14 @@ mod tests { } // Create quadratic expression: x^2 + 2*x*y + y^2 (which is (x+y)^2) - let mut quadratic_obj = QuadraticAffineExpression::new(); + let mut quadratic_obj = Expression::default(); quadratic_obj.add_quadratic_term(x, x, 1.0); quadratic_obj.add_quadratic_term(x, y, 2.0); quadratic_obj.add_quadratic_term(y, y, 1.0); let solution = vars - .minimise_quadratic(quadratic_obj) - .using(clarabel_quadratic) + .minimise(quadratic_obj) + .using(clarabel) .with(x + y >> 2.0) // x + y >= 2 .solve() .expect("Should solve successfully"); diff --git a/src/solvers/coin_cbc.rs b/src/solvers/coin_cbc.rs index 75a8ef7..2224efd 100644 --- a/src/solvers/coin_cbc.rs +++ b/src/solvers/coin_cbc.rs @@ -27,6 +27,11 @@ pub fn coin_cbc(to_solve: UnsolvedProblem) -> CoinCbcProblem { direction, variables, } = to_solve; + + #[cfg(feature = "enable_quadratic")] + if !objective.is_affine() { + panic!("coin_cbc does not support quadratic objectives"); + } let mut model = Model::default(); let mut initial_solution = Vec::with_capacity(variables.initial_solution_len()); let columns: Vec = variables diff --git a/src/solvers/lpsolve.rs b/src/solvers/lpsolve.rs index 33deb96..addb1a2 100644 --- a/src/solvers/lpsolve.rs +++ b/src/solvers/lpsolve.rs @@ -33,6 +33,9 @@ fn col_num(var: Variable) -> c_int { /// The [lp_solve](http://lpsolve.sourceforge.net/5.5/) open-source solver library. /// lp_solve is released under the LGPL license. pub fn lp_solve(to_solve: UnsolvedProblem) -> LpSolveProblem { + #[cfg(feature = "enable_quadratic")] + panic!("lp_solve does not support quadratic objectives"); + let UnsolvedProblem { objective, direction, diff --git a/src/solvers/microlp.rs b/src/solvers/microlp.rs index 56e6788..21e2787 100644 --- a/src/solvers/microlp.rs +++ b/src/solvers/microlp.rs @@ -12,6 +12,9 @@ use crate::{Constraint, Variable}; /// The [microlp](https://docs.rs/microlp) solver, /// to be used with [UnsolvedProblem::using]. pub fn microlp(to_solve: UnsolvedProblem) -> MicroLpProblem { + #[cfg(feature = "enable_quadratic")] + panic!("microlp does not support quadratic objectives"); + let UnsolvedProblem { objective, direction, diff --git a/src/variable.rs b/src/variable.rs index d389c4e..937ef93 100644 --- a/src/variable.rs +++ b/src/variable.rs @@ -8,10 +8,16 @@ use std::hash::Hash; use std::iter::{repeat, FromIterator}; use std::ops::{Div, Mul, Neg, Not, RangeBounds}; +#[cfg(feature = "enable_quadratic")] +use crate::expression::QuadraticExpression; +#[cfg(feature = "enable_quadratic")] +use crate::IntoQuadraticExpression; + use fnv::FnvHashMap as HashMap; use crate::affine_expression_trait::IntoAffineExpression; use crate::expression::{Expression, LinearExpression}; + use crate::solvers::{ObjectiveDirection, Solver}; /// A variable in a problem. Use variables to create [expressions](Expression), @@ -497,6 +503,35 @@ impl ProblemVariables { } } + /// Creates an optimization problem with the given quadratic objective. Don't solve it immediately. + #[cfg(feature = "enable_quadratic")] + pub fn optimise_quadratic( + self, + direction: ObjectiveDirection, + objective: E + ) -> UnsolvedProblem { + let objective = Expression::from_other_quadratic(objective); + // todo, correct assert check + UnsolvedProblem { + objective, + direction, + variables: self, + } + } + + /// Creates an maximization problem with the given objective. Don't solve it immediately + /// + /// ``` + /// use good_lp::{variables, variable, default_solver, SolverModel, Solution}; + /// variables!{problem: x <= 7;} + /// let solution = problem.maximise(x).using(default_solver).solve().unwrap(); + /// # use float_eq::assert_float_eq; + /// assert_float_eq!(solution.value(x), 7., abs <= 1e-8); + /// ``` + #[cfg(feature = "enable_quadratic")] + pub fn maximise(self, objective: E) -> UnsolvedProblem { + self.optimise_quadratic(ObjectiveDirection::Maximisation, objective) + } /// Creates an maximization problem with the given objective. Don't solve it immediately /// /// ``` @@ -506,6 +541,7 @@ impl ProblemVariables { /// # use float_eq::assert_float_eq; /// assert_float_eq!(solution.value(x), 7., abs <= 1e-8); /// ``` + #[cfg(not(feature = "enable_quadratic"))] pub fn maximise(self, objective: E) -> UnsolvedProblem { self.optimise(ObjectiveDirection::Maximisation, objective) } @@ -518,6 +554,20 @@ impl ProblemVariables { /// # use float_eq::assert_float_eq; /// assert_float_eq!(solution.value(x), -8., abs <= 1e-8); /// ``` + #[cfg(feature = "enable_quadratic")] + pub fn minimise(self, objective: E) -> UnsolvedProblem { + self.optimise_quadratic(ObjectiveDirection::Minimisation, objective) + } + + /// Creates an minimization problem with the given objective. Don't solve it immediately + /// ``` + /// use good_lp::{variables, variable, default_solver, SolverModel, Solution}; + /// variables!{problem: x >= -8;} + /// let solution = problem.minimise(x).using(default_solver).solve().unwrap(); + /// # use float_eq::assert_float_eq; + /// assert_float_eq!(solution.value(x), -8., abs <= 1e-8); + /// ``` + #[cfg(not(feature = "enable_quadratic"))] pub fn minimise(self, objective: E) -> UnsolvedProblem { self.optimise(ObjectiveDirection::Minimisation, objective) } @@ -612,6 +662,12 @@ impl UnsolvedProblem { pub fn using(self, mut solver: S) -> S::Model { solver.create_model(self) } + + /// Check if this problem has quadratic terms + #[cfg(feature = "enable_quadratic")] + pub fn is_quadratic(&self) -> bool { + !self.objective.is_affine() + } } impl> Mul for Variable { @@ -621,6 +677,8 @@ impl> Mul for Variable { let mut coefficients = HashMap::with_capacity_and_hasher(1, Default::default()); coefficients.insert(self, rhs.into()); Expression { + #[cfg(feature = "enable_quadratic")] + quadratic: QuadraticExpression::new(), linear: LinearExpression { coefficients }, constant: 0.0, } @@ -634,6 +692,8 @@ impl Mul for f64 { let mut coefficients = HashMap::with_capacity_and_hasher(1, Default::default()); coefficients.insert(rhs, self); Expression { + #[cfg(feature = "enable_quadratic")] + quadratic: QuadraticExpression::new(), linear: LinearExpression { coefficients }, constant: 0.0, } diff --git a/tests/quadratic_solver_tests.rs b/tests/quadratic_solver_tests.rs index b7703b1..cf1975f 100644 --- a/tests/quadratic_solver_tests.rs +++ b/tests/quadratic_solver_tests.rs @@ -3,10 +3,7 @@ #[cfg(all(feature = "clarabel", feature = "enable_quadratic"))] mod quadratic_integration_tests { - use good_lp::{ - clarabel_quadratic, variables, QuadraticAffineExpression, ResolutionError, Solution, - SolverModel, - }; + use good_lp::{clarabel, variables, Expression, ResolutionError, Solution, SolverModel}; #[test] fn test_constrained_quadratic_with_bounds() { @@ -16,7 +13,7 @@ mod quadratic_integration_tests { } // Minimize x^2 - x - let mut objective = QuadraticAffineExpression::new(); + let mut objective = Expression::default(); objective.add_quadratic_term(x, x, 1.0); objective.add_linear_term(x, -1.0); @@ -25,8 +22,8 @@ mod quadratic_integration_tests { println!("Objective: {:?}", objective); let solution = vars - .minimise_quadratic(objective) - .using(clarabel_quadratic) + .minimise(objective) + .using(clarabel) .solve() .expect("Bounded quadratic should solve"); @@ -46,14 +43,14 @@ mod quadratic_integration_tests { } // Minimize x^2 + y^2 + z^2 subject to x + y + z = 6 and x - y + 2z = 4 - let mut objective = QuadraticAffineExpression::new(); + let mut objective = Expression::default(); objective.add_quadratic_term(x, x, 1.0); objective.add_quadratic_term(y, y, 1.0); objective.add_quadratic_term(z, z, 1.0); let solution = vars - .minimise_quadratic(objective) - .using(clarabel_quadratic) + .minimise(objective) + .using(clarabel) .with((x + y + z).eq(6.0)) // First equality constraint .with((x - y + 2.0 * z).eq(4.0)) // Second equality constraint .solve() @@ -87,6 +84,34 @@ mod quadratic_integration_tests { ); } + #[test] + fn test_with_many_variables() { + let mut objective = Expression::default(); + let mut problem_vars = variables!(); + let mut vars = vec![]; + + for i in 0..1_000 { + let var = problem_vars.add_variable(); + vars.push(var); + + // Build (var - i)^2 = var^2 - 2*i*var + i^2 manually + objective.add_quadratic_term(var, var, 1.0); // var^2 + objective.add_linear_term(var, -2.0 * i as f64); // -2*i*var + objective.add_constant((i as f64) * (i as f64)); // i^2 + } + + let solution = problem_vars + .minimise(objective) + .using(clarabel) + .solve() + .expect("Quadratic with many variables should solve"); + + for i in 0..1_000 { + let var = vars[i]; + assert!((solution.value(var) - i as f64).abs() < 1e-3); + } + } + #[test] fn test_quadratic_with_inequality_constraints() { variables! { @@ -96,13 +121,13 @@ mod quadratic_integration_tests { } // Minimize x^2 + y^2 subject to x + y ≥ 4 - let mut objective = QuadraticAffineExpression::new(); + let mut objective = Expression::default(); objective.add_quadratic_term(x, x, 1.0); objective.add_quadratic_term(y, y, 1.0); let solution = vars - .minimise_quadratic(objective) - .using(clarabel_quadratic) + .minimise(objective) + .using(clarabel) .with((x + y).geq(4.0)) // inequality constraint .solve() .expect("Quadratic with inequality should solve"); @@ -130,7 +155,7 @@ mod quadratic_integration_tests { } // Minimize (x-1)^2 + (y-3)^2 subject to x + y = 5, x ≥ 0 - let mut objective = QuadraticAffineExpression::new(); + let mut objective = Expression::default(); objective.add_quadratic_term(x, x, 1.0); objective.add_quadratic_term(y, y, 1.0); objective.add_linear_term(x, -2.0); @@ -138,8 +163,8 @@ mod quadratic_integration_tests { objective.add_constant(10.0); let solution = vars - .minimise_quadratic(objective) - .using(clarabel_quadratic) + .minimise(objective) + .using(clarabel) .with((x + y).eq(5.0)) // equality .with((1.0 * x).geq(0.0)) // inequality .solve() @@ -164,15 +189,15 @@ mod quadratic_integration_tests { } // Maximize a concave quadratic: -x^2 - y^2 + 8x + 6y - let mut objective = QuadraticAffineExpression::new(); + let mut objective = Expression::default(); objective.add_quadratic_term(x, x, -1.0); // -x^2 (concave) objective.add_quadratic_term(y, y, -1.0); // -y^2 (concave) objective.add_linear_term(x, 8.0); // +8x objective.add_linear_term(y, 6.0); // +6y let solution = vars - .maximise_quadratic(objective) - .using(clarabel_quadratic) + .maximise(objective) + .using(clarabel) .solve() .expect("Quadratic maximization should solve"); @@ -195,14 +220,14 @@ mod quadratic_integration_tests { y; } - let mut objective = QuadraticAffineExpression::new(); + let mut objective = Expression::default(); objective.add_quadratic_term(x, x, 1.0); objective.add_quadratic_term(y, y, 1.0); // Create conflicting constraints: x + y >= 5 and x + y <= 2 let result = vars - .minimise_quadratic(objective) - .using(clarabel_quadratic) + .minimise(objective) + .using(clarabel) .with(x + y >> 5.0) // x + y >= 5 .with(x + y << 2.0) // x + y <= 2 (conflicts with above) .solve(); @@ -226,13 +251,13 @@ mod quadratic_integration_tests { } // Large-scale problem: minimize 1000000*x^2 + 1000000*y^2 - let mut objective = QuadraticAffineExpression::new(); + let mut objective = Expression::default(); objective.add_quadratic_term(x, x, 1_000_000.0); objective.add_quadratic_term(y, y, 1_000_000.0); let solution = vars - .minimise_quadratic(objective) - .using(clarabel_quadratic) + .minimise(objective) + .using(clarabel) .with(x + y >> 1e-3) // Very small constraint .solve() .expect("Scaled quadratic should solve"); diff --git a/tests/variables.rs b/tests/variables.rs index bbb8c0f..beaa2d3 100644 --- a/tests/variables.rs +++ b/tests/variables.rs @@ -33,7 +33,7 @@ fn debug_format() { let mut vars = variables!(); let a = vars.add_variable(); let b = vars.add_variable(); - let expr_str = format!("{:?}", (9 * (1. + a + b / 3)).leq(a + 1)); + let expr_str = format!("{:?}", (9_i32 * (1.0_f64 + a + b / 3.0_f64)).leq(a + 1)); let possibilities = vec!["3 v1 + 8 v0 <= -8", "8 v0 + 3 v1 <= -8"]; assert!( possibilities.contains(&expr_str.as_str()), From 217fffd980e4d7832e235b9b0233dde1560b42e0 Mon Sep 17 00:00:00 2001 From: Shurikal Date: Tue, 30 Sep 2025 07:47:36 +0200 Subject: [PATCH 8/8] Refactor imports and clean up formatting in expression and solver modules --- src/expression.rs | 21 +++++++++------------ src/solvers/clarabel.rs | 16 +++++++--------- src/variable.rs | 2 +- 3 files changed, 17 insertions(+), 22 deletions(-) diff --git a/src/expression.rs b/src/expression.rs index c9915ce..97ce22e 100644 --- a/src/expression.rs +++ b/src/expression.rs @@ -5,9 +5,9 @@ use fnv::FnvHashMap as HashMap; use crate::affine_expression_trait::IntoAffineExpression; use crate::constraint; +use crate::variable::{FormatWithVars, Variable}; #[cfg(feature = "enable_quadratic")] use crate::IntoQuadraticExpression; -use crate::variable::{FormatWithVars, Variable}; use crate::{Constraint, Solution}; /// Represents a pair of variables in a quadratic term @@ -495,7 +495,7 @@ impl FormatWithVars for Expression { FUN: FnMut(&mut Formatter<'_>, Variable) -> std::fmt::Result, { let mut has_terms = false; - + // Format linear terms first for (&var, &coeff) in &self.linear.coefficients { if coeff != 0f64 { @@ -509,7 +509,7 @@ impl FormatWithVars for Expression { has_terms = true; } } - + // Format quadratic terms #[cfg(feature = "enable_quadratic")] { @@ -528,7 +528,7 @@ impl FormatWithVars for Expression { } } } - + // Format constant if self.constant.abs() >= f64::EPSILON { if has_terms { @@ -539,7 +539,7 @@ impl FormatWithVars for Expression { } else if !has_terms { write!(f, "0")?; } - + Ok(()) } } @@ -594,14 +594,14 @@ macro_rules! impl_extended_num_ops { Expression::from(self) + rhs } } - + impl Sub for $num { type Output = Expression; fn sub(self, rhs: Variable) -> Self::Output { Expression::from(self) - rhs } } - + impl Sub for $num { type Output = Expression; fn sub(self, rhs: Expression) -> Self::Output { @@ -656,7 +656,6 @@ impl Sub for Expression { } } - macro_rules! impl_var_num_ops { ($($num:ty),*) => {$( impl Add<$num> for Variable { @@ -870,7 +869,7 @@ mod tests { let mut values = HashMap::new(); values.insert(a, 100); values.insert(b, -1); - assert_eq!((a + 3.0_f64 * (b + 3.0_f64 )).eval_with(&values), 106.) + assert_eq!((a + 3.0_f64 * (b + 3.0_f64)).eval_with(&values), 106.) } #[cfg(feature = "enable_quadratic")] @@ -884,9 +883,7 @@ mod tests { assert_eq!(expression_0_str, "-1 v0 + -1 v1*v1 + 3"); println!("{:?}", expression_0); println!("{:?}", expression_1); - assert_eq!( - expression_0, expression_1 - ); + assert_eq!(expression_0, expression_1); let values = HashMap::from([(v0, 100.), (v1, -1.)]); assert_eq!(expression_0.eval_with(&values), -98.0) } diff --git a/src/solvers/clarabel.rs b/src/solvers/clarabel.rs index 81d2ead..474472b 100644 --- a/src/solvers/clarabel.rs +++ b/src/solvers/clarabel.rs @@ -2,10 +2,10 @@ use crate::affine_expression_trait::IntoAffineExpression; #[cfg(feature = "enable_quadratic")] -use crate::expression::VariablePair; +use crate::expression::Expression; use crate::expression::LinearExpression; #[cfg(feature = "enable_quadratic")] -use crate::expression::Expression; +use crate::expression::VariablePair; use crate::variable::UnsolvedProblem; use crate::{ constraint::ConstraintReference, @@ -257,7 +257,7 @@ impl CscMatrixBuilder { is_quadratic: false, } } - + /// Create a new builder for a square matrix (used for quadratic matrices) #[cfg(feature = "enable_quadratic")] fn new_square(n_vars: usize) -> Self { @@ -269,7 +269,7 @@ impl CscMatrixBuilder { is_quadratic: true, } } - + fn add_row(&mut self, row: LinearExpression) { for (var, value) in row.linear_coefficients() { self.rowval[var.index()].push(self.n_rows); @@ -277,7 +277,7 @@ impl CscMatrixBuilder { } self.n_rows += 1; } - + /// Add a quadratic term to the matrix (for quadratic objectives) #[cfg(feature = "enable_quadratic")] fn add_quadratic_term(&mut self, pair: VariablePair, coeff: f64) { @@ -300,7 +300,7 @@ impl CscMatrixBuilder { self.nzval[i].push(coeff); } } - + fn build(mut self) -> clarabel::algebra::CscMatrix { // For quadratic matrices, sort columns to maintain proper CSC format #[cfg(feature = "enable_quadratic")] @@ -317,7 +317,7 @@ impl CscMatrixBuilder { self.nzval[col] = pairs.iter().map(|&(_, val)| val).collect(); } } - + let mut colptr = Vec::with_capacity(self.n_cols + 1); colptr.push(0); for col in &self.rowval { @@ -351,8 +351,6 @@ fn fast_flatten_vecs(vecs: Vec>) -> Vec { result } - - #[cfg(test)] mod tests { diff --git a/src/variable.rs b/src/variable.rs index 937ef93..02f9b31 100644 --- a/src/variable.rs +++ b/src/variable.rs @@ -508,7 +508,7 @@ impl ProblemVariables { pub fn optimise_quadratic( self, direction: ObjectiveDirection, - objective: E + objective: E, ) -> UnsolvedProblem { let objective = Expression::from_other_quadratic(objective); // todo, correct assert check