pounce_algorithm/kkt/
perturbation_handler.rs

1//! Primal-dual perturbation handler — port of
2//! `Algorithm/IpPDPerturbationHandler.{hpp,cpp}`.
3//!
4//! Owns the four perturbations `(δ_x, δ_s, δ_c, δ_d)` that the
5//! `PDFullSpaceSolver` adds to the augmented system to recover correct
6//! inertia / non-singularity. Implements upstream's full state
7//! machine:
8//!
9//! * [`Self::consider_new_system`] — first call per new aug-system.
10//!   Finalizes the previous trial's degeneracy probe, decides whether
11//!   to start a new degeneracy test, and seeds `δ_c` / `δ_d` if the
12//!   Jacobian is already known to be degenerate (or `perturb_always_cd`
13//!   is on).
14//! * [`Self::perturb_for_singular`] — escalation step taken when MA57
15//!   reports `Singular`.
16//! * [`Self::perturb_for_wrong_inertia`] — escalation step taken when
17//!   the factor's negative-eigenvalue count disagrees with what the
18//!   KKT structure requires.
19//! * [`Self::current_perturbation`] — read the most recently committed
20//!   `(δ_x, δ_s, δ_c, δ_d)`.
21//!
22//! Returns `false` when no further escalation is possible (caller must
23//! enter the restoration phase). The `info_string`-mutation calls in
24//! upstream are emitted via the `IpoptData` handle the caller passes
25//! in; if `None` is passed, the strings are simply dropped.
26
27use crate::ipopt_data::IpoptDataHandle;
28use pounce_common::types::{Index, Number};
29
30/// Trial state — port of upstream `TrialStatus` enum.
31#[derive(Debug, Clone, Copy, PartialEq, Eq)]
32pub enum TrialStatus {
33    NoTest,
34    DcEq0DxEq0,
35    DcGt0DxEq0,
36    DcEq0DxGt0,
37    DcGt0DxGt0,
38}
39
40/// Degeneracy state — port of `DegenType`.
41#[derive(Debug, Clone, Copy, PartialEq, Eq)]
42pub enum DegenType {
43    NotDegenerate,
44    Degenerate,
45    NotYetDetermined,
46}
47
48/// State + algorithmic parameters. Defaults mirror
49/// `IpPDPerturbationHandler.cpp::RegisterOptions`.
50#[derive(Debug, Clone)]
51pub struct PdPerturbationHandler {
52    // ---- algorithmic parameters (read from options) ----
53    pub delta_xs_max: Number,
54    pub delta_xs_min: Number,
55    pub delta_xs_first_inc_fact: Number,
56    pub delta_xs_inc_fact: Number,
57    pub delta_xs_dec_fact: Number,
58    pub delta_xs_init: Number,
59    pub delta_cd_val: Number,
60    pub delta_cd_exp: Number,
61    pub perturb_always_cd: bool,
62    pub reset_last: bool,
63    pub degen_iters_max: Index,
64
65    // ---- live state ----
66    pub delta_x_curr: Number,
67    pub delta_s_curr: Number,
68    pub delta_c_curr: Number,
69    pub delta_d_curr: Number,
70    pub delta_x_last: Number,
71    pub delta_s_last: Number,
72    pub delta_c_last: Number,
73    pub delta_d_last: Number,
74    pub get_deltas_for_wrong_inertia_called: bool,
75    pub hess_degenerate: DegenType,
76    pub jac_degenerate: DegenType,
77    pub degen_iters: Index,
78    pub test_status: TrialStatus,
79}
80
81impl Default for PdPerturbationHandler {
82    fn default() -> Self {
83        Self {
84            delta_xs_max: 1e20,
85            delta_xs_min: 1e-20,
86            delta_xs_first_inc_fact: 100.0,
87            delta_xs_inc_fact: 8.0,
88            delta_xs_dec_fact: 1.0 / 3.0,
89            delta_xs_init: 1e-4,
90            delta_cd_val: 1e-8,
91            delta_cd_exp: 0.25,
92            perturb_always_cd: false,
93            reset_last: false,
94            degen_iters_max: 3,
95            delta_x_curr: 0.0,
96            delta_s_curr: 0.0,
97            delta_c_curr: 0.0,
98            delta_d_curr: 0.0,
99            delta_x_last: 0.0,
100            delta_s_last: 0.0,
101            delta_c_last: 0.0,
102            delta_d_last: 0.0,
103            get_deltas_for_wrong_inertia_called: false,
104            hess_degenerate: DegenType::NotYetDetermined,
105            jac_degenerate: DegenType::NotYetDetermined,
106            degen_iters: 0,
107            test_status: TrialStatus::NoTest,
108        }
109    }
110}
111
112/// Snapshot of the four perturbations after a state-machine call.
113#[derive(Debug, Clone, Copy, PartialEq)]
114pub struct Deltas {
115    pub delta_x: Number,
116    pub delta_s: Number,
117    pub delta_c: Number,
118    pub delta_d: Number,
119}
120
121impl PdPerturbationHandler {
122    pub fn new() -> Self {
123        Self::default()
124    }
125
126    /// Configure `perturb_always_cd_` and rebuild the initial `jac`
127    /// state. Mirrors upstream's `InitializeImpl`.
128    pub fn set_perturb_always_cd(&mut self, on: bool) {
129        self.perturb_always_cd = on;
130        self.jac_degenerate = if on {
131            DegenType::NotDegenerate
132        } else {
133            DegenType::NotYetDetermined
134        };
135    }
136
137    /// First call when starting a new aug-system. `mu` is the current
138    /// barrier parameter (used by the `δ_cd` formula).
139    /// Returns `None` if no suitable starting perturbation could be
140    /// found (the caller bails).
141    pub fn consider_new_system(
142        &mut self,
143        mu: Number,
144        ip_data: Option<&IpoptDataHandle>,
145    ) -> Option<Deltas> {
146        self.finalize_test(ip_data);
147
148        // Bookkeeping: roll the previous trial's `_curr` values into
149        // `_last` (matches upstream cpp:158-183).
150        if self.reset_last {
151            self.delta_x_last = self.delta_x_curr;
152            self.delta_s_last = self.delta_s_curr;
153            self.delta_c_last = self.delta_c_curr;
154            self.delta_d_last = self.delta_d_curr;
155        } else {
156            if self.delta_x_curr > 0.0 {
157                self.delta_x_last = self.delta_x_curr;
158            }
159            if self.delta_s_curr > 0.0 {
160                self.delta_s_last = self.delta_s_curr;
161            }
162            if self.delta_c_curr > 0.0 {
163                self.delta_c_last = self.delta_c_curr;
164            }
165            if self.delta_d_curr > 0.0 {
166                self.delta_d_last = self.delta_d_curr;
167            }
168        }
169
170        let undet = matches!(self.hess_degenerate, DegenType::NotYetDetermined)
171            || matches!(self.jac_degenerate, DegenType::NotYetDetermined);
172        self.test_status = if undet {
173            if self.perturb_always_cd {
174                TrialStatus::DcGt0DxEq0
175            } else {
176                TrialStatus::DcEq0DxEq0
177            }
178        } else {
179            TrialStatus::NoTest
180        };
181
182        let mut delta_c = if matches!(self.jac_degenerate, DegenType::Degenerate) {
183            let v = self.delta_cd(mu);
184            self.delta_c_curr = v;
185            append_info(ip_data, "l");
186            v
187        } else if self.perturb_always_cd {
188            let v = self.delta_cd(mu);
189            self.delta_c_curr = v;
190            v
191        } else {
192            self.delta_c_curr = 0.0;
193            0.0
194        };
195        let mut delta_d = delta_c;
196        self.delta_d_curr = delta_d;
197
198        let mut delta_x = 0.0;
199        let mut delta_s = 0.0;
200
201        if matches!(self.hess_degenerate, DegenType::Degenerate) {
202            self.delta_x_curr = 0.0;
203            self.delta_s_curr = 0.0;
204            if !self.get_deltas_for_wrong_inertia(
205                &mut delta_x,
206                &mut delta_s,
207                &mut delta_c,
208                &mut delta_d,
209                ip_data,
210            ) {
211                return None;
212            }
213        }
214
215        self.delta_x_curr = delta_x;
216        self.delta_s_curr = delta_s;
217        self.delta_c_curr = delta_c;
218        self.delta_d_curr = delta_d;
219        set_info_regu_x(ip_data, delta_x);
220        self.get_deltas_for_wrong_inertia_called = false;
221
222        Some(Deltas {
223            delta_x,
224            delta_s,
225            delta_c,
226            delta_d,
227        })
228    }
229
230    /// Escalation after `Singular` factorization status. Mirrors
231    /// `PerturbForSingularity` (cpp:245-364).
232    pub fn perturb_for_singular(
233        &mut self,
234        mu: Number,
235        ip_data: Option<&IpoptDataHandle>,
236    ) -> Option<Deltas> {
237        let mut delta_x = 0.0;
238        let mut delta_s = 0.0;
239        let mut delta_c = 0.0;
240        let mut delta_d = 0.0;
241
242        let undet = matches!(self.hess_degenerate, DegenType::NotYetDetermined)
243            || matches!(self.jac_degenerate, DegenType::NotYetDetermined);
244        if undet {
245            match self.test_status {
246                TrialStatus::DcEq0DxEq0 => {
247                    debug_assert!(self.delta_x_curr == 0.0 && self.delta_c_curr == 0.0);
248                    if matches!(self.jac_degenerate, DegenType::NotYetDetermined) {
249                        let v = self.delta_cd(mu);
250                        self.delta_c_curr = v;
251                        self.delta_d_curr = v;
252                        self.test_status = TrialStatus::DcGt0DxEq0;
253                    } else {
254                        debug_assert!(matches!(self.hess_degenerate, DegenType::NotYetDetermined));
255                        if !self.get_deltas_for_wrong_inertia(
256                            &mut delta_x,
257                            &mut delta_s,
258                            &mut delta_c,
259                            &mut delta_d,
260                            ip_data,
261                        ) {
262                            return None;
263                        }
264                        self.test_status = TrialStatus::DcEq0DxGt0;
265                    }
266                }
267                TrialStatus::DcGt0DxEq0 => {
268                    debug_assert!(self.delta_x_curr == 0.0 && self.delta_c_curr > 0.0);
269                    debug_assert!(matches!(self.jac_degenerate, DegenType::NotYetDetermined));
270                    if !self.perturb_always_cd {
271                        self.delta_c_curr = 0.0;
272                        self.delta_d_curr = 0.0;
273                        if !self.get_deltas_for_wrong_inertia(
274                            &mut delta_x,
275                            &mut delta_s,
276                            &mut delta_c,
277                            &mut delta_d,
278                            ip_data,
279                        ) {
280                            return None;
281                        }
282                        self.test_status = TrialStatus::DcEq0DxGt0;
283                    } else if !self.get_deltas_for_wrong_inertia(
284                        &mut delta_x,
285                        &mut delta_s,
286                        &mut delta_c,
287                        &mut delta_d,
288                        ip_data,
289                    ) {
290                        return None;
291                    } else {
292                        self.test_status = TrialStatus::DcGt0DxGt0;
293                    }
294                }
295                TrialStatus::DcEq0DxGt0 => {
296                    debug_assert!(self.delta_x_curr > 0.0 && self.delta_c_curr == 0.0);
297                    let v = self.delta_cd(mu);
298                    self.delta_c_curr = v;
299                    self.delta_d_curr = v;
300                    if !self.get_deltas_for_wrong_inertia(
301                        &mut delta_x,
302                        &mut delta_s,
303                        &mut delta_c,
304                        &mut delta_d,
305                        ip_data,
306                    ) {
307                        return None;
308                    }
309                    self.test_status = TrialStatus::DcGt0DxGt0;
310                }
311                TrialStatus::DcGt0DxGt0 => {
312                    if !self.get_deltas_for_wrong_inertia(
313                        &mut delta_x,
314                        &mut delta_s,
315                        &mut delta_c,
316                        &mut delta_d,
317                        ip_data,
318                    ) {
319                        return None;
320                    }
321                }
322                TrialStatus::NoTest => {
323                    debug_assert!(false, "perturb_for_singular: NoTest in undetermined branch");
324                }
325            }
326        } else if self.delta_c_curr > 0.0 {
327            // Already perturbed C; treat as wrong-inertia.
328            if !self.get_deltas_for_wrong_inertia(
329                &mut delta_x,
330                &mut delta_s,
331                &mut delta_c,
332                &mut delta_d,
333                ip_data,
334            ) {
335                return None;
336            }
337        } else {
338            let v = self.delta_cd(mu);
339            self.delta_c_curr = v;
340            self.delta_d_curr = v;
341            append_info(ip_data, "L");
342        }
343
344        let out = Deltas {
345            delta_x: self.delta_x_curr,
346            delta_s: self.delta_s_curr,
347            delta_c: self.delta_c_curr,
348            delta_d: self.delta_d_curr,
349        };
350        set_info_regu_x(ip_data, out.delta_x);
351        Some(out)
352    }
353
354    /// Escalation after `WrongInertia` factorization status. Mirrors
355    /// `PerturbForWrongInertia` (cpp:419-450).
356    pub fn perturb_for_wrong_inertia(
357        &mut self,
358        mu: Number,
359        ip_data: Option<&IpoptDataHandle>,
360    ) -> Option<Deltas> {
361        self.finalize_test(ip_data);
362
363        let mut delta_x = 0.0;
364        let mut delta_s = 0.0;
365        let mut delta_c = 0.0;
366        let mut delta_d = 0.0;
367        let mut ok = self.get_deltas_for_wrong_inertia(
368            &mut delta_x,
369            &mut delta_s,
370            &mut delta_c,
371            &mut delta_d,
372            ip_data,
373        );
374        // Upstream "no progress on δ_x but δ_c == 0" recovery: bring
375        // up the C/D perturbation, reset Hessian degeneracy, retry.
376        // Upstream peeks at the OUT-parameter `delta_c`, but
377        // `get_deltas_for_wrong_inertia` only writes that on success;
378        // we look at the handler's own δ_c_curr instead, which
379        // matches the algorithmic intent unambiguously.
380        if !ok && self.delta_c_curr == 0.0 {
381            debug_assert_eq!(self.delta_d_curr, 0.0);
382            let v = self.delta_cd(mu);
383            self.delta_c_curr = v;
384            self.delta_d_curr = v;
385            self.delta_x_curr = 0.0;
386            self.delta_s_curr = 0.0;
387            self.test_status = TrialStatus::NoTest;
388            if matches!(self.hess_degenerate, DegenType::Degenerate) {
389                self.hess_degenerate = DegenType::NotYetDetermined;
390            }
391            ok = self.get_deltas_for_wrong_inertia(
392                &mut delta_x,
393                &mut delta_s,
394                &mut delta_c,
395                &mut delta_d,
396                ip_data,
397            );
398        }
399        if !ok {
400            return None;
401        }
402        Some(Deltas {
403            delta_x,
404            delta_s,
405            delta_c,
406            delta_d,
407        })
408    }
409
410    /// Read the most recently committed perturbations.
411    pub fn current_perturbation(&self) -> Deltas {
412        Deltas {
413            delta_x: self.delta_x_curr,
414            delta_s: self.delta_s_curr,
415            delta_c: self.delta_c_curr,
416            delta_d: self.delta_d_curr,
417        }
418    }
419
420    /// Internal — pure escalation of `δ_x` / `δ_s`. Returns `false` if
421    /// `δ_x` would exceed `delta_xs_max`. Mirrors
422    /// `get_deltas_for_wrong_inertia`.
423    fn get_deltas_for_wrong_inertia(
424        &mut self,
425        delta_x: &mut Number,
426        delta_s: &mut Number,
427        delta_c: &mut Number,
428        delta_d: &mut Number,
429        ip_data: Option<&IpoptDataHandle>,
430    ) -> bool {
431        if self.delta_x_curr == 0.0 {
432            self.delta_x_curr = if self.delta_x_last == 0.0 {
433                self.delta_xs_init
434            } else {
435                self.delta_xs_min
436                    .max(self.delta_x_last * self.delta_xs_dec_fact)
437            };
438        } else if self.delta_x_last == 0.0 || 1e5 * self.delta_x_last < self.delta_x_curr {
439            self.delta_x_curr *= self.delta_xs_first_inc_fact;
440        } else {
441            self.delta_x_curr *= self.delta_xs_inc_fact;
442        }
443        if self.delta_x_curr > self.delta_xs_max {
444            self.delta_x_last = 0.0;
445            self.delta_s_last = 0.0;
446            append_info(ip_data, "dx");
447            return false;
448        }
449        self.delta_s_curr = self.delta_x_curr;
450
451        *delta_x = self.delta_x_curr;
452        *delta_s = self.delta_s_curr;
453        *delta_c = self.delta_c_curr;
454        *delta_d = self.delta_d_curr;
455        set_info_regu_x(ip_data, *delta_x);
456        self.get_deltas_for_wrong_inertia_called = true;
457        true
458    }
459
460    fn delta_cd(&self, mu: Number) -> Number {
461        self.delta_cd_val * mu.powf(self.delta_cd_exp)
462    }
463
464    /// Read the test outcome from the just-completed (non-singular)
465    /// factor and update degeneracy flags. Mirrors `finalize_test`
466    /// (cpp:470-538).
467    fn finalize_test(&mut self, ip_data: Option<&IpoptDataHandle>) {
468        match self.test_status {
469            TrialStatus::NoTest => (),
470            TrialStatus::DcEq0DxEq0 => {
471                if matches!(self.hess_degenerate, DegenType::NotYetDetermined)
472                    && matches!(self.jac_degenerate, DegenType::NotYetDetermined)
473                {
474                    self.hess_degenerate = DegenType::NotDegenerate;
475                    self.jac_degenerate = DegenType::NotDegenerate;
476                    append_info(ip_data, "Nhj ");
477                } else if matches!(self.hess_degenerate, DegenType::NotYetDetermined) {
478                    self.hess_degenerate = DegenType::NotDegenerate;
479                    append_info(ip_data, "Nh ");
480                } else if matches!(self.jac_degenerate, DegenType::NotYetDetermined) {
481                    self.jac_degenerate = DegenType::NotDegenerate;
482                    append_info(ip_data, "Nj ");
483                }
484            }
485            TrialStatus::DcGt0DxEq0 => {
486                if matches!(self.hess_degenerate, DegenType::NotYetDetermined) {
487                    self.hess_degenerate = DegenType::NotDegenerate;
488                    append_info(ip_data, "Nh ");
489                }
490                if matches!(self.jac_degenerate, DegenType::NotYetDetermined) {
491                    self.degen_iters += 1;
492                    if self.degen_iters >= self.degen_iters_max {
493                        self.jac_degenerate = DegenType::Degenerate;
494                        append_info(ip_data, "Dj ");
495                    }
496                    append_info(ip_data, "L");
497                }
498            }
499            TrialStatus::DcEq0DxGt0 => {
500                if matches!(self.jac_degenerate, DegenType::NotYetDetermined) {
501                    self.jac_degenerate = DegenType::NotDegenerate;
502                    append_info(ip_data, "Nj ");
503                }
504                if matches!(self.hess_degenerate, DegenType::NotYetDetermined) {
505                    self.degen_iters += 1;
506                    if self.degen_iters >= self.degen_iters_max {
507                        self.hess_degenerate = DegenType::Degenerate;
508                        append_info(ip_data, "Dh ");
509                    }
510                }
511            }
512            TrialStatus::DcGt0DxGt0 => {
513                self.degen_iters += 1;
514                if self.degen_iters >= self.degen_iters_max {
515                    self.hess_degenerate = DegenType::Degenerate;
516                    self.jac_degenerate = DegenType::Degenerate;
517                    append_info(ip_data, "Dhj ");
518                }
519                append_info(ip_data, "L");
520            }
521        }
522    }
523}
524
525fn append_info(ip_data: Option<&IpoptDataHandle>, s: &str) {
526    if let Some(h) = ip_data {
527        h.borrow_mut().append_info_string(s);
528    }
529}
530
531fn set_info_regu_x(ip_data: Option<&IpoptDataHandle>, v: Number) {
532    if let Some(h) = ip_data {
533        h.borrow_mut().info_regu_x = v;
534    }
535}
536
537#[cfg(test)]
538mod tests {
539    use super::*;
540
541    #[test]
542    fn first_wrong_inertia_perturbation_is_delta_xs_init() {
543        let mut h = PdPerturbationHandler::new();
544        let d = h.perturb_for_wrong_inertia(0.1, None).unwrap();
545        // delta_xs_init = first_hessian_perturbation = 1e-4
546        assert!((d.delta_x - 1e-4).abs() < 1e-20);
547        assert_eq!(d.delta_x, d.delta_s);
548        assert_eq!(d.delta_c, 0.0);
549        assert_eq!(d.delta_d, 0.0);
550    }
551
552    #[test]
553    fn second_perturbation_uses_first_inc_fact() {
554        // After the *first* nonzero δ_x, with δ_x_last == 0, the
555        // doubling uses `delta_xs_first_inc_fact = 100` per upstream
556        // (cpp:386-389: "if delta_x_last_ == 0 ...").
557        let mut h = PdPerturbationHandler::new();
558        let d1 = h.perturb_for_wrong_inertia(0.1, None).unwrap();
559        let d2 = h.perturb_for_wrong_inertia(0.1, None).unwrap();
560        assert!((d2.delta_x - d1.delta_x * 100.0).abs() < 1e-15);
561    }
562
563    #[test]
564    fn third_perturbation_uses_inc_fact() {
565        // After delta_x_last has been set (via consider_new_system or
566        // first inc), continued growth uses `delta_xs_inc_fact = 8`.
567        let mut h = PdPerturbationHandler::new();
568        h.delta_x_curr = 1e-2;
569        h.delta_x_last = 1e-2;
570        let d = h.perturb_for_wrong_inertia(0.1, None).unwrap();
571        assert!((d.delta_x - 1e-2 * 8.0).abs() < 1e-15);
572    }
573
574    #[test]
575    fn perturbation_caps_at_max_when_dcd_already_active() {
576        // When δ_c is already > 0 (e.g., perturb_always_cd, or after a
577        // singular-recovery), the fallback path inside
578        // `perturb_for_wrong_inertia` is skipped and the
579        // δ_x-overflow surfaces as `None`.
580        let mut h = PdPerturbationHandler::new();
581        h.delta_x_curr = h.delta_xs_max;
582        h.delta_c_curr = 1e-4;
583        h.delta_d_curr = 1e-4;
584        assert!(h.perturb_for_wrong_inertia(0.1, None).is_none());
585    }
586
587    #[test]
588    fn consider_new_system_with_perturb_always_cd_seeds_dcd() {
589        let mut h = PdPerturbationHandler::new();
590        h.set_perturb_always_cd(true);
591        let mu = 0.1;
592        let d = h.consider_new_system(mu, None).unwrap();
593        let expected = h.delta_cd_val * mu.powf(h.delta_cd_exp);
594        assert!((d.delta_c - expected).abs() < 1e-15);
595        assert!((d.delta_d - expected).abs() < 1e-15);
596        assert_eq!(d.delta_x, 0.0);
597        assert_eq!(d.delta_s, 0.0);
598    }
599
600    #[test]
601    fn consider_new_system_default_zeros_dcd() {
602        let mut h = PdPerturbationHandler::new();
603        let d = h.consider_new_system(0.1, None).unwrap();
604        assert_eq!(
605            d,
606            Deltas {
607                delta_x: 0.0,
608                delta_s: 0.0,
609                delta_c: 0.0,
610                delta_d: 0.0
611            }
612        );
613    }
614
615    #[test]
616    fn singular_in_test_dc_eq0_dx_eq0_seeds_dcd() {
617        let mut h = PdPerturbationHandler::new();
618        let _ = h.consider_new_system(0.1, None).unwrap();
619        // After consider_new_system on a fresh handler, test_status
620        // should be DcEq0DxEq0 (since both flags are NotYetDetermined,
621        // and perturb_always_cd is false).
622        assert_eq!(h.test_status, TrialStatus::DcEq0DxEq0);
623        let d = h.perturb_for_singular(0.1, None).unwrap();
624        let expected = h.delta_cd_val * (0.1_f64).powf(h.delta_cd_exp);
625        assert!((d.delta_c - expected).abs() < 1e-15);
626        assert!((d.delta_d - expected).abs() < 1e-15);
627        assert_eq!(d.delta_x, 0.0);
628        assert_eq!(h.test_status, TrialStatus::DcGt0DxEq0);
629    }
630
631    #[test]
632    fn singular_when_determined_with_dc_zero_seeds_dcd() {
633        let mut h = PdPerturbationHandler::new();
634        h.hess_degenerate = DegenType::NotDegenerate;
635        h.jac_degenerate = DegenType::NotDegenerate;
636        h.test_status = TrialStatus::NoTest;
637        let d = h.perturb_for_singular(0.1, None).unwrap();
638        let expected = h.delta_cd_val * (0.1_f64).powf(h.delta_cd_exp);
639        assert!((d.delta_c - expected).abs() < 1e-15);
640    }
641
642    #[test]
643    fn finalize_test_sets_not_degenerate_after_dc_eq0_dx_eq0_pass() {
644        let mut h = PdPerturbationHandler::new();
645        let _ = h.consider_new_system(0.1, None).unwrap();
646        // Simulate the next call recognizing that the previous trial
647        // factor was non-singular: a fresh consider_new_system runs
648        // finalize_test first.
649        let _ = h.consider_new_system(0.1, None).unwrap();
650        assert_eq!(h.hess_degenerate, DegenType::NotDegenerate);
651        assert_eq!(h.jac_degenerate, DegenType::NotDegenerate);
652    }
653
654    #[test]
655    fn current_perturbation_returns_committed_values() {
656        let mut h = PdPerturbationHandler::new();
657        h.delta_x_curr = 1.0;
658        h.delta_s_curr = 2.0;
659        h.delta_c_curr = 3.0;
660        h.delta_d_curr = 4.0;
661        let d = h.current_perturbation();
662        assert_eq!(d.delta_x, 1.0);
663        assert_eq!(d.delta_s, 2.0);
664        assert_eq!(d.delta_c, 3.0);
665        assert_eq!(d.delta_d, 4.0);
666    }
667}
pounce_algorithm/kkt/perturbation_handler.rs

pounce_algorithm/kkt/
perturbation_handler.rs
