disjCut.cpp

Go to the documentation of this file.

/* $Id: disjCut.cpp 694 2011-06-18 20:13:17Z stefan $
 *
 * Name:    disjCut.cpp
 * Author:  Pietro Belotti
 * Purpose: generate one disjunctive cut based on a single disjunction
 *
 * (C) Carnegie-Mellon University, 2008. 
 * This file is licensed under the Eclipse Public License (EPL)
 */

#include "CoinTime.hpp"
#include "OsiClpSolverInterface.hpp"
#include "CouennePrecisions.hpp"
#include "CouenneDisjCuts.hpp"
//#include "CouenneCutGenerator.hpp"
//#include "CouenneProblem.hpp"
#include "CouenneInfeasCut.hpp"

using namespace Ipopt;
using namespace Couenne;

#define COEFF_BOUNDS 1.e10

#define MIN_NUM_COEFF 1.e-9
#define MAX_NUM_COEFF 1.e+9
#define MAX_NUM_RATIO 1.e+6


// print content of sparse matrix
void printMatrix (int nrows, int ncols, int nel, 
                  const int *start, const int *len, 
                  const int *ind, const double *el);

// same with CoinPackedMatrix
void printMatrix   (const CoinPackedMatrix *A);

// same with si.GetMatrixByRow()
void printLPMatrix (const OsiSolverInterface &si);

// add columns specified by start/len/ind/el to matrix Astd
void addSubMatr (int *start, int *len, int *ind, double *el, 
                 CoinPackedMatrix &Astd, CoinPackedVector &rstd, 
                 int &cur, int &curCol, int dispM, int dispVec, int nrows);


