CouenneLPtightenBoundsCLP.cpp
Go to the documentation of this file.
1 /* $Id: CouenneLPtightenBoundsCLP.cpp 490 2011-01-14 16:07:12Z pbelotti $
2  *
3  * Name: CouenneLPtightenBoundsCLP.cpp
4  * Authors: Pietro Belotti, Carnegie Mellon University
5  * Purpose: adaptation of OsiClpSolverInterface::tightenBounds()
6  *
7  * (C) Carnegie-Mellon University, 2009.
8  * This file is licensed under the Eclipse Public License (EPL)
9  */
10 
11 #include "CouennePrecisions.hpp"
12 #include "CouenneProblem.hpp"
13 #include "CouenneCutGenerator.hpp"
14 #include "CouenneExprVar.hpp"
15 
16 namespace Couenne {
17 
18 //#define COIN_DEVELOP 4
19 
20 // Tighten bounds. Returns -1 if infeasible, otherwise number of
21 // variables tightened.
22 template <class T>
23 int CouenneSolverInterface<T>::tightenBoundsCLP (int lightweight) {
24 
25  // Copied from OsiClpSolverInterface::tightenBounds
26 
27  int
28  numberRows = T::getNumRows(),
29  numberColumns = T::getNumCols(),
30  iRow, iColumn;
31 
32  const double * columnUpper = T::getColUpper();
33  const double * columnLower = T::getColLower();
34  const double * rowUpper = T::getRowUpper();
35  const double * rowLower = T::getRowLower();
36 
37  // Column copy of matrix
38  const double * element = T::getMatrixByCol()->getElements();
39  const int * row = T::getMatrixByCol()->getIndices();
40  const CoinBigIndex * columnStart = T::getMatrixByCol()->getVectorStarts();
41  const int * columnLength = T::getMatrixByCol()->getVectorLengths();
42  const double *objective = T::getObjCoefficients() ;
43 
44  double direction = T::getObjSense();
45  double * down = new double [numberRows];
46 
47  if (lightweight)
48  return tightenBoundsCLP_Light (lightweight);
49 
50  // NOT LIGHTWEIGHT /////////////////////////////////////////////////////
51 
52  double * up = new double [numberRows];
53  double * sum = new double [numberRows];
54  int * type = new int [numberRows];
55  CoinZeroN(down,numberRows);
56  CoinZeroN(up,numberRows);
57  CoinZeroN(sum,numberRows);
58  CoinZeroN(type,numberRows);
59  double infinity = T::getInfinity();
60 
61  for (iColumn=0;iColumn<numberColumns;iColumn++) {
62  CoinBigIndex start = columnStart[iColumn];
63  CoinBigIndex end = start + columnLength[iColumn];
64  double lower = columnLower[iColumn];
65  double upper = columnUpper[iColumn];
66  if (lower==upper) {
67  for (CoinBigIndex j=start;j<end;j++) {
68  int iRow = row[j];
69  double value = element[j];
70  sum[iRow]+=2.0*fabs(value*lower);
71  if ((type[iRow]&1)==0)
72  down[iRow] += value*lower;
73  if ((type[iRow]&2)==0)
74  up[iRow] += value*lower;
75  }
76  } else {
77  for (CoinBigIndex j=start;j<end;j++) {
78  int iRow = row[j];
79  double value = element[j];
80  if (value>0.0) {
81  if ((type[iRow]&1)==0) {
82  if (lower!=-infinity) {
83  down[iRow] += value*lower;
84  sum[iRow]+=fabs(value*lower);
85  } else {
86  type[iRow] |= 1;
87  }
88  }
89  if ((type[iRow]&2)==0) {
90  if (upper!=infinity) {
91  up[iRow] += value*upper;
92  sum[iRow]+=fabs(value*upper);
93  } else {
94  type[iRow] |= 2;
95  }
96  }
97  } else {
98  if ((type[iRow]&1)==0) {
99  if (upper!=infinity) {
100  down[iRow] += value*upper;
101  sum[iRow]+=fabs(value*upper);
102  } else {
103  type[iRow] |= 1;
104  }
105  }
106  if ((type[iRow]&2)==0) {
107  if (lower!=-infinity) {
108  up[iRow] += value*lower;
109  sum[iRow]+=fabs(value*lower);
110  } else {
111  type[iRow] |= 2;
112  }
113  }
114  }
115  }
116  }
117  }
118 
119  int nTightened = 0;
120  double tolerance = 1.0e-6;
121 
122  for (iRow=0;iRow<numberRows;iRow++) {
123  if ((type[iRow]&1)!=0)
124  down[iRow]=-infinity;
125  if (down[iRow]>rowUpper[iRow]) {
126  if (down[iRow]>rowUpper[iRow]+tolerance+1.0e-8*sum[iRow]) {
127  // infeasible
128 #ifdef COIN_DEVELOP
129  printf("infeasible on row %d\n",iRow);
130 #endif
131  nTightened=-1;
132  break;
133  } else {
134  down[iRow]=rowUpper[iRow];
135  }
136  }
137  if ((type[iRow]&2)!=0)
138  up[iRow]=infinity;
139  if (up[iRow]<rowLower[iRow]) {
140  if (up[iRow]<rowLower[iRow]-tolerance-1.0e-8*sum[iRow]) {
141  // infeasible
142 #ifdef COIN_DEVELOP
143  printf("infeasible on row %d\n",iRow);
144 #endif
145  nTightened=-1;
146  break;
147  } else {
148  up[iRow]=rowLower[iRow];
149  }
150  }
151  }
152 
153  if (nTightened)
154  numberColumns = 0; // so will skip
155 
156  for (iColumn=0;iColumn<numberColumns;iColumn++) {
157  double lower = columnLower[iColumn];
158  double upper = columnUpper[iColumn];
159  double gap = upper-lower;
160 
161  if (!gap)
162  continue;
163 
164  int canGo=0;
165 
166  CoinBigIndex
167  start = columnStart [iColumn],
168  end = start + columnLength [iColumn];
169 
170  if (lower < -1.0e8 && upper > 1.0e8)
171  continue; // Could do severe damage to accuracy
172 
173 
174  // there was an ifInteger condition here. We do like tightened
175  // bounds for continuous variables too, so we don't test for
176  // integrality.
177 
178  std::vector <exprVar *> &vars = cutgen_ -> Problem () -> Variables ();
179 
180  {
181  if (vars [iColumn] -> isInteger ()) {
182 
183  if (lower < ceil (lower - COUENNE_EPS) - COUENNE_EPS) {
184 #ifdef COIN_DEVELOP
185  printf("increasing lower bound on %d from %e to %e\n",iColumn,
186  lower,ceil(lower - COUENNE_EPS));
187 #endif
188  lower=ceil(lower - COUENNE_EPS);
189  gap=upper-lower;
190  T::setColLower(iColumn,lower);
191  }
192 
193  if (upper > floor(upper + COUENNE_EPS) + COUENNE_EPS) {
194 #ifdef COIN_DEVELOP
195  printf("decreasing upper bound on %d from %e to %e\n",iColumn,
196  upper,floor(upper + COUENNE_EPS));
197 #endif
198  upper=floor(upper + COUENNE_EPS);
199  gap=upper-lower;
200  T::setColUpper(iColumn,upper);
201  }
202  }
203 
204  double newLower=lower;
205  double newUpper=upper;
206 
207  for (CoinBigIndex j=start;j<end;j++) {
208  int iRow = row[j];
209  double value = element[j];
210  if (value>0.0) {
211  if ((type[iRow]&1)==0) {
212  // has to be at most something
213  if (down[iRow] + value*gap > rowUpper[iRow]+tolerance) {
214  double newGap = (rowUpper[iRow]-down[iRow])/value;
215  // adjust
216  newGap += 1.0e-10*sum[iRow];
217  if (vars [iColumn] -> isInteger ())
218  newGap = floor(newGap);
219  if (lower+newGap<newUpper)
220  newUpper=lower+newGap;
221  }
222  }
223  if (down[iRow]<rowLower[iRow])
224  canGo |=1; // can't go down without affecting result
225  if ((type[iRow]&2)==0) {
226  // has to be at least something
227  if (up[iRow] - value*gap < rowLower[iRow]-tolerance) {
228  double newGap = (up[iRow]-rowLower[iRow])/value;
229  // adjust
230  newGap += 1.0e-10*sum[iRow];
231  if (vars [iColumn] -> isInteger ())
232  newGap = floor(newGap);
233  if (upper-newGap>newLower)
234  newLower=upper-newGap;
235  }
236  }
237  if (up[iRow]>rowUpper[iRow])
238  canGo |=2; // can't go up without affecting result
239  } else {
240  if ((type[iRow]&1)==0) {
241  // has to be at least something
242  if (down[iRow] - value*gap > rowUpper[iRow]+tolerance) {
243  double newGap = -(rowUpper[iRow]-down[iRow])/value;
244  // adjust
245  newGap += 1.0e-10*sum[iRow];
246  if (vars [iColumn] -> isInteger ())
247  newGap = floor(newGap);
248  if (upper-newGap>newLower)
249  newLower=upper-newGap;
250  }
251  }
252  if (up[iRow]>rowUpper[iRow])
253  canGo |=1; // can't go down without affecting result
254  if ((type[iRow]&2)==0) {
255  // has to be at most something
256  if (up[iRow] + value*gap < rowLower[iRow]-tolerance) {
257  double newGap = -(up[iRow]-rowLower[iRow])/value;
258  // adjust
259  newGap += 1.0e-10*sum[iRow];
260  if (vars [iColumn] -> isInteger ())
261  newGap = floor(newGap);
262  if (lower+newGap<newUpper)
263  newUpper=lower+newGap;
264  }
265  }
266  if (down[iRow]<rowLower[iRow])
267  canGo |=2; // can't go up without affecting result
268  }
269  }
270 
271  if (newUpper<upper || newLower>lower) {
272  nTightened++;
273  if (newLower>newUpper) {
274  // infeasible
275 #if COIN_DEVELOP>1
276  printf("infeasible on column %d\n",iColumn);
277 #endif
278  nTightened=-1;
279  break;
280  } else {
281  T::setColLower(iColumn,newLower);
282  T::setColUpper(iColumn,newUpper);
283  }
284  for (CoinBigIndex j=start;j<end;j++) {
285  int iRow = row[j];
286  double value = element[j];
287  if (value>0.0) {
288  if ((type[iRow]&1)==0) down [iRow] += value*(newLower-lower);
289  if ((type[iRow]&2)==0) up [iRow] += value*(newUpper-upper);
290  } else {
291  if ((type[iRow]&1)==0) down [iRow] += value*(newUpper-upper);
292  if ((type[iRow]&2)==0) up [iRow] += value*(newLower-lower);
293  }
294  }
295  } else {
296 
297  if (canGo!=3) {
298 
299  double objValue = direction*objective[iColumn];
300 
301  if (objValue>=0.0&&(canGo&1)==0) {
302 #if COIN_DEVELOP>2
303  printf("dual fix down on column %d\n",iColumn);
304 #endif
305  nTightened++;
306  T::setColUpper(iColumn,lower);
307  } else if (objValue<=0.0 && (canGo&2)==0) {
308 #if COIN_DEVELOP>2
309  printf("dual fix up on column %d\n",iColumn);
310 #endif
311  nTightened++;
312  T::setColLower(iColumn,upper);
313  }
314  }
315  }
316  }
317 
318 // else {
319 
320 // // CONTINUOUS //////////////////////////////////////////
321 
322 // // just do dual tests
323 // for (CoinBigIndex j=start;j<end;j++) {
324 // int iRow = row[j];
325 // double value = element[j];
326 // if (value>0.0) {
327 // if (down [iRow] < rowLower [iRow]) canGo |=1; // can't go down without affecting result
328 // if (up [iRow] > rowUpper [iRow]) canGo |=2; // can't go up without affecting result
329 // } else {
330 // if (up [iRow] > rowUpper [iRow]) canGo |=1; // can't go down without affecting result
331 // if (down [iRow] < rowLower [iRow]) canGo |=2; // can't go up without affecting result
332 // }
333 // }
334 
335 // if (canGo!=3) {
336 // double objValue = direction*objective[iColumn];
337 // if (objValue>=0.0&&(canGo&1)==0) {
338 // #if COIN_DEVELOP>2
339 // printf("dual fix down on continuous column %d lower %g\n",
340 // iColumn,lower);
341 // #endif
342 // // Only if won't cause numerical problems
343 // if (lower>-1.0e10) {
344 // nTightened++;;
345 // setColUpper(iColumn,lower);
346 // }
347 // } else if (objValue<=0.0&&(canGo&2)==0) {
348 // #if COIN_DEVELOP>2
349 // printf("dual fix up on continuous column %d upper %g\n",
350 // iColumn,upper);
351 // #endif
352 // // Only if won't cause numerical problems
353 // if (upper<1.0e10) {
354 // nTightened++;;
355 // setColLower(iColumn,upper);
356 // }
357 // }
358 // }
359 // }
360 
361  }
362 
363  delete [] type;
364  delete [] down;
365  delete [] up;
366  delete [] sum;
367 
368  return nTightened;
369 }
370 
371 }
Couenne::CouenneSolverInterface::tightenBoundsCLP
virtual int tightenBoundsCLP(int lightweight)
Copy of the Clp version — not light version.
Definition: CouenneLPtightenBoundsCLP.cpp:23
j
static char * j
Definition: OSdtoa.cpp:3622
up
int up
Definition: OSdtoa.cpp:1817
e
void fint fint fint real fint real real real real real real real real real * e
Definition: BonBqpdSolver.cpp:27
end
fint end
Definition: BonFilterSolver.cpp:135
CouenneCutGenerator.hpp
COUENNE_EPS
#define COUENNE_EPS
Definition: CouennePrecisions.hpp:19
CouenneExprVar.hpp
Variables
The in-memory representation of the variables element.
Definition: OSInstance.h:83
CouenneProblem.hpp
CouennePrecisions.hpp
Couenne::isInteger
bool isInteger(CouNumber x)
is this number integer?
Definition: CouenneExpression.hpp:309
Generated by   doxygen 1.8.5