00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "lib/config.h"
00012
00013 #ifdef USE_CPLEX
00014 #include <unistd.h>
00015
00016 #include "lib/Cplex.h"
00017 #include "lib/io.h"
00018 #include "lib/Signal.h"
00019 #include "lib/Mathematics.h"
00020
00021 CCplex::CCplex()
00022 : CSGObject(), env(NULL), lp(NULL), lp_initialized(false)
00023 {
00024 }
00025
00026 CCplex::~CCplex()
00027 {
00028 cleanup();
00029 }
00030
00031 bool CCplex::init(E_PROB_TYPE typ, int32_t timeout)
00032 {
00033 problem_type=typ;
00034
00035 while (env==NULL)
00036 {
00037 int status = 1;
00038 env = CPXopenCPLEX (&status);
00039
00040 if ( env == NULL )
00041 {
00042 char errmsg[1024];
00043 SG_WARNING("Could not open CPLEX environment.\n");
00044 CPXgeterrorstring (env, status, errmsg);
00045 SG_WARNING("%s", errmsg);
00046 SG_WARNING("retrying in %d seconds\n", timeout);
00047 sleep(timeout);
00048 }
00049 else
00050 {
00051
00052
00053 status = CPXsetintparam (env, CPX_PARAM_SCRIND, CPX_OFF);
00054 if (status)
00055 SG_ERROR( "Failure to turn off screen indicator, error %d.\n", status);
00056
00057 {
00058 status = CPXsetintparam (env, CPX_PARAM_DATACHECK, CPX_ON);
00059 if (status)
00060 SG_ERROR( "Failure to turn on data checking, error %d.\n", status);
00061 else
00062 {
00063 lp = CPXcreateprob (env, &status, "shogun");
00064
00065 if ( lp == NULL )
00066 SG_ERROR( "Failed to create optimization problem.\n");
00067 else
00068 CPXchgobjsen (env, lp, CPX_MIN);
00069
00070 if (problem_type==E_QP)
00071 status = CPXsetintparam (env, CPX_PARAM_QPMETHOD, 0);
00072 else if (problem_type == E_LINEAR)
00073 status = CPXsetintparam (env, CPX_PARAM_LPMETHOD, 0);
00074 if (status)
00075 SG_ERROR( "Failure to select dual lp/qp optimization, error %d.\n", status);
00076
00077 }
00078 }
00079 }
00080 }
00081
00082 return (lp != NULL) && (env != NULL);
00083 }
00084
00085 bool CCplex::setup_subgradientlpm_QP(
00086 float64_t C, CLabels* labels, CSparseFeatures<float64_t>* features,
00087 int32_t* idx_bound, int32_t num_bound, int32_t* w_zero, int32_t num_zero,
00088 float64_t* vee, int32_t num_dim, bool use_bias)
00089 {
00090 const int cmatsize=10*1024*1024;
00091 bool result=false;
00092 int32_t num_variables=num_dim + num_bound+num_zero;
00093
00094 ASSERT(num_dim>0);
00095 ASSERT(num_dim>0);
00096
00097
00098 float64_t* lb=new float64_t[num_variables];
00099 float64_t* ub=new float64_t[num_variables];
00100 float64_t* obj=new float64_t[num_variables];
00101 char* sense = new char[num_dim];
00102 int* cmatbeg=new int[num_variables];
00103 int* cmatcnt=new int[num_variables];
00104 int* cmatind=new int[cmatsize];
00105 double* cmatval=new double[cmatsize];
00106
00107 for (int32_t i=0; i<num_variables; i++)
00108 {
00109 obj[i]=0;
00110
00111 if (i<num_dim)
00112 {
00113 lb[i]=-CPX_INFBOUND;
00114 ub[i]=+CPX_INFBOUND;
00115 }
00116 else
00117 {
00118 lb[i]=0.0;
00119 ub[i]=1.0;
00120 }
00121 }
00122
00123 int32_t offs=0;
00124 for (int32_t i=0; i<num_variables; i++)
00125 {
00126 if (i<num_dim)
00127 {
00128 sense[i]='E';
00129 cmatbeg[i]=offs;
00130 cmatcnt[i]=1;
00131 cmatind[offs]=offs;
00132 cmatval[offs]=1.0;
00133 offs++;
00134 ASSERT(offs<cmatsize);
00135 }
00136 else if (i<num_dim+num_zero)
00137 {
00138 cmatbeg[i]=offs;
00139 cmatcnt[i]=1;
00140 cmatind[offs]=w_zero[i-num_dim];
00141 cmatval[offs]=-1.0;
00142 offs++;
00143 ASSERT(offs<cmatsize);
00144 }
00145 else
00146 {
00147 int32_t idx=idx_bound[i-num_dim-num_zero];
00148 int32_t vlen=0;
00149 bool vfree=false;
00150
00151 TSparseEntry<float64_t>* vec=features->get_sparse_feature_vector(idx, vlen, vfree);
00152
00153
00154 cmatbeg[i]=offs;
00155 cmatcnt[i]=vlen;
00156
00157 float64_t val= -C*labels->get_label(idx);
00158
00159 if (vlen>0)
00160 {
00161 for (int32_t j=0; j<vlen; j++)
00162 {
00163 cmatind[offs]=vec[j].feat_index;
00164 cmatval[offs]=-val*vec[j].entry;
00165 offs++;
00166 ASSERT(offs<cmatsize);
00167
00168 }
00169
00170 if (use_bias)
00171 {
00172 cmatcnt[i]++;
00173 cmatind[offs]=num_dim-1;
00174 cmatval[offs]=-val;
00175 offs++;
00176 ASSERT(offs<cmatsize);
00177 }
00178 }
00179 else
00180 {
00181 if (use_bias)
00182 {
00183 cmatcnt[i]++;
00184 cmatind[offs]=num_dim-1;
00185 cmatval[offs]=-val;
00186 }
00187 else
00188 cmatval[offs]=0.0;
00189 offs++;
00190 ASSERT(offs<cmatsize);
00191 }
00192
00193 features->free_feature_vector(vec, idx, vfree);
00194 }
00195 }
00196
00197 result = CPXcopylp(env, lp, num_variables, num_dim, CPX_MIN,
00198 obj, vee, sense, cmatbeg, cmatcnt, cmatind, cmatval, lb, ub, NULL) == 0;
00199
00200 if (!result)
00201 SG_ERROR("CPXcopylp failed.\n");
00202
00203
00204
00205 delete[] sense;
00206 delete[] lb;
00207 delete[] ub;
00208 delete[] obj;
00209 delete[] cmatbeg;
00210 delete[] cmatcnt;
00211 delete[] cmatind;
00212 delete[] cmatval;
00213
00215 int* qmatbeg=new int[num_variables];
00216 int* qmatcnt=new int[num_variables];
00217 int* qmatind=new int[num_variables];
00218 double* qmatval=new double[num_variables];
00219
00220 float64_t diag=2.0;
00221
00222 for (int32_t i=0; i<num_variables; i++)
00223 {
00224 if (i<num_dim)
00225
00226 {
00227 qmatbeg[i]=i;
00228 qmatcnt[i]=1;
00229 qmatind[i]=i;
00230 qmatval[i]=diag;
00231 }
00232 else
00233 {
00234
00235 qmatbeg[i]= num_dim;
00236 qmatcnt[i]=0;
00237 qmatind[i]=0;
00238 qmatval[i]=0;
00239 }
00240 }
00241
00242 if (result)
00243 result = CPXcopyquad(env, lp, qmatbeg, qmatcnt, qmatind, qmatval) == 0;
00244
00245 delete[] qmatbeg;
00246 delete[] qmatcnt;
00247 delete[] qmatind;
00248 delete[] qmatval;
00249
00250 if (!result)
00251 SG_ERROR("CPXcopyquad failed.\n");
00252
00253
00254
00255
00256 return result;
00257 }
00258
00259 bool CCplex::setup_lpboost(float64_t C, int32_t num_cols)
00260 {
00261 init(E_LINEAR);
00262 int32_t status = CPXsetintparam (env, CPX_PARAM_LPMETHOD, 1);
00263 if (status)
00264 SG_ERROR( "Failure to select dual lp optimization, error %d.\n", status);
00265
00266 double* obj=new double[num_cols];
00267 double* lb=new double[num_cols];
00268 double* ub=new double[num_cols];
00269
00270 for (int32_t i=0; i<num_cols; i++)
00271 {
00272 obj[i]=-1;
00273 lb[i]=0;
00274 ub[i]=C;
00275 }
00276
00277 status = CPXnewcols(env, lp, num_cols, obj, lb, ub, NULL, NULL);
00278 if ( status )
00279 {
00280 char errmsg[1024];
00281 CPXgeterrorstring (env, status, errmsg);
00282 SG_ERROR( "%s", errmsg);
00283 }
00284 delete[] obj;
00285 delete[] lb;
00286 delete[] ub;
00287 return status==0;
00288 }
00289
00290 bool CCplex::add_lpboost_constraint(
00291 float64_t factor, TSparseEntry<float64_t>* h, int32_t len, int32_t ulen,
00292 CLabels* label)
00293 {
00294 int amatbeg[1];
00295 int amatind[len+1];
00296 double amatval[len+1];
00297 double rhs[1];
00298 char sense[1];
00299
00300 amatbeg[0]=0;
00301 rhs[0]=1;
00302 sense[0]='L';
00303
00304 for (int32_t i=0; i<len; i++)
00305 {
00306 int32_t idx=h[i].feat_index;
00307 float64_t val=factor*h[i].entry;
00308 amatind[i]=idx;
00309 amatval[i]=label->get_label(idx)*val;
00310 }
00311
00312 int32_t status = CPXaddrows (env, lp, 0, 1, len, rhs, sense, amatbeg, amatind, amatval, NULL, NULL);
00313
00314 if ( status )
00315 SG_ERROR( "Failed to add the new row.\n");
00316
00317 return status == 0;
00318 }
00319
00320 bool CCplex::setup_lpm(
00321 float64_t C, CSparseFeatures<float64_t>* x, CLabels* y, bool use_bias)
00322 {
00323 ASSERT(x);
00324 ASSERT(y);
00325 int32_t num_vec=y->get_num_labels();
00326 int32_t num_feat=x->get_num_features();
00327 int64_t nnz=x->get_num_nonzero_entries();
00328
00329
00330 int32_t num_dims=1+2*num_feat+num_vec;
00331 int32_t num_constraints=num_vec;
00332
00333 float64_t* lb=new float64_t[num_dims];
00334 float64_t* ub=new float64_t[num_dims];
00335 float64_t* f=new float64_t[num_dims];
00336 float64_t* b=new float64_t[num_dims];
00337
00338
00339 int64_t amatsize=((int64_t) num_vec)+nnz+nnz+num_vec;
00340
00341 int* amatbeg=new int[num_dims];
00342 int* amatcnt=new int[num_dims];
00343 int* amatind=new int[amatsize];
00344 double* amatval= new double[amatsize];
00345
00346 for (int32_t i=0; i<num_dims; i++)
00347 {
00348 if (i==0)
00349 {
00350 if (use_bias)
00351 {
00352 lb[i]=-CPX_INFBOUND;
00353 ub[i]=+CPX_INFBOUND;
00354 }
00355 else
00356 {
00357 lb[i]=0;
00358 ub[i]=0;
00359 }
00360 f[i]=0;
00361 }
00362 else if (i<2*num_feat+1)
00363 {
00364 lb[i]=0;
00365 ub[i]=CPX_INFBOUND;
00366 f[i]=1;
00367 }
00368 else
00369 {
00370 lb[i]=0;
00371 ub[i]=CPX_INFBOUND;
00372 f[i]=C;
00373 }
00374 }
00375
00376 for (int32_t i=0; i<num_constraints; i++)
00377 b[i]=-1;
00378
00379 char* sense=new char[num_constraints];
00380 memset(sense,'L',sizeof(char)*num_constraints);
00381
00382
00383 int64_t offs=0;
00384
00385
00386 amatbeg[0]=offs;
00387 amatcnt[0]=num_vec;
00388
00389 for (int32_t i=0; i<num_vec; i++)
00390 {
00391 amatind[offs]=i;
00392 amatval[offs]=-y->get_label(i);
00393 offs++;
00394 }
00395
00396
00397 int32_t num_sfeat=0;
00398 int32_t num_svec=0;
00399 TSparse<float64_t>* sfeat= x->get_transposed(num_sfeat, num_svec);
00400 ASSERT(sfeat);
00401
00402 for (int32_t i=0; i<num_svec; i++)
00403 {
00404 amatbeg[i+1]=offs;
00405 amatcnt[i+1]=sfeat[i].num_feat_entries;
00406
00407 for (int32_t j=0; j<sfeat[i].num_feat_entries; j++)
00408 {
00409 int32_t row=sfeat[i].features[j].feat_index;
00410 float64_t val=sfeat[i].features[j].entry;
00411
00412 amatind[offs]=row;
00413 amatval[offs]=-y->get_label(row)*val;
00414 offs++;
00415 }
00416 }
00417
00418 for (int32_t i=0; i<num_svec; i++)
00419 {
00420 amatbeg[num_svec+i+1]=offs;
00421 amatcnt[num_svec+i+1]=sfeat[i].num_feat_entries;
00422
00423 for (int32_t j=0; j<sfeat[i].num_feat_entries; j++)
00424 {
00425 int32_t row=sfeat[i].features[j].feat_index;
00426 float64_t val=sfeat[i].features[j].entry;
00427
00428 amatind[offs]=row;
00429 amatval[offs]=y->get_label(row)*val;
00430 offs++;
00431 }
00432 }
00433
00434 x->clean_tsparse(sfeat, num_svec);
00435
00436
00437 for (int32_t i=0; i<num_vec; i++)
00438 {
00439 amatbeg[1+2*num_feat+i]=offs;
00440 amatcnt[1+2*num_feat+i]=1;
00441 amatind[offs]=i;
00442 amatval[offs]=-1;
00443 offs++;
00444 }
00445
00446 int32_t status = CPXsetintparam (env, CPX_PARAM_LPMETHOD, 1);
00447 if (status)
00448 SG_ERROR( "Failure to select barrier optimization, error %d.\n", status);
00449 CPXsetintparam (env, CPX_PARAM_SCRIND, CPX_ON);
00450
00451 bool result = CPXcopylp(env, lp, num_dims, num_constraints, CPX_MIN,
00452 f, b, sense, amatbeg, amatcnt, amatind, amatval, lb, ub, NULL) == 0;
00453
00454
00455 delete[] amatval;
00456 delete[] amatcnt;
00457 delete[] amatind;
00458 delete[] amatbeg;
00459 delete[] b;
00460 delete[] f;
00461 delete[] ub;
00462 delete[] lb;
00463
00464 return result;
00465 }
00466
00467 bool CCplex::cleanup()
00468 {
00469 bool result=false;
00470
00471 if (lp)
00472 {
00473 int32_t status = CPXfreeprob(env, &lp);
00474 lp = NULL;
00475 lp_initialized = false;
00476
00477 if (status)
00478 SG_WARNING( "CPXfreeprob failed, error code %d.\n", status);
00479 else
00480 result = true;
00481 }
00482
00483 if (env)
00484 {
00485 int32_t status = CPXcloseCPLEX (&env);
00486 env=NULL;
00487
00488 if (status)
00489 {
00490 char errmsg[1024];
00491 SG_WARNING( "Could not close CPLEX environment.\n");
00492 CPXgeterrorstring (env, status, errmsg);
00493 SG_WARNING( "%s", errmsg);
00494 }
00495 else
00496 result = true;
00497 }
00498 return result;
00499 }
00500
00501 bool CCplex::dense_to_cplex_sparse(
00502 float64_t* H, int32_t rows, int32_t cols, int* &qmatbeg, int* &qmatcnt,
00503 int* &qmatind, double* &qmatval)
00504 {
00505 qmatbeg=new int[cols];
00506 qmatcnt=new int[cols];
00507 qmatind=new int[cols*rows];
00508 qmatval = H;
00509
00510 if (!(qmatbeg && qmatcnt && qmatind))
00511 {
00512 delete[] qmatbeg;
00513 delete[] qmatcnt;
00514 delete[] qmatind;
00515 return false;
00516 }
00517
00518 for (int32_t i=0; i<cols; i++)
00519 {
00520 qmatcnt[i]=rows;
00521 qmatbeg[i]=i*rows;
00522 for (int32_t j=0; j<rows; j++)
00523 qmatind[i*rows+j]=j;
00524 }
00525
00526 return true;
00527 }
00528
00529 bool CCplex::setup_lp(
00530 float64_t* objective, float64_t* constraints_mat, int32_t rows,
00531 int32_t cols, float64_t* rhs, float64_t* lb, float64_t* ub)
00532 {
00533 bool result=false;
00534
00535 int* qmatbeg=NULL;
00536 int* qmatcnt=NULL;
00537 int* qmatind=NULL;
00538 double* qmatval=NULL;
00539
00540 char* sense = NULL;
00541
00542 if (constraints_mat==0)
00543 {
00544 ASSERT(rows==0);
00545 rows=1;
00546 float64_t dummy=0;
00547 rhs=&dummy;
00548 sense=new char[rows];
00549 memset(sense,'L',sizeof(char)*rows);
00550 constraints_mat=new float64_t[cols];
00551 memset(constraints_mat, 0, sizeof(float64_t)*cols);
00552
00553 result=dense_to_cplex_sparse(constraints_mat, 0, cols, qmatbeg, qmatcnt, qmatind, qmatval);
00554 ASSERT(result);
00555 result = CPXcopylp(env, lp, cols, rows, CPX_MIN,
00556 objective, rhs, sense, qmatbeg, qmatcnt, qmatind, qmatval, lb, ub, NULL) == 0;
00557 delete[] constraints_mat;
00558 }
00559 else
00560 {
00561 sense=new char[rows];
00562 memset(sense,'L',sizeof(char)*rows);
00563 result=dense_to_cplex_sparse(constraints_mat, rows, cols, qmatbeg, qmatcnt, qmatind, qmatval);
00564 result = CPXcopylp(env, lp, cols, rows, CPX_MIN,
00565 objective, rhs, sense, qmatbeg, qmatcnt, qmatind, qmatval, lb, ub, NULL) == 0;
00566 }
00567
00568 delete[] sense;
00569 delete[] qmatbeg;
00570 delete[] qmatcnt;
00571 delete[] qmatind;
00572
00573 if (!result)
00574 SG_WARNING("CPXcopylp failed.\n");
00575
00576 return result;
00577 }
00578
00579 bool CCplex::setup_qp(float64_t* H, int32_t dim)
00580 {
00581 int* qmatbeg=NULL;
00582 int* qmatcnt=NULL;
00583 int* qmatind=NULL;
00584 double* qmatval=NULL;
00585 bool result=dense_to_cplex_sparse(H, dim, dim, qmatbeg, qmatcnt, qmatind, qmatval);
00586 if (result)
00587 result = CPXcopyquad(env, lp, qmatbeg, qmatcnt, qmatind, qmatval) == 0;
00588
00589 delete[] qmatbeg;
00590 delete[] qmatcnt;
00591 delete[] qmatind;
00592
00593 if (!result)
00594 SG_WARNING("CPXcopyquad failed.\n");
00595
00596 return result;
00597 }
00598
00599 bool CCplex::optimize(float64_t* sol, float64_t* lambda)
00600 {
00601 int solnstat;
00602 double objval;
00603 int status=1;
00604
00605 if (problem_type==E_QP)
00606 status = CPXqpopt (env, lp);
00607 else if (problem_type == E_LINEAR)
00608 status = CPXlpopt (env, lp);
00609
00610 if (status)
00611 SG_WARNING( "Failed to optimize QP.\n");
00612
00613 status = CPXsolution (env, lp, &solnstat, &objval, sol, lambda, NULL, NULL);
00614
00615
00616
00617 return (status==0);
00618 }
00619 #endif