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