mpc bug test (#280)

* mpc bug test

* explicit memset to 0

* not sure what is failing

* fix n

* add some prints

* better print

* test disable res infeas

* disable other tests

* Revert "disable other tests"

This reverts commit b1a41ba.

* cleanup

* add more tests

* rename var

---------

Co-authored-by: Brendan O'Donoghue <bodonoghue@google.com>

include/glbopts.h

@@ -39,7 +39,8 @@ extern "C" {
 #define WRITE_DATA_FILENAME (0)
 #define LOG_CSV_FILENAME (0)
 #define TIME_LIMIT_SECS (0.)
-
+/* Tolerance to check negativity condition for infeasibility */
+#define INFEAS_NEGATIVITY_TOL (1e-9)
 /* redefine printfs as needed */
 #if NO_PRINTING > 0     /* Disable all printing */
 #define scs_printf(...) /* No-op */

src/scs.c

@@ -207,13 +207,13 @@ static void compute_residuals(ScsResiduals *r, scs_int m, scs_int n) {
   r->res_unbdd_a = NAN;
   r->res_unbdd_p = NAN;
   r->res_infeas = NAN;
-  if (r->ctx_tau < 0) {
+  if (r->ctx_tau < -INFEAS_NEGATIVITY_TOL) {
     nm_ax_s = NORM(r->ax_s, m);
     nm_px = NORM(r->px, n);
     r->res_unbdd_a = SAFEDIV_POS(nm_ax_s, -r->ctx_tau);
     r->res_unbdd_p = SAFEDIV_POS(nm_px, -r->ctx_tau);
   }
-  if (r->bty_tau < 0) {
+  if (r->bty_tau < -INFEAS_NEGATIVITY_TOL) {
     nm_aty = NORM(r->aty, n);
     r->res_infeas = SAFEDIV_POS(nm_aty, -r->bty_tau);
   }
@@ -587,22 +587,22 @@ static void print_summary(ScsWork *w, scs_int i, SCS(timer) * solve_timer) {
   scs_printf("\n");
 
 #if VERBOSITY > 0
-  scs_printf("Norm u = %4f, ", SCS(norm_2)(w->u, w->d->n + w->d->m + 1));
-  scs_printf("Norm u_t = %4f, ", SCS(norm_2)(w->u_t, w->d->n + w->d->m + 1));
-  scs_printf("Norm v = %4f, ", SCS(norm_2)(w->v, w->d->n + w->d->m + 1));
-  scs_printf("Norm rsk = %4f, ", SCS(norm_2)(w->rsk, w->d->n + w->d->m + 1));
-  scs_printf("Norm x = %4f, ", SCS(norm_2)(w->xys_orig->x, w->d->n));
-  scs_printf("Norm y = %4f, ", SCS(norm_2)(w->xys_orig->y, w->d->m));
-  scs_printf("Norm s = %4f, ", SCS(norm_2)(w->xys_orig->s, w->d->m));
-  scs_printf("Norm |Ax + s| = %1.2e, ", SCS(norm_2)(r->ax_s, w->d->m));
-  scs_printf("tau = %4f, ", w->u[w->d->n + w->d->m]);
-  scs_printf("kappa = %4f, ", w->rsk[w->d->n + w->d->m]);
-  scs_printf("|u - u_t| = %1.2e, ",
+  scs_printf("Norm u = %1.6e, ", SCS(norm_2)(w->u, w->d->n + w->d->m + 1));
+  scs_printf("Norm u_t = %1.6e, ", SCS(norm_2)(w->u_t, w->d->n + w->d->m + 1));
+  scs_printf("Norm v = %1.6e, ", SCS(norm_2)(w->v, w->d->n + w->d->m + 1));
+  scs_printf("Norm rsk = %1.6e, ", SCS(norm_2)(w->rsk, w->d->n + w->d->m + 1));
+  scs_printf("Norm x = %1.6e, ", SCS(norm_2)(w->xys_orig->x, w->d->n));
+  scs_printf("Norm y = %1.6e, ", SCS(norm_2)(w->xys_orig->y, w->d->m));
+  scs_printf("Norm s = %1.6e, ", SCS(norm_2)(w->xys_orig->s, w->d->m));
+  scs_printf("Norm |Ax + s| = %1.6e, ", SCS(norm_2)(r->ax_s, w->d->m));
+  scs_printf("tau = %1.6e, ", w->u[w->d->n + w->d->m]);
+  scs_printf("kappa = %1.6e, ", w->rsk[w->d->n + w->d->m]);
+  scs_printf("|u - u_t| = %1.6e, ",
              SCS(norm_diff)(w->u, w->u_t, w->d->n + w->d->m + 1));
-  scs_printf("res_infeas = %1.2e, ", r->res_infeas);
-  scs_printf("res_unbdd_a = %1.2e, ", r->res_unbdd_a);
-  scs_printf("res_unbdd_p = %1.2e, ", r->res_unbdd_p);
-  scs_printf("ctx_tau = %1.2e, ", r->ctx_tau);
+  scs_printf("res_infeas = %1.6e, ", r->res_infeas);
+  scs_printf("res_unbdd_a = %1.6e, ", r->res_unbdd_a);
+  scs_printf("res_unbdd_p = %1.6e, ", r->res_unbdd_p);
+  scs_printf("ctx_tau = %1.6e, ", r->ctx_tau);
   scs_printf("bty_tau = %1.2e\n", r->bty_tau);
 #endif
 

test/problems/mpc_bug.h

@@ -0,0 +1,19 @@
+#include "glbopts.h"
+#include "problems/test_prob_from_data_file.h"
+#include "scs.h"
+
+static const char *mpc_bug(void) {
+  const char *fail;
+  scs_float OPT1 = -0.473957794500; /* from scs */
+  scs_float OPT2 = -0.029336830816; /* from scs */
+  scs_float OPT3 = -0.002215217478; /* from scs */
+  fail = _test_prob_from_data("test/problems/mpc_bug1", OPT1);
+  if (fail) {
+    return fail;
+  }
+  fail = _test_prob_from_data("test/problems/mpc_bug2", OPT2);
+  if (fail) {
+    return fail;
+  }
+  return _test_prob_from_data("test/problems/mpc_bug3", OPT3);
+}

test/problems/test_prob_from_data_file.h

@@ -19,6 +19,8 @@ static const char *_test_prob_from_data(const char *file, scs_float OPT) {
   scs_float perr, derr;
   scs_int success;
   const char *fail;
+  scs_float xt_p_x;
+  scs_float *px = SCS_NULL;
 
   read_status = SCS(read_data)(file, &d, &k, &stgs);
 
@@ -28,12 +30,24 @@ static const char *_test_prob_from_data(const char *file, scs_float OPT) {
 
   stgs->eps_abs = 1e-6;
   stgs->eps_rel = 1e-6;
+  /* Force verbosity for the test */
+  stgs->verbose = 1;
 
   sol = (ScsSolution *)scs_calloc(1, sizeof(ScsSolution));
   exitflag = scs(d, k, stgs, sol, &info);
 
-  perr = SCS(dot)(d->c, sol->x, d->n) - OPT;
-  derr = -SCS(dot)(d->b, sol->y, d->m) - OPT;
+  if (d->P) {
+    /* px = Px */
+    px = (scs_float *)scs_calloc(d->n, sizeof(scs_float));
+    memset(px, 0, d->n * sizeof(scs_float));
+    SCS(accum_by_p)(d->P, sol->x, px);
+    xt_p_x = SCS(dot)(px, sol->x, d->n);
+  } else {
+    xt_p_x = 0.;
+  }
+
+  perr = 0.5 * xt_p_x + SCS(dot)(d->c, sol->x, d->n) - OPT;
+  derr = -0.5 * xt_p_x - SCS(dot)(d->b, sol->y, d->m) - OPT;
   scs_printf("primal obj error %4e\n", perr);
   scs_printf("dual obj error %4e\n", derr);
 
@@ -47,7 +61,9 @@ static const char *_test_prob_from_data(const char *file, scs_float OPT) {
   SCS(free_cone)(k);
   SCS(free_sol)(sol);
   scs_free(stgs);
-
+  if (px) {
+    scs_free(px);
+  }
   if (fail) {
     scs_printf("%s: FAILED\n", file);
   }

test/run_tests.c

@@ -43,9 +43,11 @@ _SKIP(random_prob)
 #if NO_READ_WRITE == 0 /* reads / writes */
 #include "problems/hs21_tiny_qp_rw.h"
 #include "problems/max_ent.h"
+#include "problems/mpc_bug.h"
 #else
 _SKIP(hs21_tiny_qp_rw)
 _SKIP(max_ent)
+_SKIP(mpc_bug)
 #endif
 
 static const char *all_tests(void) {
@@ -61,6 +63,7 @@ static const char *all_tests(void) {
   mu_run_test(unbounded_tiny_qp);
   mu_run_test(random_prob);
   mu_run_test(max_ent);
+  mu_run_test(mpc_bug);
   mu_run_test(test_exp_cone);
   return 0;
 }