1 | |
---|
2 | /* |
---|
3 | Copyright (C) 2010, International Business Machines Corporation and others. |
---|
4 | All Rights Reserved. |
---|
5 | |
---|
6 | This code is licensed under the terms of the Eclipse Public License (EPL). |
---|
7 | |
---|
8 | $Id$ |
---|
9 | */ |
---|
10 | |
---|
11 | #if defined(_MSC_VER) |
---|
12 | // Turn off compiler warning about long names |
---|
13 | # pragma warning(disable:4786) |
---|
14 | #endif |
---|
15 | |
---|
16 | #define PROBING100 0 |
---|
17 | |
---|
18 | // #define PRINT_DEBUG |
---|
19 | |
---|
20 | #include "CoinPackedMatrix.hpp" |
---|
21 | #include "CglProbing.hpp" |
---|
22 | #include "CglProbingRowCut.hpp" |
---|
23 | |
---|
24 | // To enable debugging, set CGL_DEBUG in CglProbingDebug.hpp |
---|
25 | #include "CglProbingDebug.hpp" |
---|
26 | |
---|
27 | namespace { |
---|
28 | |
---|
29 | /* |
---|
30 | Allocating one big block for data arrays is efficient but makes it |
---|
31 | difficult for debugging tools to detect accesses outside the boundary of an |
---|
32 | array. The assert checks that a double is larger than an int and that the |
---|
33 | ratio is a power of two. |
---|
34 | */ |
---|
35 | #if CGL_DEBUG == 0 |
---|
36 | |
---|
37 | const unsigned int DIratio = sizeof(double)/sizeof(int) ; |
---|
38 | |
---|
39 | #define ONE_ARRAY |
---|
40 | |
---|
41 | #endif |
---|
42 | |
---|
43 | /* |
---|
44 | Bit constants used to specify bound information. |
---|
45 | |
---|
46 | 0x0: no bounds tightened |
---|
47 | 0x1: u<j> tightened |
---|
48 | 0x2: l<j> tightened |
---|
49 | 0x4: l<j> < -1e10 (-infty) |
---|
50 | 0x8: u<j> > 1e10 (+infty) |
---|
51 | */ |
---|
52 | const unsigned int tightenUpper = 0x1 ; |
---|
53 | const unsigned int tightenLower = 0x2 ; |
---|
54 | const unsigned int infLower = 0x4 ; |
---|
55 | const unsigned int infUpper = 0x8 ; |
---|
56 | |
---|
57 | /* |
---|
58 | Integer and bit constants used to specify probing. |
---|
59 | |
---|
60 | The second down pass is used to restore the down probe state when down |
---|
61 | is feasible but up is not. |
---|
62 | |
---|
63 | The first set (downProbe, upProbe, downProbe2, oneProbeTooMany) *must* |
---|
64 | match the iteration number for the corresponding pass in the PROBE |
---|
65 | loop. The symbolic names make the code a bit easier to read. |
---|
66 | */ |
---|
67 | const unsigned int downProbe = 0 ; |
---|
68 | const unsigned int upProbe = 1 ; |
---|
69 | const unsigned int downProbe2 = 2 ; |
---|
70 | const unsigned int oneProbeTooMany = 3 ; |
---|
71 | |
---|
72 | const unsigned int probeDown = tightenUpper ; |
---|
73 | const unsigned int probeUp = tightenLower ; |
---|
74 | |
---|
75 | const unsigned int downProbeFeas = 0x1 ; |
---|
76 | const unsigned int upProbeFeas = 0x2 ; |
---|
77 | const unsigned int downProbe2Feas = 0x4 ; |
---|
78 | |
---|
79 | } // end file-local namespace |
---|
80 | |
---|
81 | // =============================================== |
---|
82 | |
---|
83 | |
---|
84 | /* |
---|
85 | Run through stackC and create disaggregation cuts. See the typeset |
---|
86 | documentation for the full derivation. Suppose the probe variable is x<p> |
---|
87 | and we're trying to generate a disaggregation cut for target x<t>. |
---|
88 | |
---|
89 | For a down probe, we have orig_u<p> reduced to u<p> = down = colUpper[p]. |
---|
90 | It may be that this causes orig_u<t> to be reduced to u<t>. The upper |
---|
91 | disaggregation cut is |
---|
92 | -(orig_u<t>-u<t>)x<p> + x<t> <= u<t> - u<p>(orig_u<t>-u<t>) |
---|
93 | |
---|
94 | It may be that this causes orig_l<t> to be increased to l<t>. The lower |
---|
95 | disaggregation cut is |
---|
96 | -(orig_l<t>-l<t>)x<p> + x<t> >= l<t> - u<p>(orig_l<t>-l<t>) |
---|
97 | |
---|
98 | For an up probe, we have orig_l<p> increased to l<p> = up = colLower[p]. |
---|
99 | It may be that this causes orig_u<t> to be reduced to u<t>. The upper |
---|
100 | disaggregation cut is |
---|
101 | (orig_u<t>-u<t>)x<p> + x<t> <= u<t> + l<p>(orig_u<t>-u<t>) |
---|
102 | |
---|
103 | It may be that this causes orig_l<t> to be increased to l<t>. The lower |
---|
104 | disaggregation cut is |
---|
105 | (orig_l<t>-l<t>)x<p> + x<t> >= l<t> + l<p>(orig_l<t>-l<t>) |
---|
106 | |
---|
107 | Note that stackC[0] = pndx, the index of the probed variable, so we |
---|
108 | don't want to process stackC[0] in the main loop. |
---|
109 | |
---|
110 | It must be true that goingToTrueBound == 2 and justReplace == false when |
---|
111 | this method is called. |
---|
112 | |
---|
113 | TODO: If you look at the math, the critical thing to get the cut form |
---|
114 | used here is that the old and new bounds on x<p> be separated by |
---|
115 | one. E.g., for a down probe, (orig_u<p>-u<p>) = 1. It looks to |
---|
116 | me like goingToTrueBound might have been an attempt to guarantee this |
---|
117 | condition, but the test is wrong, hence it's patched up with tweak |
---|
118 | that essentially cuts this cut back to binary probing variables. |
---|
119 | -- lh, 101214 -- |
---|
120 | |
---|
121 | probeDir should be coded using probeDown, probeUp. |
---|
122 | */ |
---|
123 | void disaggCuts (int nstackC, unsigned int probeDir, |
---|
124 | double primalTolerance_, double disaggEffectiveness, |
---|
125 | const OsiSolverInterface &si, |
---|
126 | CglProbingRowCut &rowCut, const int *const stackC, |
---|
127 | const double *const colsol, |
---|
128 | const double *const colUpper, const double *const colLower, |
---|
129 | const double *const saveU, const double *const saveL, |
---|
130 | int *const index, double *const element) |
---|
131 | { |
---|
132 | int pndx = stackC[0] ; |
---|
133 | |
---|
134 | # if CGL_DEBUG > 0 |
---|
135 | std::cout |
---|
136 | << "Entering disaggCuts, probed variable " |
---|
137 | << si.getColName(pndx) << "(" << pndx << ")." << std::endl ; |
---|
138 | const OsiRowCutDebugger *debugger = si.getRowCutDebugger() ; |
---|
139 | // down or up are the only possibilities |
---|
140 | assert((probeDir&(~(probeDown|probeUp))) == 0) ; |
---|
141 | # endif |
---|
142 | |
---|
143 | int plusMinus = 0 ; |
---|
144 | double lu_p = 0.0 ; |
---|
145 | double origsol_p = colsol[pndx] ; |
---|
146 | double origdelta_p = 0.0 ; |
---|
147 | /* |
---|
148 | Set up the quantities that vary depending on whether this is an up or down |
---|
149 | probe: |
---|
150 | * lu_p is the new upper (down probe) or lower (up probe) bound on x<p> |
---|
151 | * plusMinus will be -1 (down probe) or +1 (up probe). |
---|
152 | * origdelta_p will be x*<p>-u<p> (down probe) or l<p>-x*<p> (up probe) |
---|
153 | */ |
---|
154 | if (probeDir == probeDown) { |
---|
155 | lu_p = colUpper[pndx] ; |
---|
156 | plusMinus = -1 ; |
---|
157 | origdelta_p = origsol_p-lu_p ; |
---|
158 | } else { |
---|
159 | lu_p = colLower[pndx] ; |
---|
160 | plusMinus = 1 ; |
---|
161 | origdelta_p = lu_p-origsol_p ; |
---|
162 | } |
---|
163 | |
---|
164 | for (int istackC = nstackC-1 ; istackC > 0 ; istackC--) { |
---|
165 | int icol = stackC[istackC] ; |
---|
166 | double origu_t = saveU[istackC] ; |
---|
167 | double origl_t = saveL[istackC] ; |
---|
168 | double u_t = colUpper[icol] ; |
---|
169 | double sol_t = colsol[icol] ; |
---|
170 | double l_t = colLower[icol] ; |
---|
171 | double boundChange = origu_t-u_t ; |
---|
172 | double coeff = plusMinus*boundChange ; |
---|
173 | /* |
---|
174 | Generate the upper disaggregation cut, if it's violated and the coefficients |
---|
175 | will be reasonable. The effectiveness is the improvement (change) in |
---|
176 | x<pndx>. |
---|
177 | */ |
---|
178 | if (boundChange > 0.0 && origu_t < 1.0e10 && |
---|
179 | (sol_t > u_t+boundChange*origdelta_p+primalTolerance_)) { |
---|
180 | OsiRowCut rc ; |
---|
181 | rc.setLb(-DBL_MAX) ; |
---|
182 | rc.setUb(u_t+coeff*lu_p) ; |
---|
183 | index[0] = icol ; |
---|
184 | element[0] = 1.0 ; |
---|
185 | index[1] = pndx ; |
---|
186 | element[1] = coeff ; |
---|
187 | /* |
---|
188 | Effectiveness is how far x<p> moves. We can save effort by transforming |
---|
189 | sol_p = lu_p + ((u_t-sol_t)/coeff) to delta_p = ((u_t-sol_t)/coeff). |
---|
190 | */ |
---|
191 | double delta_p = ((u_t-sol_t)/coeff) ; |
---|
192 | assert(delta_p > origdelta_p) ; |
---|
193 | rc.setEffectiveness(delta_p-origdelta_p) ; |
---|
194 | if (rc.effectiveness() > disaggEffectiveness) { |
---|
195 | rc.setRow(2,index,element,false) ; |
---|
196 | # if CGL_DEBUG > 0 |
---|
197 | if (debugger) assert(!debugger->invalidCut(rc)); |
---|
198 | # endif |
---|
199 | rowCut.addCutIfNotDuplicate(rc) ; |
---|
200 | } |
---|
201 | } |
---|
202 | /* |
---|
203 | Generate the upper disaggregation cut. |
---|
204 | */ |
---|
205 | boundChange = origl_t-l_t ; |
---|
206 | coeff = plusMinus*boundChange ; |
---|
207 | if (boundChange < 0.0 && origl_t > -1.0e10 && |
---|
208 | (sol_t < l_t+boundChange*origdelta_p-primalTolerance_)) { |
---|
209 | OsiRowCut rc ; |
---|
210 | rc.setLb(l_t+coeff*lu_p) ; |
---|
211 | rc.setUb(DBL_MAX) ; |
---|
212 | index[0] = icol ; |
---|
213 | element[0] = 1.0 ; |
---|
214 | index[1] = pndx ; |
---|
215 | element[1] = coeff ; |
---|
216 | double delta_p = ((sol_t-l_t)/coeff) ; |
---|
217 | assert(delta_p > origdelta_p) ; |
---|
218 | rc.setEffectiveness(delta_p-origdelta_p) ; |
---|
219 | if (rc.effectiveness() > disaggEffectiveness) { |
---|
220 | rc.setRow(2,index,element,false) ; |
---|
221 | # if CGL_DEBUG > 0 |
---|
222 | if (debugger) assert(!debugger->invalidCut(rc)); |
---|
223 | # endif |
---|
224 | rowCut.addCutIfNotDuplicate(rc) ; |
---|
225 | #if 0 |
---|
226 | printf("%d original bounds %g, %g new Lo %g sol= %g int %d sol= %g\n",icol,origl_t,origu_t,colLower[icol],colsol[icol], pndx, colsol[pndx]) ; |
---|
227 | printf("-1.0 * x(%d) + %g * y(%d) <= %g\n", |
---|
228 | icol,boundChange,pndx,rc.ub()) ; |
---|
229 | #endif |
---|
230 | } |
---|
231 | } |
---|
232 | } |
---|
233 | |
---|
234 | # if CGL_DEBUG > 0 |
---|
235 | std::cout |
---|
236 | << "Leaving disaggCuts, " |
---|
237 | << rowCut.numberCuts() << " cuts." << std::endl ; |
---|
238 | # endif |
---|
239 | return ; |
---|
240 | } |
---|
241 | |
---|
242 | |
---|
243 | |
---|
244 | |
---|
245 | /* |
---|
246 | Walk the column for this variable and propagate a bound change on x<j> to |
---|
247 | the row bounds for rows referencing x<j>. Recall that the row bounds are |
---|
248 | |
---|
249 | L<i> = SUM{P}(l<j>) + SUM{M}(u<j>) |
---|
250 | U<i> = SUM{P}(u<j>) + SUM{M}(l<j>) |
---|
251 | |
---|
252 | for negative coefficients M and positive coefficients P. |
---|
253 | The update cases are: |
---|
254 | |
---|
255 | column bound a<ij> row bound |
---|
256 | l<j> >0 L<i> |
---|
257 | l<j> <0 U<i> |
---|
258 | u<j> >0 U<i> |
---|
259 | u<j> <0 L<i> |
---|
260 | |
---|
261 | Movement is used to determine the column bound that's being tightened. It |
---|
262 | should be > 0 to tighten l<j>, movement < 0 to tighten u<j>. Once we've |
---|
263 | selected the row bound to update, the update is always += value*movement. |
---|
264 | */ |
---|
265 | |
---|
266 | bool updateRowBounds (int j, double movement, |
---|
267 | const int *const columnStart, |
---|
268 | const int *const columnLength, |
---|
269 | const int *const row, |
---|
270 | const double *const columnElements, |
---|
271 | const double *const rowUpper, |
---|
272 | const double *const rowLower, |
---|
273 | int &nstackR, int *const stackR, int *const markR, |
---|
274 | double *const minR, double *const maxR, |
---|
275 | double *const saveMin, double *const saveMax) |
---|
276 | { |
---|
277 | bool feasible = true ; |
---|
278 | |
---|
279 | for (int k = columnStart[j] ; k < columnStart[j]+columnLength[j] ; k++) { |
---|
280 | int irow = row[k] ; |
---|
281 | double value = columnElements[k] ; |
---|
282 | double delta = value*movement ; |
---|
283 | /* |
---|
284 | markR is initialised by tighten2. A code of -1 indicates this row should be |
---|
285 | processed; -2 says ignore. Other codes should not happen. |
---|
286 | |
---|
287 | TODO: This assert, and the disabled continue, go back to the change in |
---|
288 | tighten2 where rows with infinite bounds are handed arbitrary |
---|
289 | bounds of 1e60 and labelled suitable for further processing. The |
---|
290 | assert should be removed and the continue clause reinstated. |
---|
291 | -- lh, 101127 -- |
---|
292 | |
---|
293 | TODO: Replace the original assert (!= -2) with a more effective one. |
---|
294 | But don't remove the rest of the code structure. I still think I'll |
---|
295 | end up reintroducing the `ignore' code. -- lh, 101203 -- |
---|
296 | */ |
---|
297 | # if CGL_DEBUG > 0 |
---|
298 | if (markR[irow] < -1 || markR[irow]%10000 >= nstackR) { |
---|
299 | std::cout |
---|
300 | << "Row " << irow << " marked " << markR[irow] |
---|
301 | << "; expected between -1 and " << nstackR << "." << std::endl ; |
---|
302 | |
---|
303 | dump_row(irow,rowLower[irow],rowUpper[irow],minR[irow],maxR[irow], |
---|
304 | 0,false,false,0,0,0,0,0,0,0) ; |
---|
305 | std::cout << std::endl ; |
---|
306 | } |
---|
307 | assert (markR[irow] >= -1 && markR[irow]%10000 < nstackR) ; |
---|
308 | # endif |
---|
309 | if (markR[irow] == -1) { |
---|
310 | stackR[nstackR] = irow ; |
---|
311 | markR[irow] = nstackR ; |
---|
312 | saveMin[nstackR] = minR[irow] ; |
---|
313 | saveMax[nstackR] = maxR[irow] ; |
---|
314 | nstackR++ ; |
---|
315 | # if 0 |
---|
316 | } else if (markR[irow] == -2) { |
---|
317 | continue ; |
---|
318 | # endif |
---|
319 | } |
---|
320 | /* |
---|
321 | LOU_DEBUG: clear count as row bound will change. |
---|
322 | */ |
---|
323 | markR[irow] = markR[irow]%10000 ; |
---|
324 | /* |
---|
325 | Case analysis to tighten the appropriate bound. |
---|
326 | */ |
---|
327 | if (movement > 0.0) { |
---|
328 | if (value > 0.0) { |
---|
329 | if (minR[irow] > -1.0e10) |
---|
330 | minR[irow] += delta ; |
---|
331 | if (minR[irow] > rowUpper[irow]+1.0e-5) { |
---|
332 | feasible = false ; |
---|
333 | break ; |
---|
334 | } |
---|
335 | } else { |
---|
336 | if (maxR[irow] < 1.0e10) |
---|
337 | maxR[irow] += delta ; |
---|
338 | if (maxR[irow] < rowLower[irow]-1.0e-5) { |
---|
339 | feasible = false ; |
---|
340 | break ; |
---|
341 | } |
---|
342 | } |
---|
343 | } else { |
---|
344 | if (value > 0.0) { |
---|
345 | if (maxR[irow] < 1.0e10) |
---|
346 | maxR[irow] += delta ; |
---|
347 | if (maxR[irow] < rowLower[irow]-1.0e-5) { |
---|
348 | feasible = false ; |
---|
349 | break ; |
---|
350 | } |
---|
351 | } else { |
---|
352 | if (minR[irow] > -1.0e10) |
---|
353 | minR[irow] += delta ; |
---|
354 | if (minR[irow] > rowUpper[irow]+1.0e-5) { |
---|
355 | feasible = false ; |
---|
356 | break ; |
---|
357 | } |
---|
358 | } |
---|
359 | } |
---|
360 | } |
---|
361 | return (feasible) ; |
---|
362 | } |
---|
363 | |
---|
364 | |
---|
365 | /* |
---|
366 | jjf: clean up djs and solution |
---|
367 | |
---|
368 | Calculate domains (u<j>-l<j>) for variables, and potential swings in rows. Do |
---|
369 | some cleanup of reduced costs. |
---|
370 | |
---|
371 | If CGL_DEBUG > 0, we'll also do some matrix consistency checks. |
---|
372 | */ |
---|
373 | void groomSoln (double direction, double primalTolerance, double *const djs, |
---|
374 | const double *const colLower, double *const colsol, |
---|
375 | const double *const colUpper, double *const columnGap, |
---|
376 | const CoinPackedMatrix *const rowCopy, |
---|
377 | const int *const rowStartPos, |
---|
378 | double *const largestNegativeInRow, |
---|
379 | double *const largestPositiveInRow) |
---|
380 | { |
---|
381 | int nRows = rowCopy->getNumRows() ; |
---|
382 | int nCols = rowCopy->getNumCols() ; |
---|
383 | const int *rowStart = rowCopy->getVectorStarts() ; |
---|
384 | const int *column = rowCopy->getIndices() ; |
---|
385 | const double *rowElements = rowCopy->getElements() ; |
---|
386 | |
---|
387 | # ifndef NDEBUG |
---|
388 | // Note that rowLength is only used in the assert |
---|
389 | const int *rowLength = rowCopy->getVectorLengths() ; |
---|
390 | # endif |
---|
391 | |
---|
392 | /* |
---|
393 | Find the largest potential postive and negative swing in a row, where |
---|
394 | swing is defined as a<ij>(u<j>-l<j>). Check correctness of rowStartPos if |
---|
395 | we're debugging. |
---|
396 | */ |
---|
397 | for (int i = 0 ; i < nRows ; i++) { |
---|
398 | |
---|
399 | assert (rowStart[i]+rowLength[i] == rowStart[i+1]) ; |
---|
400 | |
---|
401 | int kk ; |
---|
402 | double value ; |
---|
403 | |
---|
404 | #if CGL_DEBUG > 0 |
---|
405 | // Check correctness of rowStartPos |
---|
406 | for (kk = rowStart[i] ; kk < rowStart[i+1] ; kk++) { |
---|
407 | value = rowElements[kk] ; |
---|
408 | if (value > 0.0) |
---|
409 | break ; |
---|
410 | } |
---|
411 | assert (rowStartPos[i] == kk) ; |
---|
412 | #endif |
---|
413 | |
---|
414 | value = 0.0 ; |
---|
415 | for (kk = rowStart[i] ; kk < rowStartPos[i] ; kk++) { |
---|
416 | int iColumn = column[kk] ; |
---|
417 | double gap = CoinMin(1.0e100,colUpper[iColumn]-colLower[iColumn]) ; |
---|
418 | value = CoinMin(value,gap*rowElements[kk]) ; |
---|
419 | } |
---|
420 | largestNegativeInRow[i] = value ; |
---|
421 | value = 0.0 ; |
---|
422 | for ( ; kk < rowStart[i+1] ; kk++) { |
---|
423 | int iColumn = column[kk] ; |
---|
424 | double gap = CoinMin(1.0e100,colUpper[iColumn]-colLower[iColumn]) ; |
---|
425 | value = CoinMax(value,gap*rowElements[kk]) ; |
---|
426 | } |
---|
427 | largestPositiveInRow[i] = value ; |
---|
428 | } |
---|
429 | /* |
---|
430 | Convert reduced costs to minimisation and clean up any that are on the wrong |
---|
431 | side of zero. |
---|
432 | |
---|
433 | TODO: Seems like we could move this ahead of the previous loops and then |
---|
434 | use columnGap[i] to calculate largest[Negative,Positive]InRow. |
---|
435 | -- lh, 101125 -- |
---|
436 | */ |
---|
437 | for (int i = 0 ; i < nCols ; ++i) { |
---|
438 | double djValue = djs[i]*direction ; |
---|
439 | double gap = colUpper[i]-colLower[i] ; |
---|
440 | if (gap > 1.0e-8) { |
---|
441 | if (colsol[i] < colLower[i]+primalTolerance) { |
---|
442 | colsol[i] = colLower[i] ; |
---|
443 | djs[i] = CoinMax(0.0,djValue) ; |
---|
444 | } else if (colsol[i] > colUpper[i]-primalTolerance) { |
---|
445 | colsol[i] = colUpper[i] ; |
---|
446 | djs[i] = CoinMin(0.0,djValue) ; |
---|
447 | } else { |
---|
448 | djs[i] = 0.0 ; |
---|
449 | } |
---|
450 | } |
---|
451 | columnGap[i] = gap-primalTolerance ; |
---|
452 | } |
---|
453 | |
---|
454 | return ; |
---|
455 | } |
---|
456 | |
---|
457 | /* |
---|
458 | Given that the variable being probed has been discovered to be monotone, |
---|
459 | this method will save the resulting bound changes for integer variables as |
---|
460 | column cuts. |
---|
461 | |
---|
462 | One could do the same for continuous variables, but at present we don't. |
---|
463 | Presumably we need to test against origColLower, origColUpper (the |
---|
464 | bounds held in the si) to decide if we have real cuts on the original |
---|
465 | problem. One could argue that we could test against saveL and saveU and |
---|
466 | save the trouble of repeatedly recording the same column cut. |
---|
467 | |
---|
468 | Index and elements are scratch arrays. |
---|
469 | |
---|
470 | Returns true if any of the bound improvements are actually cuts, false |
---|
471 | otherwise. |
---|
472 | */ |
---|
473 | bool monotoneActions (double primalTolerance_, |
---|
474 | const OsiSolverInterface &si, |
---|
475 | OsiCuts &cs, |
---|
476 | int nstackC, const int *const stackC, |
---|
477 | const char *const intVar, |
---|
478 | const double *const colLower, |
---|
479 | const double *const colsol, |
---|
480 | const double *const colUpper, |
---|
481 | int *const index, double *const element) |
---|
482 | { |
---|
483 | OsiColCut cc ; |
---|
484 | bool anyColumnCuts = false ; |
---|
485 | int nTot = 0 ; |
---|
486 | int nFix = 0 ; |
---|
487 | const double *origColLower = si.getColLower() ; |
---|
488 | const double *origColUpper = si.getColUpper() ; |
---|
489 | # if CGL_DEBUG > 0 |
---|
490 | int probedVar = stackC[0] ; |
---|
491 | std::cout |
---|
492 | << " " << si.getColName(probedVar) << " (" << probedVar |
---|
493 | << ") monotone [" << colLower[probedVar] << ", " |
---|
494 | << colUpper[probedVar] << "] from " << colsol[probedVar] |
---|
495 | << "." << std::endl ; |
---|
496 | # endif |
---|
497 | /* |
---|
498 | Stash the improved column bounds, u<j> first, then l<j>. If any of them |
---|
499 | actually cut off the current solution, indicate that we have cuts. |
---|
500 | */ |
---|
501 | bool ifCut = false ; |
---|
502 | for (int istackC = 0 ; istackC < nstackC ; istackC++) { |
---|
503 | int icol = stackC[istackC] ; |
---|
504 | if (intVar[icol]) { |
---|
505 | if (colUpper[icol] < origColUpper[icol]-1.0e-4) { |
---|
506 | element[nFix] = colUpper[icol] ; |
---|
507 | index[nFix++] = icol ; |
---|
508 | if (colsol[icol] > colUpper[icol]+primalTolerance_) { |
---|
509 | ifCut = true ; |
---|
510 | anyColumnCuts = true ; |
---|
511 | } |
---|
512 | } |
---|
513 | } |
---|
514 | } |
---|
515 | if (nFix) { |
---|
516 | nTot = nFix ; |
---|
517 | cc.setUbs(nFix,index,element) ; |
---|
518 | nFix = 0 ; |
---|
519 | } |
---|
520 | for (int istackC = 0 ; istackC < nstackC ; istackC++) { |
---|
521 | int icol = stackC[istackC] ; |
---|
522 | if (intVar[icol]) { |
---|
523 | if (colLower[icol] > origColLower[icol]+1.0e-4) { |
---|
524 | element[nFix] = colLower[icol] ; |
---|
525 | index[nFix++] = icol ; |
---|
526 | if (colsol[icol] < colLower[icol]-primalTolerance_) { |
---|
527 | ifCut = true ; |
---|
528 | anyColumnCuts = true ; |
---|
529 | } |
---|
530 | } |
---|
531 | } |
---|
532 | } |
---|
533 | if (nFix) { |
---|
534 | nTot += nFix ; |
---|
535 | cc.setLbs(nFix,index,element) ; |
---|
536 | } |
---|
537 | if (nTot) { |
---|
538 | if (ifCut) { |
---|
539 | cc.setEffectiveness(100.0) ; |
---|
540 | } else { |
---|
541 | cc.setEffectiveness(1.0e-5) ; |
---|
542 | } |
---|
543 | #if CGL_DEBUG > 0 |
---|
544 | checkBounds(si,cc) ; |
---|
545 | #endif |
---|
546 | cs.insert(cc) ; |
---|
547 | } |
---|
548 | |
---|
549 | return (anyColumnCuts) ; |
---|
550 | } |
---|
551 | |
---|
552 | void clearStacks (double primalTolerance_, |
---|
553 | int nstackC, int *const stackC, int *const markC, |
---|
554 | const double *const colLower, const double *const colUpper, |
---|
555 | int nstackR, int *const stackR, int *const markR) |
---|
556 | { |
---|
557 | /* |
---|
558 | Clear off the bound change markers, except for fixed variables. We'll be |
---|
559 | moving on to another variable. |
---|
560 | */ |
---|
561 | for (int istackC = 0 ; istackC < nstackC ; istackC++) { |
---|
562 | int icol = stackC[istackC] ; |
---|
563 | if (colUpper[icol]-colLower[icol] > primalTolerance_) { |
---|
564 | markC[icol] &= ~(tightenLower|tightenUpper) ; |
---|
565 | } else { |
---|
566 | markC[icol] = tightenLower|tightenUpper ; |
---|
567 | } |
---|
568 | } |
---|
569 | for (int istackR = 0 ; istackR < nstackR ; istackR++) { |
---|
570 | int irow = stackR[istackR] ; |
---|
571 | markR[irow] = -1 ; |
---|
572 | } |
---|
573 | } |
---|
574 | |
---|
575 | |
---|
576 | |
---|
577 | /* |
---|
578 | This method examines the rows on stackR looking for redundant rows that |
---|
579 | consist solely of singleton variables (i.e., variables that appear in just |
---|
580 | this constraint). It then looks to see if any of the variables in the row |
---|
581 | can be fixed at bound. |
---|
582 | |
---|
583 | Consider the case where U<i> < b<i> and all unfixed variables in the row |
---|
584 | are singletons (occur in no other constraint). Given x<j> is a singleton, |
---|
585 | the reduced cost is c<j> - ya<j> = c<j> - ye<i> = c<j> - y<i>. But if |
---|
586 | U<i> < b<i>, the constraint is loose by definition, hence y<i> = 0 and |
---|
587 | the reduced cost is c<j>. If a<ij> > 0, x<j> should be u<j> in U<i>. If |
---|
588 | c<j> < 0, minimising will tend to drive x<j> to u<j>. This cannot violate |
---|
589 | constraint i (U<i> < b<i>) and (since x<j> is a singleton) will have no |
---|
590 | effect elsewhere. Hence we can fix x<j> at u<j>. |
---|
591 | |
---|
592 | Do the case analysis, and what you find is that against U<i> we want |
---|
593 | a<ij> > 0 ==> c<j> < 0 ==> push x<j> to u<j> |
---|
594 | a<ij> < 0 ==> c<j> > 0 ==> push x<j> to l<j> |
---|
595 | and against L<i> we want |
---|
596 | a<ij> > 0 ==> c<j> > 0 ==> push x<j> to l<j> |
---|
597 | a<ij> < 0 ==> c<j> < 0 ==> push x<j> to u<j> |
---|
598 | |
---|
599 | Extend it one more time by taking the objective direction (max/min) into |
---|
600 | account (dir = 1.0 for min, -1.0 for max) and you have |
---|
601 | against U<i> ==> a<ij>*c<j>*dir < 0 |
---|
602 | against L<i> ==> a<ij>*c<j>*dir > 0 |
---|
603 | |
---|
604 | John's original comment for this code was |
---|
605 | // also see if singletons can go to good objective |
---|
606 | // Taken out as should be found elsewhere |
---|
607 | // and has to be original column length |
---|
608 | but he reinstated it. Apparently there were cases where fixing the probing |
---|
609 | variable was required to satisfy the condition that all unfixed variables be |
---|
610 | singletons. Enough of them to justify reinstating this code. |
---|
611 | */ |
---|
612 | |
---|
613 | void singletonRows (int jProbe, double primalTolerance_, |
---|
614 | const OsiSolverInterface &si, |
---|
615 | const CoinPackedMatrix *rowCopy, |
---|
616 | int *const markC, |
---|
617 | int &nstackC, int *const stackC, |
---|
618 | double *const saveL, double *const saveU, |
---|
619 | double *const colUpper, double *const colLower, |
---|
620 | double *const colGap, |
---|
621 | const int nstackR, const int *const stackR, |
---|
622 | const double *const rowUpper, |
---|
623 | const double *const rowLower, |
---|
624 | const double *const maxR, const double *const minR) |
---|
625 | { |
---|
626 | /* |
---|
627 | Unpack a few vectors from the row-major matrix. |
---|
628 | */ |
---|
629 | const CoinBigIndex *rowStart = rowCopy->getVectorStarts() ; |
---|
630 | const int *column = rowCopy->getIndices() ; |
---|
631 | const double *rowElements = rowCopy->getElements() ; |
---|
632 | const int nCols = rowCopy->getNumCols() ; |
---|
633 | /* |
---|
634 | `Singleton' must be based on the column length in the original system. |
---|
635 | */ |
---|
636 | const double *objective = si.getObjCoefficients() ; |
---|
637 | const int *columnLengths = si.getMatrixByCol()->getVectorLengths() ; |
---|
638 | const double objSense = si.getObjSense() ; |
---|
639 | /* |
---|
640 | Open a loop to work through the rows on stackR. |
---|
641 | */ |
---|
642 | for (int istackR = 0 ; istackR < nstackR ; istackR++) { |
---|
643 | int i = stackR[istackR] ; |
---|
644 | /* |
---|
645 | Check the gaps. If the constraint is potentially tight in both directions, |
---|
646 | there's nothing more to do here. |
---|
647 | */ |
---|
648 | const double uGap = rowUpper[i]-maxR[i] ; |
---|
649 | const double lGap = minR[i]-rowLower[i] ; |
---|
650 | if (uGap < primalTolerance_ && lGap < primalTolerance_) continue ; |
---|
651 | /* |
---|
652 | Scan the row to see if it meets the `all singletons' condition. Again, |
---|
653 | if this fails, there's nothing more to be done. |
---|
654 | |
---|
655 | Note that the original code didn't check the probing variable x<p>, |
---|
656 | because if you're probing a binary variable it's fixed. But for general |
---|
657 | integer variables a probe does not in general fix the variable. So we |
---|
658 | check all variables. |
---|
659 | |
---|
660 | We should not be executing this method if we have prima facie infeasibility. |
---|
661 | */ |
---|
662 | bool canFix = true ; |
---|
663 | for (int kk = rowStart[i] ; kk < rowStart[i+1] ; kk++) { |
---|
664 | int j = column[kk] ; |
---|
665 | assert(colUpper[j]-colLower[j] > -primalTolerance_) ; |
---|
666 | if (colUpper[j] > colLower[j] && columnLengths[j] != 1) { |
---|
667 | canFix = false ; |
---|
668 | break ; |
---|
669 | } |
---|
670 | } |
---|
671 | if (!canFix) continue ; |
---|
672 | /* |
---|
673 | If we've passed the tests, we can look for variables suitable to drive to |
---|
674 | bound. Work against U<i> first. We're looking for variables with a<ij> > 0 |
---|
675 | that will be naturally driven to u<j>, or variables with a<ij> < 0 that will |
---|
676 | be naturally driven to l<j>. |
---|
677 | |
---|
678 | Don't overwrite the saved bounds if we've |
---|
679 | tightened this variable already! |
---|
680 | */ |
---|
681 | if (uGap > primalTolerance_) { |
---|
682 | for (int kk = rowStart[i] ; kk < rowStart[i+1] ; kk++) { |
---|
683 | int j = column[kk] ; |
---|
684 | const double lj = colLower[j] ; |
---|
685 | const double uj = colUpper[j] ; |
---|
686 | if (uj > lj) { |
---|
687 | double value = rowElements[kk] ; |
---|
688 | if (objSense*objective[j]*value < 0.0) |
---|
689 | { assert(jProbe != j) ; |
---|
690 | if (!(markC[j]&(tightenLower|tightenUpper))) { |
---|
691 | stackC[nstackC] = j ; |
---|
692 | saveL[nstackC] = lj ; |
---|
693 | saveU[nstackC] = uj ; |
---|
694 | } |
---|
695 | assert(saveU[nstackC] > saveL[nstackC]) ; |
---|
696 | if (value > 0.0) { |
---|
697 | colLower[j] = uj ; |
---|
698 | } else { |
---|
699 | colUpper[j] = lj ; |
---|
700 | } |
---|
701 | markC[j] |= tightenLower|tightenUpper ; |
---|
702 | colGap[j] = -primalTolerance_ ; |
---|
703 | assert(nstackC < nCols) ; |
---|
704 | nstackC++ ; |
---|
705 | } |
---|
706 | } |
---|
707 | } |
---|
708 | } |
---|
709 | /* |
---|
710 | And now the identical code, except that we're working against L<i>, hence |
---|
711 | the sense of the elibigility test is reversed and we want variables with |
---|
712 | a<ij> > 0 that will be naturally driven to l<j>, or variables with |
---|
713 | a<ij> < 0 that will be naturally driven to u<j>. |
---|
714 | */ |
---|
715 | if (lGap > primalTolerance_) { |
---|
716 | for (int kk = rowStart[i] ; kk < rowStart[i+1] ; kk++) { |
---|
717 | int j = column[kk] ; |
---|
718 | const double lj = colLower[j] ; |
---|
719 | const double uj = colUpper[j] ; |
---|
720 | if (uj > lj) { |
---|
721 | double value = rowElements[kk] ; |
---|
722 | if (objSense*objective[j]*value > 0.0) |
---|
723 | { assert(jProbe != j) ; |
---|
724 | if (!(markC[j]&(tightenLower|tightenUpper))) { |
---|
725 | stackC[nstackC] = j ; |
---|
726 | saveL[nstackC] = lj ; |
---|
727 | saveU[nstackC] = uj ; |
---|
728 | } |
---|
729 | assert(saveU[nstackC] > saveL[nstackC]) ; |
---|
730 | if (value < 0.0) { |
---|
731 | colLower[j] = uj ; |
---|
732 | } else { |
---|
733 | colUpper[j] = lj ; |
---|
734 | } |
---|
735 | markC[j] |= tightenLower|tightenUpper ; |
---|
736 | colGap[j] = -primalTolerance_ ; |
---|
737 | assert(nstackC < nCols) ; |
---|
738 | nstackC++ ; |
---|
739 | } |
---|
740 | } |
---|
741 | } |
---|
742 | } |
---|
743 | } |
---|
744 | return ; |
---|
745 | } |
---|
746 | |
---|
747 | |
---|
748 | /* |
---|
749 | Generate coefficient strengthening cuts. |
---|
750 | |
---|
751 | Assume we have a binary probe variable x<p> with a<ip> > 0. Assume a |
---|
752 | down probe. Assume that U<i> > b<i> before the probe, but now U'<i> < b<i> |
---|
753 | (i.e., the constraint is redundant against b<i> after the down probe forces |
---|
754 | x<p> to 0 and reduces U<i> to U'<i>). We would like to have the following in |
---|
755 | a strengthened constraint a'<ip>x<p> + (stuff) <= b'<i>: |
---|
756 | * When x<p> = 0, b'<i> = U'<i> (the lhs can't be larger than U'<i>) |
---|
757 | * When x<p> = 1, b'<i> = b<i> (the original value) |
---|
758 | Define delta = b<i> - U'<i>, b'<i> = b<i>-delta, and a'<ip> = a<ip>-delta. |
---|
759 | When x<p> = 1, the delta term on each side cancels and we're left with the |
---|
760 | original constraint. When x<p> = 0, the rhs tightens to |
---|
761 | b'<i> = b<i>+delta<i> = U'<i>. |
---|
762 | |
---|
763 | For an honest derivation that works for binary and general integer |
---|
764 | variables, see the typeset documentation. As you'd expect, there are four |
---|
765 | cases (a<ip> < 0, a<ip> > 0) X (up probe, down probe). Assume that a down |
---|
766 | probe bounds x<p> as [l<p>, u'<p>] and an up probe bounds x<p> as [l'<p>, |
---|
767 | u<p>]. Then the summary (omitting subscripts to fit the line) is: |
---|
768 | |
---|
769 | a<ip> dir delta blow' a'<ip> b' |
---|
770 | |
---|
771 | >0 d b - U' a-delta b-(u'+1)delta |
---|
772 | >0 u blow - L' blow+(l'-1)delta a+delta |
---|
773 | <0 d b - L' blow-(u'+1)delta a-delta |
---|
774 | <0 u blow - U' a+delta b+(l'-1)delta |
---|
775 | |
---|
776 | In the code, negating delta for down probes unifies the up and down math. |
---|
777 | |
---|
778 | TODO: In the original code, this method will not execute if column cuts |
---|
779 | are present (i.e., bounds were tightened so that some x*<j> may not |
---|
780 | be feasible). The current code retains those guards. Clearly it's |
---|
781 | a problem if x*<p> isn't feasible --- we can't generate a useful |
---|
782 | cut over floor(x*<p>), ceil(x*<p>) if the interval is outside the |
---|
783 | polytope. It's not so clear this is an obstacle for other variables, |
---|
784 | except that the lhs activity may be inaccurate. |
---|
785 | */ |
---|
786 | |
---|
787 | #define STRENGTHEN_PRINT |
---|
788 | void strengthenCoeff ( |
---|
789 | int jProbe, unsigned int probeDir, |
---|
790 | double primalTolerance_, |
---|
791 | bool justReplace, bool canReplace, |
---|
792 | double needEffectiveness, |
---|
793 | const OsiSolverInterface &si, |
---|
794 | CglProbingRowCut &rowCut, |
---|
795 | const CoinPackedMatrix *rowCopy, |
---|
796 | double *const colUpper, double *const colLower, |
---|
797 | const double *const colsol, |
---|
798 | const int nstackR, const int *const stackR, |
---|
799 | const double *const rowUpper, |
---|
800 | const double *const rowLower, |
---|
801 | const double *const maxR, |
---|
802 | const double *const minR, |
---|
803 | const int *const realRows, |
---|
804 | double *const element, |
---|
805 | int *const index, |
---|
806 | CglTreeInfo *const info |
---|
807 | ) |
---|
808 | { |
---|
809 | # if CGL_DEBUG > 0 |
---|
810 | std::cout |
---|
811 | << "Entering strengthenCoeff, probing " |
---|
812 | << si.getColName(jProbe) << "(" << jProbe << ") " |
---|
813 | << ((probeDir == probeDown)?"down":"up") << "." << std::endl ; |
---|
814 | const OsiRowCutDebugger *debugger = si.getRowCutDebugger() ; |
---|
815 | # endif |
---|
816 | /* |
---|
817 | Define magic numbers up front. |
---|
818 | * Limits on `reasonable' coefficients. We don't want coefficients |
---|
819 | outside this range. |
---|
820 | * `Infinity' for the rhs. This'll get fixed later. |
---|
821 | * Effectiveness in the case that a<ip> = 0 and a'<ip> != 0. |
---|
822 | All copied from the original code. |
---|
823 | */ |
---|
824 | const double bigCoeff = 1.0e8 ; |
---|
825 | const double tinyCoeff = 1.0e-12 ; |
---|
826 | const double rhsInf = 1.0e20 ; |
---|
827 | const double aipZeroEffectiveness = 1.0e-7 ; |
---|
828 | /* |
---|
829 | Unpack a few vectors from the row-major matrix. |
---|
830 | */ |
---|
831 | const CoinBigIndex *rowStart = rowCopy->getVectorStarts() ; |
---|
832 | const int *column = rowCopy->getIndices() ; |
---|
833 | const double *rowElements = rowCopy->getElements() ; |
---|
834 | /* |
---|
835 | Open up an outer loop to walk stackR and look for interesting constraints. |
---|
836 | */ |
---|
837 | for (int istackR = 0 ; istackR < nstackR ; istackR++) { |
---|
838 | int irow = stackR[istackR] ; |
---|
839 | double bi = rowUpper[irow] ; |
---|
840 | double blowi = rowLower[irow] ; |
---|
841 | /* |
---|
842 | We can't get anywhere unless probing has made one end or the other of the |
---|
843 | constraint redundant. If there's no gap, we can't do anything. |
---|
844 | */ |
---|
845 | double uGap = bi-maxR[irow] ; |
---|
846 | double lGap = blowi-minR[irow] ; |
---|
847 | if (uGap < primalTolerance_ && -lGap < primalTolerance_) continue ; |
---|
848 | bool isRangeCon = ((blowi > -rhsInf) && (bi < rhsInf)) ; |
---|
849 | /* |
---|
850 | We'll need the lhs activity, excluding the probed variable, to evaluate the |
---|
851 | effectiveness of the cut. Extract the coefficient of the probe variable |
---|
852 | while we're here; we need it to determine which case we're working. |
---|
853 | |
---|
854 | TODO: As noted at the head of the routine, we could correct for column cuts |
---|
855 | by checking bounds when calculating the sum below, if it was worth the |
---|
856 | effort. -- lh, 110113 -- |
---|
857 | */ |
---|
858 | double sum = 0.0 ; |
---|
859 | double aip = 0.0 ; |
---|
860 | bool aipNonZero = false ; |
---|
861 | for (CoinBigIndex kk = rowStart[irow] ; kk < rowStart[irow+1] ; kk++) { |
---|
862 | int k = column[kk] ; |
---|
863 | double aik = rowElements[kk] ; |
---|
864 | if (k == jProbe) { |
---|
865 | aipNonZero = true ; |
---|
866 | aip = aik ; |
---|
867 | } else { |
---|
868 | sum += aik*colsol[k] ; |
---|
869 | } |
---|
870 | } |
---|
871 | /* |
---|
872 | Now figure out which case we're in and do some setup. At the end, see if the |
---|
873 | coefficient will have reasonable magnitude. If not, move on to the next |
---|
874 | iteration. |
---|
875 | */ |
---|
876 | double delta ; |
---|
877 | double deltaMul ; |
---|
878 | bool revisebi = true ; |
---|
879 | if (probeDir == probeDown) { |
---|
880 | deltaMul = colUpper[jProbe]+1.0 ; |
---|
881 | if (aip >= 0) { |
---|
882 | delta = -uGap ; |
---|
883 | } else { |
---|
884 | delta = -lGap ; |
---|
885 | revisebi = false ; |
---|
886 | } |
---|
887 | } else { |
---|
888 | deltaMul = colLower[jProbe]-1.0 ; |
---|
889 | if (aip >= 0) { |
---|
890 | delta = lGap ; |
---|
891 | revisebi = false ; |
---|
892 | } else { |
---|
893 | delta = uGap ; |
---|
894 | } |
---|
895 | } |
---|
896 | if (CoinAbs(delta) < primalTolerance_ || CoinAbs(delta) > bigCoeff) |
---|
897 | continue ; |
---|
898 | /* |
---|
899 | Now decide if we have something that cuts off the current solution. Augment |
---|
900 | the lhs activity with the contribution for a'<ip>x*<p>, calculate the new |
---|
901 | rhs value, and compare. |
---|
902 | As an alternate measure of effectiveness, consider the gap between the |
---|
903 | current activity and the revised lhs bound, normalised by the gap between |
---|
904 | the original rhs and the revised lhs bound. |
---|
905 | |
---|
906 | We'll also generate the cut if we can strengthen it in place and we're not |
---|
907 | dealing with a range constraint. (The reason for excluding range constraints |
---|
908 | is that we might try to alter a<ip> against both b<i> and blow<i>. That way |
---|
909 | lies madness.) |
---|
910 | |
---|
911 | TODO: It's not clear to me how the first two effectiveness calculations |
---|
912 | relate to one another. Think more on this -- lh, 110115 -- |
---|
913 | */ |
---|
914 | double aipPrime = aip+delta ; |
---|
915 | bool aipPrimeNonZero = true ; |
---|
916 | if (CoinAbs(aipPrime) < tinyCoeff) { |
---|
917 | aipPrime = 0.0 ; |
---|
918 | aipPrimeNonZero = false ; |
---|
919 | } |
---|
920 | double xp = colsol[jProbe] ; |
---|
921 | sum += aipPrime*xp ; |
---|
922 | double biPrime = DBL_MAX ; |
---|
923 | double blowiPrime = -DBL_MAX ; |
---|
924 | double effectiveness = 0.0 ; |
---|
925 | bool genCut = false ; |
---|
926 | if (revisebi) { |
---|
927 | biPrime = rowUpper[irow]+deltaMul*delta ; |
---|
928 | effectiveness = sum-biPrime ; |
---|
929 | effectiveness = CoinMax(effectiveness,(sum-maxR[irow])/uGap) ; |
---|
930 | } else { |
---|
931 | blowiPrime = rowLower[irow]+deltaMul*delta ; |
---|
932 | effectiveness = blowiPrime-sum ; |
---|
933 | effectiveness = CoinMax(effectiveness,(minR[irow]-sum)/lGap) ; |
---|
934 | } |
---|
935 | if (!aipNonZero && aipPrimeNonZero) |
---|
936 | effectiveness = CoinMax(effectiveness,aipZeroEffectiveness) ; |
---|
937 | if (effectiveness > needEffectiveness) genCut = true ; |
---|
938 | if (info->strengthenRow && !isRangeCon) genCut = true ; |
---|
939 | /* |
---|
940 | Are we going to generate the cut? If not, on to the next iteration. |
---|
941 | */ |
---|
942 | if (!genCut) continue ; |
---|
943 | /* |
---|
944 | Generate the coefficients. Copy whatever's present and plug in aipPrime at |
---|
945 | the end. |
---|
946 | */ |
---|
947 | int n = 0 ; |
---|
948 | int aip_n = -1 ; |
---|
949 | for (CoinBigIndex kk = rowStart[irow] ; kk < rowStart[irow+1] ; kk++) { |
---|
950 | int k = column[kk] ; |
---|
951 | double aik = rowElements[kk] ; |
---|
952 | if (k == jProbe) |
---|
953 | aip_n = k ; |
---|
954 | index[n] = k ; |
---|
955 | element[n++] = aik ; |
---|
956 | } |
---|
957 | if (aip_n < 0) { |
---|
958 | aip_n = n ; |
---|
959 | n++ ; |
---|
960 | } |
---|
961 | index[aip_n] = jProbe ; |
---|
962 | element[aip_n] = aipPrime ; |
---|
963 | /* |
---|
964 | Fill in a cut structure with the cut. |
---|
965 | */ |
---|
966 | OsiRowCut rc ; |
---|
967 | rc.setLb(blowiPrime) ; |
---|
968 | rc.setUb(biPrime) ; |
---|
969 | rc.setEffectiveness(effectiveness) ; |
---|
970 | rc.setRow(n,index,element,false) ; |
---|
971 | # if CGL_DEBUG > 0 |
---|
972 | if (debugger) assert(!debugger->invalidCut(rc)); |
---|
973 | # endif |
---|
974 | # ifdef STRENGTHEN_PRINT |
---|
975 | if (canReplace && !isRangeCon) { |
---|
976 | printf("Strengthen Cut (1):\n") ; |
---|
977 | dump_row(irow,rc.lb(),rc.ub(), |
---|
978 | nan(""),nan(""),&si,true,false,4, |
---|
979 | n,index,element, |
---|
980 | 1.0e-10,colLower,colUpper) ; |
---|
981 | printf("Original Row:\n") ; |
---|
982 | int k = rowStart[irow] ; |
---|
983 | dump_row(irow,rowLower[irow],rowUpper[irow], |
---|
984 | minR[irow],maxR[irow],&si,true,false,4, |
---|
985 | rowStart[irow+1]-k,&column[k],&rowElements[k], |
---|
986 | 1.0e-10,colLower,colUpper) ; |
---|
987 | } |
---|
988 | # endif |
---|
989 | /* |
---|
990 | If we're in preprocessing, we might try to simply replace the existing |
---|
991 | constraint (justReplace = true; preprocessing is the only place this will |
---|
992 | happen). Otherwise, drop the cut into the cut set. |
---|
993 | |
---|
994 | realRows comes in as a parameter. This is the translation array created if |
---|
995 | we modified the constraint system during preprocessing in gutsOfGenerateCuts. |
---|
996 | |
---|
997 | Effectiveness, if we're strengthening in place, seems to be absolute size of |
---|
998 | coefficients; smaller is better. (Why?) |
---|
999 | |
---|
1000 | TODO: It seems to me that we can get here with |
---|
1001 | justReplace = true and canReplace = false, and this condition is |
---|
1002 | constant over an iteration of the way loop. In other words, we've done |
---|
1003 | all the work to generate this cut and we knew before we started that |
---|
1004 | we would discard it here. -- lh, 110107 -- |
---|
1005 | Put in an assert to see if we ever do all the work, only to reject |
---|
1006 | the cut. -- lh, 110115 -- |
---|
1007 | */ |
---|
1008 | int realRow = (canReplace && !isRangeCon)?(irow):(-1) ; |
---|
1009 | if (realRows && realRow >= 0) |
---|
1010 | realRow = realRows[realRow] ; |
---|
1011 | if (!justReplace) { |
---|
1012 | rowCut.addCutIfNotDuplicate(rc,realRow) ; |
---|
1013 | } else if (realRow >= 0) { |
---|
1014 | double effectiveness = 0.0 ; |
---|
1015 | for (int i = 0 ; i < n ; i++) |
---|
1016 | effectiveness += CoinAbs(element[i]) ; |
---|
1017 | if (!info->strengthenRow[realRow] || |
---|
1018 | info->strengthenRow[realRow]->effectiveness() > effectiveness) { |
---|
1019 | delete info->strengthenRow[realRow] ; |
---|
1020 | rc.setEffectiveness(effectiveness) ; |
---|
1021 | info->strengthenRow[realRow] = rc.clone() ; |
---|
1022 | } else { |
---|
1023 | # if CGL_DEBUG > 0 |
---|
1024 | std::cout |
---|
1025 | << "INPLACE: rejected on cut coeff magnitude." << std::endl ; |
---|
1026 | # endif |
---|
1027 | } |
---|
1028 | } else { |
---|
1029 | # if CGL_DEBUG > 0 |
---|
1030 | std::cout |
---|
1031 | << "INPLACE: rejected because no real row." << std::endl ; |
---|
1032 | # endif |
---|
1033 | } |
---|
1034 | } |
---|
1035 | |
---|
1036 | # if CGL_DEBUG > 0 |
---|
1037 | std::cout |
---|
1038 | << "Leaving strengthenCoeff, " |
---|
1039 | << rowCut.numberCuts() << " cuts." << std::endl ; |
---|
1040 | # endif |
---|
1041 | |
---|
1042 | return ; |
---|
1043 | } |
---|
1044 | |
---|
1045 | |
---|
1046 | |
---|
1047 | |
---|
1048 | // ========================================================= |
---|
1049 | |
---|
1050 | |
---|
1051 | |
---|
1052 | /* |
---|
1053 | jjf: Does probing and adding cuts |
---|
1054 | |
---|
1055 | Note that this method is heavily commented and instrumented in CbcAnn. |
---|
1056 | |
---|
1057 | It would appear that probe() has received more attention that probeClique or |
---|
1058 | probeSlack. Neither of them has the ONE_ARRAY option. |
---|
1059 | */ |
---|
1060 | |
---|
1061 | /* |
---|
1062 | We're going to probe integer variables, and we're interested in propagating |
---|
1063 | the effect of each probe (force a variable up or down). |
---|
1064 | |
---|
1065 | The bulk of the method consists of three nested loops: |
---|
1066 | * PASS LOOP: makes multiple passes, probing all variables in lookedAt_ |
---|
1067 | * LOOKEDAT LOOP: iterates over the variables in lookedAt_ |
---|
1068 | * PROBE LOOP: probes down/up/down for each variable. |
---|
1069 | |
---|
1070 | The body of the probe loop runs a bit over 2600 lines and propagates |
---|
1071 | the effect of forcing the probed variable down or up. |
---|
1072 | |
---|
1073 | * The start of the body updates the row bounds of affected rows, then |
---|
1074 | initialises stackC with the probed variable. |
---|
1075 | |
---|
1076 | * The middle is a 1,000 line loop that does the bulk of propagation. |
---|
1077 | At the start there's some dead code (PROBING100) associated with |
---|
1078 | CglTreeProbingInfo. The purpose remains to be determined. This is |
---|
1079 | followed by six specialised replications of the same functionality: |
---|
1080 | walk the column of a variable from stackC, and for each affected |
---|
1081 | row, try to tighten the bounds of other variables in the row. If |
---|
1082 | we successfully tighten bounds on a variable, walk the column of |
---|
1083 | that variable, tighten the bounds of affected rows, and place the |
---|
1084 | variable on stackC. |
---|
1085 | |
---|
1086 | * The end of the body deals with the result of the probe. At the end |
---|
1087 | of each iteration, there's a decision: have we proven infeasibility |
---|
1088 | (INFEASIBILITY) or monotonicity (MONOTONE), or do we need another |
---|
1089 | iteration (ITERATE). Monotonicity and iteration are large blocks |
---|
1090 | in themselves. |
---|
1091 | |
---|
1092 | There was a fair bit of preamble and postamble around the PASS loop. |
---|
1093 | The preamble code is protected by PROBING_EXTRA_STUFF and has been disabled |
---|
1094 | since 2006. It appears to find disaggregation cuts. The postamble code |
---|
1095 | purports to find bigM constraints. It's clearly a work in progress. I've |
---|
1096 | moved both out to tmp files. -- lh, 101203 -- |
---|
1097 | |
---|
1098 | Comments on the data structures: |
---|
1099 | |
---|
1100 | markC 0: not fixed |
---|
1101 | 1: tightened upper bound |
---|
1102 | 2: tightened lower bound |
---|
1103 | 3: fixed |
---|
1104 | |
---|
1105 | or perhaps the above should be interpreted as bit codes, along with the |
---|
1106 | bits used to indicate infinite bounds: |
---|
1107 | |
---|
1108 | 0x0: no bounds tightened |
---|
1109 | 0x1: u<j> tightened |
---|
1110 | 0x2: l<j> tightened |
---|
1111 | 0x4: l<j> < -1e10 (-infty) |
---|
1112 | 0x8: u<j> > 1e10 (+infty) |
---|
1113 | |
---|
1114 | Note that we'll only tighten a bound once --- this is enforced during |
---|
1115 | bounds propagation. |
---|
1116 | |
---|
1117 | stackC Records columns to be processed. This record is started anew |
---|
1118 | each time we probe, pushing a variable in a given direction. |
---|
1119 | The variable being probed is always entry 0 on stackC. |
---|
1120 | saveL and saveU are correlated with stackC. |
---|
1121 | |
---|
1122 | This isn't a stack --- it's a record of every variable |
---|
1123 | processed. nstackC is the tail of the queue, and istackC is |
---|
1124 | the head. At the end of a probing pass it's scanned to recover |
---|
1125 | the actual implications that were generated. It's not clear |
---|
1126 | to me why it's allocated at 2*ncols. Perhaps, at one time, |
---|
1127 | there were separate entries for upper and lower bound changes? |
---|
1128 | |
---|
1129 | No, just as I'm thinking I'm almost done, I notice there's some |
---|
1130 | fancy footwork in the portion of the code that handles the case |
---|
1131 | where both the up and down probes were feasible. Seems likely |
---|
1132 | that we're allowing for each direction to process all variables |
---|
1133 | at most once. |
---|
1134 | |
---|
1135 | And as I'm working through the code that handles the case |
---|
1136 | where down and up probes were both feasible, there in the |
---|
1137 | move singleton code is an assert that nstackC < nCols. |
---|
1138 | |
---|
1139 | And by the time I finally reach the end, I'm pretty well |
---|
1140 | convinced that it's an error if we exceed nCols on stackC. |
---|
1141 | The only thing to worry about is whether singleton processing |
---|
1142 | could add a variable that's already been added due to |
---|
1143 | propagation, and a check of the code says the answer is 'no'. |
---|
1144 | |
---|
1145 | TODO: When I rewrite this, allocate stackC as nCols and put in some |
---|
1146 | asserts to catch any violations. There are already asserts |
---|
1147 | in place at the point where a variable is placed on stackC |
---|
1148 | during propagation. (Look for `Is there any change to |
---|
1149 | propagate'.) -- lh, 101128 -- |
---|
1150 | |
---|
1151 | |
---|
1152 | stackR Stack rows to be processed. This stack is started anew each |
---|
1153 | time we probe pushing a variable in a given direction. |
---|
1154 | |
---|
1155 | minR, maxR, markR Initialised externally by tighten2. For row i, minR[i] = |
---|
1156 | LB(i), maxR[i] = UB(i) (row lhs lower and upper bounds, |
---|
1157 | respectively). markR is set to -1 if there's at least |
---|
1158 | one finite lhs bound, -2 if we shouldn't process this |
---|
1159 | row (no finite bounds, or too many coefficients). |
---|
1160 | |
---|
1161 | Then, as we work, markR[i] will be set to point to the entry in stackR for |
---|
1162 | row i. saveMin and saveMax are correlated with stackR, and hold the |
---|
1163 | original values for LB(i) and UB(i) (from minR and maxR) while we're |
---|
1164 | working. As it turns out, however, the the back pointer from markR is never |
---|
1165 | used. |
---|
1166 | |
---|
1167 | largestPositiveInRow for row i, max{j in P} a<ij>(u<i>-l<j>), where P is |
---|
1168 | the set of positive coefficients a<ij> |
---|
1169 | largestNegativeInRow for row i, min{j in M} a<ij>(u<i>-l<j>), where M is |
---|
1170 | the set of negative coefficients a<ij> |
---|
1171 | |
---|
1172 | element and index are scratch arrays used to construct cuts. |
---|
1173 | |
---|
1174 | realRows is the translation array created if we modified the constraint |
---|
1175 | system during preprocessing in gutsOfGenerateCuts. |
---|
1176 | |
---|
1177 | */ |
---|
1178 | int CglProbing::probe( const OsiSolverInterface &si, |
---|
1179 | OsiCuts &cs, |
---|
1180 | double *colLower, double *colUpper, |
---|
1181 | CoinPackedMatrix *rowCopy, |
---|
1182 | CoinPackedMatrix *columnCopy, |
---|
1183 | const CoinBigIndex *rowStartPos, const int *realRows, |
---|
1184 | const double *rowLower, const double *rowUpper, |
---|
1185 | const char *intVar, double *minR, double *maxR, |
---|
1186 | int *markR, |
---|
1187 | CglTreeInfo *info) const |
---|
1188 | |
---|
1189 | { |
---|
1190 | # if CGL_DEBUG > 0 |
---|
1191 | std::cout << "Entering CglProbing::probe." << std::endl ; |
---|
1192 | # endif |
---|
1193 | /* |
---|
1194 | PREPARATION |
---|
1195 | |
---|
1196 | All the code down through `PASS LOOP: HEAD' is preparation. Do all of the |
---|
1197 | setup that will not change over the nested loops that do the work of probing. |
---|
1198 | Note that rowCopy may have considerably fewer rows than the original system. |
---|
1199 | */ |
---|
1200 | |
---|
1201 | int nRows = rowCopy->getNumRows() ; |
---|
1202 | int nCols = rowCopy->getNumCols() ; |
---|
1203 | int maxStack = info->inTree ? maxStack_ : maxStackRoot_ ; |
---|
1204 | |
---|
1205 | /* |
---|
1206 | Fiddle maxStack. Odd sort of function: |
---|
1207 | 1 < maxStack <= 25 => 2*maxStack |
---|
1208 | 26 < maxStack <= 50 => 50 |
---|
1209 | 51 < maxStack => maxStack |
---|
1210 | |
---|
1211 | TODO: Grepping through the code says there's no way that totalTimesCalled_ |
---|
1212 | can become negative. Perhaps in some derived class? More likely a |
---|
1213 | magic number for something, or else a quick way to disable this |
---|
1214 | code. -- lh, 101021 -- |
---|
1215 | */ |
---|
1216 | if ((totalTimesCalled_%10) == -1) { |
---|
1217 | int newMax = CoinMin(2*maxStack,50) ; |
---|
1218 | maxStack = CoinMax(newMax,maxStack) ; |
---|
1219 | } |
---|
1220 | totalTimesCalled_++ ; |
---|
1221 | |
---|
1222 | # ifdef ONE_ARRAY |
---|
1223 | assert((DIratio != 0) && ((DIratio&(DIratio-1)) == 0)) ; |
---|
1224 | int nSpace = 8*nCols+4*nRows+2*maxStack ; |
---|
1225 | nSpace += (4*nCols+nRows+maxStack+DIratio-1)>>(DIratio-1) ; |
---|
1226 | double * colsol = new double[nSpace] ; |
---|
1227 | double * djs = colsol + nCols ; |
---|
1228 | double * columnGap = djs + nCols ; |
---|
1229 | double * saveL = columnGap + nCols ; |
---|
1230 | double * saveU = saveL + 2*nCols ; |
---|
1231 | double * saveMin = saveU + 2*nCols ; |
---|
1232 | double * saveMax = saveMin + nRows ; |
---|
1233 | double * largestPositiveInRow = saveMax + nRows ; |
---|
1234 | double * largestNegativeInRow = largestPositiveInRow + nRows ; |
---|
1235 | double * element = largestNegativeInRow + nRows ; |
---|
1236 | double * lo0 = element + nCols ; |
---|
1237 | double * up0 = lo0 + maxStack ; |
---|
1238 | int * markC = reinterpret_cast<int *> (up0+maxStack) ; |
---|
1239 | int * stackC = markC + nCols ; |
---|
1240 | int * stackR = stackC + 2*nCols ; |
---|
1241 | int * index = stackR + nRows ; |
---|
1242 | int * stackC0 = index + nCols ; |
---|
1243 | # else |
---|
1244 | double * colsol = new double[nCols] ; |
---|
1245 | double * djs = new double[nCols] ; |
---|
1246 | double * columnGap = new double [nCols] ; |
---|
1247 | double * saveL = new double [2*nCols] ; |
---|
1248 | double * saveU = new double [2*nCols] ; |
---|
1249 | double * saveMin = new double [nRows] ; |
---|
1250 | double * saveMax = new double [nRows] ; |
---|
1251 | double * largestPositiveInRow = new double [nRows] ; |
---|
1252 | double * largestNegativeInRow = new double [nRows] ; |
---|
1253 | double * element = new double[nCols] ; |
---|
1254 | double * lo0 = new double[maxStack] ; |
---|
1255 | double * up0 = new double[maxStack] ; |
---|
1256 | int * markC = new int [nCols] ; |
---|
1257 | int * stackC = new int [2*nCols] ; |
---|
1258 | int * stackR = new int [nRows] ; |
---|
1259 | int * index = new int[nCols] ; |
---|
1260 | int * stackC0 = new int[maxStack] ; |
---|
1261 | # endif |
---|
1262 | |
---|
1263 | /* |
---|
1264 | Create a local container to hold the cuts we generate. At the end we'll |
---|
1265 | transfer these to the container passed as a parameter. |
---|
1266 | |
---|
1267 | Typically, rowCopy (nRows) will be smaller than the original system in the |
---|
1268 | si, but it could be one larger (all original constraints plus the objective). |
---|
1269 | At the root, allow lots of room; even more if this is the first call at the |
---|
1270 | root. In the tree, much less. |
---|
1271 | |
---|
1272 | If all we're doing is strengthening rows in place (option 0x40 indicates |
---|
1273 | we're in preprocessing, a good place to do this), we can only work with |
---|
1274 | what we have in rowCopy (and we need a translation array). |
---|
1275 | |
---|
1276 | TODO: There's a problem with realRows here --- if we didn't delete any |
---|
1277 | rows during preprocessing in gutsOfGenerateCuts, realRows is not |
---|
1278 | allocated. So we can't strengthen in place in that case. Not likely |
---|
1279 | the intent. -- lh, 101209 -- |
---|
1280 | */ |
---|
1281 | bool justReplace = ((info->options&0x40) != 0) && (realRows != NULL) ; |
---|
1282 | int nRowsFake = 0 ; |
---|
1283 | if (justReplace) { |
---|
1284 | nRowsFake = nRows ; |
---|
1285 | } else { |
---|
1286 | int nRowsSafe = CoinMin(nRows,si.getNumRows()) ; |
---|
1287 | nRowsFake = info->inTree ? nRowsSafe/3 : nRowsSafe*10 ; |
---|
1288 | if (!info->inTree && info->pass == 0) nRowsFake *= 10 ; |
---|
1289 | } |
---|
1290 | CglProbingRowCut rowCut(nRowsFake,!info->inTree) ; |
---|
1291 | |
---|
1292 | |
---|
1293 | // Unpack matrices |
---|
1294 | const int * column = rowCopy->getIndices() ; |
---|
1295 | const CoinBigIndex * rowStart = rowCopy->getVectorStarts() ; |
---|
1296 | const double * rowElements = rowCopy->getElements() ; |
---|
1297 | |
---|
1298 | const int * row = columnCopy->getIndices() ; |
---|
1299 | const CoinBigIndex * columnStart = columnCopy->getVectorStarts() ; |
---|
1300 | const int * columnLength = columnCopy->getVectorLengths(); |
---|
1301 | const double * columnElements = columnCopy->getElements() ; |
---|
1302 | |
---|
1303 | /* |
---|
1304 | Grab the column solution and reduced costs and groom them. |
---|
1305 | */ |
---|
1306 | CoinMemcpyN(si.getReducedCost(),nCols,djs) ; |
---|
1307 | CoinMemcpyN(si.getColSolution(),nCols,colsol) ; |
---|
1308 | double direction = si.getObjSense() ; |
---|
1309 | groomSoln(direction,primalTolerance_,djs,colLower,colsol,colUpper,columnGap, |
---|
1310 | rowCopy,rowStartPos,largestNegativeInRow,largestPositiveInRow) ; |
---|
1311 | /* |
---|
1312 | Determine objective cutoff (minimisation convention). |
---|
1313 | |
---|
1314 | TODO: Seems to me that this bit of code is a guaranteed fail for a |
---|
1315 | maximisation problem with no cutoff. It also repeats in a number |
---|
1316 | of places. -- lh, 101125 -- |
---|
1317 | */ |
---|
1318 | double cutoff ; |
---|
1319 | bool cutoff_available = si.getDblParam(OsiDualObjectiveLimit,cutoff) ; |
---|
1320 | if (!cutoff_available || usingObjective_ < 0) { |
---|
1321 | cutoff = si.getInfinity() ; |
---|
1322 | } |
---|
1323 | cutoff *= direction ; |
---|
1324 | if (CoinAbs(cutoff) > 1.0e30) |
---|
1325 | assert (cutoff > 1.0e30) ; |
---|
1326 | double current = si.getObjValue() ; |
---|
1327 | current *= direction ; |
---|
1328 | /* |
---|
1329 | Scan the variables, noting the ones that are fixed or have a bound too large |
---|
1330 | to be useful. |
---|
1331 | */ |
---|
1332 | //int nFix = 0 ; |
---|
1333 | for (int i = 0 ; i < nCols ; i++) { |
---|
1334 | if (colUpper[i]-colLower[i] < 1.0e-8) { |
---|
1335 | markC[i] = tightenLower|tightenUpper ; |
---|
1336 | //nFix++ ; |
---|
1337 | } else { |
---|
1338 | markC[i] = 0 ; |
---|
1339 | if (colUpper[i] > 1.0e10) |
---|
1340 | markC[i] |= infUpper ; |
---|
1341 | if (colLower[i] < -1.0e10) |
---|
1342 | markC[i] |= infLower ; |
---|
1343 | } |
---|
1344 | } |
---|
1345 | //printf("PROBE %d fixed out of %d\n",nFix,nCols) ; |
---|
1346 | |
---|
1347 | /* |
---|
1348 | jjf: If we are going to replace coefficient then we don't need to be |
---|
1349 | effective |
---|
1350 | |
---|
1351 | Seems like this comment does not agree with CglProbingRowCut.addCuts(), |
---|
1352 | which checks effectiveness before adding a cut to the OsiCuts collection or |
---|
1353 | entering it into the strenghtenRow array. |
---|
1354 | |
---|
1355 | But it does agree with the way coefficient strengthening cuts are handled |
---|
1356 | down in the end-of-iteration processing for the down/up/down probe loop. |
---|
1357 | |
---|
1358 | It looks like the idea of strengthenRow is that a cut that's simply |
---|
1359 | strengthening an existing row i is entered in slot i of strengthenRow. It's |
---|
1360 | left to the client to deal with it on return. I need to get a better idea of |
---|
1361 | how justReplace is handled, here and in the client. -- lh, 101125 -- |
---|
1362 | */ |
---|
1363 | //double needEffectiveness = info->strengthenRow ? -1.0e10 : 1.0e-3 ; |
---|
1364 | double needEffectiveness = info->strengthenRow ? 1.0e-8 : 1.0e-3 ; |
---|
1365 | if (justReplace && ((info->pass&1) != 0)) |
---|
1366 | needEffectiveness = -1.0e10 ; |
---|
1367 | |
---|
1368 | double tolerance = 1.0e1*primalTolerance_ ; |
---|
1369 | |
---|
1370 | /* for both way coding */ |
---|
1371 | int nstackC0 = -1 ; |
---|
1372 | int nstackR,nstackC ; |
---|
1373 | /* |
---|
1374 | So, let's assume that the real business starts here, and the previous block |
---|
1375 | is irrelevant. |
---|
1376 | */ |
---|
1377 | /* |
---|
1378 | All the initialisation that occurs here is specific to CglTreeProbingInfo. |
---|
1379 | Sort of begs the question `Why is this handled as a derived class of |
---|
1380 | CglTreeInfo?' And why (what?) is CglTreeInfo, in the grand scheme of things? |
---|
1381 | |
---|
1382 | Some grepping about in Cbc says that we might see a CglTreeProbingInfo |
---|
1383 | object during processing of the root node. My tentative hypothesis is that |
---|
1384 | it's a vehicle to pass implications back to cbc, where they will be stashed |
---|
1385 | in a CglImplication object for further use. It's worth mentioning that the |
---|
1386 | existence of a CglTreeProbingInfo object depends on an independent symbol |
---|
1387 | (CLIQUE_ANALYSIS) being defined in CbcModel::branchAndBound. It's also worth |
---|
1388 | mentioning that the code here will core dump if the info object is not a |
---|
1389 | CglTreeProbingInfo object. |
---|
1390 | |
---|
1391 | And finally, it's worth mentioning that this code is disabled --- PROBING100 |
---|
1392 | is defined to 0 at the head of CglProbing --- since 080722. |
---|
1393 | |
---|
1394 | The trivial implementation CglTreeInfo::initializeFixing() always returns 0, |
---|
1395 | so in effect this block ensures saveFixingInfo is false. |
---|
1396 | |
---|
1397 | -- lh, 101126 -- |
---|
1398 | */ |
---|
1399 | bool saveFixingInfo = false ; |
---|
1400 | #if PROBING100 |
---|
1401 | CglTreeProbingInfo *probingInfo = dynamic_cast<CglTreeProbingInfo *> (info) ; |
---|
1402 | const int *backward = NULL ; |
---|
1403 | const int *integerVariable = NULL ; |
---|
1404 | const int *toZero = NULL ; |
---|
1405 | const int *toOne = NULL ; |
---|
1406 | const fixEntry *fixEntries = NULL ; |
---|
1407 | #endif |
---|
1408 | if (info->inTree) { |
---|
1409 | #if PROBING100 |
---|
1410 | backward = probingInfo->backward() ; |
---|
1411 | integerVariable = probingInfo->integerVariable() ; |
---|
1412 | toZero = probingInfo->toZero() ; |
---|
1413 | toOne = probingInfo->toOne() ; |
---|
1414 | fixEntries = probingInfo->fixEntries() ; |
---|
1415 | #endif |
---|
1416 | } else { |
---|
1417 | saveFixingInfo = (info->initializeFixing(&si) > 0) ; |
---|
1418 | } |
---|
1419 | /* |
---|
1420 | PASS LOOP: HEAD |
---|
1421 | |
---|
1422 | Major block #2: Main pass loop. |
---|
1423 | |
---|
1424 | From CglProbingAnn: anyColumnCuts is set only in the case that we've |
---|
1425 | fixed a variable by probing (i.e., one of the up or down probe resulted |
---|
1426 | in infeasibility) and that probe entailed column cuts. Once set, it is |
---|
1427 | never rescinded. In the reworked code, it's set as the return value of |
---|
1428 | monotoneActions(). |
---|
1429 | */ |
---|
1430 | bool anyColumnCuts = false ; |
---|
1431 | int ninfeas = 0 ; |
---|
1432 | int rowCuts = rowCuts_ ; |
---|
1433 | double disaggEffectiveness = 1.0e-3 ; |
---|
1434 | int ipass = 0 ; |
---|
1435 | int nfixed = -1 ; |
---|
1436 | int maxPass = info->inTree ? maxPass_ : maxPassRoot_ ; |
---|
1437 | |
---|
1438 | while (ipass < maxPass && nfixed) { |
---|
1439 | int iLook ; |
---|
1440 | ipass++ ; |
---|
1441 | //printf("pass %d\n",ipass) ; |
---|
1442 | nfixed = 0 ; |
---|
1443 | /* |
---|
1444 | We went through a fair bit of trouble in gutsOfGenerateCut to determine |
---|
1445 | the set of variables to be probed and loaded the indices into lookedAt_; |
---|
1446 | numberThisTime_ reflects that. |
---|
1447 | |
---|
1448 | The net effect here is that the root gets special treatment |
---|
1449 | (maxProbeRoot_) and the first pass at the root gets extra special treatment |
---|
1450 | (numberThisTime_). |
---|
1451 | |
---|
1452 | See CbcCutGenerator::generateCuts for the origin of magic number 123. |
---|
1453 | Comments there indicate it's intended to take effect deep in the tree, but |
---|
1454 | the code here will also work at the root. |
---|
1455 | |
---|
1456 | Special cases aside, the net effect of the loops is to walk lookedAt_, |
---|
1457 | promoting every nth variable (n = cutDown) to the front of lookedAt_. When |
---|
1458 | the loop finishes, the rest of the variables (which were copied off to |
---|
1459 | stackC) are appended to lookedAt_. If XXXXXX is not defined, we have an |
---|
1460 | expensive noop. -- lh, 101126 -- |
---|
1461 | */ |
---|
1462 | int justFix = (!info->inTree && !info->pass) ? -1 : 0 ; |
---|
1463 | int maxProbe = info->inTree ? maxProbe_ : maxProbeRoot_ ; |
---|
1464 | if (justFix < 0) |
---|
1465 | maxProbe = numberThisTime_ ; |
---|
1466 | if (maxProbe == 123) { |
---|
1467 | maxProbe = 0 ; |
---|
1468 | if (!info->inTree) { |
---|
1469 | if (!info->pass || numberThisTime_ < -100) { |
---|
1470 | maxProbe = numberThisTime_ ; |
---|
1471 | } else { |
---|
1472 | int cutDown = 4 ; |
---|
1473 | int offset = info->pass%cutDown ; |
---|
1474 | int i ; |
---|
1475 | int k = 0 ; |
---|
1476 | int kk = offset ; |
---|
1477 | for (i = 0 ; i < numberThisTime_ ; i++) { |
---|
1478 | if (!kk) { |
---|
1479 | #define XXXXXX |
---|
1480 | #ifdef XXXXXX |
---|
1481 | lookedAt_[maxProbe] = lookedAt_[i] ; |
---|
1482 | #endif |
---|
1483 | maxProbe++ ; |
---|
1484 | kk = cutDown-1 ; |
---|
1485 | } else { |
---|
1486 | stackC[k++] = lookedAt_[i] ; |
---|
1487 | kk-- ; |
---|
1488 | } |
---|
1489 | } |
---|
1490 | #ifdef XXXXXX |
---|
1491 | memcpy(lookedAt_+maxProbe,stackC,k*sizeof(int)) ; |
---|
1492 | #endif |
---|
1493 | } |
---|
1494 | } else { |
---|
1495 | if (numberThisTime_ < 200) { |
---|
1496 | maxProbe = numberThisTime_ ; |
---|
1497 | } else { |
---|
1498 | int cutDown = CoinMax(numberThisTime_/100,4) ; |
---|
1499 | int offset = info->pass%cutDown ; |
---|
1500 | int i ; |
---|
1501 | int k = 0 ; |
---|
1502 | int kk = offset ; |
---|
1503 | for (i = 0 ; i < numberThisTime_ ; i++) { |
---|
1504 | if (!kk) { |
---|
1505 | #ifdef XXXXXX |
---|
1506 | lookedAt_[maxProbe] = lookedAt_[i] ; |
---|
1507 | #endif |
---|
1508 | maxProbe++ ; |
---|
1509 | kk = cutDown-1 ; |
---|
1510 | } else { |
---|
1511 | stackC[k++] = lookedAt_[i] ; |
---|
1512 | kk-- ; |
---|
1513 | } |
---|
1514 | } |
---|
1515 | #ifdef XXXXXX |
---|
1516 | memcpy(lookedAt_+maxProbe,stackC,k*sizeof(int)) ; |
---|
1517 | #endif |
---|
1518 | } |
---|
1519 | } |
---|
1520 | } // end maxProbe == 123 |
---|
1521 | |
---|
1522 | /* |
---|
1523 | This looks to be an overall limit on probing. It's decremented every time a |
---|
1524 | variable is popped off stackC for processing. |
---|
1525 | |
---|
1526 | TODO: PROBING5 would disable this limit; currently inactive. -- lh, 101127 -- |
---|
1527 | */ |
---|
1528 | int leftTotalStack = maxStack*CoinMax(200,maxProbe) ; |
---|
1529 | #ifdef PROBING5 |
---|
1530 | if (!info->inTree&&!info->pass) |
---|
1531 | leftTotalStack = 1234567890 ; |
---|
1532 | #endif |
---|
1533 | //printf("maxStack %d maxPass %d numberThisTime %d info pass %d\n", |
---|
1534 | // maxStack,maxPass,numberThisTime_,info->pass) ; |
---|
1535 | /* |
---|
1536 | LOOKEDAT LOOP: HEAD |
---|
1537 | |
---|
1538 | Main loop to probe each variable in lookedAt_ |
---|
1539 | |
---|
1540 | We're finally ready to get down to the business of probing. Open a loop to |
---|
1541 | probe each variable in lookedAt_. |
---|
1542 | */ |
---|
1543 | for (iLook = 0 ; iLook < numberThisTime_ ; iLook++) { |
---|
1544 | double solval ; |
---|
1545 | double down ; |
---|
1546 | double up ; |
---|
1547 | /* |
---|
1548 | Too successful? Consider bailing out. |
---|
1549 | |
---|
1550 | If we're generating row cuts, but haven't fixed any variables or we're in |
---|
1551 | the tree, break from probing. |
---|
1552 | |
---|
1553 | JustFix < 0 says this is the first pass at the root; for any other |
---|
1554 | iteration of the pass loop, it'll be initialised to 0. If it is the first |
---|
1555 | pass at the root, turn off row cuts, keep on fixing variables, and stop |
---|
1556 | at the end of this pass. Otherwise, if we haven't fixed any variables, break. |
---|
1557 | |
---|
1558 | TODO: Otherwise keep going with no changes? That doesn't seem right. The |
---|
1559 | logic here does not cover all situations. -- lh, 101126 -- |
---|
1560 | */ |
---|
1561 | if (rowCut.outOfSpace() || leftTotalStack <= 0) { |
---|
1562 | if (!justFix && (!nfixed || info->inTree)) { |
---|
1563 | #ifdef COIN_DEVELOP |
---|
1564 | if (!info->inTree) |
---|
1565 | printf("Exiting a on pass %d, maxProbe %d\n", |
---|
1566 | ipass,maxProbe) ; |
---|
1567 | #endif |
---|
1568 | break ; |
---|
1569 | } else if (justFix <= 0) { |
---|
1570 | if (!info->inTree) { |
---|
1571 | rowCuts = 0 ; |
---|
1572 | justFix = 1 ; |
---|
1573 | disaggEffectiveness = COIN_DBL_MAX ; |
---|
1574 | needEffectiveness = COIN_DBL_MAX ; |
---|
1575 | //maxStack=10 ; |
---|
1576 | maxPass = 1 ; |
---|
1577 | } else if (!nfixed) { |
---|
1578 | #ifdef COIN_DEVELOP |
---|
1579 | printf("Exiting b on pass %d, maxProbe %d\n", |
---|
1580 | ipass,maxProbe) ; |
---|
1581 | #endif |
---|
1582 | break ; |
---|
1583 | } |
---|
1584 | } |
---|
1585 | } |
---|
1586 | /* |
---|
1587 | awayFromBound isn't quite accurate --- we'll also set awayFromBound = 0 |
---|
1588 | for a general integer variable at an integer value strictly within bounds. |
---|
1589 | */ |
---|
1590 | int awayFromBound = 1 ; |
---|
1591 | /* |
---|
1592 | Have a look at the current variable. We're not interested in processing it |
---|
1593 | if either or both bounds have been improved (0x02|0x01). If the variable has |
---|
1594 | become fixed, claim both bounds have improved. |
---|
1595 | |
---|
1596 | TODO: if x<j> is not an integer variable we have big problems. This should |
---|
1597 | be an assert. -- lh, 101126 -- |
---|
1598 | */ |
---|
1599 | int j = lookedAt_[iLook] ; |
---|
1600 | solval = colsol[j] ; |
---|
1601 | down = floor(solval+tolerance) ; |
---|
1602 | up = ceil(solval-tolerance) ; |
---|
1603 | if (columnGap[j] < 0.0) markC[j] = tightenUpper|tightenLower ; |
---|
1604 | if ((markC[j]&(tightenUpper|tightenLower)) != 0 || !intVar[j]) continue ; |
---|
1605 | /* |
---|
1606 | `Normalize' variables that are near (or past) their bounds. |
---|
1607 | |
---|
1608 | In normal circumstances, u<j> - l<j> is at least 1, while tol will be |
---|
1609 | down around 10e-5 (given a primal tolerance of 10e-6, about as large as |
---|
1610 | you'd normally see). Then |
---|
1611 | |
---|
1612 | l<j> < l<j>+tol < l<j>+2tol << u<j>-2tol < u<j>-tol < u<j> |
---|
1613 | |
---|
1614 | For x<j> > u<j>-tol, (x<j>,u<j>,u<j>) => (u<j>-.1, l<j>, u<j>) |
---|
1615 | For x<j> < l<j>+tol, (x<j>,l<j>,l<j>) => (l<j>+.1, l<j>, u<j>) |
---|
1616 | For up == down, (x<j>,v,v) => (v+.1,v,v+1) |
---|
1617 | |
---|
1618 | So we're forcing a spread of 1 between up and down, and moving x<j> |
---|
1619 | to something far enough (.1) from integer to avoid accuracy problems when |
---|
1620 | calculating movement. A side effect is that primal infeasibility is |
---|
1621 | corrected. |
---|
1622 | |
---|
1623 | TODO: Why this structure? Another approach would be to test using up and |
---|
1624 | down calculated just above. It works out pretty much the same if |
---|
1625 | you stick with variables of type double. Moving to int would likely |
---|
1626 | make it faster. -- lh, 101127 -- |
---|
1627 | |
---|
1628 | TODO: Things will get seriously wierd if (u<j> - l<j>) <= 3tol. For |
---|
1629 | u<j>-tol < x<j> < l<j>+2tol, values will change. But the symmetric |
---|
1630 | case u<j>-2tol < x<j> < l<j>+tol will never change values because |
---|
1631 | of the order of the if statements. Clearly the code is not designed |
---|
1632 | to cope with this level of pathology. -- lh, 101127 -- |
---|
1633 | |
---|
1634 | TODO: It's worth noting that awayFromBound is never used. -- lh, 101127 -- |
---|
1635 | |
---|
1636 | TODO: It's worth noting that the values set for solval are never used except |
---|
1637 | in a debug print. -- lh, 101213 -- |
---|
1638 | |
---|
1639 | TODO: The final case below corresponds to l<j>+1 <= down == up < u<j>-1, |
---|
1640 | e.g., an interior integer value. Seems like we'd want to probe +1 and |
---|
1641 | -1. Presumably the code isn't set up to do this, so we arbitrarily |
---|
1642 | pick the up probe. Could the code be augmented to do the stronger |
---|
1643 | probe? -- lh, 101127 -- |
---|
1644 | */ |
---|
1645 | if (solval >= colUpper[j]-tolerance || |
---|
1646 | solval <= colLower[j]+tolerance || up == down) { |
---|
1647 | awayFromBound = 0 ; |
---|
1648 | if (solval <= colLower[j]+2.0*tolerance) { |
---|
1649 | solval = colLower[j]+1.0e-1 ; |
---|
1650 | down = colLower[j] ; |
---|
1651 | up = down+1.0 ; |
---|
1652 | } else if (solval >= colUpper[j]-2.0*tolerance) { |
---|
1653 | solval = colUpper[j]-1.0e-1 ; |
---|
1654 | up = colUpper[j] ; |
---|
1655 | down = up-1 ; |
---|
1656 | } else { |
---|
1657 | up = down+1.0 ; |
---|
1658 | solval = down+1.0e-1 ; |
---|
1659 | } |
---|
1660 | } |
---|
1661 | assert (up <= colUpper[j]) ; |
---|
1662 | assert (down >= colLower[j]) ; |
---|
1663 | assert (up > down) ; |
---|
1664 | # if CGL_DEBUG > 0 |
---|
1665 | const double lj = colLower[j] ; |
---|
1666 | const double uj = colUpper[j] ; |
---|
1667 | const bool downIsLower = (CoinAbs(down-colLower[j]) < 1.0e-7) ; |
---|
1668 | const bool upIsUpper = (CoinAbs(up-colUpper[j]) < 1.0e-7) ; |
---|
1669 | # endif |
---|
1670 | /* |
---|
1671 | Set up to probe each variable down (1), up (2), down (1). |
---|
1672 | |
---|
1673 | The notion is that we'll set a bit (1, 2, 4) in the result record if we're |
---|
1674 | feasible for a particular way. Way is defined as: |
---|
1675 | 1: we're looking at movement between floor(x*<j>) and x*<j>, u<j> |
---|
1676 | 2: we're looking at movement between ceil(x*<j>) and x*<j>, l<j> |
---|
1677 | As defined here, movement will be negative for way = 1. |
---|
1678 | |
---|
1679 | Why do we have a second pass with way = 1? Glad you asked. About 1200 |
---|
1680 | lines farther down, it finally becomes apparent that we execute the |
---|
1681 | final pass only when the initial down probe was feasible, but the up probe |
---|
1682 | was not. In a nutshell, we're restoring the state produced by the down |
---|
1683 | probe, which we'll then nail down in the section of code labelled `keep', |
---|
1684 | a mere 600 lines farther on. |
---|
1685 | */ |
---|
1686 | unsigned int iway ; |
---|
1687 | unsigned int way[] = { probeDown, probeUp, probeDown } ; |
---|
1688 | unsigned int feasValue[] = |
---|
1689 | { downProbeFeas, upProbeFeas, downProbe2Feas } ; |
---|
1690 | unsigned int feasRecord = 0 ; |
---|
1691 | bool notFeasible = false ; |
---|
1692 | int istackC = 0 ; |
---|
1693 | int istackR = 0 ; |
---|
1694 | /* |
---|
1695 | PROBE LOOP: HEAD |
---|
1696 | |
---|
1697 | Open a loop to probe up and down for the current variable. Get things |
---|
1698 | started by placing the variable on the propagation queue (stackC). |
---|
1699 | |
---|
1700 | As with the previous loops (variables in lookedAt_, passes) this loop |
---|
1701 | extends to the end of major block #2). |
---|
1702 | */ |
---|
1703 | for (iway = downProbe ; iway < oneProbeTooMany ; iway++) { |
---|
1704 | stackC[0] = j ; |
---|
1705 | markC[j] = way[iway] ; |
---|
1706 | /* |
---|
1707 | Calculate movement given the current direction. Given the setup above, |
---|
1708 | movement should be at least +/-1. |
---|
1709 | |
---|
1710 | TODO: Notice we reset up or down yet again. Really, this shouldn't be |
---|
1711 | necessary. For that matter, why did we bother with awayFromBound |
---|
1712 | immediately above? -- lh, 101127 -- |
---|
1713 | |
---|
1714 | TODO: The more I think about this, the less happy I am. Seems to me that |
---|
1715 | neither colUpper or colLower can change during iterations of this |
---|
1716 | loop. If we discover monotone, we save bounds and break. If we |
---|
1717 | don't discover monotone, we restore. Let's test that. |
---|
1718 | -- lh, 110105 -- |
---|
1719 | */ |
---|
1720 | double solMovement ; |
---|
1721 | double movement ; |
---|
1722 | int goingToTrueBound = 0 ; |
---|
1723 | |
---|
1724 | # if CGL_DEBUG > 0 |
---|
1725 | assert(lj == colLower[j]) ; |
---|
1726 | assert(uj == colUpper[j]) ; |
---|
1727 | assert(downIsLower == (CoinAbs(down-colLower[j]) < 1.0e-7)) ; |
---|
1728 | assert(upIsUpper == (CoinAbs(up-colUpper[j]) < 1.0e-7)) ; |
---|
1729 | # endif |
---|
1730 | |
---|
1731 | if (way[iway] == probeDown) { |
---|
1732 | movement = down-colUpper[j] ; |
---|
1733 | solMovement = down-colsol[j] ; |
---|
1734 | assert(movement < -0.99999) ; |
---|
1735 | if (CoinAbs(down-colLower[j]) < 1.0e-7) { |
---|
1736 | goingToTrueBound = 2 ; |
---|
1737 | down = colLower[j] ; |
---|
1738 | } |
---|
1739 | } else { |
---|
1740 | movement = up-colLower[j] ; |
---|
1741 | solMovement = up-colsol[j] ; |
---|
1742 | assert(movement > 0.99999) ; |
---|
1743 | if (CoinAbs(up-colUpper[j]) < 1.0e-7) { |
---|
1744 | goingToTrueBound = 2 ; |
---|
1745 | up = colUpper[j] ; |
---|
1746 | } |
---|
1747 | } |
---|
1748 | /* |
---|
1749 | Coding for goingToTrueBound is: |
---|
1750 | 0: We're somewhere in the interior of a general integer. |
---|
1751 | 1: We're heading for one of the bounds on a general integer. |
---|
1752 | 2: We're processing a binary variable (u<j>-l<j> = 1 and l<j> = 0). |
---|
1753 | You can view this next test as correcting for the fact that the previous |
---|
1754 | code assumes binary variables are all there is. |
---|
1755 | |
---|
1756 | TODO: Now that I've figured out the math for the constraints generated |
---|
1757 | by disaggCuts (see typeset documentation), I see what's wrong |
---|
1758 | (?) here. The disagg cut that's generated relies on the new probe |
---|
1759 | bound being distance one from the original bound. But (for example) |
---|
1760 | that's (colUpper-down) = 1 for a down probe. Not quite what's |
---|
1761 | being checked for here. But this next test ensures binary variables, |
---|
1762 | and then the tests are equivalent. Before I willy-nilly change this, I |
---|
1763 | need to sort out a couple of other places where goingToTrueBound |
---|
1764 | controls behaviour. -- lh, 101214 -- |
---|
1765 | |
---|
1766 | TODO: Apropos the above, for example, the test for coefficient strengthening |
---|
1767 | cuts includes goingToTrueBound != 0 -- lh, 110105 -- |
---|
1768 | */ |
---|
1769 | if (goingToTrueBound && (colUpper[j]-colLower[j] > 1.5 || colLower[j])) |
---|
1770 | goingToTrueBound = 1 ; |
---|
1771 | /* |
---|
1772 | If the user doesn't want disaggregation cuts, pretend we're in the interior |
---|
1773 | of a general integer variable. |
---|
1774 | |
---|
1775 | TODO: Why is row strengthening limited to binary variables? If we can |
---|
1776 | figure out what to do, the limitation seems artificial. |
---|
1777 | -- lh, 101127 -- |
---|
1778 | The test here also says that we can't have coefficient strengthening |
---|
1779 | without disaggregation. I don't see any reason for this dominance. |
---|
1780 | -- lh, 110105 -- |
---|
1781 | */ |
---|
1782 | if ((rowCuts&1) == 0) |
---|
1783 | goingToTrueBound = 0 ; |
---|
1784 | bool canReplace = info->strengthenRow && (goingToTrueBound == 2) ; |
---|
1785 | |
---|
1786 | #ifdef PRINT_DEBUG |
---|
1787 | // Alert for general integer with nontrivial range. |
---|
1788 | // Also last use of solval |
---|
1789 | if (CoinAbs(movement)>1.01) { |
---|
1790 | printf("big %d %g %g %g\n",j,colLower[j],solval,colUpper[j]) ; |
---|
1791 | } |
---|
1792 | #endif |
---|
1793 | /* |
---|
1794 | Recall that we adjusted the reduced costs (djs) and current objective |
---|
1795 | (current) to the minimisation sign convention. objVal will accumulate |
---|
1796 | objective change due to forced bound changes. |
---|
1797 | */ |
---|
1798 | double objVal = current ; |
---|
1799 | if (solMovement*djs[j] > 0.0) |
---|
1800 | objVal += solMovement*djs[j] ; |
---|
1801 | nstackC = 1 ; |
---|
1802 | nstackR = 0 ; |
---|
1803 | saveL[0] = colLower[j] ; |
---|
1804 | saveU[0] = colUpper[j] ; |
---|
1805 | assert (saveU[0] > saveL[0]) ; |
---|
1806 | notFeasible = false ; |
---|
1807 | /* |
---|
1808 | Recalculate the upper (down probe) or lower (up probe) bound. You can see |
---|
1809 | from the debug print that arbitrary integers are an afterthought in this |
---|
1810 | code. |
---|
1811 | |
---|
1812 | TODO: And why not just set the bounds to down (colUpper) or up (colLower)? |
---|
1813 | We set movement = down-colUpper just above, and we've taken a lot of |
---|
1814 | care to groom up and down. This is pointless effort. -- lh, 101214 -- |
---|
1815 | */ |
---|
1816 | if (movement < 0.0) { |
---|
1817 | colUpper[j] += movement ; |
---|
1818 | colUpper[j] = floor(colUpper[j]+0.5) ; |
---|
1819 | columnGap[j] = colUpper[j]-colLower[j]-primalTolerance_ ; |
---|
1820 | #ifdef PRINT_DEBUG |
---|
1821 | printf("** Trying %d down to 0\n",j) ; |
---|
1822 | #endif |
---|
1823 | } else { |
---|
1824 | colLower[j] += movement ; |
---|
1825 | colLower[j] = floor(colLower[j]+0.5) ; |
---|
1826 | columnGap[j] = colUpper[j]-colLower[j]-primalTolerance_ ; |
---|
1827 | #ifdef PRINT_DEBUG |
---|
1828 | printf("** Trying %d up to 1\n",j) ; |
---|
1829 | #endif |
---|
1830 | } |
---|
1831 | /* |
---|
1832 | Is this variable now fixed? |
---|
1833 | */ |
---|
1834 | if (CoinAbs(colUpper[j]-colLower[j]) < 1.0e-6) |
---|
1835 | markC[j] = tightenUpper|tightenLower ; |
---|
1836 | /* |
---|
1837 | Clear the infinite bound bits (infUpper|infLower) and check again. Bounds |
---|
1838 | may have improved. |
---|
1839 | */ |
---|
1840 | markC[j] &= ~(infUpper|infLower) ; |
---|
1841 | if (colUpper[j] > 1.0e10) |
---|
1842 | markC[j] |= infUpper ; |
---|
1843 | if (colLower[j] < -1.0e10) |
---|
1844 | markC[j] |= infLower ; |
---|
1845 | istackC = 0 ; |
---|
1846 | /* |
---|
1847 | Update row bounds to reflect the change in variable bound. |
---|
1848 | |
---|
1849 | TODO: Replacing code to update row bounds with updateRowBounds(). We'll |
---|
1850 | see how it looks. -- lh, 101203 -- |
---|
1851 | */ |
---|
1852 | |
---|
1853 | if (!updateRowBounds(j,movement,columnStart,columnLength,row,columnElements, |
---|
1854 | rowUpper,rowLower,nstackR,stackR,markR, |
---|
1855 | minR,maxR,saveMin,saveMax)) { |
---|
1856 | notFeasible = true ; |
---|
1857 | istackC = 1 ; |
---|
1858 | } |
---|
1859 | /* |
---|
1860 | PROBE LOOP: BEGIN PROPAGATION |
---|
1861 | |
---|
1862 | Row bounds are now adjusted for all rows with a<ij> != 0. Time to consider |
---|
1863 | the effects. nstackC is incremented each time we add a variable, istackC is |
---|
1864 | incremented each time we finish processing it. |
---|
1865 | |
---|
1866 | If we lost feasibility above, istackC = nstackC = 1 and this loop will not |
---|
1867 | execute. |
---|
1868 | |
---|
1869 | TODO: Is stackC really a stack? Or a queue? And what is the role of |
---|
1870 | maxStack? stackC is allocated at 2*ncols, but maxStack defaults to |
---|
1871 | 50. -- lh, 101127 -- |
---|
1872 | |
---|
1873 | TODO: stackC is clearly a queue. Allocation strategy is still unclear. |
---|
1874 | -- lh, 101128 -- |
---|
1875 | |
---|
1876 | jjf: could be istackC<maxStack? |
---|
1877 | */ |
---|
1878 | while (istackC < nstackC && nstackC < maxStack) { |
---|
1879 | leftTotalStack-- ; |
---|
1880 | int jway ; |
---|
1881 | int jcol = stackC[istackC] ; |
---|
1882 | jway = markC[jcol] ; |
---|
1883 | /* |
---|
1884 | jjf: If not first and fixed then skip |
---|
1885 | |
---|
1886 | It's clear that x<j>, the probing variable, can be marked as fixed at this |
---|
1887 | point (happens just above, when the upper or lower bound is adjusted). And |
---|
1888 | there's an earlier check that avoids probing a variable that's fixed when |
---|
1889 | plucked from lookedAt_. But ... why leave the if with an empty body? |
---|
1890 | |
---|
1891 | TODO: Disabled since r118 050225, but seems like the right thing to do. Is it |
---|
1892 | somehow guaranteed that fixed variables are not stacked? |
---|
1893 | -- lh, 101127 -- |
---|
1894 | |
---|
1895 | TODO: Based on what's happening down below, it looks like we can encounter |
---|
1896 | variables here which have had both bounds tightened and are then |
---|
1897 | pushed onto istackC to propagate the effect of that. This bit of code |
---|
1898 | should simply go away. -- lh, 101127 -- |
---|
1899 | */ |
---|
1900 | if (((jway&(tightenLower|tightenUpper)) == |
---|
1901 | (tightenLower|tightenUpper)) && |
---|
1902 | istackC) { |
---|
1903 | //istackC++ ; |
---|
1904 | //continue ; |
---|
1905 | //printf("fixed %d on stack\n",jcol) ; |
---|
1906 | } |
---|
1907 | #if PROBING100 |
---|
1908 | /* |
---|
1909 | This block of code has been disabled since r661 080722. I have no notes on |
---|
1910 | it from my old annotated CglProbing. Given the use of backward, it's tied to |
---|
1911 | CglTreeProbingInfo. -- lh, 101127 -- |
---|
1912 | */ |
---|
1913 | if (backward) { |
---|
1914 | int jColumn = backward[jcol] ; |
---|
1915 | if (jColumn>=0) { |
---|
1916 | int nFix=0 ; |
---|
1917 | // 0-1 see what else could be fixed |
---|
1918 | if (jway==tightenUpper) { |
---|
1919 | // fixed to 0 |
---|
1920 | int j ; |
---|
1921 | for (j = toZero_[jColumn];j<toOne_[jColumn];j++) { |
---|
1922 | int kColumn=fixEntry_[j].sequence ; |
---|
1923 | kColumn = integerVariable_[kColumn] ; |
---|
1924 | bool fixToOne = fixEntry_[j].oneFixed ; |
---|
1925 | if (fixToOne) { |
---|
1926 | if (colLower[kColumn]==0.0) { |
---|
1927 | if (colUpper[kColumn]==1.0) { |
---|
1928 | // See if on list |
---|
1929 | if (!(markC[kColumn]&(tightenLower|tightenUpper))) { |
---|
1930 | if(nStackC<nCols) { |
---|
1931 | stackC[nstackC]=kColumn ; |
---|
1932 | saveL[nstackC]=colLower[kColumn] ; |
---|
1933 | saveU[nstackC]=colUpper[kColumn] ; |
---|
1934 | assert (saveU[nstackC]>saveL[nstackC]) ; |
---|
1935 | assert (nstackC<nCols) ; |
---|
1936 | nstackC++ ; |
---|
1937 | markC[kColumn]|=tightenLower ; |
---|
1938 | nFix++ ; |
---|
1939 | } |
---|
1940 | } else if ((markC[kColumn]&(tightenLower|tightenUpper))==tightenUpper) { |
---|
1941 | notFeasible = true ; |
---|
1942 | } |
---|
1943 | } else { |
---|
1944 | // infeasible! |
---|
1945 | notFeasible = true ; |
---|
1946 | } |
---|
1947 | } |
---|
1948 | } else { |
---|
1949 | if (colUpper[kColumn]==1.0) { |
---|
1950 | if (colLower[kColumn]==0.0) { |
---|
1951 | // See if on list |
---|
1952 | if (!(markC[kColumn]&(tightenLower|tightenUpper))) { |
---|
1953 | if(nStackC<nCols) { |
---|
1954 | stackC[nstackC]=kColumn ; |
---|
1955 | saveL[nstackC]=colLower[kColumn] ; |
---|
1956 | saveU[nstackC]=colUpper[kColumn] ; |
---|
1957 | assert (saveU[nstackC]>saveL[nstackC]) ; |
---|
1958 | assert (nstackC<nCols) ; |
---|
1959 | nstackC++ ; |
---|
1960 | markC[kColumn]|=tightenUpper ; |
---|
1961 | nFix++ ; |
---|
1962 | } |
---|
1963 | } else if ((markC[kColumn]&(tightenLower|tightenUpper))==tightenLower) { |
---|
1964 | notFeasible = true ; |
---|
1965 | } |
---|
1966 | } else { |
---|
1967 | // infeasible! |
---|
1968 | notFeasible = true ; |
---|
1969 | } |
---|
1970 | } |
---|
1971 | } |
---|
1972 | } |
---|
1973 | } else if (jway==tightenLower) { |
---|
1974 | int j ; |
---|
1975 | for ( j=toOne_[jColumn];j<toZero_[jColumn+1];j++) { |
---|
1976 | int kColumn=fixEntry_[j].sequence ; |
---|
1977 | kColumn = integerVariable_[kColumn] ; |
---|
1978 | bool fixToOne = fixEntry_[j].oneFixed ; |
---|
1979 | if (fixToOne) { |
---|
1980 | if (colLower[kColumn]==0.0) { |
---|
1981 | if (colUpper[kColumn]==1.0) { |
---|
1982 | // See if on list |
---|
1983 | if (!(markC[kColumn]&(tightenLower|tightenUpper))) { |
---|
1984 | if(nStackC<nCols) { |
---|
1985 | stackC[nstackC]=kColumn ; |
---|
1986 | saveL[nstackC]=colLower[kColumn] ; |
---|
1987 | saveU[nstackC]=colUpper[kColumn] ; |
---|
1988 | assert (saveU[nstackC]>saveL[nstackC]) ; |
---|
1989 | assert (nstackC<nCols) ; |
---|
1990 | nstackC++ ; |
---|
1991 | markC[kColumn]|=tightenLower ; |
---|
1992 | nFix++ ; |
---|
1993 | } |
---|
1994 | } else if ((markC[kColumn]&(tightenLower|tightenUpper))==tightenUpper) { |
---|
1995 | notFeasible = true ; |
---|
1996 | } |
---|
1997 | } else { |
---|
1998 | // infeasible! |
---|
1999 | notFeasible = true ; |
---|
2000 | } |
---|
2001 | } |
---|
2002 | } else { |
---|
2003 | if (colUpper[kColumn]==1.0) { |
---|
2004 | if (colLower[kColumn]==0.0) { |
---|
2005 | // See if on list |
---|
2006 | if (!(markC[kColumn]&(tightenLower|tightenUpper))) { |
---|
2007 | if(nStackC<nCols) { |
---|
2008 | stackC[nstackC]=kColumn ; |
---|
2009 | saveL[nstackC]=colLower[kColumn] ; |
---|
2010 | saveU[nstackC]=colUpper[kColumn] ; |
---|
2011 | assert (saveU[nstackC]>saveL[nstackC]) ; |
---|
2012 | assert (nstackC<nCols) ; |
---|
2013 | nstackC++ ; |
---|
2014 | markC[kColumn]|=tightenUpper ; |
---|
2015 | nFix++ ; |
---|
2016 | } |
---|
2017 | } else if ((markC[kColumn]&(tightenLower|tightenUpper))==tightenLower) { |
---|
2018 | notFeasible = true ; |
---|
2019 | } |
---|
2020 | } else { |
---|
2021 | // infeasible! |
---|
2022 | notFeasible = true ; |
---|
2023 | } |
---|
2024 | } |
---|
2025 | } |
---|
2026 | } |
---|
2027 | } |
---|
2028 | } |
---|
2029 | } |
---|
2030 | #endif // PROBING100 (disabled) |
---|
2031 | /* |
---|
2032 | PROBE LOOP: WALK COLUMN |
---|
2033 | |
---|
2034 | Loop to walk the column of a variable popped off stackC. |
---|
2035 | |
---|
2036 | We've pulled x<jcol> off stackC. We're going to walk the column and process |
---|
2037 | each row where a<i,jcol> != 0. Note that we've already updated the row |
---|
2038 | bounds to reflect the change in bound for x<jcol>. In effect, we've stacked |
---|
2039 | this variable as a surrogate for stacking the rows whose bounds were |
---|
2040 | changed by the change to bounds on this variable. Now we're going to examine |
---|
2041 | those rows and see if any look promising for further propagation. |
---|
2042 | |
---|
2043 | */ |
---|
2044 | for (int k = columnStart[jcol] ; |
---|
2045 | k < columnStart[jcol]+columnLength[jcol] ; k++) { |
---|
2046 | if (notFeasible) |
---|
2047 | break ; |
---|
2048 | int irow = row[k] ; |
---|
2049 | // Row already processed with these bounds. |
---|
2050 | if (markR[irow]/10000 > 0) continue ; |
---|
2051 | |
---|
2052 | int rStart = rowStart[irow] ; |
---|
2053 | int rEnd = rowStartPos[irow] ; |
---|
2054 | double rowUp = rowUpper[irow] ; |
---|
2055 | double rowUp2 = 0.0 ; |
---|
2056 | bool doRowUpN ; |
---|
2057 | bool doRowUpP ; |
---|
2058 | /* |
---|
2059 | So, is this row promising? |
---|
2060 | |
---|
2061 | If our current L<i> (minR) is larger than the original b<i> (rowUp), we're |
---|
2062 | infeasible. |
---|
2063 | |
---|
2064 | Otherwise, a bit of linear algebra brings the conclusion that if the |
---|
2065 | largest negative gap (min{j in N} a<ij>(u<j>-l<j>)) in the row is less than |
---|
2066 | -(b<i>-L<i>), we aren't going to be able to make any progress shrinking the |
---|
2067 | domains of variables with negative coefficients. Similarly, if the largest |
---|
2068 | positive gap (max{j in P} a<ij>(u<j>-l<j>) is greater than b<i>-L<i>, we |
---|
2069 | can't shrink the domain of any variable with a positive coefficient. |
---|
2070 | |
---|
2071 | There are analogous conditions using U<i> and blow<i>. |
---|
2072 | |
---|
2073 | Note that we don't update the largest negative and positive gaps, so as |
---|
2074 | propagation grinds on, this simple test becomes less accurate. On the other |
---|
2075 | hand, the gaps (b<i>-L<i>) and (U<i>-blow<i>) are also shrinking. Still, it's |
---|
2076 | another justification for the strict limits on propagation. |
---|
2077 | |
---|
2078 | Summary: |
---|
2079 | |
---|
2080 | minR = SUM{P}(a<ij>l<j>) + SUM{N}(a<ij>u<j>) |
---|
2081 | maxR = SUM{P}(a<ij>u<j>) + SUM{N}(a<ij>l<j>) |
---|
2082 | |
---|
2083 | doRowUpN: a<ij> < 0, minR can violate rowUp => increase l<j> (more negative) |
---|
2084 | doRowLoN: a<ij> < 0, maxR can violate rowLo => decrease u<j> (more positive) |
---|
2085 | doRowUpP: a<ij> > 0, minR can violate rowUp => decrease u<j> (less positive) |
---|
2086 | doRowLoP: a<ij> > 0, maxR can violate rowLo => increase l<j> (less negative) |
---|
2087 | |
---|
2088 | Note that the bound (l<j>, u<j>) that can be tightened in any given |
---|
2089 | situation is *not* the bound involved in minR or maxR. |
---|
2090 | |
---|
2091 | For binary variables, of course, `shrink the domain' corresponds to fixing |
---|
2092 | the variable. |
---|
2093 | */ |
---|
2094 | |
---|
2095 | if (rowUp < 1.0e10) { |
---|
2096 | doRowUpN = true ; |
---|
2097 | doRowUpP = true ; |
---|
2098 | rowUp2 = rowUp-minR[irow] ; |
---|
2099 | if (rowUp2 < -primalTolerance_) { |
---|
2100 | notFeasible = true ; |
---|
2101 | break ; |
---|
2102 | } else { |
---|
2103 | if (rowUp2+largestNegativeInRow[irow] > 0) |
---|
2104 | doRowUpN = false ; |
---|
2105 | if (rowUp2-largestPositiveInRow[irow] > 0) |
---|
2106 | doRowUpP = false ; |
---|
2107 | } |
---|
2108 | } else { |
---|
2109 | doRowUpN = false ; |
---|
2110 | doRowUpP = false ; |
---|
2111 | rowUp2 = COIN_DBL_MAX ; |
---|
2112 | } |
---|
2113 | double rowLo = rowLower[irow] ; |
---|
2114 | double rowLo2 = 0.0 ; |
---|
2115 | bool doRowLoN ; |
---|
2116 | bool doRowLoP ; |
---|
2117 | if (rowLo > -1.0e10) { |
---|
2118 | doRowLoN = true ; |
---|
2119 | doRowLoP = true ; |
---|
2120 | rowLo2 = rowLo-maxR[irow] ; |
---|
2121 | if (rowLo2 > primalTolerance_) { |
---|
2122 | notFeasible = true ; |
---|
2123 | break ; |
---|
2124 | } else { |
---|
2125 | if (rowLo2-largestNegativeInRow[irow] < 0) |
---|
2126 | doRowLoN = false ; |
---|
2127 | if (rowLo2+largestPositiveInRow[irow] < 0) |
---|
2128 | doRowLoP = false ; |
---|
2129 | } |
---|
2130 | } else { |
---|
2131 | doRowLoN = false ; |
---|
2132 | doRowLoP = false ; |
---|
2133 | rowLo2 = -COIN_DBL_MAX ; |
---|
2134 | } |
---|
2135 | markR[irow] += 10000 ; |
---|
2136 | /* |
---|
2137 | LOU_DEBUG: see if we're processing again without a bound change. |
---|
2138 | if (markR[irow]/10000 > 1) |
---|
2139 | std::cout |
---|
2140 | << "REDUNDANT: processing " << irow |
---|
2141 | << " for the " << markR[irow]/10000 << " time." << std::endl ; |
---|
2142 | */ |
---|
2143 | /* |
---|
2144 | PROBE LOOP: WALK ROW |
---|
2145 | |
---|
2146 | Given the above analysis, work on the columns with negative coefficients. |
---|
2147 | |
---|
2148 | doRowUpN: a<ij> < 0, minR can violate rowUp => increase l<j> (more negative) |
---|
2149 | doRowLoN: a<ij> < 0, maxR can violate rowLo => decrease u<j> (more positive) |
---|
2150 | |
---|
2151 | TODO: The asserts here should be compiled out unless we're debugging. |
---|
2152 | -- lh, 101210 -- |
---|
2153 | */ |
---|
2154 | if (doRowUpN && doRowLoN) { |
---|
2155 | //doRowUpN=doRowLoN=false ; |
---|
2156 | // Start neg values loop |
---|
2157 | for (int kk = rStart ; kk < rEnd ; kk++) { |
---|
2158 | int kcol = column[kk] ; |
---|
2159 | int markIt = markC[kcol] ; |
---|
2160 | // Skip columns with fixed variables. |
---|
2161 | if ((markIt&(tightenLower|tightenUpper)) != (tightenLower|tightenUpper)) { |
---|
2162 | double value2 = rowElements[kk] ; |
---|
2163 | if (colUpper[kcol] <= 1e10) |
---|
2164 | assert ((markIt&infUpper) == 0) ; |
---|
2165 | else |
---|
2166 | assert ((markIt&infUpper) != 0) ; |
---|
2167 | if (colLower[kcol] >= -1e10) |
---|
2168 | assert ((markIt&infLower) == 0) ; |
---|
2169 | else |
---|
2170 | assert ((markIt&infLower) != 0) ; |
---|
2171 | assert (value2 < 0.0) ; |
---|
2172 | /* |
---|
2173 | Not every column will be productive. Can we do anything with this one? |
---|
2174 | */ |
---|
2175 | double gap = columnGap[kcol]*value2 ; |
---|
2176 | bool doUp = (rowUp2+gap < 0.0) ; |
---|
2177 | bool doDown = (rowLo2-gap > 0.0) ; |
---|
2178 | if (doUp || doDown) { |
---|
2179 | double moveUp = 0.0 ; |
---|
2180 | double moveDown = 0.0 ; |
---|
2181 | double newUpper = -1.0 ; |
---|
2182 | double newLower = 1.0 ; |
---|
2183 | /* |
---|
2184 | doUp attempts to increase the lower bound. The math looks like this: |
---|
2185 | a<irow,kcol>x<kcol> + (minR[irow]-a<irow,kcol>u<kcol>) <= b<irow> |
---|
2186 | value2*x<kcol> - value2*u<kcol> <= (b<irow>-minR[irow]) = rowUp2 |
---|
2187 | x<kcol> - u<kcol> >= rowUp2/value2 |
---|
2188 | l<kcol> = u<kcol> + rowUp2/value2 |
---|
2189 | Hence we need u<kcol> to be finite (0x8 not set). We also refuse to tighten |
---|
2190 | l<kcol> a second time (0x2 not set). |
---|
2191 | |
---|
2192 | TODO: So, why don't we allow tightening a second time? Clearly we might |
---|
2193 | process a row on some other iteration that imposes a tighter bound. |
---|
2194 | Is this an optimisation for binary variables, or is there some |
---|
2195 | deeper issue at work in terms of correctness? -- lh, 101127 -- |
---|
2196 | |
---|
2197 | TODO: Also, it's clear that this bit of code would like to handle |
---|
2198 | continuous variables, but it's also clear that it does it incorrectly. |
---|
2199 | When the new lower bound straddles the existing upper bound, that |
---|
2200 | looks to me to be sufficient justification to fix the variable. But |
---|
2201 | the code here does just the opposite: merely tightening the bound is |
---|
2202 | flagged as fixing the variable, with no consideration that we may have |
---|
2203 | gone infeasible. But if we prove we straddle the upper bound, all we |
---|
2204 | claim is we've tightened the upper bound! -- lh, 101127 -- |
---|
2205 | */ |
---|
2206 | if (doUp && ((markIt&(tightenLower|infUpper)) == 0)) { |
---|
2207 | double dbound = colUpper[kcol]+rowUp2/value2 ; |
---|
2208 | if (intVar[kcol]) { |
---|
2209 | markIt |= tightenLower ; |
---|
2210 | newLower = ceil(dbound-primalTolerance_) ; |
---|
2211 | } else { |
---|
2212 | newLower = dbound ; |
---|
2213 | if (newLower+primalTolerance_ > colUpper[kcol] && |
---|
2214 | newLower-primalTolerance_ <= colUpper[kcol]) { |
---|
2215 | newLower = colUpper[kcol] ; |
---|
2216 | markIt |= tightenLower ; |
---|
2217 | //markIt = (tightenUpper|tightenLower) ; |
---|
2218 | } else { |
---|
2219 | // avoid problems - fix later ? |
---|
2220 | markIt = (tightenUpper|tightenLower) ; |
---|
2221 | } |
---|
2222 | } |
---|
2223 | moveUp = newLower-colLower[kcol] ; |
---|
2224 | } |
---|
2225 | /* |
---|
2226 | The same, but attempt to decrease the upper bound by comparison with the |
---|
2227 | blow for the row. |
---|
2228 | */ |
---|
2229 | if (doDown && ((markIt&(tightenUpper|infLower)) == 0)) { |
---|
2230 | double dbound = colLower[kcol] + rowLo2/value2 ; |
---|
2231 | if (intVar[kcol]) { |
---|
2232 | markIt |= tightenUpper ; |
---|
2233 | newUpper = floor(dbound+primalTolerance_) ; |
---|
2234 | } else { |
---|
2235 | newUpper = dbound ; |
---|
2236 | if (newUpper-primalTolerance_ < colLower[kcol] && |
---|
2237 | newUpper+primalTolerance_ >= colLower[kcol]) { |
---|
2238 | newUpper = colLower[kcol] ; |
---|
2239 | markIt |= tightenUpper ; |
---|
2240 | //markIt = (tightenUpper|tightenLower) ; |
---|
2241 | } else { |
---|
2242 | // avoid problems - fix later ? |
---|
2243 | markIt = (tightenUpper|tightenLower) ; |
---|
2244 | } |
---|
2245 | } |
---|
2246 | moveDown = newUpper-colUpper[kcol] ; |
---|
2247 | } |
---|
2248 | /* |
---|
2249 | Is there any change to propagate? Assuming we haven't exceeded the limit |
---|
2250 | for propagating, queue the variable. (If the variable's already on the list, |
---|
2251 | that's not a problem.) |
---|
2252 | |
---|
2253 | Note that markIt is initially loaded from markC and whatever bounds are |
---|
2254 | changed are then or'd in. So changes accumulate. |
---|
2255 | */ |
---|
2256 | if (!moveUp && !moveDown) |
---|
2257 | continue ; |
---|
2258 | bool onList = ((markC[kcol]&(tightenLower|tightenUpper)) != 0) ; |
---|
2259 | if (nstackC < 2*maxStack) { |
---|
2260 | markC[kcol] = markIt ; |
---|
2261 | } |
---|
2262 | if (moveUp && (nstackC < 2*maxStack)) { |
---|
2263 | if (!onList) { |
---|
2264 | stackC[nstackC] = kcol ; |
---|
2265 | saveL[nstackC] = colLower[kcol] ; |
---|
2266 | saveU[nstackC] = colUpper[kcol] ; |
---|
2267 | assert (saveU[nstackC] > saveL[nstackC]) ; |
---|
2268 | assert (nstackC < nCols) ; |
---|
2269 | nstackC++ ; |
---|
2270 | onList = true ; |
---|
2271 | } |
---|
2272 | /* |
---|
2273 | The first assert here says `If we're minimising, and d<j> < 0, then |
---|
2274 | x<j> must be NBUB, so if the new l<j> > x*<j> = u<j>, we're infeasible.' |
---|
2275 | |
---|
2276 | TODO: Why haven't we detected infeasibility? Put another way, why isn't |
---|
2277 | there an assert(notFeasible)? -- lh, 101127 -- |
---|
2278 | |
---|
2279 | TODO: Seems like the change in obj should be ((new l<j>) - x*<j>)*d<j>, but |
---|
2280 | moveUp is (newLower-colLower). Was there some guarantee that x*<j> |
---|
2281 | equals l<j>? -- lh, 101127 -- |
---|
2282 | TODO: I think this is another symptom of incomplete conversion for general |
---|
2283 | integers. For binary variables, it's not an issue. -- lh, 101203 -- |
---|
2284 | */ |
---|
2285 | if (newLower > colsol[kcol]) { |
---|
2286 | if (djs[kcol] < 0.0) { |
---|
2287 | assert (newLower > colUpper[kcol]+primalTolerance_) ; |
---|
2288 | } else { |
---|
2289 | objVal += moveUp*djs[kcol] ; |
---|
2290 | } |
---|
2291 | } |
---|
2292 | /* |
---|
2293 | TODO: Up above we've already set newLower = ceil(dbound-primalTolerance_). |
---|
2294 | It's highly unlikely that ceil(newLower-1.0e4) will be different. |
---|
2295 | -- lh, 101127 -- |
---|
2296 | TODO: My first inclination is to claim that newLower should never be worse |
---|
2297 | than the existing bound. But I'd better extend the benefit of the |
---|
2298 | doubt for a bit. Some rows will produce stronger bounds than others. |
---|
2299 | -- lh, 101127 -- |
---|
2300 | */ |
---|
2301 | if (intVar[kcol]) |
---|
2302 | newLower = CoinMax(colLower[kcol], |
---|
2303 | ceil(newLower-1.0e-4)) ; |
---|
2304 | colLower[kcol] = newLower ; |
---|
2305 | columnGap[kcol] = colUpper[kcol]-newLower-primalTolerance_ ; |
---|
2306 | /* |
---|
2307 | Propagate the column bound change. |
---|
2308 | |
---|
2309 | TODO: Notice that we're not setting a variable here when we discover |
---|
2310 | infeasibility; rather, we're setting the column bound to a ridiculous |
---|
2311 | value which will be discovered farther down. -- lh, 101127 -- |
---|
2312 | TODO: Replacing code to update row bounds with updateRowBounds(). We'll |
---|
2313 | see how it looks. -- lh, 101203 -- |
---|
2314 | */ |
---|
2315 | if (CoinAbs(colUpper[kcol]-colLower[kcol]) < 1.0e-6) { |
---|
2316 | markC[kcol] = (tightenLower|tightenUpper) ; |
---|
2317 | } |
---|
2318 | markC[kcol] &= ~(infLower|infUpper) ; |
---|
2319 | if (colUpper[kcol] > 1.0e10) |
---|
2320 | markC[kcol] |= infUpper ; |
---|
2321 | if (colLower[kcol] < -1.0e10) |
---|
2322 | markC[kcol] |= infLower ; |
---|
2323 | |
---|
2324 | if (!updateRowBounds(kcol,moveUp, |
---|
2325 | columnStart,columnLength,row,columnElements, |
---|
2326 | rowUpper,rowLower,nstackR,stackR,markR, |
---|
2327 | minR,maxR,saveMin,saveMax)) { |
---|
2328 | colLower[kcol] = 1.0e50 ; |
---|
2329 | } |
---|
2330 | } |
---|
2331 | /* |
---|
2332 | Repeat the whole business, in the down direction. Stack the variable, if |
---|
2333 | it's not already there, do the infeasibility check / objective update, and |
---|
2334 | update minR / maxR for affected constraints. |
---|
2335 | */ |
---|
2336 | if (moveDown&&nstackC<2*maxStack) { |
---|
2337 | if (!onList) { |
---|
2338 | stackC[nstackC]=kcol ; |
---|
2339 | saveL[nstackC]=colLower[kcol] ; |
---|
2340 | saveU[nstackC]=colUpper[kcol] ; |
---|
2341 | assert (saveU[nstackC]>saveL[nstackC]) ; |
---|
2342 | assert (nstackC<nCols) ; |
---|
2343 | nstackC++ ; |
---|
2344 | onList=true ; |
---|
2345 | } |
---|
2346 | if (newUpper<colsol[kcol]) { |
---|
2347 | if (djs[kcol]>0.0) { |
---|
2348 | /* should be infeasible */ |
---|
2349 | assert (colLower[kcol]>newUpper+primalTolerance_) ; |
---|
2350 | } else { |
---|
2351 | objVal += moveDown*djs[kcol] ; |
---|
2352 | } |
---|
2353 | } |
---|
2354 | if (intVar[kcol]) |
---|
2355 | newUpper = CoinMin(colUpper[kcol],floor(newUpper+1.0e-4)) ; |
---|
2356 | colUpper[kcol]=newUpper ; |
---|
2357 | columnGap[kcol] = newUpper-colLower[kcol]-primalTolerance_ ; |
---|
2358 | if (CoinAbs(colUpper[kcol]-colLower[kcol])<1.0e-6) { |
---|
2359 | markC[kcol] = (tightenLower|tightenUpper); // say fixed |
---|
2360 | } |
---|
2361 | markC[kcol] &= ~(infLower|infUpper) ; |
---|
2362 | if (colUpper[kcol]>1.0e10) |
---|
2363 | markC[kcol] |= infUpper ; |
---|
2364 | if (colLower[kcol]<-1.0e10) |
---|
2365 | markC[kcol] |= infLower ; |
---|
2366 | |
---|
2367 | if (!updateRowBounds(kcol,moveDown, |
---|
2368 | columnStart,columnLength,row,columnElements, |
---|
2369 | rowUpper,rowLower,nstackR,stackR,markR, |
---|
2370 | minR,maxR,saveMin,saveMax)) { |
---|
2371 | colUpper[kcol] = -1.0e50 ; |
---|
2372 | } |
---|
2373 | } |
---|
2374 | /* |
---|
2375 | We've propagated the bound changes for this column. Check to see if we |
---|
2376 | discovered infeasibility. Abort by claiming we've cleared the propagation |
---|
2377 | stack. |
---|
2378 | */ |
---|
2379 | if (colLower[kcol] > colUpper[kcol]+primalTolerance_) { |
---|
2380 | notFeasible = true ; |
---|
2381 | k = columnStart[jcol]+columnLength[jcol] ; |
---|
2382 | istackC = nstackC+1 ; |
---|
2383 | break ; |
---|
2384 | } |
---|
2385 | } // end if (doUp || doDown) (productive column) |
---|
2386 | } // end column not fixed |
---|
2387 | } // end loop on negative coefficients of this row |
---|
2388 | } else if (doRowUpN) { |
---|
2389 | /* |
---|
2390 | The above block propagated change for a variable with a negative coefficient, |
---|
2391 | where we could work against both the row upper and lower bounds. Now we're |
---|
2392 | going to do it again, but specialised for the case where we can only work |
---|
2393 | against the upper bound. |
---|
2394 | */ |
---|
2395 | // Start neg values loop |
---|
2396 | for (int kk = rStart ; kk < rEnd ; kk++) { |
---|
2397 | int kcol =column[kk] ; |
---|
2398 | int markIt=markC[kcol] ; |
---|
2399 | if ((markIt&(tightenLower|tightenUpper))!= (tightenLower|tightenUpper)) { |
---|
2400 | double value2=rowElements[kk] ; |
---|
2401 | double gap = columnGap[kcol]*value2 ; |
---|
2402 | if (!(rowUp2 + gap < 0.0)) |
---|
2403 | continue ; |
---|
2404 | double moveUp=0.0 ; |
---|
2405 | double newLower=1.0 ; |
---|
2406 | if ((markIt&(tightenLower|infUpper))==0) { |
---|
2407 | double dbound = colUpper[kcol]+rowUp2/value2 ; |
---|
2408 | if (intVar[kcol]) { |
---|
2409 | markIt |= tightenLower ; |
---|
2410 | newLower = ceil(dbound-primalTolerance_) ; |
---|
2411 | } else { |
---|
2412 | newLower=dbound ; |
---|
2413 | if (newLower+primalTolerance_>colUpper[kcol]&& |
---|
2414 | newLower-primalTolerance_<=colUpper[kcol]) { |
---|
2415 | newLower=colUpper[kcol] ; |
---|
2416 | markIt |= tightenLower ; |
---|
2417 | //markIt= (tightenLower|tightenUpper) ; |
---|
2418 | } else { |
---|
2419 | // avoid problems - fix later ? |
---|
2420 | markIt= (tightenLower|tightenUpper) ; |
---|
2421 | } |
---|
2422 | } |
---|
2423 | moveUp = newLower-colLower[kcol] ; |
---|
2424 | if (!moveUp) |
---|
2425 | continue ; |
---|
2426 | bool onList = ((markC[kcol]&(tightenLower|tightenUpper))!=0) ; |
---|
2427 | if (nstackC<2*maxStack) { |
---|
2428 | markC[kcol] = markIt ; |
---|
2429 | } |
---|
2430 | if (moveUp&&nstackC<2*maxStack) { |
---|
2431 | if (!onList) { |
---|
2432 | stackC[nstackC]=kcol ; |
---|
2433 | saveL[nstackC]=colLower[kcol] ; |
---|
2434 | saveU[nstackC]=colUpper[kcol] ; |
---|
2435 | assert (saveU[nstackC]>saveL[nstackC]) ; |
---|
2436 | assert (nstackC<nCols) ; |
---|
2437 | nstackC++ ; |
---|
2438 | onList=true ; |
---|
2439 | } |
---|
2440 | if (newLower>colsol[kcol]) { |
---|
2441 | if (djs[kcol]<0.0) { |
---|
2442 | /* should be infeasible */ |
---|
2443 | assert (newLower>colUpper[kcol]+primalTolerance_) ; |
---|
2444 | } else { |
---|
2445 | objVal += moveUp*djs[kcol] ; |
---|
2446 | } |
---|
2447 | } |
---|
2448 | if (intVar[kcol]) |
---|
2449 | newLower = CoinMax(colLower[kcol],ceil(newLower-1.0e-4)) ; |
---|
2450 | colLower[kcol]=newLower ; |
---|
2451 | columnGap[kcol] = colUpper[kcol]-newLower-primalTolerance_ ; |
---|
2452 | if (CoinAbs(colUpper[kcol]-colLower[kcol])<1.0e-6) { |
---|
2453 | markC[kcol]= (tightenLower|tightenUpper); // say fixed |
---|
2454 | } |
---|
2455 | markC[kcol] &= ~(infLower|infUpper) ; |
---|
2456 | if (colUpper[kcol]>1.0e10) |
---|
2457 | markC[kcol] |= infUpper ; |
---|
2458 | if (colLower[kcol]<-1.0e10) |
---|
2459 | markC[kcol] |= infLower ; |
---|
2460 | |
---|
2461 | if (!updateRowBounds(kcol,moveUp, |
---|
2462 | columnStart,columnLength,row,columnElements, |
---|
2463 | rowUpper,rowLower,nstackR,stackR,markR, |
---|
2464 | minR,maxR,saveMin,saveMax)) { |
---|
2465 | colLower[kcol] = 1.0e50 ; |
---|
2466 | } |
---|
2467 | } |
---|
2468 | if (colLower[kcol] > colUpper[kcol]+primalTolerance_) { |
---|
2469 | notFeasible = true ; |
---|
2470 | k = columnStart[jcol]+columnLength[jcol] ; |
---|
2471 | istackC = nstackC+1 ; |
---|
2472 | break ; |
---|
2473 | } |
---|
2474 | } |
---|
2475 | } |
---|
2476 | } // end big loop rStart->rPos |
---|
2477 | } else if (doRowLoN) { |
---|
2478 | /* |
---|
2479 | And yet again, for the case where we can only work against the lower bound. |
---|
2480 | */ |
---|
2481 | // Start neg values loop |
---|
2482 | for (int kk = rStart ; kk < rEnd ; kk++) { |
---|
2483 | int kcol =column[kk] ; |
---|
2484 | if ((markC[kcol]&(tightenLower|tightenUpper))!= (tightenLower|tightenUpper)) { |
---|
2485 | double moveDown=0.0 ; |
---|
2486 | double newUpper=-1.0 ; |
---|
2487 | double value2=rowElements[kk] ; |
---|
2488 | int markIt=markC[kcol] ; |
---|
2489 | assert (value2<0.0) ; |
---|
2490 | double gap = columnGap[kcol]*value2 ; |
---|
2491 | bool doDown = (rowLo2 -gap > 0.0) ; |
---|
2492 | if (doDown && ((markIt&(tightenUpper|infLower)) == 0)) { |
---|
2493 | double dbound = colLower[kcol] + rowLo2/value2 ; |
---|
2494 | if (intVar[kcol]) { |
---|
2495 | markIt |= tightenUpper ; |
---|
2496 | newUpper = floor(dbound+primalTolerance_) ; |
---|
2497 | } else { |
---|
2498 | newUpper=dbound ; |
---|
2499 | if (newUpper-primalTolerance_<colLower[kcol]&& |
---|
2500 | newUpper+primalTolerance_>=colLower[kcol]) { |
---|
2501 | newUpper=colLower[kcol] ; |
---|
2502 | markIt |= tightenUpper ; |
---|
2503 | //markIt= (tightenLower|tightenUpper) ; |
---|
2504 | } else { |
---|
2505 | // avoid problems - fix later ? |
---|
2506 | markIt= (tightenLower|tightenUpper) ; |
---|
2507 | } |
---|
2508 | } |
---|
2509 | moveDown = newUpper-colUpper[kcol] ; |
---|
2510 | if (!moveDown) |
---|
2511 | continue ; |
---|
2512 | bool onList = ((markC[kcol]&(tightenLower|tightenUpper))!=0) ; |
---|
2513 | if (nstackC<2*maxStack) { |
---|
2514 | markC[kcol] = markIt ; |
---|
2515 | } |
---|
2516 | if (moveDown&&nstackC<2*maxStack) { |
---|
2517 | if (!onList) { |
---|
2518 | stackC[nstackC]=kcol ; |
---|
2519 | saveL[nstackC]=colLower[kcol] ; |
---|
2520 | saveU[nstackC]=colUpper[kcol] ; |
---|
2521 | assert (saveU[nstackC]>saveL[nstackC]) ; |
---|
2522 | assert (nstackC<nCols) ; |
---|
2523 | nstackC++ ; |
---|
2524 | onList=true ; |
---|
2525 | } |
---|
2526 | if (newUpper<colsol[kcol]) { |
---|
2527 | if (djs[kcol]>0.0) { |
---|
2528 | /* should be infeasible */ |
---|
2529 | assert (colLower[kcol]>newUpper+primalTolerance_) ; |
---|
2530 | } else { |
---|
2531 | objVal += moveDown*djs[kcol] ; |
---|
2532 | } |
---|
2533 | } |
---|
2534 | if (intVar[kcol]) |
---|
2535 | newUpper = CoinMin(colUpper[kcol],floor(newUpper+1.0e-4)) ; |
---|
2536 | colUpper[kcol]=newUpper ; |
---|
2537 | columnGap[kcol] = newUpper-colLower[kcol]-primalTolerance_ ; |
---|
2538 | if (CoinAbs(colUpper[kcol]-colLower[kcol])<1.0e-6) { |
---|
2539 | markC[kcol]= (tightenLower|tightenUpper); // say fixed |
---|
2540 | } |
---|
2541 | markC[kcol] &= ~(infLower|infUpper) ; |
---|
2542 | if (colUpper[kcol]>1.0e10) |
---|
2543 | markC[kcol] |= infUpper ; |
---|
2544 | if (colLower[kcol]<-1.0e10) |
---|
2545 | markC[kcol] |= infLower ; |
---|
2546 | |
---|
2547 | if (!updateRowBounds(kcol,moveDown, |
---|
2548 | columnStart,columnLength,row,columnElements, |
---|
2549 | rowUpper,rowLower,nstackR,stackR,markR, |
---|
2550 | minR,maxR,saveMin,saveMax)) { |
---|
2551 | colUpper[kcol] = -1.0e50 ; |
---|
2552 | } |
---|
2553 | } |
---|
2554 | if (colLower[kcol] > colUpper[kcol]+primalTolerance_) { |
---|
2555 | notFeasible = true ; |
---|
2556 | k = columnStart[jcol]+columnLength[jcol] ; |
---|
2557 | istackC = nstackC+1 ; |
---|
2558 | break ; |
---|
2559 | } |
---|
2560 | } |
---|
2561 | } |
---|
2562 | } // end big loop rStart->rPos |
---|
2563 | } |
---|
2564 | /* |
---|
2565 | We've finished working on the negative coefficients of the row. Advance |
---|
2566 | rStart and rEnd to cover the positive coefficients and repeat the previous |
---|
2567 | 500 lines. |
---|
2568 | */ |
---|
2569 | rStart = rEnd ; |
---|
2570 | rEnd = rowStart[irow+1] ; |
---|
2571 | if (doRowUpP&&doRowLoP) { |
---|
2572 | //doRowUpP=doRowLoP=false ; |
---|
2573 | // Start pos values loop |
---|
2574 | for (int kk =rStart;kk<rEnd;kk++) { |
---|
2575 | int kcol=column[kk] ; |
---|
2576 | int markIt=markC[kcol] ; |
---|
2577 | if ((markIt&(tightenLower|tightenUpper))!= (tightenLower|tightenUpper)) { |
---|
2578 | double value2=rowElements[kk] ; |
---|
2579 | assert (value2 > 0.0) ; |
---|
2580 | /* positive element */ |
---|
2581 | double gap = columnGap[kcol]*value2 ; |
---|
2582 | bool doDown = (rowLo2 + gap > 0.0) ; |
---|
2583 | bool doUp = (rowUp2 - gap < 0.0) ; |
---|
2584 | if (doDown||doUp) { |
---|
2585 | double moveUp=0.0 ; |
---|
2586 | double moveDown=0.0 ; |
---|
2587 | double newUpper=-1.0 ; |
---|
2588 | double newLower=1.0 ; |
---|
2589 | if (doDown && ((markIt&(tightenLower|infUpper)) == 0)) { |
---|
2590 | double dbound = colUpper[kcol] + rowLo2/value2 ; |
---|
2591 | if (intVar[kcol]) { |
---|
2592 | markIt |= tightenLower ; |
---|
2593 | newLower = ceil(dbound-primalTolerance_) ; |
---|
2594 | } else { |
---|
2595 | newLower=dbound ; |
---|
2596 | if (newLower+primalTolerance_>colUpper[kcol]&& |
---|
2597 | newLower-primalTolerance_<=colUpper[kcol]) { |
---|
2598 | newLower=colUpper[kcol] ; |
---|
2599 | markIt |= tightenLower ; |
---|
2600 | //markIt= (tightenLower|tightenUpper) ; |
---|
2601 | } else { |
---|
2602 | // avoid problems - fix later ? |
---|
2603 | markIt= (tightenLower|tightenUpper) ; |
---|
2604 | } |
---|
2605 | } |
---|
2606 | moveUp = newLower-colLower[kcol] ; |
---|
2607 | } |
---|
2608 | if (doUp && ((markIt&(tightenUpper|infLower)) == 0)) { |
---|
2609 | double dbound = colLower[kcol] + rowUp2/value2 ; |
---|
2610 | if (intVar[kcol]) { |
---|
2611 | markIt |= tightenUpper ; |
---|
2612 | newUpper = floor(dbound+primalTolerance_) ; |
---|
2613 | } else { |
---|
2614 | newUpper=dbound ; |
---|
2615 | if (newUpper-primalTolerance_<colLower[kcol]&& |
---|
2616 | newUpper+primalTolerance_>=colLower[kcol]) { |
---|
2617 | newUpper=colLower[kcol] ; |
---|
2618 | markIt |= tightenUpper ; |
---|
2619 | //markIt= (tightenLower|tightenUpper) ; |
---|
2620 | } else { |
---|
2621 | // avoid problems - fix later ? |
---|
2622 | markIt= (tightenLower|tightenUpper) ; |
---|
2623 | } |
---|
2624 | } |
---|
2625 | moveDown = newUpper-colUpper[kcol] ; |
---|
2626 | } |
---|
2627 | if (!moveUp&&!moveDown) |
---|
2628 | continue ; |
---|
2629 | bool onList = ((markC[kcol]&(tightenLower|tightenUpper))!=0) ; |
---|
2630 | if (nstackC<2*maxStack) { |
---|
2631 | markC[kcol] = markIt ; |
---|
2632 | } |
---|
2633 | if (moveUp&&nstackC<2*maxStack) { |
---|
2634 | if (!onList) { |
---|
2635 | stackC[nstackC]=kcol ; |
---|
2636 | saveL[nstackC]=colLower[kcol] ; |
---|
2637 | saveU[nstackC]=colUpper[kcol] ; |
---|
2638 | assert (saveU[nstackC]>saveL[nstackC]) ; |
---|
2639 | assert (nstackC<nCols) ; |
---|
2640 | nstackC++ ; |
---|
2641 | onList=true ; |
---|
2642 | } |
---|
2643 | if (newLower>colsol[kcol]) { |
---|
2644 | if (djs[kcol]<0.0) { |
---|
2645 | /* should be infeasible */ |
---|
2646 | assert (newLower>colUpper[kcol]+primalTolerance_) ; |
---|
2647 | } else { |
---|
2648 | objVal += moveUp*djs[kcol] ; |
---|
2649 | } |
---|
2650 | } |
---|
2651 | if (intVar[kcol]) |
---|
2652 | newLower = CoinMax(colLower[kcol],ceil(newLower-1.0e-4)) ; |
---|
2653 | colLower[kcol]=newLower ; |
---|
2654 | columnGap[kcol] = colUpper[kcol]-newLower-primalTolerance_ ; |
---|
2655 | if (CoinAbs(colUpper[kcol]-colLower[kcol])<1.0e-6) { |
---|
2656 | markC[kcol]= (tightenLower|tightenUpper); // say fixed |
---|
2657 | } |
---|
2658 | markC[kcol] &= ~(infLower|infUpper) ; |
---|
2659 | if (colUpper[kcol]>1.0e10) |
---|
2660 | markC[kcol] |= infUpper ; |
---|
2661 | if (colLower[kcol]<-1.0e10) |
---|
2662 | markC[kcol] |= infLower ; |
---|
2663 | |
---|
2664 | if (!updateRowBounds(kcol,moveUp, |
---|
2665 | columnStart,columnLength,row,columnElements, |
---|
2666 | rowUpper,rowLower,nstackR,stackR,markR, |
---|
2667 | minR,maxR,saveMin,saveMax)) { |
---|
2668 | colLower[kcol] = 1.0e50 ; |
---|
2669 | } |
---|
2670 | } |
---|
2671 | if (moveDown&&nstackC<2*maxStack) { |
---|
2672 | if (!onList) { |
---|
2673 | stackC[nstackC]=kcol ; |
---|
2674 | saveL[nstackC]=colLower[kcol] ; |
---|
2675 | saveU[nstackC]=colUpper[kcol] ; |
---|
2676 | assert (saveU[nstackC]>saveL[nstackC]) ; |
---|
2677 | assert (nstackC<nCols) ; |
---|
2678 | nstackC++ ; |
---|
2679 | onList=true ; |
---|
2680 | } |
---|
2681 | if (newUpper<colsol[kcol]) { |
---|
2682 | if (djs[kcol]>0.0) { |
---|
2683 | /* should be infeasible */ |
---|
2684 | assert (colLower[kcol]>newUpper+primalTolerance_) ; |
---|
2685 | } else { |
---|
2686 | objVal += moveDown*djs[kcol] ; |
---|
2687 | } |
---|
2688 | } |
---|
2689 | if (intVar[kcol]) |
---|
2690 | newUpper = CoinMin(colUpper[kcol],floor(newUpper+1.0e-4)) ; |
---|
2691 | colUpper[kcol]=newUpper ; |
---|
2692 | columnGap[kcol] = newUpper-colLower[kcol]-primalTolerance_ ; |
---|
2693 | if (CoinAbs(colUpper[kcol]-colLower[kcol])<1.0e-6) { |
---|
2694 | markC[kcol]= (tightenLower|tightenUpper); // say fixed |
---|
2695 | } |
---|
2696 | markC[kcol] &= ~(infLower|infUpper) ; |
---|
2697 | if (colUpper[kcol]>1.0e10) |
---|
2698 | markC[kcol] |= infUpper ; |
---|
2699 | if (colLower[kcol]<-1.0e10) |
---|
2700 | markC[kcol] |= infLower ; |
---|
2701 | |
---|
2702 | if (!updateRowBounds(kcol,moveDown, |
---|
2703 | columnStart,columnLength,row,columnElements, |
---|
2704 | rowUpper,rowLower,nstackR,stackR,markR, |
---|
2705 | minR,maxR,saveMin,saveMax)) { |
---|
2706 | colUpper[kcol] = -1.0e50 ; |
---|
2707 | } |
---|
2708 | } |
---|
2709 | if (colLower[kcol]>colUpper[kcol]+primalTolerance_) { |
---|
2710 | notFeasible = true; ; |
---|
2711 | k=columnStart[jcol]+columnLength[jcol] ; |
---|
2712 | istackC=nstackC+1 ; |
---|
2713 | break ; |
---|
2714 | } |
---|
2715 | } |
---|
2716 | } |
---|
2717 | } // end big loop rPos->rEnd |
---|
2718 | } else if (doRowUpP) { |
---|
2719 | // Start pos values loop |
---|
2720 | for (int kk =rStart;kk<rEnd;kk++) { |
---|
2721 | int kcol =column[kk] ; |
---|
2722 | int markIt=markC[kcol] ; |
---|
2723 | if ((markIt&(tightenLower|tightenUpper))!= (tightenLower|tightenUpper)) { |
---|
2724 | double value2=rowElements[kk] ; |
---|
2725 | assert (value2 > 0.0) ; |
---|
2726 | /* positive element */ |
---|
2727 | double gap = columnGap[kcol]*value2 ; |
---|
2728 | bool doUp = (rowUp2 - gap < 0.0) ; |
---|
2729 | if (doUp && ((markIt&(tightenUpper|infLower)) == 0)) { |
---|
2730 | double newUpper=-1.0 ; |
---|
2731 | double dbound = colLower[kcol] + rowUp2/value2 ; |
---|
2732 | if (intVar[kcol]) { |
---|
2733 | markIt |= tightenUpper ; |
---|
2734 | newUpper = floor(dbound+primalTolerance_) ; |
---|
2735 | } else { |
---|
2736 | newUpper=dbound ; |
---|
2737 | if (newUpper-primalTolerance_<colLower[kcol]&& |
---|
2738 | newUpper+primalTolerance_>=colLower[kcol]) { |
---|
2739 | newUpper=colLower[kcol] ; |
---|
2740 | markIt |= tightenUpper ; |
---|
2741 | //markIt= (tightenLower|tightenUpper) ; |
---|
2742 | } else { |
---|
2743 | // avoid problems - fix later ? |
---|
2744 | markIt= (tightenLower|tightenUpper) ; |
---|
2745 | } |
---|
2746 | } |
---|
2747 | double moveDown = newUpper-colUpper[kcol] ; |
---|
2748 | if (!moveDown) |
---|
2749 | continue ; |
---|
2750 | bool onList = ((markC[kcol]&(tightenLower|tightenUpper))!=0) ; |
---|
2751 | if (nstackC<2*maxStack) { |
---|
2752 | markC[kcol] = markIt ; |
---|
2753 | } |
---|
2754 | if (nstackC<2*maxStack) { |
---|
2755 | if (!onList) { |
---|
2756 | stackC[nstackC]=kcol ; |
---|
2757 | saveL[nstackC]=colLower[kcol] ; |
---|
2758 | saveU[nstackC]=colUpper[kcol] ; |
---|
2759 | assert (saveU[nstackC]>saveL[nstackC]) ; |
---|
2760 | assert (nstackC<nCols) ; |
---|
2761 | nstackC++ ; |
---|
2762 | onList=true ; |
---|
2763 | } |
---|
2764 | if (newUpper<colsol[kcol]) { |
---|
2765 | if (djs[kcol]>0.0) { |
---|
2766 | /* should be infeasible */ |
---|
2767 | assert (colLower[kcol]>newUpper+primalTolerance_) ; |
---|
2768 | } else { |
---|
2769 | objVal += moveDown*djs[kcol] ; |
---|
2770 | } |
---|
2771 | } |
---|
2772 | if (intVar[kcol]) |
---|
2773 | newUpper = CoinMin(colUpper[kcol],floor(newUpper+1.0e-4)) ; |
---|
2774 | colUpper[kcol]=newUpper ; |
---|
2775 | columnGap[kcol] = newUpper-colLower[kcol]-primalTolerance_ ; |
---|
2776 | if (CoinAbs(colUpper[kcol]-colLower[kcol])<1.0e-6) { |
---|
2777 | markC[kcol]= (tightenLower|tightenUpper); // say fixed |
---|
2778 | } |
---|
2779 | markC[kcol] &= ~(infLower|infUpper) ; |
---|
2780 | if (colUpper[kcol]>1.0e10) |
---|
2781 | markC[kcol] |= infUpper ; |
---|
2782 | if (colLower[kcol]<-1.0e10) |
---|
2783 | markC[kcol] |= infLower ; |
---|
2784 | |
---|
2785 | if (!updateRowBounds(kcol,moveDown, |
---|
2786 | columnStart,columnLength,row,columnElements, |
---|
2787 | rowUpper,rowLower,nstackR,stackR,markR, |
---|
2788 | minR,maxR,saveMin,saveMax)) { |
---|
2789 | colUpper[kcol] = -1.0e50 ; |
---|
2790 | } |
---|
2791 | } |
---|
2792 | if (colLower[kcol]>colUpper[kcol]+primalTolerance_) { |
---|
2793 | notFeasible = true ; |
---|
2794 | k=columnStart[jcol]+columnLength[jcol] ; |
---|
2795 | istackC=nstackC+1 ; |
---|
2796 | break ; |
---|
2797 | } |
---|
2798 | } |
---|
2799 | } |
---|
2800 | } // end big loop rPos->rEnd |
---|
2801 | } else if (doRowLoP) { |
---|
2802 | // Start pos values loop |
---|
2803 | for (int kk =rStart;kk<rEnd;kk++) { |
---|
2804 | int kcol =column[kk] ; |
---|
2805 | if ((markC[kcol]&(tightenLower|tightenUpper))!= (tightenLower|tightenUpper)) { |
---|
2806 | double value2=rowElements[kk] ; |
---|
2807 | int markIt=markC[kcol] ; |
---|
2808 | assert (value2 > 0.0) ; |
---|
2809 | /* positive element */ |
---|
2810 | double gap = columnGap[kcol]*value2 ; |
---|
2811 | bool doDown = (rowLo2 +gap > 0.0) ; |
---|
2812 | if (doDown&&(markIt&(tightenLower|infUpper))==0) { |
---|
2813 | double newLower=1.0 ; |
---|
2814 | double dbound = colUpper[kcol] + rowLo2/value2 ; |
---|
2815 | if (intVar[kcol]) { |
---|
2816 | markIt |= tightenLower ; |
---|
2817 | newLower = ceil(dbound-primalTolerance_) ; |
---|
2818 | } else { |
---|
2819 | newLower=dbound ; |
---|
2820 | if (newLower+primalTolerance_>colUpper[kcol]&& |
---|
2821 | newLower-primalTolerance_<=colUpper[kcol]) { |
---|
2822 | newLower=colUpper[kcol] ; |
---|
2823 | markIt |= tightenLower ; |
---|
2824 | //markIt= (tightenLower|tightenUpper) ; |
---|
2825 | } else { |
---|
2826 | // avoid problems - fix later ? |
---|
2827 | markIt= (tightenLower|tightenUpper) ; |
---|
2828 | } |
---|
2829 | } |
---|
2830 | double moveUp = newLower-colLower[kcol] ; |
---|
2831 | if (!moveUp) |
---|
2832 | continue ; |
---|
2833 | bool onList = ((markC[kcol]&(tightenLower|tightenUpper))!=0) ; |
---|
2834 | if (nstackC<2*maxStack) { |
---|
2835 | markC[kcol] = markIt ; |
---|
2836 | } |
---|
2837 | if (nstackC<2*maxStack) { |
---|
2838 | if (!onList) { |
---|
2839 | stackC[nstackC]=kcol ; |
---|
2840 | saveL[nstackC]=colLower[kcol] ; |
---|
2841 | saveU[nstackC]=colUpper[kcol] ; |
---|
2842 | assert (saveU[nstackC]>saveL[nstackC]) ; |
---|
2843 | assert (nstackC<nCols) ; |
---|
2844 | nstackC++ ; |
---|
2845 | onList=true ; |
---|
2846 | } |
---|
2847 | if (newLower>colsol[kcol]) { |
---|
2848 | if (djs[kcol]<0.0) { |
---|
2849 | /* should be infeasible */ |
---|
2850 | assert (newLower>colUpper[kcol]+primalTolerance_) ; |
---|
2851 | } else { |
---|
2852 | objVal += moveUp*djs[kcol] ; |
---|
2853 | } |
---|
2854 | } |
---|
2855 | if (intVar[kcol]) |
---|
2856 | newLower = CoinMax(colLower[kcol],ceil(newLower-1.0e-4)) ; |
---|
2857 | colLower[kcol]=newLower ; |
---|
2858 | columnGap[kcol] = colUpper[kcol]-newLower-primalTolerance_ ; |
---|
2859 | if (CoinAbs(colUpper[kcol]-colLower[kcol])<1.0e-6) { |
---|
2860 | markC[kcol]= (tightenLower|tightenUpper); // say fixed |
---|
2861 | } |
---|
2862 | markC[kcol] &= ~(infLower|infUpper) ; |
---|
2863 | if (colUpper[kcol]>1.0e10) |
---|
2864 | markC[kcol] |= infUpper ; |
---|
2865 | if (colLower[kcol]<-1.0e10) |
---|
2866 | markC[kcol] |= infLower ; |
---|
2867 | |
---|
2868 | if (!updateRowBounds(kcol,moveUp, |
---|
2869 | columnStart,columnLength,row,columnElements, |
---|
2870 | rowUpper,rowLower,nstackR,stackR,markR, |
---|
2871 | minR,maxR,saveMin,saveMax)) { |
---|
2872 | colLower[kcol] = 1.0e50 ; |
---|
2873 | } |
---|
2874 | } |
---|
2875 | if (colLower[kcol]>colUpper[kcol]+primalTolerance_) { |
---|
2876 | notFeasible = true ; |
---|
2877 | k=columnStart[jcol]+columnLength[jcol] ; |
---|
2878 | istackC=nstackC+1 ; |
---|
2879 | break ; |
---|
2880 | } |
---|
2881 | } |
---|
2882 | } |
---|
2883 | } // end big loop rPos->rEnd |
---|
2884 | } // end processing of positive coefficients of row. |
---|
2885 | } // end loop to walk the column of a variable popped off stackC |
---|
2886 | istackC++ ; |
---|
2887 | } // end stackC processing loop |
---|
2888 | /* |
---|
2889 | PROBE LOOP: END PROPAGATION |
---|
2890 | |
---|
2891 | End propagation of probe in one direction for a single variable. Hard to |
---|
2892 | believe 1,000 lines of code can achieve so little. |
---|
2893 | */ |
---|
2894 | /* |
---|
2895 | PROBE LOOP: INFEASIBILITY |
---|
2896 | |
---|
2897 | Feasibility check. Primal infeasibility is noted in notFeasible; we still |
---|
2898 | need to test against the objective bound. If we have shown up and down |
---|
2899 | infeasibility, we're infeasible, period. Break from the probe loop and |
---|
2900 | terminate the lookedAt and pass loops by forcing their iteration counts past |
---|
2901 | the maximum. |
---|
2902 | |
---|
2903 | If we're not feasible, then we surely can't drive the variable to bound, |
---|
2904 | so reset goingToTrueBound. |
---|
2905 | |
---|
2906 | TODO: I'm coming to think forcing goingToTrueBound = 0 is really an overload |
---|
2907 | to suppress various post-probe activities (cut generation, singleton |
---|
2908 | motion, etc) that should not be performed when the probe shows |
---|
2909 | infeasible. It might be that this is not the best approach. |
---|
2910 | -- lh, 101216 -- |
---|
2911 | */ |
---|
2912 | if (notFeasible || (objVal > cutoff)) { |
---|
2913 | # if CGL_DEBUG > 1 |
---|
2914 | std::cout << " Column " << j << " infeasible: " ; |
---|
2915 | if (notFeasible && (objVal > cutoff)) |
---|
2916 | std::cout << "primal bounds and objective cutoff" ; |
---|
2917 | else if (notFeasible) |
---|
2918 | std::cout << "primal bounds" ; |
---|
2919 | else |
---|
2920 | std::cout << "objective cutoff" ; |
---|
2921 | std::cout << "." << std::endl ; |
---|
2922 | # endif |
---|
2923 | notFeasible = true ; |
---|
2924 | if (iway == upProbe && feasRecord == 0) { |
---|
2925 | ninfeas = 1 ; |
---|
2926 | j = nCols-1 ; |
---|
2927 | iLook = numberThisTime_ ; |
---|
2928 | ipass = maxPass ; |
---|
2929 | break ; |
---|
2930 | } |
---|
2931 | goingToTrueBound = 0 ; |
---|
2932 | } else { |
---|
2933 | feasRecord |= feasValue[iway] ; |
---|
2934 | } |
---|
2935 | /* |
---|
2936 | Save the implications? Currently restricted to binary variables. If info |
---|
2937 | were a CglTreeProbing object, fixes() would save the implication in the |
---|
2938 | arrays kept there, returning false only when space is exceeded. |
---|
2939 | |
---|
2940 | TODO: This would be one of the places that will fail if fixes() and |
---|
2941 | initializeFixing() are removed from CglTreeInfo. -- lh, 101127 -- |
---|
2942 | */ |
---|
2943 | if (!notFeasible && saveFixingInfo) { |
---|
2944 | // save fixing info |
---|
2945 | assert (j == stackC[0]) ; |
---|
2946 | int toValue = (way[iway] == probeDown) ? -1 : +1 ; |
---|
2947 | for (istackC = 1 ; istackC < nstackC ; istackC++) { |
---|
2948 | int icol = stackC[istackC] ; |
---|
2949 | // for now back to just 0-1 |
---|
2950 | if (colUpper[icol]-colLower[icol]<1.0e-12 && |
---|
2951 | !saveL[istackC] && saveU[istackC]==1.0) { |
---|
2952 | assert(saveL[istackC] == colLower[icol] || |
---|
2953 | saveU[istackC] == colUpper[icol]) ; |
---|
2954 | saveFixingInfo = |
---|
2955 | info->fixes(j,toValue,icol, |
---|
2956 | (colLower[icol] == saveL[istackC])) ; |
---|
2957 | } |
---|
2958 | } |
---|
2959 | } |
---|
2960 | /* |
---|
2961 | PROBE LOOP: MONOTONE (BREAK AND KEEP) |
---|
2962 | |
---|
2963 | We can't prove infeasibility. What's left? If we're at the end of the up |
---|
2964 | probe, we may know enough to prove monotonicity: |
---|
2965 | |
---|
2966 | * If this is the up probe and we were feasible on this probe but not on |
---|
2967 | the down pass, we're monotone. |
---|
2968 | * If this is the second down probe, we were feasible on down, infeasible on |
---|
2969 | up (monotone) and repeated the down probe to recreate the bounds. |
---|
2970 | |
---|
2971 | Either way, execute monotoneActions to capture the state and move on to the |
---|
2972 | next variable. Terminate the probe loop by forcing iway = oneProbeTooMany. |
---|
2973 | |
---|
2974 | This if-then-else extends to the end of the down/up/down probe loop. |
---|
2975 | |
---|
2976 | TODO: There's something not right about this `second down pass' business. |
---|
2977 | Why do we copy stackC to stackC0 on a feasible first pass? Surely that |
---|
2978 | could be used to restore the state of the first down pass. |
---|
2979 | -- lh, 101128 -- |
---|
2980 | */ |
---|
2981 | if (iway == downProbe2 || |
---|
2982 | (iway == upProbe && feasRecord == upProbeFeas)) { |
---|
2983 | nfixed++ ; |
---|
2984 | bool retVal = monotoneActions(primalTolerance_,si,cs, |
---|
2985 | nstackC,stackC,intVar, |
---|
2986 | colLower,colsol,colUpper, |
---|
2987 | index,element) ; |
---|
2988 | if (retVal) anyColumnCuts = true ; |
---|
2989 | clearStacks(primalTolerance_,nstackC,stackC,markC,colLower,colUpper, |
---|
2990 | nstackR,stackR,markR) ; |
---|
2991 | break ; |
---|
2992 | } |
---|
2993 | /* |
---|
2994 | PROBE LOOP: ITERATE |
---|
2995 | |
---|
2996 | If we reach here, we may iterate the probe loop. |
---|
2997 | |
---|
2998 | Cases remaining include |
---|
2999 | |
---|
3000 | a) End of the first down probe; in this case we will iterate regardless of |
---|
3001 | feasibility. |
---|
3002 | b) End of the up probe, up and down feasible; in this case we will not |
---|
3003 | iterate (we're done, cannot prove monotonicity). |
---|
3004 | c) End of the up probe, down feasible, up infeasible; in this case we will |
---|
3005 | iterate to restore the down feasible state. |
---|
3006 | |
---|
3007 | TODO: Unless I miss my guess, large chunks of this block of code will be |
---|
3008 | replicated in the code that handles case b), as we generate cuts |
---|
3009 | related to forcing the variable up. Further inspection confirms |
---|
3010 | this. -- lh, 101128 -- |
---|
3011 | |
---|
3012 | TODO: I'm becoming more convinced that the second down probe iteration is |
---|
3013 | unnecessary. Given that we have stackC0, lo0, and up0, seems like |
---|
3014 | calling monotoneActions with lo0 = colLower, up0 = colUpper, and |
---|
3015 | stackC0 = stackC should do the trick. -- lh, 101216 -- |
---|
3016 | |
---|
3017 | TODO: As I get the control flow sorted out a bit better, I should try to |
---|
3018 | separate b) (which does not iterate) from a) and c). |
---|
3019 | |
---|
3020 | The next block deals with the end of the first down probe. If the probe |
---|
3021 | is feasible, attempt to tighten column bounds with singletonRows, process |
---|
3022 | stackC to generate disaggregation cuts, and copy the tightened bounds |
---|
3023 | to stackC0, lo0, and up0 for use after the up probe. Feasible or not, |
---|
3024 | we'll need to go on to the up probe, so restore the original bounds. |
---|
3025 | |
---|
3026 | Note that GoingToTrueBound was set to 0 way back in initial setup unless |
---|
3027 | the user requested disaggregation cuts (rowCuts&0x01). |
---|
3028 | |
---|
3029 | I've rearranged the control flow here, moving the call to singletonRows |
---|
3030 | ahead of the code that captures column bound changes to stackC0, et al. |
---|
3031 | This makes it possible to use the same code at the end of case b) (down and |
---|
3032 | up probe both feasible). Also, singletonRows is just another way to tighten |
---|
3033 | column bounds. In its original position, it was subject to the same |
---|
3034 | restrictions as coefficient strengthening (coefficient strengthening enabled |
---|
3035 | and goingToTrueBound for the probing variable). I don't see any reason to |
---|
3036 | keep them (but it might be good to add a separate control bit precisely for |
---|
3037 | singletonRows). -- lh, 110113 -- |
---|
3038 | */ |
---|
3039 | if (iway == downProbe) { |
---|
3040 | if (!notFeasible) { |
---|
3041 | /* |
---|
3042 | Attempt to tighten singleton bounds and generate disaggregation cuts. |
---|
3043 | */ |
---|
3044 | singletonRows(j,primalTolerance_,si,rowCopy,markC, |
---|
3045 | nstackC,stackC,saveL,saveU, |
---|
3046 | colUpper,colLower,columnGap, |
---|
3047 | nstackR,stackR,rowUpper,rowLower,maxR,minR) ; |
---|
3048 | if (goingToTrueBound == 2 && !justReplace) { |
---|
3049 | disaggCuts(nstackC,probeDown,primalTolerance_, |
---|
3050 | disaggEffectiveness,si,rowCut,stackC,colsol, |
---|
3051 | colUpper,colLower,saveU,saveL,index,element) ; |
---|
3052 | } |
---|
3053 | /* |
---|
3054 | Copy off the bound changes from the down probe -- we'll need them after the |
---|
3055 | up probe. |
---|
3056 | */ |
---|
3057 | nstackC0 = CoinMin(nstackC,maxStack) ; |
---|
3058 | for (istackC = 0 ; istackC < nstackC0 ; istackC++) { |
---|
3059 | int icol = stackC[istackC] ; |
---|
3060 | stackC0[istackC] = icol ; |
---|
3061 | lo0[istackC] = colLower[icol] ; |
---|
3062 | up0[istackC] = colUpper[icol] ; |
---|
3063 | } |
---|
3064 | } else { |
---|
3065 | nstackC0 = 0 ; |
---|
3066 | } |
---|
3067 | /* |
---|
3068 | Now reset column bounds to their original values in preparation for the |
---|
3069 | up probe. |
---|
3070 | */ |
---|
3071 | for (istackC = nstackC-1 ; istackC >= 0 ; istackC--) { |
---|
3072 | int icol = stackC[istackC] ; |
---|
3073 | double oldU = saveU[istackC] ; |
---|
3074 | double oldL = saveL[istackC] ; |
---|
3075 | colUpper[icol] = oldU ; |
---|
3076 | colLower[icol] = oldL ; |
---|
3077 | columnGap[icol] = oldU-oldL-primalTolerance_ ; |
---|
3078 | markC[icol] = 0 ; |
---|
3079 | if (oldU > 1.0e10) |
---|
3080 | markC[icol] |= infUpper ; |
---|
3081 | if (oldL < -1.0e10) |
---|
3082 | markC[icol] |= infLower ; |
---|
3083 | } |
---|
3084 | /* |
---|
3085 | And now for the rows. Remember, we're still finishing off the first down |
---|
3086 | probe and readying for the next iteration (up probe). The final three lines |
---|
3087 | (of about 300) actually restore the row bounds. The rest is cut generation. |
---|
3088 | */ |
---|
3089 | if ((rowCuts&0x02) != 0 && goingToTrueBound && !anyColumnCuts) |
---|
3090 | strengthenCoeff(j,iway,primalTolerance_,justReplace,canReplace, |
---|
3091 | needEffectiveness,si,rowCut, |
---|
3092 | rowCopy,colUpper,colLower,colsol,nstackR,stackR, |
---|
3093 | rowUpper,rowLower,maxR,minR,realRows, |
---|
3094 | element,index,info) ; |
---|
3095 | /* |
---|
3096 | Restore row bounds. |
---|
3097 | */ |
---|
3098 | for (istackR = 0 ; istackR < nstackR ; istackR++) { |
---|
3099 | int irow = stackR[istackR] ; |
---|
3100 | minR[irow] = saveMin[istackR] ; |
---|
3101 | maxR[irow] = saveMax[istackR] ; |
---|
3102 | markR[irow] = -1 ; |
---|
3103 | } |
---|
3104 | } // end of code to iterate after first down pass |
---|
3105 | /* |
---|
3106 | PROBE LOOP: DOWN AND UP FEASIBLE |
---|
3107 | |
---|
3108 | Again, not quite the full story. Cases remaining: |
---|
3109 | b) End of up probe, up and down feasible; in this case we will not |
---|
3110 | iterate. |
---|
3111 | c) End of up probe, down feasible, up infeasible; in this case we will |
---|
3112 | iterate to restore the down feasible state. |
---|
3113 | |
---|
3114 | The next block deals with case b). We don't want to iterate the |
---|
3115 | down/up/down loop, so force iway to oneProbeTooMany. As usual, the move |
---|
3116 | singleton code consists of two nearly identical blocks, one working off |
---|
3117 | U(i), the next off L(i). |
---|
3118 | |
---|
3119 | TODO: The conditions guarding the move singleton code here are different than |
---|
3120 | those for pass 0: here there's an added guard using strengthenRow. |
---|
3121 | Just a failure to change all instances? Nor do I see how John's |
---|
3122 | comment (`just at root') is enforced here, unless it's from outside |
---|
3123 | via rowCuts or strengthenRow. |
---|
3124 | |
---|
3125 | The check for any column cuts is not present, so presumably we can |
---|
3126 | do singleton motion in the presence of column cuts. The check for gap |
---|
3127 | does not have the high end (gap < 1.0e8) test. |
---|
3128 | |
---|
3129 | Note also that here the test for cut generation is moved outside the |
---|
3130 | loop that iterates over rows. Presumably that's because in the |
---|
3131 | previous case, the loop also restored markR, etc. |
---|
3132 | |
---|
3133 | Note that the code that handles the bound change is considerably |
---|
3134 | different than the previous instance in pass 0. We're adding the |
---|
3135 | changes to stackC rather than stackC0. |
---|
3136 | -- lh, 101128 -- |
---|
3137 | |
---|
3138 | Moved the call to disaggCuts up from the common code for b) and c), because |
---|
3139 | c) (up infeasible) sets goingToTrueBound to 0 and we can't execute |
---|
3140 | disaggCuts. Similarly for strengthenCoeff. |
---|
3141 | */ |
---|
3142 | else { |
---|
3143 | if (iway == upProbe && |
---|
3144 | feasRecord == (downProbeFeas|upProbeFeas)) { |
---|
3145 | iway = oneProbeTooMany ; |
---|
3146 | |
---|
3147 | singletonRows(j,primalTolerance_,si,rowCopy,markC, |
---|
3148 | nstackC,stackC,saveL,saveU, |
---|
3149 | colUpper,colLower,columnGap, |
---|
3150 | nstackR,stackR,rowUpper,rowLower,maxR,minR) ; |
---|
3151 | if (goingToTrueBound == 2 && !justReplace) { |
---|
3152 | disaggCuts(nstackC,probeUp,primalTolerance_, |
---|
3153 | disaggEffectiveness,si, |
---|
3154 | rowCut,stackC,colsol,colUpper,colLower,saveU,saveL, |
---|
3155 | index,element) ; |
---|
3156 | if ((rowCuts&0x02) != 0 && goingToTrueBound && !anyColumnCuts) |
---|
3157 | strengthenCoeff(j,iway,primalTolerance_,justReplace,canReplace, |
---|
3158 | needEffectiveness,si,rowCut, |
---|
3159 | rowCopy,colUpper,colLower,colsol,nstackR,stackR, |
---|
3160 | rowUpper,rowLower,maxR,minR,realRows, |
---|
3161 | element,index,info) ; |
---|
3162 | } |
---|
3163 | /* |
---|
3164 | jjf: point back to stack |
---|
3165 | |
---|
3166 | We're done playing with singletons. Get to the real work here. We don't need |
---|
3167 | markC any more; coopt it for a cross-reference, var index -> stackC index. |
---|
3168 | The +1000 offset allows us to distinguish invalid entries (not all variables |
---|
3169 | are on stackC). |
---|
3170 | */ |
---|
3171 | for (istackC = nstackC-1 ; istackC >= 0 ; istackC--) { |
---|
3172 | int icol = stackC[istackC] ; |
---|
3173 | markC[icol] = istackC+1000 ; |
---|
3174 | } |
---|
3175 | OsiColCut cc ; |
---|
3176 | int nTot = 0, nFix = 0, nInt = 0 ; |
---|
3177 | bool ifCut = false ; |
---|
3178 | /* |
---|
3179 | See if we have bounds improvement. Given a lower bound ld<j> from the |
---|
3180 | down probe, a bound lu<j> from the up probe, and the original bound lo<j>, |
---|
3181 | the bound we want is max(min(ld<j>,lu<j>),lo<j>). Use a conservative |
---|
3182 | tolerance. |
---|
3183 | */ |
---|
3184 | for (istackC = 1 ; istackC < nstackC0 ; istackC++) { |
---|
3185 | int icol = stackC0[istackC] ; |
---|
3186 | int istackC1 = markC[icol]-1000 ; |
---|
3187 | if (istackC1 >= 0) { |
---|
3188 | if (CoinMin(lo0[istackC],colLower[icol]) > |
---|
3189 | saveL[istackC1]+1.0e-4) { |
---|
3190 | saveL[istackC1] = CoinMin(lo0[istackC],colLower[icol]) ; |
---|
3191 | if (intVar[icol] /*||!info->inTree*/ ) { |
---|
3192 | element[nFix] = saveL[istackC1] ; |
---|
3193 | index[nFix++] = icol ; |
---|
3194 | nInt++ ; |
---|
3195 | if (colsol[icol] < saveL[istackC1]-primalTolerance_) |
---|
3196 | ifCut = true ; |
---|
3197 | } |
---|
3198 | nfixed++ ; |
---|
3199 | } |
---|
3200 | } |
---|
3201 | } |
---|
3202 | if (nFix) { |
---|
3203 | nTot = nFix ; |
---|
3204 | cc.setLbs(nFix,index,element) ; |
---|
3205 | nFix = 0 ; |
---|
3206 | } |
---|
3207 | /* |
---|
3208 | Repeat for upper bounds. There's an odd bit of asymmetry here; this loop |
---|
3209 | tries to generate a cut if there's no bound improvement. |
---|
3210 | |
---|
3211 | TODO: What, pray tell, is a `two cut'? Binary variables only, it seems. |
---|
3212 | Judging from the conditions, we're looking for a variable that's |
---|
3213 | fixed at 0 on a down probe and fixed at 1 on an up probe, or vice |
---|
3214 | versa. |
---|
3215 | -- lh, 101128 -- |
---|
3216 | */ |
---|
3217 | for (istackC = 1 ; istackC < nstackC0 ; istackC++) { |
---|
3218 | int icol = stackC0[istackC] ; |
---|
3219 | int istackC1 = markC[icol]-1000 ; |
---|
3220 | if (istackC1 >= 0) { |
---|
3221 | if (CoinMax(up0[istackC],colUpper[icol]) < |
---|
3222 | saveU[istackC1]-1.0e-4) { |
---|
3223 | saveU[istackC1] = CoinMax(up0[istackC],colUpper[icol]) ; |
---|
3224 | if (intVar[icol] /*||!info->inTree*/ ) { |
---|
3225 | element[nFix] = saveU[istackC1] ; |
---|
3226 | index[nFix++] = icol ; |
---|
3227 | nInt++ ; |
---|
3228 | if (colsol[icol] > saveU[istackC1]+primalTolerance_) |
---|
3229 | ifCut = true ; |
---|
3230 | } |
---|
3231 | nfixed++ ; |
---|
3232 | } else if (!info->inTree && |
---|
3233 | saveL[0] == 0.0 && saveU[0] == 1.0) { |
---|
3234 | // See if can do two cut |
---|
3235 | double upperWhenDown = up0[istackC] ; |
---|
3236 | double lowerWhenDown = lo0[istackC] ; |
---|
3237 | double upperWhenUp = colUpper[icol] ; |
---|
3238 | double lowerWhenUp = colLower[icol] ; |
---|
3239 | double upperOriginal = saveU[istackC1] ; |
---|
3240 | double lowerOriginal = saveL[istackC1] ; |
---|
3241 | if (upperWhenDown < lowerOriginal+1.0e-12 && |
---|
3242 | lowerWhenUp > upperOriginal-1.0e-12) { |
---|
3243 | OsiRowCut rc ; |
---|
3244 | rc.setLb(lowerOriginal) ; |
---|
3245 | rc.setUb(lowerOriginal) ; |
---|
3246 | rc.setEffectiveness(1.0e-5) ; |
---|
3247 | int index[2] ; |
---|
3248 | double element[2] ; |
---|
3249 | index[0] = j ; |
---|
3250 | index[1] = icol ; |
---|
3251 | element[0] = -(upperOriginal-lowerOriginal) ; |
---|
3252 | /* |
---|
3253 | jjf: If zero then - must have been fixed without noticing! |
---|
3254 | |
---|
3255 | TODO: Uh, fixed without noticing?! That's an algorithm problem. |
---|
3256 | -- lh, 101128 -- |
---|
3257 | */ |
---|
3258 | if (CoinAbs(element[0]) > 1.0e-8) { |
---|
3259 | element[1] = 1.0 ; |
---|
3260 | rc.setRow(2,index,element,false) ; |
---|
3261 | cs.insert(rc) ; |
---|
3262 | } |
---|
3263 | } else if (upperWhenUp < lowerOriginal+1.0e-12 && |
---|
3264 | lowerWhenDown > upperOriginal-1.0e-12) { |
---|
3265 | OsiRowCut rc ; |
---|
3266 | rc.setLb(upperOriginal) ; |
---|
3267 | rc.setUb(upperOriginal) ; |
---|
3268 | rc.setEffectiveness(1.0e-5) ; |
---|
3269 | int index[2] ; |
---|
3270 | double element[2] ; |
---|
3271 | index[0] = j ; |
---|
3272 | index[1] = icol ; |
---|
3273 | element[0] = upperOriginal-lowerOriginal ; |
---|
3274 | element[1] = 1.0 ; |
---|
3275 | rc.setRow(2,index,element,false) ; |
---|
3276 | cs.insert(rc) ; |
---|
3277 | } |
---|
3278 | } |
---|
3279 | } |
---|
3280 | } // end loop to update upper bounds (w/ two cuts) |
---|
3281 | if (nFix) { |
---|
3282 | nTot+=nFix ; |
---|
3283 | cc.setUbs(nFix,index,element) ; |
---|
3284 | } |
---|
3285 | // could tighten continuous as well |
---|
3286 | if (nInt) { |
---|
3287 | if (ifCut) { |
---|
3288 | cc.setEffectiveness(100.0) ; |
---|
3289 | } else { |
---|
3290 | cc.setEffectiveness(1.0e-5) ; |
---|
3291 | } |
---|
3292 | #if CGL_DEBUG > 0 |
---|
3293 | checkBounds(si,cc) ; |
---|
3294 | #endif |
---|
3295 | cs.insert(cc) ; |
---|
3296 | } |
---|
3297 | } // end of code to handle up and down feasible. |
---|
3298 | /* |
---|
3299 | PROBE LOOP: DOWN FEASIBLE, UP INFEASIBLE |
---|
3300 | |
---|
3301 | The only remaining case is down feasible, up infeasible, in which case we |
---|
3302 | need to reset to the initial state and run the down probe again. Surely |
---|
3303 | there must be a better way? |
---|
3304 | */ |
---|
3305 | else { |
---|
3306 | goingToTrueBound = 0 ; |
---|
3307 | } |
---|
3308 | /* |
---|
3309 | PROBE LOOP: RESTORE |
---|
3310 | |
---|
3311 | And now the code to reset to the initial state. |
---|
3312 | */ |
---|
3313 | /* restore all */ |
---|
3314 | for (istackC = nstackC-1 ; istackC >= 0 ; istackC--) { |
---|
3315 | int icol = stackC[istackC] ; |
---|
3316 | double oldU = saveU[istackC] ; |
---|
3317 | double oldL = saveL[istackC] ; |
---|
3318 | /* |
---|
3319 | TODO: This next bit differs in subtle ways from the restore at the end |
---|
3320 | of the down probe. Here we force markC to show fixed if the bounds are |
---|
3321 | within 1.0e-4 of one another, a fairly broad tolerance; code for |
---|
3322 | the pass 0 restore does not restore fixed status. Also, |
---|
3323 | we're not distinguishing here between actions for down feasible, |
---|
3324 | up feasible, vs. down feasible, up infeasible. -- lh, 101128 -- |
---|
3325 | |
---|
3326 | I still don't see any reason for the restore here to differ from |
---|
3327 | the restore after the down probe. -- lh, 110113 -- |
---|
3328 | */ |
---|
3329 | colUpper[icol] = oldU ; |
---|
3330 | colLower[icol] = oldL ; |
---|
3331 | columnGap[icol] = oldU-oldL-primalTolerance_ ; |
---|
3332 | if (oldU > oldL+1.0e-4) { |
---|
3333 | markC[icol] = 0 ; |
---|
3334 | if (oldU > 1.0e10) |
---|
3335 | markC[icol] |= infUpper ; |
---|
3336 | if (oldL < -1.0e10) |
---|
3337 | markC[icol] |= infLower ; |
---|
3338 | } else { |
---|
3339 | markC[icol] = (tightenUpper|tightenLower) ; |
---|
3340 | } |
---|
3341 | } |
---|
3342 | /* |
---|
3343 | End of column restore. On to rows. |
---|
3344 | */ |
---|
3345 | for (istackR = 0 ; istackR < nstackR ; istackR++) { |
---|
3346 | int irow = stackR[istackR] ; |
---|
3347 | minR[irow] = saveMin[istackR] ; |
---|
3348 | maxR[irow] = saveMax[istackR] ; |
---|
3349 | markR[irow] = -1 ; |
---|
3350 | } |
---|
3351 | } // end of PROBE: DOWN AND UP FEASIBLE |
---|
3352 | } // PROBE LOOP: END |
---|
3353 | } // LOOKEDAT LOOP: END |
---|
3354 | } // PASS LOOP: END |
---|
3355 | /* |
---|
3356 | Free up the space we've been using. |
---|
3357 | */ |
---|
3358 | #ifndef ONE_ARRAY |
---|
3359 | delete [] stackC0 ; |
---|
3360 | delete [] lo0 ; |
---|
3361 | delete [] up0 ; |
---|
3362 | delete [] columnGap ; |
---|
3363 | delete [] markC ; |
---|
3364 | delete [] stackC ; |
---|
3365 | delete [] stackR ; |
---|
3366 | delete [] saveL ; |
---|
3367 | delete [] saveU ; |
---|
3368 | delete [] saveMin ; |
---|
3369 | delete [] saveMax ; |
---|
3370 | delete [] index ; |
---|
3371 | delete [] element ; |
---|
3372 | delete [] djs ; |
---|
3373 | delete [] largestPositiveInRow ; |
---|
3374 | delete [] largestNegativeInRow ; |
---|
3375 | #endif |
---|
3376 | delete [] colsol ; |
---|
3377 | /* |
---|
3378 | jjf: Add in row cuts |
---|
3379 | |
---|
3380 | Transfer row cuts from the local container into something that can escape |
---|
3381 | this method. If we're not in `just replace' mode, simply copy them over to |
---|
3382 | the container (cs) passed as a parameter using addCuts. If we're in `just |
---|
3383 | replace' mode, ? need to explore further ? -- lh, 101125 -- |
---|
3384 | |
---|
3385 | TODO: If there's any possibility of realRows[i] < 0, this code will break |
---|
3386 | trying to read strengthenRow[]! |
---|
3387 | -- lh, 101128 -- |
---|
3388 | */ |
---|
3389 | if (!ninfeas) { |
---|
3390 | if (!justReplace) { |
---|
3391 | rowCut.addCuts(cs,info->strengthenRow,info->pass) ; |
---|
3392 | } else { |
---|
3393 | for (int i = 0 ; i < nRows ; i++) { |
---|
3394 | int realRow = realRows[i] ; |
---|
3395 | OsiRowCut *cut = info->strengthenRow[realRow] ; |
---|
3396 | if (realRow >= 0 && cut) { |
---|
3397 | #ifdef CLP_INVESTIGATE |
---|
3398 | printf("Row %d, real row %d effectiveness %g\n", |
---|
3399 | i,realRow,cut->effectiveness()) ; |
---|
3400 | #endif |
---|
3401 | cs.insert(cut) ; |
---|
3402 | } |
---|
3403 | } |
---|
3404 | } |
---|
3405 | } |
---|
3406 | #if 0 |
---|
3407 | /* |
---|
3408 | Debug print. |
---|
3409 | */ |
---|
3410 | { |
---|
3411 | int numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
3412 | int k ; |
---|
3413 | for (k = 0;k<numberRowCutsAfter;k++) { |
---|
3414 | OsiRowCut thisCut = cs.rowCut(k) ; |
---|
3415 | printf("Cut %d is %g <=",k,thisCut.lb()) ; |
---|
3416 | int n=thisCut.row().getNumElements() ; |
---|
3417 | const int * column = thisCut.row().getIndices() ; |
---|
3418 | const double * element = thisCut.row().getElements() ; |
---|
3419 | assert (n) ; |
---|
3420 | for (int i=0;i<n;i++) { |
---|
3421 | printf(" %g*x%d",element[i],column[i]) ; |
---|
3422 | } |
---|
3423 | printf(" <= %g\n",thisCut.ub()) ; |
---|
3424 | } |
---|
3425 | } |
---|
3426 | #endif |
---|
3427 | |
---|
3428 | # if CGL_DEBUG > 0 |
---|
3429 | std::cout |
---|
3430 | << "Leaving CglProbing::probe, ninfeas " << ninfeas |
---|
3431 | << "." << std::endl ; |
---|
3432 | # endif |
---|
3433 | |
---|
3434 | return (ninfeas) ; |
---|
3435 | } |
---|