source: trunk/Couenne/src/branch/CouenneChooseStrong.cpp @ 577

Last change on this file since 577 was 577, checked in by pbelotti, 9 years ago

restoring FP for tests. Changing branching point on integer variables with integer value. Fixing large branching points for variables with large (upper/lower) bounds and large LP values. Added debug code to twoImplBounds for check against optimal solutions. Added fictitious objective function of 0 when none is defined. Redeclared size_t in invmap.cpp for compatibility with MSVS.

  • Property svn:keywords set to Author Date Id Revision
File size: 17.1 KB
Line 
1/* $Id: CouenneChooseStrong.cpp 577 2011-05-21 20:38:48Z pbelotti $
2 *
3 * Name:    CouenneChooseStrong.cpp
4 * Authors: Andreas Waechter, IBM Corp.
5 *          Pietro Belotti, Lehigh University
6 * Purpose: Strong branching objects for Couenne
7 *
8 * (C) Carnegie-Mellon University, 2006-10.
9 * This file is licensed under the Eclipse Public License (EPL)
10 */
11
12#include "CouenneObject.hpp"
13#include "BonChooseVariable.hpp"
14#include "CouenneChooseStrong.hpp"
15#include "CouenneProblem.hpp"
16#include "CouenneProblemElem.hpp"
17#include "CouenneBranchingObject.hpp"
18
19#include "CouenneRecordBestSol.hpp"
20
21//#define TRACE_STRONG
22
23using namespace Couenne;
24
25const CouNumber estProdEps = 1e-6;
26
27
28  /// constructor
29  CouenneChooseStrong::CouenneChooseStrong (Bonmin::BabSetupBase &b, CouenneProblem* p, JnlstPtr jnlst) :
30
31  Bonmin::BonChooseVariable (b, b.continuousSolver()),
32  problem_          (p),
33  jnlst_            (jnlst),
34  branchtime_       (0.) {
35
36    std::string s;
37
38    b.options () -> GetStringValue ("pseudocost_mult_lp", s, "couenne.");
39    pseudoUpdateLP_ = (s == "yes");     
40
41    b.options () -> GetStringValue ("estimate_select", s, "couenne.");
42    estimateProduct_ = (s == "product");
43
44    b.options () -> GetStringValue ("trust_strong", s, "couenne.");
45
46    // trust solution from strong branching to provide right LB
47
48    setTrustStrongForSolution (s == "yes");
49    setTrustStrongForBound    (s == "yes");
50  }
51
52  /// copy constructor
53  CouenneChooseStrong::CouenneChooseStrong (const CouenneChooseStrong& rhs) :
54    Bonmin::BonChooseVariable (rhs),
55    problem_          (rhs.problem_),
56    pseudoUpdateLP_   (rhs.pseudoUpdateLP_),
57    estimateProduct_  (rhs.estimateProduct_),
58    jnlst_            (rhs.jnlst_),
59    branchtime_       (rhs.branchtime_)
60  {}
61
62  /// destructor
63  CouenneChooseStrong::~CouenneChooseStrong()
64  {if (branchtime_ > 1e-9) jnlst_ -> Printf (J_ERROR, J_BRANCHING, "Strong branching: total time %g\n", branchtime_);}
65
66  /// cloning method
67  OsiChooseVariable *
68  CouenneChooseStrong::clone() const
69  {return new CouenneChooseStrong(*this);}
70
71  /// assignment operator
72  CouenneChooseStrong&
73  CouenneChooseStrong::operator=(const CouenneChooseStrong & rhs)
74  {
75    if (this != &rhs) {
76      Bonmin::BonChooseVariable::operator=(rhs);
77      problem_         = rhs.problem_;
78      pseudoUpdateLP_  = rhs.pseudoUpdateLP_;
79      estimateProduct_ = rhs.estimateProduct_;
80      jnlst_           = rhs.jnlst_;
81      branchtime_      = rhs.branchtime_;
82    }
83    return *this;
84  }
85
86
87  /* Choose a variable. Returns:
88
89    -1  Node is infeasible
90     0  Normal termination - we have a candidate
91     1  All look satisfied - no candidate
92     2  We can change the bound on a variable - but we also have a strong branching candidate
93     3  We can change the bound on a variable - but we have a non-strong branching candidate
94     4  We can change the bound on a variable - no other candidates
95
96     We can pick up branch from whichObject() and whichWay()
97     We can pick up a forced branch (can change bound) from whichForcedObject() and whichForcedWay()
98     If we have a solution then we can pick up from goodObjectiveValue() and goodSolution()
99  */
100  int CouenneChooseStrong::chooseVariable (OsiSolverInterface * solver,
101                                           OsiBranchingInformation *info,
102                                           bool fixVariables) {
103
104    /// Note: had to copy code from
105    /// BonChooseVariable::chooseVariable() in order to test product
106    /// thing
107
108    problem_ -> domain () -> push
109      (problem_ -> nVars (),
110       info -> solution_,
111       info -> lower_,
112       info -> upper_);
113
114    int retval;
115
116    //int retval = BonChooseVariable::chooseVariable (solver, info, fixVariables);
117
118    // COPY of Bonmin starts here ////////////////////////////////////////////
119
120    // We assume here that chooseVariable is called once at the very
121    // beginning with fixVariables set to true.  This is then the root
122    // node.
123
124    bool isRoot = isRootNode(info);
125    int numberStrong = numberStrong_;
126    if (isRoot) {
127      numberStrong = CoinMax(numberStrong_, numberStrongRoot_);
128    }
129    if (numberUnsatisfied_) {
130      const double* upTotalChange = pseudoCosts_.upTotalChange();
131      const double* downTotalChange = pseudoCosts_.downTotalChange();
132      const int* upNumber = pseudoCosts_.upNumber();
133      const int* downNumber = pseudoCosts_.downNumber();
134      int numberBeforeTrusted = pseudoCosts_.numberBeforeTrusted();
135      int numberLeft = CoinMin(numberStrong -numberStrongDone_,numberUnsatisfied_);
136      results_.clear();
137      int returnCode=0;
138      bestObjectIndex_ = -1;
139      bestWhichWay_ = -1;
140      firstForcedObjectIndex_ = -1;
141      firstForcedWhichWay_ =-1;
142      double bestTrusted=-COIN_DBL_MAX;
143      for (int i=0;i<numberLeft;i++) {
144        int iObject = list_[i];
145        if (numberBeforeTrusted == 0||
146            i < minNumberStrongBranch_ ||
147            (
148              only_pseudo_when_trusted_ && number_not_trusted_>0 ) ||
149              ( !isRoot && (upNumber[iObject]<numberBeforeTrusted ||
150                          downNumber[iObject]<numberBeforeTrusted ))||
151              ( isRoot && (!upNumber[iObject] && !downNumber[iObject])) ) {
152         
153#ifdef TRACE_STRONG
154          printf("Push object %d for strong branch\n", iObject);
155#endif
156
157          results_.push_back(Bonmin::HotInfo(solver, info,
158                                solver->objects(), iObject));
159        }
160        else {
161          const OsiObject * obj = solver->object(iObject);
162
163#ifdef TRACE_STRONG
164          printf("Use pseudo cost for object %d\n", iObject);
165#endif
166
167          double
168            upEstimate       = (upTotalChange[iObject]*obj->upEstimate())/upNumber[iObject],
169            downEstimate     = (downTotalChange[iObject]*obj->downEstimate())/downNumber[iObject],
170            MAXMIN_CRITERION = maxminCrit (info),
171            minVal, maxVal, value;
172
173          if (downEstimate < upEstimate) {
174            minVal = downEstimate;
175            maxVal = upEstimate;
176          } else {
177            minVal = upEstimate;
178            maxVal = downEstimate;
179          }
180
181          value = 
182            estimateProduct_ ? 
183            ((estProdEps + minVal) * maxVal) :
184            (       MAXMIN_CRITERION  * minVal + 
185             (1.0 - MAXMIN_CRITERION) * maxVal);
186
187          if (value > bestTrusted) {
188            bestObjectIndex_=iObject;
189            bestWhichWay_ = upEstimate>downEstimate ? 0 : 1;
190            bestTrusted = value;
191          }
192        }
193      }
194      int numberFixed=0;
195      if (results_.size() > 0) {
196
197        // do strong branching //////////////////////////////////////////////////
198        returnCode = doStrongBranching (solver, info, results_.size(), 1);
199
200        if (bb_log_level_>=3) {
201          const char* stat_msg[] = {"NOTDON", "FEAS", "INFEAS", "NOFINI"};
202          message(SB_HEADER)<<CoinMessageEol;
203          for (unsigned int i = 0; i< results_.size(); i++) {
204            double up_change = results_[i].upChange();
205            double down_change = results_[i].downChange();
206            int up_status = results_[i].upStatus();
207            int down_status = results_[i].downStatus();
208            message(SB_RES)<<(int) i<<stat_msg[down_status+1]<<down_change
209            <<stat_msg[up_status+1]<< up_change<< CoinMessageEol;
210          }
211        }
212        if (returnCode >= 0 && 
213            returnCode <= 2) {
214          if (returnCode)
215            returnCode = (bestObjectIndex_>=0) ? 3 : 4;
216          for (unsigned int i=0;i < results_.size();i++) {
217            int iObject = results_[i].whichObject();
218
219            // UP ESTIMATE ////////////////////////////////////////
220
221            double upEstimate;
222            if (results_[i].upStatus()!=1) {
223              assert (results_[i].upStatus()>=0);
224              upEstimate = results_[i].upChange();
225            }
226            else {
227              // infeasible - just say expensive
228              if (info->cutoff_<1.0e50)
229                upEstimate = 2.0*(info->cutoff_-info->objectiveValue_);
230              else
231                upEstimate = 2.0*fabs(info->objectiveValue_);
232              if (firstForcedObjectIndex_ <0) {
233                // first fixed variable
234                firstForcedObjectIndex_ = iObject;
235                firstForcedWhichWay_ =0;
236              }
237              numberFixed++;
238              if (fixVariables) {
239                const OsiObject * obj = solver->object(iObject);
240                OsiBranchingObject * branch = obj->createBranch(solver,info,0);
241                branch->branch(solver);
242                delete branch;
243              }
244            }
245
246            // DOWN ESTIMATE /////////////////////////////////////
247
248            double downEstimate;
249            if (results_[i].downStatus()!=1) {
250              assert (results_[i].downStatus()>=0);
251              downEstimate = results_[i].downChange();
252            }
253            else {
254              // infeasible - just say expensive
255              if (info->cutoff_<1.0e50)
256                downEstimate = 2.0*(info->cutoff_-info->objectiveValue_);
257              else
258                downEstimate = 2.0*fabs(info->objectiveValue_);
259              if (firstForcedObjectIndex_ <0) {
260                firstForcedObjectIndex_ = iObject;
261                firstForcedWhichWay_ =1;
262              }
263              numberFixed++;
264              if (fixVariables) {
265                const OsiObject * obj = solver->object(iObject);
266                OsiBranchingObject * branch = obj->createBranch(solver,info,1);
267                branch->branch(solver);
268                delete branch;
269              }
270            }
271
272            double
273              MAXMIN_CRITERION = maxminCrit(info),           
274              minVal, maxVal, value;
275
276            if (downEstimate < upEstimate) {
277              minVal = downEstimate;
278              maxVal = upEstimate;
279            } else {
280              minVal = upEstimate;
281              maxVal = downEstimate;
282            }
283
284            value = 
285              estimateProduct_ ? 
286              ((estProdEps + minVal) * maxVal) :
287              (       MAXMIN_CRITERION  * minVal + 
288               (1.0 - MAXMIN_CRITERION) * maxVal);
289
290            if (value>bestTrusted) {
291              bestTrusted = value;
292              bestObjectIndex_ = iObject;
293              bestWhichWay_ = upEstimate>downEstimate ? 0 : 1;
294              // but override if there is a preferred way
295              const OsiObject * obj = solver->object(iObject);
296              if (obj->preferredWay()>=0&&obj->infeasibility())
297                bestWhichWay_ = obj->preferredWay();
298              if (returnCode)
299                returnCode = 2;
300            }
301          }
302        }
303        else if (returnCode==3) {
304          // max time - just choose one
305          bestObjectIndex_ = list_[0];
306          bestWhichWay_ = 0;
307          returnCode=0;
308        }
309      }
310      else {
311        bestObjectIndex_=list_[0];
312      }
313
314      if ( bestObjectIndex_ >=0 ) {
315        OsiObject * obj = solver->objects()[bestObjectIndex_];
316        obj->setWhichWay(bestWhichWay_);
317        message(BRANCH_VAR)<<obj->columnNumber()<< bestWhichWay_ <<CoinMessageEol;
318      }
319
320      message(CHOSEN_VAR)<<bestObjectIndex_<<CoinMessageEol;
321
322      if (numberFixed==numberUnsatisfied_&&numberFixed)
323        returnCode=4;
324      retval = returnCode;
325    }
326    else {
327      retval = 1;
328    }
329
330    // COPY of Bonmin ends here //////////////////////////////////////////////
331
332    problem_ -> domain () -> pop ();
333
334    return retval;
335  }
336
337
338  // Sets up strong list and clears all if initialize is true.
339  // Returns number of infeasibilities.
340  int CouenneChooseStrong::setupList (OsiBranchingInformation *info, bool initialize) {
341
342    static bool warned = false;
343
344    initialize = true; // to avoid failed assert in BonChooseVariable::setupList()
345
346    problem_ -> domain () -> push
347      (problem_ -> nVars (),
348       info -> solution_, 
349       info -> lower_, 
350       info -> upper_); // have to alloc+copy
351
352#ifdef COIN_HAS_NTY
353
354    if (problem_ -> orbitalBranching ()) {
355
356      problem_ -> ChangeBounds (info -> lower_,
357                                info -> upper_, 
358                                problem_ -> nVars ());
359   
360      problem_ -> Compute_Symmetry();
361    }
362
363#endif
364
365    jnlst_ -> Printf (J_ITERSUMMARY, J_BRANCHING, 
366                      "----------------- (strong) setup list\n");
367
368    if (jnlst_ -> ProduceOutput (J_DETAILED, J_BRANCHING)) {
369      for (int i=0; i<problem_ -> domain () -> current () -> Dimension (); i++)
370        printf ("%4d %20.4g [%20.4g %20.4g]\n", i,
371                info -> solution_ [i], info -> lower_ [i], info -> upper_ [i]);
372    }
373
374#ifdef TRACE_STRONG
375    printObjViol(info);
376#endif
377
378    // call Bonmin's setuplist
379    int retval = Bonmin::BonChooseVariable::setupList (info, initialize);
380
381    if (!retval) { // No branching is possible
382
383      double ckObj = info->objectiveValue_;
384
385      if (!(problem_ -> checkNLP (info -> solution_, ckObj, true))) {
386
387        if (!warned) {
388          printf ("CouenneChooseStrong::setupList() -- Warning: checkNLP returns infeasible, no branching object selected\n");
389          warned = true;
390        }
391        //exit (1);
392      }
393    }
394
395    jnlst_ -> Printf (J_ITERSUMMARY, J_BRANCHING, 
396                      "----------------- (strong) setup list done - %d infeasibilities\n", retval);
397
398
399    if(retval == 0) { // No branching is possible
400      double ckObj = 0;
401
402#ifdef FM_CHECKNLP2
403      if(!(problem_->checkNLP2(info->solution_, 
404                               info->objectiveValue_, true, // care about obj
405                               false, // do not stop at first viol
406                               true, // checkAll
407                               problem_->getFeasTol()))) {
408                                // false for NOT stopping at first violation
409        printf("CouenneChooseStrong::setupList(): ### WARNING: checkNLP2() returns infeasible, no branching object selected\n");
410      }
411#else /* not FM_CHECKNLP2 */
412      if(!(problem_->checkNLP(info->solution_, ckObj, true))) {
413        printf("CouenneChooseStrong::setupList(): ### WARNING: checkNLP() returns infeasible, no branching object selected\n");
414      }
415#endif /* not FM_CHECKNLP2 */
416       
417#ifdef FM_TRACE_OPTSOL
418#ifdef FM_CHECKNLP2
419      problem_->getRecordBestSol()->update();
420#else /* not FM_CHECKNLP2 */
421      problem_->getRecordBestSol()->update(info->solution_, problem_->nVars(),
422                                           ckObj, problem_->getFeasTol());
423#endif /* not FM_CHECKNLP2 */
424#endif
425    }
426
427#ifdef TRACE_STRONG
428    printf("Strong list: (obj_ind var_ind priority useful)\n");
429    printf("numberStrong: %d  numberStrongRoot: %d  retval: %d\n", 
430           numberStrong_, numberStrongRoot_, retval);
431
432    OsiObject ** object = info->solver_->objects();
433
434    for(int i=0; i<retval; i++) {
435      CouenneObject *co =  dynamic_cast <CouenneObject *> (object[list_[i]]);
436      int objectInd = -1;
437      if(co) {
438        objectInd = co->Reference()->Index();
439      }
440      else {
441        objectInd = object[list_[i]]->columnNumber();
442      }
443      printf(" (%d %d %d %6.4f)", list_[i], objectInd, 
444             object[list_[i]]->priority(), useful_[i]);
445    }
446    printf("\n");
447#endif
448
449
450    problem_ -> domain () -> pop ();
451    return retval;
452  }
453
454
455  /// Add list of options to be read from file ////////////////////////////////////////
456  void CouenneChooseStrong::registerOptions (Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions) {
457
458    roptions -> AddStringOption6
459      ("pseudocost_mult",
460       "Multipliers of pseudocosts for estimating and update estimation of bound",
461       "interval_br_rev",
462
463       "infeasibility", "infeasibility returned by object",
464
465       "projectDist",   "distance between current LP point and resulting branches' LP points",
466
467       "interval_lp",   "width of the interval between bound and current lp point",
468       "interval_lp_rev",   "similar to interval_lp, reversed",
469
470       "interval_br",   "width of the interval between bound and branching point",
471       "interval_br_rev",   "similar to interval_br, reversed");
472
473    roptions -> AddStringOption2
474      ("pseudocost_mult_lp",
475       "Use distance between LP points to update multipliers of pseudocosts " 
476       "after simulating branching",
477       "no",
478       "yes", "",
479       "no",  "");
480
481    roptions -> AddStringOption2
482      ("estimate_select",
483       "How the min/max estimates of the subproblems' bounds are used in strong branching",
484       "normal",
485       "normal",   "as usual in literature",
486       "product",  "use their product");
487
488    roptions -> AddStringOption2
489      ("trust_strong",
490       "Fathom strong branching LPs when their bound is above the cutoff",
491       "yes",
492       "yes", "",
493       "no",  "");
494  }
495
496
497  // Returns true if solution looks feasible against given objects
498  bool CouenneChooseStrong::feasibleSolution (const OsiBranchingInformation * info,
499                                              const double * solution,
500                                              int numberObjects,
501                                              const OsiObject ** objects) {
502
503    double obj = solution [problem_ -> Obj (0) -> Body () -> Index ()];
504
505#ifdef FM_CHECKNLP2
506    bool isFeas = problem_->checkNLP2(solution, 0, false, true, true, 
507                                      problem_->getFeasTol());
508    return isFeas;
509#else
510    return problem_ -> checkNLP (solution, obj);
511#endif
512  }
513/****************************************************************************/
514  void CouenneChooseStrong::printObjViol(OsiBranchingInformation *info) {
515
516    OsiObject ** object = info->solver_->objects();
517    int numberObjects = info->solver_->numberObjects();
518
519    printf("CouenneChooseStrong::printObjViol(): Object violations: (obj_ind  var_ind  violation)");
520    double maxViol = 0;
521    double minPosViol = 1e50;
522    for(int i=0; i<numberObjects; i++) {
523      int way;
524      double value = object[i]->infeasibility(info,way);
525      int indVar = object[i]->columnNumber();
526      maxViol = (value > maxViol ? value : maxViol);
527      if(value > 0.0) {
528        printf("(%d %d %f)", i, indVar, value);
529        minPosViol = (value < minPosViol ? value : minPosViol);
530      }
531    }
532    printf("\nmaxViol: %g  minPosViol: %g\n", maxViol, minPosViol);
533
534  } /* printObjViol */
535
536
537//}
Note: See TracBrowser for help on using the repository browser.