#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. |
1.3.5