SHOGUN
v3.2.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2012 Fernando José Iglesias García 00008 * Copyright (C) 2012 Fernando José Iglesias García 00009 */ 00010 00011 #ifdef USE_MOSEK 00012 00013 //#define DEBUG_MOSEK 00014 //#define DEBUG_SOLUTION 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 // Create MOSEK environment 00031 #if (MSK_VERSION_MAJOR == 6) 00032 m_rescode = MSK_makeenv(&m_env, NULL, NULL, NULL, NULL); 00033 #elif (MSK_VERSION_MAJOR == 7) 00034 m_rescode = MSK_makeenv(&m_env, NULL); 00035 #else 00036 #error "Unsupported Mosek version" 00037 #endif 00038 00039 #ifdef DEBUG_MOSEK 00040 // Direct the environment's log stream to SG_PRINT 00041 if ( m_rescode == MSK_RES_OK ) 00042 { 00043 m_rescode = MSK_linkfunctoenvstream(m_env, MSK_STREAM_LOG, 00044 NULL, CMosek::print); 00045 } 00046 #endif 00047 00048 // Initialize the environment 00049 if ( m_rescode == MSK_RES_OK ) 00050 { 00051 m_rescode = MSK_initenv(m_env); 00052 } 00053 00054 // Create an optimization task 00055 if ( m_rescode == MSK_RES_OK ) 00056 { 00057 m_rescode = MSK_maketask(m_env, num_con, num_var, &m_task); 00058 } 00059 00060 #ifdef DEBUG_MOSEK 00061 // Direct the task's log stream to SG_PRINT 00062 if ( m_rescode == MSK_RES_OK ) 00063 { 00064 m_rescode = MSK_linkfunctotaskstream(m_task, MSK_STREAM_LOG, 00065 NULL, CMosek::print); 00066 } 00067 #endif 00068 } 00069 00070 CMosek::~CMosek() 00071 { 00072 delete_problem(); 00073 } 00074 00075 void MSKAPI CMosek::print(void* handle, char str[]) 00076 { 00077 SG_SPRINT("%s", str) 00078 } 00079 00080 MSKrescodee CMosek::init_sosvm(int32_t M, int32_t N, 00081 int32_t num_aux, int32_t num_aux_con, 00082 SGMatrix< float64_t > C, SGVector< float64_t > lb, 00083 SGVector< float64_t > ub, SGMatrix< float64_t > A, 00084 SGVector< float64_t > b) 00085 { 00086 // Give an estimate of the size of the input data to increase the 00087 // speed of inputting 00088 int32_t num_var = M+N+num_aux; 00089 int32_t num_con = N*N+num_aux_con; 00090 // NOTE: However, to input this step is completely optional and MOSEK 00091 // will automatically allocate more resources if necessary 00092 m_rescode = MSK_putmaxnumvar(m_task, num_var); 00093 // Neither the number of constraints nor the number of non-zero elements 00094 // is known a priori, rough estimates are given here 00095 m_rescode = MSK_putmaxnumcon(m_task, num_con); 00096 // A = [-dPsi(y) | -I_N ] with M+N columns => max. M+1 nnz per row 00097 m_rescode = MSK_putmaxnumanz(m_task, (M+1)*N*N); 00098 00099 // Append optimization variables initialized to zero 00100 #if (MSK_VERSION_MAJOR == 6) 00101 m_rescode = MSK_append(m_task, MSK_ACC_VAR, num_var); 00102 #elif (MSK_VERSION_MAJOR == 7) 00103 m_rescode = MSK_appendvars(m_task, num_var); 00104 #else 00105 #error "Unsupported Mosek version" 00106 #endif 00107 // Append empty constraints initialized with no bounds 00108 #if (MSK_VERSION_MAJOR == 6) 00109 m_rescode = MSK_append(m_task, MSK_ACC_CON, num_con); 00110 #elif (MSK_VERSION_MAJOR == 7) 00111 m_rescode = MSK_appendcons(m_task, num_con); 00112 #else 00113 #error "Unsupported Mosek version" 00114 #endif 00115 // Set the constant term in the objective equal to zero 00116 m_rescode = MSK_putcfix(m_task, 0.0); 00117 00118 for ( int32_t j = 0 ; j < num_var && m_rescode == MSK_RES_OK ; ++j ) 00119 { 00120 // Set the linear term c_j in the objective 00121 if ( j < M+num_aux ) 00122 m_rescode = MSK_putcj(m_task, j, 0.0); 00123 else 00124 m_rescode = MSK_putcj(m_task, j, 1.0); 00125 00126 // Set the bounds on x_j: blx[j] <= x_j <= bux[j] 00127 // TODO set bounds lb and ub given by init_opt for w 00128 if ( j < M ) 00129 { 00130 m_rescode = MSK_putbound(m_task, MSK_ACC_VAR, j, 00131 MSK_BK_FR, -MSK_INFINITY, +MSK_INFINITY); 00132 } 00133 00134 // The slack and the auxiliary variables are required to be positive 00135 if ( j >= M ) 00136 { 00137 m_rescode = MSK_putbound(m_task, MSK_ACC_VAR, j, 00138 MSK_BK_LO, 0.0, +MSK_INFINITY); 00139 } 00140 } 00141 00142 // Input the matrix Q^0 for the objective 00143 // 00144 // NOTE: In MOSEK we minimize x'*Q^0*x. C != Q0 but Q0 is 00145 // just an extended version of C with zeros that make no 00146 // difference in MOSEK's sparse representation 00147 m_rescode = wrapper_putqobj(C); 00148 00149 // Input the matrix A and the vector b for the contraints A*x <= b 00150 m_rescode = wrapper_putaveclist(m_task, A); 00151 m_rescode = wrapper_putboundlist(m_task, b); 00152 00153 REQUIRE(m_rescode == MSK_RES_OK, "MOSEK Error in CMosek::init_sosvm(). " 00154 "Enable DEBUG_MOSEK for details.\n"); 00155 00156 return m_rescode; 00157 } 00158 00159 MSKrescodee CMosek::add_constraint_sosvm( 00160 SGVector< float64_t > dPsi, 00161 index_t con_idx, 00162 index_t train_idx, 00163 int32_t num_aux, 00164 float64_t bi) 00165 { 00166 // Count the number of non-zero element in dPsi 00167 int32_t nnz = CMath::get_num_nonzero(dPsi.vector, dPsi.vlen); 00168 // Indices of the non-zero elements in the row of A to add 00169 SGVector< index_t > asub(nnz+1); // +1 because of the -1 element 00170 // Values of the non-zero elements 00171 SGVector< float64_t > aval(nnz+1); 00172 // Next element to add in asub and aval 00173 index_t idx = 0; 00174 00175 for ( int32_t i = 0 ; i < dPsi.vlen ; ++i ) 00176 { 00177 if ( dPsi[i] != 0 ) 00178 { 00179 asub[idx] = i; 00180 aval[idx] = dPsi[i]; 00181 ++idx; 00182 } 00183 } 00184 00185 ASSERT(idx == nnz) 00186 00187 asub[idx] = dPsi.vlen + num_aux + train_idx; 00188 aval[idx] = -1; 00189 00190 #if (MSK_VERSION_MAJOR == 6) 00191 m_rescode = MSK_putavec(m_task, MSK_ACC_CON, con_idx, nnz+1, 00192 asub.vector, aval.vector); 00193 #elif (MSK_VERSION_MAJOR == 7) 00194 m_rescode = MSK_putarow(m_task, con_idx, nnz+1, asub.vector, aval.vector); 00195 #else 00196 #error "Unsupported Mosek version" 00197 #endif 00198 00199 if ( m_rescode == MSK_RES_OK ) 00200 { 00201 m_rescode = MSK_putbound(m_task, MSK_ACC_CON, con_idx, 00202 MSK_BK_UP, -MSK_INFINITY, bi); 00203 } 00204 00205 return m_rescode; 00206 } 00207 00208 MSKrescodee CMosek::wrapper_putaveclist( 00209 MSKtask_t & task, 00210 SGMatrix< float64_t > A) 00211 { 00212 // Indices to the rows of A to replace, all the rows 00213 SGVector< index_t > sub(A.num_rows); 00214 for ( index_t i = 0 ; i < A.num_rows ; ++i ) 00215 sub[i] = i; 00216 00217 // Non-zero elements of A 00218 int32_t nnza = CMath::get_num_nonzero(A.matrix, A.num_rows*A.num_cols); 00219 SGVector< float64_t > aval(nnza); 00220 // For each of the rows, indices to non-zero elements 00221 SGVector< index_t > asub(nnza); 00222 // For each row, pointer to the first non-zero element 00223 // in aval 00224 SGVector< int32_t > ptrb(A.num_rows); 00225 // Next position to write in aval and asub 00226 index_t idx = 0; 00227 // Switch if the first non-zero element in each row 00228 // has been found 00229 bool first_nnz_found = false; 00230 00231 for ( index_t i = 0 ; i < A.num_rows ; ++i ) 00232 { 00233 first_nnz_found = false; 00234 00235 for ( index_t j = 0 ; j < A.num_cols ; ++j ) 00236 { 00237 if ( A(i,j) ) 00238 { 00239 aval[idx] = A(i,j); 00240 asub[idx] = j; 00241 00242 if ( !first_nnz_found ) 00243 { 00244 ptrb[i] = idx; 00245 first_nnz_found = true; 00246 } 00247 00248 ++idx; 00249 } 00250 } 00251 00252 // Handle rows whose elements are all zero 00253 // TODO does it make sense that a row in A has all its elements 00254 // equal to zero? 00255 if ( !first_nnz_found ) 00256 ptrb[i] = ( i ? ptrb[i-1] : 0 ); 00257 } 00258 00259 // For each row, pointer to the last+1 non-zero element 00260 // in aval 00261 SGVector< int32_t > ptre(A.num_rows); 00262 for ( index_t i = 0 ; i < A.num_rows-1 ; ++i ) 00263 ptre[i] = ptrb[i+1]; 00264 00265 if ( A.num_rows > 0 ) 00266 ptre[A.num_rows-1] = nnza; 00267 00268 MSKrescodee ret; 00269 #if (MSK_VERSION_MAJOR == 6) 00270 ret = MSK_putaveclist(task, MSK_ACC_CON, A.num_rows, sub.vector, 00271 ptrb.vector, ptre.vector, 00272 asub.vector, aval.vector); 00273 #elif (MSK_VERSION_MAJOR == 7) 00274 ret = MSK_putarowlist(task, A.num_rows, sub.vector, ptrb.vector, ptre.vector, 00275 asub.vector, aval.vector); 00276 #else 00277 #error "Unsupported Mosek version" 00278 #endif 00279 00280 REQUIRE(ret == MSK_RES_OK, "MOSEK Error in CMosek::wrapper_putaveclist(). " 00281 "Enable DEBUG_MOSEK for details.\n"); 00282 00283 return ret; 00284 } 00285 00286 MSKrescodee CMosek::wrapper_putboundlist(MSKtask_t & task, SGVector< float64_t > b) 00287 { 00288 // Indices to the bounds that should be replaced, b.vlen bounds starting 00289 // from zero 00290 SGVector< index_t > sub(b.vlen); 00291 for ( index_t i = 0 ; i < b.vlen ; ++i ) 00292 sub[i] = i; 00293 00294 // Type of the bounds and lower bound values 00295 MSKboundkeye* bk = SG_MALLOC(MSKboundkeye, b.vlen); 00296 SGVector< float64_t > bl(b.vlen); 00297 for ( index_t i = 0 ; i < b.vlen ; ++i ) 00298 { 00299 bk[i] = MSK_BK_UP; 00300 bl[i] = -MSK_INFINITY; 00301 } 00302 00303 MSKrescodee ret = MSK_putboundlist(task, MSK_ACC_CON, b.vlen, sub.vector, 00304 bk, bl.vector, b.vector); 00305 00306 SG_FREE(bk); 00307 00308 REQUIRE(ret == MSK_RES_OK, "MOSEK Error in CMosek::wrapper_putboundlist(). " 00309 "Enable DEBUG_MOSEK for details.\n"); 00310 00311 return ret; 00312 } 00313 00314 MSKrescodee CMosek::wrapper_putqobj(SGMatrix< float64_t > Q0) const 00315 { 00316 // Shorthands for the dimensions of the matrix 00317 index_t N = Q0.num_rows; 00318 index_t M = Q0.num_cols; 00319 00320 // Count the number of non-zero elements in the lower 00321 // triangular part of the matrix 00322 int32_t nnz = 0; 00323 for ( index_t i = 0 ; i < N ; ++i ) 00324 for ( index_t j = i ; j < M ; ++j ) 00325 nnz += ( Q0[j + i*M] ? 1 : 0 ); 00326 00327 // i subscript for non-zero elements of Q0 00328 SGVector< index_t > qosubi(nnz); 00329 // j subscript for non-zero elements of Q0 00330 SGVector< index_t > qosubj(nnz); 00331 // Values of non-zero elements of Q0 00332 SGVector< float64_t > qoval(nnz); 00333 // Next position to write in the vectors 00334 index_t idx = 0; 00335 00336 for ( index_t i = 0 ; i < N ; ++i ) 00337 for ( index_t j = i ; j < M ; ++j ) 00338 { 00339 if ( Q0[j + i*M] ) 00340 { 00341 qosubi[idx] = i; 00342 qosubj[idx] = j; 00343 qoval[idx] = Q0[j + i*M]; 00344 00345 ++idx; 00346 } 00347 } 00348 00349 return MSK_putqobj(m_task, nnz, qosubi.vector, 00350 qosubj.vector, qoval.vector); 00351 } 00352 00353 MSKrescodee CMosek::optimize(SGVector< float64_t > sol) 00354 { 00355 m_rescode = MSK_optimize(m_task); 00356 00357 #ifdef DEBUG_MOSEK 00358 // Print a summary containing information about the solution 00359 MSK_solutionsummary(m_task, MSK_STREAM_LOG); 00360 #endif 00361 00362 // Read the solution 00363 if ( m_rescode == MSK_RES_OK ) 00364 { 00365 // Solution status 00366 MSKsolstae solsta; 00367 // FIXME posible solutions are: 00368 // MSK_SOL_ITR: the interior solution 00369 // MSK_SOL_BAS: the basic solution 00370 // MSK_SOL_ITG: the integer solution 00371 #if (MSK_VERSION_MAJOR == 6) 00372 MSK_getsolutionstatus(m_task, MSK_SOL_ITR, NULL, &solsta); 00373 #elif (MSK_VERSION_MAJOR == 7) 00374 MSK_getsolsta(m_task, MSK_SOL_ITR, &solsta); 00375 #else 00376 #error "Unsupported Mosek Version" 00377 #endif 00378 switch (solsta) 00379 { 00380 case MSK_SOL_STA_OPTIMAL: 00381 case MSK_SOL_STA_NEAR_OPTIMAL: 00382 MSK_getsolutionslice(m_task, 00383 // Request the interior solution 00384 MSK_SOL_ITR, 00385 // of the optimization vector 00386 MSK_SOL_ITEM_XX, 00387 0, 00388 sol.vlen, 00389 sol.vector); 00390 #ifdef DEBUG_SOLUTION 00391 sol.display_vector("Solution"); 00392 #endif 00393 break; 00394 case MSK_SOL_STA_DUAL_INFEAS_CER: 00395 case MSK_SOL_STA_PRIM_INFEAS_CER: 00396 case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER: 00397 case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER: 00398 #ifdef DEBUG_MOSEK 00399 SG_PRINT("Primal or dual infeasibility " 00400 "certificate found\n"); 00401 #endif 00402 break; 00403 case MSK_SOL_STA_UNKNOWN: 00404 #ifdef DEBUG_MOSEK 00405 SG_PRINT("Undetermined solution status\n") 00406 #endif 00407 break; 00408 default: 00409 #ifdef DEBUG_MOSEK 00410 SG_PRINT("Other solution status\n") 00411 #endif 00412 break; // to avoid compile error when DEBUG_MOSEK 00413 // is not defined 00414 } 00415 } 00416 00417 // In case any error occurred, print the appropriate error message 00418 if ( m_rescode != MSK_RES_OK ) 00419 { 00420 char symname[MSK_MAX_STR_LEN]; 00421 char desc[MSK_MAX_STR_LEN]; 00422 00423 MSK_getcodedesc(m_rescode, symname, desc); 00424 00425 SG_PRINT("An error occurred optimizing with MOSEK\n") 00426 SG_PRINT("ERROR %s - '%s'\n", symname, desc) 00427 } 00428 00429 return m_rescode; 00430 } 00431 00432 void CMosek::delete_problem() 00433 { 00434 MSK_deletetask(&m_task); 00435 MSK_deleteenv(&m_env); 00436 } 00437 00438 void CMosek::display_problem() 00439 { 00440 int32_t num_var, num_con; 00441 m_rescode = MSK_getnumvar(m_task, &num_var); 00442 m_rescode = MSK_getnumcon(m_task, &num_con); 00443 00444 SG_PRINT("\nMatrix Q^0:\n") 00445 for ( int32_t i = 0 ; i < num_var ; ++i ) 00446 { 00447 for ( int32_t j = 0 ; j < num_var ; ++j ) 00448 { 00449 float64_t qij; 00450 m_rescode = MSK_getqobjij(m_task, i, j, &qij); 00451 if ( qij != 0.0 ) 00452 SG_PRINT("(%d,%d)\t%.2f\n", i, j, qij) 00453 } 00454 } 00455 SG_PRINT("\n") 00456 00457 SG_PRINT("\nVector c:\n") 00458 SGVector< float64_t > c(num_var); 00459 m_rescode = MSK_getc(m_task, c.vector); 00460 c.display_vector(); 00461 00462 SG_PRINT("\n\nMatrix A:\n") 00463 for ( int32_t i = 0 ; i < num_con ; ++i ) 00464 { 00465 for ( int32_t j = 0 ; j < num_var ; ++j ) 00466 { 00467 float64_t aij; 00468 m_rescode = MSK_getaij(m_task, i, j, &aij); 00469 if ( aij != 0.0 ) 00470 SG_PRINT("(%d,%d)\t%.2f\n", i, j, aij) 00471 } 00472 } 00473 SG_PRINT("\n") 00474 00475 SG_PRINT("\nConstraint Bounds, vector b:\n") 00476 for ( int32_t i = 0 ; i < num_con ; ++i ) 00477 { 00478 MSKboundkeye bk; 00479 float64_t bl, bu; 00480 m_rescode = MSK_getbound(m_task, MSK_ACC_CON, i, &bk, &bl, &bu); 00481 00482 SG_PRINT("%6.2f %6.2f\n", bl, bu) 00483 } 00484 00485 SG_PRINT("\nVariable Bounds, vectors lb and ub:\n") 00486 for ( int32_t i = 0 ; i < num_var ; ++i ) 00487 { 00488 MSKboundkeye bk; 00489 float64_t bl, bu; 00490 m_rescode = MSK_getbound(m_task, MSK_ACC_VAR, i, &bk, &bl, &bu); 00491 00492 SG_PRINT("%6.2f %6.2f\n", bl, bu) 00493 } 00494 SG_PRINT("\n") 00495 } 00496 00497 00498 float64_t CMosek::get_primal_objective_value() const 00499 { 00500 float64_t po = 0.0; 00501 MSK_getprimalobj(m_task, MSK_SOL_ITR, &po); 00502 00503 return po; 00504 } 00505 00506 #endif /* USE_MOSEK */