source: trunk/Couenne/src/branch/infeasibilityVT.cpp @ 43

Last change on this file since 43 was 43, checked in by pbelotti, 11 years ago

restored previous fix

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