CouenneSdpCuts.cpp

Go to the documentation of this file.

/* $Id: CouenneSdpCuts.cpp 1259 2018-08-28 10:34:22Z pbelotti $
 *
 * Name:    CouenneSdpCuts.cpp
 * Author:  Pietro Belotti
 * Purpose: wrapper for Couenne to insert sdpcuts
 *
 * This file is licensed under the Eclipse Public License (EPL)
 */

#include "IpOptionsList.hpp"

#include "CouenneJournalist.hpp"
#include "CouenneMatrix.hpp"
#include "CouennePSDcon.hpp"
#include "CouenneSdpCuts.hpp"
#include "CouenneProblem.hpp"
#include "CouenneExprVar.hpp"
#include "CouenneExprAux.hpp"
#include "operators/CouenneExprPow.hpp"
#include "operators/CouenneExprMul.hpp"

#include "CoinTime.hpp"

#include <set>

using namespace Couenne;

//#define DEBUG

CouenneSdpCuts::CouenneSdpCuts (CouenneProblem *p,
                                JnlstPtr jnlst,
                                const Ipopt::SmartPtr <Ipopt::OptionsList> options):
  problem_  (p),
  doNotUse_ (false) {

  std::string s;

  options -> GetIntegerValue ("sdp_cuts_num_ev",      numEigVec_, "couenne.");
  options -> GetStringValue  ("sdp_cuts_neg_ev",      s,          "couenne."); onlyNegEV_        = (s == "yes");
  options -> GetStringValue  ("sdp_cuts_sparsify",    s,          "couenne."); useSparsity_      = (s == "yes");
  options -> GetStringValue  ("sdp_cuts_fillmissing", s,          "couenne."); fillMissingTerms_ = (s == "yes");

  CouenneExprMatrix *cauldron = new CouenneExprMatrix;

  // 1) Construct matrix with entries x_i, x_j

  for (std::vector <exprVar *>::iterator 
         i  = p -> Variables (). begin (); 
       i   != p -> Variables (). end   (); ++i)

    if ((*i) -> Type () == AUX) {

      expression *image = (*i) -> Image ();

      if ((image -> code () == COU_EXPRMUL) && 
          (image -> ArgList () [0] -> Type () != CONST) && 
          (image -> ArgList () [1] -> Type () != CONST)) {

        int 
          index0 = image -> ArgList () [0] -> Index (),
          index1 = image -> ArgList () [1] -> Index ();

        if ((index0 >= 0) &&
            (index1 >= 0) &&
            ((*i) -> Index () >= 0)) {

          // cauldron -> add_element (CoinMin (index0, index1),
          //                          CoinMax (index0, index1), (*i));

          cauldron -> add_element (index0, index1, (*i));
          cauldron -> add_element (index1, index0, (*i));
        }
      }

      if ((image -> code () == COU_EXPRPOW) && 
          (image -> ArgList () [0] -> Type  () != CONST) && 
          (image -> ArgList () [1] -> Value () == 2.)) {

        int index0 = image -> ArgList () [0] -> Index ();

        if (        (index0   >= 0) &&
            ((*i) -> Index () >= 0))

          cauldron -> add_element (index0, index0, (*i));
      }
    }

#ifdef DEBUG
  printf ("Cauldron so far\n");
  cauldron -> print ();
#endif

  // TODO

  // 2) Block-partition it (optional), obtain matrices. Replace line
  //    below to decompose cauldron

  minors_ . push_back (cauldron);

  // 2.5) If option says so, add fictitious auxiliary variables if not there

  if (fillMissingTerms_) {

    for (std::vector <CouenneExprMatrix *>::iterator 
           i  = minors_ . begin ();
         i   != minors_ . end   (); ++i) {

      // First: construct (possibly sparse) index set, to check
      // against each row (whose index set is a SUBSET, if non-proper,
      // of varNumIndices). Do not use varIndices, as it is empty for
      // now.

      std::set <int> varNumIndices;

      for (std::set <std::pair <int, CouenneSparseVector *>, CouenneExprMatrix::compare_pair_ind>::const_iterator 
             rowIt  = (*i) -> getRows (). begin ();
           rowIt   != (*i) -> getRows (). end   (); ++rowIt) {

        varNumIndices. insert (rowIt -> first);

        for (std::set <CouenneScalar *, CouenneSparseVector::compare_scalars>::const_iterator 
               elemIt  = rowIt -> second -> getElements () . begin ();
             elemIt   != rowIt -> second -> getElements () . end   (); ++elemIt)

          varNumIndices. insert ((*elemIt) -> getIndex ());
      }

      // Second: check every row for elements (i,j) not in this row by
      // parallel scanning of varNumINdices

      for (std::set <std::pair <int, CouenneSparseVector *>, CouenneExprMatrix::compare_pair_ind>::const_iterator 
             rowIt  = (*i) -> getRows (). begin (); 
           rowIt   != (*i) -> getRows (). end   (); ++rowIt) {

        int rowInd = rowIt -> first;

        std::set <int>::iterator vniIt = varNumIndices . begin ();

        for (std::set <CouenneScalar *, CouenneSparseVector::compare_scalars>::const_iterator 
               elemIt  = rowIt -> second -> getElements () . begin ();
             elemIt   != rowIt -> second -> getElements () . end   (); ++elemIt) {

          int colInd = (*elemIt) -> getIndex ();

          while ((vniIt != varNumIndices . end ()) && (*vniIt < colInd)) {

            if (rowInd <= *vniIt) {
              //printf ("missing term: %d, %d\n", rowInd, *vniIt);

              expression *image;

              if (rowInd == *vniIt) image = new exprPow (new exprClone (problem_ -> Var (rowInd)), new exprConst (2.));
              else                  image = new exprMul (new exprClone (problem_ -> Var (rowInd)), new exprClone (problem_ -> Var (*vniIt)));

              exprAux *yIJ = problem_ -> addAuxiliary (image);

              // seek expression in the set
              if (problem_ -> AuxSet () -> find (yIJ) == 
                  problem_ -> AuxSet () -> end ()) {

                // no such expression found in the set, create entry therein
                problem_ -> Variables () . push_back (yIJ);
                problem_ -> AuxSet () -> insert (yIJ); // insert into repetition checking structure
              }

              (*i) -> add_element (rowInd, *vniIt, yIJ);
              (*i) -> add_element (*vniIt, rowInd, yIJ);
            }

            ++vniIt;
          }

          if (vniIt == varNumIndices . end ())
            break;
          else ++vniIt;
        }
      }
    }

    // post-rescan: update 
    //
    // numbering_
    // domain_
    // commuted_
    // optimum_
    // integerRank_
    // unusedOriginalsIndices_
    //
    // since there are new variables
  }

  // 3) Bottom-right border each block with a row vector, a column vector,
  // and the constant 1

  for (std::vector <CouenneExprMatrix *>::iterator 
         i  = minors_ . begin ();
       i   != minors_ . end   (); ++i) {

    int size = problem_ -> nVars ();

#ifdef DEBUG
    printf ("minor has %ld rows and %ld columns\n",
            (*i) -> getRows () . size (),
            (*i) -> getCols () . size ());
#endif

    for (std::set <std::pair <int, CouenneSparseVector *>, CouenneExprMatrix::compare_pair_ind>::const_iterator 
           j  = (*i) -> getCols () . begin (); 
         j   != (*i) -> getCols () . end   (); ++j)

      (*i) -> varIndices () . push_back (problem_ -> Var (j -> first));

    for (std::vector <expression *>::iterator 
           j  = (*i) -> varIndices () . begin ();
         j   != (*i) -> varIndices () . end   (); ++j) {

      int indexVar = (*j) -> Index ();
#ifdef DEBUG
      printf ("adding at [%d,%d] and viceversa\n", indexVar, size);
#endif
      (*i) -> add_element (indexVar, size, *j); // note: problem_ -> Var (indexVar) == (*j)
      (*i) -> add_element (size, indexVar, *j);
    }

    (*i) -> add_element (size, size, new exprConst (1.));

#ifdef DEBUG
    (*i) -> print ();
#endif
  }

  // 0) Search for X \succeq 0 constraints, if any, then add matrix to
  //    minors for each such constraint

  if (p -> ConstraintClass ("PSDcon"))
    for (std::vector <CouenneConstraint *>::iterator 
           i  = p -> ConstraintClass ("PSDcon") -> begin (); 
         i   != p -> ConstraintClass ("PSDcon") -> end   (); ++i) {

      CouennePSDcon *con = dynamic_cast <CouennePSDcon *> (*i);

      if (!con) 
        continue;

      minors_ . push_back (con -> getX ());
    }
}


