#include <CglLiftAndProject.hpp>
Inheritance diagram for CglLiftAndProject:
Public Member Functions | |
Generate Cuts | |
virtual void | generateCuts (const OsiSolverInterface &si, OsiCuts &cs) const |
double | getBeta () const |
void | setBeta (int oneOrMinusOne) |
Constructors and destructors | |
CglLiftAndProject () | |
Default constructor. | |
CglLiftAndProject (const CglLiftAndProject &) | |
Copy constructor. | |
CglLiftAndProject & | operator= (const CglLiftAndProject &rhs) |
Assignment operator. | |
virtual | ~CglLiftAndProject () |
Destructor. | |
Private Attributes | |
Private member data | |
double | beta_ |
The normalization is beta_=1 or beta_=-1. | |
double | epsilon_ |
epsilon | |
double | onetol_ |
1-epsilon | |
Friends | |
void | CglLiftAndProjectUnitTest (const OsiSolverInterface *siP, const std::string mpdDir) |
Definition at line 11 of file CglLiftAndProject.hpp.
|
Generate lift-and-project cuts for the model of the solver interface, si. Insert the generated cuts into OsiCut, cs. Implements CglCutGenerator. Definition at line 23 of file CglLiftAndProject.cpp. References beta_.
00025 { 00026 // Assumes the mixed 0-1 problem 00027 // 00028 // min {cx: <Atilde,x> >= btilde} 00029 // 00030 // is in cannonical form with all bounds, 00031 // including x_t>=0, -x_t>=-1 for x_t binary, 00032 // explicitly stated in the constraint matrix. 00033 // See ~/COIN/Examples/Cgl2/cgl2.cpp 00034 // for a general purpose "convert" function. 00035 00036 // Reference [BCC]: Balas, Ceria, and Corneujols, 00037 // "A lift-and-project cutting plane algorithm 00038 // for mixed 0-1 program", Math Prog 58, (1993) 00039 // 295-324. 00040 00041 // This implementation uses Normalization 1. 00042 00043 // Given cannonical problem and 00044 // the lp-relaxation solution, x, 00045 // the LAP cut generator attempts to construct 00046 // a cut for every x_j such that 0<x_j<1 00047 // [BCC:307] 00048 00049 00050 // x_j is the strictly fractional binary variable 00051 // the cut is generated from 00052 int j = 0; 00053 00054 // Get basic problem information 00055 // let Atilde be an m by n matrix 00056 const int m = si.getNumRows(); 00057 const int n = si.getNumCols(); 00058 const double * x = si.getColSolution(); 00059 00060 // Remember - Atildes may have gaps.. 00061 const CoinPackedMatrix * Atilde = si.getMatrixByRow(); 00062 const double * AtildeElements = Atilde->getElements(); 00063 const int * AtildeIndices = Atilde->getIndices(); 00064 const int * AtildeStarts = Atilde->getVectorStarts(); 00065 const int * AtildeLengths = Atilde->getVectorLengths(); 00066 const int AtildeFullSize = AtildeStarts[m]; 00067 const double * btilde = si.getRowLower(); 00068 00069 // Set up memory for system (10) [BCC:307] 00070 // (the problem over the norm intersected 00071 // with the polar cone) 00072 // 00073 // min <<x^T,Atilde^T>,u> + x_ju_0 00074 // s.t. 00075 // <B,w> = (0,...,0,beta_,beta)^T 00076 // w is nonneg for all but the 00077 // last two entries, which are free. 00078 // where 00079 // w = (u,v,v_0,u_0)in BCC notation 00080 // u and v are m-vectors; u,v >=0 00081 // v_0 and u_0 are free-scalars, and 00082 // 00083 // B = Atilde^T -Atilde^T -e_j e_j 00084 // btilde^T e_0^T 0 0 00085 // e_0^T btilde^T 1 0 00086 00087 // ^T indicates Transpose 00088 // e_0 is a (AtildeNCols x 1) vector of all zeros 00089 // e_j is e_0 with a 1 in the jth position 00090 00091 // Storing B in column order. B is a (n+2 x 2m+2) matrix 00092 // But need to allow for possible gaps in Atilde. 00093 // At each iteration, only need to change 2 cols and objfunc 00094 // Sane design of OsiSolverInterface does not permit mucking 00095 // with matrix. 00096 // Because we must delete and add cols to alter matrix, 00097 // and we can only add columns on the end of the matrix 00098 // put the v_0 and u_0 columns on the end. 00099 // rather than as described in [BCC] 00100 00101 // Initially allocating B with space for v_0 and u_O cols 00102 // but not populating, for efficiency. 00103 00104 // B without u_0 and v_0 is a (n+2 x 2m) size matrix. 00105 00106 int twoM = 2*m; 00107 int BNumRows = n+2; 00108 int BNumCols = twoM+2; 00109 int BFullSize = 2*AtildeFullSize+twoM+3; 00110 double * BElements = new double[BFullSize]; 00111 int * BIndices = new int[BFullSize]; 00112 int * BStarts = new int[BNumCols+1]; 00113 int * BLengths = new int[BNumCols]; 00114 00115 00116 int i, ij, k=0; 00117 int nPlus1=n+1; 00118 int offset = AtildeStarts[m]+m; 00119 for (i=0; i<m; i++){ 00120 for (ij=AtildeStarts[i];ij<AtildeStarts[i]+AtildeLengths[i];ij++){ 00121 BElements[k]=AtildeElements[ij]; 00122 BElements[k+offset]=-AtildeElements[ij]; 00123 BIndices[k]= AtildeIndices[ij]; 00124 BIndices[k+offset]= AtildeIndices[ij]; 00125 00126 k++; 00127 } 00128 BElements[k]=btilde[i]; 00129 BElements[k+offset]=btilde[i]; 00130 BIndices[k]=n; 00131 BIndices[k+offset]=nPlus1; 00132 BStarts[i]= AtildeStarts[i]+i; 00133 BStarts[i+m]=offset+BStarts[i];// = AtildeStarts[m]+m+AtildeStarts[i]+i 00134 BLengths[i]= AtildeLengths[i]+1; 00135 BLengths[i+m]= AtildeLengths[i]+1; 00136 k++; 00137 } 00138 00139 BStarts[twoM]=BStarts[twoM-1]+BLengths[twoM-1]; 00140 00141 // Cols that will be deleted each iteration 00142 int BNumColsLessOne=BNumCols-1; 00143 int BNumColsLessTwo=BNumCols-2; 00144 const int delCols[2] = {BNumColsLessOne, BNumColsLessTwo}; 00145 00146 // Set lower bound on u and v 00147 // u_0, v_0 will be reset as free 00148 const double solverINFINITY = si.getInfinity(); 00149 double * BColLowers = new double[BNumCols]; 00150 double * BColUppers = new double[BNumCols]; 00151 CoinFillN(BColLowers,BNumCols,0.0); 00152 CoinFillN(BColUppers,BNumCols,solverINFINITY); 00153 00154 // Set row lowers and uppers. 00155 // The rhs is zero, for but the last two rows. 00156 // For these the rhs is beta_ 00157 double * BRowLowers = new double[BNumRows]; 00158 double * BRowUppers = new double[BNumRows]; 00159 CoinFillN(BRowLowers,BNumRows,0.0); 00160 CoinFillN(BRowUppers,BNumRows,0.0); 00161 BRowLowers[BNumRows-2]=beta_; 00162 BRowUppers[BNumRows-2]=beta_; 00163 BRowLowers[BNumRows-1]=beta_; 00164 BRowUppers[BNumRows-1]=beta_; 00165 00166 00167 // Calculate base objective <<x^T,Atilde^T>,u> 00168 // Note: at each iteration coefficient u_0 00169 // changes to <x^T,e_j> 00170 // w=(u,v,beta,v_0,u_0) size 2m+3 00171 // So, BOjective[2m+2]=x[j] 00172 double * BObjective= new double[BNumCols]; 00173 double * Atildex = new double[m]; 00174 CoinFillN(BObjective,BNumCols,0.0); 00175 Atilde->times(x,Atildex); // Atildex is size m, x is size n 00176 CoinDisjointCopyN(Atildex,m,BObjective); 00177 00178 // Number of cols and size of Elements vector 00179 // in B without the v_0 and u_0 cols 00180 int BFullSizeLessThree = BFullSize-3; 00181 00182 // Load B matrix into a column orders CoinPackedMatrix 00183 CoinPackedMatrix * BMatrix = new CoinPackedMatrix(true, BNumRows, 00184 BNumColsLessTwo, 00185 BFullSizeLessThree, 00186 BElements,BIndices, 00187 BStarts,BLengths); 00188 // Assign problem into a solver interface 00189 // Note: coneSi will cleanup the memory itself 00190 OsiSolverInterface * coneSi = si.clone(false); 00191 coneSi->assignProblem (BMatrix, BColLowers, BColUppers, 00192 BObjective, 00193 BRowLowers, BRowUppers); 00194 00195 // Problem sense should defalut to "min" by default, 00196 // but just to be virtuous... 00197 coneSi->setObjSense(1.0); 00198 00199 // The plot outline from here on down: 00200 // coneSi has been assigned B without the u_0 and v_0 columns 00201 // Calculate base objective <<x^T,Atilde^T>,u> 00202 // bool haveWarmStart = false; 00203 // For (j=0; j<n, j++) 00204 // if (!isBinary(x_j) || x_j<=0 || x_j>=1) continue; 00205 // // IMPROVEME: if(haveWarmStart) check if j attractive 00206 // add {-e_j,0,-1} matrix column for v_0 00207 // add {e_j,0,0} matrix column for u_0 00208 // objective coefficient for u_0 is x_j 00209 // if (haveWarmStart) 00210 // set warmstart info 00211 // solve min{objw:Bw=0; w>=0,except v_0, u_0 free} 00212 // if (bounded) 00213 // get warmstart info 00214 // haveWarmStart=true; 00215 // ustar = optimal u solution 00216 // ustar_0 = optimal u_0 solution 00217 // alpha^T= <ustar^T,Atilde> -ustar_0e_j^T 00218 // (double check <alpha^T,x> >= beta_ should be violated) 00219 // add <alpha^T,x> >= beta_ to cutset 00220 // endif 00221 // delete column for u_0 // this deletes all column info. 00222 // delete column for v_0 00223 // endFor 00224 // clean up memory 00225 // return 0; 00226 00227 int * nVectorIndices = new int[n]; 00228 CoinIotaN(nVectorIndices, n, 0); 00229 00230 bool haveWarmStart = false; 00231 bool equalObj1, equalObj2; 00232 CoinRelFltEq eq; 00233 00234 double v_0Elements[2] = {-1,1}; 00235 double u_0Elements[1] = {1}; 00236 00237 CoinWarmStart * warmStart = 0; 00238 00239 double * ustar = new double[m]; 00240 CoinFillN(ustar, m, 0.0); 00241 00242 double* alpha = new double[n]; 00243 CoinFillN(alpha, n, 0.0); 00244 00245 for (j=0;j<n;j++){ 00246 if (!si.isBinary(j)) continue; // Better to ask coneSi? No! 00247 // coneSi has no binInfo. 00248 equalObj1=eq(x[j],0); 00249 equalObj2=eq(x[j],1); 00250 if (equalObj1 || equalObj2) continue; 00251 // IMPROVEME: if (haveWarmStart) check if j attractive; 00252 00253 // AskLL:wanted to declare u_0 and v_0 packedVec outside loop 00254 // and setIndices, but didn't see a method to do that(?) 00255 // (Could "insert". Seems inefficient) 00256 int v_0Indices[2]={j,nPlus1}; 00257 int u_0Indices[1]={j}; 00258 // 00259 CoinPackedVector v_0(2,v_0Indices,v_0Elements,false); 00260 CoinPackedVector u_0(1,u_0Indices,u_0Elements,false); 00261 00262 #if CGL_DEBUG 00263 const CoinPackedMatrix *see1 = coneSi->getMatrixByRow(); 00264 #endif 00265 00266 coneSi->addCol(v_0,-solverINFINITY,solverINFINITY,0); 00267 coneSi->addCol(u_0,-solverINFINITY,solverINFINITY,x[j]); 00268 if(haveWarmStart) { 00269 coneSi->setWarmStart(warmStart); 00270 coneSi->resolve(); 00271 } 00272 else { 00273 00274 #if CGL_DEBUG 00275 const CoinPackedMatrix *see2 = coneSi->getMatrixByRow(); 00276 #endif 00277 00278 coneSi->initialSolve(); 00279 } 00280 if(coneSi->isProvenOptimal()){ 00281 warmStart = coneSi->getWarmStart(); 00282 haveWarmStart=true; 00283 const double * wstar = coneSi->getColSolution(); 00284 CoinDisjointCopyN(wstar, m, ustar); 00285 Atilde->transposeTimes(ustar,alpha); 00286 alpha[j]+=wstar[BNumCols-1]; 00287 00288 #if debug 00289 int p; 00290 double sum; 00291 for(p=0;p<n;p++)sum+=alpha[p]*x[p]; 00292 if (sum<=beta_){ 00293 throw CoinError("Cut not violated", 00294 "cutGeneration", 00295 "CglLiftAndProject"); 00296 } 00297 #endif 00298 00299 // add <alpha^T,x> >= beta_ to cutset 00300 OsiRowCut rc; 00301 rc.setRow(n,nVectorIndices,alpha); 00302 rc.setLb(beta_); 00303 rc.setUb(solverINFINITY); 00304 cs.insert(rc); 00305 } 00306 // delete col for u_o and v_0 00307 coneSi->deleteCols(2,delCols); 00308 00309 // clean up memory 00310 } 00311 // clean up 00312 delete [] alpha; 00313 delete [] ustar; 00314 delete [] nVectorIndices; 00315 // BMatrix, BColLowers,BColUppers, BObjective, BRowLowers, BRowUppers 00316 // are all freed by OsiSolverInterface destructor (?) 00317 delete [] BLengths; 00318 delete [] BStarts; 00319 delete [] BIndices; 00320 delete [] BElements; 00321 } |
|
Get the normalization : Either beta=+1 or beta=-1. Definition at line 27 of file CglLiftAndProject.hpp. References beta_.
00027 { 00028 return beta_; 00029 } |
|
Set the normalization : Either beta=+1 or beta=-1. Default value is 1. Definition at line 34 of file CglLiftAndProject.hpp. References beta_.
|
|
A function that tests the methods in the CglLiftAndProject class. The only reason for it not to be a member method is that this way it doesn't have to be compiled into the library. And that's a gain, because the library should be compiled with optimization on, but this method should be compiled with debugging. |