sensitivity_factor_reuse_bench/
sensitivity_factor_reuse_bench.rs

1//! Tracking-only timing comparison: factor-once + N sensitivity steps
2//! vs N fresh cold-start solves. The point of `pounce_sensitivity::Solver`
3//! is that the converged KKT factor survives between operations, so
4//! each follow-up parametric step / reduced-Hessian / KKT-solve is a
5//! back-solve rather than a full IPM re-run.
6//!
7//! Not a regression test — numbers vary across machines and builds.
8//! Run with:
9//!   `cargo run --release --example sensitivity_factor_reuse_bench -p pounce-sensitivity`
10
11use std::cell::RefCell;
12use std::rc::Rc;
13use std::time::Instant;
14
15use pounce_algorithm::application::IpoptApplication;
16use pounce_common::types::{Index, Number};
17use pounce_nlp::tnlp::{
18    BoundsInfo, IndexStyle, IpoptCq, IpoptData, NlpInfo, Solution, SparsityRequest, StartingPoint,
19    TNLP,
20};
21use pounce_sensitivity::{SensSolve, Solver};
22
23struct ParametricTNLP {
24    eta1: Number,
25    eta2: Number,
26}
27
28impl TNLP for ParametricTNLP {
29    fn get_nlp_info(&mut self) -> Option<NlpInfo> {
30        Some(NlpInfo {
31            n: 5,
32            m: 4,
33            nnz_jac_g: 10,
34            nnz_h_lag: 5,
35            index_style: IndexStyle::C,
36        })
37    }
38    fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
39        for k in 0..3 {
40            b.x_l[k] = 0.0;
41            b.x_u[k] = 1.0e19;
42        }
43        b.x_l[3] = -1.0e19;
44        b.x_u[3] = 1.0e19;
45        b.x_l[4] = -1.0e19;
46        b.x_u[4] = 1.0e19;
47        b.g_l[0] = 0.0;
48        b.g_u[0] = 0.0;
49        b.g_l[1] = 0.0;
50        b.g_u[1] = 0.0;
51        b.g_l[2] = self.eta1;
52        b.g_u[2] = self.eta1;
53        b.g_l[3] = self.eta2;
54        b.g_u[3] = self.eta2;
55        true
56    }
57    fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
58        sp.x[0] = 0.15;
59        sp.x[1] = 0.15;
60        sp.x[2] = 0.0;
61        sp.x[3] = 0.0;
62        sp.x[4] = 0.0;
63        true
64    }
65    fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
66        Some(x[0] * x[0] + x[1] * x[1] + x[2] * x[2])
67    }
68    fn eval_grad_f(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
69        g[0] = 2.0 * x[0];
70        g[1] = 2.0 * x[1];
71        g[2] = 2.0 * x[2];
72        g[3] = 0.0;
73        g[4] = 0.0;
74        true
75    }
76    fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
77        let (x1, x2, x3, e1, e2) = (x[0], x[1], x[2], x[3], x[4]);
78        g[0] = 6.0 * x1 + 3.0 * x2 + 2.0 * x3 - e1;
79        g[1] = e2 * x1 + x2 - x3 - 1.0;
80        g[2] = e1;
81        g[3] = e2;
82        true
83    }
84    fn eval_jac_g(
85        &mut self,
86        x: Option<&[Number]>,
87        _new_x: bool,
88        mode: SparsityRequest<'_>,
89    ) -> bool {
90        match mode {
91            SparsityRequest::Structure { irow, jcol } => {
92                let rs: [Index; 10] = [0, 0, 0, 0, 1, 1, 1, 1, 2, 3];
93                let cs: [Index; 10] = [0, 1, 2, 3, 0, 1, 2, 4, 3, 4];
94                irow.copy_from_slice(&rs);
95                jcol.copy_from_slice(&cs);
96            }
97            SparsityRequest::Values { values } => {
98                let x = x.expect("Values without x");
99                values[0] = 6.0;
100                values[1] = 3.0;
101                values[2] = 2.0;
102                values[3] = -1.0;
103                values[4] = x[4];
104                values[5] = 1.0;
105                values[6] = -1.0;
106                values[7] = x[0];
107                values[8] = 1.0;
108                values[9] = 1.0;
109            }
110        }
111        true
112    }
113    fn eval_h(
114        &mut self,
115        _x: Option<&[Number]>,
116        _new_x: bool,
117        obj_factor: Number,
118        lambda: Option<&[Number]>,
119        _new_lambda: bool,
120        mode: SparsityRequest<'_>,
121    ) -> bool {
122        match mode {
123            SparsityRequest::Structure { irow, jcol } => {
124                let rs: [Index; 5] = [0, 1, 2, 4, 0];
125                let cs: [Index; 5] = [0, 1, 2, 0, 0];
126                irow.copy_from_slice(&rs);
127                jcol.copy_from_slice(&cs);
128            }
129            SparsityRequest::Values { values } => {
130                let lam = lambda.expect("Values without lambda");
131                values[0] = 2.0 * obj_factor;
132                values[1] = 2.0 * obj_factor;
133                values[2] = 2.0 * obj_factor;
134                values[3] = lam[1];
135                values[4] = 0.0;
136            }
137        }
138        true
139    }
140    fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
141}
142
143fn make_app() -> IpoptApplication {
144    let mut app = IpoptApplication::new();
145    app.options_mut()
146        .set_integer_value("print_level", 0, true, false)
147        .unwrap();
148    app.options_mut()
149        .set_string_value("sb", "yes", true, false)
150        .unwrap();
151    app.initialize().unwrap();
152    app
153}
154
155const N: usize = 20;
156
157fn main() {
158    let deltas: Vec<Vec<Number>> = (1..=N).map(|k| vec![0.0, 0.01 * k as Number]).collect();
159    let pins = vec![2 as Index, 3];
160
161    // (1) Held-factor path: 1 solve + N parametric steps.
162    let t = Instant::now();
163    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
164        eta1: 5.0,
165        eta2: 1.0,
166    }));
167    let mut solver = Solver::new(make_app(), tnlp);
168    solver.solve();
169    let solve_dt = t.elapsed();
170    let mut held_steps = Vec::with_capacity(N);
171    let t = Instant::now();
172    for d in &deltas {
173        let dx = solver
174            .parametric_step(&pins, d)
175            .expect("parametric_step ok");
176        held_steps.push(dx);
177    }
178    let held_step_dt = t.elapsed();
179    let held_total = solve_dt + held_step_dt;
180
181    // (2) Cold path: N fresh SensSolve runs (one full IPM each).
182    let mut cold_steps = Vec::with_capacity(N);
183    let t = Instant::now();
184    for d in &deltas {
185        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
186            eta1: 5.0,
187            eta2: 1.0,
188        }));
189        let mut app = make_app();
190        let r = SensSolve::new(pins.clone())
191            .with_deltas(d.clone())
192            .run(&mut app, tnlp);
193        cold_steps.push(r.dx.expect("dx"));
194    }
195    let cold_total = t.elapsed();
196
197    // Cross-check: held vs cold dx should agree.
198    let mut max_err = 0.0_f64;
199    for (h, c) in held_steps.iter().zip(cold_steps.iter()) {
200        for (a, b) in h.iter().zip(c.iter()) {
201            max_err = max_err.max((a - b).abs());
202        }
203    }
204
205    println!("Workload: 1 IPM solve + {N} parametric steps (Δeta2 sweep).");
206    println!();
207    println!(
208        "Held-factor path : solve {:.2} ms + {N} steps {:.2} ms = {:.2} ms total",
209        solve_dt.as_secs_f64() * 1e3,
210        held_step_dt.as_secs_f64() * 1e3,
211        held_total.as_secs_f64() * 1e3,
212    );
213    println!(
214        "Cold-restart path:                          {N} fresh solves = {:.2} ms total",
215        cold_total.as_secs_f64() * 1e3,
216    );
217    println!();
218    println!(
219        "Speedup of held-factor over cold-restart: {:.1}x",
220        cold_total.as_secs_f64() / held_total.as_secs_f64()
221    );
222    println!("Numerical agreement (held vs cold dx): max |err| = {max_err:.2e}");
223}
sensitivity_factor_reuse_bench/sensitivity_factor_reuse_bench.rs

sensitivity_factor_reuse_bench/
sensitivity_factor_reuse_bench.rs
