00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #ifdef USE_MOSEK
00012
00013
00014
00015
00016 #include <shogun/mathematics/Math.h>
00017 #include <shogun/mathematics/Mosek.h>
00018 #include <shogun/lib/SGVector.h>
00019
00020 using namespace shogun;
00021
00022 CMosek::CMosek()
00023 : CSGObject()
00024 {
00025 }
00026
00027 CMosek::CMosek(int32_t num_con, int32_t num_var)
00028 : CSGObject()
00029 {
00030
00031 m_rescode = MSK_makeenv(&m_env, NULL, NULL, NULL, NULL);
00032
00033 #ifdef DEBUG_MOSEK
00034
00035 if ( m_rescode == MSK_RES_OK )
00036 {
00037 m_rescode = MSK_linkfunctoenvstream(m_env, MSK_STREAM_LOG,
00038 NULL, CMosek::print);
00039 }
00040 #endif
00041
00042
00043 if ( m_rescode == MSK_RES_OK )
00044 {
00045 m_rescode = MSK_initenv(m_env);
00046 }
00047
00048
00049 if ( m_rescode == MSK_RES_OK )
00050 {
00051 m_rescode = MSK_maketask(m_env, num_con, num_var, &m_task);
00052 }
00053
00054 #ifdef DEBUG_MOSEK
00055
00056 if ( m_rescode == MSK_RES_OK )
00057 {
00058 m_rescode = MSK_linkfunctotaskstream(m_task, MSK_STREAM_LOG,
00059 NULL, CMosek::print);
00060 }
00061 #endif
00062 }
00063
00064 CMosek::~CMosek()
00065 {
00066 delete_problem();
00067 }
00068
00069 void MSKAPI CMosek::print(void* handle, char str[])
00070 {
00071 SG_SPRINT("%s", str);
00072 }
00073
00074 MSKrescodee CMosek::init_sosvm(int32_t M, int32_t N,
00075 int32_t num_aux, int32_t num_aux_con,
00076 SGMatrix< float64_t > C, SGVector< float64_t > lb,
00077 SGVector< float64_t > ub, SGMatrix< float64_t > A,
00078 SGVector< float64_t > b)
00079 {
00080
00081
00082 int32_t num_var = M+N+num_aux;
00083 int32_t num_con = N*N+num_aux_con;
00084
00085
00086 m_rescode = MSK_putmaxnumvar(m_task, num_var);
00087
00088
00089 m_rescode = MSK_putmaxnumcon(m_task, num_con);
00090
00091 m_rescode = MSK_putmaxnumanz(m_task, (M+1)*N*N);
00092
00093
00094 m_rescode = MSK_append(m_task, MSK_ACC_VAR, num_var);
00095
00096 m_rescode = MSK_append(m_task, MSK_ACC_CON, num_con);
00097
00098 m_rescode = MSK_putcfix(m_task, 0.0);
00099
00100 for ( int32_t j = 0 ; j < num_var && m_rescode == MSK_RES_OK ; ++j )
00101 {
00102
00103 if ( j < M+num_aux )
00104 m_rescode = MSK_putcj(m_task, j, 0.0);
00105 else
00106 m_rescode = MSK_putcj(m_task, j, 1.0);
00107
00108
00109
00110 if ( j < M )
00111 {
00112 m_rescode = MSK_putbound(m_task, MSK_ACC_VAR, j,
00113 MSK_BK_FR, -MSK_INFINITY, +MSK_INFINITY);
00114 }
00115
00116
00117 if ( j >= M )
00118 {
00119 m_rescode = MSK_putbound(m_task, MSK_ACC_VAR, j,
00120 MSK_BK_LO, 0.0, +MSK_INFINITY);
00121 }
00122 }
00123
00124
00125
00126
00127
00128
00129 m_rescode = wrapper_putqobj(C);
00130
00131
00132 m_rescode = wrapper_putaveclist(m_task, A);
00133 m_rescode = wrapper_putboundlist(m_task, b);
00134
00135 REQUIRE(m_rescode == MSK_RES_OK, "MOSEK Error in CMosek::init_sosvm(). "
00136 "Enable DEBUG_MOSEK for details.\n");
00137
00138 return m_rescode;
00139 }
00140
00141 MSKrescodee CMosek::add_constraint_sosvm(
00142 SGVector< float64_t > dPsi,
00143 index_t con_idx,
00144 index_t train_idx,
00145 int32_t num_aux,
00146 float64_t bi)
00147 {
00148
00149 int32_t nnz = CMath::get_num_nonzero(dPsi.vector, dPsi.vlen);
00150
00151 SGVector< index_t > asub(nnz+1);
00152
00153 SGVector< float64_t > aval(nnz+1);
00154
00155 index_t idx = 0;
00156
00157 for ( int32_t i = 0 ; i < dPsi.vlen ; ++i )
00158 {
00159 if ( dPsi[i] != 0 )
00160 {
00161 asub[idx] = i;
00162 aval[idx] = dPsi[i];
00163 ++idx;
00164 }
00165 }
00166
00167 ASSERT(idx == nnz);
00168
00169 asub[idx] = dPsi.vlen + num_aux + train_idx;
00170 aval[idx] = -1;
00171
00172 m_rescode = MSK_putavec(m_task, MSK_ACC_CON, con_idx, nnz+1,
00173 asub.vector, aval.vector);
00174
00175 if ( m_rescode == MSK_RES_OK )
00176 {
00177 m_rescode = MSK_putbound(m_task, MSK_ACC_CON, con_idx,
00178 MSK_BK_UP, -MSK_INFINITY, bi);
00179 }
00180
00181 return m_rescode;
00182 }
00183
00184 MSKrescodee CMosek::wrapper_putaveclist(
00185 MSKtask_t & task,
00186 SGMatrix< float64_t > A)
00187 {
00188
00189 SGVector< index_t > sub(A.num_rows);
00190 for ( index_t i = 0 ; i < A.num_rows ; ++i )
00191 sub[i] = i;
00192
00193
00194 int32_t nnza = CMath::get_num_nonzero(A.matrix, A.num_rows*A.num_cols);
00195 SGVector< float64_t > aval(nnza);
00196
00197 SGVector< index_t > asub(nnza);
00198
00199
00200 SGVector< int32_t > ptrb(A.num_rows);
00201
00202 index_t idx = 0;
00203
00204
00205 bool first_nnz_found = false;
00206
00207 for ( index_t i = 0 ; i < A.num_rows ; ++i )
00208 {
00209 first_nnz_found = false;
00210
00211 for ( index_t j = 0 ; j < A.num_cols ; ++j )
00212 {
00213 if ( A(i,j) )
00214 {
00215 aval[idx] = A(i,j);
00216 asub[idx] = j;
00217
00218 if ( !first_nnz_found )
00219 {
00220 ptrb[i] = idx;
00221 first_nnz_found = true;
00222 }
00223
00224 ++idx;
00225 }
00226 }
00227
00228
00229
00230
00231 if ( !first_nnz_found )
00232 ptrb[i] = ( i ? ptrb[i-1] : 0 );
00233 }
00234
00235
00236
00237 SGVector< int32_t > ptre(A.num_rows);
00238 for ( index_t i = 0 ; i < A.num_rows-1 ; ++i )
00239 ptre[i] = ptrb[i+1];
00240
00241 if ( A.num_rows > 0 )
00242 ptre[A.num_rows-1] = nnza;
00243
00244 MSKrescodee ret = MSK_putaveclist(task, MSK_ACC_CON, A.num_rows, sub.vector,
00245 ptrb.vector, ptre.vector,
00246 asub.vector, aval.vector);
00247
00248 REQUIRE(ret == MSK_RES_OK, "MOSEK Error in CMosek::wrapper_putaveclist(). "
00249 "Enable DEBUG_MOSEK for details.\n");
00250
00251 return ret;
00252 }
00253
00254 MSKrescodee CMosek::wrapper_putboundlist(MSKtask_t & task, SGVector< float64_t > b)
00255 {
00256
00257
00258 SGVector< index_t > sub(b.vlen);
00259 for ( index_t i = 0 ; i < b.vlen ; ++i )
00260 sub[i] = i;
00261
00262
00263 MSKboundkeye* bk = SG_MALLOC(MSKboundkeye, b.vlen);
00264 SGVector< float64_t > bl(b.vlen);
00265 for ( index_t i = 0 ; i < b.vlen ; ++i )
00266 {
00267 bk[i] = MSK_BK_UP;
00268 bl[i] = -MSK_INFINITY;
00269 }
00270
00271 MSKrescodee ret = MSK_putboundlist(task, MSK_ACC_CON, b.vlen, sub.vector,
00272 bk, bl.vector, b.vector);
00273
00274 SG_FREE(bk);
00275
00276 REQUIRE(ret == MSK_RES_OK, "MOSEK Error in CMosek::wrapper_putboundlist(). "
00277 "Enable DEBUG_MOSEK for details.\n");
00278
00279 return ret;
00280 }
00281
00282 MSKrescodee CMosek::wrapper_putqobj(SGMatrix< float64_t > Q0) const
00283 {
00284
00285 index_t N = Q0.num_rows;
00286 index_t M = Q0.num_cols;
00287
00288
00289
00290 int32_t nnz = 0;
00291 for ( index_t i = 0 ; i < N ; ++i )
00292 for ( index_t j = i ; j < M ; ++j )
00293 nnz += ( Q0[j + i*M] ? 1 : 0 );
00294
00295
00296 SGVector< index_t > qosubi(nnz);
00297
00298 SGVector< index_t > qosubj(nnz);
00299
00300 SGVector< float64_t > qoval(nnz);
00301
00302 index_t idx = 0;
00303
00304 for ( index_t i = 0 ; i < N ; ++i )
00305 for ( index_t j = i ; j < M ; ++j )
00306 {
00307 if ( Q0[j + i*M] )
00308 {
00309 qosubi[idx] = i;
00310 qosubj[idx] = j;
00311 qoval[idx] = Q0[j + i*M];
00312
00313 ++idx;
00314 }
00315 }
00316
00317 return MSK_putqobj(m_task, nnz, qosubi.vector,
00318 qosubj.vector, qoval.vector);
00319 }
00320
00321 MSKrescodee CMosek::optimize(SGVector< float64_t > sol)
00322 {
00323 m_rescode = MSK_optimize(m_task);
00324
00325 #ifdef DEBUG_MOSEK
00326
00327 MSK_solutionsummary(m_task, MSK_STREAM_LOG);
00328 #endif
00329
00330
00331 if ( m_rescode == MSK_RES_OK )
00332 {
00333
00334 MSKsolstae solsta;
00335
00336
00337
00338
00339 MSK_getsolutionstatus(m_task, MSK_SOL_ITR, NULL, &solsta);
00340
00341 switch (solsta)
00342 {
00343 case MSK_SOL_STA_OPTIMAL:
00344 case MSK_SOL_STA_NEAR_OPTIMAL:
00345 MSK_getsolutionslice(m_task,
00346
00347 MSK_SOL_ITR,
00348
00349 MSK_SOL_ITEM_XX,
00350 0,
00351 sol.vlen,
00352 sol.vector);
00353 #ifdef DEBUG_SOLUTION
00354 sol.display_vector("Solution");
00355 #endif
00356 break;
00357 case MSK_SOL_STA_DUAL_INFEAS_CER:
00358 case MSK_SOL_STA_PRIM_INFEAS_CER:
00359 case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER:
00360 case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER:
00361 #ifdef DEBUG_MOSEK
00362 SG_PRINT("Primal or dual infeasibility "
00363 "certificate found\n");
00364 #endif
00365 break;
00366 case MSK_SOL_STA_UNKNOWN:
00367 #ifdef DEBUG_MOSEK
00368 SG_PRINT("Undetermined solution status\n");
00369 #endif
00370 break;
00371 default:
00372 #ifdef DEBUG_MOSEK
00373 SG_PRINT("Other solution status\n");
00374 #endif
00375 break;
00376
00377 }
00378 }
00379
00380
00381 if ( m_rescode != MSK_RES_OK )
00382 {
00383 char symname[MSK_MAX_STR_LEN];
00384 char desc[MSK_MAX_STR_LEN];
00385
00386 MSK_getcodedesc(m_rescode, symname, desc);
00387
00388 SG_PRINT("An error occurred optimizing with MOSEK\n");
00389 SG_PRINT("ERROR %s - '%s'\n", symname, desc);
00390 }
00391
00392 return m_rescode;
00393 }
00394
00395 void CMosek::delete_problem()
00396 {
00397 MSK_deletetask(&m_task);
00398 MSK_deleteenv(&m_env);
00399 }
00400
00401 void CMosek::display_problem()
00402 {
00403 int32_t num_var, num_con;
00404 m_rescode = MSK_getnumvar(m_task, &num_var);
00405 m_rescode = MSK_getnumcon(m_task, &num_con);
00406
00407 SG_PRINT("\nMatrix Q^0:\n");
00408 for ( int32_t i = 0 ; i < num_var ; ++i )
00409 {
00410 for ( int32_t j = 0 ; j < num_var ; ++j )
00411 {
00412 float64_t qij;
00413 m_rescode = MSK_getqobjij(m_task, i, j, &qij);
00414 if ( qij != 0.0 )
00415 SG_PRINT("(%d,%d)\t%.2f\n", i, j, qij);
00416 }
00417 }
00418 SG_PRINT("\n");
00419
00420 SG_PRINT("\nVector c:\n");
00421 SGVector< float64_t > c(num_var);
00422 m_rescode = MSK_getc(m_task, c.vector);
00423 c.display_vector();
00424
00425 SG_PRINT("\n\nMatrix A:\n");
00426 for ( int32_t i = 0 ; i < num_con ; ++i )
00427 {
00428 for ( int32_t j = 0 ; j < num_var ; ++j )
00429 {
00430 float64_t aij;
00431 m_rescode = MSK_getaij(m_task, i, j, &aij);
00432 if ( aij != 0.0 )
00433 SG_PRINT("(%d,%d)\t%.2f\n", i, j, aij);
00434 }
00435 }
00436 SG_PRINT("\n");
00437
00438 SG_PRINT("\nConstraint Bounds, vector b:\n");
00439 for ( int32_t i = 0 ; i < num_con ; ++i )
00440 {
00441 MSKboundkeye bk;
00442 float64_t bl, bu;
00443 m_rescode = MSK_getbound(m_task, MSK_ACC_CON, i, &bk, &bl, &bu);
00444
00445 SG_PRINT("%6.2f %6.2f\n", bl, bu);
00446 }
00447
00448 SG_PRINT("\nVariable Bounds, vectors lb and ub:\n");
00449 for ( int32_t i = 0 ; i < num_var ; ++i )
00450 {
00451 MSKboundkeye bk;
00452 float64_t bl, bu;
00453 m_rescode = MSK_getbound(m_task, MSK_ACC_VAR, i, &bk, &bl, &bu);
00454
00455 SG_PRINT("%6.2f %6.2f\n", bl, bu);
00456 }
00457 SG_PRINT("\n");
00458 }
00459
00460
00461 float64_t CMosek::get_primal_objective_value() const
00462 {
00463 float64_t po = 0.0;
00464 MSK_getprimalobj(m_task, MSK_SOL_ITR, &po);
00465
00466 return po;
00467 }
00468
00469 #endif