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, int32_t timeout)
00032 {
00033     problem_type=typ;
00034 
00035     while (env==NULL)
00036     {
00037         int status = 1; /* for calling external lib */
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(
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; //no more than 10mil. elements
00091     bool result=false;
00092     int32_t 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     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) //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     int32_t offs=0;
00124     for (int32_t 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             int32_t idx=idx_bound[i-num_dim-num_zero];
00148             int32_t vlen=0;
00149             bool vfree=false;
00150             //SG_PRINT("idx=%d\n", idx);
00151             TSparseEntry<float64_t>* 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             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                     //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     float64_t diag=2.0;
00221 
00222     for (int32_t 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(float64_t C, int32_t num_cols)
00260 {
00261     init(E_LINEAR);
00262     int32_t 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 (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]; /* for calling external lib */
00295     int amatind[len+1]; /* for calling external lib */
00296     double amatval[len+1]; /* for calling external lib */
00297     double rhs[1]; /* for calling external lib */
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     //number of variables: b,w+,w-,xi concatenated
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     //number of non zero entries in A (b,w+,w-,xi)
00339     int64_t amatsize=((int64_t) num_vec)+nnz+nnz+num_vec; 
00340 
00341     int* amatbeg=new int[num_dims]; /* for calling external lib */
00342     int* amatcnt=new int[num_dims]; /* for calling external lib */
00343     int* amatind=new int[amatsize]; /* for calling external lib */
00344     double* amatval= new double[amatsize]; /* for calling external lib */
00345 
00346     for (int32_t i=0; i<num_dims; i++)
00347     {
00348         if (i==0) //b
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) //w+,w-
00363         {
00364             lb[i]=0;
00365             ub[i]=CPX_INFBOUND;
00366             f[i]=1;
00367         }
00368         else //xi
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     //construct A
00383     int64_t offs=0;
00384 
00385     //b part of A
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     //w+ and w- part of A
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     //xi part of A
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); //barrier
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; /* for calling external lib */
00602     double   objval;
00603     int status=1; /* for calling external lib */
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     //SG_PRINT("obj:%f\n", objval);
00616 
00617     return (status==0);
00618 }
00619 #endif

SHOGUN Machine Learning Toolbox - Documentation