Changeset 936 for trunk

Timestamp:

Dec 30, 2012 8:25:14 PM (7 years ago)

Author:

pbelotti

Message:

added option to fill matrices for SDP cuts. removed delete_ field from CouenneScalar?

Location:

trunk/Couenne/src

Files:

• cut/sdpcuts/CouenneMatrix.cpp (modified) (1 diff)
• cut/sdpcuts/CouenneMatrix.hpp (modified) (4 diffs)
• cut/sdpcuts/CouenneSdpCuts.cpp (modified) (8 diffs)
• cut/sdpcuts/CouenneSdpCuts.hpp (modified) (1 diff)
• cut/sdpcuts/CutGenSparse.cpp (modified) (6 diffs)
• main/BonCouenneSetup.cpp (modified) (8 diffs)
• problem/CouenneProblem.cpp (modified) (1 diff)
• problem/CouenneProblem.hpp (modified) (2 diffs)
• problem/CouenneProblemConstructors.cpp (modified) (2 diffs)
• problem/CouenneRecordBestSol.cpp (modified) (2 diffs)
• standardize/standardize.cpp (modified) (1 diff)

trunk/Couenne/src/cut/sdpcuts/CouenneMatrix.cpp

@@ -219 +219 @@
 
 
-CouenneScalar::~CouenneScalar () {
-
-  if (delete_) // only delete constants
-    delete elem_;
-} 
+CouenneScalar::~CouenneScalar () 
+{delete elem_;} 
 
 #define WRAP 20

trunk/Couenne/src/cut/sdpcuts/CouenneMatrix.hpp

@@ -29 +29 @@
     int         index_;  ///< index of element in vector
     expression *elem_;   ///< element
-    bool        delete_; ///< destructor should delete pointer to expression
 
   public:
@@ -35 +34 @@
     CouenneScalar (int index, expression *elem):
       index_  (index),
-      elem_   (elem),
-      delete_ (elem_ -> code () == COU_EXPRCONST || (elem_ != elem_ -> Original ())) {}
+      elem_   (elem -> code () == COU_EXPRCONST ? elem : new exprClone (elem)) {}
 
     ~CouenneScalar ();
@@ -42 +40 @@
     CouenneScalar (const CouenneScalar &rhs):
       index_  (rhs.index_),
-      elem_   (new exprClone (rhs.elem_)),
-      delete_ (true) {}
+      elem_   (new exprClone (rhs.elem_)) {}
 
     CouenneScalar &operator= (const CouenneScalar &rhs) {
       index_  = rhs.index_;
       elem_   = new exprClone (rhs.elem_);
-      delete_ = true;
       return *this;
     }
@@ -84 +80 @@
    ~CouenneSparseVector ();
 
-    CouenneSparseVector            (const CouenneSparseVector &rhs);
-    CouenneSparseVector &operator= (const CouenneSparseVector &rhs);
+    CouenneSparseVector                (const CouenneSparseVector &rhs);
+    CouenneSparseVector &operator=     (const CouenneSparseVector &rhs);
     CouenneSparseVector *clone () {return new CouenneSparseVector (*this);}
 

trunk/Couenne/src/cut/sdpcuts/CouenneSdpCuts.cpp

@@ -16 +16 @@
 #include "CouenneProblem.hpp"
 #include "CouenneExprVar.hpp"
+#include "CouenneExprAux.hpp"
+#include "operators/CouenneExprPow.hpp"
+#include "operators/CouenneExprMul.hpp"
 
 #include "CoinTime.hpp"
@@ -32 +35 @@
   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 -> 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;
@@ -75 +79 @@
         int index0 = image -> ArgList () [0] -> Index ();
 
