sensitivity_session/
sensitivity_session.rs

1//! Sensitivity workflow without callbacks.
2//!
3//! Solves the parametric NLP from
4//! `tests/parametric_cpp.rs::ParametricTNLP` once, then issues several
5//! cheap follow-up operations against the held KKT factor:
6//!
7//!   * `parametric_step` for two different parameter perturbations
8//!   * `compute_reduced_hessian` over the pinned-row set
9//!   * `kkt_solve` for a raw back-solve against a synthetic RHS
10//!
11//! Each follow-up reuses the symbolic + numeric factor cached inside
12//! the linear-solver backend — no IPM iteration, no refactorization.
13//!
14//! Run with:
15//!   `cargo run --example sensitivity_session -p pounce-sensitivity`
16
17use std::cell::RefCell;
18use std::rc::Rc;
19
20use pounce_algorithm::application::IpoptApplication;
21use pounce_common::types::{Index, Number};
22use pounce_nlp::tnlp::{
23    BoundsInfo, IndexStyle, IpoptCq, IpoptData, NlpInfo, Solution, SparsityRequest, StartingPoint,
24    TNLP,
25};
26use pounce_sensitivity::Solver;
27
28struct ParametricTNLP {
29    eta1: Number,
30    eta2: Number,
31}
32
33impl TNLP for ParametricTNLP {
34    fn get_nlp_info(&mut self) -> Option<NlpInfo> {
35        Some(NlpInfo {
36            n: 5,
37            m: 4,
38            nnz_jac_g: 10,
39            nnz_h_lag: 5,
40            index_style: IndexStyle::C,
41        })
42    }
43    fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
44        for k in 0..3 {
45            b.x_l[k] = 0.0;
46            b.x_u[k] = 1.0e19;
47        }
48        b.x_l[3] = -1.0e19;
49        b.x_u[3] = 1.0e19;
50        b.x_l[4] = -1.0e19;
51        b.x_u[4] = 1.0e19;
52        b.g_l[0] = 0.0;
53        b.g_u[0] = 0.0;
54        b.g_l[1] = 0.0;
55        b.g_u[1] = 0.0;
56        b.g_l[2] = self.eta1;
57        b.g_u[2] = self.eta1;
58        b.g_l[3] = self.eta2;
59        b.g_u[3] = self.eta2;
60        true
61    }
62    fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
63        sp.x[0] = 0.15;
64        sp.x[1] = 0.15;
65        sp.x[2] = 0.0;
66        sp.x[3] = 0.0;
67        sp.x[4] = 0.0;
68        true
69    }
70    fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
71        Some(x[0] * x[0] + x[1] * x[1] + x[2] * x[2])
72    }
73    fn eval_grad_f(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
74        g[0] = 2.0 * x[0];
75        g[1] = 2.0 * x[1];
76        g[2] = 2.0 * x[2];
77        g[3] = 0.0;
78        g[4] = 0.0;
79        true
80    }
81    fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
82        let (x1, x2, x3, e1, e2) = (x[0], x[1], x[2], x[3], x[4]);
83        g[0] = 6.0 * x1 + 3.0 * x2 + 2.0 * x3 - e1;
84        g[1] = e2 * x1 + x2 - x3 - 1.0;
85        g[2] = e1;
86        g[3] = e2;
87        true
88    }
89    fn eval_jac_g(
90        &mut self,
91        x: Option<&[Number]>,
92        _new_x: bool,
93        mode: SparsityRequest<'_>,
94    ) -> bool {
95        match mode {
96            SparsityRequest::Structure { irow, jcol } => {
97                let rs: [Index; 10] = [0, 0, 0, 0, 1, 1, 1, 1, 2, 3];
98                let cs: [Index; 10] = [0, 1, 2, 3, 0, 1, 2, 4, 3, 4];
99                irow.copy_from_slice(&rs);
100                jcol.copy_from_slice(&cs);
101            }
102            SparsityRequest::Values { values } => {
103                let x = x.expect("Values without x");
104                values[0] = 6.0;
105                values[1] = 3.0;
106                values[2] = 2.0;
107                values[3] = -1.0;
108                values[4] = x[4];
109                values[5] = 1.0;
110                values[6] = -1.0;
111                values[7] = x[0];
112                values[8] = 1.0;
113                values[9] = 1.0;
114            }
115        }
116        true
117    }
118    fn eval_h(
119        &mut self,
120        _x: Option<&[Number]>,
121        _new_x: bool,
122        obj_factor: Number,
123        lambda: Option<&[Number]>,
124        _new_lambda: bool,
125        mode: SparsityRequest<'_>,
126    ) -> bool {
127        match mode {
128            SparsityRequest::Structure { irow, jcol } => {
129                let rs: [Index; 5] = [0, 1, 2, 4, 0];
130                let cs: [Index; 5] = [0, 1, 2, 0, 0];
131                irow.copy_from_slice(&rs);
132                jcol.copy_from_slice(&cs);
133            }
134            SparsityRequest::Values { values } => {
135                let lam = lambda.expect("Values without lambda");
136                values[0] = 2.0 * obj_factor;
137                values[1] = 2.0 * obj_factor;
138                values[2] = 2.0 * obj_factor;
139                values[3] = lam[1];
140                values[4] = 0.0;
141            }
142        }
143        true
144    }
145    fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
146}
147
148fn make_app() -> IpoptApplication {
149    let mut app = IpoptApplication::new();
150    app.options_mut()
151        .set_integer_value("print_level", 0, true, false)
152        .unwrap();
153    app.options_mut()
154        .set_string_value("sb", "yes", true, false)
155        .unwrap();
156    app.initialize().unwrap();
157    app
158}
159
160fn main() {
161    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
162        eta1: 5.0,
163        eta2: 1.0,
164    }));
165
166    let mut solver = Solver::new(make_app(), tnlp);
167    let status = solver.solve();
168    println!("solve status: {status:?}");
169    assert!(solver.converged().is_some(), "solver did not converge");
170
171    let pins = vec![2 as Index, 3];
172
173    // Two cheap parametric steps against the same factor.
174    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
175        let dx = solver
176            .parametric_step(&pins, deltas)
177            .expect("parametric_step ok");
178        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
179    }
180
181    // Reduced Hessian over the same pinned-row set.
182    let hr = solver
183        .compute_reduced_hessian(&pins, 1.0)
184        .expect("reduced Hessian ok");
185    println!("reduced Hessian (2x2, column-major) = {hr:?}");
186
187    // Raw back-solve against a zero RHS — must come back zero.
188    let dim = solver.kkt_dim().expect("kkt_dim available");
189    let rhs = vec![0.0; dim];
190    let mut lhs = vec![1.0; dim];
191    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
192    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
193    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
194    assert!(max_abs < 1e-10);
195}
sensitivity_session/sensitivity_session.rs

sensitivity_session/
sensitivity_session.rs
