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) 1999-2009 Soeren Sonnenburg 00008 * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society 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; /* for calling external lib */ 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 /* Turn on output to the screen */ 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); /* Problem is minimization */ 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; //no more than 10mil. elements 00093 bool result=false; 00094 int32_t num_variables=num_dim + num_bound+num_zero; // xi, beta 00095 00096 ASSERT(num_dim>0) 00097 ASSERT(num_dim>0) 00098 00099 // setup LP part 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) //xi part 00114 { 00115 lb[i]=-CPX_INFBOUND; 00116 ub[i]=+CPX_INFBOUND; 00117 } 00118 else //beta part 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) //sum_xi 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) // Z_w*beta_w 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 // Z_x*beta_x 00148 { 00149 int32_t idx=idx_bound[i-num_dim-num_zero]; 00150 int32_t vlen=0; 00151 bool vfree=false; 00152 //SG_PRINT("idx=%d\n", idx) 00153 SGSparseVector<float64_t> vec=features->get_sparse_feature_vector(idx); 00154 //SG_PRINT("vlen=%d\n", vlen) 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 //SG_PRINT("vec[%d]=%10.10f\n", j, vec.features[j].entry) 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 //write_problem("problem.lp"); 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) //|| (!use_bias && i<num_dim)) //xi 00227 //if ((i<num_dim-1) || (!use_bias && i<num_dim)) //xi 00228 { 00229 qmatbeg[i]=i; 00230 qmatcnt[i]=1; 00231 qmatind[i]=i; 00232 qmatval[i]=diag; 00233 } 00234 else 00235 { 00236 //qmatbeg[i]= (use_bias) ? (num_dim-1) : (num_dim); 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 //write_problem("problem.lp"); 00256 //write_Q("problem.qp"); 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); //primal simplex 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]; /* for calling external lib */ 00297 int amatind[len+1]; /* for calling external lib */ 00298 double amatval[len+1]; /* for calling external lib */ 00299 double rhs[1]; /* for calling external lib */ 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 //number of variables: b,w+,w-,xi concatenated 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 //number of non zero entries in A (b,w+,w-,xi) 00341 int64_t amatsize=((int64_t) num_vec)+nnz+nnz+num_vec; 00342 00343 int* amatbeg=SG_MALLOC(int, num_dims); /* for calling external lib */ 00344 int* amatcnt=SG_MALLOC(int, num_dims); /* for calling external lib */ 00345 int* amatind=SG_MALLOC(int, amatsize); /* for calling external lib */ 00346 double* amatval= SG_MALLOC(double, amatsize); /* for calling external lib */ 00347 00348 for (int32_t i=0; i<num_dims; i++) 00349 { 00350 if (i==0) //b 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) //w+,w- 00365 { 00366 lb[i]=0; 00367 ub[i]=CPX_INFBOUND; 00368 f[i]=1; 00369 } 00370 else //xi 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 //construct A 00385 int64_t offs=0; 00386 00387 //b part of A 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 //w+ and w- part of A 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 //xi part of A 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); //barrier 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; /* for calling external lib */ 00604 double objval; 00605 int status=1; /* for calling external lib */ 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 //SG_PRINT("obj:%f\n", objval) 00618 00619 return (status==0); 00620 } 00621 #endif