pounce_cli/
builtin.rs

1//! Built-in TNLP test problems for the CLI. Each problem is a
2//! self-contained `impl TNLP` so the CLI can run end-to-end without
3//! parsing an `.nl` file.
4//!
5//! Currently shipped:
6//!
7//! * `quadratic` — `min (x[0]-3)^2 + (x[1]-4)^2`, unconstrained,
8//!   exact Hessian = `2I`. Optimum at `(3, 4)`, `f* = 0`.
9//! * `rosenbrock` — `min 100*(x[1]-x[0]^2)^2 + (1-x[0])^2`,
10//!   unconstrained, exact Hessian. Optimum at `(1, 1)`, `f* = 0`.
11
12use pounce_common::types::{Index, Number};
13use pounce_nlp::tnlp::{
14    BoundsInfo, IndexStyle, IpoptCq, IpoptData, NlpInfo, Solution, SparsityRequest, StartingPoint,
15    TNLP,
16};
17use std::cell::RefCell;
18use std::rc::Rc;
19
20pub fn list() -> Vec<&'static str> {
21    vec![
22        "quadratic",
23        "rosenbrock",
24        "bounded-quadratic",
25        "eq-quadratic",
26        "circle",
27        "infeasible-eq",
28        "wachter-biegler",
29    ]
30}
31
32pub fn lookup(name: &str) -> Option<Rc<RefCell<dyn TNLP>>> {
33    match name {
34        "quadratic" => Some(Rc::new(RefCell::new(Quadratic::default()))),
35        "rosenbrock" => Some(Rc::new(RefCell::new(Rosenbrock::default()))),
36        "bounded-quadratic" => Some(Rc::new(RefCell::new(BoundedQuadratic::default()))),
37        "eq-quadratic" => Some(Rc::new(RefCell::new(EqQuadratic::default()))),
38        "circle" => Some(Rc::new(RefCell::new(Circle::default()))),
39        "infeasible-eq" => Some(Rc::new(RefCell::new(InfeasibleEq::default()))),
40        "wachter-biegler" => Some(Rc::new(RefCell::new(WachterBiegler::default()))),
41        _ => None,
42    }
43}
44
45// --------------------------------------------------------------------
46// Quadratic: min (x0 - 3)^2 + (x1 - 4)^2
47// --------------------------------------------------------------------
48
49#[derive(Debug, Default)]
50pub struct Quadratic {
51    pub final_x: Option<[Number; 2]>,
52    pub final_obj: Number,
53}
54
55impl TNLP for Quadratic {
56    fn get_nlp_info(&mut self) -> Option<NlpInfo> {
57        Some(NlpInfo {
58            n: 2,
59            m: 0,
60            nnz_jac_g: 0,
61            nnz_h_lag: 2, // diagonal Hessian, lower triangle
62            index_style: IndexStyle::C,
63        })
64    }
65
66    fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
67        b.x_l.iter_mut().for_each(|v| *v = -2e19);
68        b.x_u.iter_mut().for_each(|v| *v = 2e19);
69        true
70    }
71
72    fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
73        sp.x.copy_from_slice(&[0.0, 0.0]);
74        true
75    }
76
77    fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
78        Some((x[0] - 3.0).powi(2) + (x[1] - 4.0).powi(2))
79    }
80
81    fn eval_grad_f(&mut self, x: &[Number], _new_x: bool, grad: &mut [Number]) -> bool {
82        grad[0] = 2.0 * (x[0] - 3.0);
83        grad[1] = 2.0 * (x[1] - 4.0);
84        true
85    }
86
87    fn eval_g(&mut self, _x: &[Number], _new_x: bool, _g: &mut [Number]) -> bool {
88        true
89    }
90
91    fn eval_jac_g(
92        &mut self,
93        _x: Option<&[Number]>,
94        _new_x: bool,
95        _mode: SparsityRequest<'_>,
96    ) -> bool {
97        true
98    }
99
100    fn eval_h(
101        &mut self,
102        _x: Option<&[Number]>,
103        _new_x: bool,
104        obj_factor: Number,
105        _lambda: Option<&[Number]>,
106        _new_lambda: bool,
107        mode: SparsityRequest<'_>,
108    ) -> bool {
109        match mode {
110            SparsityRequest::Structure { irow, jcol } => {
111                irow.copy_from_slice(&[0, 1]);
112                jcol.copy_from_slice(&[0, 1]);
113            }
114            SparsityRequest::Values { values } => {
115                values[0] = 2.0 * obj_factor;
116                values[1] = 2.0 * obj_factor;
117            }
118        }
119        true
120    }
121
122    fn finalize_solution(&mut self, sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {
123        self.final_x = Some([sol.x[0], sol.x[1]]);
124        self.final_obj = sol.obj_value;
125    }
126}
127
128// --------------------------------------------------------------------
129// Rosenbrock: min 100 (x1 - x0^2)^2 + (1 - x0)^2
130// --------------------------------------------------------------------
131
132#[derive(Debug, Default)]
133pub struct Rosenbrock {
134    pub final_x: Option<[Number; 2]>,
135    pub final_obj: Number,
136}
137
138impl TNLP for Rosenbrock {
139    fn get_nlp_info(&mut self) -> Option<NlpInfo> {
140        Some(NlpInfo {
141            n: 2,
142            m: 0,
143            nnz_jac_g: 0,
144            nnz_h_lag: 3, // dense 2x2 lower triangle: (0,0), (1,0), (1,1)
145            index_style: IndexStyle::C,
146        })
147    }
148
149    fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
150        b.x_l.iter_mut().for_each(|v| *v = -2e19);
151        b.x_u.iter_mut().for_each(|v| *v = 2e19);
152        true
153    }
154
155    fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
156        sp.x.copy_from_slice(&[-1.2, 1.0]);
157        true
158    }
159
160    fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
161        let a = x[1] - x[0] * x[0];
162        let b = 1.0 - x[0];
163        Some(100.0 * a * a + b * b)
164    }
165
166    fn eval_grad_f(&mut self, x: &[Number], _new_x: bool, grad: &mut [Number]) -> bool {
167        // d/dx0 = -400 x0 (x1 - x0^2) - 2 (1 - x0)
168        // d/dx1 =  200 (x1 - x0^2)
169        grad[0] = -400.0 * x[0] * (x[1] - x[0] * x[0]) - 2.0 * (1.0 - x[0]);
170        grad[1] = 200.0 * (x[1] - x[0] * x[0]);
171        true
172    }
173
174    fn eval_g(&mut self, _x: &[Number], _new_x: bool, _g: &mut [Number]) -> bool {
175        true
176    }
177
178    fn eval_jac_g(
179        &mut self,
180        _x: Option<&[Number]>,
181        _new_x: bool,
182        _mode: SparsityRequest<'_>,
183    ) -> bool {
184        true
185    }
186
187    fn eval_h(
188        &mut self,
189        x: Option<&[Number]>,
190        _new_x: bool,
191        obj_factor: Number,
192        _lambda: Option<&[Number]>,
193        _new_lambda: bool,
194        mode: SparsityRequest<'_>,
195    ) -> bool {
196        match mode {
197            SparsityRequest::Structure { irow, jcol } => {
198                // Lower triangle (row >= col): (0,0), (1,0), (1,1).
199                irow.copy_from_slice(&[0, 1, 1]);
200                jcol.copy_from_slice(&[0, 0, 1]);
201            }
202            SparsityRequest::Values { values } => {
203                let x = x.unwrap_or(&[0.0, 0.0]);
204                let h00 = -400.0 * (x[1] - 3.0 * x[0] * x[0]) + 2.0;
205                let h10 = -400.0 * x[0];
206                let h11 = 200.0;
207                values[0] = obj_factor * h00;
208                values[1] = obj_factor * h10;
209                values[2] = obj_factor * h11;
210            }
211        }
212        true
213    }
214
215    fn finalize_solution(&mut self, sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {
216        self.final_x = Some([sol.x[0], sol.x[1]]);
217        self.final_obj = sol.obj_value;
218    }
219}
220
221// Small helper for doc references
222#[allow(dead_code)]
223fn _ix(_: Index) {}
224
225// --------------------------------------------------------------------
226// BoundedQuadratic: min (x0-3)^2 + (x1-4)^2 s.t. 0 <= x0 <= 2, 0 <= x1 <= 2
227// Optimum is at the corner (2, 2) with f* = 1 + 4 = 5.
228// --------------------------------------------------------------------
229
230#[derive(Debug, Default)]
231pub struct BoundedQuadratic {
232    pub final_x: Option<[Number; 2]>,
233    pub final_obj: Number,
234}
235
236impl TNLP for BoundedQuadratic {
237    fn get_nlp_info(&mut self) -> Option<NlpInfo> {
238        Some(NlpInfo {
239            n: 2,
240            m: 0,
241            nnz_jac_g: 0,
242            nnz_h_lag: 2,
243            index_style: IndexStyle::C,
244        })
245    }
246
247    fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
248        b.x_l.copy_from_slice(&[0.0, 0.0]);
249        b.x_u.copy_from_slice(&[2.0, 2.0]);
250        true
251    }
252
253    fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
254        sp.x.copy_from_slice(&[1.0, 1.0]);
255        true
256    }
257
258    fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
259        Some((x[0] - 3.0).powi(2) + (x[1] - 4.0).powi(2))
260    }
261
262    fn eval_grad_f(&mut self, x: &[Number], _new_x: bool, grad: &mut [Number]) -> bool {
263        grad[0] = 2.0 * (x[0] - 3.0);
264        grad[1] = 2.0 * (x[1] - 4.0);
265        true
266    }
267
268    fn eval_g(&mut self, _x: &[Number], _new_x: bool, _g: &mut [Number]) -> bool {
269        true
270    }
271
272    fn eval_jac_g(
273        &mut self,
274        _x: Option<&[Number]>,
275        _new_x: bool,
276        _mode: SparsityRequest<'_>,
277    ) -> bool {
278        true
279    }
280
281    fn eval_h(
282        &mut self,
283        _x: Option<&[Number]>,
284        _new_x: bool,
285        obj_factor: Number,
286        _lambda: Option<&[Number]>,
287        _new_lambda: bool,
288        mode: SparsityRequest<'_>,
289    ) -> bool {
290        match mode {
291            SparsityRequest::Structure { irow, jcol } => {
292                irow.copy_from_slice(&[0, 1]);
293                jcol.copy_from_slice(&[0, 1]);
294            }
295            SparsityRequest::Values { values } => {
296                values[0] = 2.0 * obj_factor;
297                values[1] = 2.0 * obj_factor;
298            }
299        }
300        true
301    }
302
303    fn finalize_solution(&mut self, sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {
304        self.final_x = Some([sol.x[0], sol.x[1]]);
305        self.final_obj = sol.obj_value;
306    }
307}
308
309// --------------------------------------------------------------------
310// EqQuadratic: min x0^2 + x1^2  s.t.  x0 + x1 = 1
311// Optimum at (1/2, 1/2), f* = 1/2, multiplier y = -1.
312// --------------------------------------------------------------------
313
314#[derive(Debug, Default)]
315pub struct EqQuadratic {
316    pub final_x: Option<[Number; 2]>,
317    pub final_obj: Number,
318}
319
320impl TNLP for EqQuadratic {
321    fn get_nlp_info(&mut self) -> Option<NlpInfo> {
322        Some(NlpInfo {
323            n: 2,
324            m: 1,
325            nnz_jac_g: 2,
326            nnz_h_lag: 2,
327            index_style: IndexStyle::C,
328        })
329    }
330
331    fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
332        b.x_l.iter_mut().for_each(|v| *v = -2e19);
333        b.x_u.iter_mut().for_each(|v| *v = 2e19);
334        b.g_l[0] = 1.0;
335        b.g_u[0] = 1.0;
336        true
337    }
338
339    fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
340        sp.x.copy_from_slice(&[0.0, 0.0]);
341        true
342    }
343
344    fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
345        Some(x[0] * x[0] + x[1] * x[1])
346    }
347
348    fn eval_grad_f(&mut self, x: &[Number], _new_x: bool, grad: &mut [Number]) -> bool {
349        grad[0] = 2.0 * x[0];
350        grad[1] = 2.0 * x[1];
351        true
352    }
353
354    fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
355        g[0] = x[0] + x[1];
356        true
357    }
358
359    fn eval_jac_g(
360        &mut self,
361        _x: Option<&[Number]>,
362        _new_x: bool,
363        mode: SparsityRequest<'_>,
364    ) -> bool {
365        match mode {
366            SparsityRequest::Structure { irow, jcol } => {
367                irow.copy_from_slice(&[0, 0]);
368                jcol.copy_from_slice(&[0, 1]);
369            }
370            SparsityRequest::Values { values } => {
371                values[0] = 1.0;
372                values[1] = 1.0;
373            }
374        }
375        true
376    }
377
378    fn eval_h(
379        &mut self,
380        _x: Option<&[Number]>,
381        _new_x: bool,
382        obj_factor: Number,
383        _lambda: Option<&[Number]>,
384        _new_lambda: bool,
385        mode: SparsityRequest<'_>,
386    ) -> bool {
387        match mode {
388            SparsityRequest::Structure { irow, jcol } => {
389                irow.copy_from_slice(&[0, 1]);
390                jcol.copy_from_slice(&[0, 1]);
391            }
392            SparsityRequest::Values { values } => {
393                values[0] = 2.0 * obj_factor;
394                values[1] = 2.0 * obj_factor;
395            }
396        }
397        true
398    }
399
400    fn finalize_solution(&mut self, sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {
401        self.final_x = Some([sol.x[0], sol.x[1]]);
402        self.final_obj = sol.obj_value;
403    }
404}
405
406// --------------------------------------------------------------------
407// Circle: min  x0  s.t.  x0^2 + x1^2 = 1
408// Optimum at (-1, 0), f* = -1, multiplier y = 1/2.
409// Tests nonlinear equality constraint with non-trivial Hessian
410// contribution from the constraint (∇²g_0 = 2I) into the Lagrangian.
411// --------------------------------------------------------------------
412
413#[derive(Debug, Default)]
414pub struct Circle {
415    pub final_x: Option<[Number; 2]>,
416    pub final_obj: Number,
417}
418
419impl TNLP for Circle {
420    fn get_nlp_info(&mut self) -> Option<NlpInfo> {
421        Some(NlpInfo {
422            n: 2,
423            m: 1,
424            nnz_jac_g: 2,
425            nnz_h_lag: 2,
426            index_style: IndexStyle::C,
427        })
428    }
429
430    fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
431        b.x_l.iter_mut().for_each(|v| *v = -2e19);
432        b.x_u.iter_mut().for_each(|v| *v = 2e19);
433        b.g_l[0] = 1.0;
434        b.g_u[0] = 1.0;
435        true
436    }
437
438    fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
439        sp.x.copy_from_slice(&[-0.5, 0.5]);
440        true
441    }
442
443    fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
444        Some(x[0])
445    }
446
447    fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, grad: &mut [Number]) -> bool {
448        grad[0] = 1.0;
449        grad[1] = 0.0;
450        true
451    }
452
453    fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
454        g[0] = x[0] * x[0] + x[1] * x[1];
455        true
456    }
457
458    fn eval_jac_g(
459        &mut self,
460        x: Option<&[Number]>,
461        _new_x: bool,
462        mode: SparsityRequest<'_>,
463    ) -> bool {
464        match mode {
465            SparsityRequest::Structure { irow, jcol } => {
466                irow.copy_from_slice(&[0, 0]);
467                jcol.copy_from_slice(&[0, 1]);
468            }
469            SparsityRequest::Values { values } => {
470                let x = x.unwrap_or(&[0.0, 0.0]);
471                values[0] = 2.0 * x[0];
472                values[1] = 2.0 * x[1];
473            }
474        }
475        true
476    }
477
478    fn eval_h(
479        &mut self,
480        _x: Option<&[Number]>,
481        _new_x: bool,
482        _obj_factor: Number,
483        lambda: Option<&[Number]>,
484        _new_lambda: bool,
485        mode: SparsityRequest<'_>,
486    ) -> bool {
487        match mode {
488            SparsityRequest::Structure { irow, jcol } => {
489                irow.copy_from_slice(&[0, 1]);
490                jcol.copy_from_slice(&[0, 1]);
491            }
492            SparsityRequest::Values { values } => {
493                // ∇²L = obj_factor * ∇²f + λ * ∇²g_0 = 0 + λ * 2I
494                let lam = lambda.map(|l| l[0]).unwrap_or(0.0);
495                values[0] = 2.0 * lam;
496                values[1] = 2.0 * lam;
497            }
498        }
499        true
500    }
501
502    fn finalize_solution(&mut self, sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {
503        self.final_x = Some([sol.x[0], sol.x[1]]);
504        self.final_obj = sol.obj_value;
505    }
506}
507
508// --------------------------------------------------------------------
509// InfeasibleEq: min x0^2 + x1^2
510//   s.t.  x0 + x1 = 1   (g_0)
511//         x0 + x1 = 2   (g_1)
512// The two equalities are mutually contradictory, so no feasible point
513// exists. The standard solve drives the restoration phase, which also
514// cannot achieve feasibility, returning Restoration_Failed. With
515// `l1_fallback_on_restoration_failure=yes` (or
516// `l1_exact_penalty_barrier=yes`), the CLI then performs a second
517// inner solve via the ℓ₁-exact penalty-barrier wrapper. That second
518// pass is what exercises the multi-pass restoration factory provider
519// path — the very path that previously panicked with
520// "restoration factory invoked more than once".
521// --------------------------------------------------------------------
522
523#[derive(Debug, Default)]
524pub struct InfeasibleEq {
525    pub final_x: Option<[Number; 2]>,
526    pub final_obj: Number,
527}
528
529impl TNLP for InfeasibleEq {
530    fn get_nlp_info(&mut self) -> Option<NlpInfo> {
531        Some(NlpInfo {
532            n: 2,
533            m: 2,
534            nnz_jac_g: 4,
535            nnz_h_lag: 2,
536            index_style: IndexStyle::C,
537        })
538    }
539
540    fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
541        b.x_l.iter_mut().for_each(|v| *v = -2e19);
542        b.x_u.iter_mut().for_each(|v| *v = 2e19);
543        b.g_l[0] = 1.0;
544        b.g_u[0] = 1.0;
545        b.g_l[1] = 2.0;
546        b.g_u[1] = 2.0;
547        true
548    }
549
550    fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
551        sp.x.copy_from_slice(&[0.0, 0.0]);
552        true
553    }
554
555    fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
556        Some(x[0] * x[0] + x[1] * x[1])
557    }
558
559    fn eval_grad_f(&mut self, x: &[Number], _new_x: bool, grad: &mut [Number]) -> bool {
560        grad[0] = 2.0 * x[0];
561        grad[1] = 2.0 * x[1];
562        true
563    }
564
565    fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
566        g[0] = x[0] + x[1];
567        g[1] = x[0] + x[1];
568        true
569    }
570
571    fn eval_jac_g(
572        &mut self,
573        _x: Option<&[Number]>,
574        _new_x: bool,
575        mode: SparsityRequest<'_>,
576    ) -> bool {
577        match mode {
578            SparsityRequest::Structure { irow, jcol } => {
579                irow.copy_from_slice(&[0, 0, 1, 1]);
580                jcol.copy_from_slice(&[0, 1, 0, 1]);
581            }
582            SparsityRequest::Values { values } => {
583                values.copy_from_slice(&[1.0, 1.0, 1.0, 1.0]);
584            }
585        }
586        true
587    }
588
589    fn eval_h(
590        &mut self,
591        _x: Option<&[Number]>,
592        _new_x: bool,
593        obj_factor: Number,
594        _lambda: Option<&[Number]>,
595        _new_lambda: bool,
596        mode: SparsityRequest<'_>,
597    ) -> bool {
598        match mode {
599            SparsityRequest::Structure { irow, jcol } => {
600                irow.copy_from_slice(&[0, 1]);
601                jcol.copy_from_slice(&[0, 1]);
602            }
603            SparsityRequest::Values { values } => {
604                values[0] = 2.0 * obj_factor;
605                values[1] = 2.0 * obj_factor;
606            }
607        }
608        true
609    }
610
611    fn finalize_solution(&mut self, sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {
612        self.final_x = Some([sol.x[0], sol.x[1]]);
613        self.final_obj = sol.obj_value;
614    }
615}
616
617// --------------------------------------------------------------------
618// WachterBiegler: the Wächter–Biegler counterexample (scaled).
619//
620//   min  x0
621//   s.t. x0^2 - 10*x1 - 1 = 0     (g_0)
622//        x0    -  3*x2 - 0.5 = 0  (g_1)
623//        x1 >= 0,  x2 >= 0
624//   start (-2, 3, 1)
625//
626// Feasibility forces x0 >= 1 (x1 = (x0^2-1)/10 >= 0 needs |x0| >= 1;
627// x2 = (x0-0.5)/3 >= 0 needs x0 >= 0.5), so the optimum is
628// (1, 0, 1/6) with f* = 1.
629//
630// This is the canonical hard case for line-search interior-point /
631// filter methods (Wächter & Biegler 2000, "Failure of global
632// convergence for a class of interior point methods for nonlinear
633// programming"). From this start the default `monotone` barrier takes a
634// first step that drives x0 strongly negative; the dual infeasibility
635// blows up and the solve enters restoration and converges to a point of
636// local infeasibility near x0 ~ -1, never reaching the true optimum. A
637// cold `adaptive` (mu_strategy=adaptive) start solves it. It is shipped
638// as a built-in so the debugger screencast can demonstrate diagnosing
639// and steering past a default-strategy failure.
640// --------------------------------------------------------------------
641
642#[derive(Debug, Default)]
643pub struct WachterBiegler {
644    pub final_x: Option<[Number; 3]>,
645    pub final_obj: Number,
646}
647
648impl TNLP for WachterBiegler {
649    fn get_nlp_info(&mut self) -> Option<NlpInfo> {
650        Some(NlpInfo {
651            n: 3,
652            m: 2,
653            // g_0: d/dx0 = 2*x0, d/dx1 = -10  -> 2 nz
654            // g_1: d/dx0 = 1,    d/dx2 = -3   -> 2 nz
655            nnz_jac_g: 4,
656            // Only g_0 has a nonzero second derivative: d2/dx0^2 = 2.
657            nnz_h_lag: 1,
658            index_style: IndexStyle::C,
659        })
660    }
661
662    fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
663        b.x_l.copy_from_slice(&[-2e19, 0.0, 0.0]);
664        b.x_u.copy_from_slice(&[2e19, 2e19, 2e19]);
665        // Both constraints are equalities (= 0).
666        b.g_l.copy_from_slice(&[0.0, 0.0]);
667        b.g_u.copy_from_slice(&[0.0, 0.0]);
668        true
669    }
670
671    fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
672        sp.x.copy_from_slice(&[-2.0, 3.0, 1.0]);
673        true
674    }
675
676    fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
677        Some(x[0])
678    }
679
680    fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, grad: &mut [Number]) -> bool {
681        grad[0] = 1.0;
682        grad[1] = 0.0;
683        grad[2] = 0.0;
684        true
685    }
686
687    fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
688        g[0] = x[0] * x[0] - 10.0 * x[1] - 1.0;
689        g[1] = x[0] - 3.0 * x[2] - 0.5;
690        true
691    }
692
693    fn eval_jac_g(
694        &mut self,
695        x: Option<&[Number]>,
696        _new_x: bool,
697        mode: SparsityRequest<'_>,
698    ) -> bool {
699        match mode {
700            SparsityRequest::Structure { irow, jcol } => {
701                // (0,0), (0,1), (1,0), (1,2)
702                irow.copy_from_slice(&[0, 0, 1, 1]);
703                jcol.copy_from_slice(&[0, 1, 0, 2]);
704            }
705            SparsityRequest::Values { values } => {
706                let x = x.unwrap_or(&[0.0, 0.0, 0.0]);
707                values[0] = 2.0 * x[0]; // d g_0 / d x0
708                values[1] = -10.0; // d g_0 / d x1
709                values[2] = 1.0; // d g_1 / d x0
710                values[3] = -3.0; // d g_1 / d x2
711            }
712        }
713        true
714    }
715
716    fn eval_h(
717        &mut self,
718        _x: Option<&[Number]>,
719        _new_x: bool,
720        _obj_factor: Number,
721        lambda: Option<&[Number]>,
722        _new_lambda: bool,
723        mode: SparsityRequest<'_>,
724    ) -> bool {
725        match mode {
726            SparsityRequest::Structure { irow, jcol } => {
727                // Only (0,0) is ever nonzero (f and g_1 are linear; g_0
728                // is quadratic only in x0).
729                irow.copy_from_slice(&[0]);
730                jcol.copy_from_slice(&[0]);
731            }
732            SparsityRequest::Values { values } => {
733                // ∇²L = obj_factor*0 + λ_0 * ∇²g_0 + λ_1 * ∇²g_1
734                //     = λ_0 * 2 (at (0,0)), λ_1 term is 0.
735                let lam0 = lambda.map(|l| l[0]).unwrap_or(0.0);
736                values[0] = 2.0 * lam0;
737            }
738        }
739        true
740    }
741
742    fn finalize_solution(&mut self, sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {
743        self.final_x = Some([sol.x[0], sol.x[1], sol.x[2]]);
744        self.final_obj = sol.obj_value;
745    }
746}
747
748#[cfg(test)]
749mod tests {
750    use super::*;
751
752    #[test]
753    fn list_contains_known_problems() {
754        let l = list();
755        assert!(l.contains(&"quadratic"));
756        assert!(l.contains(&"rosenbrock"));
757    }
758
759    #[test]
760    fn quadratic_evaluates_correctly() {
761        let mut q = Quadratic::default();
762        let f = q.eval_f(&[3.0, 4.0], true).unwrap();
763        assert_eq!(f, 0.0);
764        let mut g = [0.0; 2];
765        q.eval_grad_f(&[0.0, 0.0], true, &mut g);
766        assert_eq!(g, [-6.0, -8.0]);
767    }
768
769    #[test]
770    fn rosenbrock_grad_zero_at_optimum() {
771        let mut r = Rosenbrock::default();
772        let f = r.eval_f(&[1.0, 1.0], true).unwrap();
773        assert!(f.abs() < 1e-15);
774        let mut g = [0.0; 2];
775        r.eval_grad_f(&[1.0, 1.0], true, &mut g);
776        assert!(g[0].abs() < 1e-12);
777        assert!(g[1].abs() < 1e-12);
778    }
779
780    #[test]
781    fn lookup_returns_known_and_rejects_unknown() {
782        assert!(lookup("quadratic").is_some());
783        assert!(lookup("rosenbrock").is_some());
784        assert!(lookup("wachter-biegler").is_some());
785        assert!(lookup("nonsense").is_none());
786    }
787
788    #[test]
789    fn wachter_biegler_optimum_is_feasible() {
790        // True optimum (1, 0, 1/6): both equalities hold and f = x0 = 1.
791        let mut w = WachterBiegler::default();
792        let xopt = [1.0, 0.0, 1.0 / 6.0];
793        let mut g = [0.0; 2];
794        w.eval_g(&xopt, true, &mut g);
795        assert!(g[0].abs() < 1e-12, "g0 = {}", g[0]);
796        assert!(g[1].abs() < 1e-12, "g1 = {}", g[1]);
797        assert_eq!(w.eval_f(&xopt, true).unwrap(), 1.0);
798    }
799
800    #[test]
801    fn wachter_biegler_jacobian_matches_finite_difference() {
802        // Spot-check the analytic Jacobian against a central difference
803        // at the start point (-2, 3, 1).
804        let mut w = WachterBiegler::default();
805        let x = [-2.0, 3.0, 1.0];
806        let mut vals = [0.0; 4];
807        w.eval_jac_g(
808            Some(&x),
809            true,
810            SparsityRequest::Values { values: &mut vals },
811        );
812        // (0,0)=2*x0=-4, (0,1)=-10, (1,0)=1, (1,2)=-3.
813        assert!((vals[0] - (-4.0)).abs() < 1e-12);
814        assert!((vals[1] - (-10.0)).abs() < 1e-12);
815        assert!((vals[2] - 1.0).abs() < 1e-12);
816        assert!((vals[3] - (-3.0)).abs() < 1e-12);
817    }
818}
pounce_cli/builtin.rs

pounce_cli/
builtin.rs
