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        if std::env::var_os("POUNCE_DBG_PERT").is_some() {
362            let it = ip_data.map(|d| d.borrow().iter_count).unwrap_or(-1);
363            tracing::debug!(target: "pounce::linsol",
364                "[PERT] iter={} WRONG_INERTIA mu={:.2e} dx_last={:.2e} dx_curr={:.2e}",
365                it, mu, self.delta_x_last, self.delta_x_curr
366            );
367        }
368        self.finalize_test(ip_data);
369
370        let mut delta_x = 0.0;
371        let mut delta_s = 0.0;
372        let mut delta_c = 0.0;
373        let mut delta_d = 0.0;
374        let mut ok = self.get_deltas_for_wrong_inertia(
375            &mut delta_x,
376            &mut delta_s,
377            &mut delta_c,
378            &mut delta_d,
379            ip_data,
380        );
381        // Upstream "no progress on δ_x but δ_c == 0" recovery: bring
382        // up the C/D perturbation, reset Hessian degeneracy, retry.
383        // Upstream peeks at the OUT-parameter `delta_c`, but
384        // `get_deltas_for_wrong_inertia` only writes that on success;
385        // we look at the handler's own δ_c_curr instead, which
386        // matches the algorithmic intent unambiguously.
387        if !ok && self.delta_c_curr == 0.0 {
388            debug_assert_eq!(self.delta_d_curr, 0.0);
389            let v = self.delta_cd(mu);
390            self.delta_c_curr = v;
391            self.delta_d_curr = v;
392            self.delta_x_curr = 0.0;
393            self.delta_s_curr = 0.0;
394            self.test_status = TrialStatus::NoTest;
395            if matches!(self.hess_degenerate, DegenType::Degenerate) {
396                self.hess_degenerate = DegenType::NotYetDetermined;
397            }
398            ok = self.get_deltas_for_wrong_inertia(
399                &mut delta_x,
400                &mut delta_s,
401                &mut delta_c,
402                &mut delta_d,
403                ip_data,
404            );
405        }
406        if !ok {
407            return None;
408        }
409        Some(Deltas {
410            delta_x,
411            delta_s,
412            delta_c,
413            delta_d,
414        })
415    }
416
417    /// Read the most recently committed perturbations.
418    pub fn current_perturbation(&self) -> Deltas {
419        Deltas {
420            delta_x: self.delta_x_curr,
421            delta_s: self.delta_s_curr,
422            delta_c: self.delta_c_curr,
423            delta_d: self.delta_d_curr,
424        }
425    }
426
427    /// Internal — pure escalation of `δ_x` / `δ_s`. Returns `false` if
428    /// `δ_x` would exceed `delta_xs_max`. Mirrors
429    /// `get_deltas_for_wrong_inertia`.
430    fn get_deltas_for_wrong_inertia(
431        &mut self,
432        delta_x: &mut Number,
433        delta_s: &mut Number,
434        delta_c: &mut Number,
435        delta_d: &mut Number,
436        ip_data: Option<&IpoptDataHandle>,
437    ) -> bool {
438        if self.delta_x_curr == 0.0 {
439            self.delta_x_curr = if self.delta_x_last == 0.0 {
440                self.delta_xs_init
441            } else {
442                self.delta_xs_min
443                    .max(self.delta_x_last * self.delta_xs_dec_fact)
444            };
445        } else if self.delta_x_last == 0.0 || 1e5 * self.delta_x_last < self.delta_x_curr {
446            self.delta_x_curr *= self.delta_xs_first_inc_fact;
447        } else {
448            self.delta_x_curr *= self.delta_xs_inc_fact;
449        }
450        if self.delta_x_curr > self.delta_xs_max {
451            self.delta_x_last = 0.0;
452            self.delta_s_last = 0.0;
453            append_info(ip_data, "dx");
454            return false;
455        }
456        self.delta_s_curr = self.delta_x_curr;
457
458        *delta_x = self.delta_x_curr;
459        *delta_s = self.delta_s_curr;
460        *delta_c = self.delta_c_curr;
461        *delta_d = self.delta_d_curr;
462        set_info_regu_x(ip_data, *delta_x);
463        self.get_deltas_for_wrong_inertia_called = true;
464        true
465    }
466
467    fn delta_cd(&self, mu: Number) -> Number {
468        self.delta_cd_val * mu.powf(self.delta_cd_exp)
469    }
470
471    /// Read the test outcome from the just-completed (non-singular)
472    /// factor and update degeneracy flags. Mirrors `finalize_test`
473    /// (cpp:470-538).
474    fn finalize_test(&mut self, ip_data: Option<&IpoptDataHandle>) {
475        match self.test_status {
476            TrialStatus::NoTest => (),
477            TrialStatus::DcEq0DxEq0 => {
478                if matches!(self.hess_degenerate, DegenType::NotYetDetermined)
479                    && matches!(self.jac_degenerate, DegenType::NotYetDetermined)
480                {
481                    self.hess_degenerate = DegenType::NotDegenerate;
482                    self.jac_degenerate = DegenType::NotDegenerate;
483                    append_info(ip_data, "Nhj ");
484                } else if matches!(self.hess_degenerate, DegenType::NotYetDetermined) {
485                    self.hess_degenerate = DegenType::NotDegenerate;
486                    append_info(ip_data, "Nh ");
487                } else if matches!(self.jac_degenerate, DegenType::NotYetDetermined) {
488                    self.jac_degenerate = DegenType::NotDegenerate;
489                    append_info(ip_data, "Nj ");
490                }
491            }
492            TrialStatus::DcGt0DxEq0 => {
493                if matches!(self.hess_degenerate, DegenType::NotYetDetermined) {
494                    self.hess_degenerate = DegenType::NotDegenerate;
495                    append_info(ip_data, "Nh ");
496                }
497                if matches!(self.jac_degenerate, DegenType::NotYetDetermined) {
498                    self.degen_iters += 1;
499                    if self.degen_iters >= self.degen_iters_max {
500                        self.jac_degenerate = DegenType::Degenerate;
501                        append_info(ip_data, "Dj ");
502                    }
503                    append_info(ip_data, "L");
504                }
505            }
506            TrialStatus::DcEq0DxGt0 => {
507                if matches!(self.jac_degenerate, DegenType::NotYetDetermined) {
508                    self.jac_degenerate = DegenType::NotDegenerate;
509                    append_info(ip_data, "Nj ");
510                }
511                if matches!(self.hess_degenerate, DegenType::NotYetDetermined) {
512                    self.degen_iters += 1;
513                    if self.degen_iters >= self.degen_iters_max {
514                        self.hess_degenerate = DegenType::Degenerate;
515                        append_info(ip_data, "Dh ");
516                    }
517                }
518            }
519            TrialStatus::DcGt0DxGt0 => {
520                self.degen_iters += 1;
521                if self.degen_iters >= self.degen_iters_max {
522                    self.hess_degenerate = DegenType::Degenerate;
523                    self.jac_degenerate = DegenType::Degenerate;
524                    append_info(ip_data, "Dhj ");
525                }
526                append_info(ip_data, "L");
527            }
528        }
529    }
530}
531
532fn append_info(ip_data: Option<&IpoptDataHandle>, s: &str) {
533    if let Some(h) = ip_data {
534        h.borrow_mut().append_info_string(s);
535    }
536}
537
538fn set_info_regu_x(ip_data: Option<&IpoptDataHandle>, v: Number) {
539    if let Some(h) = ip_data {
540        h.borrow_mut().info_regu_x = v;
541    }
542}
543
544#[cfg(test)]
545mod tests {
546    use super::*;
547
548    #[test]
549    fn first_wrong_inertia_perturbation_is_delta_xs_init() {
550        let mut h = PdPerturbationHandler::new();
551        let d = h.perturb_for_wrong_inertia(0.1, None).unwrap();
552        // delta_xs_init = first_hessian_perturbation = 1e-4
553        assert!((d.delta_x - 1e-4).abs() < 1e-20);
554        assert_eq!(d.delta_x, d.delta_s);
555        assert_eq!(d.delta_c, 0.0);
556        assert_eq!(d.delta_d, 0.0);
557    }
558
559    #[test]
560    fn second_perturbation_uses_first_inc_fact() {
561        // After the *first* nonzero δ_x, with δ_x_last == 0, the
562        // doubling uses `delta_xs_first_inc_fact = 100` per upstream
563        // (cpp:386-389: "if delta_x_last_ == 0 ...").
564        let mut h = PdPerturbationHandler::new();
565        let d1 = h.perturb_for_wrong_inertia(0.1, None).unwrap();
566        let d2 = h.perturb_for_wrong_inertia(0.1, None).unwrap();
567        assert!((d2.delta_x - d1.delta_x * 100.0).abs() < 1e-15);
568    }
569
570    #[test]
571    fn third_perturbation_uses_inc_fact() {
572        // After delta_x_last has been set (via consider_new_system or
573        // first inc), continued growth uses `delta_xs_inc_fact = 8`.
574        let mut h = PdPerturbationHandler::new();
575        h.delta_x_curr = 1e-2;
576        h.delta_x_last = 1e-2;
577        let d = h.perturb_for_wrong_inertia(0.1, None).unwrap();
578        assert!((d.delta_x - 1e-2 * 8.0).abs() < 1e-15);
579    }
580
581    #[test]
582    fn perturbation_caps_at_max_when_dcd_already_active() {
583        // When δ_c is already > 0 (e.g., perturb_always_cd, or after a
584        // singular-recovery), the fallback path inside
585        // `perturb_for_wrong_inertia` is skipped and the
586        // δ_x-overflow surfaces as `None`.
587        let mut h = PdPerturbationHandler::new();
588        h.delta_x_curr = h.delta_xs_max;
589        h.delta_c_curr = 1e-4;
590        h.delta_d_curr = 1e-4;
591        assert!(h.perturb_for_wrong_inertia(0.1, None).is_none());
592    }
593
594    #[test]
595    fn consider_new_system_with_perturb_always_cd_seeds_dcd() {
596        let mut h = PdPerturbationHandler::new();
597        h.set_perturb_always_cd(true);
598        let mu = 0.1;
599        let d = h.consider_new_system(mu, None).unwrap();
600        let expected = h.delta_cd_val * mu.powf(h.delta_cd_exp);
601        assert!((d.delta_c - expected).abs() < 1e-15);
602        assert!((d.delta_d - expected).abs() < 1e-15);
603        assert_eq!(d.delta_x, 0.0);
604        assert_eq!(d.delta_s, 0.0);
605    }
606
607    #[test]
608    fn consider_new_system_default_zeros_dcd() {
609        let mut h = PdPerturbationHandler::new();
610        let d = h.consider_new_system(0.1, None).unwrap();
611        assert_eq!(
612            d,
613            Deltas {
614                delta_x: 0.0,
615                delta_s: 0.0,
616                delta_c: 0.0,
617                delta_d: 0.0
618            }
619        );
620    }
621
622    #[test]
623    fn singular_in_test_dc_eq0_dx_eq0_seeds_dcd() {
624        let mut h = PdPerturbationHandler::new();
625        let _ = h.consider_new_system(0.1, None).unwrap();
626        // After consider_new_system on a fresh handler, test_status
627        // should be DcEq0DxEq0 (since both flags are NotYetDetermined,
628        // and perturb_always_cd is false).
629        assert_eq!(h.test_status, TrialStatus::DcEq0DxEq0);
630        let d = h.perturb_for_singular(0.1, None).unwrap();
631        let expected = h.delta_cd_val * (0.1_f64).powf(h.delta_cd_exp);
632        assert!((d.delta_c - expected).abs() < 1e-15);
633        assert!((d.delta_d - expected).abs() < 1e-15);
634        assert_eq!(d.delta_x, 0.0);
635        assert_eq!(h.test_status, TrialStatus::DcGt0DxEq0);
636    }
637
638    #[test]
639    fn singular_when_determined_with_dc_zero_seeds_dcd() {
640        let mut h = PdPerturbationHandler::new();
641        h.hess_degenerate = DegenType::NotDegenerate;
642        h.jac_degenerate = DegenType::NotDegenerate;
643        h.test_status = TrialStatus::NoTest;
644        let d = h.perturb_for_singular(0.1, None).unwrap();
645        let expected = h.delta_cd_val * (0.1_f64).powf(h.delta_cd_exp);
646        assert!((d.delta_c - expected).abs() < 1e-15);
647    }
648
649    #[test]
650    fn finalize_test_sets_not_degenerate_after_dc_eq0_dx_eq0_pass() {
651        let mut h = PdPerturbationHandler::new();
652        let _ = h.consider_new_system(0.1, None).unwrap();
653        // Simulate the next call recognizing that the previous trial
654        // factor was non-singular: a fresh consider_new_system runs
655        // finalize_test first.
656        let _ = h.consider_new_system(0.1, None).unwrap();
657        assert_eq!(h.hess_degenerate, DegenType::NotDegenerate);
658        assert_eq!(h.jac_degenerate, DegenType::NotDegenerate);
659    }
660
661    #[test]
662    fn current_perturbation_returns_committed_values() {
663        let mut h = PdPerturbationHandler::new();
664        h.delta_x_curr = 1.0;
665        h.delta_s_curr = 2.0;
666        h.delta_c_curr = 3.0;
667        h.delta_d_curr = 4.0;
668        let d = h.current_perturbation();
669        assert_eq!(d.delta_x, 1.0);
670        assert_eq!(d.delta_s, 2.0);
671        assert_eq!(d.delta_c, 3.0);
672        assert_eq!(d.delta_d, 4.0);
673    }
674}
pounce_algorithm/kkt/perturbation_handler.rs

pounce_algorithm/kkt/
perturbation_handler.rs