CouenneSdpCuts::~CouenneSdpCuts () {

  for (std::vector <CouenneExprMatrix *>::iterator 
         i  = minors_ . begin ();
       i   != minors_ . end   (); ++i)

    delete (*i);
}


CouenneSdpCuts::CouenneSdpCuts (const CouenneSdpCuts &rhs):

  problem_          (rhs. problem_),
  doNotUse_         (rhs. doNotUse_),
  numEigVec_        (rhs. numEigVec_),
  onlyNegEV_        (rhs. onlyNegEV_),
  useSparsity_      (rhs. useSparsity_),
  fillMissingTerms_ (rhs. fillMissingTerms_) {

  for (std::vector <CouenneExprMatrix *>::const_iterator 
         i  = rhs.minors_ . begin ();
       i   != rhs.minors_ . end   (); ++i)

    minors_ . push_back (new CouenneExprMatrix (**i));
}


CouenneSdpCuts &CouenneSdpCuts::operator= (const CouenneSdpCuts &rhs) {

  problem_          = rhs. problem_;
  doNotUse_         = rhs. doNotUse_;
  numEigVec_        = rhs. numEigVec_;
  onlyNegEV_        = rhs. onlyNegEV_;
  useSparsity_      = rhs. useSparsity_;
  fillMissingTerms_ = rhs. fillMissingTerms_;

  for (std::vector <CouenneExprMatrix *>::const_iterator 
         i  = rhs.minors_ . begin ();
       i   != rhs.minors_ . end   (); ++i)

    minors_ . push_back (new CouenneExprMatrix (**i));

  return *this;
}