int CouenneDisjCuts::generateDisjCuts (std::vector <std::pair <OsiCuts *, OsiCuts *> > &disjunctions, 
                                       OsiSolverInterface &si, 
                                       OsiCuts &cs, 
                                       const CglTreeInfo &info) const {

  // create CGLP with si+left and si+right, for each (left,right) in
  // the vector of disjunctions
  //
  // Here's the CGLP (row and column order as shown)
  //
  //        {n}          {1}    {m}    {m}   {mL} {mR}
  //
  // max    alpha xbar - beta                                        --- maximize so can copy xbar
  // s.t.  -alpha              + u A       + u'C           =  0      --- n rows
  //       -alpha                    + v A       + v'D     =  0      --- n
  //                    -beta  + u b       + u'c          <=  0      --- 1
  //                    -beta        + v b       + v'd    <=  0      --- 1
  //                           |(u,    v,    u',   v')|_1  =  1      --- normalized multipliers
  //
  //                             u,    v,    u',   v'     >=  0      --- non-negativity 
  //
  //
  //  And here are the submatrices (M' is M transposed)
  //
  //     First        Second  Third  Fourth Fifth Sixth   --- Columns
  //
  //       xbar          -1       0     0     0     0
  //      ---------------------------------------------
  //       -I_n           .       A'    .     C'    .     = 0
  //       -I_n           .       .     A'    .     D'    = 0
  //        .            -1       b'    .     c'    .    <= 0
  //        .            -1       .     b'    .     d'   <= 0
  //        .             .       e     e     e     e     = 1
  //
  //
  // build a different A such that
  // 
  // - only active rows (resp. columns) of A are included if active_rows_
  //   (resp. active_columns_) are set
  // - equality constraints are duplicated in a <= and a >= constraint
  // - >= constraints are inverted
  //
  // Also, add disjunctive cut to CGLP for use with next disjunction

  // put matrix from base problem in canonical form //////////////////////////////
  CoinPackedMatrix Astd;
  CoinPackedVector rstd;
  OsiSI2MatrVec (Astd, rstd, si); 

  int
    n   = si.   getNumCols (),      //mC   = 2*n + 3,
    m   = Astd. getMajorDim (),     nC   = 1 + n + 2 * m,
    nnz = Astd. getNumElements (),  nnzC = 2 * (n + 1 + nnz + 2*m);

  if (jnlst_ -> ProduceOutput (J_DETAILED, J_DISJCUTS))
    printf ("canonical form has %d cols, %d rows, %d nonzeros --> cglp has %d,%d,%d\n", 
            n, m, nnz, nC, 2*n + 3, nnzC);

  double 
    *elements = new double [nnzC];

  int 
    *indices = new int [nnzC],
    *start   = new int [nC + 1],
    *length  = new int [nC],
    cur      = 0,
    curCol   = 0;

  // first column: two identity matrices
  for (int i=0, i2 = n; i<n;) {
    start   [curCol]   = cur;
    length  [curCol++] = 2;
    indices [cur] =    i++; elements [cur++] = -1.;
    indices [cur] =   i2++; elements [cur++] = -1.;
  }

  // second column: two "-1" at position 2n and 2n+1
  start   [curCol]   = cur;
  length  [curCol++] = 2;
  indices [cur] = 2*n;   elements [cur++] = -1.;
  indices [cur] = 2*n+1; elements [cur++] = -1.;

  // third...
  addSubMatr (start + curCol, length + curCol, 
              indices + cur, elements + cur, 
              Astd, rstd,
              cur, curCol, 
              0, 2*n,   2*n+2);

  if (jnlst_ -> ProduceOutput (J_MATRIX, J_DISJCUTS)) {
    printf ("with third column\n");
    printMatrix (curCol, 2*n+3, cur, start, length, indices, elements);
  }

  // ... and fourth column: get single rows from Astd
  addSubMatr (start + curCol, length + curCol, 
              indices + cur, elements + cur, 
              Astd, rstd, 
              cur, curCol, 
              n, 2*n+1, 2*n+2);

  if (jnlst_ -> ProduceOutput (J_MATRIX, J_DISJCUTS)) {
    printf ("with 4th column\n");
    printMatrix (curCol, 2*n+3, cur, start, length, indices, elements);
  }

  CoinPackedMatrix *baseA = new CoinPackedMatrix;

  baseA -> assignMatrix (true,     // column ordered
                         2*n+3,    // minor dimension
                         curCol,   // major dimension
                         cur,      // number of elements
                         elements, // elements
                         indices,  // indices
                         start,    // starting positions
                         length);  // length

  //printf ("should be copy of above\n");
  //printMatrix (baseA);

  OsiClpSolverInterface cglp;

  cglp. messageHandler () -> setLogLevel (0);

  int 
    N = baseA -> getMajorDim (),        // # cols in base problem
    M = baseA -> getMinorDim ();        // # rows in base problem

  assert (M == 2 * n + 3);

  // vectors of the problem
  double
    *collb  = new double [N], // variable lower bounds
    *colub  = new double [N], // variable upper bounds
    *obj    = new double [N], // objective coefficients
    *rowrhs = new double [M], // right hand sides (all zero except the last, 1)
    *rowrng = new double [M]; // row range (empty)

  // bounds
  CoinFillN (collb,       n+1,       -COEFF_BOUNDS);
  CoinFillN (collb + n+1, N - (n+1),  0.);
  CoinFillN (colub,       n+1,        COEFF_BOUNDS);
  CoinFillN (colub + n+1, N - (n+1),  1.);

  // objective coefficients
  CoinCopyN (si.getColSolution (), n, obj);
  obj [n] = -1.;
  CoinFillN (obj + (n+1), N-(n+1),    0.);

  // rhs
  CoinFillN (rowrhs,      M-1,        0.);
  rowrhs [M-1] = 1.;

  // rhs range
  CoinFillN (rowrng,      M,          COIN_DBL_MAX);

  // signs of the inequalities
  char *rowsen = new char [M];
  CoinFillN (rowsen, M, 'E');
  rowsen [M-3] = rowsen [M-2] = 'L';

  cglp.assignProblem (baseA,   // matrix
                      collb,   // lower bounds 
                      colub,   // upper bounds
                      obj,     // obj coefficients
                      rowsen,  // row sense
                      rowrhs,  // right hand sides
                      rowrng); // no row range

  // this is a maximization problem
  cglp. setObjSense (-1);


  // generate and solve one CGLP for each disjunction

  bool first = true;

  for (std::vector <std::pair <OsiCuts *, OsiCuts *> >::iterator disjI = disjunctions.begin ();
       (disjI != disjunctions.end ()) && (CoinCpuTime () < cpuTime_); ++disjI) {

    OsiCuts
      *left  = disjI -> first,
      *right = disjI -> second;

    int 
      ncolLeft  = OsiCuts2MatrVec (&cglp,  left, 0, 2*n),
      ncolRight = OsiCuts2MatrVec (&cglp, right, n, 2*n+1);    

    /*char filename [30];
    static int iter = 0;
    sprintf (filename, "cglp-%04d-%04d-%04d", info.level, iter++, info.pass);
    cglp.writeLp (filename);*/

    if (jnlst_ -> ProduceOutput (J_MOREMATRIX, J_DISJCUTS)) {
      printf ("current CGLP:\n");
      printLPMatrix (cglp);
    }

    //cglp.options () -> SetNumericValue ("max_cpu_time", couenneCG_ -> Problem () -> getMaxCpuTime () - CoinCpuTime ());
    if (first) {cglp.initialSolve (); first = false;}
    else        cglp.resolve (); // segfault in ex1244

    if (cglp. isProvenOptimal () && (cglp.getObjValue () > COUENNE_EPS)) {

      const double *AlphaBeta = cglp. getColSolution ();

      int    *colInd = NULL, nnzCut = 0;
      double *colCoe = NULL;

      // count nonzero entries, compute ratio max/min coefficient
      double mincoeff = COIN_DBL_MAX, maxcoeff = 0.;

      for (register int i=n+1; i--;) {
        double value = fabs (AlphaBeta [i]);
        if (value == 0.) continue;
        if (value > maxcoeff) maxcoeff = value;
        if (value < mincoeff) mincoeff = value;
        if ((maxcoeff            > MAX_NUM_COEFF) ||
            (maxcoeff            < MIN_NUM_COEFF) ||
            (maxcoeff / mincoeff > MAX_NUM_RATIO)) 
          break;
        nnzCut++;
      }

      if (nnzCut &&
          (maxcoeff            < MAX_NUM_COEFF) &&
          (maxcoeff            > MIN_NUM_COEFF) &&
          (maxcoeff / mincoeff < MAX_NUM_RATIO)) {

        // cut data
        double *nzcoeff = new double [nnzCut];
        int    *indices = new int    [nnzCut];

        // fill in indices and coefficient
        for (int i = nnzCut = 0; i<n; i++)
          if (fabs (AlphaBeta [i]) > MIN_NUM_COEFF) {
            indices [nnzCut]   = i;
            nzcoeff [nnzCut++] = AlphaBeta [i];
          }

        OsiRowCut *cut = new OsiRowCut;
        cut -> setRow (nnzCut, indices, nzcoeff);
        cut -> setUb  (AlphaBeta [n]);

        /*if (1) {

          printf ("---- RESOLVING\n");
          si.applyRowCuts (1, cut);
          si.writeLp ("added");
          si.resolve ();
          printf ("---- RESOLVED\n");

          double *obj    = new double [N]; // objective coefficients

          // objective coefficients
          CoinCopyN (si.getColSolution (), n, obj);
          obj [n] = -1.;
          CoinFillN (obj + (n+1), N-(n+1),    0.);

          cglp.setObjective (obj);
          }*/
 
        // add it to CGLP
        if (addPreviousCut_) {

          colInd = new int    [2 * (nnzCut + 2)];
          colCoe = new double [2 * (nnzCut + 2)];

          // first column
          CoinCopyN    (nzcoeff, nnzCut, colCoe);
          CoinCopyN    (indices, nnzCut, colInd); 
          colInd [nnzCut]       = 2*n;   colCoe [nnzCut]   = AlphaBeta [n];
          colInd [nnzCut+1]     = 2*n+2; colCoe [nnzCut+1] = 1; // entry in norm constraint

          // second column
          CoinCopyN    (nzcoeff, nnzCut, colCoe + nnzCut + 2);
          CoinCopyDisp (indices, nnzCut, colInd + nnzCut + 2, n); 
          colInd [2*nnzCut + 2] = 2*n+1; colCoe [2*nnzCut+2] = AlphaBeta [n];
          colInd [2*nnzCut + 3] = 2*n+2; colCoe [2*nnzCut+3] = 1.; // entry in norm constraint

          // extra vectors
          double lb  [2] = {0., 0.};
          double ub  [2] = {1., 1.};
          double obj [2] = {0., 0.};

          int start [3];
          *start = 0;
          start [2] = 2 * (start [1] = nnzCut + 2);

          cglp. addCols (2,        // const int numcols, 
                         start,    // const int* columnStarts,
                         colInd,   // const int* rows, 
                         colCoe,   // const double* elements,
                         lb,       // const double* collb, 
                         ub,       // const double* colub,   
                         obj);     // const double* obj

          delete [] colCoe;
          delete [] colInd;
        }

        delete [] nzcoeff;
        delete [] indices;

        if (jnlst_ -> ProduceOutput (J_DETAILED, J_DISJCUTS)) {
          printf ("====== disjunctive cut: "); 
          cut -> print ();
        }

        // add cut to cs
        cs. insert (cut);
      }
    }

    // remove last ncolLeft + ncolRight columns from cglp
    int *delIndices = new int [ncolLeft + ncolRight];
    for (register int nc = ncolLeft + ncolRight, j = N + nc; nc--;)
      *delIndices++ = --j;
    delIndices -= (ncolLeft + ncolRight);
    cglp.deleteCols (ncolLeft + ncolRight, delIndices);
    delete [] delIndices;
  }

  return COUENNE_FEASIBLE;
}