source: trunk/Couenne/src/branch/doStrongBranching.cpp @ 39

Last change on this file since 39 was 39, checked in by pbelotti, 12 years ago

Merging some of the latest changes from Bonmin/branches/Couenne: added Complementarity objects, disjunctive cuts, SOS (these don't work yet). Fixed some output. Improved cut-optimum checking. Sparser initial problem by replacing constraints ax-b>=0, w=ax-b with equivalent w>=0, w=ax-b. Fixed some bugs in new getBounds.

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