source: stable/0.3/Couenne/src/branch/infeasibilityVT.cpp @ 581

Last change on this file since 581 was 581, checked in by pbelotti, 9 years ago

intInfeasibility() checks lower/upper bounds to avoid non-null infeasibility upon out-of-bounds values, which triggers infeasibility in Couenne but not in Cbc, and this in turn makes the BB loop

  • Property svn:keywords set to Id
File size: 8.9 KB
Line 
1/* $Id: infeasibilityVT.cpp 581 2011-05-22 11:37:33Z pbelotti $
2 *
3 * Name:    infeasibilityVT.cpp
4 * Authors: Pietro Belotti, Carnegie Mellon University
5 * Purpose: Compute violation transfer of a variable. See Tawarmalani
6 *          and Sahinidis' book or article on MathProg 2002
7 *
8 * (C) Carnegie-Mellon University, 2008-09.
9 * This file is licensed under the Common Public License (CPL)
10 */
11
12#include "CoinHelperFunctions.hpp"
13#include "CouenneProblem.hpp"
14#include "CouenneVTObject.hpp"
15
16
17/// compute infeasibility of this variable, |w - f(x)| (where w is the
18/// auxiliary variable defined as w = f(x))
19double CouenneVTObject::infeasibility (const OsiBranchingInformation *info, int &way) const {
20
21  int indexVar = reference_ -> Index ();
22  assert (indexVar >= 0);
23
24  CouNumber tol = CoinMin (COUENNE_EPS, feas_tolerance_);
25
26  // feasible variable
27  if (info -> upper_ [indexVar] - 
28      info -> lower_ [indexVar] < tol) {
29    if (reference_ -> isInteger ()) {
30      double point = info -> solution_ [reference_ -> Index ()];
31      if (downEstimate_ <       point  - floor (point)) downEstimate_ =       point  - floor (point);
32      if (upEstimate_   < ceil (point) -        point)  upEstimate_   = ceil (point) -        point;
33      return intInfeasibility (point,
34                               info -> lower_ [indexVar],
35                               info -> upper_ [indexVar]);
36    } else return (upEstimate_ = downEstimate_ = 0.);
37  }
38
39  // let Couenne know of current status (variables, bounds)
40  problem_ -> domain () -> push
41    (problem_ -> nVars (),
42     info -> solution_, 
43     info -> lower_, 
44     info -> upper_);
45
46  // debug output ////////////////////////////////////////////////////////////////
47
48  if (jnlst_ -> ProduceOutput (J_DETAILED, J_BRANCHING)) {
49    printf ("VT infeas on ");
50    reference_ -> print ();
51    if (reference_ -> Image ()) { // if no list, print image
52      printf (" := ");
53      reference_ -> Image () -> print ();
54    }
55    const std::set <int> &dependence = problem_ -> Dependence () [indexVar];
56    if (dependence.size () > 0) {
57      printf (" -- ");
58      for (std::set <int>::const_iterator i = dependence.begin (); 
59           i != dependence.end (); ++i) {
60        problem_ -> Var (*i) -> print ();
61        printf (" ");
62      }
63    }
64    printf ("\n");
65  }
66
67  // get set of variable indices that depend on reference_
68  const std::set <int> &dependence = problem_ -> Dependence () [indexVar];
69
70  CouNumber
71    xcurr    = info -> solution_ [indexVar], // current value of variable
72    retval   = 0.,                           // will be returned
73    lFeas    = xcurr,                        // left-most feasible point with all but x_i constant
74    rFeas    = xcurr,                        // right-most ...
75    leanLeft = 0.,                           // [0,1], 
76    fx       = xcurr,                        // value of expression associated with variable (if aux)
77    maxInf   = 0.;                           // max infeasibility over expressions depending on this
78
79  if (reference_ -> Type () == AUX)
80    fx = (*(reference_ -> Image ())) ();
81
82  if (dependence.size () == 0) { // this is a top level auxiliary,
83                                 // nowhere an independent
84
85    // that means, for VT, that letting w vary and keeping everything
86    // else constant we return the difference between current value
87    // and value of function at this point
88
89    // these variables may also be disabled in BonCouenneSetup.cpp
90
91    retval = (reference_ -> Type () == AUX) ? // if this is an isolated variable
92      // check if this w=f(x) is used nowhere and is feasible
93      (upEstimate_ = downEstimate_ = maxInf = checkInfeasibility (info)) :
94      0.;
95
96    // check if this w=f(x) is used nowhere and is feasible
97    retval = upEstimate_ = downEstimate_ = maxInf = checkInfeasibility (info);
98
99  } else {
100
101    // this appears as independent in all auxs of the "dependence" list
102    // check all non-linear objects containing this variable
103
104    for (std::set <int>::const_iterator i = dependence.begin ();
105         i != dependence.end (); ++i) {
106
107      const CouenneObject *obj = problem_ -> Objects () [*i];
108      assert (obj -> Reference ());
109
110      CouNumber
111        left   = xcurr,
112        right  = xcurr,
113        infeas = obj -> checkInfeasibility (info);
114
115      if (infeas > maxInf)
116        maxInf = infeas;
117
118      // check feasibility of nonlinear object
119      if (infeas > 0.) {
120
121        // measure how far have to go, left and right, to restore
122        // feasibility (see page 582, Tawarmalani and Sahinidis,
123        // MathProg A: 99, pp. 563-591)
124        //if (obj. Reference () -> Image ()) // not needed! obj only has nonlinear objects
125        obj -> Reference () -> Image () -> closestFeasible
126          (reference_, obj -> Reference (), left, right);
127
128        if (left  < lFeas) lFeas = left;
129        if (right > rFeas) rFeas = right;
130      }
131
132      if (jnlst_ -> ProduceOutput (J_MATRIX, J_BRANCHING)) { // debug output
133        expression *ref = obj -> Reference ();
134        jnlst_ -> Printf (J_MATRIX, J_BRANCHING, "[%g,%g] --> %g - %g = %g (diff = %g - %g = %g): ", 
135                          left, right, rFeas, lFeas, rFeas - lFeas,
136                          (ref -> Image ()) ? 
137                          (*(ref-> Image ())) () : 0., 
138                          (*(ref)) (),
139                          (ref -> Image ()) ? 
140                          (*(ref -> Image ())) () - (*(ref)) () : 0.);
141        ref -> print (); 
142        if (ref -> Image ()) 
143          {printf (" := "); ref -> Image() -> print();}
144        printf ("\n");
145      }
146    }
147
148    if (lFeas < info -> lower_ [indexVar]) lFeas = info -> lower_ [indexVar];
149    if (rFeas > info -> upper_ [indexVar]) rFeas = info -> upper_ [indexVar];
150
151    retval = rFeas - lFeas;
152
153    upEstimate_   = rFeas - xcurr;
154    downEstimate_ = xcurr - lFeas;
155
156    // need a nonzero value for the estimate. Although maxInf is
157    // related to the expressions that depend on this variable, it is
158    // still an estimate of how the point will change
159    if (upEstimate_   <= tol) upEstimate_   = maxInf;
160    if (downEstimate_ <= tol) downEstimate_ = maxInf;
161  }
162
163  //////////////////////////////////////////////////////////////////////
164
165  // object is infeasible.
166  // check how in the middle of the interval we are
167
168  const CouNumber threshold = .5; // can be in [0,0.5]
169
170  if (retval > COUENNE_EPS) {
171
172    leanLeft  = (xcurr - lFeas) / retval;
173
174    if      (leanLeft <     threshold) way = 0;
175    else if (leanLeft > 1 - threshold) way = 1;
176    else way = (CoinDrand48 () < 0.5) ? 0 : 1;
177  }
178
179  // done computing delta. Now transfer violation on LP relaxation ///////////////
180
181  // coefficient in objective is in {0,1} and there is only one
182  // variable with coefficient 1
183  CouNumber vt_delta =
184    (indexVar == problem_ -> Obj (0) -> Body () -> Index ()) ? 1. : 0.;
185
186  for (int i=0, n_el = info -> columnLength_ [indexVar]; i < n_el; i++) {
187
188    int indRow = info -> columnStart_ [indexVar] + i;
189
190    vt_delta += 
191      fabs (info -> pi_ [info -> row_ [indRow]] * 
192            info -> elementByColumn_  [indRow]);
193
194    //     vt_delta +=
195    //       info -> pi_ [info -> row_ [indRow]] *
196    //       fabs (info -> elementByColumn_  [indRow]);
197
198    jnlst_ -> Printf (J_MATRIX, J_BRANCHING, "+ (pi[%d]=%g) * (el[%d]=%g) [=%g] --> vtd = %g\n",
199                      info -> row_ [indRow],
200                      info -> pi_ [info -> row_ [indRow]], 
201                      indRow,
202                      info -> elementByColumn_  [indRow],
203                      info -> pi_ [info -> row_ [indRow]] * 
204                      info -> elementByColumn_  [indRow],
205                      vt_delta);
206  }
207
208  // weights for VT itself and width of interval
209  const CouNumber
210    alpha = 1.0,
211    beta  = 0.0;
212
213  jnlst_ -> Printf (J_MATRIX, J_BRANCHING, "return %g * %g + %g * %g + %g * %g --> ", 
214                    alpha, fabs (retval * vt_delta), beta, retval,
215                    1-alpha-beta, leanLeft * (1-leanLeft));
216
217  retval = 
218    alpha          * fabs (retval*vt_delta) +  // violation transfer itself
219    beta           * retval +                  // width of feasibility interval
220    (1-alpha-beta) * leanLeft * (1-leanLeft);  // how in the middle of the interval x is
221
222  if (jnlst_ -> ProduceOutput (J_MATRIX, J_BRANCHING)) {
223    if (retval > CoinMin (COUENNE_EPS, feas_tolerance_)) {
224      printf ("vt-delta is %-10g [", retval); 
225      reference_ -> print (); 
226      if (reference_ -> Image ()) { // if no list, print image
227        printf (" := ");
228        reference_ -> Image () -> print ();
229      }
230      if (dependence.size () > 0) {
231        printf (" -- ");
232        for (std::set <int>::const_iterator i = dependence.begin (); 
233             i != dependence.end (); ++i) {
234          problem_ -> Var (*i) -> print ();
235          printf (" ");
236        }
237      }
238      printf ("]\n");
239    } else printf ("feasible...\n");
240  }
241
242  problem_ -> domain () -> pop ();
243
244  if ((retval < tol) &&
245      (maxInf > tol)) {
246
247    // no improvement with this variable, but need to return nonzero
248    // if this is an auxiliary whose value is different from relative
249    // expression
250    retval = maxInf; 
251  }
252
253  if (retval < CoinMin (COUENNE_EPS, feas_tolerance_)) 
254    retval = 0.;
255
256#define ALMOST_ZERO 1e-8
257
258  assert ((retval < ALMOST_ZERO && upEstimate_ < ALMOST_ZERO && downEstimate_ < ALMOST_ZERO) ||
259          (retval > ALMOST_ZERO && upEstimate_ > ALMOST_ZERO && downEstimate_ > ALMOST_ZERO));
260
261  return (reference_ -> isInteger ()) ? 
262    CoinMax (retval, intInfeasibility (info -> solution_ [indexVar],
263                                       info -> lower_    [indexVar],
264                                       info -> upper_    [indexVar])) :
265    retval;
266}
Note: See TracBrowser for help on using the repository browser.