SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Mosek.cpp
Go to the documentation of this file.
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 */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation