Made Python and C implementations match again. Fixed double summary p…

…rinting in some cases
buduo15 · May 5, 2017 · a67f02a · a67f02a
1 parent bcfe8bf
commit a67f02a
Show file tree

Hide file tree

Showing 11 changed files with 80 additions and 47 deletions.
diff --git a/include/types.h b/include/types.h
@@ -243,6 +243,10 @@ typedef struct {
         /// flag indicating whether the solve function has been run before
         c_int first_run;
         #endif
+
+        #ifdef PRINTING
+        c_int summary_printed;  ///< Has last summary been printed? (true/false)
+        #endif
 } OSQPWorkspace;
 
 

diff --git a/include/util.h b/include/util.h
@@ -65,15 +65,15 @@ void print_header(void);
 
 /**
  * Print iteration summary
- * @param info current information
+ * @param work current workspace
  */
-void print_summary(OSQPInfo * info);
+void print_summary(OSQPWorkspace * work);
 
 /**
  * Print information after polish
- * @param info information structure
+ * @param work current workspace
  */
-void print_polish(OSQPInfo * info);
+void print_polish(OSQPWorkspace * work);
 
 /**
  * Print footer when algorithm terminates

diff --git a/interfaces/python/modulepurepy/_osqp.py b/interfaces/python/modulepurepy/_osqp.py
@@ -19,7 +19,7 @@
 OSQP_UNSOLVED = -10
 
 # Printing interval
-PRINT_INTERVAL = 1
+PRINT_INTERVAL = 100
 
 # OSQP Infinity
 OSQP_INFTY = 1e+20
@@ -246,6 +246,8 @@ def __init__(self, work):
                             spspa.eye(work.data.n), work.data.A.T]),
               spspa.hstack([work.data.A,
                            -1./work.settings.rho * spspa.eye(work.data.m)])])
+        # print("KKT")
+        # print(KKT.todense())
 
         # Initialize private structure
         self.kkt_factor = spla.splu(KKT.tocsc())
@@ -360,7 +362,7 @@ def print_setup_header(self, data, settings):
             self.version)
         print("              Pure Python Implementation")
         print("     (c) Bartolomeo Stellato, Goran Banjac")
-        print("   University of Oxford  -  Stanford University 2016")
+        print("   University of Oxford  -  Stanford University 2017")
         print("-------------------------------------------------------")
 
         print("Problem:  variables n = %d, constraints m = %d" % \
@@ -482,7 +484,9 @@ def compute_pri_res(self, polish):
             return la.norm(np.maximum(self.work.data.l - self.work.pol.z, 0) +
                            np.maximum(self.work.pol.z - self.work.data.u, 0), np.inf)
         else:
-            return la.norm(self.work.data.A.dot(self.work.x) - self.work.z, np.inf)
+
+            return la.norm(self.work.data.A.dot(self.work.x) - self.work.z,
+                           np.inf)
 
     def compute_dua_res(self, polish):
         """
@@ -835,6 +839,7 @@ def solve(self):
 
             # Second step: update x and z
             self.update_x()
+
             self.update_z()
 
             # Third step: update y

diff --git a/interfaces/python/tests/general/small_prob.py b/interfaces/python/tests/general/small_prob.py
@@ -1,26 +1,26 @@
 import osqp
-# import osqppurepy as osqp
+import osqppurepy as osqppurepy
 import scipy.sparse as sparse
 import scipy as sp
 import numpy as np
 import mathprogbasepy as mpbpy
 sp.random.seed(2)
 
-n = 10
-m = 50
-# A = sparse.random(m, n, density=0.9, format='csc')
-A = sparse.eye(n)
+n = 2
+m = 3
+A = sparse.random(m, n, density=0.6, format='csc')
+# A = sparse.eye(m)
 
 # l = -sp.rand(m) * 2.
 # u = sp.rand(m) * 2.
 
-l = -sp.rand(n) * 2. + 1e6
-u = sp.rand(n) * 2. + 1e6
+l = -sp.rand(m)
+u = sp.rand(m)
 
 
 # norm_l = np.linalg.norm(l)
 # norm_u = np.linalg.norm(u)
-Escal = 1e3
+Escal = 1
 
 A /= Escal
 l /= Escal
@@ -29,8 +29,8 @@
 
 
 
-P = sparse.random(n, n, density=0.9, format='csc')
-P = P.dot(P.T)
+P = sparse.random(n, n, density=0.6, format='csc')
+P = P.dot(P.T).tocsc()
 q = sp.randn(n)
 
 
@@ -39,33 +39,40 @@
 # norm_q = 1
 
 
-q = q/norm_q
+# q = q/norm_q
 
 # print("old P ")
 # print(P.todense())
-P = P/norm_q
+# P = P/norm_q
 # print("new P ")
 # print(P.todense())
 
-osqp_opts = {'rho': 0.001,
-             'auto_rho': True,
+osqp_opts = {'rho': 0.5,
+             'sigma': 0.001,
+             'auto_rho': False,
              'polish': False,
              'eps_abs': 1e-03,
              'eps_rel': 1e-03,
-            #  'early_terminate_interval': 1,
-            #  'max_iter': 10,
-             'scaling': True}
+             'early_terminate_interval': 1,
+             'max_iter': 174,
+             'scaling': False}
 
 
-qp = mpbpy.QuadprogProblem(P, q, A, l, u)
-res_gurobi = qp.solve(solver=mpbpy.GUROBI, verbose=True)
+
+
+
+# qp = mpbpy.QuadprogProblem(P, q, A, l, u)
+# res_gurobi = qp.solve(solver=mpbpy.GUROBI, verbose=True)
 # res_purepy = qp.solve(solver=OSQP_PUREPY, **osqp_opts)
 # res_osqp = qp.solve(solver=mpbpy.OSQP, **osqp_opts)
 
 model = osqp.OSQP()
 model.setup(P=P, q=q, A=A, l=l, u=u, **osqp_opts)
+res = model.solve()
 
 
+model = osqppurepy.OSQP()
+model.setup(P=P, q=q, A=A, l=l, u=u, **osqp_opts)
 res = model.solve()
 
 

diff --git a/lin_sys/direct/suitesparse/private.c b/lin_sys/direct/suitesparse/private.c
@@ -262,6 +262,7 @@ c_int solve_lin_sys(const OSQPSettings *settings, Priv *p, c_float *b) {
     /* returns solution to linear system */
     /* Ax = b with solution stored in b */
     LDLSolve(b, b, p->L, p->Dinv, p->P, p->bp);
