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

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

restored objects_ as vector of POINTERS to object rather than objects -- was not applying complementarity branching

  • Property svn:keywords set to Id
File size: 8.8 KB
Line 
1/* $Id: infeasibilityVT.cpp 259 2009-10-04 14:20:13Z 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    } 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        expression *ref = obj -> Reference ();
132        jnlst_ -> Printf (J_MATRIX, J_BRANCHING, "[%g,%g] --> %g - %g = %g (diff = %g - %g = %g): ", 
133                          left, right, rFeas, lFeas, rFeas - lFeas,
134                          (ref -> Image ()) ? 
135                          (*(ref-> Image ())) () : 0., 
136                          (*(ref)) (),
137                          (ref -> Image ()) ? 
138                          (*(ref -> Image ())) () - (*(ref)) () : 0.);
139        ref -> print (); 
140        if (ref -> Image ()) 
141          {printf (" := "); ref -> Image() -> print();}
142        printf ("\n");
143      }
144    }
145
146    if (lFeas < info -> lower_ [indexVar]) lFeas = info -> lower_ [indexVar];
147    if (rFeas > info -> upper_ [indexVar]) rFeas = info -> upper_ [indexVar];
148
149    retval = rFeas - lFeas;
150
151    upEstimate_   = rFeas - xcurr;
152    downEstimate_ = xcurr - lFeas;
153
154    // need a nonzero value for the estimate. Although maxInf is
155    // related to the expressions that depend on this variable, it is
156    // still an estimate of how the point will change
157    if (upEstimate_   <= tol) upEstimate_   = maxInf;
158    if (downEstimate_ <= tol) downEstimate_ = maxInf;
159  }
160
161  //////////////////////////////////////////////////////////////////////
162
163  // object is infeasible.
164  // check how in the middle of the interval we are
165
166  const CouNumber threshold = .5; // can be in [0,0.5]
167
168  if (retval > COUENNE_EPS) {
169
170    leanLeft  = (xcurr - lFeas) / retval;
171
172    if      (leanLeft <     threshold) way = 0;
173    else if (leanLeft > 1 - threshold) way = 1;
174    else way = (CoinDrand48 () < 0.5) ? 0 : 1;
175  }
176
177  // done computing delta. Now transfer violation on LP relaxation ///////////////
178
179  // coefficient in objective is in {0,1} and there is only one
180  // variable with coefficient 1
181  CouNumber vt_delta =
182    (indexVar == problem_ -> Obj (0) -> Body () -> Index ()) ? 1. : 0.;
183
184  for (int i=0, n_el = info -> columnLength_ [indexVar]; i < n_el; i++) {
185
186    int indRow = info -> columnStart_ [indexVar] + i;
187
188    vt_delta += 
189      fabs (info -> pi_ [info -> row_ [indRow]] * 
190            info -> elementByColumn_  [indRow]);
191
192    //     vt_delta +=
193    //       info -> pi_ [info -> row_ [indRow]] *
194    //       fabs (info -> elementByColumn_  [indRow]);
195
196    jnlst_ -> Printf (J_MATRIX, J_BRANCHING, "+ (pi[%d]=%g) * (el[%d]=%g) [=%g] --> vtd = %g\n",
197                      info -> row_ [indRow],
198                      info -> pi_ [info -> row_ [indRow]], 
199                      indRow,
200                      info -> elementByColumn_  [indRow],
201                      info -> pi_ [info -> row_ [indRow]] * 
202                      info -> elementByColumn_  [indRow],
203                      vt_delta);
204  }
205
206  // weights for VT itself and width of interval
207  const CouNumber
208    alpha = 1.0,
209    beta  = 0.0;
210
211  jnlst_ -> Printf (J_MATRIX, J_BRANCHING, "return %g * %g + %g * %g + %g * %g --> ", 
212                    alpha, fabs (retval * vt_delta), beta, retval,
213                    1-alpha-beta, leanLeft * (1-leanLeft));
214
215  retval = 
216    alpha          * fabs (retval*vt_delta) +  // violation transfer itself
217    beta           * retval +                  // width of feasibility interval
218    (1-alpha-beta) * leanLeft * (1-leanLeft);  // how in the middle of the interval x is
219
220  if (jnlst_ -> ProduceOutput (J_MATRIX, J_BRANCHING)) {
221    if (retval > CoinMin (COUENNE_EPS, feas_tolerance_)) {
222      printf ("vt-delta is %-10g [", retval); 
223      reference_ -> print (); 
224      if (reference_ -> Image ()) { // if no list, print image
225        printf (" := ");
226        reference_ -> Image () -> print ();
227      }
228      if (dependence.size () > 0) {
229        printf (" -- ");
230        for (std::set <int>::const_iterator i = dependence.begin (); 
231             i != dependence.end (); ++i) {
232          problem_ -> Var (*i) -> print ();
233          printf (" ");
234        }
235      }
236      printf ("]\n");
237    } else printf ("feasible...\n");
238  }
239
240  problem_ -> domain () -> pop ();
241
242  if ((retval < tol) &&
243      (maxInf > tol)) {
244
245    // no improvement with this variable, but need to return nonzero
246    // if this is an auxiliary whose value is different from relative
247    // expression
248    retval = maxInf; 
249  }
250
251  if (retval < CoinMin (COUENNE_EPS, feas_tolerance_)) 
252    retval = 0.;
253
254#define ALMOST_ZERO 1e-8
255
256  assert ((retval < ALMOST_ZERO && upEstimate_ < ALMOST_ZERO && downEstimate_ < ALMOST_ZERO) ||
257          (retval > ALMOST_ZERO && upEstimate_ > ALMOST_ZERO && downEstimate_ > ALMOST_ZERO));
258
259  return (reference_ -> isInteger ()) ? 
260    CoinMax (retval, intInfeasibility (info -> solution_ [reference_ -> Index ()])) :
261    retval;
262}
Note: See TracBrowser for help on using the repository browser.