pounce_sensitivity/
barrier_iterate.rs

1//! pounce's own converged iterate behind [`BarrierIterate`], the trait the
2//! corrector and the activity classification in `pounce-sens-core` read.
3//!
4//! Everything here reads `IpoptData`, the calculated quantities and the NLP
5//! of the solve the backsolver holds, in the scaled frame the algorithm
6//! iterated in. The arithmetic of both analyses lives in `pounce-sens-core`.
7
8use std::rc::Rc;
9
10use pounce_algorithm::ipopt_data::IpoptDataHandle;
11use pounce_algorithm::iterates_vector::IteratesVector;
12use pounce_common::types::{Index, Number};
13use pounce_linalg::Matrix;
14use pounce_linalg::expansion_matrix::ExpansionMatrix;
15use pounce_sens_core::barrier::{
16    ActivityData, BarrierIterate, BoundSide, EngineScaling, KktOperator, ModelMaps,
17};
18
19use crate::algorithm_backsolver::{CorrectorSigma, PdSensBacksolver};
20use crate::backsolver::SensBacksolver;
21use crate::vec_util::dense_to_vec;
22
23impl BarrierIterate for PdSensBacksolver {
24    fn block_dims(&self) -> [usize; 8] {
25        PdSensBacksolver::block_dims(self)
26    }
27
28    fn offsets(&self) -> [usize; 9] {
29        self.offsets_public()
30    }
31
32    fn barrier_mu(&self) -> Number {
33        PdSensBacksolver::barrier_mu(self)
34    }
35
36    fn converged_iterate(&self, out: &mut [Number]) -> Result<(), String> {
37        self.curr_flat(out)
38            .map_err(|_| "the solve holds no current iterate of that length".to_string())
39    }
40
41    fn scaled_rhs_factor(&self) -> Option<&[Number]> {
42        PdSensBacksolver::scaled_rhs_factor(self)
43    }
44
45    /// Installs `point` as the trial iterate and reads the residual off the
46    /// calculated quantities there, so scaling, fixed variables and the bound
47    /// expansions are handled exactly as the solve handles them. Setting a
48    /// trial point leaves `curr` alone.
49    fn barrier_residual(&self, point: &[Number], out: &mut [Number]) -> Result<(), String> {
50        let dim = self.dim();
51        if point.len() != dim || out.len() != dim {
52            return Err(format!(
53                "barrier residual: point of length {} and out of length {} for a KKT system of {dim}",
54                point.len(),
55                out.len()
56            ));
57        }
58        let iv = self
59            .pack_public(point)
60            .map_err(|_| "corrector: pack failed".to_string())?;
61        let (data, cq, _) = self.activity_handles();
62        data.borrow_mut().set_trial(iv.freeze());
63
64        let off = self.offsets_public();
65        let cqb = cq.borrow();
66        let blocks: [Rc<dyn pounce_linalg::Vector>; 8] = [
67            cqb.trial_grad_lag_x(),
68            cqb.trial_grad_lag_s(),
69            cqb.trial_c(),
70            cqb.trial_d_minus_s(),
71            cqb.trial_compl_x_l(),
72            cqb.trial_compl_x_u(),
73            cqb.trial_compl_s_l(),
74            cqb.trial_compl_s_u(),
75        ];
76        for (i, v) in blocks.into_iter().enumerate() {
77            let vals = dense_to_vec(&*v);
78            let (a, b) = (off[i], off[i + 1]);
79            if vals.len() != b - a {
80                return Err(format!(
81                    "corrector: block {i} is {} long, expected {}",
82                    vals.len(),
83                    b - a
84                ));
85            }
86            out[a..b].copy_from_slice(&vals);
87        }
88        Ok(())
89    }
90
91    /// Swaps the predicted point in as the current iterate for the life of
92    /// the returned operator, so every block of the operator evaluates there,
93    /// and builds the barrier diagonal after the swap. The stored base-point
94    /// diagonal, which the ceiling (gh#737) or crossover (gh#654) would
95    /// otherwise freeze, is never consulted: the frame rule and the ceiling
96    /// are re-derived at the predicted point. See `corrector_sigma`.
97    ///
98    /// The solve takes its Hessian from `IpoptData::w`, the matrix the
99    /// algorithm last stored, so moving the iterate alone moves only the
100    /// Jacobians and slacks. Under an exact-Hessian solve the matrix is
101    /// re-evaluated at the predicted point, with the step's own multipliers.
102    /// A limited-memory solve keeps its quasi-Newton matrix, since there is
103    /// no exact Hessian to evaluate anywhere.
104    fn corrector_operator<'a>(
105        &'a self,
106        predicted: &[Number],
107        pinned: &[(usize, Number)],
108        exact_hessian: bool,
109    ) -> Result<Box<dyn KktOperator + 'a>, String> {
110        let (data_h, cq_h, _) = self.activity_handles();
111        let predicted = self
112            .pack_natural(predicted)
113            .ok_or_else(|| "corrector: packing the predicted point failed".to_string())?;
114        let mut restore = RestoreCurr {
115            data: data_h,
116            saved: data_h.borrow().curr.clone(),
117            saved_w: None,
118            restore_w: false,
119        };
120        data_h.borrow_mut().set_curr(predicted);
121        if exact_hessian {
122            restore.saved_w = data_h.borrow().w.clone();
123            restore.restore_w = true;
124            let w = cq_h.borrow().curr_exact_hessian();
125            data_h.borrow_mut().w = Some(w);
126        }
127        let op = self
128            .corrector_sigma(pinned)
129            .ok_or_else(|| "corrector: operator diagonal unavailable".to_string())?;
130        Ok(Box::new(PdCorrectorOperator {
131            bs: self,
132            op,
133            _restore: restore,
134        }))
135    }
136
137    fn activity_data(&self) -> Result<ActivityData, String> {
138        let (data, cq, nlp) = self.activity_handles();
139
140        // Scoped borrows: the calculated-quantity getters below re-borrow
141        // the NLP (mutably, for lazy evaluation) and the data, so nothing
142        // here may hold either across one of those calls.
143        let (z_l, z_u, v_l, v_u) = {
144            let d = data.borrow();
145            let curr = d
146                .curr
147                .as_ref()
148                .ok_or("activity: the converged state has no iterate")?;
149            (
150                dense_to_vec(curr.z_l.as_ref()),
151                dense_to_vec(curr.z_u.as_ref()),
152                dense_to_vec(curr.v_l.as_ref()),
153                dense_to_vec(curr.v_u.as_ref()),
154            )
155        };
156        let (px_l, px_u, pd_l, pd_u) = {
157            let nl = nlp.borrow();
158            (
159                positions(&nl.px_l())?,
160                positions(&nl.px_u())?,
161                positions(&nl.pd_l())?,
162                positions(&nl.pd_u())?,
163            )
164        };
165        let (s_x_l, s_x_u, s_d_l, s_d_u, hessian, jac_c, jac_d) = {
166            let c = cq.borrow();
167            (
168                dense_to_vec(c.curr_slack_x_l().as_ref()),
169                dense_to_vec(c.curr_slack_x_u().as_ref()),
170                dense_to_vec(c.curr_slack_s_l().as_ref()),
171                dense_to_vec(c.curr_slack_s_u().as_ref()),
172                c.curr_exact_hessian(),
173                c.curr_jac_c(),
174                c.curr_jac_d(),
175            )
176        };
177        let sigma_x = dense_to_vec(self.barrier_sigma_x().as_ref());
178        let sigma_s = dense_to_vec(self.barrier_sigma_s().as_ref());
179
180        let side =
181            |pos: Vec<usize>, mult: Vec<Number>, slack: Vec<Number>| BoundSide { pos, mult, slack };
182        Ok(ActivityData {
183            x_l: side(px_l, z_l, s_x_l),
184            x_u: side(px_u, z_u, s_x_u),
185            d_l: side(pd_l, v_l, s_d_l),
186            d_u: side(pd_u, v_u, s_d_u),
187            sigma_x,
188            sigma_s,
189            hessian,
190            jac_c,
191            jac_d,
192        })
193    }
194
195    /// The NLP's own maps. The inequality positions are its
196    /// `full_g_to_d_block`, the exact complement of `full_g_to_c_block`:
197    /// one ascending scan over the model's rows assigns each to the next
198    /// position of its block.
199    fn model_maps(&self) -> Result<ModelMaps, String> {
200        let (_, _, nlp) = self.activity_handles();
201        let nl = nlp.borrow();
202        let n_x = self.block_dims()[0];
203        let n_full_g = nl.n_full_g();
204        let block = |p: Option<Index>| p.map(|p| p as usize);
205        Ok(ModelMaps {
206            n_full_x: nl.n_full_x() as usize,
207            var_x_to_full_x: (0..n_x)
208                .map(|i| nl.var_x_to_full_x(i as Index) as usize)
209                .collect(),
210            full_g_to_c: (0..n_full_g)
211                .map(|g| block(nl.full_g_to_c_block(g)))
212                .collect(),
213            full_g_to_d: (0..n_full_g)
214                .map(|g| block(nl.full_g_to_d_block(g)))
215                .collect(),
216        })
217    }
218
219    fn scaling(&self) -> EngineScaling {
220        let (obj, c_scale, d_scale) = self.nlp_scaling();
221        EngineScaling {
222            objective: obj,
223            variables: self.variable_scaling().map(|v| v.to_vec()),
224            equality_rows: c_scale,
225            inequality_rows: d_scale,
226        }
227    }
228
229    /// The bounds the NLP holds, relaxed by `bound_relax_factor` when the
230    /// solve relaxed them, expanded from its compressed bound vectors.
231    fn primal_bounds(&self) -> Result<(Vec<Number>, Vec<Number>), String> {
232        let (_, _, nlp) = self.activity_handles();
233        let nl = nlp.borrow();
234        let [n_x, n_s, ..] = self.block_dims();
235        let (mut lo, mut hi) = pounce_sens_core::boundcheck::expand_bounds(
236            n_x,
237            &nl.px_l(),
238            &nl.px_u(),
239            nl.x_l(),
240            nl.x_u(),
241        );
242        let (s_lo, s_hi) = pounce_sens_core::boundcheck::expand_bounds(
243            n_s,
244            &nl.pd_l(),
245            &nl.pd_u(),
246            nl.d_l(),
247            nl.d_u(),
248        );
249        lo.extend_from_slice(&s_lo);
250        hi.extend_from_slice(&s_hi);
251        Ok((lo, hi))
252    }
253
254    fn kkt_perturbations(&self) -> Option<[Number; 4]> {
255        Some(PdSensBacksolver::kkt_perturbations(self))
256    }
257
258    /// The released diagonals built once and held, so every solve passes
259    /// the same objects and the factorization cache, which keys on the
260    /// diagonal's tag, factors the released system once.
261    fn released_operator<'a>(
262        &'a self,
263        released: &[usize],
264    ) -> Result<Box<dyn KktOperator + 'a>, String> {
265        let (sigma, sigma_s) = self
266            .released_sigmas(released)
267            .ok_or_else(|| "released sigma unavailable".to_string())?;
268        Ok(Box::new(PdReleasedOperator {
269            bs: self,
270            released: released.to_vec(),
271            sigma,
272            sigma_s,
273        }))
274    }
275}
276
277/// The operator [`PdSensBacksolver`] hands a calculation that solves
278/// against one released system many times.
279struct PdReleasedOperator<'a> {
280    bs: &'a PdSensBacksolver,
281    released: Vec<usize>,
282    sigma: Rc<dyn pounce_linalg::Vector>,
283    sigma_s: Option<Rc<dyn pounce_linalg::Vector>>,
284}
285
286impl KktOperator for PdReleasedOperator<'_> {
287    /// [`SensBacksolver::solve_released`] against the held diagonals: it
288    /// does nothing to the right-hand side.
289    fn solve(&self, rhs: &[Number], lhs: &mut [Number]) -> bool {
290        self.bs.solve_released_prebuilt(
291            &self.released,
292            Rc::clone(&self.sigma),
293            self.sigma_s.as_ref().map(Rc::clone),
294            None,
295            rhs,
296            lhs,
297            false,
298        )
299    }
300}
301
302/// The entries a bound expansion places its compressed vector at.
303fn positions(m: &Rc<dyn Matrix>) -> Result<Vec<usize>, String> {
304    m.as_any()
305        .downcast_ref::<ExpansionMatrix>()
306        .map(|e| {
307            e.expanded_pos_indices()
308                .iter()
309                .map(|&i| i as usize)
310                .collect()
311        })
312        .ok_or_else(|| "activity: a bound projection is not an expansion matrix".to_string())
313}
314
315/// The operator [`PdSensBacksolver`] hands the corrector: a solve against
316/// the held factorization with the predicted point's diagonal, and the
317/// iterate swap that diagonal was built under, undone on drop.
318struct PdCorrectorOperator<'a> {
319    bs: &'a PdSensBacksolver,
320    op: CorrectorSigma,
321    _restore: RestoreCurr<'a>,
322}
323
324impl KktOperator for PdCorrectorOperator<'_> {
325    /// The same `Rc` every call: the factorization cache keys on its tag,
326    /// so the predicted point's operator is factored once and every later
327    /// solve is a back-solve.
328    ///
329    /// The ratios go with the diagonals: the operator eliminates each bound
330    /// row into a diagonal the gh#737 ceiling may have softened, so the row
331    /// is read back through the same cap (gh#828). A no-op wherever the
332    /// ceiling did not bind.
333    ///
334    /// No bound is released here. A bound the corrector released has a
335    /// zero multiplier at the point the diagonal was built at, and a zero
336    /// right-hand side, so it adds nothing to the system.
337    fn solve(&self, rhs: &[Number], lhs: &mut [Number]) -> bool {
338        self.bs.solve_released_prebuilt(
339            &[],
340            Rc::clone(&self.op.sigma_x),
341            Some(Rc::clone(&self.op.sigma_s)),
342            Some((&self.op.ratio_x, &self.op.ratio_s)),
343            rhs,
344            lhs,
345            false,
346        )
347    }
348}
349
350/// Puts back the current iterate, and the stored Hessian when it was
351/// replaced, that [`BarrierIterate::corrector_operator`] swapped out.
352struct RestoreCurr<'a> {
353    data: &'a IpoptDataHandle,
354    saved: Option<IteratesVector>,
355    saved_w: Option<Rc<dyn pounce_linalg::SymMatrix>>,
356    restore_w: bool,
357}
358
359impl Drop for RestoreCurr<'_> {
360    fn drop(&mut self) {
361        if let Some(c) = self.saved.take() {
362            self.data.borrow_mut().set_curr(c);
363        }
364        if self.restore_w {
365            self.data.borrow_mut().w = self.saved_w.take();
366        }
367    }
368}
pounce_sensitivity/barrier_iterate.rs

pounce_sensitivity/
barrier_iterate.rs
