# source:stable/0.2/Couenne/src/branch/doStrongBranching.cpp@159

Last change on this file since 159 was 159, checked in by pbelotti, 11 years ago

created new stable branch 0.2 from trunk (rev. 157)

File size: 11.9 KB
Line
1/* \$Id: doStrongBranching.cpp 141 2009-06-03 04:19:19Z pbelotti \$ */
2/*
3 * Name:    doStrongBranching.cpp
4 * Authors: Andreas Waechter, IBM Corp.
5 *          Pietro Belotti, CMU
6 * Purpose: actual strong branching method
7 *
8 * (C) Carnegie-Mellon University, 2008.
10 */
11
12#include "CoinTime.hpp"
13#include "CouenneChooseStrong.hpp"
14#include "CouenneProblem.hpp"
15#include "CouenneBranchingObject.hpp"
16#include "CouenneSolverInterface.hpp"
17
18/// compute Euclidean distance between two points (most likely LP solutions)
19/// l_2 norm by default, but can change it by fourth parameter
20double distance (const double *p1, const double *p2, int size, double k=2.) {
21
22  double
23    result = 0.,
24    element;
25
26  if (k == 2.) // a bit faster, probably
27
28    while (size--) {
29      element = *p1++ - *p2++;
30      result += element * element;
31    }
32
33  else
34
35    while (size--) {
36      element = *p1++ - *p2++;
37      result += pow (element, k);
38    }
39
40  return pow (result, 1. / k);
41}
42
43
44namespace Bonmin {
45
46  /**  This is a utility function which does strong branching on
47       a list of objects and stores the results in OsiHotInfo.objects.
48       On entry the object sequence is stored in the OsiHotInfo object
49       and maybe more.
50       It returns -
51
52      -1 - one branch was infeasible both ways
53       0 - all inspected - nothing can be fixed
54       1 - all inspected - some can be fixed (returnCriterion==0)
55       2 - may be returning early - one can be fixed (last one done) (returnCriterion==1)
56       3 - returning because max time
57  */
58  int CouenneChooseStrong::doStrongBranching (CouenneSolverInterface * solver,
59                                              OsiBranchingInformation *info,
60                                              int numberToDo, int returnCriterion)
61  {
62
63    jnlst_ -> Printf (J_ITERSUMMARY, J_BRANCHING,
64                      "\n-\n------- CCS: trying %d objects:\n", numberToDo);
65
66    solver -> doingResolve () = false; // turns off setCutoff and restoreUnused
67
68    int numberColumns = solver -> getNumCols ();
69
70    solver -> markHotStart (); // save current LP point
71
72    const double
73      *lower = info -> lower_,
74      *upper = info -> upper_;
75
76    double
77      *saveLower = CoinCopyOfArray (info -> lower_, numberColumns),
78      *saveUpper = CoinCopyOfArray (info -> upper_, numberColumns),
79      *Lower0    = CoinCopyOfArray (info -> lower_, numberColumns), // delete afterwards
80      *Upper0    = CoinCopyOfArray (info -> upper_, numberColumns),
81      *oldLower  = new double [numberColumns],
82      *oldUpper  = new double [numberColumns],
83      *lpSol     = NULL,
84       timeStart = CoinCpuTime ();
85
86    // LP solution for distance
87    if (pseudoUpdateLP_)
88      lpSol = CoinCopyOfArray (info -> solution_, numberColumns);
89
90    // provide Couenne problem with point/bounds contained in info
91    problem_ -> domain () -> push
92      (problem_ -> nVars (),
93       info -> solution_,
94       info -> lower_,
95       info -> upper_);
96
97    int returnCode = 0, iDo = 0;
98
99    for (;iDo < numberToDo; iDo++) {
100
101      HotInfo * result = results_ () + iDo; // retrieve i-th object to test
102
103      CouenneObject *CouObj = dynamic_cast <CouenneObject *>
104        (solver_ -> objects () [result -> whichObject ()]);
105
106      // For now just 2 way
107      OsiBranchingObject * branch = result -> branchingObject ();
108      assert (branch->numberBranches()==2);
109
110      CouenneBranchingObject *cb = dynamic_cast <CouenneBranchingObject *> (branch);
111      if (cb) cb -> setSimulate (true);
112
113      /* Try the first direction.  Each subsequent call to branch()
114         performs the specified branch and advances the branch object
115         state to the next branch alternative.) */
116
117      int
118        status0 = -1,
119        status1 = -1;
120
121      OsiSolverInterface * thisSolver = solver;
122
123      // DOWN DIRECTION ///////////////////////////////////////////////////////
124
125      if (branch -> boundBranch ()) { // a (variable) bound branch
126
127        if (branch -> branch (solver) > COUENNE_INFINITY) // branch is infeasible
128          result -> setDownStatus (status0 = 1);
129
130        else { // branch is feasible, solve and compare
131
132          solver -> solveFromHotStart ();
133          if (pseudoUpdateLP_ && CouObj && solver -> isProvenOptimal ()) {
134            CouNumber dist = distance (lpSol, solver -> getColSolution (), numberColumns);
135            if (dist > COUENNE_EPS)
136              CouObj -> setEstimate (dist, 0);
137          }
138        }
139
140      } else {                       // some more complex branch, have to clone solver
141
142        // adding cuts or something
143        thisSolver = solver -> clone ();
144
145        if (branch -> branch (thisSolver) > COUENNE_INFINITY)
146          result -> setDownStatus (status0 = 1);
147
148        else { // set hot start iterations
149          int limit;
150          thisSolver -> getIntParam (OsiMaxNumIterationHotStart, limit);
151          thisSolver -> setIntParam (OsiMaxNumIteration,         limit);
152
153          thisSolver -> resolve ();
154          if (pseudoUpdateLP_ && CouObj && thisSolver -> isProvenOptimal ()) {
155            CouNumber dist = distance (lpSol, thisSolver -> getColSolution (), numberColumns);
156            if (dist > COUENNE_EPS)
157              CouObj -> setEstimate (dist, 0);
158            //CouObj -> setEstimate (distance (lpSol, thisSolver->getColSolution (),numberColumns), 0);
159          }
160        }
161      }
162
163      // can check if we got solution
164      // status is 0 finished, 1 infeasible and 2 unfinished and 3 is solution
165
166      // only update information if this branch is feasible
167      if (status0 < 0)
168        status0 = result -> updateInformation (thisSolver, info, this);
169
170      numberStrongIterations_ += thisSolver -> getIterationCount();
171
172      if ((status0 == 3) && (trustStrongForSolution_)) {
173        // new solution already saved
174        info -> cutoff_ = goodObjectiveValue_;
175        problem_ -> setCutOff (goodObjectiveValue_);
176        status0 = 0;
177      }
178
179      if (solver != thisSolver)
180        delete thisSolver;
181
182      // save current bounds as tightened by the down branch; will be
183      // used below to update global bounding box in solver
184      CoinCopyN (problem_ -> Lb (), numberColumns, oldLower);
185      CoinCopyN (problem_ -> Ub (), numberColumns, oldUpper);
186
187      // Restore pre-left-branch bounds in solver
188      for (int j=0; j<numberColumns; j++) {
189
190        if (saveLower [j] != lower [j]) solver -> setColLower (j, saveLower [j]);
191        if (saveUpper [j] != upper [j]) solver -> setColUpper (j, saveUpper [j]);
192      }
193
194      // UP DIRECTION ///////////////////////////////////////////////////////
195
196      thisSolver = solver;
197
198      if (branch -> boundBranch ()) { // (variable) bound branch
199
200        if (branch -> branch (solver) > COUENNE_INFINITY)
201          result -> setUpStatus (status1 = 1);
202
203        else {
204          solver -> solveFromHotStart ();
205          if (pseudoUpdateLP_ && CouObj && solver -> isProvenOptimal ()) {
206            CouNumber dist = distance (lpSol, solver -> getColSolution (), numberColumns);
207            if (dist > COUENNE_EPS)
208              CouObj -> setEstimate (dist, 1);
209            //CouObj -> setEstimate (distance (lpSol, solver -> getColSolution (), numberColumns), 1);
210          }
211        }
212      } else {                     // some more complex branch, have to clone solver
213        // adding cuts or something
214        thisSolver = solver -> clone ();
215
216        if (branch -> branch (thisSolver) > COUENNE_INFINITY)
217          result -> setUpStatus (status1 = 1);
218
219        else {
220          // set hot start iterations
221          int limit;
222          thisSolver -> getIntParam (OsiMaxNumIterationHotStart, limit);
223          thisSolver -> setIntParam (OsiMaxNumIteration,         limit);
224
225          thisSolver -> resolve();
226          if (pseudoUpdateLP_ && CouObj && thisSolver -> isProvenOptimal ()) {
227            CouNumber dist = distance (lpSol, thisSolver -> getColSolution (), numberColumns);
228            if (dist > COUENNE_EPS)
229              CouObj -> setEstimate (dist, 1);
230            //CouObj -> setEstimate (distance (lpSol, thisSolver->getColSolution (),numberColumns), 1);
231          }
232        }
233      }
234
235      // can check if we got solution
236      // status is 0 finished, 1 infeasible and 2 unfinished and 3 is solution
237
238      // only update information if this branch is feasible
239      if (status1 < 0)
240        status1 = result -> updateInformation (thisSolver, info, this);
241
242      numberStrongDone_++;
243      numberStrongIterations_ += thisSolver->getIterationCount();
244
245      if ((status1 == 3) && (trustStrongForSolution_)) {
246        // new solution already saved
247        info -> cutoff_ = goodObjectiveValue_;
248        problem_ -> setCutOff (goodObjectiveValue_);
249        status1 = 0;
250      }
251
252      jnlst_ -> Printf (J_ITERSUMMARY, J_BRANCHING, "-------\n");
253
254      if (cb) cb -> setSimulate (false);
255
256      /////////////////////////////////////////////////////////////////////////////
257
258      if (solver != thisSolver)
259        delete thisSolver;
260
261      bool tightened = false;
262
263      t_chg_bounds *chg_bds = new t_chg_bounds [numberColumns];
264
265      // extend problem_'s bounding box to include downbranch's tightened
266      for (int j=0; j<numberColumns; j++) {
267
268        if (oldLower [j] < problem_ -> Lb (j)) problem_ -> Lb (j) = oldLower [j];
269        if (oldUpper [j] > problem_ -> Ub (j)) problem_ -> Ub (j) = oldUpper [j];
270
271        if (problem_ -> Lb (j) > lower [j] + COUENNE_EPS) {
272          chg_bds [j].setLower (t_chg_bounds::CHANGED);
273          tightened = true;
274        }
275
276        if (problem_ -> Ub (j) < upper [j] - COUENNE_EPS) {
277          chg_bds [j].setUpper (t_chg_bounds::CHANGED);
278          tightened = true;
279        }
280      }
281
282      if (tightened &&                     // have tighter bounds
283          (problem_ -> doFBBT ()) &&       // selected FBBT
284          !(problem_ -> btCore (chg_bds))) // tighten again on root
285
286        status0 = status1 = 1;             // if returns false, problem is infeasible
287
288      delete [] chg_bds;
289
290      // create union of bounding box from both branching directions
291      for (int j=0; j<numberColumns; j++) {
292
293        if (oldLower [j] < problem_ -> Lb (j)) problem_ -> Lb (j) = oldLower [j];
294        if (oldUpper [j] > problem_ -> Ub (j)) problem_ -> Ub (j) = oldUpper [j];
295      }
296
297      // set new bounding box as the possibly tightened one (a subset
298      // of the initial)
299      for (int j=0; j<numberColumns; j++) {
300
301        solver -> setColLower (j, saveLower [j] = problem_ -> Lb (j));
302        solver -> setColUpper (j, saveUpper [j] = problem_ -> Ub (j));
303      }
304
305      /*
306        End of evaluation for this candidate object. Possibilities are:
307
308        * Both sides below cutoff; this variable is a candidate for
309          branching.
310
311        * Both sides infeasible or above the objective cutoff: no
312          further action here. Break from the evaluation loop and
313          assume the node will be purged by the caller.
314
315        * One side feasible and below cutoff: Install the branch
316          (i.e., fix the variable). Possibly break from the evaluation
317          loop and assume the node will be reoptimised by the caller.
318      */
319
320      if (status0 == 1 &&
321          status1 == 1) { // infeasible
322        returnCode=-1;
323        break; // exit loop
324      } else if (status0==1 || status1==1) {
325        numberStrongFixed_++;
326        if (!returnCriterion) {
327          returnCode=1;
328        } else {
329          returnCode=2;
330          break;
331        }
332      }
333
334      bool hitMaxTime = ( CoinCpuTime()-timeStart > info->timeRemaining_);
335      if (hitMaxTime) {
336        returnCode=3;
337        break;
338      }
339    } // end loop /***********************************/
340
341
342    if (jnlst_ -> ProduceOutput (J_DETAILED, J_BRANCHING)) {
343      printf ("tightened bounds: ");
344      // create union of bounding box from both branching directions
345      for (int j=0; j<numberColumns; j++) {
346
347        if (problem_ -> Lb (j) > Lower0 [j]) printf ("l%d (%g-->%g) ", j,Lower0[j], problem_->Lb (j));
348        if (problem_ -> Ub (j) < Upper0 [j]) printf ("u%d (%g-->%g) ", j,Upper0[j], problem_->Ub (j));
349      }
350    }
351
352    delete [] Lower0;
353    delete [] Upper0;
354
355    problem_ -> domain () -> pop (); // discard current point/bounds from problem
356
357    delete [] lpSol;
358
359    jnlst_ -> Printf (J_ITERSUMMARY, J_BRANCHING, "----------------------done\n\n\n");
360
361    if (iDo < numberToDo) iDo++; // exited due to infeasibility
362    assert (iDo <= (int) results_.size());
363    results_.resize (iDo);
364
365    delete [] oldLower;
366    delete [] oldUpper;
367
368    delete [] saveLower;
369    delete [] saveUpper;
370
371    solver -> unmarkHotStart ();     // Delete the snapshot
372
373    solver -> doingResolve () = true;
374
375    return returnCode;
376  }
377}
Note: See TracBrowser for help on using the repository browser.