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

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

fixing constant-objective glitches (merge from trunk)

  • Property svn:keywords set to Author Date Id Revision
File size: 10.2 KB
Line 
1/* $Id: fake_tightening.cpp 802 2012-01-30 23:11:34Z 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 ? -sqrt (-lb) : sqrt (ub));
40
41      //if (fabs (x) < COUENNE_EPS) return (direction ? AGGR_MUL : - AGGR_MUL);
42      //else                        return (direction ? AGGR_MUL : - AGGR_MUL) * fabs (x);
43
44    } else { // ]-inf,u]
45
46      if (!direction)
47        return -sqrt (-lb); // try to tighten interval from a very large value
48
49      if      (x < -COUENNE_EPS) return (CoinMin (0., (x+ub)/2));
50      else if (x >  COUENNE_EPS) return ((x + (ub-x)/AGGR_DIV));
51      else                       return ((ub/AGGR_DIV));
52
53      // if      (x < -COUENNE_EPS) return (direction ? CoinMin (0., (x+ub)/2) : AGGR_MUL * x);
54      // else if (x >  COUENNE_EPS) return (direction ? (x + (ub-x)/AGGR_DIV) : 0);
55      // else                       return (direction ? (ub/AGGR_DIV) : -AGGR_MUL);
56    }
57  }
58  else {
59    if (ub >  LARGE_BOUND) { // [l,+inf[
60
61      if (direction)
62        return sqrt (ub);
63
64      if      (x < -COUENNE_EPS) return ((x - (x-lb) / AGGR_DIV));
65      else if (x >  COUENNE_EPS) return (CoinMax (0.,(x+lb)/2));
66      else                       return (lb/AGGR_DIV);
67
68      // if      (x < -COUENNE_EPS) return (direction ? 0 : (x - (x-lb) / AGGR_DIV));
69      // else if (x >  COUENNE_EPS) return (direction ? (AGGR_MUL * x) : CoinMax (0.,(x+lb)/2));
70      // else                       return (direction ? AGGR_MUL : lb/AGGR_DIV);
71
72    } else // [l,u]
73      return (direction ? 
74              (x + (ub-x) / AGGR_DIV) : 
75              (x - (x-lb) / AGGR_DIV));
76  }
77}
78
79
80// Single fake tightening. Return
81//
82// -1   if infeasible
83//  0   if no improvement
84// +1   if improved
85int CouenneProblem::
86fake_tighten (char direction,  // 0: left, 1: right
87              int index,       // index of the variable tested
88              const double *X, // point round which tightening is done
89              CouNumber *olb,  // cur. lower bound
90              CouNumber *oub,  //      upper
91              t_chg_bounds *chg_bds,
92              t_chg_bounds *f_chg) const {
93
94  int
95    ncols    = nVars (),
96    objind   = Obj (0) -> Body  () -> Index ();
97
98  //assert (objind >= 0);
99
100  bool 
101    tightened = false,
102    intvar    = variables_ [index] -> isInteger ();
103
104  CouNumber
105    xcur      = X [index],
106    inner     = xcur,                                                 // point closest to current x
107    outer     = (direction ? oub : olb) [index],                      // point closest to bound
108    fb        = fictBounds (direction, xcur, Lb (index), Ub (index)); // starting point
109
110  // This is a one-dimensional optimization problem between inner and
111  // outer, on a monotone function of which we can compute the value
112  // (with relative expense) but not the derivative.
113
114  jnlst_ -> Printf (Ipopt::J_ERROR, J_BOUNDTIGHTENING, 
115                    "  x_%d.  x = %10g, lb = %g, cutoff = %g-----------------\n", 
116                    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,
180                      " [%10g,%10g] lb = %g {fea=%d,btr=%d} ",
181                      Lb (index), Ub (index), objind >= 0 ? Lb (objind) : 0., feasible,betterbds);
182
183    if (feasible && !betterbds) {
184
185      // case 1: too tight, move inner out
186      inner = fb;
187
188      // restore initial bound
189      CoinCopyN (chg_bds, ncols, f_chg);
190      CoinCopyN (olb, ncols, Lb ());
191      CoinCopyN (oub, ncols, Ub ());
192
193    } else {
194
195      // Here, !feasible || betterbds
196      //
197      // If !feasible OR
198      //    (betterbds and the new lb is above the cutoff)
199      //
200      // then there is a tightening
201
202      // case 2: tightening valid, apply and move outer in
203
204      //printf (" --> %cbound [x_%d]: %g --> %g",direction?'U':'L',index,(direction?oub:olb)[index],fb);
205
206      if (optimum_ && 
207          ((!direction &&
208            (optimum_ [index] >= olb [index]) && 
209            (optimum_ [index] <= Lb (index) - COUENNE_EPS)) ||
210           (direction &&
211            (optimum_ [index] <= oub [index]) && 
212            (optimum_ [index] >= Ub (index) + COUENNE_EPS)))) {
213
214        jnlst_ -> Printf (Ipopt::J_ERROR, J_BOUNDTIGHTENING, 
215                          "fake tightening CUTS optimum: x%d=%g in [%g,%g] but not in [%g,%g]\n",
216                          index, olb [index], oub [index], Lb (index), Ub (index));
217      }
218
219      /*bool do_not_tighten = false;
220
221      // check if cut best known solution
222      if (optimum_) {
223        if (direction) {
224          if ((oub [index] > optimum_ [index]) &&
225              (fb          < optimum_ [index])) {
226            printf ("aggressive bt cuts optimum ub %d: %g < %g < %g\n",
227                    index, fb, optimum_ [index], oub [index]);
228            do_not_tighten = true;
229          }
230        } else {
231          if ((olb [index] < optimum_ [index]) &&
232              (fb          > optimum_ [index])) {
233            printf ("aggressive bt cuts optimum lb %d: %g < %g < %g\n",
234                    index, olb [index], optimum_ [index], fb);
235            do_not_tighten = true;
236          }
237        }
238        }*/
239
240      //if (!do_not_tighten) {
241
242      // apply bound
243      if (direction) {
244
245        oub [index] = Ub (index) = intvar ? floor (fb + COUENNE_EPS) : fb; 
246        chg_bds [index]. setUpper (t_chg_bounds::CHANGED);
247
248      } else {
249
250        olb [index] = Lb (index) = intvar ? ceil  (fb - COUENNE_EPS) : fb; 
251        chg_bds [index]. setLower (t_chg_bounds::CHANGED);
252      }
253
254      outer = fb; // we have at least a tightened bound, save it
255
256      tightened = true;
257      //}
258
259      // restore initial bound
260      CoinCopyN (chg_bds, ncols, f_chg);
261      CoinCopyN (olb, ncols, Lb ());
262      CoinCopyN (oub, ncols, Ub ());
263
264      //#if BR_TEST_LOG < 0 // for fair testing
265      // check tightened problem for feasibility
266      if (!(btCore (chg_bds))) {
267
268        jnlst_ -> Printf (Ipopt::J_ERROR, J_BOUNDTIGHTENING, 
269                          "\n    pruned by Probing\n");
270        return -1;
271
272      } else {
273
274        // bounds further tightened should be saved
275       
276        CoinCopyN (Lb (), ncols, olb);
277        CoinCopyN (Ub (), ncols, oub);
278      }
279      //#endif
280    }
281
282    // TODO: compute golden section
283    //fb = (inner + outer) / 2;
284
285    //fb = fictBounds (direction, fb, CoinMin (inner, outer), CoinMax (inner, outer));
286
287    // inner and outer might have to change. Update
288
289    CouNumber
290      lb = Lb (index),
291      ub = Ub (index);
292
293    if ((!direction && ((inner > ub) || (outer < lb))) ||
294        ( direction && ((inner < lb) || (outer > ub)))) {
295
296      // keep it simple
297
298      inner = direction ? lb : ub;
299      outer = direction ? ub : lb;
300    }
301
302    CouNumber diff = fabs (inner - outer);
303
304    if (diff <= COUENNE_EPS) break;
305
306    if (diff > 1.) {
307
308      CouNumber L = log (diff) / log (10.);
309
310      if (direction) fb = inner + exp (log (10.) * L/2);
311      else           fb = inner - exp (log (10.) * L/2);
312
313    } else fb = (inner + outer)/2;
314
315    //    if () fb = (          inner + (phi-1) * outer) / phi;
316    //    else  fb = ((phi-1) * inner +           outer) / phi;
317
318    //  if (!feasible)       
319    //    fb = fictBounds (direction, xcur,
320    //       direction ? lb [index] : outer,
321    //       direction ? outer      : ub [index]);
322
323    jnlst_ -> Printf (Ipopt::J_ERROR, J_BOUNDTIGHTENING, "\n");
324  }
325
326  Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_BOUNDTIGHTENING, "\n");
327  if (tightened) 
328    Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_BOUNDTIGHTENING, 
329                    "  [x%2d] pruned %s [%g, %g] -- lb = %g cutoff = %g\n", 
330                    index,direction?"right":"left",
331                    olb[index],oub[index], objind >= 0 ? Lb (objind) : 0., getCutOff ());
332
333  return tightened ? 1 : 0;
334}
Note: See TracBrowser for help on using the repository browser.