+
     return 0;
 }
 

diff --git a/src/auxil.c b/src/auxil.c
@@ -157,6 +157,7 @@ c_float compute_pri_res(OSQPWorkspace * work, c_int polish){
         return prim_resid;
     } else {
     #endif
+
         // Called from ADMM algorithm (store temporary vector in z_prev)
         mat_vec(work->data->A, work->x, work->z_prev, 0);
         return vec_norm_inf_diff(work->z_prev, work->z, work->data->m);
@@ -379,6 +380,11 @@ void update_info(OSQPWorkspace *work, c_int iter, c_int compute_objective, c_int
     #ifndef EMBEDDED
     }
     #endif
+
+
+    #ifdef PRINTING
+    work->summary_printed = 0;  // The just updated info have not been printed
+    #endif
 }
 
 

diff --git a/src/lin_alg.c b/src/lin_alg.c
@@ -20,7 +20,7 @@ c_float vec_norm_inf(const c_float *v, c_int l) {
     c_float max = 0.0;
     for (i = 0; i < l; i++) {
         abs_v_i = c_absval(v[i]);
-        if ( abs_v_i > max) max = abs_v_i;
+        if ( abs_v_i > max ) max = abs_v_i;
     }
     return max;
 }
@@ -30,7 +30,7 @@ c_float vec_norm_inf_diff(const c_float *a, const c_float *b, c_int l){
     c_int i;
     for (i = 0; i < l; i++) {
         tmp = c_absval(a[i] - b[i]);
-        if ( tmp > nmDiff) nmDiff = tmp;
+        if ( tmp > nmDiff ) nmDiff = tmp;
     }
     return nmDiff;
 }

