/home/coin/SVN-release/OS-2.4.0/Couenne/src/disjunctive/disjCut.cpp

Go to the documentation of this file.
00001 /* $Id: disjCut.cpp 694 2011-06-18 20:13:17Z stefan $
00002  *
00003  * Name:    disjCut.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: generate one disjunctive cut based on a single disjunction
00006  *
00007  * (C) Carnegie-Mellon University, 2008. 
00008  * This file is licensed under the Eclipse Public License (EPL)
00009  */
00010 
00011 #include "CoinTime.hpp"
00012 #include "OsiClpSolverInterface.hpp"
00013 #include "CouennePrecisions.hpp"
00014 #include "CouenneDisjCuts.hpp"
00015 //#include "CouenneCutGenerator.hpp"
00016 //#include "CouenneProblem.hpp"
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 // print content of sparse matrix
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 // same with CoinPackedMatrix
00035 void printMatrix   (const CoinPackedMatrix *A);
00036 
00037 // same with si.GetMatrixByRow()
00038 void printLPMatrix (const OsiSolverInterface &si);
00039 
00040 // add columns specified by start/len/ind/el to matrix Astd
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   // create CGLP with si+left and si+right, for each (left,right) in
00053   // the vector of disjunctions
00054   //
00055   // Here's the CGLP (row and column order as shown)
00056   //
00057   //        {n}          {1}    {m}    {m}   {mL} {mR}
00058   //
00059   // max    alpha xbar - beta                                        --- maximize so can copy xbar
00060   // s.t.  -alpha              + u A       + u'C           =  0      --- n rows
00061   //       -alpha                    + v A       + v'D     =  0      --- n
00062   //                    -beta  + u b       + u'c          <=  0      --- 1
00063   //                    -beta        + v b       + v'd    <=  0      --- 1
00064   //                           |(u,    v,    u',   v')|_1  =  1      --- normalized multipliers
00065   //
00066   //                             u,    v,    u',   v'     >=  0      --- non-negativity 
00067   //
00068   //
00069   //  And here are the submatrices (M' is M transposed)
00070   //
00071   //     First        Second  Third  Fourth Fifth Sixth   --- Columns
00072   //
00073   //       xbar          -1       0     0     0     0
00074   //      ---------------------------------------------
00075   //       -I_n           .       A'    .     C'    .     = 0
00076   //       -I_n           .       .     A'    .     D'    = 0
00077   //        .            -1       b'    .     c'    .    <= 0
00078   //        .            -1       .     b'    .     d'   <= 0
00079   //        .             .       e     e     e     e     = 1
00080   //
00081   //
00082   // build a different A such that
00083   // 
00084   // - only active rows (resp. columns) of A are included if active_rows_
00085   //   (resp. active_columns_) are set
00086   // - equality constraints are duplicated in a <= and a >= constraint
00087   // - >= constraints are inverted
00088   //
00089   // Also, add disjunctive cut to CGLP for use with next disjunction
00090 
00091   // put matrix from base problem in canonical form //////////////////////////////
00092   CoinPackedMatrix Astd;
00093   CoinPackedVector rstd;
00094   OsiSI2MatrVec (Astd, rstd, si); 
00095 
00096   int
00097     n   = si.   getNumCols (),      //mC   = 2*n + 3,
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   // first column: two identity matrices
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   // second column: two "-1" at position 2n and 2n+1
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   // third...
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   // ... and fourth column: get single rows from Astd
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,     // column ordered
00156                          2*n+3,    // minor dimension
00157                          curCol,   // major dimension
00158                          cur,      // number of elements
00159                          elements, // elements
00160                          indices,  // indices
00161                          start,    // starting positions
00162                          length);  // length
00163 
00164   //printf ("should be copy of above\n");
00165   //printMatrix (baseA);
00166 
00167   OsiClpSolverInterface cglp;
00168 
00169   cglp. messageHandler () -> setLogLevel (0);
00170 
00171   int 
00172     N = baseA -> getMajorDim (),        // # cols in base problem
00173     M = baseA -> getMinorDim ();        // # rows in base problem
00174 
00175   assert (M == 2 * n + 3);
00176 
00177   // vectors of the problem
00178   double
00179     *collb  = new double [N], // variable lower bounds
00180     *colub  = new double [N], // variable upper bounds
00181     *obj    = new double [N], // objective coefficients
00182     *rowrhs = new double [M], // right hand sides (all zero except the last, 1)
00183     *rowrng = new double [M]; // row range (empty)
00184 
00185   // bounds
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   // objective coefficients
00192   CoinCopyN (si.getColSolution (), n, obj);
00193   obj [n] = -1.;
00194   CoinFillN (obj + (n+1), N-(n+1),    0.);
00195 
00196   // rhs
00197   CoinFillN (rowrhs,      M-1,        0.);
00198   rowrhs [M-1] = 1.;
00199 
00200   // rhs range
00201   CoinFillN (rowrng,      M,          COIN_DBL_MAX);
00202 
00203   // signs of the inequalities
00204   char *rowsen = new char [M];
00205   CoinFillN (rowsen, M, 'E');
00206   rowsen [M-3] = rowsen [M-2] = 'L';
00207 
00208   cglp.assignProblem (baseA,   // matrix
00209                       collb,   // lower bounds 
00210                       colub,   // upper bounds
00211                       obj,     // obj coefficients
00212                       rowsen,  // row sense
00213                       rowrhs,  // right hand sides
00214                       rowrng); // no row range
00215 
00216   // this is a maximization problem
00217   cglp. setObjSense (-1);
00218 
00220 
00221   // generate and solve one CGLP for each disjunction
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     /*char filename [30];
00237     static int iter = 0;
00238     sprintf (filename, "cglp-%04d-%04d-%04d", info.level, iter++, info.pass);
00239     cglp.writeLp (filename);*/
00240 
00241     if (jnlst_ -> ProduceOutput (J_MOREMATRIX, J_DISJCUTS)) {
00242       printf ("current CGLP:\n");
00243       printLPMatrix (cglp);
00244     }
00245 
00246     //cglp.options () -> SetNumericValue ("max_cpu_time", couenneCG_ -> Problem () -> getMaxCpuTime () - CoinCpuTime ());
00247     if (first) {cglp.initialSolve (); first = false;}
00248     else        cglp.resolve (); // segfault in ex1244
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       // count nonzero entries, compute ratio max/min coefficient
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         // cut data
00278         double *nzcoeff = new double [nnzCut];
00279         int    *indices = new int    [nnzCut];
00280 
00281         // fill in indices and coefficient
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         /*if (1) {
00293 
00294           printf ("---- RESOLVING\n");
00295           si.applyRowCuts (1, cut);
00296           si.writeLp ("added");
00297           si.resolve ();
00298           printf ("---- RESOLVED\n");
00299 
00300           double *obj    = new double [N]; // objective coefficients
00301 
00302           // objective coefficients
00303           CoinCopyN (si.getColSolution (), n, obj);
00304           obj [n] = -1.;
00305           CoinFillN (obj + (n+1), N-(n+1),    0.);
00306 
00307           cglp.setObjective (obj);
00308           }*/
00309  
00310         // add it to CGLP
00311         if (addPreviousCut_) {
00312 
00313           colInd = new int    [2 * (nnzCut + 2)];
00314           colCoe = new double [2 * (nnzCut + 2)];
00315 
00316           // first column
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; // entry in norm constraint
00321 
00322           // second column
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.; // entry in norm constraint
00327 
00328           // extra vectors
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,        // const int numcols, 
00338                          start,    // const int* columnStarts,
00339                          colInd,   // const int* rows, 
00340                          colCoe,   // const double* elements,
00341                          lb,       // const double* collb, 
00342                          ub,       // const double* colub,   
00343                          obj);     // const double* 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         // add cut to cs
00358         cs. insert (cut);
00359       }
00360     }
00361 
00362     // remove last ncolLeft + ncolRight columns from cglp
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 }

Generated on Thu Sep 22 03:05:57 2011 by  doxygen 1.4.7