1use crate::problem::{is_valid_bound_pair, ConstraintType, SolveStatus};
4use crate::sparse::CscMatrix;
5
6#[derive(Debug, Clone, Default)]
11pub struct QcqpMatrix {
12 pub n: usize,
14 pub triplets: Vec<(usize, usize, f64)>,
16}
17
18impl QcqpMatrix {
19 pub fn new(n: usize) -> Self {
21 Self {
22 n,
23 triplets: Vec::new(),
24 }
25 }
26
27 pub fn nnz(&self) -> usize {
29 self.triplets.len()
30 }
31}
32
33#[non_exhaustive]
34#[derive(Debug)]
35pub enum QpProblemError {
36 DimensionMismatch(String),
37 NonFiniteCoefficient {
39 field: &'static str,
40 index: usize,
41 },
42 InvalidBounds {
44 index: usize,
45 lb: f64,
46 ub: f64,
47 },
48 TripletIndexOutOfBounds {
50 constraint: usize,
51 row: usize,
52 col: usize,
53 n: usize,
54 },
55}
56
57impl std::fmt::Display for QpProblemError {
58 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
59 match self {
60 QpProblemError::DimensionMismatch(msg) => write!(f, "dimension mismatch: {}", msg),
61 QpProblemError::NonFiniteCoefficient { field, index } => {
62 write!(f, "non-finite coefficient in {}: index {}", field, index)
63 }
64 QpProblemError::InvalidBounds { index, lb, ub } => {
65 write!(
66 f,
67 "invalid bounds at index {}: lb={} > ub={} or NaN",
68 index, lb, ub
69 )
70 }
71 QpProblemError::TripletIndexOutOfBounds {
72 constraint,
73 row,
74 col,
75 n,
76 } => {
77 write!(
78 f,
79 "triplet index out of bounds in constraint {}: ({},{}) >= n={}",
80 constraint, row, col, n
81 )
82 }
83 }
84 }
85}
86
87impl std::error::Error for QpProblemError {}
88
89#[derive(Debug, Clone)]
96pub struct QpProblem {
97 pub q: CscMatrix,
98 pub c: Vec<f64>,
99 pub a: CscMatrix,
100 pub b: Vec<f64>,
101 pub bounds: Vec<(f64, f64)>,
102 pub num_vars: usize,
103 pub num_constraints: usize,
104 pub constraint_types: Vec<ConstraintType>,
105 pub quadratic_constraints: Vec<QcqpMatrix>,
110 pub obj_offset: f64,
112}
113
114impl QpProblem {
115 pub fn new(
116 q: CscMatrix,
117 c: Vec<f64>,
118 a: CscMatrix,
119 b: Vec<f64>,
120 bounds: Vec<(f64, f64)>,
121 constraint_types: Vec<ConstraintType>,
122 ) -> Result<Self, QpProblemError> {
123 let n = c.len();
124 let m = b.len();
125 if q.nrows != n || q.ncols != n {
126 return Err(QpProblemError::DimensionMismatch(format!(
127 "Q must be {}x{}, got {}x{}",
128 n, n, q.nrows, q.ncols
129 )));
130 }
131 if a.nrows != m || a.ncols != n {
132 return Err(QpProblemError::DimensionMismatch(format!(
133 "A must be {}x{}, got {}x{}",
134 m, n, a.nrows, a.ncols
135 )));
136 }
137 if bounds.len() != n {
138 return Err(QpProblemError::DimensionMismatch(format!(
139 "bounds length must be {}, got {}",
140 n,
141 bounds.len()
142 )));
143 }
144 if constraint_types.len() != m {
145 return Err(QpProblemError::DimensionMismatch(format!(
146 "constraint_types length must be {}, got {}",
147 m,
148 constraint_types.len()
149 )));
150 }
151 for (i, &v) in c.iter().enumerate() {
152 if !v.is_finite() {
153 return Err(QpProblemError::NonFiniteCoefficient {
154 field: "c",
155 index: i,
156 });
157 }
158 }
159 for (i, &v) in b.iter().enumerate() {
160 if !v.is_finite() {
161 return Err(QpProblemError::NonFiniteCoefficient {
162 field: "b",
163 index: i,
164 });
165 }
166 }
167 for (i, &v) in q.values.iter().enumerate() {
168 if !v.is_finite() {
169 return Err(QpProblemError::NonFiniteCoefficient {
170 field: "Q",
171 index: i,
172 });
173 }
174 }
175 for (i, &v) in a.values.iter().enumerate() {
176 if !v.is_finite() {
177 return Err(QpProblemError::NonFiniteCoefficient {
178 field: "A",
179 index: i,
180 });
181 }
182 }
183 for (i, &(lb, ub)) in bounds.iter().enumerate() {
184 if !is_valid_bound_pair(lb, ub) {
185 return Err(QpProblemError::InvalidBounds { index: i, lb, ub });
186 }
187 }
188 Ok(QpProblem {
189 q,
190 c,
191 a,
192 b,
193 bounds,
194 num_vars: n,
195 num_constraints: m,
196 constraint_types,
197 quadratic_constraints: vec![],
198 obj_offset: 0.0,
199 })
200 }
201
202 #[doc(hidden)]
204 pub fn new_all_le(
205 q: CscMatrix,
206 c: Vec<f64>,
207 a: CscMatrix,
208 b: Vec<f64>,
209 bounds: Vec<(f64, f64)>,
210 ) -> Result<Self, QpProblemError> {
211 let m = b.len();
212 Self::new(q, c, a, b, bounds, vec![ConstraintType::Le; m])
213 }
214
215 pub fn is_zero_q(&self) -> bool {
223 self.q.values.iter().all(|&v| v == 0.0)
224 }
225
226 pub fn has_qcqp_constraints(&self) -> bool {
230 self.quadratic_constraints
231 .iter()
232 .any(|q| !q.triplets.is_empty())
233 }
234
235 pub fn set_quadratic_constraints(
241 &mut self,
242 qcs: Vec<QcqpMatrix>,
243 ) -> Result<(), QpProblemError> {
244 if !qcs.is_empty() && qcs.len() != self.num_constraints {
245 return Err(QpProblemError::DimensionMismatch(format!(
246 "quadratic_constraints length must be 0 or {}, got {}",
247 self.num_constraints,
248 qcs.len()
249 )));
250 }
251 for (k, qc) in qcs.iter().enumerate() {
252 if qc.n != self.num_vars {
253 return Err(QpProblemError::DimensionMismatch(format!(
254 "quadratic_constraints[{}].n must be {}, got {}",
255 k, self.num_vars, qc.n
256 )));
257 }
258 for &(row, col, v) in &qc.triplets {
259 if !v.is_finite() {
260 return Err(QpProblemError::NonFiniteCoefficient {
261 field: "quadratic_constraints",
262 index: k,
263 });
264 }
265 if row >= qc.n || col >= qc.n {
266 return Err(QpProblemError::TripletIndexOutOfBounds {
267 constraint: k,
268 row,
269 col,
270 n: qc.n,
271 });
272 }
273 }
274 }
275 self.quadratic_constraints = qcs;
276 Ok(())
277 }
278
279 pub fn is_diagonal_q(&self) -> bool {
281 for col in 0..self.num_vars {
282 let start = self.q.col_ptr[col];
283 let end = self.q.col_ptr[col + 1];
284 for k in start..end {
285 let row = self.q.row_ind[k];
286 if row != col && self.q.values[k].abs() > 1e-12 {
287 return false;
288 }
289 }
290 }
291 true
292 }
293}
294
295impl crate::problem::SolverResult {
296 pub fn infeasible() -> Self {
297 crate::problem::SolverResult {
298 status: SolveStatus::Infeasible,
299 objective: f64::INFINITY,
300 solution: vec![],
301 dual_solution: vec![],
302 bound_duals: vec![],
303
304 iterations: 0,
305 ..Default::default()
306 }
307 }
308
309 pub fn unbounded() -> Self {
310 crate::problem::SolverResult {
311 status: SolveStatus::Unbounded,
312 objective: f64::NEG_INFINITY,
313 solution: vec![],
314 dual_solution: vec![],
315 bound_duals: vec![],
316 iterations: 0,
317 ..Default::default()
318 }
319 }
320
321 pub fn numerical_error() -> Self {
322 crate::problem::SolverResult {
323 status: SolveStatus::NumericalError,
324 objective: f64::INFINITY,
325 solution: vec![],
326 dual_solution: vec![],
327 bound_duals: vec![],
328
329 iterations: 0,
330 ..Default::default()
331 }
332 }
333
334 pub fn not_supported(msg: impl Into<String>) -> Self {
335 crate::problem::SolverResult {
336 status: SolveStatus::NotSupported(msg.into()),
337 objective: f64::INFINITY,
338 solution: vec![],
339 dual_solution: vec![],
340 bound_duals: vec![],
341 iterations: 0,
342 ..Default::default()
343 }
344 }
345}
346
347pub use crate::options::QpWarmStart;
348
349#[cfg(test)]
350mod tests {
351 use super::*;
352 use crate::problem::ConstraintType;
353 use crate::sparse::CscMatrix;
354
355 fn make_qp(
356 c: Vec<f64>,
357 b: Vec<f64>,
358 q_vals: Vec<f64>,
359 a_vals: Vec<f64>,
360 bounds: Vec<(f64, f64)>,
361 ) -> Result<QpProblem, QpProblemError> {
362 let n = c.len();
363 let m = b.len();
364 let q = if q_vals.is_empty() {
365 CscMatrix::new(n, n)
366 } else {
367 let idx: Vec<usize> = (0..n).collect();
368 crate::sparse::CscMatrix::from_triplets(&idx, &idx, &q_vals, n, n).unwrap()
369 };
370 let a = if a_vals.is_empty() {
371 CscMatrix::new(m, n)
372 } else {
373 let rows = vec![0usize; n];
374 let cols: Vec<usize> = (0..n).collect();
375 crate::sparse::CscMatrix::from_triplets(&rows, &cols, &a_vals, m, n).unwrap()
376 };
377 let ct = vec![ConstraintType::Le; m];
378 QpProblem::new(q, c, a, b, bounds, ct)
379 }
380
381 #[test]
382 fn valid_qp_accepted() {
383 let res = make_qp(
384 vec![1.0, 2.0],
385 vec![5.0],
386 vec![],
387 vec![1.0, 1.0],
388 vec![(0.0, f64::INFINITY), (0.0, 10.0)],
389 );
390 assert!(res.is_ok());
391 }
392
393 #[test]
395 fn nan_in_c_rejected() {
396 let cases = [f64::NAN, f64::INFINITY, f64::NEG_INFINITY];
397 for bad in cases {
398 let res = make_qp(
399 vec![bad, 1.0],
400 vec![5.0],
401 vec![],
402 vec![1.0, 1.0],
403 vec![(0.0, f64::INFINITY); 2],
404 );
405 assert!(
406 matches!(
407 res,
408 Err(QpProblemError::NonFiniteCoefficient { field: "c", .. })
409 ),
410 "expected NonFiniteCoefficient for c={bad}"
411 );
412 }
413 }
414
415 #[test]
417 fn nan_in_b_rejected() {
418 let cases = [f64::NAN, f64::INFINITY, f64::NEG_INFINITY];
419 for bad in cases {
420 let res = make_qp(
421 vec![1.0, 2.0],
422 vec![bad],
423 vec![],
424 vec![1.0, 1.0],
425 vec![(0.0, f64::INFINITY); 2],
426 );
427 assert!(
428 matches!(
429 res,
430 Err(QpProblemError::NonFiniteCoefficient { field: "b", .. })
431 ),
432 "expected NonFiniteCoefficient for b={bad}"
433 );
434 }
435 }
436
437 #[test]
441 fn nan_in_q_rejected() {
442 let n = 2;
443 let bad_vals = [f64::NAN, f64::INFINITY, f64::NEG_INFINITY];
444 for bad in bad_vals {
445 let mut q = CscMatrix::from_triplets(&[0], &[0], &[1.0], n, n).unwrap();
446 q.values[0] = bad; let a = CscMatrix::new(1, n);
448 let c = vec![1.0, 2.0];
449 let b = vec![5.0];
450 let bounds = vec![(0.0, f64::INFINITY); n];
451 let ct = vec![ConstraintType::Le];
452 let res = QpProblem::new(q, c, a, b, bounds, ct);
453 assert!(
454 matches!(
455 res,
456 Err(QpProblemError::NonFiniteCoefficient { field: "Q", .. })
457 ),
458 "expected NonFiniteCoefficient for Q val={bad}"
459 );
460 }
461 }
462
463 #[test]
465 fn nan_in_a_rejected() {
466 let n = 2;
467 let bad_vals = [f64::NAN, f64::INFINITY, f64::NEG_INFINITY];
468 for bad in bad_vals {
469 let q = CscMatrix::new(n, n);
470 let mut a = CscMatrix::from_triplets(&[0], &[0], &[1.0], 1, n).unwrap();
471 a.values[0] = bad; let c = vec![1.0, 2.0];
473 let b = vec![5.0];
474 let bounds = vec![(0.0, f64::INFINITY); n];
475 let ct = vec![ConstraintType::Le];
476 let res = QpProblem::new(q, c, a, b, bounds, ct);
477 assert!(
478 matches!(
479 res,
480 Err(QpProblemError::NonFiniteCoefficient { field: "A", .. })
481 ),
482 "expected NonFiniteCoefficient for A val={bad}"
483 );
484 }
485 }
486
487 #[test]
489 fn nan_in_bounds_rejected() {
490 let cases: Vec<(f64, f64)> = vec![(f64::NAN, 1.0), (0.0, f64::NAN), (f64::NAN, f64::NAN)];
491 for (lb, ub) in cases {
492 let res = make_qp(
493 vec![1.0, 2.0],
494 vec![5.0],
495 vec![],
496 vec![1.0, 1.0],
497 vec![(lb, ub), (0.0, f64::INFINITY)],
498 );
499 assert!(
500 matches!(res, Err(QpProblemError::InvalidBounds { index: 0, .. })),
501 "expected InvalidBounds for ({lb},{ub})"
502 );
503 }
504 }
505
506 #[test]
508 fn lb_gt_ub_rejected() {
509 let cases: Vec<(f64, f64)> = vec![
510 (5.0, 1.0),
511 (1.0, 0.0),
512 (f64::INFINITY, f64::NEG_INFINITY),
513 (0.1, 0.0),
514 ];
515 for (lb, ub) in cases {
516 let res = make_qp(
517 vec![1.0, 2.0],
518 vec![5.0],
519 vec![],
520 vec![1.0, 1.0],
521 vec![(lb, ub), (0.0, f64::INFINITY)],
522 );
523 assert!(
524 matches!(res, Err(QpProblemError::InvalidBounds { .. })),
525 "expected InvalidBounds for lb={lb} ub={ub}"
526 );
527 }
528 }
529
530 #[test]
532 fn inf_bounds_accepted() {
533 let res = make_qp(
534 vec![1.0, 2.0],
535 vec![5.0],
536 vec![],
537 vec![1.0, 1.0],
538 vec![(f64::NEG_INFINITY, f64::INFINITY), (0.0, f64::INFINITY)],
539 );
540 assert!(res.is_ok(), "±inf bounds should be valid");
541 }
542
543 #[test]
546 fn qp_empty_interval_same_inf_rejected() {
547 let cases: Vec<(f64, f64)> = vec![
548 (f64::INFINITY, f64::INFINITY),
549 (f64::NEG_INFINITY, f64::NEG_INFINITY),
550 ];
551 for (lb, ub) in cases {
552 let res = make_qp(
553 vec![1.0, 2.0],
554 vec![5.0],
555 vec![],
556 vec![1.0, 1.0],
557 vec![(lb, ub), (0.0, f64::INFINITY)],
558 );
559 assert!(
560 matches!(res, Err(QpProblemError::InvalidBounds { index: 0, .. })),
561 "expected InvalidBounds for ({lb},{ub})"
562 );
563 }
564 }
565
566 #[test]
568 fn qp_valid_bound_variants_accepted() {
569 let valid_cases: Vec<Vec<(f64, f64)>> = vec![
570 vec![(f64::NEG_INFINITY, 5.0), (0.0, f64::INFINITY)],
571 vec![(0.0, f64::INFINITY), (f64::NEG_INFINITY, f64::INFINITY)],
572 vec![(3.0, 3.0), (0.0, 10.0)],
573 vec![(0.0, 0.0), (0.0, 0.0)],
574 ];
575 for bounds in valid_cases {
576 let res = make_qp(
577 vec![1.0, 2.0],
578 vec![5.0],
579 vec![],
580 vec![1.0, 1.0],
581 bounds.clone(),
582 );
583 assert!(res.is_ok(), "expected Ok for bounds={bounds:?}");
584 }
585 }
586
587 #[test]
589 fn qp_lone_inf_bound_rejected() {
590 let cases: Vec<(f64, f64)> = vec![
591 (f64::INFINITY, 5.0),
592 (f64::INFINITY, f64::INFINITY),
593 (0.0, f64::NEG_INFINITY),
594 (f64::NEG_INFINITY, f64::NEG_INFINITY),
595 ];
596 for (lb, ub) in cases {
597 let res = make_qp(
598 vec![1.0, 2.0],
599 vec![5.0],
600 vec![],
601 vec![1.0, 1.0],
602 vec![(lb, ub), (0.0, f64::INFINITY)],
603 );
604 assert!(
605 matches!(res, Err(QpProblemError::InvalidBounds { index: 0, .. })),
606 "expected InvalidBounds for ({lb},{ub})"
607 );
608 }
609 }
610
611 #[test]
613 fn dimension_mismatch_still_caught() {
614 let n = 2;
615 let q = CscMatrix::new(n, n);
616 let a = CscMatrix::new(1, n);
617 let c = vec![1.0, 2.0, 3.0]; let b = vec![5.0];
619 let bounds = vec![(0.0, f64::INFINITY); n];
620 let ct = vec![ConstraintType::Le];
621 let res = QpProblem::new(q, c, a, b, bounds, ct);
622 assert!(matches!(res, Err(QpProblemError::DimensionMismatch(_))));
623 }
624
625 #[test]
628 fn set_quadratic_constraints_length_mismatch_rejected() {
629 let mut prob = make_qp(
630 vec![1.0, 2.0],
631 vec![5.0],
632 vec![],
633 vec![1.0, 1.0],
634 vec![(0.0, f64::INFINITY); 2],
635 )
636 .unwrap();
637 let qcs = vec![QcqpMatrix::new(2), QcqpMatrix::new(2)];
639 let res = prob.set_quadratic_constraints(qcs);
640 assert!(matches!(res, Err(QpProblemError::DimensionMismatch(_))));
641 }
642
643 #[test]
644 fn set_quadratic_constraints_nan_rejected() {
645 let mut prob = make_qp(
646 vec![1.0, 2.0],
647 vec![5.0],
648 vec![],
649 vec![1.0, 1.0],
650 vec![(0.0, f64::INFINITY); 2],
651 )
652 .unwrap();
653 let bad_vals = [f64::NAN, f64::INFINITY, f64::NEG_INFINITY];
654 for bad in bad_vals {
655 let mut qc = QcqpMatrix::new(2);
656 qc.triplets.push((0, 0, bad));
657 let res = prob.set_quadratic_constraints(vec![qc]);
658 assert!(
659 matches!(
660 res,
661 Err(QpProblemError::NonFiniteCoefficient {
662 field: "quadratic_constraints",
663 ..
664 })
665 ),
666 "expected NonFiniteCoefficient for val={bad}"
667 );
668 }
669 }
670
671 #[test]
672 fn set_quadratic_constraints_empty_accepted() {
673 let mut prob = make_qp(
674 vec![1.0, 2.0],
675 vec![5.0],
676 vec![],
677 vec![1.0, 1.0],
678 vec![(0.0, f64::INFINITY); 2],
679 )
680 .unwrap();
681 assert!(prob.set_quadratic_constraints(vec![]).is_ok());
682 }
683
684 #[test]
685 fn set_quadratic_constraints_valid_accepted() {
686 let mut prob = make_qp(
687 vec![1.0, 2.0],
688 vec![5.0],
689 vec![],
690 vec![1.0, 1.0],
691 vec![(0.0, f64::INFINITY); 2],
692 )
693 .unwrap();
694 let mut qc = QcqpMatrix::new(2);
695 qc.triplets.push((0, 0, 2.0));
696 qc.triplets.push((1, 1, 3.0));
697 assert!(prob.set_quadratic_constraints(vec![qc]).is_ok());
698 }
699
700 #[test]
702 fn set_quadratic_constraints_n_mismatch_rejected() {
703 let mut prob = make_qp(
704 vec![1.0, 2.0],
705 vec![5.0],
706 vec![],
707 vec![1.0, 1.0],
708 vec![(0.0, f64::INFINITY); 2],
709 )
710 .unwrap();
711 let mut qc = QcqpMatrix::new(3);
713 qc.triplets.push((0, 0, 1.0));
714 let res = prob.set_quadratic_constraints(vec![qc]);
715 assert!(
716 matches!(res, Err(QpProblemError::DimensionMismatch(_))),
717 "n mismatch must be rejected"
718 );
719 }
720
721 #[test]
723 fn set_quadratic_constraints_triplet_row_oob_rejected() {
724 let mut prob = make_qp(
725 vec![1.0, 2.0],
726 vec![5.0],
727 vec![],
728 vec![1.0, 1.0],
729 vec![(0.0, f64::INFINITY); 2],
730 )
731 .unwrap();
732 let mut qc = QcqpMatrix::new(2);
733 qc.triplets.push((2, 0, 1.0)); let res = prob.set_quadratic_constraints(vec![qc]);
735 assert!(
736 matches!(
737 res,
738 Err(QpProblemError::TripletIndexOutOfBounds {
739 row: 2,
740 col: 0,
741 n: 2,
742 ..
743 })
744 ),
745 "row OOB must be rejected"
746 );
747 }
748
749 #[test]
751 fn set_quadratic_constraints_triplet_col_oob_rejected() {
752 let mut prob = make_qp(
753 vec![1.0, 2.0],
754 vec![5.0],
755 vec![],
756 vec![1.0, 1.0],
757 vec![(0.0, f64::INFINITY); 2],
758 )
759 .unwrap();
760 let mut qc = QcqpMatrix::new(2);
761 qc.triplets.push((0, 5, 1.0)); let res = prob.set_quadratic_constraints(vec![qc]);
763 assert!(
764 matches!(
765 res,
766 Err(QpProblemError::TripletIndexOutOfBounds { col: 5, n: 2, .. })
767 ),
768 "col OOB must be rejected"
769 );
770 }
771
772 #[test]
774 fn csc_accessor_returns_correct_data() {
775 let rows = vec![0usize, 1];
776 let cols = vec![0usize, 1];
777 let vals = vec![3.0f64, 7.0];
778 let m = CscMatrix::from_triplets(&rows, &cols, &vals, 2, 2).unwrap();
779 assert_eq!(m.nrows(), 2);
780 assert_eq!(m.ncols(), 2);
781 assert_eq!(m.nnz(), 2);
782 assert_eq!(m.col_ptr().len(), 3);
783 assert_eq!(m.row_ind().len(), 2);
784 assert_eq!(m.values().len(), 2);
785 let (ri, vs) = m.get_column(0).unwrap();
786 assert_eq!(ri, &[0]);
787 assert_eq!(vs, &[3.0]);
788 let (ri, vs) = m.get_column(1).unwrap();
789 assert_eq!(ri, &[1]);
790 assert_eq!(vs, &[7.0]);
791 }
792}