-        if ((index0 >= 0) &&
+        if (        (index0   >= 0) &&
             ((*i) -> Index () >= 0))
 
@@ -96 +100 @@
   // 2.5) If option says so, add fictitious auxiliary variables if not there
 
-  for (std::vector <CouenneExprMatrix *>::iterator 
-         i  = minors_ . begin ();
-       i   != minors_ . end   (); ++i) {
-
-    
+  if (fillMissingTerms_) {
+
+    int nOld = problem_ -> nVars ();
+
+    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>::iterator 
+             rowIt  = (*i) -> getRows (). begin (); 
+           rowIt   != (*i) -> getRows (). end   (); ++rowIt) {
+
+        varNumIndices. insert (rowIt -> first);
+
+        for (std::set <CouenneScalar *, CouenneSparseVector::compare_scalars>::iterator 
+               elemIt  = rowIt -> second -> getElements () . begin ();
+             elemIt   != rowIt -> second -> getElements () . end   (); ++elemIt)
+
+          varNumIndices. insert ((*elemIt) -> getIndex ());
+      }
+
+      printf ("varNumIndices: ");
+
+      for (std::set <int>::iterator 
+             j  = varNumIndices. begin ();
+           j   != varNumIndices. end   (); ++j)
+        printf ("%d ", *j);
+
+      printf ("\n");
+
+      // 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>::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>::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 = new exprAux (image,
+                                          problem_ -> nVars (),
+                                          1 + image -> rank (), 
+                                          exprAux::Unset, 
+                                          problem_ -> domain ());
+
+              // 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
+
+    problem_ -> resizeAuxs (nOld, problem_ -> nVars ());
   }
 
@@ -132 +236 @@
       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 (indexVar, size, *j); // note: problem_ -> Var (indexVar) == (*j)
       (*i) -> add_element (size, indexVar, *j);
     }
@@ -175 +279 @@
 CouenneSdpCuts::CouenneSdpCuts (const CouenneSdpCuts &rhs):
 
