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-2009 Soeren Sonnenburg
00008  * Copyright (C) 1999-2009 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 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, CLabels* 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=new float64_t[num_variables];
00101     float64_t* ub=new float64_t[num_variables];
00102     float64_t* obj=new float64_t[num_variables];
00103     char* sense = new char[num_dim];
00104     int* cmatbeg=new int[num_variables];
00105     int* cmatcnt=new int[num_variables];
00106     int* cmatind=new int[cmatsize];
00107     double* cmatval=new 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             TSparseEntry<float64_t>* vec=features->get_sparse_feature_vector(idx, vlen, vfree);
00154             //SG_PRINT("vlen=%d\n", vlen);
00155 
00156             cmatbeg[i]=offs;
00157             cmatcnt[i]=vlen;
00158 
00159             float64_t val= -C*labels->get_label(idx);
00160 
00161             if (vlen>0)
00162             {
00163                 for (int32_t j=0; j<vlen; j++)
00164                 {
00165                     cmatind[offs]=vec[j].feat_index;
00166                     cmatval[offs]=-val*vec[j].entry;
00167                     offs++;
00168                     ASSERT(offs<cmatsize);
00169                     //SG_PRINT("vec[%d]=%10.10f\n", j, vec[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, vfree);
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     delete[] sense;
00208     delete[] lb;
00209     delete[] ub;
00210     delete[] obj;
00211     delete[] cmatbeg;
00212     delete[] cmatcnt;
00213     delete[] cmatind;
00214     delete[] cmatval;
00215 
00217     int* qmatbeg=new int[num_variables];
00218     int* qmatcnt=new int[num_variables];
00219     int* qmatind=new int[num_variables];
00220     double* qmatval=new 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     delete[] qmatbeg;
00248     delete[] qmatcnt;
00249     delete[] qmatind;
00250     delete[] 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=new double[num_cols];
00269     double* lb=new double[num_cols];
00270     double* ub=new 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     delete[] obj;
00287     delete[] lb;
00288     delete[] ub;
00289     return status==0;
00290 }
00291 
00292 bool CCplex::add_lpboost_constraint(
00293     float64_t factor, TSparseEntry<float64_t>* h, int32_t len, int32_t ulen,
00294     CLabels* 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_label(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, CLabels* 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=new float64_t[num_dims];
00336     float64_t* ub=new float64_t[num_dims];
00337     float64_t* f=new float64_t[num_dims];
00338     float64_t* b=new 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=new int[num_dims]; /* for calling external lib */
00344     int* amatcnt=new int[num_dims]; /* for calling external lib */
00345     int* amatind=new int[amatsize]; /* for calling external lib */
00346     double* amatval= new 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=new 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_label(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     TSparse<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_label(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_label(row)*val;
00432             offs++;
00433         }
00434     }
00435 
00436     x->clean_tsparse(sfeat, num_svec);
00437 
00438     //xi part of A
00439     for (int32_t i=0; i<num_vec; i++)
00440     {
00441         amatbeg[1+2*num_feat+i]=offs;
00442         amatcnt[1+2*num_feat+i]=1;
00443         amatind[offs]=i;
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     delete[] amatval;
00458     delete[] amatcnt;
00459     delete[] amatind;
00460     delete[] amatbeg;
00461     delete[] b;
00462     delete[] f;
00463     delete[] ub;
00464     delete[] 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=new int[cols];
00508     qmatcnt=new int[cols];
00509     qmatind=new int[cols*rows];
00510     qmatval = H;
00511 
00512     if (!(qmatbeg && qmatcnt && qmatind))
00513     {
00514         delete[] qmatbeg;
00515         delete[] qmatcnt;
00516         delete[] 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=new char[rows];
00551         memset(sense,'L',sizeof(char)*rows);
00552         constraints_mat=new 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         delete[] constraints_mat;
00560     }
00561     else
00562     {
00563         sense=new 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     delete[] sense;
00571     delete[] qmatbeg;
00572     delete[] qmatcnt;
00573     delete[] 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     delete[] qmatbeg;
00592     delete[] qmatcnt;
00593     delete[] 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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation