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