-  problem_     (rhs. problem_),
-  doNotUse_    (rhs. doNotUse_),
-  numEigVec_   (rhs. numEigVec_),
-  onlyNegEV_   (rhs. onlyNegEV_),
-  useSparsity_ (rhs. useSparsity_) {
+  problem_          (rhs. problem_),
+  doNotUse_         (rhs. doNotUse_),
+  numEigVec_        (rhs. numEigVec_),
+  onlyNegEV_        (rhs. onlyNegEV_),
+  useSparsity_      (rhs. useSparsity_),
+  fillMissingTerms_ (rhs. fillMissingTerms_) {
 
   for (std::vector <CouenneExprMatrix *>::const_iterator 
@@ -192 +297 @@
 CouenneSdpCuts &CouenneSdpCuts::operator= (const CouenneSdpCuts &rhs) {
 
-  problem_     = rhs. problem_;
-  doNotUse_    = rhs. doNotUse_;
-  numEigVec_   = rhs. numEigVec_;
-  onlyNegEV_   = rhs. onlyNegEV_;
-  useSparsity_ = rhs. useSparsity_;
+  problem_          = rhs. problem_;
+  doNotUse_         = rhs. doNotUse_;
+  numEigVec_        = rhs. numEigVec_;
+  onlyNegEV_        = rhs. onlyNegEV_;
+  useSparsity_      = rhs. useSparsity_;
+  fillMissingTerms_ = rhs. fillMissingTerms_;
 
   for (std::vector <CouenneExprMatrix *>::const_iterator 
@@ -251 +357 @@
     );
 
+  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

trunk/Couenne/src/cut/sdpcuts/CouenneSdpCuts.hpp

@@ -57 +57 @@
 
     bool useSparsity_; ///< Sparsify eigenvalues before writing inequality (default: no)
+
+    bool fillMissingTerms_; ///< If minor not fully dense, create 
+                            ///< fictitious auxiliary variables that 
+                            ///< will be used in sdp cuts only (tighter
+                            ///< than sdp cuts without)
 
   public:

trunk/Couenne/src/cut/sdpcuts/CutGenSparse.cpp

@@ -111 +111 @@
       (*evdec_num)++;
 
-      //printf ("calling dsyevx SPARSIFY 2\n");
-
       //--------------------------------------------------------------------------------------------------------------------------------
       dsyevx_interface (running_n - 1, T, card_ev, w, z, EV_TOL, -COIN_DBL_MAX, 0., 1, (running_n - 1 == min_nz) ? (running_n - 1) : 1);
@@ -125 +123 @@
 
         std::memcpy (Tbest, Tcopy, rnsq            * sizeof (double));
-        CoinCopyN (z, rnsq, zbest); //std::memcpy (zbest, z,     rnsq            * sizeof (double));
+        std::memcpy (zbest, z,     rnsq            * sizeof (double));
         std::memcpy (wbest, w,     (running_n - 1) * sizeof (double));
 
@@ -533 +531 @@
                                int *card_v_mat, int *evdec_num) const {
 
-  int i, j, nchanged = 0,
+  int nchanged = 0,
     min_number_new_per_cut = 1,
     loc_card_new_selected  = 0,
@@ -558 +556 @@
   *card_v_mat = 0;
 
-  for (i=0; i<n; i++) {
+  for (int i=0; i<n; i++) {
 
     order [i] = i;
@@ -605 +603 @@
   while (card_selected < n) {
 
-    for (i=0; i<n; i++)
+    for (int i=0; i<n; i++)
 
       if (selected [order [i]] == 0) {
@@ -632 +630 @@
     // printf (":::::::::::::::::::::::::::::::::\n");
 
-    for(i=0; i<n; i++) {
+    for(int i=0; i<n; i++) {
       if(selected[i] == -1) {
         loc_selected[i] = 0;

trunk/Couenne/src/main/BonCouenneSetup.cpp

@@ -108 +108 @@
                                       CouenneInterface *ci,
                                       Bonmin::Bab *bb) {
+  int freq;
+
   bool retval = true;
 
@@ -120 +122 @@
   /* Get the basic options. */
   readOptionsFile();
- 
+
   // Suppress iteration output from nonlinear solver
   options () -> SetStringValue ("sb", "yes", false, true);
@@ -271 +273 @@
   continuousSolver_ -> messageHandler () -> setLogLevel (lpLogLevel);
 
-  //////////////////////////////////////////////////////////////
-
-  couenneCg -> Problem () -> setMaxCpuTime (getDoubleParameter (BabSetupBase::MaxTime));
-
-  std::string doHeuristic;
-  options () -> GetStringValue ("local_optimization_heuristic", doHeuristic, "couenne.");
-
-  ci -> extractLinearRelaxation (*continuousSolver_, *couenneCg, true, doHeuristic == "yes");
-
-  // In case there are no discrete variables, we have already a
-  // heuristic solution for which create a initialization heuristic
-  if (!(extraStuff -> infeasibleNode ()) &&
-      ci -> isProvenOptimal () && 
-      ci -> haveNlpSolution ()) {
-
-    /// setup initial heuristic (in principle it should only run once...)
-    InitHeuristic* initHeuristic = new InitHeuristic 
-      (ci -> getObjValue (), ci -> getColSolution (), *couenneProb_);
-    HeuristicMethod h;
-    h.id = "Couenne Rounding NLP"; // same name as the real rounding one
-    h.heuristic = initHeuristic;
-    heuristics_.push_back(h);
-  }
-
-  if (extraStuff -> infeasibleNode ()){
-    journalist() -> Printf (J_SUMMARY, J_PROBLEM, "Linear relaxation infeasible, the problem is infeasible.\n");
-    retval = false;
-  }
-
-  //continuousSolver_ -> findIntegersAndSOS (false);
-  //addSos (); // only adds embedded SOS objects
-
   // Add Couenne SOS ///////////////////////////////////////////////////////////////
 
@@ -479 +449 @@
   delete [] objects;
 
-  int freq;
-
   // Setup Fix Point bound tightener /////////////////////////////////////////////
 
@@ -519 +487 @@
 
   CouennePtr_ = couenneCg;
-
-  // Add PSD cut handler ///////////////////////////////////////////////////////
-
-  options () -> GetIntegerValue ("sdp_cuts", freq, "couenne.");
-
-  if (freq != 0) {
-
-    CouenneSdpCuts * couenneSDP = 
-      new CouenneSdpCuts (couenneProb_,
-                          journalist (),
-                          options    ());
-    CuttingMethod cg;
-    cg.frequency = freq;
-    cg.cgl = couenneSDP;
-    cg.id = "Couenne SDP cuts";
-    cutGenerators (). push_back (cg);
-  }
 
   // Add two-inequalities based bound tightening ///////////////////////////////////////////////////////
@@ -557 +508 @@
 
   // Setup heuristic to solve MINLP problems. /////////////////////////////////
+
+  std::string doHeuristic;
+  options () -> GetStringValue ("local_optimization_heuristic", doHeuristic, "couenne.");
+
+  // ----------------------------------------------------------------
+
+  //////////////////////////////////////////////////////////////
+
+  couenneCg -> Problem () -> setMaxCpuTime (getDoubleParameter (BabSetupBase::MaxTime));
+
+  ci -> extractLinearRelaxation (*continuousSolver_, *couenneCg, true, doHeuristic == "yes");
+
+  // In case there are no discrete variables, we have already a
+  // heuristic solution for which create a initialization heuristic
+  if (!(extraStuff -> infeasibleNode ()) &&
+      ci -> isProvenOptimal () && 
+      ci -> haveNlpSolution ()) {
+
+    /// setup initial heuristic (in principle it should only run once...)
+    InitHeuristic* initHeuristic = new InitHeuristic 
+      (ci -> getObjValue (), ci -> getColSolution (), *couenneProb_);
+    HeuristicMethod h;
+    h.id = "Couenne Rounding NLP"; // same name as the real rounding one
+    h.heuristic = initHeuristic;
+    heuristics_.push_back(h);
+  }
+
+  if (extraStuff -> infeasibleNode ()){
+    journalist() -> Printf (J_SUMMARY, J_PROBLEM, "Linear relaxation infeasible, the problem is infeasible.\n");
+    retval = false;
+  }
+
+  // ----------------------------------------------------------------
+
+  //continuousSolver_ -> findIntegersAndSOS (false);
+  //addSos (); // only adds embedded SOS objects
 
   if (doHeuristic == "yes") {
@@ -765 +752 @@
   intParam_ [BabSetupBase::SpecialOption] = 16 | 4;
 
+  // Add PSD cut handler ///////////////////////////////////////////////////////
+
+  options () -> GetIntegerValue ("sdp_cuts", freq, "couenne.");
+
+  if (freq != 0) {
+
+    CouenneSdpCuts * couenneSDP = 
+      new CouenneSdpCuts (couenneProb_,
+                          journalist (),
+                          options    ());
+    CuttingMethod cg;
+    cg.frequency = freq;
+    cg.cgl = couenneSDP;
+    cg.id = "Couenne SDP cuts";
+    cutGenerators (). push_back (cg);
+  }
+
   // Add disjunctive cuts ///////////////////////////////////////////////////////
 
@@ -803 +807 @@
     cutGenerators (). push_back(cg);
   }
-
-  // // Add disjunctive cuts ///////////////////////////////////////////////////////
-
-  // options () -> GetIntegerValue ("minlp_disj_cuts", freq, "couenne.");
-
-  // if (freq != 0) {
-
-  //   CouenneDisjCuts * couenneDisj = 
-  //     new CouenneDisjCuts (ci, this, 
-  //                       couenneCg, 
-  //                       branchingMethod_, 
-  //                       varSelection == OSI_STRONG, // if true, use strong branching candidates
-  //                       journalist (),
-  //                       options ());
-
-  //   CuttingMethod cg;
-  //   cg.frequency = freq;
-  //   cg.cgl = couenneDisj;
-  //   cg.id = "Couenne disjunctive cuts";
-  //   cutGenerators (). push_back(cg);
-  // }
 
   return retval;

trunk/Couenne/src/problem/CouenneProblem.cpp

@@ -238 +238 @@
 }
 
+/// rescans problem after adding new auxiliaries
+void CouenneProblem::resizeAuxs (int nOld, int nNew) {
+
+#define resizeOld(typeD,oldV,oldN,newN) { \
+  if (oldV) {                             \
+    typeD *newV = new typeD [newN];       \
+    CoinCopyN (oldV, oldN, newV);         \
+    delete [] oldV;                       \
+    oldV = newV;                          \
+  }                                       \
+}
+
+  resizeOld (int,    numbering_,   nOld, nNew);
+  resizeOld (bool,   commuted_,    nOld, nNew);
+  resizeOld (int,    integerRank_, nOld, nNew);
+  //resizeOld (double, optimum_,     nOld, nNew);
+
+  optimum_ = (double *) realloc (optimum_, nNew * sizeof (double));
+
+  if (numbering_)   for (int i=nOld; i<nNew; ++i) numbering_   [i] = i;
+  if (commuted_)    for (int i=nOld; i<nNew; ++i) commuted_    [i] = false;
+  if (integerRank_) for (int i=nOld; i<nNew; ++i) integerRank_ [i] = nNew;
+  if (optimum_)     for (int i=nOld; i<nNew; ++i) optimum_     [i] = 0.;
+
+  domain_ . current () -> resize (nNew);
+
+  // post-rescan: update 
+  //
+  // numbering_
+  // domain_
+  // commuted_
+  // optimum_
+  // integerRank_
+  //
+  // since there are new variables
+}
 
 /// Get cutoff

trunk/Couenne/src/problem/CouenneProblem.hpp

@@ -217 +217 @@
   CouNumber bestObj_;
 
-  /// Indices of variables appearing in products (used for SDP cuts)
-  int *quadIndex_;
-
   /// Variables that have commuted to auxiliary
   bool *commuted_;
@@ -725 +722 @@
   inline double constObjVal () const {return constObjVal_;}
 
+  /// refills auxiliaries vectors after adding auxiliaries
+  void resizeAuxs (int nOld, int nNew);
+
 protected:
 

trunk/Couenne/src/problem/CouenneProblemConstructors.cpp

@@ -55 +55 @@
   optimum_   (NULL),
   bestObj_   (COIN_DBL_MAX),
-  quadIndex_ (NULL),
   commuted_  (NULL),
   numbering_ (NULL),
@@ -219 +218 @@
 CouenneProblem::~CouenneProblem () {
 
+  delete auxSet_;
+
   if (perfIndicator_)
     delete perfIndicator_;

trunk/Couenne/src/problem/CouenneRecordBestSol.cpp

@@ -207 +207 @@
   else {
     if(givenCard != cardSol) {
-      printf("CouenneRecordBestSol::setSol(): ### ERROR: givenCard: %d  cardSol: %d", givenCard, cardSol);
-      exit(1);
+      //printf("CouenneRecordBestSol::setSol(): ### ERROR: givenCard: %d  cardSol: %d", givenCard, cardSol);
+      //exit(1);
+
+        double *newSol = new double [givenCard];
+        CoinCopyN (givenSol, givenCard, newSol);
+        delete [] modSol;
+        modSol = newSol;
+        cardSol = givenCard;
+
     }
   }
@@ -341 +348 @@
     else {
       if(givenModCard != cardModSol) {
-        printf("CouenneRecordBestSol::setModSol(): ### ERROR: givenModCard: %d  cardModSol: %d", givenModCard, cardModSol);
-        exit(1);
+        // printf("CouenneRecordBestSol::setModSol(): ### ERROR: givenModCard: %d  cardModSol: %d", givenModCard, cardModSol);
+        // exit(1);
+
+        double *newModSol = new double [givenModCard];
+        CoinCopyN (givenModSol, givenModCard, newModSol);
+        delete [] modSol;
+        modSol = newModSol;
+        cardModSol = givenModCard;
       }
     }

trunk/Couenne/src/standardize/standardize.cpp

@@ -476 +476 @@
 #endif
 
-  delete auxSet_;
+  //  delete auxSet_;
 
   // Create evaluation order ////////////////////////////////////////////////////