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