otspot_core/problem/
mod.rs

1//! LP問題定義モジュール
2//!
3//! 線形計画問題（LP）の構造定義・制約種別・ソルバー結果の表現を提供する。
4//! 問題は標準形 `min c^T x  s.t.  Ax {<=,>=,=} b,  x in [lb, ub]` で定義される。
5
6pub mod certificate;
7use certificate::{BoundGapCertificate, OptimalCertificate};
8
9use crate::error::SolverError;
10use crate::options::WarmStartBasis;
11use crate::sparse::CscMatrix;
12use std::fmt;
13
14/// Returns `true` if `(lb, ub)` is a valid bound pair.
15///
16/// Valid means: neither endpoint is NaN, `lb <= ub`, `lb < +∞`, and `ub > -∞`.
17/// The last two conditions reject empty intervals that `lb > ub` misses when both
18/// endpoints are the same infinity — e.g., `(+∞, +∞)` or `(-∞, -∞)`.
19pub(crate) fn is_valid_bound_pair(lb: f64, ub: f64) -> bool {
20    !lb.is_nan() && !ub.is_nan() && lb <= ub && lb < f64::INFINITY && ub > f64::NEG_INFINITY
21}
22
23/// LP問題における制約条件の種別
24#[non_exhaustive]
25#[derive(Debug, Clone, Copy, PartialEq)]
26pub enum ConstraintType {
27    /// 以下（<=）
28    Le,
29    /// 以上（>=）
30    Ge,
31    /// 等式（==）
32    Eq,
33}
34
35/// Route taken by a solve call (populated per-result, race-free).
36#[non_exhaustive]
37#[derive(Debug, Clone, Copy, PartialEq, Default)]
38pub enum SolveRoute {
39    /// Route not yet set (default for uninitialized results).
40    #[default]
41    Unknown,
42    /// Called directly via `crate::lp::solve_lp_with`.
43    LpDirect,
44    /// LP forwarded from `solve_qp_with(Q=0)`.
45    LpForwardedFromQp,
46    /// QP solved via IPM (Q≠0).
47    QpIpm,
48}
49
50/// Per-solve routing and warm-start statistics (race-free, per-result).
51///
52/// Replaces process-global `AtomicU64` counters so parallel tests observe
53/// independent stats without reset/race issues.
54#[derive(Debug, Clone, Default)]
55pub struct SolveStats {
56    /// Route taken for this solve.
57    pub route: SolveRoute,
58    /// Whether the solver stopped because the deadline (timeout_secs / deadline) was reached.
59    ///
60    /// `true` iff `result.status == SolveStatus::Timeout`. Deterministic sentinel for
61    /// deadline-enforcement tests: assert this field instead of measuring wall time.
62    pub deadline_triggered: bool,
63    /// Whether the postsolve saddle-point Krylov IR was skipped because the solution
64    /// already met the user tolerance (`kkt_already_pass`). Deterministic sentinel for
65    /// the gate: removing the gate (always refine) flips this to `false`.
66    pub postsolve_krylov_ir_skipped: bool,
67    /// For LP solves via the QP dispatch (`LpForwardedFromQp`): `true` if the
68    /// returned result came from the IPM, `false` if from simplex. Lets callers
69    /// (e.g. benchmarks) label the actual route instead of a static size guess.
70    pub lp_ipm_path: bool,
71}
72
73/// ソルバーの求解結果ステータス
74#[non_exhaustive]
75#[derive(Debug, Clone, PartialEq)]
76pub enum SolveStatus {
77    /// 最適解が求まった
78    Optimal,
79    /// 局所的最適解（非凸QP: 慣性修正付きIPMが収束したKKT点）
80    ///
81    /// Q行列が不定（indefinite）の場合、大域的最適性は保証されないが、
82    /// KKT条件を満たす局所最適解またはサドル点が返される。
83    /// 慣性修正（Gershgorin 評価から導出した δI 加算）により IPM を収束させた。
84    LocallyOptimal,
85    /// 問題が実行不可能（infeasible）
86    Infeasible,
87    /// 問題が非有界（unbounded）
88    Unbounded,
89    /// 反復回数上限に到達した（最適性未確認）
90    MaxIterations,
91    /// 解は見つかったが精度基準未達（偽Optimal検出: スケール解除後の残差超過）
92    SuboptimalSolution,
93    /// タイムアウト（timeout_secs を超過した）
94    Timeout,
95    /// 数値エラー（LDL分解失敗等、問題が数値的に解けない）
96    NumericalError,
97    /// Q行列が不定（非凸QP）。IPMはQ正半定値を前提とする。
98    NonConvex(String),
99    /// 非凸 QP の局所最適解 (= `solve_qp_global` 経由で incumbent あり、ε-global 証明なし)。
100    ///
101    /// BB driver が deadline / max_nodes / max_depth で打ち切られ、incumbent ある状態。
102    /// `LocallyOptimal` (= IPM inertia 補正後の単発解) と区別して、caller が「探索打切」
103    /// vs「単発 KKT 収束」を識別できる。`Optimal` には**含めない** (= global proof なし)。
104    NonconvexLocal,
105    /// 非凸 QP の大域 ε-最適解 (= `solve_qp_global` で gap_tol まで証明済み + Q が indefinite)。
106    ///
107    /// `Optimal` は「Q が PSD で IPM/BB が global 達成」専用に維持し、indefinite Q の場合は
108    /// 本 variant で明示分離する (caller が「global 証明済」かを fact で判別)。
109    NonconvexGlobal,
110    /// Problem type not supported by this solver.
111    ///
112    /// Returned when the caller passes a problem that cannot be handled, e.g.
113    /// a QCQP (quadratic constraints present) submitted to the QP/LP entry.
114    NotSupported(String),
115}
116
117impl fmt::Display for SolveStatus {
118    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
119        match self {
120            SolveStatus::Optimal => write!(f, "Optimal"),
121            SolveStatus::LocallyOptimal => write!(f, "LocallyOptimal"),
122            SolveStatus::Infeasible => write!(f, "Infeasible"),
123            SolveStatus::Unbounded => write!(f, "Unbounded"),
124            SolveStatus::MaxIterations => write!(f, "MaxIterations"),
125            SolveStatus::SuboptimalSolution => write!(f, "SuboptimalSolution"),
126            SolveStatus::Timeout => write!(f, "Timeout"),
127            SolveStatus::NumericalError => write!(f, "NumericalError"),
128            SolveStatus::NonConvex(msg) => write!(f, "NonConvex({})", msg),
129            SolveStatus::NonconvexLocal => write!(f, "NonconvexLocal"),
130            SolveStatus::NonconvexGlobal => write!(f, "NonconvexGlobal"),
131            SolveStatus::NotSupported(msg) => write!(f, "NotSupported({})", msg),
132        }
133    }
134}
135
136/// LP/QP共通求解結果型
137///
138/// LP求解（Simplex等）と QP求解（AS/IPM/Concurrent）の両方で使用できる統一結果型。
139/// LP固有フィールド（`reduced_costs`, `slack`, `warm_start_basis`）は QP求解時は空/None。
140/// QP固有フィールド（`bound_duals`, `iterations`）は LP求解時は空/0。
141#[derive(Debug, Clone)]
142pub struct SolverResult {
143    /// 求解ステータス
144    pub status: SolveStatus,
145    /// 最適目的関数値（最適解が存在する場合）
146    pub objective: f64,
147    /// 解ベクトル（最適解が存在する場合）
148    pub solution: Vec<f64>,
149    /// 双対変数ベクトル（各制約の影価格、最適解が存在する場合）
150    pub dual_solution: Vec<f64>,
151    // --- LP固有フィールド ---
152    /// 被縮小費用ベクトル（各決定変数に対して、最適解が存在する場合）
153    pub reduced_costs: Vec<f64>,
154    /// スラック変数ベクトル（各制約のスラック b_i - a_i^T x、最適解が存在する場合）
155    pub slack: Vec<f64>,
156    /// warm-start用の基底情報（Optimal時のみ Some）
157    pub warm_start_basis: Option<WarmStartBasis>,
158    // --- QP固有フィールド ---
159    /// Bound dual values (shadow prices for variable bounds).
160    ///
161    /// Maps to original variable indices via col_map.
162    /// Empty if no bound constraints are active.
163    ///
164    /// - 除去変数 (presolveで固定された変数) の bound_dual = 0.0 (近似)
165    /// - presolve tightening で追加された境界の dual は報告しない（元問題基準）
166    /// - 配列順: `[lb_dual(j0), ..., lb_dual(j_{n_lb-1}), ub_dual(j0), ..., ub_dual(j_{n_ub-1})]`
167    pub bound_duals: Vec<f64>,
168    /// 反復回数（WSR実績回数）
169    pub iterations: usize,
170    /// 最終反復の残差実値 (pfeas, dfeas, duality_gap)。Optimal/MaxIterations時のみ Some。
171    pub final_residuals: Option<(f64, f64, f64)>,
172    /// 相対双対ギャップ (|p_obj - d_obj| / max(|p|,|d|,1))。
173    /// IPPMM 内部の best-so-far に紐づく値。unscale_ipm_result の Suboptimal→Optimal 昇格ゲート用。
174    /// None = 未計測（LP simplex 等 gap を持たない経路）。
175    pub duality_gap_rel: Option<f64>,
176    /// 各 phase の所要時間 (LP simplex 経路のみ、None なら未計測)。
177    /// 「どこに時間が掛かっているか」事実観測用 (CLAUDE.md「順調に収束に向けて探索」)。
178    pub timing_breakdown: Option<TimingBreakdown>,
179    /// Postsolve が最終的に採用した y_orig の dfeas violation (bound-aware sup ノルム).
180    /// LP simplex + presolve 経路のみ Some。caller (solve_with) が値が `PIVOT_TOL` を
181    /// 超えるとき presolve=off で再解する fallback gate に使う (greenbea-class 問題対策)。
182    pub postsolve_dfeas: Option<f64>,
183    /// Per-solve routing and warm-start statistics (race-free).
184    pub stats: SolveStats,
185    /// Branch-and-bound gap certificate.
186    ///
187    /// Present iff the solver completed a B&B search with a fully authenticated gap
188    /// (no `proof_uncertain` region, `within_gap` satisfied). `None` for direct LP/QP solves.
189    pub bound_gap_cert: Option<BoundGapCertificate>,
190    /// KKT optimality certificate — minted by `prove_optimal` on the B&B incumbent.
191    ///
192    /// Set when `finalize_proven` verifies all KKT conditions on the returned point.
193    /// `None` for demoted (LocallyOptimal/NonconvexLocal) or non-B&B results.
194    pub opt_cert: Option<OptimalCertificate>,
195}
196
197/// 各 phase 所要時間 (μs精度)。LP simplex と QP IPM の両経路で共用。
198#[derive(Debug, Clone, Copy, Default, PartialEq)]
199pub struct TimingBreakdown {
200    // ── LP simplex 経路 ───────────────────────────────────────────────────────
201    /// Presolve 全体 (run_presolve)
202    pub presolve_us: u64,
203    /// 縮約後 simplex 本体 (solve_without_presolve)
204    pub solve_us: u64,
205    /// Postsolve 総計 (run_postsolve / QP 経路では下記 4 field の合計)
206    pub postsolve_us: u64,
207
208    // ── QP IPM 経路: IPM 反復内部 ────────────────────────────────────────────
209    /// KKT 行列 LDL 数値因子化の累計 (全 iteration 合計)
210    pub ipm_factorize_us: u64,
211    /// KKT solve (predictor/corrector/Gondzio) の累計
212    pub ipm_solve_us: u64,
213    /// LDL regularization retry の累計回数 (健全性プローブ失敗含む)
214    pub ipm_reg_retries: u32,
215    /// MINRES (iterative) backend が 1 回以上使われたか
216    pub ipm_used_iterative: bool,
217
218    // ── QP IPM 経路: postsolve 段階別 ────────────────────────────────────────
219    /// postsolve_qp_with_dual_recovery (reduced → orig 空間写像)
220    pub postsolve_map_us: u64,
221    /// refine_postsolve_dual_lsq (元空間 y LSQ refine)
222    pub postsolve_lsq_us: u64,
223    /// refine_postsolve_recovery (Stage 0: SingletonRow 後退代入)
224    pub postsolve_recovery_us: u64,
225    /// refine_post_processing (Stage 1+2: primal projection / y-z refit)
226    pub postsolve_refine_us: u64,
227    /// refine_krylov_and_projection (saddle-point Krylov IR)
228    pub postsolve_krylov_ir_us: u64,
229}
230
231impl Default for SolverResult {
232    fn default() -> Self {
233        SolverResult {
234            status: SolveStatus::NumericalError,
235            objective: 0.0,
236            solution: vec![],
237            dual_solution: vec![],
238            reduced_costs: vec![],
239            slack: vec![],
240            warm_start_basis: None,
241            bound_duals: vec![],
242            iterations: 0,
243            final_residuals: None,
244            duality_gap_rel: None,
245            timing_breakdown: None,
246            postsolve_dfeas: None,
247            stats: SolveStats::default(),
248            bound_gap_cert: None,
249            opt_cert: None,
250        }
251    }
252}
253
254impl fmt::Display for SolverResult {
255    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
256        write!(f, "Status: {}, Objective: {}", self.status, self.objective)
257    }
258}
259
260/// 線形計画問題: min c^T x + obj_offset  s.t.  Ax {op} b,  x in [lb, ub]
261///
262/// 目的関数・制約行列・右辺ベクトル・変数上下限をまとめて保持する。
263/// 制約種別（`<=`, `>=`, `=`）と変数ごとの上下限を個別に指定できる。
264///
265/// `obj_offset` は MPS/QPS の N 行 RHS (objective constant) に対応する。
266/// `LpProblem::new` / `new_general` では 0.0 に初期化される。
267#[derive(Debug, Clone)]
268pub struct LpProblem {
269    /// 目的関数係数ベクトル（長さ: `num_vars`）
270    pub c: Vec<f64>,
271    /// 制約行列（CSC形式、サイズ: `num_constraints` x `num_vars`）
272    pub a: CscMatrix,
273    /// 制約右辺ベクトル（長さ: `num_constraints`）
274    pub b: Vec<f64>,
275    /// 決定変数の数
276    pub num_vars: usize,
277    /// 制約式の数
278    pub num_constraints: usize,
279    /// 各制約の種別（長さ: `num_constraints`）
280    pub constraint_types: Vec<ConstraintType>,
281    /// 各変数の上下限 `(lower, upper)`（長さ: `num_vars`）
282    pub bounds: Vec<(f64, f64)>,
283    /// 問題名（オプション）
284    pub name: Option<String>,
285    /// 目的定数項 (MPS N-row RHS); 最適値 = c^T x* + obj_offset
286    pub obj_offset: f64,
287}
288
289impl LpProblem {
290    /// 新しいLP問題を検証付きで生成する（後方互換版）
291    ///
292    /// 標準形 `min c^T x  s.t.  Ax <= b,  x >= 0` を作成する。
293    /// 全制約を `<=`、全変数の下限を 0・上限を `+∞` とする。
294    ///
295    /// # 引数
296    /// * `c` - 目的関数係数ベクトル
297    /// * `a` - 制約行列（CSC形式）
298    /// * `b` - 制約右辺ベクトル
299    ///
300    /// # 戻り値
301    /// * `Ok(LpProblem)` - 次元が有効な場合
302    /// * `Err(String)` - 次元不一致などの検証エラー時
303    pub fn new(c: Vec<f64>, a: CscMatrix, b: Vec<f64>) -> Result<Self, SolverError> {
304        let num_vars = c.len();
305        let num_constraints = b.len();
306
307        // Set defaults for backward compatibility
308        let constraint_types = vec![ConstraintType::Le; num_constraints];
309        let bounds = vec![(0.0, f64::INFINITY); num_vars];
310        let name = None;
311
312        Self::new_general(c, a, b, constraint_types, bounds, name)
313    }
314
315    /// 制約種別と変数上下限を完全指定して新しいLP問題を生成する
316    ///
317    /// # 引数
318    /// * `c` - 目的関数係数ベクトル
319    /// * `a` - 制約行列（CSC形式）
320    /// * `b` - 制約右辺ベクトル
321    /// * `constraint_types` - 各制約の種別（`Le` / `Ge` / `Eq`）
322    /// * `bounds` - 各変数の上下限 `(lower, upper)`
323    /// * `name` - 問題名（オプション）
324    ///
325    /// # 戻り値
326    /// * `Ok(LpProblem)` - 次元が有効な場合
327    /// * `Err(String)` - 次元不一致などの検証エラー時
328    pub fn new_general(
329        c: Vec<f64>,
330        a: CscMatrix,
331        b: Vec<f64>,
332        constraint_types: Vec<ConstraintType>,
333        bounds: Vec<(f64, f64)>,
334        name: Option<String>,
335    ) -> Result<Self, SolverError> {
336        // Validate dimensions
337        if c.len() != a.ncols {
338            return Err(SolverError::DimensionMismatch {
339                field: "c",
340                expected: a.ncols,
341                got: c.len(),
342            });
343        }
344        if b.len() != a.nrows {
345            return Err(SolverError::DimensionMismatch {
346                field: "b",
347                expected: a.nrows,
348                got: b.len(),
349            });
350        }
351        if constraint_types.len() != b.len() {
352            return Err(SolverError::DimensionMismatch {
353                field: "constraint_types",
354                expected: b.len(),
355                got: constraint_types.len(),
356            });
357        }
358        if bounds.len() != c.len() {
359            return Err(SolverError::DimensionMismatch {
360                field: "bounds",
361                expected: c.len(),
362                got: bounds.len(),
363            });
364        }
365        for (i, &v) in c.iter().enumerate() {
366            if !v.is_finite() {
367                return Err(SolverError::NonFiniteCoefficient {
368                    field: "c",
369                    index: i,
370                });
371            }
372        }
373        for (i, &v) in b.iter().enumerate() {
374            if !v.is_finite() {
375                return Err(SolverError::NonFiniteCoefficient {
376                    field: "b",
377                    index: i,
378                });
379            }
380        }
381        for (i, &v) in a.values.iter().enumerate() {
382            if !v.is_finite() {
383                return Err(SolverError::NonFiniteCoefficient {
384                    field: "A",
385                    index: i,
386                });
387            }
388        }
389        for (i, &(lb, ub)) in bounds.iter().enumerate() {
390            if !is_valid_bound_pair(lb, ub) {
391                return Err(SolverError::InvalidBounds { index: i, lb, ub });
392            }
393        }
394
395        Ok(LpProblem {
396            num_vars: c.len(),
397            num_constraints: b.len(),
398            c,
399            a,
400            b,
401            constraint_types,
402            bounds,
403            name,
404            obj_offset: 0.0,
405        })
406    }
407}
408
409impl fmt::Display for LpProblem {
410    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
411        write!(
412            f,
413            "LP: min c^T x, {} vars, {} constraints",
414            self.num_vars, self.num_constraints
415        )
416    }
417}
418
419#[cfg(test)]
420mod tests {
421    use super::*;
422    use crate::error::SolverError;
423
424    #[test]
425    fn test_lp_problem_new_valid() {
426        // 2 variables, 2 constraints
427        let c = vec![1.0, 2.0];
428        let a = CscMatrix::new(2, 2);
429        let b = vec![5.0, 6.0];
430
431        let lp = LpProblem::new(c, a, b).unwrap();
432        assert_eq!(lp.num_vars, 2);
433        assert_eq!(lp.num_constraints, 2);
434    }
435
436    #[test]
437    fn test_lp_problem_new_invalid_c_dimension() {
438        // c.len() = 3, but a.ncols = 2
439        let c = vec![1.0, 2.0, 3.0];
440        let a = CscMatrix::new(2, 2);
441        let b = vec![5.0, 6.0];
442
443        let result = LpProblem::new(c, a, b);
444        assert!(result.is_err());
445        assert!(matches!(
446            result.unwrap_err(),
447            SolverError::DimensionMismatch { field: "c", .. }
448        ));
449    }
450
451    #[test]
452    fn test_lp_problem_new_invalid_b_dimension() {
453        // b.len() = 3, but a.nrows = 2
454        let c = vec![1.0, 2.0];
455        let a = CscMatrix::new(2, 2);
456        let b = vec![5.0, 6.0, 7.0];
457
458        let result = LpProblem::new(c, a, b);
459        assert!(result.is_err());
460        assert!(matches!(
461            result.unwrap_err(),
462            SolverError::DimensionMismatch { field: "b", .. }
463        ));
464    }
465
466    #[test]
467    fn test_lp_problem_display() {
468        let c = vec![1.0, 2.0];
469        let a = CscMatrix::new(2, 2);
470        let b = vec![5.0, 6.0];
471        let lp = LpProblem::new(c, a, b).unwrap();
472
473        let display = format!("{}", lp);
474        assert_eq!(display, "LP: min c^T x, 2 vars, 2 constraints");
475    }
476
477    #[test]
478    fn test_solve_status_display() {
479        assert_eq!(format!("{}", SolveStatus::Optimal), "Optimal");
480        assert_eq!(format!("{}", SolveStatus::Infeasible), "Infeasible");
481        assert_eq!(format!("{}", SolveStatus::Unbounded), "Unbounded");
482    }
483
484    #[test]
485    fn test_solver_result_display() {
486        let result = SolverResult {
487            status: SolveStatus::Optimal,
488            objective: 42.5,
489            solution: vec![1.0, 2.0],
490            dual_solution: vec![],
491            reduced_costs: vec![],
492            slack: vec![],
493            warm_start_basis: None,
494            ..Default::default()
495        };
496        let display = format!("{}", result);
497        assert_eq!(display, "Status: Optimal, Objective: 42.5");
498    }
499
500    #[test]
501    fn solver_result_default_is_not_success() {
502        let result = SolverResult::default();
503        assert_eq!(result.status, SolveStatus::NumericalError);
504        assert!(result.solution.is_empty());
505    }
506
507    fn make_lp(
508        c: Vec<f64>,
509        b: Vec<f64>,
510        a_vals: Vec<f64>,
511        bounds: Vec<(f64, f64)>,
512    ) -> Result<LpProblem, SolverError> {
513        let n = c.len();
514        let m = b.len();
515        let a = if a_vals.is_empty() {
516            CscMatrix::new(m, n)
517        } else {
518            let rows = vec![0usize; n];
519            let cols: Vec<usize> = (0..n).collect();
520            CscMatrix::from_triplets(&rows, &cols, &a_vals, m, n).unwrap()
521        };
522        let ct = vec![ConstraintType::Le; m];
523        LpProblem::new_general(c, a, b, ct, bounds, None)
524    }
525
526    #[test]
527    fn lp_valid_accepted() {
528        let res = make_lp(
529            vec![1.0, 2.0],
530            vec![5.0],
531            vec![1.0, 1.0],
532            vec![(0.0, f64::INFINITY), (0.0, 10.0)],
533        );
534        assert!(res.is_ok());
535    }
536
537    #[test]
538    fn lp_nan_in_c_rejected() {
539        let bad_vals = [f64::NAN, f64::INFINITY, f64::NEG_INFINITY];
540        for bad in bad_vals {
541            let res = make_lp(
542                vec![bad, 1.0],
543                vec![5.0],
544                vec![1.0, 1.0],
545                vec![(0.0, f64::INFINITY); 2],
546            );
547            assert!(
548                matches!(
549                    res,
550                    Err(SolverError::NonFiniteCoefficient { field: "c", .. })
551                ),
552                "expected NonFiniteCoefficient for c={bad}"
553            );
554        }
555    }
556
557    #[test]
558    fn lp_nan_in_b_rejected() {
559        let bad_vals = [f64::NAN, f64::INFINITY, f64::NEG_INFINITY];
560        for bad in bad_vals {
561            let res = make_lp(
562                vec![1.0, 2.0],
563                vec![bad],
564                vec![1.0, 1.0],
565                vec![(0.0, f64::INFINITY); 2],
566            );
567            assert!(
568                matches!(
569                    res,
570                    Err(SolverError::NonFiniteCoefficient { field: "b", .. })
571                ),
572                "expected NonFiniteCoefficient for b={bad}"
573            );
574        }
575    }
576
577    #[test]
578    fn lp_nan_in_a_rejected() {
579        let n = 2;
580        let bad_vals = [f64::NAN, f64::INFINITY, f64::NEG_INFINITY];
581        for bad in bad_vals {
582            // from_triplets drops NaN via DROP_TOL; inject bad value directly.
583            let mut a = CscMatrix::from_triplets(&[0], &[0], &[1.0], 1, n).unwrap();
584            a.values[0] = bad;
585            let res = LpProblem::new_general(
586                vec![1.0, 2.0],
587                a,
588                vec![5.0],
589                vec![ConstraintType::Le],
590                vec![(0.0, f64::INFINITY); n],
591                None,
592            );
593            assert!(
594                matches!(
595                    res,
596                    Err(SolverError::NonFiniteCoefficient { field: "A", .. })
597                ),
598                "expected NonFiniteCoefficient for A val={bad}"
599            );
600        }
601    }
602
603    #[test]
604    fn lp_nan_in_bounds_rejected() {
605        let cases: Vec<(f64, f64)> = vec![(f64::NAN, 1.0), (0.0, f64::NAN), (f64::NAN, f64::NAN)];
606        for (lb, ub) in cases {
607            let res = make_lp(
608                vec![1.0, 2.0],
609                vec![5.0],
610                vec![1.0, 1.0],
611                vec![(lb, ub), (0.0, f64::INFINITY)],
612            );
613            assert!(
614                matches!(res, Err(SolverError::InvalidBounds { index: 0, .. })),
615                "expected InvalidBounds for ({lb},{ub})"
616            );
617        }
618    }
619
620    #[test]
621    fn lp_lb_gt_ub_rejected() {
622        let cases: Vec<(f64, f64)> = vec![
623            (5.0, 1.0),
624            (1.0, 0.0),
625            (f64::INFINITY, f64::NEG_INFINITY),
626            (0.1, 0.0),
627        ];
628        for (lb, ub) in cases {
629            let res = make_lp(
630                vec![1.0, 2.0],
631                vec![5.0],
632                vec![1.0, 1.0],
633                vec![(lb, ub), (0.0, f64::INFINITY)],
634            );
635            assert!(
636                matches!(res, Err(SolverError::InvalidBounds { .. })),
637                "expected InvalidBounds for lb={lb} ub={ub}"
638            );
639        }
640    }
641
642    #[test]
643    fn lp_inf_bounds_accepted() {
644        let res = make_lp(
645            vec![1.0, 2.0],
646            vec![5.0],
647            vec![1.0, 1.0],
648            vec![(f64::NEG_INFINITY, f64::INFINITY), (0.0, f64::INFINITY)],
649        );
650        assert!(res.is_ok(), "±inf bounds should be valid");
651    }
652
653    // Bug C/D: (inf,inf) and (-inf,-inf) form empty intervals that lb>ub misses
654    // because inf==inf and -inf==-inf. These must be rejected.
655    #[test]
656    fn lp_empty_interval_same_inf_rejected() {
657        let cases: Vec<(f64, f64)> = vec![
658            (f64::INFINITY, f64::INFINITY),
659            (f64::NEG_INFINITY, f64::NEG_INFINITY),
660        ];
661        for (lb, ub) in cases {
662            let res = make_lp(
663                vec![1.0, 2.0],
664                vec![5.0],
665                vec![1.0, 1.0],
666                vec![(lb, ub), (0.0, f64::INFINITY)],
667            );
668            assert!(
669                matches!(res, Err(SolverError::InvalidBounds { index: 0, .. })),
670                "expected InvalidBounds for ({lb},{ub})"
671            );
672        }
673    }
674
675    // Sentinel: valid half-infinite and fixed-var bounds must be accepted.
676    // If bound validation is replaced by a no-op, fixed-var or half-inf cases
677    // would still be Ok (no rejection) — these are *positive* sentinels that
678    // confirm the validator doesn't over-reject.  The reject sentinels above
679    // ensure no-op → FAIL.
680    #[test]
681    fn lp_valid_bound_variants_accepted() {
682        let valid_cases: Vec<Vec<(f64, f64)>> = vec![
683            vec![(f64::NEG_INFINITY, 5.0), (0.0, f64::INFINITY)],
684            vec![(0.0, f64::INFINITY), (f64::NEG_INFINITY, f64::INFINITY)],
685            vec![(3.0, 3.0), (0.0, 10.0)],
686            vec![(0.0, 0.0), (0.0, 0.0)],
687        ];
688        for bounds in valid_cases {
689            let res = make_lp(
690                vec![1.0, 2.0],
691                vec![5.0],
692                vec![1.0, 1.0],
693                bounds.clone(),
694            );
695            assert!(res.is_ok(), "expected Ok for bounds={bounds:?}");
696        }
697    }
698
699    // Sentinel: lb=+inf (alone) and ub=-inf (alone) must each be rejected.
700    #[test]
701    fn lp_lone_inf_bound_rejected() {
702        let cases: Vec<(f64, f64)> = vec![
703            (f64::INFINITY, 5.0),
704            (f64::INFINITY, f64::INFINITY),
705            (0.0, f64::NEG_INFINITY),
706            (f64::NEG_INFINITY, f64::NEG_INFINITY),
707        ];
708        for (lb, ub) in cases {
709            let res = make_lp(
710                vec![1.0, 2.0],
711                vec![5.0],
712                vec![1.0, 1.0],
713                vec![(lb, ub), (0.0, f64::INFINITY)],
714            );
715            assert!(
716                matches!(res, Err(SolverError::InvalidBounds { index: 0, .. })),
717                "expected InvalidBounds for ({lb},{ub})"
718            );
719        }
720    }
721}
otspot_core/problem/mod.rs

otspot_core/problem/
mod.rs
