00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "CoinTime.hpp"
00012 #include "OsiClpSolverInterface.hpp"
00013 #include "CouennePrecisions.hpp"
00014 #include "CouenneDisjCuts.hpp"
00015
00016 #define COEFF_BOUNDS 1.e10
00017
00018 #define MIN_NUM_COEFF 1.e-9
00019 #define MAX_NUM_COEFF 1.e+9
00020 #define MAX_NUM_RATIO 1.e+6
00021
00022
00023
00024 void printMatrix (int nrows, int ncols, int nel,
00025 const int *start, const int *len,
00026 const int *ind, const double *el);
00027
00028
00029 void printMatrix (const CoinPackedMatrix *A);
00030
00031
00032 void printLPMatrix (const OsiSolverInterface &si);
00033
00034
00035 void addSubMatr (int *start, int *len, int *ind, double *el,
00036 CoinPackedMatrix &Astd, CoinPackedVector &rstd,
00037 int &cur, int &curCol, int dispM, int dispVec, int nrows);
00038
00039
00041 int CouenneDisjCuts::generateDisjCuts (std::vector <std::pair <OsiCuts *, OsiCuts *> > &disjunctions,
00042 OsiSolverInterface &si,
00043 OsiCuts &cs,
00044 const CglTreeInfo &info) const {
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086 CoinPackedMatrix Astd;
00087 CoinPackedVector rstd;
00088 OsiSI2MatrVec (Astd, rstd, si);
00089
00090 int
00091 n = si. getNumCols (),
00092 m = Astd. getMajorDim (), nC = 1 + n + 2 * m,
00093 nnz = Astd. getNumElements (), nnzC = 2 * (n + 1 + nnz + 2*m);
00094
00095 if (jnlst_ -> ProduceOutput (J_DETAILED, J_DISJCUTS))
00096 printf ("canonical form has %d cols, %d rows, %d nonzeros --> cglp has %d,%d,%d\n",
00097 n, m, nnz, nC, 2*n + 3, nnzC);
00098
00099 double
00100 *elements = new double [nnzC];
00101
00102 int
00103 *indices = new int [nnzC],
00104 *start = new int [nC + 1],
00105 *length = new int [nC],
00106 cur = 0,
00107 curCol = 0;
00108
00109
00110 for (int i=0, i2 = n; i<n;) {
00111 start [curCol] = cur;
00112 length [curCol++] = 2;
00113 indices [cur] = i++; elements [cur++] = -1.;
00114 indices [cur] = i2++; elements [cur++] = -1.;
00115 }
00116
00117
00118 start [curCol] = cur;
00119 length [curCol++] = 2;
00120 indices [cur] = 2*n; elements [cur++] = -1.;
00121 indices [cur] = 2*n+1; elements [cur++] = -1.;
00122
00123
00124 addSubMatr (start + curCol, length + curCol,
00125 indices + cur, elements + cur,
00126 Astd, rstd,
00127 cur, curCol,
00128 0, 2*n, 2*n+2);
00129
00130 if (jnlst_ -> ProduceOutput (J_MATRIX, J_DISJCUTS)) {
00131 printf ("with third column\n");
00132 printMatrix (curCol, 2*n+3, cur, start, length, indices, elements);
00133 }
00134
00135
00136 addSubMatr (start + curCol, length + curCol,
00137 indices + cur, elements + cur,
00138 Astd, rstd,
00139 cur, curCol,
00140 n, 2*n+1, 2*n+2);
00141
00142 if (jnlst_ -> ProduceOutput (J_MATRIX, J_DISJCUTS)) {
00143 printf ("with 4th column\n");
00144 printMatrix (curCol, 2*n+3, cur, start, length, indices, elements);
00145 }
00146
00147 CoinPackedMatrix *baseA = new CoinPackedMatrix;
00148
00149 baseA -> assignMatrix (true,
00150 2*n+3,
00151 curCol,
00152 cur,
00153 elements,
00154 indices,
00155 start,
00156 length);
00157
00158
00159
00160
00161 OsiClpSolverInterface cglp;
00162
00163 cglp. messageHandler () -> setLogLevel (0);
00164
00165 int
00166 N = baseA -> getMajorDim (),
00167 M = baseA -> getMinorDim ();
00168
00169 assert (M == 2 * n + 3);
00170
00171
00172 double
00173 *collb = new double [N],
00174 *colub = new double [N],
00175 *obj = new double [N],
00176 *rowrhs = new double [M],
00177 *rowrng = new double [M];
00178
00179
00180 CoinFillN (collb, n+1, -COEFF_BOUNDS);
00181 CoinFillN (collb + n+1, N - (n+1), 0.);
00182 CoinFillN (colub, n+1, COEFF_BOUNDS);
00183 CoinFillN (colub + n+1, N - (n+1), 1.);
00184
00185
00186 CoinCopyN (si.getColSolution (), n, obj);
00187 obj [n] = -1.;
00188 CoinFillN (obj + (n+1), N-(n+1), 0.);
00189
00190
00191 CoinFillN (rowrhs, M-1, 0.);
00192 rowrhs [M-1] = 1.;
00193
00194
00195 CoinFillN (rowrng, M, COIN_DBL_MAX);
00196
00197
00198 char *rowsen = new char [M];
00199 CoinFillN (rowsen, M, 'E');
00200 rowsen [M-3] = rowsen [M-2] = 'L';
00201
00202 cglp.assignProblem (baseA,
00203 collb,
00204 colub,
00205 obj,
00206 rowsen,
00207 rowrhs,
00208 rowrng);
00209
00210
00211 cglp. setObjSense (-1);
00212
00214
00215
00216
00217 bool first = true;
00218
00219 for (std::vector <std::pair <OsiCuts *, OsiCuts *> >::iterator disjI = disjunctions.begin ();
00220 (disjI != disjunctions.end ()) && (CoinCpuTime () < cpuTime_); ++disjI) {
00221
00222 OsiCuts
00223 *left = disjI -> first,
00224 *right = disjI -> second;
00225
00226 int
00227 ncolLeft = OsiCuts2MatrVec (&cglp, left, 0, 2*n),
00228 ncolRight = OsiCuts2MatrVec (&cglp, right, n, 2*n+1);
00229
00230
00231
00232
00233
00234
00235 if (jnlst_ -> ProduceOutput (J_MOREMATRIX, J_DISJCUTS)) {
00236 printf ("current CGLP:\n");
00237 printLPMatrix (cglp);
00238 }
00239
00240 if (first) {cglp.initialSolve (); first = false;}
00241 else cglp.resolve ();
00242
00243 if (cglp. isProvenOptimal () && (cglp.getObjValue () > COUENNE_EPS)) {
00244
00245 const double *AlphaBeta = cglp. getColSolution ();
00246
00247 int *colInd = NULL, nnzCut = 0;
00248 double *colCoe = NULL;
00249
00250
00251 double mincoeff = COIN_DBL_MAX, maxcoeff = 0.;
00252
00253 for (register int i=n+1; i--;) {
00254 double value = fabs (AlphaBeta [i]);
00255 if (value == 0.) continue;
00256 if (value > maxcoeff) maxcoeff = value;
00257 if (value < mincoeff) mincoeff = value;
00258 if ((maxcoeff > MAX_NUM_COEFF) ||
00259 (maxcoeff < MIN_NUM_COEFF) ||
00260 (maxcoeff / mincoeff > MAX_NUM_RATIO))
00261 break;
00262 nnzCut++;
00263 }
00264
00265 if (nnzCut &&
00266 (maxcoeff < MAX_NUM_COEFF) &&
00267 (maxcoeff > MIN_NUM_COEFF) &&
00268 (maxcoeff / mincoeff < MAX_NUM_RATIO)) {
00269
00270
00271 double *nzcoeff = new double [nnzCut];
00272 int *indices = new int [nnzCut];
00273
00274
00275 for (int i = nnzCut = 0; i<n; i++)
00276 if (fabs (AlphaBeta [i]) > MIN_NUM_COEFF) {
00277 indices [nnzCut] = i;
00278 nzcoeff [nnzCut++] = AlphaBeta [i];
00279 }
00280
00281 OsiRowCut *cut = new OsiRowCut;
00282 cut -> setRow (nnzCut, indices, nzcoeff);
00283 cut -> setUb (AlphaBeta [n]);
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304 if (addPreviousCut_) {
00305
00306 colInd = new int [2 * (nnzCut + 2)];
00307 colCoe = new double [2 * (nnzCut + 2)];
00308
00309
00310 CoinCopyN (nzcoeff, nnzCut, colCoe);
00311 CoinCopyN (indices, nnzCut, colInd);
00312 colInd [nnzCut] = 2*n; colCoe [nnzCut] = AlphaBeta [n];
00313 colInd [nnzCut+1] = 2*n+2; colCoe [nnzCut+1] = 1;
00314
00315
00316 CoinCopyN (nzcoeff, nnzCut, colCoe + nnzCut + 2);
00317 CoinCopyDisp (indices, nnzCut, colInd + nnzCut + 2, n);
00318 colInd [2*nnzCut + 2] = 2*n+1; colCoe [2*nnzCut+2] = AlphaBeta [n];
00319 colInd [2*nnzCut + 3] = 2*n+2; colCoe [2*nnzCut+3] = 1.;
00320
00321
00322 double lb [2] = {0., 0.};
00323 double ub [2] = {1., 1.};
00324 double obj [2] = {0., 0.};
00325
00326 int start [3];
00327 *start = 0;
00328 start [2] = 2 * (start [1] = nnzCut + 2);
00329
00330 cglp. addCols (2,
00331 start,
00332 colInd,
00333 colCoe,
00334 lb,
00335 ub,
00336 obj);
00337
00338 delete [] colCoe;
00339 delete [] colInd;
00340 }
00341
00342 delete [] nzcoeff;
00343 delete [] indices;
00344
00345 if (jnlst_ -> ProduceOutput (J_DETAILED, J_DISJCUTS)) {
00346 printf ("====== disjunctive cut: ");
00347 cut -> print ();
00348 }
00349
00350
00351 cs. insert (cut);
00352 }
00353 }
00354
00355
00356 int *delIndices = new int [ncolLeft + ncolRight];
00357 for (register int nc = ncolLeft + ncolRight, j = N + nc; nc--;)
00358 *delIndices++ = --j;
00359 delIndices -= (ncolLeft + ncolRight);
00360 cglp.deleteCols (ncolLeft + ncolRight, delIndices);
00361 delete [] delIndices;
00362 }
00363
00364 return COUENNE_FEASIBLE;
00365 }