Cplex.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) 1999-2008 Soeren Sonnenburg
00008  * Copyright (C) 1999-2008 Fraunhofer Institute FIRST and Max-Planck-Society
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, INT 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             /* Turn on output to the screen */
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);  /* Problem is minimization */
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(DREAL C, CLabels* labels, CSparseFeatures<DREAL>* features, INT* idx_bound, INT num_bound,
00086         INT* w_zero, INT num_zero,
00087         DREAL* vee, INT num_dim,
00088         bool use_bias)
00089 {
00090     const int cmatsize=10*1024*1024; //no more than 10mil. elements
00091     bool result=false;
00092     INT num_variables=num_dim + num_bound+num_zero; // xi, beta
00093 
00094     ASSERT(num_dim>0);
00095     ASSERT(num_dim>0);
00096 
00097     // setup LP part
00098     DREAL* lb=new DREAL[num_variables];
00099     DREAL* ub=new DREAL[num_variables];
00100     DREAL* obj=new DREAL[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 (INT i=0; i<num_variables; i++)
00108     {
00109         obj[i]=0;
00110 
00111         if (i<num_dim) //xi part
00112         {
00113             lb[i]=-CPX_INFBOUND;
00114             ub[i]=+CPX_INFBOUND;
00115         }
00116         else //beta part
00117         {
00118             lb[i]=0.0;
00119             ub[i]=1.0;
00120         }
00121     }
00122 
00123     INT offs=0;
00124     for (INT i=0; i<num_variables; i++)
00125     {
00126         if (i<num_dim) //sum_xi
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) // Z_w*beta_w
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 // Z_x*beta_x
00146         {
00147             INT idx=idx_bound[i-num_dim-num_zero];
00148             INT vlen=0;
00149             bool vfree=false;
00150             //SG_PRINT("idx=%d\n", idx);
00151             TSparseEntry<DREAL>* vec=features->get_sparse_feature_vector(idx, vlen, vfree);
00152             //SG_PRINT("vlen=%d\n", vlen);
00153 
00154             cmatbeg[i]=offs;
00155             cmatcnt[i]=vlen;
00156 
00157             DREAL val= -C*labels->get_label(idx);
00158 
00159             if (vlen>0)
00160             {
00161                 for (INT 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                     //SG_PRINT("vec[%d]=%10.10f\n", j, vec[j].entry);
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     //write_problem("problem.lp");
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     DREAL diag=2.0;
00221 
00222     for (INT i=0; i<num_variables; i++)
00223     {
00224         if (i<num_dim) //|| (!use_bias && i<num_dim)) //xi
00225         //if ((i<num_dim-1) || (!use_bias && i<num_dim)) //xi
00226         {
00227             qmatbeg[i]=i;
00228             qmatcnt[i]=1;
00229             qmatind[i]=i;
00230             qmatval[i]=diag;
00231         }
00232         else
00233         {
00234             //qmatbeg[i]= (use_bias) ? (num_dim-1) : (num_dim);
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     //write_problem("problem.lp");
00254     //write_Q("problem.qp");
00255 
00256     return result;
00257 }
00258 
00259 bool CCplex::setup_lpboost(DREAL C, INT num_cols)
00260 {
00261     init(E_LINEAR);
00262     INT status = CPXsetintparam (env, CPX_PARAM_LPMETHOD, 1); //primal simplex
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 (INT 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(DREAL factor, TSparseEntry<DREAL>* h, INT len, INT ulen, CLabels* label)
00291 {
00292     int amatbeg[1];
00293     int amatind[len+1];
00294     double amatval[len+1];
00295     double rhs[1];
00296     char sense[1];
00297 
00298     amatbeg[0]=0;
00299     rhs[0]=1;
00300     sense[0]='L';
00301 
00302     for (INT i=0; i<len; i++)
00303     {
00304         INT idx=h[i].feat_index;
00305         DREAL val=factor*h[i].entry;
00306         amatind[i]=idx;
00307         amatval[i]=label->get_label(idx)*val;
00308     }
00309 
00310     INT status = CPXaddrows (env, lp, 0, 1, len, rhs, sense, amatbeg, amatind, amatval, NULL, NULL);
00311 
00312     if ( status ) 
00313         SG_ERROR( "Failed to add the new row.\n");
00314 
00315     return status == 0;
00316 }
00317 
00318 bool CCplex::setup_lpm(DREAL C, CSparseFeatures<DREAL>* x, CLabels* y, bool use_bias)
00319 {
00320     ASSERT(x);
00321     ASSERT(y);
00322     INT num_vec=y->get_num_labels();
00323     INT num_feat=x->get_num_features();
00324     LONG nnz=x->get_num_nonzero_entries();
00325 
00326     //number of variables: b,w+,w-,xi concatenated
00327     INT num_dims=1+2*num_feat+num_vec;
00328     INT num_constraints=num_vec; 
00329 
00330     DREAL* lb=new DREAL[num_dims];
00331     DREAL* ub=new DREAL[num_dims];
00332     DREAL* f=new DREAL[num_dims];
00333     DREAL* b=new DREAL[num_dims];
00334 
00335     //number of non zero entries in A (b,w+,w-,xi)
00336     LONG amatsize=((LONG) num_vec)+nnz+nnz+num_vec; 
00337 
00338     int* amatbeg=new int[num_dims];
00339     int* amatcnt=new int[num_dims];
00340     int* amatind=new int[amatsize];
00341     double* amatval= new double[amatsize];
00342 
00343     for (INT i=0; i<num_dims; i++)
00344     {
00345         if (i==0) //b
00346         {
00347             if (use_bias)
00348             {
00349                 lb[i]=-CPX_INFBOUND;
00350                 ub[i]=+CPX_INFBOUND;
00351             }
00352             else
00353             {
00354                 lb[i]=0;
00355                 ub[i]=0;
00356             }
00357             f[i]=0;
00358         }
00359         else if (i<2*num_feat+1) //w+,w-
00360         {
00361             lb[i]=0;
00362             ub[i]=CPX_INFBOUND;
00363             f[i]=1;
00364         }
00365         else //xi
00366         {
00367             lb[i]=0;
00368             ub[i]=CPX_INFBOUND;
00369             f[i]=C;
00370         }
00371     }
00372 
00373     for (INT i=0; i<num_constraints; i++)
00374         b[i]=-1;
00375 
00376     char* sense=new char[num_constraints];
00377     memset(sense,'L',sizeof(char)*num_constraints);
00378 
00379     //construct A
00380     LONG offs=0;
00381 
00382     //b part of A
00383     amatbeg[0]=offs;
00384     amatcnt[0]=num_vec;
00385 
00386     for (INT i=0; i<num_vec; i++)
00387     {
00388         amatind[offs]=i;
00389         amatval[offs]=-y->get_label(i);
00390         offs++;
00391     }
00392 
00393     //w+ and w- part of A
00394     INT num_sfeat=0;
00395     INT num_svec=0;
00396     TSparse<DREAL>* sfeat= x->get_transposed(num_sfeat, num_svec);
00397     ASSERT(sfeat);
00398 
00399     for (INT i=0; i<num_svec; i++)
00400     {
00401         amatbeg[i+1]=offs;
00402         amatcnt[i+1]=sfeat[i].num_feat_entries;
00403 
00404         for (INT j=0; j<sfeat[i].num_feat_entries; j++)
00405         {
00406             INT row=sfeat[i].features[j].feat_index;
00407             DREAL val=sfeat[i].features[j].entry;
00408 
00409             amatind[offs]=row;
00410             amatval[offs]=-y->get_label(row)*val;
00411             offs++;
00412         }
00413     }
00414 
00415     for (INT i=0; i<num_svec; i++)
00416     {
00417         amatbeg[num_svec+i+1]=offs;
00418         amatcnt[num_svec+i+1]=sfeat[i].num_feat_entries;
00419 
00420         for (INT j=0; j<sfeat[i].num_feat_entries; j++)
00421         {
00422             INT row=sfeat[i].features[j].feat_index;
00423             DREAL val=sfeat[i].features[j].entry;
00424 
00425             amatind[offs]=row;
00426             amatval[offs]=y->get_label(row)*val;
00427             offs++;
00428         }
00429     }
00430 
00431     x->clean_tsparse(sfeat, num_svec);
00432 
00433     //xi part of A
00434     for (INT i=0; i<num_vec; i++)
00435     {
00436         amatbeg[1+2*num_feat+i]=offs;
00437         amatcnt[1+2*num_feat+i]=1;
00438         amatind[offs]=i;
00439         amatval[offs]=-1;
00440         offs++;
00441     }
00442 
00443     INT status = CPXsetintparam (env, CPX_PARAM_LPMETHOD, 1); //barrier
00444     if (status)
00445         SG_ERROR( "Failure to select barrier optimization, error %d.\n", status);
00446     CPXsetintparam (env, CPX_PARAM_SCRIND, CPX_ON);
00447 
00448     bool result = CPXcopylp(env, lp, num_dims, num_constraints, CPX_MIN, 
00449             f, b, sense, amatbeg, amatcnt, amatind, amatval, lb, ub, NULL) == 0;
00450     
00451 
00452     delete[] amatval;
00453     delete[] amatcnt;
00454     delete[] amatind;
00455     delete[] amatbeg;
00456     delete[] b;
00457     delete[] f;
00458     delete[] ub;
00459     delete[] lb;
00460 
00461     return result;
00462 }
00463 
00464 bool CCplex::cleanup()
00465 {
00466     bool result=false;
00467 
00468     if (lp)
00469     {
00470         INT status = CPXfreeprob(env, &lp);
00471         lp = NULL;
00472         lp_initialized = false;
00473 
00474         if (status)
00475             SG_WARNING( "CPXfreeprob failed, error code %d.\n", status);
00476         else
00477             result = true;
00478     }
00479 
00480     if (env)
00481     {
00482         INT status = CPXcloseCPLEX (&env);
00483         env=NULL;
00484         
00485         if (status)
00486         {
00487             char  errmsg[1024];
00488             SG_WARNING( "Could not close CPLEX environment.\n");
00489             CPXgeterrorstring (env, status, errmsg);
00490             SG_WARNING( "%s", errmsg);
00491         }
00492         else
00493             result = true;
00494     }
00495     return result;
00496 }
00497 
00498 bool CCplex::dense_to_cplex_sparse(DREAL* H, INT rows, INT cols, int* &qmatbeg, int* &qmatcnt, int* &qmatind, double* &qmatval)
00499 {
00500     qmatbeg=new int[cols];
00501     qmatcnt=new int[cols];
00502     qmatind=new int[cols*rows];
00503     qmatval = H;
00504 
00505     if (!(qmatbeg && qmatcnt && qmatind))
00506     {
00507         delete[] qmatbeg;
00508         delete[] qmatcnt;
00509         delete[] qmatind;
00510         return false;
00511     }
00512 
00513     for (INT i=0; i<cols; i++)
00514     {
00515         qmatcnt[i]=rows;
00516         qmatbeg[i]=i*rows;
00517         for (INT j=0; j<rows; j++)
00518             qmatind[i*rows+j]=j;
00519     }
00520 
00521     return true;
00522 }
00523 
00524 bool CCplex::setup_lp(DREAL* objective, DREAL* constraints_mat, INT rows, INT cols, DREAL* rhs, DREAL* lb, DREAL* ub)
00525 {
00526     bool result=false;
00527 
00528     int* qmatbeg=NULL;
00529     int* qmatcnt=NULL;
00530     int* qmatind=NULL;
00531     double* qmatval=NULL;
00532 
00533     char* sense = NULL;
00534 
00535     if (constraints_mat==0)
00536     {
00537         ASSERT(rows==0);
00538         rows=1;
00539         DREAL dummy=0;
00540         rhs=&dummy;
00541         sense=new char[rows];
00542         memset(sense,'L',sizeof(char)*rows);
00543         constraints_mat=new DREAL[cols];
00544         memset(constraints_mat, 0, sizeof(DREAL)*cols);
00545 
00546         result=dense_to_cplex_sparse(constraints_mat, 0, cols, qmatbeg, qmatcnt, qmatind, qmatval);
00547         ASSERT(result);
00548         result = CPXcopylp(env, lp, cols, rows, CPX_MIN, 
00549                 objective, rhs, sense, qmatbeg, qmatcnt, qmatind, qmatval, lb, ub, NULL) == 0;
00550         delete[] constraints_mat;
00551     }
00552     else
00553     {
00554         sense=new char[rows];
00555         memset(sense,'L',sizeof(char)*rows);
00556         result=dense_to_cplex_sparse(constraints_mat, rows, cols, qmatbeg, qmatcnt, qmatind, qmatval);
00557         result = CPXcopylp(env, lp, cols, rows, CPX_MIN, 
00558                 objective, rhs, sense, qmatbeg, qmatcnt, qmatind, qmatval, lb, ub, NULL) == 0;
00559     }
00560 
00561     delete[] sense;
00562     delete[] qmatbeg;
00563     delete[] qmatcnt;
00564     delete[] qmatind;
00565 
00566     if (!result)
00567         SG_WARNING("CPXcopylp failed.\n");
00568 
00569     return result;
00570 }
00571 
00572 bool CCplex::setup_qp(DREAL* H, INT dim)
00573 {
00574     int* qmatbeg=NULL;
00575     int* qmatcnt=NULL;
00576     int* qmatind=NULL;
00577     double* qmatval=NULL;
00578     bool result=dense_to_cplex_sparse(H, dim, dim, qmatbeg, qmatcnt, qmatind, qmatval);
00579     if (result)
00580         result = CPXcopyquad(env, lp, qmatbeg, qmatcnt, qmatind, qmatval) == 0;
00581 
00582     delete[] qmatbeg;
00583     delete[] qmatcnt;
00584     delete[] qmatind;
00585 
00586     if (!result)
00587         SG_WARNING("CPXcopyquad failed.\n");
00588 
00589     return result;
00590 }
00591 
00592 bool CCplex::optimize(DREAL* sol, DREAL* lambda)
00593 {
00594     int      solnstat;
00595     double   objval;
00596     int status=1;
00597 
00598     if (problem_type==E_QP)
00599         status = CPXqpopt (env, lp);
00600     else if (problem_type == E_LINEAR)
00601         status = CPXlpopt (env, lp);
00602 
00603     if (status)
00604         SG_WARNING( "Failed to optimize QP.\n");
00605 
00606     status = CPXsolution (env, lp, &solnstat, &objval, sol, lambda, NULL, NULL);
00607 
00608     //SG_PRINT("obj:%f\n", objval);
00609 
00610     return (status==0);
00611 }
00612 #endif

SHOGUN Machine Learning Toolbox - Documentation