1 // Last edit: 02/05/2013
2 //
3 // Name: CglGMI.cpp
4 // Author: Giacomo Nannicini
5 // Singapore University of Technology and Design, Singapore
6 // email: nannicini@sutd.edu.sg
7 // Date: 11/17/09
8 //-----------------------------------------------------------------------------
9 // Copyright (C) 2009, Giacomo Nannicini. All Rights Reserved.
10
11 #include <cstdlib>
12 #include <cstdio>
13 #include <cmath>
14 #include <cfloat>
15 #include <cassert>
16 #include <iostream>
17 #include <climits>
18
19 #include "CoinPragma.hpp"
20 #include "CoinHelperFunctions.hpp"
21 #include "CoinPackedVector.hpp"
22 #include "CoinPackedMatrix.hpp"
23 #include "CoinIndexedVector.hpp"
24 #include "OsiSolverInterface.hpp"
25 #include "OsiRowCutDebugger.hpp"
26 #include "CoinFactorization.hpp"
27 #include "CglGMI.hpp"
28 #include "CoinFinite.hpp"
29 #include "CoinRational.hpp"
30
31 //-------------------------------------------------------------------
32 // Generate GMI cuts
33 //-------------------------------------------------------------------
34
35 /***************************************************************************/
CglGMI()36 CglGMI::CglGMI() :
37 CglCutGenerator(),
38 param(),
39 nrow(0),
40 ncol(0),
41 colLower(NULL),
42 colUpper(NULL),
43 rowLower(NULL),
44 rowUpper(NULL),
45 rowRhs(NULL),
46 isInteger(NULL),
47 cstat(NULL),
48 rstat(NULL),
49 solver(NULL),
50 xlp(NULL),
51 rowActivity(NULL),
52 byRow(NULL),
53 byCol(NULL),
54 f0(0.0),
55 f0compl(0.0),
56 ratiof0compl(0.0)
57 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
58 ,
59 trackRejection(false),
60 fracFail(0),
61 dynFail(0),
62 violFail(0),
63 suppFail(0),
64 scaleFail(0),
65 numGeneratedCuts(0)
66 #endif
67 {
68
69 }
70
71 /***************************************************************************/
CglGMI(const CglGMIParam & parameters)72 CglGMI::CglGMI(const CglGMIParam ¶meters) :
73 CglCutGenerator(),
74 param(parameters),
75 nrow(0),
76 ncol(0),
77 colLower(NULL),
78 colUpper(NULL),
79 rowLower(NULL),
80 rowUpper(NULL),
81 rowRhs(NULL),
82 isInteger(NULL),
83 cstat(NULL),
84 rstat(NULL),
85 solver(NULL),
86 xlp(NULL),
87 rowActivity(NULL),
88 byRow(NULL),
89 byCol(NULL),
90 f0(0.0),
91 f0compl(0.0),
92 ratiof0compl(0.0)
93 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
94 ,
95 trackRejection(false),
96 fracFail(0),
97 dynFail(0),
98 violFail(0),
99 suppFail(0),
100 scaleFail(0),
101 numGeneratedCuts(0)
102 #endif
103 {
104
105 }
106
107 /***************************************************************************/
CglGMI(const CglGMI & rhs)108 CglGMI::CglGMI(const CglGMI& rhs) :
109 CglCutGenerator(rhs),
110 param(rhs.param),
111 nrow(rhs.nrow),
112 ncol(rhs.ncol),
113 colLower(rhs.colLower),
114 colUpper(rhs.colUpper),
115 rowLower(rhs.rowLower),
116 rowUpper(rhs.rowUpper),
117 rowRhs(rhs.rowRhs),
118 isInteger(rhs.isInteger),
119 cstat(rhs.cstat),
120 rstat(rhs.rstat),
121 solver(rhs.solver),
122 xlp(rhs.xlp),
123 rowActivity(rhs.rowActivity),
124 byRow(rhs.byRow),
125 byCol(rhs.byCol),
126 f0(rhs.f0),
127 f0compl(rhs.f0compl),
128 ratiof0compl(rhs.ratiof0compl)
129 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
130 ,
131 trackRejection(rhs.trackRejection),
132 fracFail(rhs.fracFail),
133 dynFail(rhs.dynFail),
134 violFail(rhs.violFail),
135 suppFail(rhs.suppFail),
136 scaleFail(rhs.scaleFail),
137 numGeneratedCuts(rhs.numGeneratedCuts)
138 #endif
139 {
140
141 }
142
143 /***************************************************************************/
operator =(const CglGMI & rhs)144 CglGMI & CglGMI::operator=(const CglGMI& rhs) {
145 if(this != &rhs){
146 CglCutGenerator::operator=(rhs);
147 param = rhs.param;
148 nrow = rhs.nrow;
149 ncol = rhs.ncol;
150 colLower = rhs.colLower;
151 colUpper = rhs.colUpper;
152 rowLower = rhs.rowLower;
153 rowUpper = rhs.rowUpper;
154 rowRhs = rhs.rowRhs;
155 isInteger = rhs.isInteger;
156 cstat = rhs.cstat;
157 rstat = rhs.rstat;
158 solver = rhs.solver;
159 xlp = rhs.xlp;
160 rowActivity = rhs.rowActivity;
161 byRow = rhs.byRow;
162 byCol = rhs.byCol;
163 f0 = rhs.f0;
164 f0compl = rhs.f0compl;
165 ratiof0compl = rhs.ratiof0compl;
166 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
167 trackRejection = rhs.trackRejection;
168 fracFail = rhs.fracFail;
169 dynFail = rhs.dynFail;
170 violFail = rhs.violFail;
171 suppFail = rhs.suppFail;
172 scaleFail = rhs.scaleFail;
173 numGeneratedCuts = rhs.numGeneratedCuts;
174 #endif
175 }
176 return *this;
177 }
178
179 /***************************************************************************/
~CglGMI()180 CglGMI::~CglGMI() {
181
182 }
183
184 /*********************************************************************/
185 CglCutGenerator *
clone() const186 CglGMI::clone() const
187 {
188 return new CglGMI(*this);
189 }
190
191 /***************************************************************************/
192
193 // Returns (value - floor)
aboveInteger(double value) const194 inline double CglGMI::aboveInteger(double value) const {
195 return (value - floor(value));
196 } /* aboveInteger */
197
198 /**********************************************************/
printvecINT(const char * vecstr,const int * x,int n) const199 void CglGMI::printvecINT(const char *vecstr, const int *x, int n) const {
200 int num, fromto, upto;
201
202 num = (n/10) + 1;
203 printf("%s :\n", vecstr);
204 for (int j = 0; j < num; ++j) {
205 fromto = 10*j;
206 upto = 10 * (j+1);
207 if(n <= upto) upto = n;
208 for (int i = fromto; i < upto; ++i)
209 printf(" %4d", x[i]);
210 printf("\n");
211 }
212 printf("\n");
213 } /* printvecINT */
214
215 /**********************************************************/
printvecDBL(const char * vecstr,const double * x,int n) const216 void CglGMI::printvecDBL(const char *vecstr, const double *x, int n) const
217 {
218 int num, fromto, upto;
219
220 num = (n/10) + 1;
221 printf("%s :\n", vecstr);
222 for (int j = 0; j < num; ++j) {
223 fromto = 10*j;
224 upto = 10 * (j+1);
225 if(n <= upto) upto = n;
226 for (int i = fromto; i < upto; ++i)
227 printf(" %7.3f", x[i]);
228 printf("\n");
229 }
230 printf("\n");
231 } /* printvecDBL */
232
233 /**********************************************************/
printvecDBL(const char * vecstr,const double * elem,const int * index,int nz) const234 void CglGMI::printvecDBL(const char *vecstr, const double *elem,
235 const int * index, int nz) const
236 {
237 printf("%s\n", vecstr);
238 int written = 0;
239 for (int j = 0; j < nz; ++j) {
240 written += printf("%d:%.3f ", index[j], elem[j]);
241 if (written > 70) {
242 printf("\n");
243 written = 0;
244 }
245 }
246 if (written > 0) {
247 printf("\n");
248 }
249
250 } /* printvecDBL */
251
252 /************************************************************************/
computeCutFractionality(double varRhs,double & cutRhs)253 inline bool CglGMI::computeCutFractionality(double varRhs,
254 double& cutRhs) {
255 f0 = aboveInteger(varRhs);
256 f0compl = 1 - f0;
257 if (f0 < param.getAway() || f0compl < param.getAway())
258 return false;
259 ratiof0compl = f0/f0compl;
260 cutRhs = -f0;
261 return true;
262 } /* computeCutFractionality */
263
264 /************************************************************************/
computeCutCoefficient(double rowElem,int index)265 inline double CglGMI::computeCutCoefficient(double rowElem, int index) {
266
267 // See Wolsey "Integer Programming" (1998), p. 130, fourth line from top
268 // after correcting typo (Proposition 8.8), flipping all signs to get <=.
269
270 if (index < ncol && isInteger[index]) {
271 double f = aboveInteger(rowElem);
272 if(f > f0) {
273 return (-((1-f) * ratiof0compl));
274 }
275 else {
276 return (-f);
277 }
278 }
279 else{
280 if(rowElem < 0) {
281 return (rowElem*ratiof0compl);
282 }
283 else {
284 return (-rowElem);
285 }
286 }
287
288 } /* computeCutCoefficient */
289
290 /************************************************************************/
eliminateSlack(double cutElem,int index,double * cut,double & cutRhs,const double * elements,const CoinBigIndex * rowStart,const int * indices,const int * rowLength,const double * rhs)291 inline void CglGMI::eliminateSlack(double cutElem, int index, double* cut,
292 double& cutRhs, const double *elements,
293 const CoinBigIndex *rowStart, const int *indices,
294 const int *rowLength, const double *rhs) {
295
296 // now i is where coefficients on slack variables begin;
297 // eliminate the slacks
298 int rowpos = index - ncol;
299 if(fabs(cutElem) > param.getEPS_ELIM()) {
300 if (areEqual(rowLower[rowpos], rowUpper[rowpos],
301 param.getEPS(), param.getEPS())) {
302 // "almost" fixed slack, we'll just skip it
303 return;
304 }
305
306 CoinBigIndex upto = rowStart[rowpos] + rowLength[rowpos];
307 for (CoinBigIndex j = rowStart[rowpos]; j < upto; ++j) {
308 cut[indices[j]] -= cutElem * elements[j];
309 }
310 cutRhs -= cutElem * rhs[rowpos];
311 }
312
313 } /* eliminateSlack */
314
315 /************************************************************************/
flip(double & rowElem,int index)316 inline void CglGMI::flip(double& rowElem, int index) {
317 if ((index < ncol && cstat[index] == 2) ||
318 (index >= ncol && rstat[index-ncol] == 2)) {
319 rowElem = -rowElem;
320 }
321 } /* flip */
322
323 /************************************************************************/
unflipOrig(double & rowElem,int index,double & rowRhs)324 inline void CglGMI::unflipOrig(double& rowElem, int index, double& rowRhs) {
325 if (cstat[index] == 2) {
326 // structural variable at upper bound
327 rowElem = -rowElem;
328 rowRhs += rowElem*colUpper[index];
329 }
330 else if (cstat[index] == 3) {
331 // structural variable at lower bound
332 rowRhs += rowElem*colLower[index];
333 }
334 } /* unflipOrig */
335
336 /************************************************************************/
unflipSlack(double & rowElem,int index,double & rowRhs,const double * slackVal)337 inline void CglGMI::unflipSlack(double& rowElem, int index, double& rowRhs,
338 const double* slackVal) {
339 if (rstat[index-ncol] == 2) {
340 // artificial variable at upper bound
341 rowElem = -rowElem;
342 rowRhs += rowElem*slackVal[index-ncol];
343 }
344 else if (rstat[index-ncol] == 3) {
345 // artificial variable at lower bound
346 rowRhs += rowElem*slackVal[index-ncol];
347 }
348
349 } /* unflipSlack */
350
351 /************************************************************************/
packRow(double * row,double * rowElem,int * rowIndex,int & rowNz)352 inline void CglGMI::packRow(double* row, double* rowElem, int* rowIndex,
353 int& rowNz) {
354 rowNz = 0;
355 for (int i = 0; i < ncol; ++i) {
356 if (!isZero(fabs(row[i]))) {
357 rowElem[rowNz] = row[i];
358 rowIndex[rowNz] = i;
359 rowNz++;
360 }
361 }
362 }
363
364 /************************************************************************/
cleanCut(double * cutElem,int * cutIndex,int & cutNz,double & cutRhs,const double * xbar)365 bool CglGMI::cleanCut(double* cutElem, int* cutIndex, int& cutNz,
366 double& cutRhs, const double* xbar) {
367 CglGMIParam::CleaningProcedure cleanProc = param.getCLEAN_PROC();
368 if (cleanProc == CglGMIParam::CP_CGLLANDP1) {
369 if (!checkViolation(cutElem, cutIndex, cutNz, cutRhs, xbar)) {
370 #if defined GMI_TRACE_CLEAN
371 printf("CglGMI::cleanCut(): cut discarded: bad violation\n");
372 #endif
373 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
374 if (trackRejection) {
375 violFail++;
376 }
377 #endif
378 return false;
379 }
380 relaxRhs(cutRhs);
381 removeSmallCoefficients(cutElem, cutIndex, cutNz, cutRhs);
382 if (!checkSupport(cutNz)) {
383 #if defined GMI_TRACE_CLEAN
384 printf("CglGMI::cleanCut(): cut discarded: too large support\n");
385 #endif
386 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
387 if (trackRejection) {
388 suppFail++;
389 }
390 #endif
391 return false;
392 }
393 if (!checkDynamism(cutElem, cutIndex, cutNz)) {
394 #if defined GMI_TRACE_CLEAN
395 printf("CglGMI::cleanCut(): cut discarded: bad dynamism\n");
396 #endif
397 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
398 if (trackRejection) {
399 dynFail++;
400 }
401 #endif
402 return false;
403 }
404 if (!checkViolation(cutElem, cutIndex, cutNz, cutRhs, xbar)) {
405 #if defined GMI_TRACE_CLEAN
406 printf("CglGMI::cleanCut(): cut discarded: bad violation (final check)\n");
407 #endif
408 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
409 if (trackRejection) {
410 violFail++;
411 }
412 #endif
413 return false;
414 }
415 } /* end of cleaning procedure CP_CGLLANDP1 */
416 else if (cleanProc == CglGMIParam::CP_CGLLANDP2) {
417 if (!checkViolation(cutElem, cutIndex, cutNz, cutRhs, xbar)) {
418 #if defined GMI_TRACE_CLEAN
419 printf("CglGMI::cleanCut(): cut discarded: bad violation\n");
420 #endif
421 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
422 if (trackRejection) {
423 violFail++;
424 }
425 #endif
426 return false;
427 }
428 relaxRhs(cutRhs);
429 if (!checkDynamism(cutElem, cutIndex, cutNz)) {
430 #if defined GMI_TRACE_CLEAN
431 printf("CglGMI::cleanCut(): cut discarded: bad dynamism\n");
432 #endif
433 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
434 if (trackRejection) {
435 dynFail++;
436 }
437 #endif
438 return false;
439 }
440 if (!scaleCut(cutElem, cutIndex, cutNz, cutRhs, 1) &&
441 param.getENFORCE_SCALING()) {
442 #if defined GMI_TRACE_CLEAN
443 printf("CglGMI::cleanCut(): cut discarded: bad scaling\n");
444 #endif
445 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
446 if (trackRejection) {
447 scaleFail++;
448 }
449 #endif
450 return false;
451 }
452 removeSmallCoefficients(cutElem, cutIndex, cutNz, cutRhs);
453 if (!checkSupport(cutNz)) {
454 #if defined GMI_TRACE_CLEAN
455 printf("CglGMI::cleanCut(): cut discarded: too large support\n");
456 #endif
457 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
458 if (trackRejection) {
459 suppFail++;
460 }
461 #endif
462 return false;
463 }
464 if (!checkViolation(cutElem, cutIndex, cutNz, cutRhs, xbar)) {
465 #if defined GMI_TRACE_CLEAN
466 printf("CglGMI::cleanCut(): cut discarded: bad violation (final check)\n");
467 #endif
468 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
469 if (trackRejection) {
470 violFail++;
471 }
472 #endif
473 return false;
474 }
475 } /* end of cleaning procedure CP_CGLLANDP2 */
476 else if (cleanProc == CglGMIParam::CP_CGLREDSPLIT) {
477 if (!scaleCut(cutElem, cutIndex, cutNz, cutRhs, 3) &&
478 param.getENFORCE_SCALING()) {
479 #if defined GMI_TRACE_CLEAN
480 printf("CglGMI::cleanCut(): cut discarded: bad scaling\n");
481 #endif
482 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
483 if (trackRejection) {
484 scaleFail++;
485 }
486 #endif
487 return false;
488 }
489 removeSmallCoefficients(cutElem, cutIndex, cutNz, cutRhs);
490 if (!checkDynamism(cutElem, cutIndex, cutNz)) {
491 #if defined GMI_TRACE_CLEAN
492 printf("CglGMI::cleanCut(): cut discarded: bad dynamism\n");
493 #endif
494 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
495 if (trackRejection) {
496 dynFail++;
497 }
498 #endif
499 return false;
500 }
501 if (!checkSupport(cutNz)) {
502 #if defined GMI_TRACE_CLEAN
503 printf("CglGMI::cleanCut(): cut discarded: too large support\n");
504 #endif
505 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
506 if (trackRejection) {
507 suppFail++;
508 }
509 #endif
510 return false;
511 }
512 if (!checkViolation(cutElem, cutIndex, cutNz, cutRhs, xbar)) {
513 #if defined GMI_TRACE_CLEAN
514 printf("CglGMI::cleanCut(): cut discarded: bad violation (final check)\n");
515 #endif
516 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
517 if (trackRejection) {
518 violFail++;
519 }
520 #endif
521 return false;
522 }
523 relaxRhs(cutRhs);
524 } /* end of cleaning procedure CP_CGLREDSPLIT */
525 else if (cleanProc == CglGMIParam::CP_INTEGRAL_CUTS) {
526 removeSmallCoefficients(cutElem, cutIndex, cutNz, cutRhs);
527 if (!checkSupport(cutNz)) {
528 #if defined GMI_TRACE_CLEAN
529 printf("CglGMI::cleanCut(): cut discarded: too large support\n");
530 #endif
531 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
532 if (trackRejection) {
533 suppFail++;
534 }
535 #endif
536 return false;
537 }
538 if (!checkDynamism(cutElem, cutIndex, cutNz)) {
539 #if defined GMI_TRACE_CLEAN
540 printf("CglGMI::cleanCut(): cut discarded: bad dynamism\n");
541 #endif
542 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
543 if (trackRejection) {
544 dynFail++;
545 }
546 #endif
547 return false;
548 }
549 if (!scaleCut(cutElem, cutIndex, cutNz, cutRhs, 0) &&
550 param.getENFORCE_SCALING()) {
551 #if defined GMI_TRACE_CLEAN
552 printf("CglGMI::cleanCut(): cut discarded: bad scaling\n");
553 #endif
554 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
555 if (trackRejection) {
556 scaleFail++;
557 }
558 #endif
559 return false;
560 }
561 if (!checkViolation(cutElem, cutIndex, cutNz, cutRhs, xbar)) {
562 #if defined GMI_TRACE_CLEAN
563 printf("CglGMI::cleanCut(): cut discarded: bad violation (final check)\n");
564 #endif
565 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
566 if (trackRejection) {
567 violFail++;
568 }
569 #endif
570 return false;
571 }
572 } /* end of cleaning procedure CP_INTEGRAL_CUTS */
573 else if (cleanProc == CglGMIParam::CP_CGLLANDP1_INT) {
574 if (!checkViolation(cutElem, cutIndex, cutNz, cutRhs, xbar)) {
575 #if defined GMI_TRACE_CLEAN
576 printf("CglGMI::cleanCut(): cut discarded: bad violation\n");
577 #endif
578 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
579 if (trackRejection) {
580 violFail++;
581 }
582 #endif
583 return false;
584 }
585 removeSmallCoefficients(cutElem, cutIndex, cutNz, cutRhs);
586 if (!checkSupport(cutNz)) {
587 #if defined GMI_TRACE_CLEAN
588 printf("CglGMI::cleanCut(): cut discarded: too large support\n");
589 #endif
590 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
591 if (trackRejection) {
592 suppFail++;
593 }
594 #endif
595 return false;
596 }
597 if (!checkDynamism(cutElem, cutIndex, cutNz)) {
598 #if defined GMI_TRACE_CLEAN
599 printf("CglGMI::cleanCut(): cut discarded: bad dynamism\n");
600 #endif
601 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
602 if (trackRejection) {
603 dynFail++;
604 }
605 #endif
606 return false;
607 }
608 // scale cut so that it becomes integral, if possible
609 if (!scaleCut(cutElem, cutIndex, cutNz, cutRhs, 0)) {
610 if (param.getENFORCE_SCALING()){
611 #if defined GMI_TRACE_CLEAN
612 printf("CglGMI::cleanCut(): cut discarded: bad scaling\n");
613 #endif
614 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
615 if (trackRejection) {
616 scaleFail++;
617 }
618 #endif
619 return false;
620 }
621 else {
622 // If cannot scale to integral and not enforcing, relax rhs
623 // (as per CglLandP cleaning procedure)
624 relaxRhs(cutRhs);
625 }
626 }
627 if (!checkViolation(cutElem, cutIndex, cutNz, cutRhs, xbar)) {
628 #if defined GMI_TRACE_CLEAN
629 printf("CglGMI::cleanCut(): cut discarded: bad violation (final check)\n");
630 #endif
631 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
632 if (trackRejection) {
633 violFail++;
634 }
635 #endif
636 return false;
637 }
638 } /* end of cleaning procedure CP_CGLLANDP1_INT */
639 else if (cleanProc == CglGMIParam::CP_CGLLANDP1_SCALEMAX ||
640 cleanProc == CglGMIParam::CP_CGLLANDP1_SCALERHS) {
641 if (!checkViolation(cutElem, cutIndex, cutNz, cutRhs, xbar)) {
642 #if defined GMI_TRACE_CLEAN
643 printf("CglGMI::cleanCut(): cut discarded: bad violation\n");
644 #endif
645 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
646 if (trackRejection) {
647 violFail++;
648 }
649 #endif
650 return false;
651 }
652 if (// Try to scale cut, but do not discard if cannot scale
653 ((cleanProc == CglGMIParam::CP_CGLLANDP1_SCALEMAX &&
654 !scaleCut(cutElem, cutIndex, cutNz, cutRhs, 1)) ||
655 (cleanProc == CglGMIParam::CP_CGLLANDP1_SCALERHS &&
656 !scaleCut(cutElem, cutIndex, cutNz, cutRhs, 2))) &&
657 param.getENFORCE_SCALING()) {
658 #if defined GMI_TRACE_CLEAN
659 printf("CglGMI::cleanCut(): cut discarded: bad scaling\n");
660 #endif
661 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
662 if (trackRejection) {
663 scaleFail++;
664 }
665 #endif
666 return false;
667 }
668 relaxRhs(cutRhs);
669 removeSmallCoefficients(cutElem, cutIndex, cutNz, cutRhs);
670 if (!checkSupport(cutNz)) {
671 #if defined GMI_TRACE_CLEAN
672 printf("CglGMI::cleanCut(): cut discarded: too large support\n");
673 #endif
674 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
675 if (trackRejection) {
676 suppFail++;
677 }
678 #endif
679 return false;
680 }
681 if (!checkDynamism(cutElem, cutIndex, cutNz)) {
682 #if defined GMI_TRACE_CLEAN
683 printf("CglGMI::cleanCut(): cut discarded: bad dynamism\n");
684 #endif
685 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
686 if (trackRejection) {
687 dynFail++;
688 }
689 #endif
690 return false;
691 }
692 if (!checkViolation(cutElem, cutIndex, cutNz, cutRhs, xbar)) {
693 #if defined GMI_TRACE_CLEAN
694 printf("CglGMI::cleanCut(): cut discarded: bad violation (final check)\n");
695 #endif
696 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
697 if (trackRejection) {
698 violFail++;
699 }
700 #endif
701 return false;
702 }
703 } /* end of cleaning procedures CP_CGLLANDP1_SCALEMAX and CG_CGLLANDP1_SCALERHS */
704 return true;
705 }
706
707 /************************************************************************/
checkViolation(const double * cutElem,const int * cutIndex,int cutNz,double cutrhs,const double * xbar)708 bool CglGMI::checkViolation(const double* cutElem, const int* cutIndex,
709 int cutNz, double cutrhs, const double* xbar) {
710 double lhs = 0.0;
711 for (int i = 0; i < cutNz; ++i) {
712 lhs += cutElem[i]*xbar[cutIndex[i]];
713 }
714 double violation = lhs - cutrhs;
715 if (fabs(cutrhs) > 1) {
716 violation /= fabs(cutrhs);
717 }
718 if (violation >= param.getMINVIOL()) {
719 return true;
720 }
721 else{
722 #if defined GMI_TRACE_CLEAN
723 printf("Cut lhs %g, rhs %g, violation %g; cut discarded\n", lhs, cutrhs, violation);
724 #endif
725 return false;
726 }
727 } /* checkViolation */
728
729 /************************************************************************/
checkDynamism(const double * cutElem,const int * cutIndex,int cutNz)730 bool CglGMI::checkDynamism(const double* cutElem, const int* cutIndex,
731 int cutNz) {
732 double min = param.getINFINIT();
733 double max = 0.0;
734 double val = 0.0;
735 for (int i = 0; i < cutNz; ++i) {
736 if (!isZero(cutElem[i])) {
737 val = fabs(cutElem[i]);
738 min = CoinMin(min, val);
739 max = CoinMax(max, val);
740 }
741 }
742 if (max > min*param.getMAXDYN()) {
743 #if defined GMI_TRACE_CLEAN
744 printf("Max elem %g, min elem %g, dyn %g; cut discarded\n", max, min, max/min);
745 #endif
746 return false;
747 }
748 else{
749 return true;
750 }
751
752 } /* checkDynamism */
753
754 /************************************************************************/
checkSupport(int cutNz)755 bool CglGMI::checkSupport(int cutNz) {
756 if (cutNz > param.getMAX_SUPPORT_ABS() + param.getMAX_SUPPORT_REL()*ncol) {
757 #if defined GMI_TRACE_CLEAN
758 printf("Support %d; cut discarded\n", cutNz);
759 #endif
760 return false;
761 }
762 else{
763 return true;
764 }
765 }
766
767 /************************************************************************/
removeSmallCoefficients(double * cutElem,int * cutIndex,int & cutNz,double & cutRhs)768 bool CglGMI::removeSmallCoefficients(double* cutElem, int* cutIndex,
769 int& cutNz, double& cutRhs) {
770 double value, absval;
771 int currPos = 0;
772 int col;
773 for (int i = 0; i < cutNz; ++i) {
774 col = cutIndex[i];
775 value = cutElem[i];
776 absval = fabs(value);
777 if (!isZero(absval) && absval <= param.getEPS_COEFF()) {
778 // small coefficient: remove and adjust rhs if possible
779 if ((value > 0.0) && (colLower[col] > -param.getINFINIT())) {
780 cutRhs -= value * colLower[col];
781 }
782 else if ((value < 0.0) && (colUpper[col] < param.getINFINIT())) {
783 cutRhs -= value * colUpper[col];
784 }
785 }
786 else if (absval > param.getEPS_COEFF()) {
787 if (currPos < i) {
788 cutElem[currPos] = cutElem[i];
789 cutIndex[currPos] = cutIndex[i];
790 }
791 currPos++;
792 }
793 }
794 cutNz = currPos;
795 return true;
796 }
797
798 /************************************************************************/
relaxRhs(double & rhs)799 void CglGMI::relaxRhs(double& rhs) {
800 if(param.getEPS_RELAX_REL() > 0.0) {
801 rhs += fabs(rhs) * param.getEPS_RELAX_REL() + param.getEPS_RELAX_ABS();
802 }
803 else{
804 rhs += param.getEPS_RELAX_ABS();
805 }
806 }
807
808 /************************************************************************/
scaleCut(double * cutElem,int * cutIndex,int cutNz,double & cutRhs,int scalingType)809 bool CglGMI::scaleCut(double* cutElem, int* cutIndex, int cutNz,
810 double& cutRhs, int scalingType) {
811 /// scalingType possible values:
812 /// 0 : scale to obtain integral cut
813 /// 1 : scale to obtain largest coefficient equal to 1
814 /// 2 : scale to obtain rhs equal to 1
815 /// 3 : scale based on norm, to obtain cut norm equal to ncol
816 /// Returns true if scaling is successful.
817 if (scalingType == 0) {
818 return scaleCutIntegral(cutElem, cutIndex, cutNz, cutRhs);
819 }
820 else if (scalingType == 1) {
821 double max = fabs(cutRhs);
822 for (int i = 0; i < cutNz; ++i) {
823 if (!isZero(cutElem[i])) {
824 max = CoinMax(max, fabs(cutElem[i]));
825 }
826 }
827 if (max < param.getEPS() || max > param.getMAXDYN()) {
828 #if defined GMI_TRACE_CLEAN
829 printf("Scale %g; %g %g cut discarded\n", max, param.getEPS(), 1/param.getMAXDYN());
830 #endif
831 return false;
832 }
833 else{
834 for (int i = 0; i < cutNz; ++i) {
835 cutElem[i] /= max;
836 }
837 cutRhs /= max;
838 return true;
839 }
840 }
841 else if (scalingType == 2) {
842 double max = fabs(cutRhs);
843 if (max < param.getEPS() || max > param.getMAXDYN()) {
844 #if defined GMI_TRACE_CLEAN
845 printf("Scale %g; %g %g cut discarded\n", max, param.getEPS(), 1/param.getMAXDYN());
846 #endif
847 return false;
848 }
849 else{
850 for (int i = 0; i < cutNz; ++i) {
851 cutElem[i] /= max;
852 }
853 cutRhs /= max;
854 return true;
855 }
856 }
857 else if (scalingType == 3) {
858 int support = 0;
859 double norm = 0.0;
860 for (int i = 0; i < cutNz; ++i) {
861 if (!isZero(fabs(cutElem[i]))) {
862 support++;
863 norm += cutElem[i]*cutElem[i];
864 }
865 }
866 double scale = sqrt(norm / support);
867 if ((scale < 0.02) || (scale > 100)) {
868 #if defined GMI_TRACE_CLEAN
869 printf("Scale %g; cut discarded\n", scale);
870 #endif
871 return false;
872 }
873 else{
874 for (int i = 0; i < cutNz; ++i) {
875 cutElem[i] /= scale;
876 }
877 cutRhs /= scale;
878 return true;
879 }
880 }
881 return false;
882 } /* scaleCut */
883
884 /************************************************************************/
scaleCutIntegral(double * cutElem,int * cutIndex,int cutNz,double & cutRhs)885 bool CglGMI::scaleCutIntegral(double* cutElem, int* cutIndex, int cutNz,
886 double& cutRhs) {
887 long gcd, lcm;
888 double maxdelta = param.getEPS();
889 double maxscale = 1000;
890 long maxdnom = 1000;
891 long numerator = 0, denominator = 0;
892 // Initialize gcd and lcm
893 CoinRational r = CoinRational(cutRhs, maxdelta, maxdnom);
894 if (r.getNumerator() != 0){
895 gcd = labs(r.getNumerator());
896 lcm = r.getDenominator();
897 }
898 else{
899 #if defined GMI_TRACE_CLEAN
900 printf("Cannot compute rational number, scaling procedure aborted\n");
901 #endif
902 return false;
903 }
904 for (int i = 0; i < cutNz; ++i) {
905 if (solver->isContinuous(cutIndex[i]) && !param.getINTEGRAL_SCALE_CONT()) {
906 continue;
907 }
908 CoinRational r = CoinRational(cutElem[i], maxdelta, maxdnom);
909 if (r.getNumerator() != 0){
910 gcd = computeGcd(gcd, r.getNumerator());
911 lcm *= r.getDenominator()/(computeGcd(lcm,r.getDenominator()));
912 }
913 else{
914 #if defined GMI_TRACE_CLEAN
915 printf("Cannot compute rational number, scaling procedure aborted\n");
916 #endif
917 return false;
918 }
919 }
920 double scale = ((double)lcm)/((double)gcd);
921 if (fabs(scale) > maxscale) {
922 #if defined GMI_TRACE_CLEAN
923 printf("Scaling factor too large, scaling procedure aborted\n");
924 #endif
925 return false;
926 }
927 // Looks like we have a good scaling factor; scale and return;
928 for (int i = 0; i < cutNz; ++i) {
929 cutElem[i] *= scale;
930 }
931 cutRhs *= scale;
932 return true;
933 } /* scaleCutIntegral */
934
935 /************************************************************************/
computeGcd(long a,long b)936 long CglGMI::computeGcd(long a, long b) {
937 // This is the standard Euclidean algorithm for gcd
938 long remainder = 1;
939 // Make sure a<=b (will always remain so)
940 if (a > b) {
941 // Swap a and b
942 long temp = a;
943 a = b;
944 b = temp;
945 }
946 // If zero then gcd is nonzero
947 if (!a) {
948 if (b) {
949 return b;
950 }
951 else {
952 printf("### WARNING: CglGMI::computeGcd() given two zeroes!\n");
953 exit(1);
954 }
955 }
956 while (remainder) {
957 remainder = b % a;
958 b = a;
959 a = remainder;
960 }
961 return b;
962 } /* computeGcd */
963
964
965 /************************************************************************/
generateCuts(const OsiSolverInterface & si,OsiCuts & cs,const CglTreeInfo)966 void CglGMI::generateCuts(const OsiSolverInterface &si, OsiCuts & cs,
967 const CglTreeInfo )
968 {
969 solver = const_cast<OsiSolverInterface *>(&si);
970 if (solver == NULL) {
971 printf("### WARNING: CglGMI::generateCuts(): no solver available.\n");
972 return;
973 }
974
975 if (!solver->optimalBasisIsAvailable()) {
976 printf("### WARNING: CglGMI::generateCuts(): no optimal basis available.\n");
977 return;
978 }
979
980 #if defined OSI_TABLEAU
981 if (!solver->canDoSimplexInterface()) {
982 printf("### WARNING: CglGMI::generateCuts(): solver does not provide simplex tableau.\n");
983 printf("### WARNING: CglGMI::generateCuts(): recompile without OSI_TABLEAU.\n");
984 return;
985 }
986 #endif
987
988
989 // Get basic problem information from solver
990 ncol = solver->getNumCols();
991 nrow = solver->getNumRows();
992 colLower = solver->getColLower();
993 colUpper = solver->getColUpper();
994 rowLower = solver->getRowLower();
995 rowUpper = solver->getRowUpper();
996 rowRhs = solver->getRightHandSide();
997
998 xlp = solver->getColSolution();
999 rowActivity = solver->getRowActivity();
1000 byRow = solver->getMatrixByRow();
1001 byCol = solver->getMatrixByCol();
1002
1003 generateCuts(cs);
1004
1005 } /* generateCuts */
1006
1007 /************************************************************************/
generateCuts(OsiCuts & cs)1008 void CglGMI::generateCuts(OsiCuts &cs)
1009 {
1010 isInteger = new bool[ncol];
1011
1012 computeIsInteger();
1013
1014 cstat = new int[ncol];
1015 rstat = new int[nrow];
1016
1017
1018 solver->getBasisStatus(cstat, rstat); // 0: free 1: basic
1019 // 2: upper 3: lower
1020
1021
1022 #if defined GMI_TRACETAB
1023 printvecINT("cstat", cstat, ncol);
1024 printvecINT("rstat", rstat, nrow);
1025 #endif
1026
1027 // list of basic integer fractional variables
1028 int *listFracBasic = new int[nrow];
1029 int numFracBasic = 0;
1030 for (int i = 0; i < ncol; ++i) {
1031 // j is the variable which is basic in row i
1032 if ((cstat[i] == 1) && (isInteger[i])) {
1033 if (CoinMin(aboveInteger(xlp[i]),
1034 1-aboveInteger(xlp[i])) > param.getAway()) {
1035 listFracBasic[numFracBasic] = i;
1036 numFracBasic++;
1037 }
1038 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
1039 else if (trackRejection) {
1040 // Say that we tried to generate a cut, but it was discarded
1041 // because of small fractionality
1042 if (!isIntegerValue(xlp[i])) {
1043 fracFail++;
1044 numGeneratedCuts++;
1045 }
1046 }
1047 #endif
1048 }
1049 }
1050
1051 #if defined GMI_TRACE
1052 printf("CglGMI::generateCuts() : %d fractional rows\n", numFracBasic);
1053 #endif
1054
1055 if (numFracBasic == 0) {
1056 delete[] listFracBasic;
1057 delete[] cstat;
1058 delete[] rstat;
1059 delete[] isInteger;
1060 return;
1061 }
1062
1063 // there are rows with basic integer fractional variables, so we can
1064 // generate cuts
1065
1066 // Basis index for columns and rows; each element is -1 if corresponding
1067 // variable is nonbasic, and contains the basis index if basic.
1068 // The basis index is the row in which the variable is basic.
1069 int* colBasisIndex = new int[ncol];
1070 int* rowBasisIndex = new int[nrow];
1071
1072 #if defined OSI_TABLEAU
1073 memset(colBasisIndex, -1, ncol*sizeof(int));
1074 memset(rowBasisIndex, -1, nrow*sizeof(int));
1075 solver->enableFactorization();
1076 int* basicVars = new int[nrow];
1077 solver->getBasics(basicVars);
1078 for (int i = 0; i < nrow; ++i) {
1079 if (basicVars[i] < ncol) {
1080 colBasisIndex[basicVars[i]] = i;
1081 }
1082 else {
1083 rowBasisIndex[basicVars[i] - ncol] = i;
1084 }
1085 }
1086 #else
1087 CoinFactorization factorization;
1088 if (factorize(factorization, colBasisIndex, rowBasisIndex)) {
1089 printf("### WARNING: CglGMI::generateCuts(): error during factorization!\n");
1090 return;
1091 }
1092 #endif
1093
1094
1095 // cut in sparse form
1096 double* cutElem = new double[ncol];
1097 int* cutIndex = new int[ncol];
1098 int cutNz = 0;
1099 double cutRhs;
1100
1101 // cut in dense form
1102 double* cut = new double[ncol];
1103
1104 double *slackVal = new double[nrow];
1105
1106 for (int i = 0; i < nrow; ++i) {
1107 slackVal[i] = rowRhs[i] - rowActivity[i];
1108 }
1109
1110 #if defined OSI_TABLEAU
1111 // Column part and row part of a row of the simplex tableau
1112 double* tableauColPart = new double[ncol];
1113 double* tableauRowPart = new double[nrow];
1114 #else
1115 // Need some more data for simplex tableau computation
1116 const int * row = byCol->getIndices();
1117 const CoinBigIndex * columnStart = byCol->getVectorStarts();
1118 const int * columnLength = byCol->getVectorLengths();
1119 const double * columnElements = byCol->getElements();
1120
1121 // Create work arrays for factorization
1122 // two vectors for updating: the first one is needed to do the computations
1123 // but we do not use it, the second one contains a row of the basis inverse
1124 CoinIndexedVector work;
1125 CoinIndexedVector array;
1126 // Make sure they large enough
1127 work.reserve(nrow);
1128 array.reserve(nrow);
1129 int * arrayRows = array.getIndices();
1130 double * arrayElements = array.denseVector();
1131 // End of code to create work arrays
1132 double one = 1.0;
1133 #endif
1134
1135 // Matrix elements by row for slack substitution
1136 const double *elements = byRow->getElements();
1137 const CoinBigIndex *rowStart = byRow->getVectorStarts();
1138 const int *indices = byRow->getIndices();
1139 const int *rowLength = byRow->getVectorLengths();
1140
1141 // Indices of basic and slack variables, and cut elements
1142 int iBasic, slackIndex;
1143 double cutCoeff;
1144 double rowElem;
1145 // Now generate the cuts: obtain a row of the simplex tableau
1146 // where an integer variable is basic and fractional, and compute the cut
1147 for (int i = 0; i < numFracBasic; ++i) {
1148 if (!computeCutFractionality(xlp[listFracBasic[i]], cutRhs)) {
1149 // cut is discarded because of the small fractionalities involved
1150 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
1151 if (trackRejection) {
1152 // Say that we tried to generate a cut, but it was discarded
1153 // because of small fractionality
1154 fracFail++;
1155 numGeneratedCuts++;
1156 }
1157 #endif
1158 continue;
1159 }
1160
1161 // the variable listFracBasic[i] is basic in row iBasic
1162 iBasic = colBasisIndex[listFracBasic[i]];
1163
1164 #if defined GMI_TRACE
1165 printf("Row %d with var %d basic, f0 = %f\n", i, listFracBasic[i], f0);
1166 #endif
1167
1168 #if defined OSI_TABLEAU
1169 solver->getBInvARow(iBasic, tableauColPart, tableauRowPart);
1170 #else
1171 array.clear();
1172 array.setVector(1, &iBasic, &one);
1173
1174 factorization.updateColumnTranspose (&work, &array);
1175
1176 int numberInArray=array.getNumElements();
1177 #endif
1178
1179 // reset the cut
1180 memset(cut, 0, ncol*sizeof(double));
1181
1182 // columns
1183 for (int j = 0; j < ncol; ++j) {
1184 if ((colBasisIndex[j] >= 0) ||
1185 (areEqual(colLower[j], colUpper[j],
1186 param.getEPS(), param.getEPS()))) {
1187 // Basic or fixed variable -- skip
1188 continue;
1189 }
1190 #ifdef OSI_TABLEAU
1191 rowElem = tableauColPart[j];
1192 #else
1193 rowElem = 0.0;
1194 // add in row of tableau
1195 for (CoinBigIndex h = columnStart[j]; h < columnStart[j]+columnLength[j]; ++h) {
1196 rowElem += columnElements[h]*arrayElements[row[h]];
1197 }
1198 #endif
1199 if (!isZero(fabs(rowElem))) {
1200 // compute cut coefficient
1201 flip(rowElem, j);
1202 cutCoeff = computeCutCoefficient(rowElem, j);
1203 if (isZero(cutCoeff)) {
1204 continue;
1205 }
1206 unflipOrig(cutCoeff, j, cutRhs);
1207 cut[j] = cutCoeff;
1208 #if defined GMI_TRACE
1209 printf("var %d, row %f, cut %f\n", j, rowElem, cutCoeff);
1210 #endif
1211 }
1212 }
1213
1214 // now do slacks part
1215 #if defined OSI_TABLEAU
1216 for (int j = 0 ; j < nrow; ++j) {
1217 // index of the row corresponding to the slack variable
1218 slackIndex = j;
1219 if (rowBasisIndex[j] >= 0) {
1220 // Basic variable -- skip it
1221 continue;
1222 }
1223 rowElem = tableauRowPart[j];
1224 #else
1225 for (int j = 0 ; j < numberInArray ; ++j) {
1226 // index of the row corresponding to the slack variable
1227 slackIndex = arrayRows[j];
1228 rowElem = arrayElements[slackIndex];
1229 #endif
1230 if (!isZero(fabs(rowElem))) {
1231 slackIndex += ncol;
1232 // compute cut coefficient
1233 flip(rowElem, slackIndex);
1234 cutCoeff = computeCutCoefficient(rowElem, slackIndex);
1235 if (isZero(fabs(cutCoeff))) {
1236 continue;
1237 }
1238 unflipSlack(cutCoeff, slackIndex, cutRhs, slackVal);
1239 eliminateSlack(cutCoeff, slackIndex, cut, cutRhs,
1240 elements, rowStart, indices, rowLength, rowRhs);
1241 #if defined GMI_TRACE
1242 printf("var %d, row %f, cut %f\n", slackIndex, rowElem, cutCoeff);
1243 #endif
1244 }
1245 }
1246
1247 packRow(cut, cutElem, cutIndex, cutNz);
1248 if (cutNz == 0)
1249 continue;
1250
1251 #if defined GMI_TRACE
1252 printvecDBL("final cut:", cutElem, cutIndex, cutNz);
1253 printf("cutRhs: %f\n", cutRhs);
1254 #endif
1255
1256 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
1257 if (trackRejection) {
1258 numGeneratedCuts++;
1259 }
1260 #endif
1261 if (cleanCut(cutElem, cutIndex, cutNz, cutRhs, xlp) && cutNz > 0) {
1262 OsiRowCut rc;
1263 rc.setRow(cutNz, cutIndex, cutElem);
1264 rc.setLb(-param.getINFINIT());
1265 rc.setUb(cutRhs);
1266 if (!param.getCHECK_DUPLICATES()) {
1267 cs.insertIfNotDuplicate(rc);
1268 }
1269 else{
1270 cs.insertIfNotDuplicate(rc, CoinAbsFltEq(param.getEPS_COEFF()));
1271 }
1272 }
1273
1274 }
1275
1276 #if defined GMI_TRACE
1277 printf("CglGMI::generateCuts() : number of cuts : %d\n", cs.sizeRowCuts());
1278 #endif
1279
1280 #if defined OSI_TABLEAU
1281 solver->disableFactorization();
1282 delete[] basicVars;
1283 delete[] tableauColPart;
1284 delete[] tableauRowPart;
1285 #endif
1286
1287 delete[] colBasisIndex;
1288 delete[] rowBasisIndex;
1289 delete[] cut;
1290 delete[] slackVal;
1291 delete[] cutElem;
1292 delete[] cutIndex;
1293 delete[] listFracBasic;
1294 delete[] cstat;
1295 delete[] rstat;
1296 delete[] isInteger;
1297
1298 } /* generateCuts */
1299
1300 /***********************************************************************/
1301 void CglGMI::setParam(const CglGMIParam &source) {
1302 param = source;
1303 } /* setParam */
1304
1305 /***********************************************************************/
1306 void CglGMI::computeIsInteger() {
1307 for (int i = 0; i < ncol; ++i) {
1308 if(solver->isInteger(i)) {
1309 isInteger[i] = true;
1310 }
1311 else {
1312 if((areEqual(colLower[i], colUpper[i],
1313 param.getEPS(), param.getEPS()))
1314 && (isIntegerValue(colUpper[i]))) {
1315 // continuous variable fixed to an integer value
1316 isInteger[i] = true;
1317 }
1318 else {
1319 isInteger[i] = false;
1320 }
1321 }
1322 }
1323 } /* computeIsInteger */
1324
1325 /***********************************************************************/
1326 void CglGMI::printOptTab(OsiSolverInterface *lclSolver) const
1327 {
1328 int *cstat = new int[ncol];
1329 int *rstat = new int[nrow];
1330
1331 lclSolver->enableFactorization();
1332 lclSolver->getBasisStatus(cstat, rstat); // 0: free 1: basic
1333 // 2: upper 3: lower
1334
1335 int *basisIndex = new int[nrow]; // basisIndex[i] =
1336 // index of pivot var in row i
1337 // (slack if number >= ncol)
1338 lclSolver->getBasics(basisIndex);
1339
1340 double *z = new double[ncol]; // workspace to get row of the tableau
1341 double *slack = new double[nrow]; // workspace to get row of the tableau
1342 double *slackVal = new double[nrow];
1343
1344 for (int i = 0; i < nrow; i++) {
1345 slackVal[i] = rowRhs[i] - rowActivity[i];
1346 }
1347
1348 const double *rc = lclSolver->getReducedCost();
1349 const double *dual = lclSolver->getRowPrice();
1350 const double *solution = lclSolver->getColSolution();
1351
1352 printvecINT("cstat", cstat, ncol);
1353 printvecINT("rstat", rstat, nrow);
1354 printvecINT("basisIndex", basisIndex, nrow);
1355
1356 printvecDBL("solution", solution, ncol);
1357 printvecDBL("slackVal", slackVal, nrow);
1358 printvecDBL("reduced_costs", rc, ncol);
1359 printvecDBL("dual solution", dual, nrow);
1360
1361 printf("Optimal Tableau:\n");
1362
1363 for (int i = 0; i < nrow; i++) {
1364 lclSolver->getBInvARow(i, z, slack);
1365 for (int ii = 0; ii < ncol; ++ii) {
1366 printf("%5.2f ", z[ii]);
1367 }
1368 printf(" | ");
1369 for (int ii = 0; ii < nrow; ++ii) {
1370 printf("%5.2f ", slack[ii]);
1371 }
1372 printf(" | ");
1373 if(basisIndex[i] < ncol) {
1374 printf("%5.2f ", solution[basisIndex[i]]);
1375 }
1376 else {
1377 printf("%5.2f ", slackVal[basisIndex[i]-ncol]);
1378 }
1379 printf("\n");
1380 }
1381 for (int ii = 0; ii < 7*(ncol+nrow+1); ++ii) {
1382 printf("-");
1383 }
1384 printf("\n");
1385
1386 for (int ii = 0; ii < ncol; ++ii) {
1387 printf("%5.2f ", rc[ii]);
1388 }
1389 printf(" | ");
1390 for (int ii = 0; ii < nrow; ++ii) {
1391 printf("%5.2f ", -dual[ii]);
1392 }
1393 printf(" | ");
1394 printf("%5.2f\n", -lclSolver->getObjValue());
1395 lclSolver->disableFactorization();
1396
1397 delete[] cstat;
1398 delete[] rstat;
1399 delete[] basisIndex;
1400 delete[] slack;
1401 delete[] z;
1402 delete[] slackVal;
1403 } /* printOptTab */
1404
1405
1406 /*********************************************************************/
1407 // Create C++ lines to get to current state
1408 std::string
1409 CglGMI::generateCpp(FILE * fp)
1410 {
1411 CglGMI other;
1412 fprintf(fp,"0#include \"CglGMI.hpp\"\n");
1413 fprintf(fp,"3 CglGMI GMI;\n");
1414 if (param.getMAX_SUPPORT()!=other.param.getMAX_SUPPORT())
1415 fprintf(fp,"3 GMI.setLimit(%d);\n",param.getMAX_SUPPORT());
1416 else
1417 fprintf(fp,"4 GMI.setLimit(%d);\n",param.getMAX_SUPPORT());
1418 if (param.getAway()!=other.param.getAway())
1419 fprintf(fp,"3 GMI.setAway(%g);\n",param.getAway());
1420 else
1421 fprintf(fp,"4 GMI.setAway(%g);\n",param.getAway());
1422 if (param.getEPS()!=other.param.getEPS())
1423 fprintf(fp,"3 GMI.setEPS(%g);\n",param.getEPS());
1424 else
1425 fprintf(fp,"4 GMI.setEPS(%g);\n",param.getEPS());
1426 if (param.getEPS_COEFF()!=other.param.getEPS_COEFF())
1427 fprintf(fp,"3 GMI.setEPS_COEFF(%g);\n",param.getEPS_COEFF());
1428 else
1429 fprintf(fp,"4 GMI.set.EPS_COEFF(%g);\n",param.getEPS_COEFF());
1430 if (param.getEPS_RELAX_ABS()!=other.param.getEPS_RELAX_ABS())
1431 fprintf(fp,"3 GMI.set.EPS_RELAX(%g);\n",param.getEPS_RELAX_ABS());
1432 else
1433 fprintf(fp,"4 GMI.set.EPS_RELAX(%g);\n",param.getEPS_RELAX_ABS());
1434 if (getAggressiveness()!=other.getAggressiveness())
1435 fprintf(fp,"3 GMI.setAggressiveness(%d);\n",getAggressiveness());
1436 else
1437 fprintf(fp,"4 GMI.setAggressiveness(%d);\n",getAggressiveness());
1438 return "GMI";
1439 }
1440
1441 /*********************************************************************/
1442 int
1443 CglGMI::factorize(CoinFactorization & factorization,
1444 int* colBasisIndex, int* rowBasisIndex) {
1445 // Start of code to create a factorization from warm start ====
1446 // Taken (with small modifications) from CglGomory
1447 int status=-100;
1448 for (int i = 0; i < nrow; ++i) {
1449 if (rstat[i] == 1) {
1450 rowBasisIndex[i]=1;
1451 } else {
1452 rowBasisIndex[i]=-1;
1453 }
1454 }
1455 for (int i = 0; i < ncol; ++i) {
1456 if (cstat[i] == 1) {
1457 colBasisIndex[i]=1;
1458 } else {
1459 colBasisIndex[i]=-1;
1460 }
1461 }
1462 // returns 0 if okay, -1 singular, -2 too many in basis, -99 memory */
1463 while (status<-98) {
1464 status=factorization.factorize(*byCol, rowBasisIndex, colBasisIndex);
1465 if (status==-99) factorization.areaFactor(factorization.areaFactor()*2.0);
1466 }
1467 if (status) {
1468 return -1;
1469 }
1470 #if defined GMI_TRACE
1471 double condition = 0.0;
1472 const CoinFactorizationDouble * pivotRegion = factorization.pivotRegion();
1473 for (int i = 0; i < nrow; ++i) {
1474 condition += log(fabs(pivotRegion[i]));
1475 }
1476 printf("CglGMI::factorize(): condition number recomputed as sum of log: %g\n", (condition));
1477 #endif
1478 return 0;
1479 }
1480
1481 /*********************************************************************/
1482 void CglGMI::setTrackRejection(bool value) {
1483 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
1484 trackRejection = value;
1485 if (trackRejection) {
1486 // reset data members
1487 resetRejectionCounters();
1488 }
1489 #endif
1490 }
1491
1492 /*********************************************************************/
1493 bool CglGMI::getTrackRejection() {
1494 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
1495 return trackRejection;
1496 #else
1497 return false;
1498 #endif
1499 }
1500
1501 /*********************************************************************/
1502 void CglGMI::resetRejectionCounters() {
1503 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
1504 fracFail = 0;
1505 dynFail = 0;
1506 violFail = 0;
1507 suppFail = 0;
1508 scaleFail = 0;
1509 numGeneratedCuts = 0;
1510 #endif
1511 }
1512
1513 /*********************************************************************/
1514 int CglGMI::getNumberRejectedCuts(RejectionType reason) {
1515 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
1516 switch (reason) {
1517 case failureFractionality:
1518 return fracFail;
1519 case failureDynamism:
1520 return dynFail;
1521 case failureViolation:
1522 return violFail;
1523 case failureSupport:
1524 return suppFail;
1525 case failureScale:
1526 return scaleFail;
1527 }
1528 return 0;
1529 #else
1530 return 0;
1531 #endif
1532 }
1533
1534 /*********************************************************************/
1535 int CglGMI::getNumberGeneratedCuts() {
1536 #if defined TRACK_REJECT || defined TRACK_REJECT_SIMPLE
1537 return numGeneratedCuts;
1538 #else
1539 return 0;
1540 #endif
1541 }
1542