diff --git a/src/osqp.c b/src/osqp.c
@@ -156,6 +156,7 @@ OSQPWorkspace * osqp_setup(const OSQPData * data, OSQPSettings *settings){
     #ifdef PRINTING
     if (work->settings->verbose)
         print_setup_header(work->data, work->settings);
+        work->summary_printed = 0; // Initialize last summary  to not printed
     #endif
 
     return work;
@@ -256,7 +257,7 @@ c_int osqp_solve(OSQPWorkspace * work){
 
             if (can_print){
                 // Print summary
-                print_summary(work->info);
+                print_summary(work);
             }
 
             // Check algorithm termination
@@ -296,7 +297,7 @@ c_int osqp_solve(OSQPWorkspace * work){
         /* Print summary */
         #ifdef PRINTING
         if (work->settings->verbose)
-            print_summary(work->info);
+            print_summary(work);
         #endif
 
         /* Check whether a termination criterion is triggered */
@@ -311,10 +312,9 @@ c_int osqp_solve(OSQPWorkspace * work){
 
     /* Print summary for last iteration */
     #ifdef PRINTING
-    if (work->settings->verbose
-        && iter % PRINT_INTERVAL != 0 && iter != 1
-        && iter != work->settings->max_iter + 1)
-        print_summary(work->info);
+    if (work->settings->verbose && !work->summary_printed){
+        print_summary(work);
+    }
     #endif
 
     /* if max iterations reached, change status accordingly */

diff --git a/src/polish.c b/src/polish.c
@@ -270,7 +270,7 @@ c_int polish(OSQPWorkspace *work) {
             // Print summary
             #ifdef PRINTING
             if (work->settings->verbose)
-                print_polish(work->info);
+                print_polish(work);
             #endif
 
     } else { // Polishing failed

diff --git a/src/util.c b/src/util.c
@@ -103,21 +103,27 @@ void print_setup_header(const OSQPData *data, const OSQPSettings *settings) {
 }
 
 
-void print_summary(OSQPInfo * info){
+void print_summary(OSQPWorkspace * work){
+    OSQPInfo * info;
+    info = work->info;
+
     c_print("%*i ", (int)strlen(HEADER[0]), (int)info->iter);
     c_print("%*.4e ", (int)HSPACE, info->obj_val);
-    // c_print("%*.4e ", (int)HSPACE, info->pri_res);
-    // c_print("%*.4e ", (int)HSPACE, info->dua_res);
-    c_print("%*.4e ", (int)HSPACE, sqrt(info->pri_res));
-    c_print("%*.4e ", (int)HSPACE, sqrt(info->dua_res));
+    c_print("%*.4e ", (int)HSPACE, info->pri_res);
+    c_print("%*.4e ", (int)HSPACE, info->dua_res);
     #ifdef PROFILING
     c_print("%*.2fs", 9, info->setup_time + info->solve_time);
     #endif
     c_print("\n");
+
+    work->summary_printed = 1; // Summary has been printed
 }
 
 
-void print_polish(OSQPInfo * info) {
+void print_polish(OSQPWorkspace * work) {
+    OSQPInfo * info;
+    info = work->info;
+
     c_print("%*s ", (int)strlen(HEADER[0]), "PLSH");
     c_print("%*.4e ", (int)HSPACE, info->obj_val);
     c_print("%*.4e ", (int)HSPACE, info->pri_res);
@@ -419,11 +425,11 @@ void print_dns_matrix(c_float * M, c_int m, c_int n, const char *name)
                 for(j=0; j<n; j++) { // Cycle over columns
                         if (j < n - 1)
                                 // c_print("% 14.12e,  ", M[j*m+i]);
-                                c_print("% .5f,  ", M[j*m+i]);
+                                c_print("% .8f,  ", M[j*m+i]);
 
                         else
                                 // c_print("% 14.12e;  ", M[j*m+i]);
-                                c_print("% .5f;  ", M[j*m+i]);
+                                c_print("% .8f;  ", M[j*m+i]);
                 }
                 if (i < m - 1) {
                         c_print("\n\t");

diff --git a/tests/solve_linsys/generate_problem.py b/tests/solve_linsys/generate_problem.py
@@ -7,10 +7,14 @@
 np.random.seed(2)
 
 # Simple case
-test_solve_KKT_n = 2
-test_solve_KKT_m = 2
-test_solve_KKT_A = spa.csc_matrix(np.array([[-3., 4.], [2., -3.]]))
-test_solve_KKT_P = spa.csc_matrix(np.array([[0.2, 0.0], [0.0, 0.6]]))
+test_solve_KKT_n = 10
+test_solve_KKT_m = 20
+
+test_solve_KKT_P = spa.random(test_solve_KKT_n, test_solve_KKT_n,
+                                 density=0.4, format='csc')
+test_solve_KKT_P = test_solve_KKT_P.dot(test_solve_KKT_P.T).tocsc()
+test_solve_KKT_A = spa.random(test_solve_KKT_m, test_solve_KKT_n,
+                                 density=0.4, format='csc')
 test_solve_KKT_Pu = spa.triu(test_solve_KKT_P).tocsc()
 
 test_solve_KKT_rho = 4.0