SHOGUN  4.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules 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  // Set default lower and upper bounds on variables
119  if (lb.vlen == 0)
120  {
121  lb = SGVector< float64_t >(M);
122  SGVector< float64_t >::fill_vector(lb.vector, lb.vlen, -MSK_INFINITY);
123  }
124 
125  if (ub.vlen == 0)
126  {
127  ub = SGVector< float64_t >(M);
128  SGVector< float64_t >::fill_vector(ub.vector, ub.vlen, +MSK_INFINITY);
129  }
130 
131  REQUIRE(lb.vlen >= M, "CMosek::init_sosvm(): lb.vlen < dimension of feature space.\n");
132  REQUIRE(ub.vlen >= M, "CMosek::init_sosvm(): ub.vlen < dimension of feature space.\n");
133 
134  for ( int32_t j = 0 ; j < num_var && m_rescode == MSK_RES_OK ; ++j )
135  {
136  // Set the linear term c_j in the objective
137  if ( j < M+num_aux )
138  m_rescode = MSK_putcj(m_task, j, 0.0);
139  else
140  m_rescode = MSK_putcj(m_task, j, 1.0);
141 
142  // Set the bounds on x_j: lb[j] <= x_j <= ub[j]
143  if ( j < M )
144  {
145  if (lb[j] > -MSK_INFINITY && ub[j] < +MSK_INFINITY && lb[j] != ub[j])
146  {
147  // ranged
148  m_rescode = MSK_putbound(m_task, MSK_ACC_VAR, j,
149  MSK_BK_RA, lb[j], ub[j]);
150  }
151  else if (lb[j] > -MSK_INFINITY && ub[j] < +MSK_INFINITY && lb[j] == ub[j])
152  {
153  // fixed
154  m_rescode = MSK_putbound(m_task, MSK_ACC_VAR, j,
155  MSK_BK_FX, lb[j], ub[j]);
156  }
157  else if (lb[j] <= -MSK_INFINITY && ub[j] < +MSK_INFINITY)
158  {
159  // upper bounded
160  m_rescode = MSK_putbound(m_task, MSK_ACC_VAR, j,
161  MSK_BK_UP, -MSK_INFINITY, ub[j]);
162  }
163  else if (lb[j] > -MSK_INFINITY && ub[j] >= +MSK_INFINITY)
164  {
165  // lower bounded
166  m_rescode = MSK_putbound(m_task, MSK_ACC_VAR, j,
167  MSK_BK_LO, lb[j], +MSK_INFINITY);
168  }
169  else
170  {
171  // free
172  m_rescode = MSK_putbound(m_task, MSK_ACC_VAR, j,
173  MSK_BK_FR, -MSK_INFINITY, +MSK_INFINITY);
174  }
175  }
176 
177  // The slack and the auxiliary variables are required to be positive
178  if ( j >= M )
179  {
180  m_rescode = MSK_putbound(m_task, MSK_ACC_VAR, j,
181  MSK_BK_LO, 0.0, +MSK_INFINITY);
182  }
183  }
184 
185  // Input the matrix Q^0 for the objective
186  //
187  // NOTE: In MOSEK we minimize x'*Q^0*x. C != Q0 but Q0 is
188  // just an extended version of C with zeros that make no
189  // difference in MOSEK's sparse representation
190  m_rescode = wrapper_putqobj(C);
191 
192  // Input the matrix A and the vector b for the contraints A*x <= b
193  m_rescode = wrapper_putaveclist(m_task, A);
194  m_rescode = wrapper_putboundlist(m_task, b);
195 
196  REQUIRE(m_rescode == MSK_RES_OK, "MOSEK Error in CMosek::init_sosvm(). "
197  "Enable DEBUG_MOSEK for details.\n");
198 
199  return m_rescode;
200 }
201 
202 MSKrescodee CMosek::add_constraint_sosvm(
204  index_t con_idx,
205  index_t train_idx,
206  int32_t num_aux,
207  float64_t bi)
208 {
209  // Count the number of non-zero element in dPsi
210  int32_t nnz = CMath::get_num_nonzero(dPsi.vector, dPsi.vlen);
211  // Indices of the non-zero elements in the row of A to add
212  SGVector< index_t > asub(nnz+1); // +1 because of the -1 element
213  // Values of the non-zero elements
214  SGVector< float64_t > aval(nnz+1);
215  // Next element to add in asub and aval
216  index_t idx = 0;
217 
218  for ( int32_t i = 0 ; i < dPsi.vlen ; ++i )
219  {
220  if ( dPsi[i] != 0 )
221  {
222  asub[idx] = i;
223  aval[idx] = dPsi[i];
224  ++idx;
225  }
226  }
227 
228  ASSERT(idx == nnz)
229 
230  asub[idx] = dPsi.vlen + num_aux + train_idx;
231  aval[idx] = -1;
232 
233 #if (MSK_VERSION_MAJOR == 6)
234  m_rescode = MSK_putavec(m_task, MSK_ACC_CON, con_idx, nnz+1,
235  asub.vector, aval.vector);
236 #elif (MSK_VERSION_MAJOR == 7)
237  m_rescode = MSK_putarow(m_task, con_idx, nnz+1, asub.vector, aval.vector);
238 #else
239  #error "Unsupported Mosek version"
240 #endif
241 
242  if ( m_rescode == MSK_RES_OK )
243  {
244  m_rescode = MSK_putbound(m_task, MSK_ACC_CON, con_idx,
245  MSK_BK_UP, -MSK_INFINITY, bi);
246  }
247 
248  return m_rescode;
249 }
250 
251 MSKrescodee CMosek::wrapper_putaveclist(
252  MSKtask_t & task,
254 {
255  // Indices to the rows of A to replace, all the rows
257  for ( index_t i = 0 ; i < A.num_rows ; ++i )
258  sub[i] = i;
259 
260  // Non-zero elements of A
261  int32_t nnza = CMath::get_num_nonzero(A.matrix, A.num_rows*A.num_cols);
262  SGVector< float64_t > aval(nnza);
263  // For each of the rows, indices to non-zero elements
264  SGVector< index_t > asub(nnza);
265  // For each row, pointer to the first non-zero element
266  // in aval
268  // Next position to write in aval and asub
269  index_t idx = 0;
270  // Switch if the first non-zero element in each row
271  // has been found
272  bool first_nnz_found = false;
273 
274  for ( index_t i = 0 ; i < A.num_rows ; ++i )
275  {
276  first_nnz_found = false;
277 
278  for ( index_t j = 0 ; j < A.num_cols ; ++j )
279  {
280  if ( A(i,j) )
281  {
282  aval[idx] = A(i,j);
283  asub[idx] = j;
284 
285  if ( !first_nnz_found )
286  {
287  ptrb[i] = idx;
288  first_nnz_found = true;
289  }
290 
291  ++idx;
292  }
293  }
294 
295  // Handle rows whose elements are all zero
296  // TODO does it make sense that a row in A has all its elements
297  // equal to zero?
298  if ( !first_nnz_found )
299  ptrb[i] = ( i ? ptrb[i-1] : 0 );
300  }
301 
302  // For each row, pointer to the last+1 non-zero element
303  // in aval
305  for ( index_t i = 0 ; i < A.num_rows-1 ; ++i )
306  ptre[i] = ptrb[i+1];
307 
308  if ( A.num_rows > 0 )
309  ptre[A.num_rows-1] = nnza;
310 
311  MSKrescodee ret;
312 #if (MSK_VERSION_MAJOR == 6)
313  ret = MSK_putaveclist(task, MSK_ACC_CON, A.num_rows, sub.vector,
314  ptrb.vector, ptre.vector,
315  asub.vector, aval.vector);
316 #elif (MSK_VERSION_MAJOR == 7)
317  ret = MSK_putarowlist(task, A.num_rows, sub.vector, ptrb.vector, ptre.vector,
318  asub.vector, aval.vector);
319 #else
320  #error "Unsupported Mosek version"
321 #endif
322 
323  REQUIRE(ret == MSK_RES_OK, "MOSEK Error in CMosek::wrapper_putaveclist(). "
324  "Enable DEBUG_MOSEK for details.\n");
325 
326  return ret;
327 }
328 
329 MSKrescodee CMosek::wrapper_putboundlist(MSKtask_t & task, SGVector< float64_t > b)
330 {
331  // Indices to the bounds that should be replaced, b.vlen bounds starting
332  // from zero
333  SGVector< index_t > sub(b.vlen);
334  for ( index_t i = 0 ; i < b.vlen ; ++i )
335  sub[i] = i;
336 
337  // Type of the bounds and lower bound values
338  MSKboundkeye* bk = SG_MALLOC(MSKboundkeye, b.vlen);
340  for ( index_t i = 0 ; i < b.vlen ; ++i )
341  {
342  bk[i] = MSK_BK_UP;
343  bl[i] = -MSK_INFINITY;
344  }
345 
346  MSKrescodee ret = MSK_putboundlist(task, MSK_ACC_CON, b.vlen, sub.vector,
347  bk, bl.vector, b.vector);
348 
349  SG_FREE(bk);
350 
351  REQUIRE(ret == MSK_RES_OK, "MOSEK Error in CMosek::wrapper_putboundlist(). "
352  "Enable DEBUG_MOSEK for details.\n");
353 
354  return ret;
355 }
356 
357 MSKrescodee CMosek::wrapper_putqobj(SGMatrix< float64_t > Q0) const
358 {
359  // Shorthands for the dimensions of the matrix
360  index_t N = Q0.num_rows;
361  index_t M = Q0.num_cols;
362 
363  // Count the number of non-zero elements in the lower
364  // triangular part of the matrix
365  int32_t nnz = 0;
366  for ( index_t i = 0 ; i < N ; ++i )
367  for ( index_t j = i ; j < M ; ++j )
368  nnz += ( Q0[j + i*M] ? 1 : 0 );
369 
370  // i subscript for non-zero elements of Q0
371  SGVector< index_t > qosubi(nnz);
372  // j subscript for non-zero elements of Q0
373  SGVector< index_t > qosubj(nnz);
374  // Values of non-zero elements of Q0
375  SGVector< float64_t > qoval(nnz);
376  // Next position to write in the vectors
377  index_t idx = 0;
378 
379  for ( index_t i = 0 ; i < N ; ++i )
380  for ( index_t j = i ; j < M ; ++j )
381  {
382  if ( Q0[j + i*M] )
383  {
384  qosubi[idx] = i;
385  qosubj[idx] = j;
386  qoval[idx] = Q0[j + i*M];
387 
388  ++idx;
389  }
390  }
391 
392  return MSK_putqobj(m_task, nnz, qosubi.vector,
393  qosubj.vector, qoval.vector);
394 }
395 
396 MSKrescodee CMosek::optimize(SGVector< float64_t > sol)
397 {
398  m_rescode = MSK_optimize(m_task);
399 
400 #ifdef DEBUG_MOSEK
401  // Print a summary containing information about the solution
402  MSK_solutionsummary(m_task, MSK_STREAM_LOG);
403 #endif
404 
405  // Read the solution
406  if ( m_rescode == MSK_RES_OK )
407  {
408  // Solution status
409  MSKsolstae solsta;
410  // FIXME posible solutions are:
411  // MSK_SOL_ITR: the interior solution
412  // MSK_SOL_BAS: the basic solution
413  // MSK_SOL_ITG: the integer solution
414 #if (MSK_VERSION_MAJOR == 6)
415  MSK_getsolutionstatus(m_task, MSK_SOL_ITR, NULL, &solsta);
416 #elif (MSK_VERSION_MAJOR == 7)
417  MSK_getsolsta(m_task, MSK_SOL_ITR, &solsta);
418 #else
419  #error "Unsupported Mosek Version"
420 #endif
421  switch (solsta)
422  {
423  case MSK_SOL_STA_OPTIMAL:
424  case MSK_SOL_STA_NEAR_OPTIMAL:
425  MSK_getsolutionslice(m_task,
426  // Request the interior solution
427  MSK_SOL_ITR,
428  // of the optimization vector
429  MSK_SOL_ITEM_XX,
430  0,
431  sol.vlen,
432  sol.vector);
433 #ifdef DEBUG_SOLUTION
434  sol.display_vector("Solution");
435 #endif
436  break;
437  case MSK_SOL_STA_DUAL_INFEAS_CER:
438  case MSK_SOL_STA_PRIM_INFEAS_CER:
439  case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER:
440  case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER:
441 #ifdef DEBUG_MOSEK
442  SG_PRINT("Primal or dual infeasibility "
443  "certificate found\n");
444 #endif
445  break;
446  case MSK_SOL_STA_UNKNOWN:
447 #ifdef DEBUG_MOSEK
448  SG_PRINT("Undetermined solution status\n")
449 #endif
450  break;
451  default:
452 #ifdef DEBUG_MOSEK
453  SG_PRINT("Other solution status\n")
454 #endif
455  break; // to avoid compile error when DEBUG_MOSEK
456  // is not defined
457  }
458  }
459 
460  // In case any error occurred, print the appropriate error message
461  if ( m_rescode != MSK_RES_OK )
462  {
463  char symname[MSK_MAX_STR_LEN];
464  char desc[MSK_MAX_STR_LEN];
465 
466  MSK_getcodedesc(m_rescode, symname, desc);
467 
468  SG_PRINT("An error occurred optimizing with MOSEK\n")
469  SG_PRINT("ERROR %s - '%s'\n", symname, desc)
470  }
471 
472  return m_rescode;
473 }
474 
475 void CMosek::delete_problem()
476 {
477  MSK_deletetask(&m_task);
478  MSK_deleteenv(&m_env);
479 }
480 
481 void CMosek::display_problem()
482 {
483  int32_t num_var, num_con;
484  m_rescode = MSK_getnumvar(m_task, &num_var);
485  m_rescode = MSK_getnumcon(m_task, &num_con);
486 
487  SG_PRINT("\nMatrix Q^0:\n")
488  for ( int32_t i = 0 ; i < num_var ; ++i )
489  {
490  for ( int32_t j = 0 ; j < num_var ; ++j )
491  {
492  float64_t qij;
493  m_rescode = MSK_getqobjij(m_task, i, j, &qij);
494  if ( qij != 0.0 )
495  SG_PRINT("(%d,%d)\t%.2f\n", i, j, qij)
496  }
497  }
498  SG_PRINT("\n")
499 
500  SG_PRINT("\nVector c:\n")
501  SGVector< float64_t > c(num_var);
502  m_rescode = MSK_getc(m_task, c.vector);
503  c.display_vector();
504 
505  SG_PRINT("\n\nMatrix A:\n")
506  for ( int32_t i = 0 ; i < num_con ; ++i )
507  {
508  for ( int32_t j = 0 ; j < num_var ; ++j )
509  {
510  float64_t aij;
511  m_rescode = MSK_getaij(m_task, i, j, &aij);
512  if ( aij != 0.0 )
513  SG_PRINT("(%d,%d)\t%.2f\n", i, j, aij)
514  }
515  }
516  SG_PRINT("\n")
517 
518  SG_PRINT("\nConstraint Bounds, vector b:\n")
519  for ( int32_t i = 0 ; i < num_con ; ++i )
520  {
521  MSKboundkeye bk;
522  float64_t bl, bu;
523  m_rescode = MSK_getbound(m_task, MSK_ACC_CON, i, &bk, &bl, &bu);
524 
525  SG_PRINT("%6.2f %6.2f\n", bl, bu)
526  }
527 
528  SG_PRINT("\nVariable Bounds, vectors lb and ub:\n")
529  for ( int32_t i = 0 ; i < num_var ; ++i )
530  {
531  MSKboundkeye bk;
532  float64_t bl, bu;
533  m_rescode = MSK_getbound(m_task, MSK_ACC_VAR, i, &bk, &bl, &bu);
534 
535  SG_PRINT("%6.2f %6.2f\n", bl, bu)
536  }
537  SG_PRINT("\n")
538 }
539 
540 
541 float64_t CMosek::get_primal_objective_value() const
542 {
543  float64_t po = 0.0;
544  MSK_getprimalobj(m_task, MSK_SOL_ITR, &po);
545 
546  return po;
547 }
548 
549 #endif /* USE_MOSEK */
static void fill_vector(T *vec, int32_t len, T value)
Definition: SGVector.cpp:221
void print(int depth, node< P > &top_node)
Definition: JLCoverTree.h:132
int32_t index_t
Definition: common.h:62
#define REQUIRE(x,...)
Definition: SGIO.h:206
index_t num_cols
Definition: SGMatrix.h:376
index_t num_rows
Definition: SGMatrix.h:374
void display_vector(const char *name="vector", const char *prefix="") const
Definition: SGVector.cpp:354
index_t vlen
Definition: SGVector.h:494
#define SG_PRINT(...)
Definition: SGIO.h:137
#define SG_SPRINT(...)
Definition: SGIO.h:180
#define ASSERT(x)
Definition: SGIO.h:201
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:115
double float64_t
Definition: common.h:50
static int32_t get_num_nonzero(T *vec, int32_t len)
Definition: Math.h:1164
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18

SHOGUN Machine Learning Toolbox - Documentation