otspot_core/
lp.rs

1//! LP-specific entry point.
2//!
3//! Splits LP from the QP `Q.is_zero` dispatch so that LP-only paths
4//! (simplex, future IPM-first / crash / postsolve) are owned by this
5//! module. `solve_qp_with(Q=0)` keeps backward compat by forwarding
6//! here; the two call sites are distinguishable via `SolverResult.stats.route`.
7
8use crate::options::SolverOptions;
9use crate::problem::{LpProblem, SolveRoute, SolveStatus, SolverResult};
10
11/// Solve an LP directly. Sets `result.stats.route = SolveRoute::LpDirect`.
12///
13/// Returns [`SolveStatus::NumericalError`] if `options` fails validation;
14/// validation is performed by the underlying `simplex::solve_with`.
15pub fn solve_lp_with(problem: &LpProblem, options: &SolverOptions) -> SolverResult {
16    let mut result = crate::simplex::solve_with(problem, options);
17    if matches!(
18        result.status,
19        SolveStatus::Optimal | SolveStatus::SuboptimalSolution | SolveStatus::Timeout
20    ) {
21        result.objective += problem.obj_offset;
22    }
23    result.stats.route = SolveRoute::LpDirect;
24    result.stats.deadline_triggered = matches!(result.status, SolveStatus::Timeout);
25    result
26}
27
28/// LP entry from `solve_qp_with(Q=0)`. Sets `result.stats.route = SolveRoute::LpForwardedFromQp`.
29pub(crate) fn solve_lp_forwarded_from_qp(
30    problem: &LpProblem,
31    options: &SolverOptions,
32) -> SolverResult {
33    let mut result = crate::simplex::solve_with(problem, options);
34    if matches!(
35        result.status,
36        SolveStatus::Optimal | SolveStatus::SuboptimalSolution | SolveStatus::Timeout
37    ) {
38        result.objective += problem.obj_offset;
39    }
40    result.stats.route = SolveRoute::LpForwardedFromQp;
41    result.stats.deadline_triggered = matches!(result.status, SolveStatus::Timeout);
42    result
43}
44
45#[cfg(test)]
46mod tests {
47    use super::*;
48    use crate::problem::ConstraintType;
49    use crate::sparse::CscMatrix;
50
51    fn make_trivial_lp() -> LpProblem {
52        // minimize x  s.t.  x <= 5,  x >= 0
53        let a = CscMatrix::from_triplets(&[0], &[0], &[1.0], 1, 1).unwrap();
54        LpProblem::new_general(
55            vec![1.0],
56            a,
57            vec![5.0],
58            vec![ConstraintType::Le],
59            vec![(0.0, f64::INFINITY)],
60            None,
61        )
62        .unwrap()
63    }
64
65    /// Timeout incumbent must include `problem.obj_offset`.
66    ///
67    /// Sentinel: removing `SolveStatus::Timeout` from the match in `solve_lp_with`
68    /// causes `result.objective == 0.0` instead of 42.5 → FAIL.
69    ///
70    /// `cancel_flag = true` with `deadline = None` bypasses the pre-simplex
71    /// INFINITY timeout (entry.rs only checks `deadline.is_some_and(...)`).
72    /// The simplex loop's first-iteration cancel check fires → Timeout with
73    /// initial BFS (x_decision = 0, c^T x = 0, sf.obj_offset = 0).
74    #[test]
75    fn test_lp_timeout_incumbent_includes_obj_offset() {
76        use std::sync::{atomic::AtomicBool, Arc};
77
78        let mut lp = make_trivial_lp();
79        lp.obj_offset = 42.5;
80
81        let opts = SolverOptions {
82            cancel_flag: Some(Arc::new(AtomicBool::new(true))),
83            presolve: false,
84            ..Default::default()
85        };
86
87        let result = solve_lp_with(&lp, &opts);
88        assert_eq!(
89            result.status,
90            SolveStatus::Timeout,
91            "cancel_flag=true must produce Timeout"
92        );
93        assert!(
94            result.objective.is_finite(),
95            "Timeout incumbent must have finite objective (not INFINITY); got {}",
96            result.objective
97        );
98        assert!(
99            (result.objective - 42.5).abs() < 1e-9,
100            "Timeout incumbent must include obj_offset 42.5; got {} \
101             (sentinel: removing Timeout from match yields 0.0 ≠ 42.5)",
102            result.objective
103        );
104    }
105
106    /// Invalid options produce NumericalError via `solve_lp_with`.
107    ///
108    /// Validation is performed by `simplex::solve_with` (the load-bearing sentinel
109    /// lives in `simplex::entry::invalid_options_rejected_at_simplex_entry`).
110    #[test]
111    fn invalid_options_rejected_at_lp_entry() {
112        let lp = make_trivial_lp();
113        let cases: &[(&str, SolverOptions)] = &[
114            (
115                "nan primal_tol",
116                SolverOptions {
117                    primal_tol: f64::NAN,
118                    ..Default::default()
119                },
120            ),
121            (
122                "inf primal_tol",
123                SolverOptions {
124                    primal_tol: f64::INFINITY,
125                    ..Default::default()
126                },
127            ),
128            (
129                "neg timeout_secs",
130                SolverOptions {
131                    timeout_secs: Some(-0.5),
132                    ..Default::default()
133                },
134            ),
135            (
136                "zero threads",
137                SolverOptions {
138                    threads: 0,
139                    ..Default::default()
140                },
141            ),
142            (
143                "nan dual_tol",
144                SolverOptions {
145                    dual_tol: f64::NAN,
146                    ..Default::default()
147                },
148            ),
149        ];
150        for (label, opts) in cases {
151            let result = solve_lp_with(&lp, opts);
152            assert_eq!(
153                result.status,
154                SolveStatus::NumericalError,
155                "solve_lp_with with {label} must return NumericalError"
156            );
157        }
158    }
159
160    fn cx(lp: &LpProblem, sol: &[f64]) -> f64 {
161        lp.c.iter().zip(sol).map(|(c, x)| c * x).sum::<f64>() + lp.obj_offset
162    }
163
164    fn solve_lp_no_presolve(lp: &LpProblem) -> SolverResult {
165        let opts = SolverOptions {
166            presolve: false,
167            ..Default::default()
168        };
169        solve_lp_with(lp, &opts)
170    }
171
172    /// Reported objective must equal `c·x` of the returned solution for an LP
173    /// with a NONZERO lower bound that routes through the Big-M Phase I path.
174    ///
175    /// `min x  s.t.  x >= 5 (Ge),  x ∈ [3, ∞)` → optimum x=5, obj=5.
176    /// The Ge constraint forces artificials ⇒ `big_m_cold_start`. The standard
177    /// form shifts x = 3 + x', so `sf.obj_offset = c·lb = 3`.
178    ///
179    /// BUG (a4200da): `phase1.rs` recomputes `obj_orig = c·solution` from the
180    /// un-shifted solution (already = c·x = 5) and then ADDS `sf.obj_offset`
181    /// again ⇒ reports 8. The solution (x=5) is correct; only the scalar is wrong.
182    /// This test FAILS until the Big-M path stops double-adding `sf.obj_offset`.
183    /// Expected after fix: reported objective == 5.
184    #[test]
185    fn bigm_nonzero_lb_objective_double_count() {
186        let a = CscMatrix::from_triplets(&[0], &[0], &[1.0], 1, 1).unwrap();
187        let lp = LpProblem::new_general(
188            vec![1.0],
189            a,
190            vec![5.0],
191            vec![ConstraintType::Ge],
192            vec![(3.0, f64::INFINITY)],
193            None,
194        )
195        .unwrap();
196        let res = solve_lp_no_presolve(&lp);
197        assert_eq!(res.status, SolveStatus::Optimal);
198        assert!(
199            (res.solution[0] - 5.0).abs() < 1e-6,
200            "solution must be x=5; got {}",
201            res.solution[0]
202        );
203        assert!(
204            (res.objective - cx(&lp, &res.solution)).abs() < 1e-6,
205            "reported objective {} must equal c·x {} (Big-M path double-counts \
206             sf.obj_offset = c·lb = 3 ⇒ reports 8 instead of 5)",
207            res.objective,
208            cx(&lp, &res.solution)
209        );
210    }
211
212    /// Control: Le-only LP with nonzero lower bound routes through the Le-only
213    /// cold-start path, which correctly adds `sf.obj_offset` to the SHIFTED
214    /// `basic_obj`. `min x s.t. x <= 10, x ∈ [3, ∞)` → x=3, obj=3. PASSES.
215    /// Sentinel for the scope of the Big-M bug (this path must stay correct).
216    #[test]
217    fn le_only_nonzero_lb_objective_correct() {
218        let a = CscMatrix::from_triplets(&[0], &[0], &[1.0], 1, 1).unwrap();
219        let lp = LpProblem::new_general(
220            vec![1.0],
221            a,
222            vec![10.0],
223            vec![ConstraintType::Le],
224            vec![(3.0, f64::INFINITY)],
225            None,
226        )
227        .unwrap();
228        let res = solve_lp_no_presolve(&lp);
229        assert_eq!(res.status, SolveStatus::Optimal);
230        assert!(
231            (res.objective - cx(&lp, &res.solution)).abs() < 1e-6,
232            "Le-only path: reported {} must equal c·x {}",
233            res.objective,
234            cx(&lp, &res.solution)
235        );
236    }
237
238    /// Control: bounded LP (finite ub) with nonzero lower bound routes through
239    /// the BFRT bounded path, which is also correct. `min x s.t. x <= 10,
240    /// x ∈ [3, 8]` → x=3, obj=3. PASSES. Sentinel for the bounded path.
241    #[test]
242    fn bounded_nonzero_lb_objective_correct() {
243        let a = CscMatrix::from_triplets(&[0], &[0], &[1.0], 1, 1).unwrap();
244        let lp = LpProblem::new_general(
245            vec![1.0],
246            a,
247            vec![10.0],
248            vec![ConstraintType::Le],
249            vec![(3.0, 8.0)],
250            None,
251        )
252        .unwrap();
253        let res = solve_lp_no_presolve(&lp);
254        assert_eq!(res.status, SolveStatus::Optimal);
255        assert!(
256            (res.objective - cx(&lp, &res.solution)).abs() < 1e-6,
257            "bounded path: reported {} must equal c·x {}",
258            res.objective,
259            cx(&lp, &res.solution)
260        );
261    }
262}
otspot_core/lp.rs

otspot_core/
lp.rs
