# source:trunk/Couenne/src/bound_tightening/fake_tightening.cpp@131

Last change on this file since 131 was 131, checked in by pbelotti, 13 years ago

looser checks on bound tightening -- should prevent good upper/lower bound on objective from fathoming a node

File size: 7.7 KB
Line
1/*
2 * Name:    fake_tightening.cpp
3 * Author:  Pietro Belotti
4 * Purpose: fake single bounds in variables to exclude parts of the solution space
5 *
6 * (C) Carnegie-Mellon University, 2007.
8 */
9
10#include "CouenneProblem.hpp"
11#include "BonBabInfos.hpp"
12
13//#define DEBUG
14
15#define MAX_ITER  2 // max # fake tightening (inner) iterations
16#define AGGR_MUL  2 // the larger,  the more conservative. Must be > 0
17#define AGGR_DIV  2 // the smaller, the more conservative. Must be > 1
18
19// golden ratio, used to find the ideal bound
20const CouNumber phi = 0.5 * (1. + sqrt (5.));
21
22// create fictitious bounds to tighten current interval
23CouNumber fictBounds (char direction,
24                      CouNumber  x,
25                      CouNumber  lb,
26                      CouNumber  ub) {
27
28  if   (lb < -COUENNE_INFINITY / 10) {
29    if (ub >  COUENNE_INFINITY / 10) { // ]-inf,+inf[
30
31      if (fabs (x) < COUENNE_EPS) return (direction ? AGGR_MUL : - AGGR_MUL);
32      else                        return (direction ? AGGR_MUL : - AGGR_MUL) * fabs (x);
33
34    } else { // ]-inf,u]
35
36      if      (x < -COUENNE_EPS) return (direction ? CoinMin (0., (x+ub)/2) : AGGR_MUL * x);
37      else if (x >  COUENNE_EPS) return (direction ? (x + (ub-x)/AGGR_DIV) : 0);
38      else                       return (direction ? (ub/AGGR_DIV) : -AGGR_MUL);
39    }
40  }
41  else {
42    if (ub >  COUENNE_INFINITY / 10) { // [l,+inf[
43
44      if      (x < -COUENNE_EPS) return (direction ? 0 : (x - (x-lb) / AGGR_DIV));
45      else if (x >  COUENNE_EPS) return (direction ? (AGGR_MUL * x) : CoinMax (0.,(x+lb)/2));
46      else                       return (direction ? AGGR_MUL : lb/AGGR_DIV);
47    } else // [l,u]
48      return (direction ? (x + (ub-x) / AGGR_DIV) : x - (x-lb) / AGGR_DIV);
49  }
50}
51
52
53// Single fake tightening. Return
54//
55// -1   if infeasible
56//  0   if no improvement
57// +1   if improved
58int CouenneProblem::
59fake_tighten (char direction,  // 0: left, 1: right
60              int index,       // index of the variable tested
61              const double *X, // point round which tightening is done
62              CouNumber *olb,  // cur. lower bound
63              CouNumber *oub,  //      upper
64              t_chg_bounds *chg_bds,
65              t_chg_bounds *f_chg) const {
66  int
67    ncols    = nVars (),
68    objind   = Obj (0) -> Body  () -> Index ();
69
70  assert (objind >= 0);
71
72  bool
73    tightened = false,
74    intvar    = variables_ [index] -> isInteger ();
75
76  CouNumber
77    xcur      = X [index],
78    inner     = xcur,                                                 // point closest to current x
79    outer     = (direction ? oub : olb) [index],                      // point closest to bound
80    fb        = fictBounds (direction, xcur, Lb (index), Ub (index)); // starting point
81
82  // This is a one-dimensional optimization problem between inner and
83  // outer, on a monotone function of which we can compute the value
84  // (with relative expense) but not the derivative.
85
86#ifdef DEBUG
87  CouNumber curdb     = Lb (objind);// : Ub (objind);  // current dual bound
88  printf ("  x_%d.  x = %10g, lb = %g, cutoff = %g-----------------\n", index,xcur,curdb,getCutOff());
89#endif
90
91  for (int iter = 0; iter < MAX_ITER; iter++) {
92
93    if (intvar) {
94
95      if (!direction) {inner = floor (inner); outer = ceil  (outer);}
96      else            {inner = ceil  (inner); outer = floor (outer);}
97
98      if ( direction && (inner > outer) ||
99          !direction && (inner < outer)) {
100
101        // apply bound
102        if (direction) {oub[index] = Ub (index) = fb; chg_bds[index].setUpper(t_chg_bounds::CHANGED);}
103        else           {olb[index] = Lb (index) = fb; chg_bds[index].setLower(t_chg_bounds::CHANGED);}
104
105        tightened = true;
106
107        // restore initial bound
108        CoinCopyN (chg_bds, ncols, f_chg);
109        CoinCopyN (olb, ncols, Lb ());
110        CoinCopyN (oub, ncols, Ub ());
111
112        break;
113      }
114
115      if (direction  && ((fb < inner) || (fb > outer)) ||
116          !direction && ((fb > inner) || (fb < outer)))
117        fb = 0.5 * (inner + outer);
118    }
119
120    if (direction) {
121      Lb (index) = intvar ? ceil (fb) : fb;
122      f_chg [index].setLower (t_chg_bounds::CHANGED);
123    } else {
124      Ub (index) = intvar ? floor (fb) : fb;
125      f_chg [index].setUpper (t_chg_bounds::CHANGED);
126    }
127
128    //    (direction ? lb_ : ub_) [index] = fb;
129
130#ifdef DEBUG
131    char c1 = direction ? '-' : '>', c2 = direction ? '<' : '-';
132    printf ("    #%3d: [%+10g -%c %+10g %c- %+10g] /\\/\\ ",iter,olb[index],c1,fb,c2, oub [index]);
133    printf (" [%10g,%10g]<%g,%g>=> ",Lb (index),Ub (index),CoinMin(inner,outer),CoinMax(inner,outer));
134#endif
135
136    bool
137      feasible  = btCore (f_chg),             // true if feasible with fake bound
138      betterbds = Lb (objind) > getCutOff (); // true if over cutoff
139
140#ifdef DEBUG
141    printf(" [%10g,%10g] lb = %g {fea=%d,btr=%d} ",
142           Lb (index), Ub (index), Lb (objind),feasible,betterbds);
143#endif
144
145    if (feasible && !betterbds) {
146
147      // case 1: too tight, move inner out
148      inner = fb;
149
150      // restore initial bound
151      CoinCopyN (chg_bds, ncols, f_chg);
152      CoinCopyN (olb, ncols, Lb ());
153      CoinCopyN (oub, ncols, Ub ());
154
155    } else {
156
157      // case 2: tightening valid, apply and move outer in
158
159#ifdef DEBUG
160      printf (" --> %cbound [x_%d]: %g --> %g",direction?'U':'L',index,(direction?oub:olb)[index],fb);
161      if (optimum_ &&
162          ((!direction &&
163            (optimum_ [index] >= olb [index]) &&
164            (optimum_ [index] <= Lb (index) - COUENNE_EPS)) ||
165           (direction &&
166            (optimum_ [index] <= oub [index]) &&
167            (optimum_ [index] >= Ub (index) + COUENNE_EPS)))) {
168        printf ("fake tightening cuts out optimum: x%d=%g in [%g,%g] but not in [%g,%g]\n",
169                index, olb [index], oub [index], Lb (index), Ub (index));
170      }
171#endif
172
173      /*bool do_not_tighten = false;
174
175      // check if cut best known solution
176      if (optimum_) {
177        if (direction) {
178          if ((oub [index] > optimum_ [index]) &&
179              (fb          < optimum_ [index])) {
180            printf ("aggressive bt cuts optimum ub %d: %g < %g < %g\n",
181                    index, fb, optimum_ [index], oub [index]);
182            do_not_tighten = true;
183          }
184        } else {
185          if ((olb [index] < optimum_ [index]) &&
186              (fb          > optimum_ [index])) {
187            printf ("aggressive bt cuts optimum lb %d: %g < %g < %g\n",
188                    index, olb [index], optimum_ [index], fb);
189            do_not_tighten = true;
190          }
191        }
192        }*/
193
194      //if (!do_not_tighten) {
195
196        // apply bound
197      if (direction) {
198        oub [index] = Ub (index) = intvar ? floor (fb) : fb;
199        chg_bds [index]. setUpper (t_chg_bounds::CHANGED);
200      }
201      else {
202        olb [index] = Lb (index) = intvar ? ceil (fb) : fb;
203        chg_bds [index].setLower (t_chg_bounds::CHANGED);
204      }
205
206      outer = fb; // we have at least a tightened bound, save it
207
208      tightened = true;
209        //}
210
211      // restore initial bound
212      CoinCopyN (chg_bds, ncols, f_chg);
213      CoinCopyN (olb, ncols, Lb ());
214      CoinCopyN (oub, ncols, Ub ());
215
216      //#if BR_TEST_LOG < 0 // for fair testing
217      // check tightened problem for feasibility
218      if (!(btCore (chg_bds))) {
219#ifdef DEBUG
220        printf ("\n    pruned by aggressive BT\n");
221#endif
222        return -1;
223      }
224      //#endif
225    }
226
227    // TODO: compute golden section
228    fb = (inner + outer) / 2;
229
230    //    if () fb = (          inner + (phi-1) * outer) / phi;
231    //    else  fb = ((phi-1) * inner +           outer) / phi;
232
233    //  if (!feasible)
234    //    fb = fictBounds (direction, xcur,
235    //       direction ? lb [index] : outer,
236    //       direction ? outer      : ub [index]);
237    //    fb = fictBounds (direction, xcur, CoinMin (inner, outer), CoinMax (inner, outer));
238
239#ifdef DEBUG
240    printf ("\n");
241#endif
242  }
243
244  Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_BOUNDTIGHTENING, "\n");
245  if (tightened)
246    Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_BOUNDTIGHTENING,
247                    "  [x%2d] pruned %s [%g, %g] -- lb = %g cutoff = %g\n",
248                    index,direction?"right":"left",
249                    olb[index],oub[index], Lb (objind), getCutOff ());
250
251  return tightened ? 1 : 0;
252}
Note: See TracBrowser for help on using the repository browser.