00001
00002
00003 #if defined(_MSC_VER)
00004
00005 # pragma warning(disable:4786)
00006 #endif
00007 #include <cstdlib>
00008 #include <cstdio>
00009 #include <cmath>
00010 #include <cassert>
00011 #include <cfloat>
00012 #include <iostream>
00013
00014 #include "CoinHelperFunctions.hpp"
00015 #include "CglLiftAndProject.hpp"
00016 #include "CoinPackedVector.hpp"
00017 #include "CoinSort.hpp"
00018 #include "CoinPackedMatrix.hpp"
00019
00020
00021
00022
00023 void CglLiftAndProject::generateCuts(const OsiSolverInterface & si,
00024 OsiCuts & cs ) const
00025 {
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052 int j = 0;
00053
00054
00055
00056 const int m = si.getNumRows();
00057 const int n = si.getNumCols();
00058 const double * x = si.getColSolution();
00059
00060
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
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
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];
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
00142 int BNumColsLessOne=BNumCols-1;
00143 int BNumColsLessTwo=BNumCols-2;
00144 const int delCols[2] = {BNumColsLessOne, BNumColsLessTwo};
00145
00146
00147
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
00155
00156
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
00168
00169
00170
00171
00172 double * BObjective= new double[BNumCols];
00173 double * Atildex = new double[m];
00174 CoinFillN(BObjective,BNumCols,0.0);
00175 Atilde->times(x,Atildex);
00176 CoinDisjointCopyN(Atildex,m,BObjective);
00177
00178
00179
00180 int BFullSizeLessThree = BFullSize-3;
00181
00182
00183 CoinPackedMatrix * BMatrix = new CoinPackedMatrix(true, BNumRows,
00184 BNumColsLessTwo,
00185 BFullSizeLessThree,
00186 BElements,BIndices,
00187 BStarts,BLengths);
00188
00189
00190 OsiSolverInterface * coneSi = si.clone(false);
00191 coneSi->assignProblem (BMatrix, BColLowers, BColUppers,
00192 BObjective,
00193 BRowLowers, BRowUppers);
00194
00195
00196
00197 coneSi->setObjSense(1.0);
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
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;
00247
00248 equalObj1=eq(x[j],0);
00249 equalObj2=eq(x[j],1);
00250 if (equalObj1 || equalObj2) continue;
00251
00252
00253
00254
00255
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
00300 OsiRowCut rc;
00301 rc.setRow(n,nVectorIndices,alpha);
00302 rc.setLb(beta_);
00303 rc.setUb(solverINFINITY);
00304 cs.insert(rc);
00305 }
00306
00307 coneSi->deleteCols(2,delCols);
00308
00309
00310 }
00311
00312 delete [] alpha;
00313 delete [] ustar;
00314 delete [] nVectorIndices;
00315
00316
00317 delete [] BLengths;
00318 delete [] BStarts;
00319 delete [] BIndices;
00320 delete [] BElements;
00321 }
00322
00323
00324
00325
00326 CglLiftAndProject::CglLiftAndProject ()
00327 :
00328 CglCutGenerator(),
00329 beta_(1),
00330 epsilon_(1.0e-08),
00331 onetol_(1-epsilon_)
00332 {
00333
00334 }
00335
00336
00337
00338
00339 CglLiftAndProject::CglLiftAndProject (const CglLiftAndProject & source) :
00340 CglCutGenerator(source),
00341 beta_(source.beta_),
00342 epsilon_(source.epsilon_),
00343 onetol_(source.onetol_)
00344 {
00345
00346 }
00347
00348
00349
00350
00351 CglLiftAndProject::~CglLiftAndProject ()
00352 {
00353
00354 }
00355
00356
00357
00358
00359 CglLiftAndProject &
00360 CglLiftAndProject::operator=(
00361 const CglLiftAndProject& rhs)
00362 {
00363 if (this != &rhs) {
00364 CglCutGenerator::operator=(rhs);
00365 beta_=rhs.beta_;
00366 epsilon_=rhs.epsilon_;
00367 onetol_=rhs.onetol_;
00368 }
00369 return *this;
00370 }