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