clarabel/solver/implementations/default/
info.rs

1use super::*;
2use crate::algebra::*;
3use crate::io::PrintTarget;
4use crate::solver::core::ffi::*;
5use crate::solver::core::kktsolvers::LinearSolverInfo;
6use crate::solver::core::{traits::Info, SolverStatus};
7use crate::solver::traits::Variables;
8use crate::timers::*;
9
10/// Standard-form solver type implementing the [`Info`](crate::solver::core::traits::Info) and [`InfoPrint`](crate::solver::core::traits::InfoPrint) traits
11#[repr(C)]
12#[derive(Default, Debug, Clone)]
13pub struct DefaultInfo<T> {
14    /// interior point path parameter μ
15    pub mu: T,
16    /// interior point path parameter reduction ratio σ
17    pub sigma: T,
18    /// step length for the current iteration
19    pub step_length: T,
20    /// number of iterations
21    pub iterations: u32,
22    /// primal objective value
23    pub cost_primal: T,
24    /// dual objective value
25    pub cost_dual: T,
26    /// primal residual
27    pub res_primal: T,
28    /// dual residual
29    pub res_dual: T,
30    /// primal infeasibility residual
31    pub res_primal_inf: T,
32    /// dual infeasibility residual
33    pub res_dual_inf: T,
34    /// absolute duality gap
35    pub gap_abs: T,
36    /// relative duality gap
37    pub gap_rel: T,
38    /// κ/τ ratio
39    pub ktratio: T,
40
41    // previous iterate
42    /// primal object value from previous iteration
43    pub(crate) prev_cost_primal: T,
44    /// dual objective value from previous iteration
45    pub(crate) prev_cost_dual: T,
46    /// primal residual from previous iteration
47    pub(crate) prev_res_primal: T,
48    /// dual residual from previous iteration
49    pub(crate) prev_res_dual: T,
50    /// absolute duality gap from previous iteration
51    pub(crate) prev_gap_abs: T,
52    /// relative duality gap from previous iteration
53    pub(crate) prev_gap_rel: T,
54    /// solve time
55    pub solve_time: f64,
56    /// solver status
57    pub status: SolverStatus,
58
59    /// linear solver information
60    pub linsolver: LinearSolverInfo,
61
62    // target stream for printing
63    pub(crate) stream: PrintTarget,
64}
65
66impl<T> DefaultInfo<T>
67where
68    T: FloatT,
69{
70    /// creates a new `DefaultInfo` object
71    pub fn new() -> Self {
72        Self::default()
73    }
74}
75
76impl<T: FloatT> ClarabelFFI<Self> for DefaultInfo<T> {
77    type FFI = super::ffi::DefaultInfoFFI<T>;
78}
79
80impl<T> Info<T> for DefaultInfo<T>
81where
82    T: FloatT,
83{
84    type V = DefaultVariables<T>;
85    type R = DefaultResiduals<T>;
86
87    fn reset(&mut self, timers: &mut Timers) {
88        self.status = SolverStatus::Unsolved;
89        self.iterations = 0;
90        self.solve_time = 0f64;
91
92        timers.reset_timer("solve");
93    }
94
95    fn post_process(&mut self, residuals: &DefaultResiduals<T>, settings: &DefaultSettings<T>) {
96        // if there was an error or we ran out of time
97        // or iterations, check for partial convergence
98
99        if self.status.is_errored()
100            || matches!(self.status, SolverStatus::MaxIterations)
101            || matches!(self.status, SolverStatus::MaxTime)
102        {
103            self.check_convergence_almost(residuals, settings);
104        }
105    }
106
107    fn finalize(&mut self, timers: &mut Timers) {
108        //final check of timers
109        self.solve_time = timers.total_time().as_secs_f64();
110    }
111
112    fn update(
113        &mut self,
114        data: &mut DefaultProblemData<T>,
115        variables: &DefaultVariables<T>,
116        residuals: &DefaultResiduals<T>,
117        timers: &Timers,
118    ) {
119        // optimality termination check should be computed w.r.t
120        // the pre-homogenization x and z variables.
121        let τinv = T::recip(variables.τ);
122
123        // unscaled linear term norms
124        let normb = data.get_normb();
125        let normq = data.get_normq();
126
127        // shortcuts for the equilibration matrices
128        let d = &data.equilibration.d;
129        let e = &data.equilibration.e;
130        let dinv = &data.equilibration.dinv;
131        let einv = &data.equilibration.einv;
132        let cinv = T::recip(data.equilibration.c);
133
134        // primal and dual costs. dot products are invariant w.r.t
135        // equilibration, but we still need to back out the overall
136        // objective scaling term c
137
138        let xPx_τinvsq_over2 = residuals.dot_xPx * τinv * τinv / (2.).as_T();
139        self.cost_primal = (residuals.dot_qx * τinv + xPx_τinvsq_over2) * cinv;
140        self.cost_dual = (-residuals.dot_bz * τinv - xPx_τinvsq_over2) * cinv;
141
142        // variables norms, undoing the equilibration.  Do not unscale
143        // by τ yet because the infeasibility residuals are ratios of
144        // terms that have no affine parts anyway
145        let mut normx = variables.x.norm_scaled(d);
146        let mut normz = variables.z.norm_scaled(e) * cinv;
147        let mut norms = variables.s.norm_scaled(einv);
148
149        // primal and dual infeasibility residuals.
150        self.res_primal_inf = (residuals.rx_inf.norm_scaled(dinv) * cinv) / T::max(T::one(), normz);
151        self.res_dual_inf = T::max(
152            residuals.Px.norm_scaled(dinv) / T::max(T::one(), normx),
153            residuals.rz_inf.norm_scaled(einv) / T::max(T::one(), normx + norms),
154        );
155
156        // now back out the τ scaling so we can normalize the unscaled primal / dual errors
157        normx *= τinv;
158        normz *= τinv;
159        norms *= τinv;
160
161        // primal and dual relative residuals.
162        self.res_primal =
163            residuals.rz.norm_scaled(einv) * τinv / T::max(T::one(), normb + normx + norms);
164        self.res_dual =
165            residuals.rx.norm_scaled(dinv) * τinv * cinv / T::max(T::one(), normq + normx + normz);
166
167        // absolute and relative gaps
168        self.gap_abs = T::abs(self.cost_primal - self.cost_dual);
169        self.gap_rel = self.gap_abs
170            / T::max(
171                T::one(),
172                T::min(T::abs(self.cost_primal), T::abs(self.cost_dual)),
173            );
174
175        // κ/τ ratio (scaled)
176        self.ktratio = variables.κ * τinv;
177
178        // solve time so far (includes setup)
179        self.solve_time = timers.total_time().as_secs_f64();
180    }
181
182    fn check_termination(
183        &mut self,
184        residuals: &DefaultResiduals<T>,
185        settings: &DefaultSettings<T>,
186        iter: u32,
187    ) -> bool {
188        //  optimality or infeasibility
189        // ---------------------
190        self.check_convergence_full(residuals, settings);
191
192        //  poor progress
193        // ----------------------
194        if self.status == SolverStatus::Unsolved
195            && iter > 1u32
196            && (self.res_dual > self.prev_res_dual || self.res_primal > self.prev_res_primal)
197        {
198            // Poor progress at high tolerance.
199            if self.ktratio < T::epsilon() * (100.).as_T()
200                && (self.prev_gap_abs < settings.tol_gap_abs
201                    || self.prev_gap_rel < settings.tol_gap_rel)
202            {
203                self.status = SolverStatus::InsufficientProgress;
204            }
205
206            // Going backwards. Stop immediately if residuals diverge out of feasibility tolerance.
207            #[allow(clippy::collapsible_if)] // nested if for readability
208            if self.ktratio < T::one() {
209                if (self.res_dual > settings.tol_feas * (100.).as_T()
210                    && self.res_dual > self.prev_res_dual * (100.).as_T())
211                    || (self.res_primal > settings.tol_feas * (100.).as_T()
212                        && self.res_primal > self.prev_res_primal * (100.).as_T())
213                {
214                    self.status = SolverStatus::InsufficientProgress;
215                }
216            }
217        }
218
219        // time or iteration limits
220        // ----------------------
221        if self.status == SolverStatus::Unsolved {
222            if settings.max_iter == self.iterations {
223                self.status = SolverStatus::MaxIterations;
224            } else if self.solve_time > settings.time_limit {
225                self.status = SolverStatus::MaxTime;
226            }
227        }
228
229        // return TRUE if we settled on a final status
230        self.status != SolverStatus::Unsolved
231    }
232
233    fn save_prev_iterate(&mut self, variables: &Self::V, prev_variables: &mut Self::V) {
234        self.prev_cost_primal = self.cost_primal;
235        self.prev_cost_dual = self.cost_dual;
236        self.prev_res_primal = self.res_primal;
237        self.prev_res_dual = self.res_dual;
238        self.prev_gap_abs = self.gap_abs;
239        self.prev_gap_rel = self.gap_rel;
240
241        prev_variables.copy_from(variables);
242    }
243
244    fn reset_to_prev_iterate(&mut self, variables: &mut Self::V, prev_variables: &Self::V) {
245        self.cost_primal = self.prev_cost_primal;
246        self.cost_dual = self.prev_cost_dual;
247        self.res_primal = self.prev_res_primal;
248        self.res_dual = self.prev_res_dual;
249        self.gap_abs = self.prev_gap_abs;
250        self.gap_rel = self.prev_gap_rel;
251
252        variables.copy_from(prev_variables);
253    }
254
255    fn save_scalars(&mut self, μ: T, α: T, σ: T, iter: u32) {
256        self.mu = μ;
257        self.step_length = α;
258        self.sigma = σ;
259        self.iterations = iter;
260    }
261
262    fn get_status(&self) -> SolverStatus {
263        self.status
264    }
265
266    fn set_status(&mut self, status: SolverStatus) {
267        self.status = status;
268    }
269}
270
271// Utility functions for convergence checkiing
272
273impl<T> DefaultInfo<T>
274where
275    T: FloatT,
276{
277    fn check_convergence_full(
278        &mut self,
279        residuals: &DefaultResiduals<T>,
280        settings: &DefaultSettings<T>,
281    ) {
282        // "full" tolerances
283        let tol_gap_abs = settings.tol_gap_abs;
284        let tol_gap_rel = settings.tol_gap_rel;
285        let tol_feas = settings.tol_feas;
286        let tol_infeas_abs = settings.tol_infeas_abs;
287        let tol_infeas_rel = settings.tol_infeas_rel;
288        let tol_ktratio = settings.tol_ktratio;
289
290        let solved_status = SolverStatus::Solved;
291        let pinf_status = SolverStatus::PrimalInfeasible;
292        let dinf_status = SolverStatus::DualInfeasible;
293
294        self.check_convergence(
295            residuals,
296            tol_gap_abs,
297            tol_gap_rel,
298            tol_feas,
299            tol_infeas_abs,
300            tol_infeas_rel,
301            tol_ktratio,
302            solved_status,
303            pinf_status,
304            dinf_status,
305        );
306    }
307
308    fn check_convergence_almost(
309        &mut self,
310        residuals: &DefaultResiduals<T>,
311        settings: &DefaultSettings<T>,
312    ) {
313        // "almost" tolerances
314        let tol_gap_abs = settings.reduced_tol_gap_abs;
315        let tol_gap_rel = settings.reduced_tol_gap_rel;
316        let tol_feas = settings.reduced_tol_feas;
317        let tol_infeas_abs = settings.reduced_tol_infeas_abs;
318        let tol_infeas_rel = settings.reduced_tol_infeas_rel;
319        let tol_ktratio = settings.reduced_tol_ktratio;
320
321        let solved_status = SolverStatus::AlmostSolved;
322        let pinf_status = SolverStatus::AlmostPrimalInfeasible;
323        let dinf_status = SolverStatus::AlmostDualInfeasible;
324
325        self.check_convergence(
326            residuals,
327            tol_gap_abs,
328            tol_gap_rel,
329            tol_feas,
330            tol_infeas_abs,
331            tol_infeas_rel,
332            tol_ktratio,
333            solved_status,
334            pinf_status,
335            dinf_status,
336        );
337    }
338
339    #[allow(clippy::too_many_arguments)]
340    fn check_convergence(
341        &mut self,
342        residuals: &DefaultResiduals<T>,
343        tol_gap_abs: T,
344        tol_gap_rel: T,
345        tol_feas: T,
346        tol_infeas_abs: T,
347        tol_infeas_rel: T,
348        tol_ktratio: T,
349        solved_status: SolverStatus,
350        pinf_status: SolverStatus,
351        dinf_status: SolverStatus,
352    ) {
353        if self.ktratio <= T::one() && self.is_solved(tol_gap_abs, tol_gap_rel, tol_feas) {
354            self.status = solved_status;
355        //PJG hardcoded factor 1000 here should be fixed
356        } else if self.ktratio > tol_ktratio.recip() * (1000.0).as_T() {
357            if self.is_primal_infeasible(residuals, tol_infeas_abs, tol_infeas_rel) {
358                self.status = pinf_status;
359            } else if self.is_dual_infeasible(residuals, tol_infeas_abs, tol_infeas_rel) {
360                self.status = dinf_status;
361            }
362        }
363    }
364
365    fn is_solved(&self, tol_gap_abs: T, tol_gap_rel: T, tol_feas: T) -> bool {
366        ((self.gap_abs < tol_gap_abs) || (self.gap_rel < tol_gap_rel))
367            && (self.res_primal < tol_feas)
368            && (self.res_dual < tol_feas)
369    }
370
371    fn is_primal_infeasible(
372        &self,
373        residuals: &DefaultResiduals<T>,
374        tol_infeas_abs: T,
375        tol_infeas_rel: T,
376    ) -> bool {
377        (residuals.dot_bz < -tol_infeas_abs)
378            && (self.res_primal_inf < -tol_infeas_rel * residuals.dot_bz)
379    }
380
381    fn is_dual_infeasible(
382        &self,
383        residuals: &DefaultResiduals<T>,
384        tol_infeas_abs: T,
385        tol_infeas_rel: T,
386    ) -> bool {
387        (residuals.dot_qx < -tol_infeas_abs)
388            && (self.res_dual_inf < -tol_infeas_rel * residuals.dot_qx)
389    }
390}
clarabel/solver/implementations/default/info.rs

clarabel/solver/implementations/default/
info.rs
