SHOGUN
v2.0.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 m_rescode = MSK_makeenv(&m_env, NULL, NULL, NULL, NULL); 00032 00033 #ifdef DEBUG_MOSEK 00034 // Direct the environment's log stream to SG_PRINT 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 // Initialize the environment 00043 if ( m_rescode == MSK_RES_OK ) 00044 { 00045 m_rescode = MSK_initenv(m_env); 00046 } 00047 00048 // Create an optimization task 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 // Direct the task's log stream to SG_PRINT 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 // Give an estimate of the size of the input data to increase the 00081 // speed of inputting 00082 int32_t num_var = M+N+num_aux; 00083 int32_t num_con = N*N+num_aux_con; 00084 // NOTE: However, to input this step is completely optional and MOSEK 00085 // will automatically allocate more resources if necessary 00086 m_rescode = MSK_putmaxnumvar(m_task, num_var); 00087 // Neither the number of constraints nor the number of non-zero elements 00088 // is known a priori, rough estimates are given here 00089 m_rescode = MSK_putmaxnumcon(m_task, num_con); 00090 // A = [-dPsi(y) | -I_N ] with M+N columns => max. M+1 nnz per row 00091 m_rescode = MSK_putmaxnumanz(m_task, (M+1)*N*N); 00092 00093 // Append optimization variables initialized to zero 00094 m_rescode = MSK_append(m_task, MSK_ACC_VAR, num_var); 00095 // Append empty constraints initialized with no bounds 00096 m_rescode = MSK_append(m_task, MSK_ACC_CON, num_con); 00097 // Set the constant term in the objective equal to zero 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 // Set the linear term c_j in the objective 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 // Set the bounds on x_j: blx[j] <= x_j <= bux[j] 00109 // TODO set bounds lb and ub given by init_opt for w 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 // The slack and the auxiliary variables are required to be positive 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 // Input the matrix Q^0 for the objective 00125 // 00126 // NOTE: In MOSEK we minimize x'*Q^0*x. C != Q0 but Q0 is 00127 // just an extended version of C with zeros that make no 00128 // difference in MOSEK's sparse representation 00129 m_rescode = wrapper_putqobj(C); 00130 00131 // Input the matrix A and the vector b for the contraints A*x <= b 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 // Count the number of non-zero element in dPsi 00149 int32_t nnz = CMath::get_num_nonzero(dPsi.vector, dPsi.vlen); 00150 // Indices of the non-zero elements in the row of A to add 00151 SGVector< index_t > asub(nnz+1); // +1 because of the -1 element 00152 // Values of the non-zero elements 00153 SGVector< float64_t > aval(nnz+1); 00154 // Next element to add in asub and aval 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 // Indices to the rows of A to replace, all the rows 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 // Non-zero elements of A 00194 int32_t nnza = CMath::get_num_nonzero(A.matrix, A.num_rows*A.num_cols); 00195 SGVector< float64_t > aval(nnza); 00196 // For each of the rows, indices to non-zero elements 00197 SGVector< index_t > asub(nnza); 00198 // For each row, pointer to the first non-zero element 00199 // in aval 00200 SGVector< int32_t > ptrb(A.num_rows); 00201 // Next position to write in aval and asub 00202 index_t idx = 0; 00203 // Switch if the first non-zero element in each row 00204 // has been found 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 // Handle rows whose elements are all zero 00229 // TODO does it make sense that a row in A has all its elements 00230 // equal to zero? 00231 if ( !first_nnz_found ) 00232 ptrb[i] = ( i ? ptrb[i-1] : 0 ); 00233 } 00234 00235 // For each row, pointer to the last+1 non-zero element 00236 // in aval 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 // Indices to the bounds that should be replaced, b.vlen bounds starting 00257 // from zero 00258 SGVector< index_t > sub(b.vlen); 00259 for ( index_t i = 0 ; i < b.vlen ; ++i ) 00260 sub[i] = i; 00261 00262 // Type of the bounds and lower bound values 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 // Shorthands for the dimensions of the matrix 00285 index_t N = Q0.num_rows; 00286 index_t M = Q0.num_cols; 00287 00288 // Count the number of non-zero elements in the lower 00289 // triangular part of the matrix 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 // i subscript for non-zero elements of Q0 00296 SGVector< index_t > qosubi(nnz); 00297 // j subscript for non-zero elements of Q0 00298 SGVector< index_t > qosubj(nnz); 00299 // Values of non-zero elements of Q0 00300 SGVector< float64_t > qoval(nnz); 00301 // Next position to write in the vectors 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 // Print a summary containing information about the solution 00327 MSK_solutionsummary(m_task, MSK_STREAM_LOG); 00328 #endif 00329 00330 // Read the solution 00331 if ( m_rescode == MSK_RES_OK ) 00332 { 00333 // Solution status 00334 MSKsolstae solsta; 00335 // FIXME posible solutions are: 00336 // MSK_SOL_ITR: the interior solution 00337 // MSK_SOL_BAS: the basic solution 00338 // MSK_SOL_ITG: the integer solution 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 // Request the interior solution 00347 MSK_SOL_ITR, 00348 // of the optimization vector 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; // to avoid compile error when DEBUG_MOSEK 00376 // is not defined 00377 } 00378 } 00379 00380 // In case any error occurred, print the appropriate error message 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 #endif /* USE_MOSEK */