SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Mosek.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2012 Fernando José Iglesias García
8  * Copyright (C) 2012 Fernando José Iglesias García
9  */
10 
11 #ifdef USE_MOSEK
12 
13 //#define DEBUG_MOSEK
14 //#define DEBUG_SOLUTION
15 
18 #include <shogun/lib/SGVector.h>
19 
20 using namespace shogun;
21 
22 CMosek::CMosek()
23 : CSGObject()
24 {
25 }
26 
27 CMosek::CMosek(int32_t num_con, int32_t num_var)
28 : CSGObject()
29 {
30  // Create MOSEK environment
31 #if (MSK_VERSION_MAJOR == 6)
32  m_rescode = MSK_makeenv(&m_env, NULL, NULL, NULL, NULL);
33 #elif (MSK_VERSION_MAJOR == 7)
34  m_rescode = MSK_makeenv(&m_env, NULL);
35 #else
36  #error "Unsupported Mosek version"
37 #endif
38 
39 #ifdef DEBUG_MOSEK
40  // Direct the environment's log stream to SG_PRINT
41  if ( m_rescode == MSK_RES_OK )
42  {
43  m_rescode = MSK_linkfunctoenvstream(m_env, MSK_STREAM_LOG,
44  NULL, CMosek::print);
45  }
46 #endif
47 
48  // Initialize the environment
49  if ( m_rescode == MSK_RES_OK )
50  {
51  m_rescode = MSK_initenv(m_env);
52  }
53 
54  // Create an optimization task
55  if ( m_rescode == MSK_RES_OK )
56  {
57  m_rescode = MSK_maketask(m_env, num_con, num_var, &m_task);
58  }
59 
60 #ifdef DEBUG_MOSEK
61  // Direct the task's log stream to SG_PRINT
62  if ( m_rescode == MSK_RES_OK )
63  {
64  m_rescode = MSK_linkfunctotaskstream(m_task, MSK_STREAM_LOG,
65  NULL, CMosek::print);
66  }
67 #endif
68 }
69 
70 CMosek::~CMosek()
71 {
72  delete_problem();
73 }
74 
75 void MSKAPI CMosek::print(void* handle, char str[])
76 {
77  SG_SPRINT("%s", str)
78 }
79 
80 MSKrescodee CMosek::init_sosvm(int32_t M, int32_t N,
81  int32_t num_aux, int32_t num_aux_con,
85 {
86  // Give an estimate of the size of the input data to increase the
87  // speed of inputting
88  int32_t num_var = M+N+num_aux;
89  int32_t num_con = N*N+num_aux_con;
90  // NOTE: However, to input this step is completely optional and MOSEK
91  // will automatically allocate more resources if necessary
92  m_rescode = MSK_putmaxnumvar(m_task, num_var);
93  // Neither the number of constraints nor the number of non-zero elements
94  // is known a priori, rough estimates are given here
95  m_rescode = MSK_putmaxnumcon(m_task, num_con);
96  // A = [-dPsi(y) | -I_N ] with M+N columns => max. M+1 nnz per row
97  m_rescode = MSK_putmaxnumanz(m_task, (M+1)*N*N);
98 
99  // Append optimization variables initialized to zero
100 #if (MSK_VERSION_MAJOR == 6)
101  m_rescode = MSK_append(m_task, MSK_ACC_VAR, num_var);
102 #elif (MSK_VERSION_MAJOR == 7)
103  m_rescode = MSK_appendvars(m_task, num_var);
104 #else
105  #error "Unsupported Mosek version"
106 #endif
107  // Append empty constraints initialized with no bounds
108 #if (MSK_VERSION_MAJOR == 6)
109  m_rescode = MSK_append(m_task, MSK_ACC_CON, num_con);
110 #elif (MSK_VERSION_MAJOR == 7)
111  m_rescode = MSK_appendcons(m_task, num_con);
112 #else
113  #error "Unsupported Mosek version"
114 #endif
115  // Set the constant term in the objective equal to zero
116  m_rescode = MSK_putcfix(m_task, 0.0);
117 
118  for ( int32_t j = 0 ; j < num_var && m_rescode == MSK_RES_OK ; ++j )
119  {
120  // Set the linear term c_j in the objective
121  if ( j < M+num_aux )
122  m_rescode = MSK_putcj(m_task, j, 0.0);
123  else
124  m_rescode = MSK_putcj(m_task, j, 1.0);
125 
126  // Set the bounds on x_j: blx[j] <= x_j <= bux[j]
127  // TODO set bounds lb and ub given by init_opt for w
128  if ( j < M )
129  {
130  m_rescode = MSK_putbound(m_task, MSK_ACC_VAR, j,
131  MSK_BK_FR, -MSK_INFINITY, +MSK_INFINITY);
132  }
133 
134  // The slack and the auxiliary variables are required to be positive
135  if ( j >= M )
136  {
137  m_rescode = MSK_putbound(m_task, MSK_ACC_VAR, j,
138  MSK_BK_LO, 0.0, +MSK_INFINITY);
139  }
140  }
141 
142  // Input the matrix Q^0 for the objective
143  //
144  // NOTE: In MOSEK we minimize x'*Q^0*x. C != Q0 but Q0 is
145  // just an extended version of C with zeros that make no
146  // difference in MOSEK's sparse representation
147  m_rescode = wrapper_putqobj(C);
148 
149  // Input the matrix A and the vector b for the contraints A*x <= b
150  m_rescode = wrapper_putaveclist(m_task, A);
151  m_rescode = wrapper_putboundlist(m_task, b);
152 
153  REQUIRE(m_rescode == MSK_RES_OK, "MOSEK Error in CMosek::init_sosvm(). "
154  "Enable DEBUG_MOSEK for details.\n");
155 
156  return m_rescode;
157 }
158 
159 MSKrescodee CMosek::add_constraint_sosvm(
161  index_t con_idx,
162  index_t train_idx,
163  int32_t num_aux,
164  float64_t bi)
165 {
166  // Count the number of non-zero element in dPsi
167  int32_t nnz = CMath::get_num_nonzero(dPsi.vector, dPsi.vlen);
168  // Indices of the non-zero elements in the row of A to add
169  SGVector< index_t > asub(nnz+1); // +1 because of the -1 element
170  // Values of the non-zero elements
171  SGVector< float64_t > aval(nnz+1);
172  // Next element to add in asub and aval
173  index_t idx = 0;
174 
175  for ( int32_t i = 0 ; i < dPsi.vlen ; ++i )
176  {
177  if ( dPsi[i] != 0 )
178  {
179  asub[idx] = i;
180  aval[idx] = dPsi[i];
181  ++idx;
182  }
183  }
184 
185  ASSERT(idx == nnz)
186 
187  asub[idx] = dPsi.vlen + num_aux + train_idx;
188  aval[idx] = -1;
189 
190 #if (MSK_VERSION_MAJOR == 6)
191  m_rescode = MSK_putavec(m_task, MSK_ACC_CON, con_idx, nnz+1,
192  asub.vector, aval.vector);
193 #elif (MSK_VERSION_MAJOR == 7)
194  m_rescode = MSK_putarow(m_task, con_idx, nnz+1, asub.vector, aval.vector);
195 #else
196  #error "Unsupported Mosek version"
197 #endif
198 
199  if ( m_rescode == MSK_RES_OK )
200  {
201  m_rescode = MSK_putbound(m_task, MSK_ACC_CON, con_idx,
202  MSK_BK_UP, -MSK_INFINITY, bi);
203  }
204 
205  return m_rescode;
206 }
207 
208 MSKrescodee CMosek::wrapper_putaveclist(
209  MSKtask_t & task,
211 {
212  // Indices to the rows of A to replace, all the rows
214  for ( index_t i = 0 ; i < A.num_rows ; ++i )
215  sub[i] = i;
216 
217  // Non-zero elements of A
218  int32_t nnza = CMath::get_num_nonzero(A.matrix, A.num_rows*A.num_cols);
219  SGVector< float64_t > aval(nnza);
220  // For each of the rows, indices to non-zero elements
221  SGVector< index_t > asub(nnza);
222  // For each row, pointer to the first non-zero element
223  // in aval
225  // Next position to write in aval and asub
226  index_t idx = 0;
227  // Switch if the first non-zero element in each row
228  // has been found
229  bool first_nnz_found = false;
230 
231  for ( index_t i = 0 ; i < A.num_rows ; ++i )
232  {
233  first_nnz_found = false;
234 
235  for ( index_t j = 0 ; j < A.num_cols ; ++j )
236  {
237  if ( A(i,j) )
238  {
239  aval[idx] = A(i,j);
240  asub[idx] = j;
241 
242  if ( !first_nnz_found )
243  {
244  ptrb[i] = idx;
245  first_nnz_found = true;
246  }
247 
248  ++idx;
249  }
250  }
251 
252  // Handle rows whose elements are all zero
253  // TODO does it make sense that a row in A has all its elements
254  // equal to zero?
255  if ( !first_nnz_found )
256  ptrb[i] = ( i ? ptrb[i-1] : 0 );
257  }
258 
259  // For each row, pointer to the last+1 non-zero element
260  // in aval
262  for ( index_t i = 0 ; i < A.num_rows-1 ; ++i )
263  ptre[i] = ptrb[i+1];
264 
265  if ( A.num_rows > 0 )
266  ptre[A.num_rows-1] = nnza;
267 
268  MSKrescodee ret;
269 #if (MSK_VERSION_MAJOR == 6)
270  ret = MSK_putaveclist(task, MSK_ACC_CON, A.num_rows, sub.vector,
271  ptrb.vector, ptre.vector,
272  asub.vector, aval.vector);
273 #elif (MSK_VERSION_MAJOR == 7)
274  ret = MSK_putarowlist(task, A.num_rows, sub.vector, ptrb.vector, ptre.vector,
275  asub.vector, aval.vector);
276 #else
277  #error "Unsupported Mosek version"
278 #endif
279 
280  REQUIRE(ret == MSK_RES_OK, "MOSEK Error in CMosek::wrapper_putaveclist(). "
281  "Enable DEBUG_MOSEK for details.\n");
282 
283  return ret;
284 }
285 
286 MSKrescodee CMosek::wrapper_putboundlist(MSKtask_t & task, SGVector< float64_t > b)
287 {
288  // Indices to the bounds that should be replaced, b.vlen bounds starting
289  // from zero
290  SGVector< index_t > sub(b.vlen);
291  for ( index_t i = 0 ; i < b.vlen ; ++i )
292  sub[i] = i;
293 
294  // Type of the bounds and lower bound values
295  MSKboundkeye* bk = SG_MALLOC(MSKboundkeye, b.vlen);
297  for ( index_t i = 0 ; i < b.vlen ; ++i )
298  {
299  bk[i] = MSK_BK_UP;
300  bl[i] = -MSK_INFINITY;
301  }
302 
303  MSKrescodee ret = MSK_putboundlist(task, MSK_ACC_CON, b.vlen, sub.vector,
304  bk, bl.vector, b.vector);
305 
306  SG_FREE(bk);
307 
308  REQUIRE(ret == MSK_RES_OK, "MOSEK Error in CMosek::wrapper_putboundlist(). "
309  "Enable DEBUG_MOSEK for details.\n");
310 
311  return ret;
312 }
313 
314 MSKrescodee CMosek::wrapper_putqobj(SGMatrix< float64_t > Q0) const
315 {
316  // Shorthands for the dimensions of the matrix
317  index_t N = Q0.num_rows;
318  index_t M = Q0.num_cols;
319 
320  // Count the number of non-zero elements in the lower
321  // triangular part of the matrix
322  int32_t nnz = 0;
323  for ( index_t i = 0 ; i < N ; ++i )
324  for ( index_t j = i ; j < M ; ++j )
325  nnz += ( Q0[j + i*M] ? 1 : 0 );
326 
327  // i subscript for non-zero elements of Q0
328  SGVector< index_t > qosubi(nnz);
329  // j subscript for non-zero elements of Q0
330  SGVector< index_t > qosubj(nnz);
331  // Values of non-zero elements of Q0
332  SGVector< float64_t > qoval(nnz);
333  // Next position to write in the vectors
334  index_t idx = 0;
335 
336  for ( index_t i = 0 ; i < N ; ++i )
337  for ( index_t j = i ; j < M ; ++j )
338  {
339  if ( Q0[j + i*M] )
340  {
341  qosubi[idx] = i;
342  qosubj[idx] = j;
343  qoval[idx] = Q0[j + i*M];
344 
345  ++idx;
346  }
347  }
348 
349  return MSK_putqobj(m_task, nnz, qosubi.vector,
350  qosubj.vector, qoval.vector);
351 }
352 
353 MSKrescodee CMosek::optimize(SGVector< float64_t > sol)
354 {
355  m_rescode = MSK_optimize(m_task);
356 
357 #ifdef DEBUG_MOSEK
358  // Print a summary containing information about the solution
359  MSK_solutionsummary(m_task, MSK_STREAM_LOG);
360 #endif
361 
362  // Read the solution
363  if ( m_rescode == MSK_RES_OK )
364  {
365  // Solution status
366  MSKsolstae solsta;
367  // FIXME posible solutions are:
368  // MSK_SOL_ITR: the interior solution
369  // MSK_SOL_BAS: the basic solution
370  // MSK_SOL_ITG: the integer solution
371 #if (MSK_VERSION_MAJOR == 6)
372  MSK_getsolutionstatus(m_task, MSK_SOL_ITR, NULL, &solsta);
373 #elif (MSK_VERSION_MAJOR == 7)
374  MSK_getsolsta(m_task, MSK_SOL_ITR, &solsta);
375 #else
376  #error "Unsupported Mosek Version"
377 #endif
378  switch (solsta)
379  {
380  case MSK_SOL_STA_OPTIMAL:
381  case MSK_SOL_STA_NEAR_OPTIMAL:
382  MSK_getsolutionslice(m_task,
383  // Request the interior solution
384  MSK_SOL_ITR,
385  // of the optimization vector
386  MSK_SOL_ITEM_XX,
387  0,
388  sol.vlen,
389  sol.vector);
390 #ifdef DEBUG_SOLUTION
391  sol.display_vector("Solution");
392 #endif
393  break;
394  case MSK_SOL_STA_DUAL_INFEAS_CER:
395  case MSK_SOL_STA_PRIM_INFEAS_CER:
396  case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER:
397  case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER:
398 #ifdef DEBUG_MOSEK
399  SG_PRINT("Primal or dual infeasibility "
400  "certificate found\n");
401 #endif
402  break;
403  case MSK_SOL_STA_UNKNOWN:
404 #ifdef DEBUG_MOSEK
405  SG_PRINT("Undetermined solution status\n")
406 #endif
407  break;
408  default:
409 #ifdef DEBUG_MOSEK
410  SG_PRINT("Other solution status\n")
411 #endif
412  break; // to avoid compile error when DEBUG_MOSEK
413  // is not defined
414  }
415  }
416 
417  // In case any error occurred, print the appropriate error message
418  if ( m_rescode != MSK_RES_OK )
419  {
420  char symname[MSK_MAX_STR_LEN];
421  char desc[MSK_MAX_STR_LEN];
422 
423  MSK_getcodedesc(m_rescode, symname, desc);
424 
425  SG_PRINT("An error occurred optimizing with MOSEK\n")
426  SG_PRINT("ERROR %s - '%s'\n", symname, desc)
427  }
428 
429  return m_rescode;
430 }
431 
432 void CMosek::delete_problem()
433 {
434  MSK_deletetask(&m_task);
435  MSK_deleteenv(&m_env);
436 }
437 
438 void CMosek::display_problem()
439 {
440  int32_t num_var, num_con;
441  m_rescode = MSK_getnumvar(m_task, &num_var);
442  m_rescode = MSK_getnumcon(m_task, &num_con);
443 
444  SG_PRINT("\nMatrix Q^0:\n")
445  for ( int32_t i = 0 ; i < num_var ; ++i )
446  {
447  for ( int32_t j = 0 ; j < num_var ; ++j )
448  {
449  float64_t qij;
450  m_rescode = MSK_getqobjij(m_task, i, j, &qij);
451  if ( qij != 0.0 )
452  SG_PRINT("(%d,%d)\t%.2f\n", i, j, qij)
453  }
454  }
455  SG_PRINT("\n")
456 
457  SG_PRINT("\nVector c:\n")
458  SGVector< float64_t > c(num_var);
459  m_rescode = MSK_getc(m_task, c.vector);
460  c.display_vector();
461 
462  SG_PRINT("\n\nMatrix A:\n")
463  for ( int32_t i = 0 ; i < num_con ; ++i )
464  {
465  for ( int32_t j = 0 ; j < num_var ; ++j )
466  {
467  float64_t aij;
468  m_rescode = MSK_getaij(m_task, i, j, &aij);
469  if ( aij != 0.0 )
470  SG_PRINT("(%d,%d)\t%.2f\n", i, j, aij)
471  }
472  }
473  SG_PRINT("\n")
474 
475  SG_PRINT("\nConstraint Bounds, vector b:\n")
476  for ( int32_t i = 0 ; i < num_con ; ++i )
477  {
478  MSKboundkeye bk;
479  float64_t bl, bu;
480  m_rescode = MSK_getbound(m_task, MSK_ACC_CON, i, &bk, &bl, &bu);
481 
482  SG_PRINT("%6.2f %6.2f\n", bl, bu)
483  }
484 
485  SG_PRINT("\nVariable Bounds, vectors lb and ub:\n")
486  for ( int32_t i = 0 ; i < num_var ; ++i )
487  {
488  MSKboundkeye bk;
489  float64_t bl, bu;
490  m_rescode = MSK_getbound(m_task, MSK_ACC_VAR, i, &bk, &bl, &bu);
491 
492  SG_PRINT("%6.2f %6.2f\n", bl, bu)
493  }
494  SG_PRINT("\n")
495 }
496 
497 
498 float64_t CMosek::get_primal_objective_value() const
499 {
500  float64_t po = 0.0;
501  MSK_getprimalobj(m_task, MSK_SOL_ITR, &po);
502 
503  return po;
504 }
505 
506 #endif /* USE_MOSEK */

SHOGUN Machine Learning Toolbox - Documentation