source: trunk/Couenne/src/bound_tightening/fake_tightening.cpp @ 39

Last change on this file since 39 was 39, checked in by pbelotti, 11 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: 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.
7 * This file is licensed under the Common Public License (CPL)
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
141#ifdef DEBUG
142    printf(" [%10g,%10g] lb = %g {fea=%d,btr=%d} ",
143           Lb (index), Ub (index), Lb (objind),feasible,betterbds);
144#endif
145
146    if (feasible && !betterbds) {
147
148      // case 1: too tight, move inner out
149      inner = fb;
150
151      // restore initial bound
152      CoinCopyN (chg_bds, ncols, f_chg);
153      CoinCopyN (olb, ncols, Lb ());
154      CoinCopyN (oub, ncols, Ub ());
155
156    } else {
157
158      // case 2: tightening valid, apply and move outer in
159
160#ifdef DEBUG
161      printf (" --> %cbound [x_%d]: %g --> %g",direction?'U':'L',index,(direction?oub:olb)[index],fb);
162      if (optimum_ && 
163          ((!direction &&
164            (optimum_ [index] >= olb [index]) && 
165            (optimum_ [index] <= Lb (index) - COUENNE_EPS)) ||
166           (direction &&
167            (optimum_ [index] <= oub [index]) && 
168            (optimum_ [index] >= Ub (index) + COUENNE_EPS)))) {
169        printf ("fake tightening cuts out optimum: x%d=%g in [%g,%g] but not in [%g,%g]\n",
170                index, olb [index], oub [index], Lb (index), Ub (index));
171      }
172#endif
173
174      /*bool do_not_tighten = false;
175
176      // check if cut best known solution
177      if (optimum_) {
178        if (direction) {
179          if ((oub [index] > optimum_ [index]) &&
180              (fb          < optimum_ [index])) {
181            printf ("aggressive bt cuts optimum ub %d: %g < %g < %g\n",
182                    index, fb, optimum_ [index], oub [index]);
183            do_not_tighten = true;
184          }
185        } else {
186          if ((olb [index] < optimum_ [index]) &&
187              (fb          > optimum_ [index])) {
188            printf ("aggressive bt cuts optimum lb %d: %g < %g < %g\n",
189                    index, olb [index], optimum_ [index], fb);
190            do_not_tighten = true;
191          }
192        }
193        }*/
194
195      //if (!do_not_tighten) {
196
197        // apply bound
198      if (direction) {
199        oub [index] = Ub (index) = intvar ? floor (fb) : fb; 
200        chg_bds [index]. setUpper (t_chg_bounds::CHANGED);
201      }
202      else {
203        olb [index] = Lb (index) = intvar ? ceil (fb) : fb; 
204        chg_bds [index].setLower (t_chg_bounds::CHANGED);
205      }
206
207      outer = fb; // we have at least a tightened bound, save it
208
209      tightened = true;
210        //}
211
212      // restore initial bound
213      CoinCopyN (chg_bds, ncols, f_chg);
214      CoinCopyN (olb, ncols, Lb ());
215      CoinCopyN (oub, ncols, Ub ());
216
217      //#if BR_TEST_LOG < 0 // for fair testing
218      // check tightened problem for feasibility
219      if (!(btCore (chg_bds))) {
220#ifdef DEBUG
221        printf ("\n    pruned by aggressive BT\n");
222#endif
223        return -1;
224      }
225      //#endif
226    }
227
228    // TODO: compute golden section
229    fb = (inner + outer) / 2;
230
231    //    if () fb = (          inner + (phi-1) * outer) / phi;
232    //    else  fb = ((phi-1) * inner +           outer) / phi;
233
234    //  if (!feasible)       
235    //    fb = fictBounds (direction, xcur,
236    //       direction ? lb [index] : outer,
237    //       direction ? outer      : ub [index]);
238    //    fb = fictBounds (direction, xcur, CoinMin (inner, outer), CoinMax (inner, outer));
239
240#ifdef DEBUG
241    printf ("\n");
242#endif
243  }
244
245
246  Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_BOUNDTIGHTENING, "\n");
247  if (tightened) 
248    Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_BOUNDTIGHTENING, 
249                    "  [x%2d] pruned %s [%g, %g] -- lb = %g cutoff = %g\n", 
250                    index,direction?"right":"left",
251                    olb[index],oub[index], Lb (objind), getCutOff ());
252
253  return tightened ? 1 : 0;
254}
Note: See TracBrowser for help on using the repository browser.