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