OpEn: good error handling

examples/panoc_ex1.rs

@@ -29,7 +29,9 @@ fn main() {
     // define the cost function and its gradient
     let df = |u: &[f64], grad: &mut [f64]| -> Result<(), SolverError> {
         if a < 0.0 || b < 0.0 {
-            Err(SolverError::Cost)
+            Err(SolverError::Cost(
+                "Rosenbrock parameters must be nonnegative",
+            ))
         } else {
             rosenbrock_grad(a, b, u, grad);
             Ok(())
@@ -38,7 +40,9 @@ fn main() {
 
     let f = |u: &[f64], c: &mut f64| -> Result<(), SolverError> {
         if a < 0.0 || b < 0.0 {
-            Err(SolverError::Cost)
+            Err(SolverError::Cost(
+                "Rosenbrock parameters must be nonnegative",
+            ))
         } else {
             *c = rosenbrock_cost(a, b, u);
             Ok(())

open-codegen/opengen/templates/c/optimizer_cinterface.rs.jinja

@@ -179,8 +179,11 @@ pub unsafe extern "C" fn {{meta.optimizer_name|lower}}_solve(
         },
         Err(e) => {{meta.optimizer_name}}SolverStatus {
             exit_status: match e {
-                SolverError::Cost => {{meta.optimizer_name}}ExitStatus::{{meta.optimizer_name}}NotConvergedCost,
-                SolverError::NotFiniteComputation => {{meta.optimizer_name}}ExitStatus::{{meta.optimizer_name}}NotConvergedNotFiniteComputation,
+                SolverError::Cost(_)
+                | SolverError::ProjectionFailed(_)
+                | SolverError::LinearAlgebraFailure(_)
+                | SolverError::InvalidProblemState(_) => {{meta.optimizer_name}}ExitStatus::{{meta.optimizer_name}}NotConvergedCost,
+                SolverError::NotFiniteComputation(_) => {{meta.optimizer_name}}ExitStatus::{{meta.optimizer_name}}NotConvergedNotFiniteComputation,
             },
             num_outer_iterations: u64::MAX as c_ulong,
             num_inner_iterations: u64::MAX as c_ulong,
@@ -209,4 +212,4 @@ pub unsafe extern "C" fn {{meta.optimizer_name|lower}}_free(instance: *mut {{met
     assert!(!instance.is_null());
     drop(Box::from_raw(instance));
 }
-{% endif %}
+{% endif %}

open-codegen/opengen/templates/tcp/tcp_server.rs

@@ -210,7 +210,7 @@ fn execution_handler(
             return_solution_to_client(ok_status, u, stream);
         }
         Err(err) => {
-            let error_message = format!("problem solution failed: {:?}", err);
+            let error_message = format!("problem solution failed: {}", err);
             write_error_message(stream, 2000, &error_message);
         }
     }

open-codegen/test/test.py

@@ -267,7 +267,7 @@ def setUpSolverError(cls):
             anchor +
             '\n'
             '    if p[0] < 0.0 {\n'
-            '        return Err(SolverError::Cost);\n'
+            '        return Err(SolverError::Cost("forced solver error for TCP test"));\n'
             '    }\n'
         )
         if anchor not in solver_lib:
@@ -544,7 +544,9 @@ def test_rust_build_solver_error_details(self):
         self.assertFalse(response.is_ok())
         status = response.get()
         self.assertEqual(2000, status.code)
-        self.assertEqual("problem solution failed: Cost", status.message)
+        self.assertEqual(
+            "problem solution failed: cost or gradient evaluation failed: forced solver error for TCP test",
+            status.message)
 
     def test_rust_build_parametric_f2(self):
         # introduced to tackle issue #123

src/alm/alm_factory.rs

@@ -155,16 +155,16 @@
     /// );
     /// ```
     ///
     pub fn new(
         f: Cost,
         df: CostGradient,
         mapping_f1: Option<MappingF1>,
         jacobian_mapping_f1_trans: Option<JacobianMappingF1Trans>,
         mapping_f2: Option<MappingF2>,
         jacobian_mapping_f2_trans: Option<JacobianMappingF2Trans>,
         set_c: Option<SetC>,
         n2: usize,
     ) -> Self {
         assert!(
             !(mapping_f2.is_none() ^ (n2 == 0)),
             "if n2 > 0 then and only then should you provide an F2"
@@ -233,7 +233,7 @@
                 .zip(y_lagrange_mult.iter())
                 .for_each(|(ti, yi)| *ti += yi / f64::max(penalty_parameter, 1.0));
             s.copy_from_slice(&f1_u_plus_y_over_c);
-            set_c.project(&mut s);
+            set_c.project(&mut s)?;
             *cost += 0.5
                 * penalty_parameter
                 * matrix_operations::norm2_squared_diff(&f1_u_plus_y_over_c, &s);
@@ -296,7 +296,7 @@
                 .zip(y_lagrange_mult.iter())
                 .for_each(|(ti, yi)| *ti += yi / c_penalty_parameter);
             s_aux_var.copy_from_slice(&f1_u_plus_y_over_c); // s = t
-            set_c.project(&mut s_aux_var); // s = Proj_C(F1(u) + y/c)
+            set_c.project(&mut s_aux_var)?; // s = Proj_C(F1(u) + y/c)
 
             // t = F1(u) + y/c - Proj_C(F1(u) + y/c)
             f1_u_plus_y_over_c
@@ -412,7 +412,9 @@
         let f2 = mapping_f2;
         let jac_f2_tr =
             |_u: &[f64], _d: &[f64], _res: &mut [f64]| -> Result<(), crate::SolverError> {
-                Err(SolverError::NotFiniteComputation)
+                Err(SolverError::NotFiniteComputation(
+                    "mock Jacobian-transpose product returned a non-finite result",
+                ))
             };
         let factory = AlmFactory::new(
             mocks::f0,

src/alm/alm_optimizer.rs

@@ -650,7 +650,7 @@ where
                 .for_each(|((y_plus_i, y_i), w_alm_aux_i)| *y_plus_i = w_alm_aux_i + y_i / c);
 
             // Step #3: y_plus := Proj_C(y_plus)
-            alm_set_c.project(y_plus);
+            alm_set_c.project(y_plus)?;
 
             // Step #4
             y_plus
@@ -667,7 +667,7 @@ where
     }
 
     /// Project y on set Y
-    fn project_on_set_y(&mut self) {
+    fn project_on_set_y(&mut self) -> FunctionCallResult {
         let problem = &self.alm_problem;
         if let Some(y_set) = &problem.alm_set_y {
             // NOTE: as_mut() converts from &mut Option<T> to Option<&mut T>
@@ -676,9 +676,10 @@ where
             // *  which can be treated as  Option<&mut [f64]>
             // * y_vec is                  &mut [f64]
             if let Some(xi_vec) = self.alm_cache.xi.as_mut() {
-                y_set.project(&mut xi_vec[1..]);
+                y_set.project(&mut xi_vec[1..])?;
             }
         }
+        Ok(())
     }
 
     /// Solve inner problem
@@ -740,7 +741,7 @@ where
         inner_solver.solve(u)
     }
 
-    fn is_exit_criterion_satisfied(&self) -> bool {
+    fn is_exit_criterion_satisfied(&self) -> Result<bool, SolverError> {
         let cache = &self.alm_cache;
         let problem = &self.alm_problem;
         // Criterion 1: ||Delta y|| <= c * delta
@@ -764,8 +765,14 @@ where
         //              This should never happen because we set the AKKT tolerance
         //              in the constructor and can never become `None` again
         let criterion_3 =
-            cache.panoc_cache.akkt_tolerance.unwrap() <= self.epsilon_tolerance + SMALL_EPSILON;
-        criterion_1 && criterion_2 && criterion_3
+            cache
+                .panoc_cache
+                .akkt_tolerance
+                .ok_or(SolverError::InvalidProblemState(
+                    "missing inner AKKT tolerance while checking the exit criterion",
+                ))?
+                <= self.epsilon_tolerance + SMALL_EPSILON;
+        Ok(criterion_1 && criterion_2 && criterion_3)
     }
 
     /// Whether the penalty parameter should not be updated
@@ -802,13 +809,21 @@ where
         }
     }
 
-    fn update_inner_akkt_tolerance(&mut self) {
+    fn update_inner_akkt_tolerance(&mut self) -> FunctionCallResult {
         let cache = &mut self.alm_cache;
         // epsilon_{nu+1} := max(epsilon, beta*epsilon_nu)
+        let akkt_tolerance =
+            cache
+                .panoc_cache
+                .akkt_tolerance
+                .ok_or(SolverError::InvalidProblemState(
+                    "missing inner AKKT tolerance while updating it",
+                ))?;
         cache.panoc_cache.set_akkt_tolerance(f64::max(
-            cache.panoc_cache.akkt_tolerance.unwrap() * self.epsilon_update_factor,
+            akkt_tolerance * self.epsilon_update_factor,
             self.epsilon_tolerance,
         ));
+        Ok(())
     }
 
     fn final_cache_update(&mut self) {
@@ -843,7 +858,7 @@ where
         let mut inner_exit_status: ExitStatus = ExitStatus::Converged;
 
         // Project y on Y
-        self.project_on_set_y();
+        self.project_on_set_y()?;
 
         // If the inner problem fails miserably, the failure should be propagated
         // upstream (using `?`). If the inner problem has not converged, that is fine,
@@ -867,7 +882,7 @@ where
         self.compute_alm_infeasibility()?; // ALM: ||y_plus - y||
 
         // Check exit criterion
-        if self.is_exit_criterion_satisfied() {
+        if self.is_exit_criterion_satisfied()? {
             // Do not continue the outer iteration
             // An (epsilon, delta)-AKKT point has been found
             return Ok(InnerProblemStatus::new(false, inner_exit_status));
@@ -876,7 +891,7 @@ where
         }
 
         // Update inner problem tolerance
-        self.update_inner_akkt_tolerance();
+        self.update_inner_akkt_tolerance()?;
 
         // conclusive step: updated iteration count, resets PANOC cache,
         // sets f2_norm = f2_norm_plus etc
@@ -979,12 +994,11 @@ where
             .with_penalty(c)
             .with_cost(cost);
         if self.alm_problem.n1 > 0 {
-            let status = status.with_lagrange_multipliers(
-                self.alm_cache
-                    .y_plus
-                    .as_ref()
-                    .expect("Although n1 > 0, there is no vector y (Lagrange multipliers)"),
-            );
+            let status = status.with_lagrange_multipliers(self.alm_cache.y_plus.as_ref().ok_or(
+                SolverError::InvalidProblemState(
+                    "missing Lagrange multipliers at the ALM solution",
+                ),
+            )?);
             Ok(status)
         } else {
             Ok(status)
@@ -1129,7 +1143,7 @@ mod tests {
             .with_initial_penalty(25.0)
             .with_initial_lagrange_multipliers(&[2., 3., 4., 10.]);
 
-        alm_optimizer.project_on_set_y();
+        alm_optimizer.project_on_set_y().unwrap();
         if let Some(xi_after_proj) = &alm_optimizer.alm_cache.xi {
             println!("xi = {:#?}", xi_after_proj);
             let y_projected_correct = [
@@ -1282,7 +1296,7 @@ mod tests {
             .with_initial_inner_tolerance(1e-1)
             .with_inner_tolerance_update_factor(0.2);
 
-        alm_optimizer.update_inner_akkt_tolerance();
+        alm_optimizer.update_inner_akkt_tolerance().unwrap();
 
         unit_test_utils::assert_nearly_equal(
             0.1,
@@ -1305,7 +1319,7 @@ mod tests {
         );
 
         for _i in 1..=5 {
-            alm_optimizer.update_inner_akkt_tolerance();
+            alm_optimizer.update_inner_akkt_tolerance().unwrap();
         }
         unit_test_utils::assert_nearly_equal(
             2e-5,
@@ -1411,20 +1425,20 @@ mod tests {
 
         // should not exit yet...
         assert!(
-            !alm_optimizer.is_exit_criterion_satisfied(),
+            !alm_optimizer.is_exit_criterion_satisfied().unwrap(),
             "exists right away"
         );
 
         let alm_optimizer = alm_optimizer
             .with_initial_inner_tolerance(1e-3)
             .with_epsilon_tolerance(1e-3);
-        assert!(!alm_optimizer.is_exit_criterion_satisfied());
+        assert!(!alm_optimizer.is_exit_criterion_satisfied().unwrap());
 
         alm_optimizer.alm_cache.delta_y_norm_plus = 1e-3;
-        assert!(!alm_optimizer.is_exit_criterion_satisfied());
+        assert!(!alm_optimizer.is_exit_criterion_satisfied().unwrap());
 
         alm_optimizer.alm_cache.f2_norm_plus = 1e-3;
-        assert!(alm_optimizer.is_exit_criterion_satisfied());
+        assert!(alm_optimizer.is_exit_criterion_satisfied().unwrap());
     }
 
     #[test]

src/constraints/affine_space.rs

@@ -1,6 +1,5 @@
 use super::Constraint;
-use crate::matrix_operations;
-use crate::CholeskyFactorizer;
+use crate::{matrix_operations, CholeskyFactorizer, FunctionCallResult, SolverError};
 
 use ndarray::{ArrayView1, ArrayView2};
 
@@ -82,26 +81,30 @@ impl Constraint for AffineSpace {
     /// ```
     ///
     /// The result is stored in `x` and it can be verified that $Ax = b$.
-    fn project(&self, x: &mut [f64]) {
+    fn project(&self, x: &mut [f64]) -> FunctionCallResult {
         let n = self.n_cols;
         assert!(x.len() == n, "x has wrong dimension");
 
         // Step 1: Compute e = Ax - b
-        let a = ArrayView2::from_shape((self.n_rows, self.n_cols), &self.a_mat)
-            .expect("invalid A shape");
+        let a = ArrayView2::from_shape((self.n_rows, self.n_cols), &self.a_mat).map_err(|_| {
+            SolverError::InvalidProblemState("failed to construct the affine-space matrix view")
+        })?;
         let x_view = ArrayView1::from(&x[..]);
         let b = ArrayView1::from(&self.b_vec[..]);
         let e = a.dot(&x_view) - b;
-        let e_slice: &[f64] = e.as_slice().unwrap();
+        let e_slice: &[f64] = e.as_slice().ok_or(SolverError::InvalidProblemState(
+            "affine-space residual vector is not stored contiguously",
+        ))?;
 
         // Step 2: Solve AA' z = e and compute z
-        let z = self.factorizer.solve(e_slice).unwrap();
+        let z = self.factorizer.solve(e_slice)?;
 
         // Step 3: Compute x = x - A'z
         let at_z = a.t().dot(&ArrayView1::from(&z[..]));
         for (xi, corr) in x.iter_mut().zip(at_z.iter()) {
             *xi -= *corr;
         }
+        Ok(())
     }
 
     /// Affine sets are convex.

src/constraints/ball1.rs

@@ -1,5 +1,6 @@
 use super::Constraint;
 use super::Simplex;
+use crate::FunctionCallResult;
 
 #[derive(Copy, Clone)]
 /// A norm-1 ball, that is, a set given by $B_1^r = \\{x \in \mathbb{R}^n {}:{} \Vert{}x{}\Vert_1 \leq r\\}$
@@ -23,24 +24,25 @@ impl<'a> Ball1<'a> {
         }
     }
 
-    fn project_on_ball1_centered_at_origin(&self, x: &mut [f64]) {
+    fn project_on_ball1_centered_at_origin(&self, x: &mut [f64]) -> FunctionCallResult {
         if crate::matrix_operations::norm1(x) > self.radius {
             // u = |x| (copied)
             let mut u = vec![0.0; x.len()];
             u.iter_mut()
                 .zip(x.iter())
                 .for_each(|(ui, &xi)| *ui = f64::abs(xi));
             // u = P_simplex(u)
-            self.simplex.project(&mut u);
+            self.simplex.project(&mut u)?;
             x.iter_mut()
                 .zip(u.iter())
                 .for_each(|(xi, &ui)| *xi = f64::signum(*xi) * ui);
         }
+        Ok(())
     }
 }
 
 impl<'a> Constraint for Ball1<'a> {
-    fn project(&self, x: &mut [f64]) {
+    fn project(&self, x: &mut [f64]) -> FunctionCallResult {
         if let Some(center) = &self.center {
             assert_eq!(
                 x.len(),
@@ -50,13 +52,14 @@ impl<'a> Constraint for Ball1<'a> {
             x.iter_mut()
                 .zip(center.iter())
                 .for_each(|(xi, &ci)| *xi -= ci);
-            self.project_on_ball1_centered_at_origin(x);
+            self.project_on_ball1_centered_at_origin(x)?;
             x.iter_mut()
                 .zip(center.iter())
                 .for_each(|(xi, &ci)| *xi += ci);
         } else {
-            self.project_on_ball1_centered_at_origin(x);
+            self.project_on_ball1_centered_at_origin(x)?;
         }
+        Ok(())
     }
 
     fn is_convex(&self) -> bool {

src/constraints/ball2.rs

@@ -1,4 +1,5 @@
 use super::Constraint;
+use crate::FunctionCallResult;
 
 #[derive(Copy, Clone)]
 /// A Euclidean ball, that is, a set given by $B_2^r = \\{x \in \mathbb{R}^n {}:{} \Vert{}x{}\Vert \leq r\\}$
@@ -19,7 +20,7 @@ impl<'a> Ball2<'a> {
 }
 
 impl<'a> Constraint for Ball2<'a> {
-    fn project(&self, x: &mut [f64]) {
+    fn project(&self, x: &mut [f64]) -> FunctionCallResult {
         if let Some(center) = &self.center {
             assert_eq!(
                 x.len(),
@@ -46,6 +47,7 @@ impl<'a> Constraint for Ball2<'a> {
                 x.iter_mut().for_each(|x_| *x_ /= norm_over_radius);
             }
         }
+        Ok(())
     }
 
     fn is_convex(&self) -> bool {