source: stable/0.4/Couenne/src/bound_tightening/fake_tightening.cpp @ 830

Last change on this file since 830 was 830, checked in by pbelotti, 8 years ago

merging change 828 from trunk (avoid probing with ridiculous bounds)

  • Property svn:keywords set to Author Date Id Revision
File size: 10.2 KB
Line 
1/* $Id: fake_tightening.cpp 830 2012-02-10 18:02:11Z 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-10.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11#include "CoinHelperFunctions.hpp"
12
13#include "CouenneProblem.hpp"
14#include "CouenneProblemElem.hpp"
15#include "CouenneExprVar.hpp"
16
17#include "BonBabInfos.hpp"
18
19using namespace Couenne;
20
21#define MAX_ITER  3 // max # fake tightening (inner) iterations
22#define AGGR_MUL  2 // the larger,  the more conservative. Must be > 0
23#define AGGR_DIV  2 // the smaller, the more conservative. Must be > 1
24
25// golden ratio, used to find the ideal bound
26const CouNumber phi = 0.5 * (1. + sqrt (5.));
27
28// create fictitious bounds to tighten current interval
29CouNumber fictBounds (char direction,
30                      CouNumber  x,
31                      CouNumber  lb,   
32                      CouNumber  ub) {
33
34#define LARGE_BOUND 1e10
35
36  if   (lb < -LARGE_BOUND) {
37    if (ub >  LARGE_BOUND) { // ]-inf,+inf[
38
39      return (!direction ? -LARGE_BOUND / AGGR_DIV : LARGE_BOUND / AGGR_DIV);
40
41      //return (!direction ? -sqrt (-lb) : sqrt (ub));
42
43      //if (fabs (x) < COUENNE_EPS) return (direction ? AGGR_MUL : - AGGR_MUL);
44      //else                        return (direction ? AGGR_MUL : - AGGR_MUL) * fabs (x);
45
46    } else { // ]-inf,u]
47
48      if (!direction)
49        return -LARGE_BOUND / AGGR_DIV; //-sqrt (-lb); // try to tighten interval from a very large value
50
51      if      (x < -COUENNE_EPS) return (CoinMin (0., (x+ub)/2));
52      else if (x >  COUENNE_EPS) return ((x + (ub-x)/AGGR_DIV));
53      else                       return ((ub/AGGR_DIV));
54
55      // if      (x < -COUENNE_EPS) return (direction ? CoinMin (0., (x+ub)/2) : AGGR_MUL * x);
56      // else if (x >  COUENNE_EPS) return (direction ? (x + (ub-x)/AGGR_DIV) : 0);
57      // else                       return (direction ? (ub/AGGR_DIV) : -AGGR_MUL);
58    }
59  }
60  else {
61    if (ub >  LARGE_BOUND) { // [l,+inf[
62
63      if (direction)
64        return LARGE_BOUND / AGGR_DIV; //sqrt (ub);
65
66      if      (x < -COUENNE_EPS) return ((x - (x-lb) / AGGR_DIV));
67      else if (x >  COUENNE_EPS) return (CoinMax (0.,(x+lb)/2));
68      else                       return (lb/AGGR_DIV);
69
70      // if      (x < -COUENNE_EPS) return (direction ? 0 : (x - (x-lb) / AGGR_DIV));
71      // else if (x >  COUENNE_EPS) return (direction ? (AGGR_MUL * x) : CoinMax (0.,(x+lb)/2));
72      // else                       return (direction ? AGGR_MUL : lb/AGGR_DIV);
73
74    } else // [l,u]
75      return (direction ? 
76              (x + (ub-x) / AGGR_DIV) : 
77              (x - (x-lb) / AGGR_DIV));
78  }
79}
80
81
82// Single fake tightening. Return
83//
84// -1   if infeasible
85//  0   if no improvement
86// +1   if improved
87int CouenneProblem::
88fake_tighten (char direction,  // 0: left, 1: right
89              int index,       // index of the variable tested
90              const double *X, // point round which tightening is done
91              CouNumber *olb,  // cur. lower bound
92              CouNumber *oub,  //      upper
93              t_chg_bounds *chg_bds,
94              t_chg_bounds *f_chg) const {
95
96  int
97    ncols    = nVars (),
98    objind   = Obj (0) -> Body  () -> Index ();
99
100  //assert (objind >= 0);
101
102  bool 
103    tightened = false,
104    intvar    = variables_ [index] -> isInteger ();
105
106  CouNumber
107    xcur      = X [index],
108    inner     = xcur,                                                 // point closest to current x
109    outer     = (direction ? oub : olb) [index],                      // point closest to bound
110    fb        = fictBounds (direction, xcur, Lb (index), Ub (index)); // starting point
111
112  // This is a one-dimensional optimization problem between inner and
113  // outer, on a monotone function of which we can compute the value
114  // (with relative expense) but not the derivative.
115
116  jnlst_ -> Printf (Ipopt::J_ERROR, J_BOUNDTIGHTENING, "  x_%d.  x = %10g, lb = %g, cutoff = %g-----------------\n", index,xcur,objind >= 0 ? Lb (objind) : 0., getCutOff());
117
118  /*if (index == objind)
119    printf ("  x_%d [%g,%g].  x = %10g, break at %g, cutoff = %g-----------------\n",
120    index, Lb (index), Ub (index), xcur, fb, getCutOff());*/
121
122  for (int iter = 0; iter < MAX_ITER; iter++) {
123
124    if (intvar) {
125
126      if (!direction) {inner = floor (inner + COUENNE_EPS); outer = ceil  (outer - COUENNE_EPS);}
127      else            {inner = ceil  (inner - COUENNE_EPS); outer = floor (outer + COUENNE_EPS);}
128
129      if ( (direction && (inner > outer + .5)) || // more robust check on integer-valued doubles
130          (!direction && (inner < outer - .5))) {
131
132        // fictitious interval is empty, hence useless to check.
133
134        // apply new (valid, tightened) bound
135        if (direction) {oub[index] = Ub (index) = fb; chg_bds[index].setUpper(t_chg_bounds::CHANGED);}
136        else           {olb[index] = Lb (index) = fb; chg_bds[index].setLower(t_chg_bounds::CHANGED);}
137
138        tightened = true;
139
140        if (!(btCore (f_chg))) 
141          return -1;
142
143        CoinCopyN (Lb (), ncols, olb);
144        CoinCopyN (Ub (), ncols, oub);
145
146        // restore initial bound.
147        CoinCopyN (chg_bds, ncols, f_chg);
148        //CoinCopyN (olb, ncols, Lb ());
149        //CoinCopyN (oub, ncols, Ub ());
150
151        break;
152      }
153
154      if ( (direction && ((fb < inner) || (fb > outer))) ||
155          (!direction && ((fb > inner) || (fb < outer))))
156        fb = 0.5 * (inner + outer);
157    }
158
159    if (direction) {
160      Lb (index) = intvar ? ceil (fb - COUENNE_EPS)  : fb; 
161      f_chg [index].setLower (t_chg_bounds::CHANGED);
162    } else {
163      Ub (index) = intvar ? floor (fb + COUENNE_EPS) : fb; 
164      f_chg [index].setUpper (t_chg_bounds::CHANGED);
165    }
166
167    //    (direction ? lb_ : ub_) [index] = fb;
168
169    if (jnlst_ -> ProduceOutput (Ipopt::J_ERROR, J_BOUNDTIGHTENING)) {
170      char c1 = direction ? '-' : '>', c2 = direction ? '<' : '-';
171      printf ("    # x%d = %g iter %3d: [%+10g -%c %+10g %c- %+10g] /\\/\\ ",index,xcur,iter,olb[index],c1,fb,c2, oub [index]);
172      printf (" [%10g,%10g]<%g,%g>=> ",Lb (index),Ub (index),CoinMin(inner,outer),CoinMax(inner,outer));
173    }
174
175    bool
176      feasible  = btCore (f_chg),                                            // true if feasible with fake bound
177      betterbds = objind >= 0 && (Lb (objind) > getCutOff () + COUENNE_EPS); // true if over cutoff
178
179    jnlst_ -> Printf (Ipopt::J_ERROR, J_BOUNDTIGHTENING, " [%10g,%10g] lb = %g {fea=%d,btr=%d} ", Lb (index), Ub (index), objind >= 0 ? Lb (objind) : 0., feasible,betterbds);
180
181    if (feasible && !betterbds) {
182
183      // case 1: too tight, move inner out
184      inner = fb;
185
186      // restore initial bound
187      CoinCopyN (chg_bds, ncols, f_chg);
188      CoinCopyN (olb, ncols, Lb ());
189      CoinCopyN (oub, ncols, Ub ());
190
191    } else {
192
193      // Here, !feasible || betterbds
194      //
195      // If !feasible OR
196      //    (betterbds and the new lb is above the cutoff)
197      //
198      // then there is a tightening
199
200      // case 2: tightening valid, apply and move outer in
201
202      //printf (" --> %cbound [x_%d]: %g --> %g",direction?'U':'L',index,(direction?oub:olb)[index],fb);
203
204      if (optimum_ && 
205          ((!direction &&
206            (optimum_ [index] >= olb [index]) && 
207            (optimum_ [index] <= Lb (index) - COUENNE_EPS)) ||
208           (direction &&
209            (optimum_ [index] <= oub [index]) && 
210            (optimum_ [index] >= Ub (index) + COUENNE_EPS)))) {
211
212        jnlst_ -> Printf (Ipopt::J_ERROR, J_BOUNDTIGHTENING, 
213                          "fake tightening CUTS optimum: x%d=%g in [%g,%g] but not in [%g,%g]\n",
214                          index, olb [index], oub [index], Lb (index), Ub (index));
215      }
216
217      /*bool do_not_tighten = false;
218
219      // check if cut best known solution
220      if (optimum_) {
221        if (direction) {
222          if ((oub [index] > optimum_ [index]) &&
223              (fb          < optimum_ [index])) {
224            printf ("aggressive bt cuts optimum ub %d: %g < %g < %g\n",
225                    index, fb, optimum_ [index], oub [index]);
226            do_not_tighten = true;
227          }
228        } else {
229          if ((olb [index] < optimum_ [index]) &&
230              (fb          > optimum_ [index])) {
231            printf ("aggressive bt cuts optimum lb %d: %g < %g < %g\n",
232                    index, olb [index], optimum_ [index], fb);
233            do_not_tighten = true;
234          }
235        }
236        }*/
237
238      //if (!do_not_tighten) {
239
240      // apply bound
241      if (direction) {
242
243        oub [index] = Ub (index) = intvar ? floor (fb + COUENNE_EPS) : fb; 
244        chg_bds [index]. setUpper (t_chg_bounds::CHANGED);
245
246      } else {
247
248        olb [index] = Lb (index) = intvar ? ceil  (fb - COUENNE_EPS) : fb; 
249        chg_bds [index]. setLower (t_chg_bounds::CHANGED);
250      }
251
252      outer = fb; // we have at least a tightened bound, save it
253
254      tightened = true;
255      //}
256
257      // restore initial bound
258      CoinCopyN (chg_bds, ncols, f_chg);
259      CoinCopyN (olb, ncols, Lb ());
260      CoinCopyN (oub, ncols, Ub ());
261
262      //#if BR_TEST_LOG < 0 // for fair testing
263      // check tightened problem for feasibility
264      if (!(btCore (chg_bds))) {
265
266        jnlst_ -> Printf (Ipopt::J_ERROR, J_BOUNDTIGHTENING, 
267                          "\n    pruned by Probing\n");
268        return -1;
269
270      } else {
271
272        // bounds further tightened should be saved
273       
274        CoinCopyN (Lb (), ncols, olb);
275        CoinCopyN (Ub (), ncols, oub);
276      }
277      //#endif
278    }
279
280    // TODO: compute golden section
281    //fb = (inner + outer) / 2;
282
283    //fb = fictBounds (direction, fb, CoinMin (inner, outer), CoinMax (inner, outer));
284
285    // inner and outer might have to change. Update
286
287    CouNumber
288      lb = Lb (index),
289      ub = Ub (index);
290
291    if ((!direction && ((inner > ub) || (outer < lb))) ||
292        ( direction && ((inner < lb) || (outer > ub)))) {
293
294      // keep it simple
295
296      inner = direction ? lb : ub;
297      outer = direction ? ub : lb;
298    }
299
300    CouNumber diff = fabs (inner - outer);
301
302    if (diff <= COUENNE_EPS) break;
303
304    if (diff > 1.) {
305
306      CouNumber L = log (diff) / log (10.);
307
308      if (direction) fb = inner + exp (log (10.) * L/2);
309      else           fb = inner - exp (log (10.) * L/2);
310
311    } else fb = (inner + outer)/2;
312
313    //    if () fb = (          inner + (phi-1) * outer) / phi;
314    //    else  fb = ((phi-1) * inner +           outer) / phi;
315
316    //  if (!feasible)       
317    //    fb = fictBounds (direction, xcur,
318    //       direction ? lb [index] : outer,
319    //       direction ? outer      : ub [index]);
320
321    jnlst_ -> Printf (Ipopt::J_ERROR, J_BOUNDTIGHTENING, "\n");
322  }
323
324  Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_BOUNDTIGHTENING, "\n");
325  if (tightened) Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_BOUNDTIGHTENING, "  [x%2d] pruned %s [%g, %g] -- lb = %g cutoff = %g\n", index,direction?"right":"left", olb[index],oub[index], objind >= 0 ? Lb (objind) : 0., getCutOff ());
326
327  return tightened ? 1 : 0;
328}
Note: See TracBrowser for help on using the repository browser.