Skip to main content

otspot_core/qp/
problem.rs

1//! QP問題のデータ構造定義。
2
3use crate::problem::{is_valid_bound_pair, ConstraintType, SolveStatus};
4use crate::sparse::CscMatrix;
5
6/// Quadratic term storage for a single QCQP constraint.
7///
8/// Uses COO (coordinate) format: symmetrized `(row, col, val)` triplets.
9/// Memory is O(nnz), unlike `CscMatrix` which requires O(n) for `col_ptr`.
10#[derive(Debug, Clone, Default)]
11pub struct QcqpMatrix {
12    /// Problem dimension (n×n symmetric matrix).
13    pub n: usize,
14    /// Symmetrized COO triplets: (row, col, value), 0-indexed.
15    pub triplets: Vec<(usize, usize, f64)>,
16}
17
18impl QcqpMatrix {
19    /// Creates an empty n×n matrix with no non-zero entries.
20    pub fn new(n: usize) -> Self {
21        Self {
22            n,
23            triplets: Vec::new(),
24        }
25    }
26
27    /// Returns the number of stored entries (after symmetrization).
28    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    /// Non-finite coefficient (NaN or ±∞) in the named field at the given index.
38    NonFiniteCoefficient {
39        field: &'static str,
40        index: usize,
41    },
42    /// Invalid variable bound: NaN or lb > ub at the given index.
43    InvalidBounds {
44        index: usize,
45        lb: f64,
46        ub: f64,
47    },
48    /// A triplet (row, col) index exceeds the matrix dimension.
49    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/// min 1/2 x^T Q x + c^T x  s.t. Ax {<=,=,>=} b, lb <= x <= ub
90///
91/// When `quadratic_constraints` is non-empty the problem is a QCQP.
92/// Entry `k` holds the symmetric `n×n` matrix `Q_k` for the quadratic part
93/// of constraint `k`: `1/2 x^T Q_k x + a_k^T x {<=,=,>=} b_k`.
94/// An empty `QcqpMatrix` at index `k` means constraint `k` has no quadratic part.
95#[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    /// Per-constraint quadratic matrices for QCQP, stored in COO format.
106    ///
107    /// Length is either 0 (pure QP/LP) or `num_constraints` (QCQP).
108    /// Entry `k` holds symmetrized triplets for `Q_k`; empty means no quadratic part.
109    pub quadratic_constraints: Vec<QcqpMatrix>,
110    /// 目的関数値 = 1/2 x^T Q x + c^T x + obj_offset
111    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    /// Convenience constructor: all constraints are `≤`.
203    #[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    /// Q が構造的にゼロ(LP)かを判定する。
216    ///
217    /// LP/QP の dispatch 判定に使う。数値閾値ではなく**構造的ゼロ**で判定する:
218    /// stored 値が一つでも非ゼロなら QP (IPM 経路)。`CscMatrix::from_triplets` は
219    /// `|v| ≤ DROP_TOL` を構築時に落とすため、ここに残る値は構造的非ゼロである。
220    /// 閾値判定 (例 `|v| < 1e-12`) は微小 Q QP を LP 化して status を変える
221    /// (例 bounded QP → false-Unbounded)。dispatch は status を変えてはならない。
222    pub fn is_zero_q(&self) -> bool {
223        self.q.values.iter().all(|&v| v == 0.0)
224    }
225
226    /// Returns `true` if the problem has at least one constraint with a non-zero quadratic term.
227    ///
228    /// Used as a guard: the QP/LP solver cannot handle QCQP and must reject such problems.
229    pub fn has_qcqp_constraints(&self) -> bool {
230        self.quadratic_constraints
231            .iter()
232            .any(|q| !q.triplets.is_empty())
233    }
234
235    /// Set per-constraint quadratic matrices for QCQP, with validation.
236    ///
237    /// `qcs` must be either empty (pure QP) or have length equal to `num_constraints`.
238    /// Each `QcqpMatrix` must have `n == num_vars`, finite values, and triplet indices
239    /// within `[0, n)`.
240    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    /// Q が対角行列かどうかを検査する
280    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    // --- NaN / ±inf in c ---
394    #[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    // --- NaN / ±inf in b ---
416    #[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    // --- NaN / ±inf in Q values ---
438    // Note: from_triplets drops NaN via DROP_TOL (NaN.abs() > tol == false).
439    // Inject bad values directly to test validation independent of that path.
440    #[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; // inject bad value directly
447            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    // --- NaN / ±inf in A values ---
464    #[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; // inject bad value directly
472            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    // --- NaN in bounds ---
488    #[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    // --- lb > ub ---
507    #[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    // --- inf bounds are valid ---
531    #[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    // Bug C/D: (inf,inf) and (-inf,-inf) form empty intervals that lb>ub misses
544    // because inf==inf and -inf==-inf. These must be rejected.
545    #[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    // Sentinel: valid half-infinite and fixed-var bounds must be accepted.
567    #[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    // Sentinel: lb=+inf (alone) and ub=-inf (alone) must each be rejected.
588    #[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    // --- dimension mismatch still caught ---
612    #[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]; // wrong length
618        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    // --- set_quadratic_constraints ---
626
627    #[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        // num_constraints = 1, so length 2 is invalid
638        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    // --- sentinel: n mismatch rejected (no-op rewrite → FAIL) ---
701    #[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        // num_vars=2, but QcqpMatrix.n=3 — must be rejected
712        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    // --- sentinel: OOB triplet row index rejected (no-op rewrite → FAIL) ---
722    #[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)); // row=2 >= n=2
734        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    // --- sentinel: OOB triplet col index rejected (no-op rewrite → FAIL) ---
750    #[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)); // col=5 >= n=2
762        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    // --- sentinel: CscMatrix accessor returns correct data (no-op rewrite → FAIL) ---
773    #[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}