CglCutGenerator *CouenneSdpCuts::clone () const
{return new CouenneSdpCuts (*this);}


void CouenneSdpCuts::registerOptions (Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions) {

  roptions -> AddLowerBoundedIntegerOption
    ("sdp_cuts",
     "The frequency (in terms of nodes) at which Couenne SDP cuts are generated.",
     -99, 0,
     "A frequency of 0 (default) means these cuts are never generated. \
Any positive number n instructs Couenne to generate them at every n nodes of the B&B tree. \
A negative number -n means that generation should be attempted at the root node, and if successful it can be repeated at every n nodes, otherwise it is stopped altogether."
    );

  roptions -> AddLowerBoundedIntegerOption
    ("sdp_cuts_num_ev",
     "The number of eigenvectors of matrix X to be used to create sdp cuts.",
     -1, -1,
     "Set to -1 to indicate that all n eigenvectors should be used. Eigenvalues are \
sorted in non-decreasing order, hence selecting 1 will provide cuts on the most negative eigenvalue."
    );

  roptions -> AddStringOption2
    ("sdp_cuts_neg_ev",
     "Only use negative eigenvalues to create sdp cuts.",
     "yes", 
     "no", "use all eigenvalues regardless of their sign.",
     "yes", "exclude all non-negative eigenvalues."
    );


  roptions -> AddStringOption2
    ("sdp_cuts_sparsify",
     "Make cuts sparse by greedily reducing X one column at a time before extracting eigenvectors.",
     "no", 
     "no", "",
     "yes", ""
    );

  roptions -> AddStringOption2
    ("sdp_cuts_fillmissing",
     "Create fictitious auxiliary variables to fill non-fully dense minors. Can make a difference when Q has at least one zero term.",
     "no", 
     "no", "Do not create auxiliaries and simply use Fourier-Motzkin to substitute a missing auxiliary y_ij with inequalities that use bounds and the definition y_ij = x_i x_j \
Advantage: limits the creation of auxiliaries, reformulation stays small. Default.",
     "yes", "Create (at the beginning) auxiliaries that are linearized (through McCormick) and used within an SDP cut. This allows tighter cuts although it increases the size \
of the reformulation and hence of the linear relaxation."
    );

#if 0
  roptions -> AddStringOption2
    ("sdp_cuts_",
     "",
     "yes", 
     "no", "",
     "yes", ""
    );
#endif
}