Changeset 936 for trunk


Ignore:
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:
11 edited

Legend:

Unmodified
Added
Removed
  • trunk/Couenne/src/cut/sdpcuts/CouenneMatrix.cpp

    r932 r936  
    219219
    220220
    221 CouenneScalar::~CouenneScalar () {
    222 
    223   if (delete_) // only delete constants
    224     delete elem_;
    225 }
     221CouenneScalar::~CouenneScalar ()
     222{delete elem_;}
    226223
    227224#define WRAP 20
  • trunk/Couenne/src/cut/sdpcuts/CouenneMatrix.hpp

    r934 r936  
    2929    int         index_;  ///< index of element in vector
    3030    expression *elem_;   ///< element
    31     bool        delete_; ///< destructor should delete pointer to expression
    3231
    3332  public:
     
    3534    CouenneScalar (int index, expression *elem):
    3635      index_  (index),
    37       elem_   (elem),
    38       delete_ (elem_ -> code () == COU_EXPRCONST || (elem_ != elem_ -> Original ())) {}
     36      elem_   (elem -> code () == COU_EXPRCONST ? elem : new exprClone (elem)) {}
    3937
    4038    ~CouenneScalar ();
     
    4240    CouenneScalar (const CouenneScalar &rhs):
    4341      index_  (rhs.index_),
    44       elem_   (new exprClone (rhs.elem_)),
    45       delete_ (true) {}
     42      elem_   (new exprClone (rhs.elem_)) {}
    4643
    4744    CouenneScalar &operator= (const CouenneScalar &rhs) {
    4845      index_  = rhs.index_;
    4946      elem_   = new exprClone (rhs.elem_);
    50       delete_ = true;
    5147      return *this;
    5248    }
     
    8480   ~CouenneSparseVector ();
    8581
    86     CouenneSparseVector            (const CouenneSparseVector &rhs);
    87     CouenneSparseVector &operator= (const CouenneSparseVector &rhs);
     82    CouenneSparseVector                (const CouenneSparseVector &rhs);
     83    CouenneSparseVector &operator=     (const CouenneSparseVector &rhs);
    8884    CouenneSparseVector *clone () {return new CouenneSparseVector (*this);}
    8985
  • trunk/Couenne/src/cut/sdpcuts/CouenneSdpCuts.cpp

    r935 r936  
    1616#include "CouenneProblem.hpp"
    1717#include "CouenneExprVar.hpp"
     18#include "CouenneExprAux.hpp"
     19#include "operators/CouenneExprPow.hpp"
     20#include "operators/CouenneExprMul.hpp"
    1821
    1922#include "CoinTime.hpp"
     
    3235  std::string s;
    3336
    34   options -> GetIntegerValue ("sdp_cuts_num_ev",   numEigVec_, "couenne.");
    35   options -> GetStringValue  ("sdp_cuts_neg_ev",   s,          "couenne."); onlyNegEV_   = (s == "yes");
    36   options -> GetStringValue  ("sdp_cuts_sparsify", s,          "couenne."); useSparsity_ = (s == "yes");
     37  options -> GetIntegerValue ("sdp_cuts_num_ev",      numEigVec_, "couenne.");
     38  options -> GetStringValue  ("sdp_cuts_neg_ev",      s,          "couenne."); onlyNegEV_        = (s == "yes");
     39  options -> GetStringValue  ("sdp_cuts_sparsify",    s,          "couenne."); useSparsity_      = (s == "yes");
     40  options -> GetStringValue  ("sdp_cuts_fillmissing", s,          "couenne."); fillMissingTerms_ = (s == "yes");
    3741
    3842  CouenneExprMatrix *cauldron = new CouenneExprMatrix;
     
    7579        int index0 = image -> ArgList () [0] -> Index ();
    7680
    77         if ((index0 >= 0) &&
     81        if (        (index0  >= 0) &&
    7882            ((*i) -> Index () >= 0))
    7983
     
    96100  // 2.5) If option says so, add fictitious auxiliary variables if not there
    97101
    98   for (std::vector <CouenneExprMatrix *>::iterator
    99          i  = minors_ . begin ();
    100        i   != minors_ . end   (); ++i) {
    101 
    102    
     102  if (fillMissingTerms_) {
     103
     104    int nOld = problem_ -> nVars ();
     105
     106    for (std::vector <CouenneExprMatrix *>::iterator
     107           i  = minors_ . begin ();
     108         i   != minors_ . end   (); ++i) {
     109
     110      // First: construct (possibly sparse) index set, to check
     111      // against each row (whose index set is a SUBSET, if non-proper,
     112      // of varNumIndices). Do not use varIndices, as it is empty for
     113      // now.
     114
     115      std::set <int> varNumIndices;
     116
     117      for (std::set <std::pair <int, CouenneSparseVector *>, CouenneExprMatrix::compare_pair_ind>::iterator
     118             rowIt  = (*i) -> getRows (). begin ();
     119           rowIt   != (*i) -> getRows (). end   (); ++rowIt) {
     120
     121        varNumIndices. insert (rowIt -> first);
     122
     123        for (std::set <CouenneScalar *, CouenneSparseVector::compare_scalars>::iterator
     124               elemIt  = rowIt -> second -> getElements () . begin ();
     125             elemIt   != rowIt -> second -> getElements () . end   (); ++elemIt)
     126
     127          varNumIndices. insert ((*elemIt) -> getIndex ());
     128      }
     129
     130      printf ("varNumIndices: ");
     131
     132      for (std::set <int>::iterator
     133             j  = varNumIndices. begin ();
     134           j   != varNumIndices. end   (); ++j)
     135        printf ("%d ", *j);
     136
     137      printf ("\n");
     138
     139      // Second: check every row for elements (i,j) not in this row by
     140      // parallel scanning of varNumINdices
     141
     142      for (std::set <std::pair <int, CouenneSparseVector *>, CouenneExprMatrix::compare_pair_ind>::iterator
     143             rowIt  = (*i) -> getRows (). begin ();
     144           rowIt   != (*i) -> getRows (). end   (); ++rowIt) {
     145
     146        int rowInd = rowIt -> first;
     147
     148        std::set <int>::iterator vniIt = varNumIndices . begin ();
     149
     150        for (std::set <CouenneScalar *, CouenneSparseVector::compare_scalars>::iterator
     151               elemIt  = rowIt -> second -> getElements () . begin ();
     152             elemIt   != rowIt -> second -> getElements () . end   (); ++elemIt) {
     153
     154          int colInd = (*elemIt) -> getIndex ();
     155
     156          while ((vniIt != varNumIndices . end ()) && (*vniIt < colInd)) {
     157
     158            if (rowInd <= *vniIt) {
     159              printf ("missing term: %d, %d\n", rowInd, *vniIt);
     160
     161              expression *image;
     162
     163              if (rowInd == *vniIt) image = new exprPow (new exprClone (problem_ -> Var (rowInd)), new exprConst (2.));
     164              else                  image = new exprMul (new exprClone (problem_ -> Var (rowInd)), new exprClone (problem_ -> Var (*vniIt)));
     165
     166              exprAux *yIJ = new exprAux (image,
     167                                          problem_ -> nVars (),
     168                                          1 + image -> rank (),
     169                                          exprAux::Unset,
     170                                          problem_ -> domain ());
     171
     172              // seek expression in the set
     173              if (problem_ -> AuxSet () -> find (yIJ) ==
     174                  problem_ -> AuxSet () -> end ()) {
     175
     176                // no such expression found in the set, create entry therein
     177                problem_ -> Variables () . push_back (yIJ);
     178                problem_ -> AuxSet () -> insert (yIJ); // insert into repetition checking structure
     179              }
     180
     181              (*i) -> add_element (rowInd, *vniIt, yIJ);
     182              (*i) -> add_element (*vniIt, rowInd, yIJ);
     183            }
     184
     185            ++vniIt;
     186          }
     187
     188          if (vniIt == varNumIndices . end ())
     189            break;
     190          else ++vniIt;
     191        }
     192      }
     193    }
     194
     195    // post-rescan: update
     196    //
     197    // numbering_
     198    // domain_
     199    // commuted_
     200    // optimum_
     201    // integerRank_
     202    // unusedOriginalsIndices_
     203    //
     204    // since there are new variables
     205
     206    problem_ -> resizeAuxs (nOld, problem_ -> nVars ());
    103207  }
    104208
     
    132236      printf ("adding at [%d,%d] and viceversa\n", indexVar, size);
    133237#endif
    134       (*i) -> add_element (indexVar, size, *j); // note: problem_ -> Var (indexVar) = (*j)
     238      (*i) -> add_element (indexVar, size, *j); // note: problem_ -> Var (indexVar) == (*j)
    135239      (*i) -> add_element (size, indexVar, *j);
    136240    }
     
    175279CouenneSdpCuts::CouenneSdpCuts (const CouenneSdpCuts &rhs):
    176280
    177   problem_     (rhs. problem_),
    178   doNotUse_    (rhs. doNotUse_),
    179   numEigVec_   (rhs. numEigVec_),
    180   onlyNegEV_   (rhs. onlyNegEV_),
    181   useSparsity_ (rhs. useSparsity_) {
     281  problem_          (rhs. problem_),
     282  doNotUse_         (rhs. doNotUse_),
     283  numEigVec_        (rhs. numEigVec_),
     284  onlyNegEV_        (rhs. onlyNegEV_),
     285  useSparsity_      (rhs. useSparsity_),
     286  fillMissingTerms_ (rhs. fillMissingTerms_) {
    182287
    183288  for (std::vector <CouenneExprMatrix *>::const_iterator
     
    192297CouenneSdpCuts &CouenneSdpCuts::operator= (const CouenneSdpCuts &rhs) {
    193298
    194   problem_     = rhs. problem_;
    195   doNotUse_    = rhs. doNotUse_;
    196   numEigVec_   = rhs. numEigVec_;
    197   onlyNegEV_   = rhs. onlyNegEV_;
    198   useSparsity_ = rhs. useSparsity_;
     299  problem_          = rhs. problem_;
     300  doNotUse_         = rhs. doNotUse_;
     301  numEigVec_        = rhs. numEigVec_;
     302  onlyNegEV_        = rhs. onlyNegEV_;
     303  useSparsity_      = rhs. useSparsity_;
     304  fillMissingTerms_ = rhs. fillMissingTerms_;
    199305
    200306  for (std::vector <CouenneExprMatrix *>::const_iterator
     
    251357    );
    252358
     359  roptions -> AddStringOption2
     360    ("sdp_cuts_fillmissing",
     361     "Create fictitious auxiliary variables to fill non-fully dense minors. Can make a difference when Q has at least one zero term.",
     362     "no",
     363     "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 \
     364Advantage: limits the creation of auxiliaries, reformulation stays small. Default.",
     365     "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 \
     366of the reformulation and hence of the linear relaxation."
     367    );
     368
    253369#if 0
    254370  roptions -> AddStringOption2
  • trunk/Couenne/src/cut/sdpcuts/CouenneSdpCuts.hpp

    r935 r936  
    5757
    5858    bool useSparsity_; ///< Sparsify eigenvalues before writing inequality (default: no)
     59
     60    bool fillMissingTerms_; ///< If minor not fully dense, create
     61                            ///< fictitious auxiliary variables that
     62                            ///< will be used in sdp cuts only (tighter
     63                            ///< than sdp cuts without)
    5964
    6065  public:
  • trunk/Couenne/src/cut/sdpcuts/CutGenSparse.cpp

    r935 r936  
    111111      (*evdec_num)++;
    112112
    113       //printf ("calling dsyevx SPARSIFY 2\n");
    114 
    115113      //--------------------------------------------------------------------------------------------------------------------------------
    116114      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);
     
    125123
    126124        std::memcpy (Tbest, Tcopy, rnsq            * sizeof (double));
    127         CoinCopyN (z, rnsq, zbest); //std::memcpy (zbest, z,     rnsq            * sizeof (double));
     125        std::memcpy (zbest, z,     rnsq            * sizeof (double));
    128126        std::memcpy (wbest, w,     (running_n - 1) * sizeof (double));
    129127
     
    533531                               int *card_v_mat, int *evdec_num) const {
    534532
    535   int i, j, nchanged = 0,
     533  int nchanged = 0,
    536534    min_number_new_per_cut = 1,
    537535    loc_card_new_selected  = 0,
     
    558556  *card_v_mat = 0;
    559557
    560   for (i=0; i<n; i++) {
     558  for (int i=0; i<n; i++) {
    561559
    562560    order [i] = i;
     
    605603  while (card_selected < n) {
    606604
    607     for (i=0; i<n; i++)
     605    for (int i=0; i<n; i++)
    608606
    609607      if (selected [order [i]] == 0) {
     
    632630    // printf (":::::::::::::::::::::::::::::::::\n");
    633631
    634     for(i=0; i<n; i++) {
     632    for(int i=0; i<n; i++) {
    635633      if(selected[i] == -1) {
    636634        loc_selected[i] = 0;
  • trunk/Couenne/src/main/BonCouenneSetup.cpp

    r933 r936  
    108108                                      CouenneInterface *ci,
    109109                                      Bonmin::Bab *bb) {
     110  int freq;
     111
    110112  bool retval = true;
    111113
     
    120122  /* Get the basic options. */
    121123  readOptionsFile();
    122  
     124
    123125  // Suppress iteration output from nonlinear solver
    124126  options () -> SetStringValue ("sb", "yes", false, true);
     
    271273  continuousSolver_ -> messageHandler () -> setLogLevel (lpLogLevel);
    272274
    273   //////////////////////////////////////////////////////////////
    274 
    275   couenneCg -> Problem () -> setMaxCpuTime (getDoubleParameter (BabSetupBase::MaxTime));
    276 
    277   std::string doHeuristic;
    278   options () -> GetStringValue ("local_optimization_heuristic", doHeuristic, "couenne.");
    279 
    280   ci -> extractLinearRelaxation (*continuousSolver_, *couenneCg, true, doHeuristic == "yes");
    281 
    282   // In case there are no discrete variables, we have already a
    283   // heuristic solution for which create a initialization heuristic
    284   if (!(extraStuff -> infeasibleNode ()) &&
    285       ci -> isProvenOptimal () &&
    286       ci -> haveNlpSolution ()) {
    287 
    288     /// setup initial heuristic (in principle it should only run once...)
    289     InitHeuristic* initHeuristic = new InitHeuristic
    290       (ci -> getObjValue (), ci -> getColSolution (), *couenneProb_);
    291     HeuristicMethod h;
    292     h.id = "Couenne Rounding NLP"; // same name as the real rounding one
    293     h.heuristic = initHeuristic;
    294     heuristics_.push_back(h);
    295   }
    296 
    297   if (extraStuff -> infeasibleNode ()){
    298     journalist() -> Printf (J_SUMMARY, J_PROBLEM, "Linear relaxation infeasible, the problem is infeasible.\n");
    299     retval = false;
    300   }
    301 
    302   //continuousSolver_ -> findIntegersAndSOS (false);
    303   //addSos (); // only adds embedded SOS objects
    304 
    305275  // Add Couenne SOS ///////////////////////////////////////////////////////////////
    306276
     
    479449  delete [] objects;
    480450
    481   int freq;
    482 
    483451  // Setup Fix Point bound tightener /////////////////////////////////////////////
    484452
     
    519487
    520488  CouennePtr_ = couenneCg;
    521 
    522   // Add PSD cut handler ///////////////////////////////////////////////////////
    523 
    524   options () -> GetIntegerValue ("sdp_cuts", freq, "couenne.");
    525 
    526   if (freq != 0) {
    527 
    528     CouenneSdpCuts * couenneSDP =
    529       new CouenneSdpCuts (couenneProb_,
    530                           journalist (),
    531                           options    ());
    532     CuttingMethod cg;
    533     cg.frequency = freq;
    534     cg.cgl = couenneSDP;
    535     cg.id = "Couenne SDP cuts";
    536     cutGenerators (). push_back (cg);
    537   }
    538489
    539490  // Add two-inequalities based bound tightening ///////////////////////////////////////////////////////
     
    557508
    558509  // Setup heuristic to solve MINLP problems. /////////////////////////////////
     510
     511  std::string doHeuristic;
     512  options () -> GetStringValue ("local_optimization_heuristic", doHeuristic, "couenne.");
     513
     514  // ----------------------------------------------------------------
     515
     516  //////////////////////////////////////////////////////////////
     517
     518  couenneCg -> Problem () -> setMaxCpuTime (getDoubleParameter (BabSetupBase::MaxTime));
     519
     520  ci -> extractLinearRelaxation (*continuousSolver_, *couenneCg, true, doHeuristic == "yes");
     521
     522  // In case there are no discrete variables, we have already a
     523  // heuristic solution for which create a initialization heuristic
     524  if (!(extraStuff -> infeasibleNode ()) &&
     525      ci -> isProvenOptimal () &&
     526      ci -> haveNlpSolution ()) {
     527
     528    /// setup initial heuristic (in principle it should only run once...)
     529    InitHeuristic* initHeuristic = new InitHeuristic
     530      (ci -> getObjValue (), ci -> getColSolution (), *couenneProb_);
     531    HeuristicMethod h;
     532    h.id = "Couenne Rounding NLP"; // same name as the real rounding one
     533    h.heuristic = initHeuristic;
     534    heuristics_.push_back(h);
     535  }
     536
     537  if (extraStuff -> infeasibleNode ()){
     538    journalist() -> Printf (J_SUMMARY, J_PROBLEM, "Linear relaxation infeasible, the problem is infeasible.\n");
     539    retval = false;
     540  }
     541
     542  // ----------------------------------------------------------------
     543
     544  //continuousSolver_ -> findIntegersAndSOS (false);
     545  //addSos (); // only adds embedded SOS objects
    559546
    560547  if (doHeuristic == "yes") {
     
    765752  intParam_ [BabSetupBase::SpecialOption] = 16 | 4;
    766753
     754  // Add PSD cut handler ///////////////////////////////////////////////////////
     755
     756  options () -> GetIntegerValue ("sdp_cuts", freq, "couenne.");
     757
     758  if (freq != 0) {
     759
     760    CouenneSdpCuts * couenneSDP =
     761      new CouenneSdpCuts (couenneProb_,
     762                          journalist (),
     763                          options    ());
     764    CuttingMethod cg;
     765    cg.frequency = freq;
     766    cg.cgl = couenneSDP;
     767    cg.id = "Couenne SDP cuts";
     768    cutGenerators (). push_back (cg);
     769  }
     770
    767771  // Add disjunctive cuts ///////////////////////////////////////////////////////
    768772
     
    803807    cutGenerators (). push_back(cg);
    804808  }
    805 
    806   // // Add disjunctive cuts ///////////////////////////////////////////////////////
    807 
    808   // options () -> GetIntegerValue ("minlp_disj_cuts", freq, "couenne.");
    809 
    810   // if (freq != 0) {
    811 
    812   //   CouenneDisjCuts * couenneDisj =
    813   //     new CouenneDisjCuts (ci, this,
    814   //                       couenneCg,
    815   //                       branchingMethod_,
    816   //                       varSelection == OSI_STRONG, // if true, use strong branching candidates
    817   //                       journalist (),
    818   //                       options ());
    819 
    820   //   CuttingMethod cg;
    821   //   cg.frequency = freq;
    822   //   cg.cgl = couenneDisj;
    823   //   cg.id = "Couenne disjunctive cuts";
    824   //   cutGenerators (). push_back(cg);
    825   // }
    826809
    827810  return retval;
  • trunk/Couenne/src/problem/CouenneProblem.cpp

    r752 r936  
    238238}
    239239
     240/// rescans problem after adding new auxiliaries
     241void CouenneProblem::resizeAuxs (int nOld, int nNew) {
     242
     243#define resizeOld(typeD,oldV,oldN,newN) { \
     244  if (oldV) {                             \
     245    typeD *newV = new typeD [newN];       \
     246    CoinCopyN (oldV, oldN, newV);         \
     247    delete [] oldV;                       \
     248    oldV = newV;                          \
     249  }                                       \
     250}
     251
     252  resizeOld (int,    numbering_,   nOld, nNew);
     253  resizeOld (bool,   commuted_,    nOld, nNew);
     254  resizeOld (int,    integerRank_, nOld, nNew);
     255  //resizeOld (double, optimum_,     nOld, nNew);
     256
     257  optimum_ = (double *) realloc (optimum_, nNew * sizeof (double));
     258
     259  if (numbering_)   for (int i=nOld; i<nNew; ++i) numbering_   [i] = i;
     260  if (commuted_)    for (int i=nOld; i<nNew; ++i) commuted_    [i] = false;
     261  if (integerRank_) for (int i=nOld; i<nNew; ++i) integerRank_ [i] = nNew;
     262  if (optimum_)     for (int i=nOld; i<nNew; ++i) optimum_     [i] = 0.;
     263
     264  domain_ . current () -> resize (nNew);
     265
     266  // post-rescan: update
     267  //
     268  // numbering_
     269  // domain_
     270  // commuted_
     271  // optimum_
     272  // integerRank_
     273  //
     274  // since there are new variables
     275}
    240276
    241277/// Get cutoff
  • trunk/Couenne/src/problem/CouenneProblem.hpp

    r910 r936  
    217217  CouNumber bestObj_;
    218218
    219   /// Indices of variables appearing in products (used for SDP cuts)
    220   int *quadIndex_;
    221 
    222219  /// Variables that have commuted to auxiliary
    223220  bool *commuted_;
     
    725722  inline double constObjVal () const {return constObjVal_;}
    726723
     724  /// refills auxiliaries vectors after adding auxiliaries
     725  void resizeAuxs (int nOld, int nNew);
     726
    727727protected:
    728728
  • trunk/Couenne/src/problem/CouenneProblemConstructors.cpp

    r821 r936  
    5555  optimum_   (NULL),
    5656  bestObj_   (COIN_DBL_MAX),
    57   quadIndex_ (NULL),
    5857  commuted_  (NULL),
    5958  numbering_ (NULL),
     
    219218CouenneProblem::~CouenneProblem () {
    220219
     220  delete auxSet_;
     221
    221222  if (perfIndicator_)
    222223    delete perfIndicator_;
  • trunk/Couenne/src/problem/CouenneRecordBestSol.cpp

    r821 r936  
    207207  else {
    208208    if(givenCard != cardSol) {
    209       printf("CouenneRecordBestSol::setSol(): ### ERROR: givenCard: %d  cardSol: %d", givenCard, cardSol);
    210       exit(1);
     209      //printf("CouenneRecordBestSol::setSol(): ### ERROR: givenCard: %d  cardSol: %d", givenCard, cardSol);
     210      //exit(1);
     211
     212        double *newSol = new double [givenCard];
     213        CoinCopyN (givenSol, givenCard, newSol);
     214        delete [] modSol;
     215        modSol = newSol;
     216        cardSol = givenCard;
     217
    211218    }
    212219  }
     
    341348    else {
    342349      if(givenModCard != cardModSol) {
    343         printf("CouenneRecordBestSol::setModSol(): ### ERROR: givenModCard: %d  cardModSol: %d", givenModCard, cardModSol);
    344         exit(1);
     350        // printf("CouenneRecordBestSol::setModSol(): ### ERROR: givenModCard: %d  cardModSol: %d", givenModCard, cardModSol);
     351        // exit(1);
     352
     353        double *newModSol = new double [givenModCard];
     354        CoinCopyN (givenModSol, givenModCard, newModSol);
     355        delete [] modSol;
     356        modSol = newModSol;
     357        cardModSol = givenModCard;
    345358      }
    346359    }
  • trunk/Couenne/src/standardize/standardize.cpp

    r811 r936  
    476476#endif
    477477
    478   delete auxSet_;
     478  //  delete auxSet_;
    479479
    480480  // Create evaluation order ////////////////////////////////////////////////////
Note: See TracChangeset for help on using the changeset viewer.