merge dev/release and fix conflicts

Author: alphaville

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

@@ -285,7 +285,7 @@ where
                 .zip(y_lagrange_mult.iter())
                 .for_each(|(ti, yi)| *ti += *yi / penalty_scale);
             s.copy_from_slice(&f1_u_plus_y_over_c);
-            set_c.project(&mut s);
+            set_c.project(&mut s)?;
             let dist_sq: T = matrix_operations::norm2_squared_diff(&f1_u_plus_y_over_c, &s);
             let scaling: T = cast::<T>(0.5) * penalty_parameter;
             *cost += scaling * dist_sq;
@@ -403,7 +403,7 @@ where
                 .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
@@ -519,7 +519,9 @@ mod tests {
         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

@@ -674,7 +674,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
@@ -691,7 +691,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>
@@ -700,9 +700,10 @@ where
             // *  which can be treated as  Option<&mut [T]>
             // * y_vec is                  &mut [T]
             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
@@ -764,7 +765,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
@@ -787,9 +788,14 @@ where
         //              This function will panic is there is no akkt_tolerance
         //              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 + T::epsilon();
-        criterion_1 && criterion_2 && criterion_3
+        let criterion_3 = cache
+            .panoc_cache
+            .akkt_tolerance
+            .ok_or(SolverError::InvalidProblemState(
+                "missing inner AKKT tolerance while checking the exit criterion",
+            ))?
+            <= self.epsilon_tolerance + T::epsilon();
+        Ok(criterion_1 && criterion_2 && criterion_3)
     }
 
     /// Whether the penalty parameter should not be updated
@@ -826,17 +832,23 @@ 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 next_tolerance = cache.panoc_cache.akkt_tolerance.unwrap() * self.epsilon_update_factor;
-        cache
-            .panoc_cache
-            .set_akkt_tolerance(if next_tolerance > self.epsilon_tolerance {
-                next_tolerance
-            } else {
-                self.epsilon_tolerance
-            });
+        let akkt_tolerance =
+            cache
+                .panoc_cache
+                .akkt_tolerance
+                .ok_or(SolverError::InvalidProblemState(
+                    "missing inner AKKT tolerance while updating it",
+                ))?;
+        let next_tolerance = akkt_tolerance * self.epsilon_update_factor;
+        cache.panoc_cache.set_akkt_tolerance(if next_tolerance > self.epsilon_tolerance {
+            next_tolerance
+        } else {
+            self.epsilon_tolerance
+        });
+        Ok(())
     }
 
     fn final_cache_update(&mut self) {
@@ -871,7 +883,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,
@@ -895,7 +907,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));
@@ -904,7 +916,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
@@ -1010,12 +1022,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)
@@ -1163,7 +1174,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 = [
@@ -1316,7 +1327,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,
@@ -1339,7 +1350,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,
@@ -1445,20 +1456,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,7 @@
 use super::Constraint;
-use crate::matrix_operations;
-use crate::{CholeskyError, CholeskyFactorizer};
+use crate::{
+    matrix_operations, CholeskyError, CholeskyFactorizer, FunctionCallResult, SolverError,
+};
 
 use ndarray::{ArrayView1, ArrayView2, LinalgScalar};
 use num::Float;
@@ -148,26 +149,30 @@ where
     /// ```
     ///
     /// The result is stored in `x` and it can be verified that $Ax = b$.
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> 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: &[T] = e.as_slice().unwrap();
+        let e_slice: &[T] = 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 = *xi - *corr;
         }
+        Ok(())
     }
 
     /// Affine sets are convex.

src/constraints/ball1.rs

@@ -2,6 +2,7 @@ use super::Constraint;
 use super::Simplex;
 use num::Float;
 use std::iter::Sum;
+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\\}$
@@ -35,7 +36,7 @@ impl<'a, T: Float> Ball1<'a, T> {
         }
     }
 
-    fn project_on_ball1_centered_at_origin(&self, x: &mut [T])
+    fn project_on_ball1_centered_at_origin(&self, x: &mut [T]) -> FunctionCallResult
     where
         T: Sum<T>,
     {
@@ -46,19 +47,20 @@ impl<'a, T: Float> Ball1<'a, T> {
                 .zip(x.iter())
                 .for_each(|(ui, &xi)| *ui = xi.abs());
             // 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 = xi.signum() * ui);
         }
+        Ok(())
     }
 }
 
 impl<'a, T> Constraint<T> for Ball1<'a, T>
 where
     T: Float + Sum<T>,
 {
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         if let Some(center) = &self.center {
             assert_eq!(
                 x.len(),
@@ -68,13 +70,14 @@ where
             x.iter_mut()
                 .zip(center.iter())
                 .for_each(|(xi, &ci)| *xi = *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 = *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 {