1 // $Id$
2 // Copyright (C) 2000, International Business Machines
3 // Corporation and others. All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5
6 #include <cstdlib>
7 #include <cstdio>
8 #include <cmath>
9 #include <cfloat>
10 #include <cassert>
11 #include <iostream>
12
13 #include "CoinPragma.hpp"
14 #include "CglTreeInfo.hpp"
15 #include "CoinHelperFunctions.hpp"
16 #include "CoinSort.hpp"
17 #include "CoinPackedMatrix.hpp"
18 #include "CglStored.hpp"
19 #include "OsiRowCut.hpp"
20
21 // Default constructor
CglTreeInfo()22 CglTreeInfo::CglTreeInfo()
23 : level(-1)
24 , pass(-1)
25 , formulation_rows(-1)
26 , options(0)
27 , inTree(false)
28 , hasParent(0)
29 , parentSolver(NULL)
30 , originalColumns(NULL)
31 , strengthenRow(NULL)
32 , randomNumberGenerator(NULL)
33 {
34 }
35
36 // Copy constructor
CglTreeInfo(const CglTreeInfo & rhs)37 CglTreeInfo::CglTreeInfo(const CglTreeInfo &rhs)
38 : level(rhs.level)
39 , pass(rhs.pass)
40 , formulation_rows(rhs.formulation_rows)
41 , options(rhs.options)
42 , inTree(rhs.inTree)
43 , hasParent(rhs.hasParent)
44 , parentSolver(rhs.parentSolver)
45 , originalColumns(rhs.originalColumns)
46 , strengthenRow(rhs.strengthenRow)
47 , randomNumberGenerator(rhs.randomNumberGenerator)
48 {
49 }
50 // Clone
51 CglTreeInfo *
clone() const52 CglTreeInfo::clone() const
53 {
54 return new CglTreeInfo(*this);
55 }
56
57 // Assignment operator
58 CglTreeInfo &
operator =(const CglTreeInfo & rhs)59 CglTreeInfo::operator=(const CglTreeInfo &rhs)
60 {
61 if (this != &rhs) {
62 //CglCutGenerator::operator=(rhs);
63 level = rhs.level;
64 pass = rhs.pass;
65 formulation_rows = rhs.formulation_rows;
66 options = rhs.options;
67 inTree = rhs.inTree;
68 hasParent = rhs.hasParent;
69 parentSolver = rhs.parentSolver;
70 originalColumns = rhs.originalColumns;
71 strengthenRow = rhs.strengthenRow;
72 randomNumberGenerator = rhs.randomNumberGenerator;
73 }
74 return *this;
75 }
76
77 // Destructor
~CglTreeInfo()78 CglTreeInfo::~CglTreeInfo()
79 {
80 }
81
82 // Default constructor
CglTreeProbingInfo()83 CglTreeProbingInfo::CglTreeProbingInfo()
84 : CglTreeInfo()
85 , fixEntry_(NULL)
86 , toZero_(NULL)
87 , toOne_(NULL)
88 , integerVariable_(NULL)
89 , backward_(NULL)
90 , fixingEntry_(NULL)
91 , numberVariables_(0)
92 , numberIntegers_(0)
93 , maximumEntries_(0)
94 , numberEntries_(-1)
95 {
96 }
97 // Constructor from model
CglTreeProbingInfo(const OsiSolverInterface * model)98 CglTreeProbingInfo::CglTreeProbingInfo(const OsiSolverInterface *model)
99 : CglTreeInfo()
100 , fixEntry_(NULL)
101 , toZero_(NULL)
102 , toOne_(NULL)
103 , integerVariable_(NULL)
104 , backward_(NULL)
105 , fixingEntry_(NULL)
106 , numberVariables_(0)
107 , numberIntegers_(0)
108 , maximumEntries_(0)
109 , numberEntries_(-1)
110 {
111 numberVariables_ = model->getNumCols();
112 // Too many ... but
113 integerVariable_ = new int[numberVariables_];
114 backward_ = new int[numberVariables_];
115 int i;
116 // Get integer types
117 const char *columnType = model->getColType(true);
118 for (i = 0; i < numberVariables_; i++) {
119 backward_[i] = -1;
120 if (columnType[i]) {
121 if (columnType[i] == 1) {
122 backward_[i] = numberIntegers_;
123 integerVariable_[numberIntegers_++] = i;
124 } else {
125 backward_[i] = -2;
126 }
127 }
128 }
129 // Set up to arrays
130 toOne_ = new int[numberIntegers_];
131 toZero_ = new int[numberIntegers_ + 1];
132 // zero out
133 CoinZeroN(toOne_, numberIntegers_);
134 CoinZeroN(toZero_, numberIntegers_ + 1);
135 }
136
137 // Copy constructor
CglTreeProbingInfo(const CglTreeProbingInfo & rhs)138 CglTreeProbingInfo::CglTreeProbingInfo(const CglTreeProbingInfo &rhs)
139 : CglTreeInfo(rhs)
140 , fixEntry_(NULL)
141 , toZero_(NULL)
142 , toOne_(NULL)
143 , integerVariable_(NULL)
144 , backward_(NULL)
145 , fixingEntry_(NULL)
146 , numberVariables_(rhs.numberVariables_)
147 , numberIntegers_(rhs.numberIntegers_)
148 , maximumEntries_(rhs.maximumEntries_)
149 , numberEntries_(rhs.numberEntries_)
150 {
151 if (numberVariables_) {
152 fixEntry_ = new CliqueEntry[maximumEntries_];
153 memcpy(fixEntry_, rhs.fixEntry_, maximumEntries_ * sizeof(CliqueEntry));
154 if (numberEntries_ < 0) {
155 // in order
156 toZero_ = CoinCopyOfArray(rhs.toZero_, numberIntegers_ + 1);
157 toOne_ = CoinCopyOfArray(rhs.toOne_, numberIntegers_);
158 } else {
159 // not in order
160 fixingEntry_ = CoinCopyOfArray(rhs.fixingEntry_, maximumEntries_);
161 }
162 integerVariable_ = CoinCopyOfArray(rhs.integerVariable_, numberIntegers_);
163 backward_ = CoinCopyOfArray(rhs.backward_, numberVariables_);
164 }
165 }
166 // Clone
167 CglTreeInfo *
clone() const168 CglTreeProbingInfo::clone() const
169 {
170 return new CglTreeProbingInfo(*this);
171 }
172
173 // Assignment operator
174 CglTreeProbingInfo &
operator =(const CglTreeProbingInfo & rhs)175 CglTreeProbingInfo::operator=(const CglTreeProbingInfo &rhs)
176 {
177 if (this != &rhs) {
178 CglTreeInfo::operator=(rhs);
179 delete[] fixEntry_;
180 delete[] toZero_;
181 delete[] toOne_;
182 delete[] integerVariable_;
183 delete[] backward_;
184 delete[] fixingEntry_;
185 numberVariables_ = rhs.numberVariables_;
186 numberIntegers_ = rhs.numberIntegers_;
187 maximumEntries_ = rhs.maximumEntries_;
188 numberEntries_ = rhs.numberEntries_;
189 if (numberVariables_) {
190 fixEntry_ = new CliqueEntry[maximumEntries_];
191 memcpy(fixEntry_, rhs.fixEntry_, maximumEntries_ * sizeof(CliqueEntry));
192 if (numberEntries_ < 0) {
193 // in order
194 toZero_ = CoinCopyOfArray(rhs.toZero_, numberIntegers_ + 1);
195 toOne_ = CoinCopyOfArray(rhs.toOne_, numberIntegers_);
196 fixingEntry_ = NULL;
197 } else {
198 // not in order
199 fixingEntry_ = CoinCopyOfArray(rhs.fixingEntry_, maximumEntries_);
200 toZero_ = NULL;
201 toOne_ = NULL;
202 }
203 toZero_ = CoinCopyOfArray(rhs.toZero_, numberIntegers_ + 1);
204 toOne_ = CoinCopyOfArray(rhs.toOne_, numberIntegers_);
205 integerVariable_ = CoinCopyOfArray(rhs.integerVariable_, numberIntegers_);
206 backward_ = CoinCopyOfArray(rhs.backward_, numberVariables_);
207 } else {
208 fixEntry_ = NULL;
209 toZero_ = NULL;
210 toOne_ = NULL;
211 integerVariable_ = NULL;
212 backward_ = NULL;
213 fixingEntry_ = NULL;
214 }
215 }
216 return *this;
217 }
218
219 // Destructor
~CglTreeProbingInfo()220 CglTreeProbingInfo::~CglTreeProbingInfo()
221 {
222 delete[] fixEntry_;
223 delete[] toZero_;
224 delete[] toOne_;
225 delete[] integerVariable_;
226 delete[] backward_;
227 delete[] fixingEntry_;
228 }
outDupsEtc(int numberIntegers,int & numberCliques,int & numberMatrixCliques,CoinBigIndex * & cliqueStart,char * & cliqueType,CliqueEntry * & entry,int numberLastTime,int printit)229 static int outDupsEtc(int numberIntegers, int &numberCliques, int &numberMatrixCliques,
230 CoinBigIndex *&cliqueStart, char *&cliqueType, CliqueEntry *&entry,
231 int numberLastTime, int printit)
232 {
233 bool allNew = false;
234 int *whichP = new int[numberIntegers];
235 int iClique;
236 assert(sizeof(int) == 4);
237 assert(sizeof(CliqueEntry) == 4);
238 // If lots then get rid of short ones
239 #define KEEP_CLIQUES 10000
240 if (numberCliques - numberMatrixCliques > KEEP_CLIQUES) {
241 int *sort = new int[numberCliques];
242 for (iClique = numberMatrixCliques; iClique < numberCliques; iClique++) {
243 CoinBigIndex j = cliqueStart[iClique];
244 int n = static_cast< int >(cliqueStart[iClique + 1] - j);
245 sort[iClique] = n;
246 }
247 std::sort(sort + numberMatrixCliques, sort + numberCliques);
248 int allow = sort[numberCliques - KEEP_CLIQUES];
249 int nEqual = 0;
250 for (iClique = numberCliques - KEEP_CLIQUES; iClique < numberCliques; iClique++) {
251 if (sort[iClique] > allow)
252 break;
253 else
254 nEqual++;
255 }
256 delete[] sort;
257 CoinBigIndex j = cliqueStart[numberMatrixCliques];
258 CoinBigIndex put = j;
259 int nClique = numberMatrixCliques;
260 for (iClique = numberMatrixCliques; iClique < numberCliques; iClique++) {
261 CoinBigIndex end = cliqueStart[iClique + 1];
262 int n = static_cast< int >(end - j);
263 bool copy = false;
264 if (n > allow) {
265 copy = true;
266 } else if (n == allow && nEqual) {
267 copy = true;
268 nEqual--;
269 }
270 if (copy) {
271 cliqueType[nClique++] = cliqueType[iClique];
272 for (; j < end; j++)
273 entry[put++] = entry[j];
274 }
275 j = cliqueStart[iClique + 1];
276 cliqueStart[nClique] = put;
277 }
278 numberCliques = nClique;
279 }
280 // sort
281 for (iClique = 0; iClique < numberCliques; iClique++) {
282 CoinBigIndex j = cliqueStart[iClique];
283 int n = static_cast< int >(cliqueStart[iClique + 1] - j);
284 for (int i = 0; i < n; i++)
285 whichP[i] = sequenceInCliqueEntry(entry[i + j]);
286 CoinSort_2(whichP, whichP + n, (reinterpret_cast< int * >(entry)) + j);
287 }
288 // lexicographic sort
289 int *which = new int[numberCliques];
290 CoinBigIndex *position = new CoinBigIndex[numberCliques];
291 int *sort = new int[numberCliques];
292 int *value = new int[numberCliques];
293 for (iClique = 0; iClique < numberCliques; iClique++) {
294 which[iClique] = iClique;
295 sort[iClique] = sequenceInCliqueEntry(entry[cliqueStart[iClique]]);
296 value[iClique] = sort[iClique];
297 position[iClique] = 0;
298 }
299 CoinSort_2(sort, sort + numberCliques, which);
300 int lastDone = -1;
301 int nDup = 0;
302 CoinBigIndex nSave = 0;
303 while (lastDone < numberCliques - 1) {
304 int jClique = lastDone + 1;
305 int jFirst = jClique;
306 int iFirst = which[jFirst];
307 int iValue = value[iFirst];
308 CoinBigIndex iPos = position[iFirst];
309 jClique++;
310 for (; jClique < numberCliques; jClique++) {
311 int kClique = which[jClique];
312 int jValue = value[kClique];
313 if (jValue > iValue || position[kClique] < iPos)
314 break;
315 }
316 if (jClique == jFirst + 1) {
317 // done that bit
318 lastDone++;
319 } else {
320 // use next bit to sort and then repeat
321 int jLast = jClique;
322 for (jClique = jFirst; jClique < jLast; jClique++) {
323 int kClique = which[jClique];
324 int iValue = value[kClique];
325 // put at end if finished
326 if (iValue < numberIntegers) {
327 CoinBigIndex kPos = position[kClique] + 1;
328 position[kClique] = kPos;
329 kPos += cliqueStart[kClique];
330 if (kPos == cliqueStart[kClique + 1]) {
331 iValue = numberIntegers;
332 } else {
333 iValue = sequenceInCliqueEntry(entry[kPos]);
334 }
335 value[kClique] = iValue;
336 }
337 sort[jClique] = iValue;
338 }
339 CoinSort_2(sort + jFirst, sort + jLast, which + jFirst);
340 // if duplicate mark and move on
341 int iLowest = numberCliques;
342 for (jClique = jFirst; jClique < jLast; jClique++) {
343 int kClique = which[jClique];
344 int iValue = value[kClique];
345 if (iValue < numberIntegers)
346 break;
347 iLowest = CoinMin(iLowest, kClique);
348 }
349 if (jClique > jFirst) {
350 // mark all apart from lowest number as duplicate and move on
351 lastDone = jClique - 1;
352 for (jClique = jFirst; jClique <= lastDone; jClique++) {
353 int kClique = which[jClique];
354 if (kClique != iLowest) {
355 value[kClique] = -2;
356 nDup++;
357 nSave += cliqueStart[kClique + 1] - cliqueStart[kClique];
358 }
359 }
360 }
361 }
362 }
363 if (printit)
364 printf("%d duplicates\n", nDup);
365 // Now see if any subset
366 int nOut = 0;
367 for (int jClique = 0; jClique < numberCliques; jClique++) {
368 if (value[jClique] != -2) {
369 position[jClique] = cliqueStart[jClique];
370 value[jClique] = sequenceInCliqueEntry(entry[cliqueStart[jClique]]);
371 }
372 }
373 nSave = 0;
374 int startLooking = 0;
375 for (int jClique = 0; jClique < numberCliques; jClique++) {
376 int kClique = which[jClique];
377 if (value[kClique] == -2) {
378 nOut++;
379 nSave += cliqueStart[kClique + 1] - cliqueStart[kClique];
380 if (jClique == startLooking)
381 startLooking++;
382 continue;
383 }
384 int kValue = value[kClique];
385 for (int iiClique = startLooking; iiClique < jClique; iiClique++) {
386 int iClique = which[iiClique];
387 int iValue = value[iClique];
388 if (iValue == -2 || iValue == numberIntegers) {
389 if (iiClique == startLooking)
390 startLooking++;
391 continue;
392 } else {
393 if (kValue > static_cast< int >(sequenceInCliqueEntry(entry[cliqueStart[iClique + 1] - 1]))) {
394 value[iClique] = numberIntegers;
395 continue;
396 }
397 }
398 if (iValue < kValue) {
399 while (iValue < kValue) {
400 CoinBigIndex iPos = position[iClique] + 1;
401 position[iClique] = iPos;
402 if (iPos == cliqueStart[iClique + 1]) {
403 iValue = numberIntegers;
404 } else {
405 iValue = sequenceInCliqueEntry(entry[iPos]);
406 }
407 value[iClique] = iValue;
408 }
409 }
410 if (iValue > kValue)
411 continue; // not a candidate
412 // See if subset (remember duplicates have gone)
413 if (cliqueStart[iClique + 1] - position[iClique] > cliqueStart[kClique + 1] - cliqueStart[kClique]) {
414 // could be subset ?
415 CoinBigIndex offset = cliqueStart[iClique] - position[kClique];
416 CoinBigIndex j;
417 bool subset = true;
418 // what about different fixes bool odd=false;
419 for (j = cliqueStart[kClique] + 1; j < cliqueStart[kClique + 1]; j++) {
420 int kColumn = sequenceInCliqueEntry(entry[j]);
421 int iColumn = sequenceInCliqueEntry(entry[j + offset]);
422 if (iColumn > kColumn) {
423 subset = false;
424 } else {
425 while (iColumn < kColumn) {
426 offset++;
427 if (j + offset < cliqueStart[iClique + 1]) {
428 iColumn = sequenceInCliqueEntry(entry[j + offset]);
429 } else {
430 subset = false;
431 break;
432 }
433 }
434 }
435 if (!subset)
436 break;
437 }
438 if (subset) {
439 value[kClique] = -2;
440 if (printit > 1)
441 printf("clique %d is subset of %d\n", kClique, iClique);
442 nOut++;
443 break;
444 }
445 }
446 }
447 }
448 if (nOut) {
449 if (printit)
450 printf("Can get rid of %d cliques\n", nOut);
451 // make new copy
452 int nNewC = numberCliques - nOut;
453 int size = static_cast< int >(cliqueStart[numberCliques] - nSave);
454 int n = 0;
455 CoinBigIndex *start = new CoinBigIndex[nNewC + 1];
456 char *type = new char[nNewC];
457 start[0] = 0;
458 CliqueEntry *entryC = new CliqueEntry[size];
459 CoinBigIndex nel = 0;
460 allNew = true;
461 for (int jClique = 0; jClique < numberCliques; jClique++) {
462 int kClique = which[jClique];
463 if (value[kClique] != -2 && kClique < numberMatrixCliques) {
464 if (kClique >= numberLastTime)
465 allNew = false;
466 int nn = static_cast< int >(cliqueStart[kClique + 1] - cliqueStart[kClique]);
467 memcpy(entryC + nel, entry + cliqueStart[kClique], nn * sizeof(CliqueEntry));
468 nel += nn;
469 type[n++] = cliqueType[kClique];
470 start[n] = nel;
471 }
472 }
473 int nM = n;
474 for (int jClique = 0; jClique < numberCliques; jClique++) {
475 int kClique = which[jClique];
476 if (value[kClique] != -2 && kClique >= numberMatrixCliques) {
477 if (kClique >= numberLastTime)
478 allNew = false;
479 int nn = static_cast< int >(cliqueStart[kClique + 1] - cliqueStart[kClique]);
480 memcpy(entryC + nel, entry + cliqueStart[kClique], nn * sizeof(CliqueEntry));
481 nel += nn;
482 type[n++] = cliqueType[kClique];
483 start[n] = nel;
484 }
485 }
486 // move
487 numberCliques = n;
488 numberMatrixCliques = nM;
489 delete[] cliqueStart;
490 cliqueStart = start;
491 delete[] entry;
492 entry = entryC;
493 delete[] cliqueType;
494 cliqueType = type;
495 if (printit > 1) {
496 for (int jClique = 0; jClique < numberCliques; jClique++) {
497 printf("%d [ ", jClique);
498 for (CoinBigIndex i = cliqueStart[jClique]; i < cliqueStart[jClique + 1]; i++)
499 printf("%d(%d) ", sequenceInCliqueEntry(entry[i]), oneFixesInCliqueEntry(entry[i]));
500 printf("]\n");
501 }
502 }
503 if (printit)
504 printf("%d matrix cliques and %d found by probing\n", numberMatrixCliques, numberCliques - numberMatrixCliques);
505 }
506 delete[] value;
507 delete[] sort;
508 delete[] which;
509 delete[] position;
510 delete[] whichP;
511 if (!allNew)
512 return nOut;
513 else
514 return -1;
515 }
516 OsiSolverInterface *
analyze(const OsiSolverInterface & si,int createSolver,int numberExtraCliques,const CoinBigIndex * starts,const CliqueEntry * entries,const char * type)517 CglTreeProbingInfo::analyze(const OsiSolverInterface &si, int createSolver,
518 int numberExtraCliques, const CoinBigIndex *starts,
519 const CliqueEntry *entries, const char *type)
520 {
521 if (!createSolver)
522 return NULL;
523 convert();
524 if (!numberIntegers_)
525 return NULL;
526 bool alwaysDo = false;
527 if (numberExtraCliques < 0) {
528 alwaysDo = true;
529 numberExtraCliques = 0;
530 }
531 bool printit = false;
532 int numberCliques = 0;
533 CoinBigIndex numberEntries = 0;
534 CoinBigIndex *cliqueStart = NULL;
535 CliqueEntry *entry = NULL;
536 char *cliqueType = NULL;
537 int *whichP = new int[numberIntegers_];
538 int *whichM = new int[numberIntegers_];
539 int *whichClique = NULL;
540 int numberRows = si.getNumRows();
541 int numberMatrixCliques = 0;
542 const CoinPackedMatrix *rowCopy = si.getMatrixByRow();
543 assert(numberRows && si.getNumCols());
544 int iRow;
545 const int *column = rowCopy->getIndices();
546 const double *elementByRow = rowCopy->getElements();
547 const CoinBigIndex *rowStart = rowCopy->getVectorStarts();
548 const int *rowLength = rowCopy->getVectorLengths();
549 const double *lower = si.getColLower();
550 const double *upper = si.getColUpper();
551 const double *rowLower = si.getRowLower();
552 const double *rowUpper = si.getRowUpper();
553 for (int iPass = 0; iPass < 2; iPass++) {
554 if (iPass) {
555 CoinBigIndex numberExtraEntries = 0;
556 if (numberExtraCliques)
557 numberExtraEntries = starts[numberExtraCliques];
558 cliqueStart = new CoinBigIndex[numberCliques + 1 + numberExtraCliques];
559 cliqueStart[0] = 0;
560 entry = new CliqueEntry[numberEntries + numberExtraEntries];
561 cliqueType = new char[numberCliques + numberExtraCliques];
562 whichClique = new int[numberEntries + numberExtraEntries];
563 numberCliques = 0;
564 numberEntries = 0;
565 }
566 #if 1
567 for (iRow = 0; iRow < numberRows; iRow++) {
568 int numberP1 = 0, numberM1 = 0;
569 int numberTotal = 0;
570 CoinBigIndex j;
571 double upperValue = rowUpper[iRow];
572 double lowerValue = rowLower[iRow];
573 bool good = true;
574 for (j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
575 int iColumn = column[j];
576 double value = elementByRow[j];
577 if (upper[iColumn] - lower[iColumn] < 1.0e-8) {
578 // fixed
579 upperValue -= lower[iColumn] * value;
580 lowerValue -= lower[iColumn] * value;
581 continue;
582 } else if (backward_[iColumn] < 0) {
583 good = false;
584 break;
585 } else {
586 iColumn = backward_[iColumn];
587 numberTotal++;
588 }
589 if (fabs(value) != 1.0) {
590 good = false;
591 } else if (value > 0.0) {
592 assert(numberP1 < numberIntegers_);
593 whichP[numberP1++] = iColumn;
594 ;
595 } else {
596 assert(numberM1 < numberIntegers_);
597 whichM[numberM1++] = iColumn;
598 }
599 }
600 int iUpper = static_cast< int >(floor(upperValue + 1.0e-5));
601 int iLower = static_cast< int >(ceil(lowerValue - 1.0e-5));
602 int state = 0;
603 if (upperValue < 1.0e6) {
604 if (iUpper == 1 - numberM1)
605 state = 1;
606 else if (iUpper == -numberM1)
607 state = 2;
608 else if (iUpper < -numberM1)
609 state = 3;
610 if (fabs(static_cast< double >(iUpper) - upperValue) > 1.0e-9)
611 state = -1;
612 }
613 if (!state && lowerValue > -1.0e6) {
614 if (-iLower == 1 - numberP1)
615 state = -1;
616 else if (-iLower == -numberP1)
617 state = -2;
618 else if (-iLower < -numberP1)
619 state = -3;
620 if (fabs(static_cast< double >(iLower) - lowerValue) > 1.0e-9)
621 state = -1;
622 }
623 if (numberP1 + numberM1 < 2)
624 state = -1;
625 if (good && state > 0) {
626 if (abs(state) == 3) {
627 // infeasible
628 printf("FFF Infeasible\n");
629 ;
630 //feasible=false;
631 break;
632 } else if (abs(state) == 2) {
633 // we can fix all
634 //numberFixed += numberP1+numberM1;
635 printf("FFF can fix %d\n", numberP1 + numberM1);
636 } else {
637 for (j = 0; j < numberP1; j++) {
638 int iColumn = whichP[j];
639 if (iPass) {
640 CliqueEntry temp;
641 setOneFixesInCliqueEntry(temp, true);
642 setSequenceInCliqueEntry(temp, iColumn);
643 entry[numberEntries] = temp;
644 }
645 numberEntries++;
646 }
647 for (j = 0; j < numberM1; j++) {
648 int iColumn = whichM[j];
649 if (iPass) {
650 CliqueEntry temp;
651 setOneFixesInCliqueEntry(temp, false);
652 setSequenceInCliqueEntry(temp, iColumn);
653 entry[numberEntries] = temp;
654 }
655 numberEntries++;
656 }
657 if (iPass) {
658 if (iLower != iUpper) {
659 // slack
660 cliqueType[numberCliques] = 'S';
661 } else {
662 cliqueType[numberCliques] = 'E';
663 }
664 cliqueStart[numberCliques + 1] = numberEntries;
665 }
666 numberCliques++;
667 }
668 }
669 }
670 #endif
671 if (numberExtraCliques) {
672 CoinBigIndex numberExtraEntries = starts[numberExtraCliques];
673 memcpy(entry + numberEntries, entries, numberExtraEntries * sizeof(CliqueEntry));
674 for (int iClique = 0; iClique < numberExtraCliques; iClique++) {
675 cliqueType[numberCliques] = type[iClique];
676 numberCliques++;
677 cliqueStart[numberCliques] = starts[iClique] + numberEntries;
678 }
679 numberEntries += numberExtraEntries;
680 }
681 numberMatrixCliques = numberCliques;
682 //int size = toZero_[numberIntegers_];
683 //char * used = new char [size];
684 //memset(used,0,size);
685 // find two cliques
686 int nFix = 0;
687 for (int iColumn = 0; iColumn < static_cast< int >(numberIntegers_); iColumn++) {
688 int j;
689 for (j = toZero_[iColumn]; j < toOne_[iColumn]; j++) {
690 int jColumn = sequenceInCliqueEntry(fixEntry_[j]);
691 // just look at ones beore (this also skips non 0-1)
692 if (jColumn < iColumn) {
693 int k;
694 for (k = toZero_[jColumn]; k < toOne_[jColumn]; k++) {
695 if (sequenceInCliqueEntry(fixEntry_[k]) == (iColumn)) {
696 if (oneFixesInCliqueEntry(fixEntry_[j])) {
697 if (oneFixesInCliqueEntry(fixEntry_[k])) {
698 if (printit && !iPass)
699 printf("%d to zero implies %d to one and %d to zero implies %d to one\n",
700 iColumn, jColumn, jColumn, iColumn);
701 //0-0 illegal
702 if (iPass) {
703 CliqueEntry temp;
704 setOneFixesInCliqueEntry(temp, false);
705 setSequenceInCliqueEntry(temp, iColumn);
706 entry[numberEntries] = temp;
707 }
708 numberEntries++;
709 if (iPass) {
710 CliqueEntry temp;
711 setOneFixesInCliqueEntry(temp, false);
712 setSequenceInCliqueEntry(temp, jColumn);
713 entry[numberEntries] = temp;
714 }
715 numberEntries++;
716 if (iPass) {
717 // slack
718 cliqueType[numberCliques] = 'S';
719 cliqueStart[numberCliques + 1] = numberEntries;
720 }
721 numberCliques++;
722 } else {
723 if (printit && !iPass)
724 printf("%d to zero implies %d to one and %d to zero implies %d to zero\n",
725 iColumn, jColumn, jColumn, iColumn);
726 nFix++; // jColumn is 1
727 }
728 } else {
729 if (oneFixesInCliqueEntry(fixEntry_[k])) {
730 if (printit && !iPass)
731 printf("%d to zero implies %d to zero and %d to zero implies %d to one\n",
732 iColumn, jColumn, jColumn, iColumn);
733 nFix++; // iColumn is 1
734 } else {
735 if (printit && !iPass)
736 printf("%d to zero implies %d to zero and %d to zero implies %d to zero\n",
737 iColumn, jColumn, jColumn, iColumn);
738 nFix++; // jColumn=iColumn
739 }
740 }
741 }
742 }
743 for (k = toOne_[jColumn]; k < toZero_[jColumn + 1]; k++) {
744 if (sequenceInCliqueEntry(fixEntry_[k]) == (iColumn)) {
745 if (oneFixesInCliqueEntry(fixEntry_[j])) {
746 if (oneFixesInCliqueEntry(fixEntry_[k])) {
747 if (printit && !iPass)
748 printf("%d to zero implies %d to one and %d to one implies %d to one\n",
749 iColumn, jColumn, jColumn, iColumn);
750 nFix++; //iColumn is 1
751 } else {
752 if (printit && !iPass)
753 printf("%d to zero implies %d to one and %d to one implies %d to zero\n",
754 iColumn, jColumn, jColumn, iColumn);
755 nFix++; // iColumn+jcolumn=1
756 }
757 } else {
758 if (oneFixesInCliqueEntry(fixEntry_[k])) {
759 if (printit && !iPass)
760 printf("%d to zero implies %d to zero and %d to one implies %d to one\n",
761 iColumn, jColumn, jColumn, iColumn);
762 // 0-1 illegal
763 if (iPass) {
764 CliqueEntry temp;
765 setOneFixesInCliqueEntry(temp, false);
766 setSequenceInCliqueEntry(temp, iColumn);
767 entry[numberEntries] = temp;
768 }
769 numberEntries++;
770 if (iPass) {
771 CliqueEntry temp;
772 setOneFixesInCliqueEntry(temp, true);
773 setSequenceInCliqueEntry(temp, jColumn);
774 entry[numberEntries] = temp;
775 }
776 numberEntries++;
777 if (iPass) {
778 // slack
779 cliqueType[numberCliques] = 'S';
780 cliqueStart[numberCliques + 1] = numberEntries;
781 }
782 numberCliques++;
783 } else {
784 if (printit && !iPass)
785 printf("%d to zero implies %d to zero and %d to one implies %d to zero\n",
786 iColumn, jColumn, jColumn, iColumn);
787 nFix++; // jColumn is 0
788 }
789 }
790 }
791 }
792 }
793 }
794 for (j = toOne_[iColumn]; j < toZero_[iColumn + 1]; j++) {
795 int jColumn = sequenceInCliqueEntry(fixEntry_[j]);
796 if (jColumn < iColumn) {
797 int k;
798 for (k = toZero_[jColumn]; k < toOne_[jColumn]; k++) {
799 if (sequenceInCliqueEntry(fixEntry_[k]) == (iColumn)) {
800 if (oneFixesInCliqueEntry(fixEntry_[j])) {
801 if (oneFixesInCliqueEntry(fixEntry_[k])) {
802 if (printit && !iPass)
803 printf("%d to one implies %d to one and %d to zero implies %d to one\n",
804 iColumn, jColumn, jColumn, iColumn);
805 nFix++; // jColumn is 1
806 } else {
807 if (printit && !iPass)
808 printf("%d to one implies %d to one and %d to zero implies %d to zero\n",
809 iColumn, jColumn, jColumn, iColumn);
810 // 1-0 illegal
811 if (iPass) {
812 CliqueEntry temp;
813 setOneFixesInCliqueEntry(temp, true);
814 setSequenceInCliqueEntry(temp, iColumn);
815 entry[numberEntries] = temp;
816 }
817 numberEntries++;
818 if (iPass) {
819 CliqueEntry temp;
820 setOneFixesInCliqueEntry(temp, false);
821 setSequenceInCliqueEntry(temp, jColumn);
822 entry[numberEntries] = temp;
823 }
824 numberEntries++;
825 if (iPass) {
826 // slack
827 cliqueType[numberCliques] = 'S';
828 cliqueStart[numberCliques + 1] = numberEntries;
829 }
830 numberCliques++;
831 }
832 } else {
833 if (oneFixesInCliqueEntry(fixEntry_[k])) {
834 if (printit && !iPass)
835 printf("%d to one implies %d to zero and %d to zero implies %d to one\n",
836 iColumn, jColumn, jColumn, iColumn);
837 nFix++; // iColumn+jColumn=1
838 } else {
839 if (printit && !iPass)
840 printf("%d to one implies %d to zero and %d to zero implies %d to zero\n",
841 iColumn, jColumn, jColumn, iColumn);
842 nFix++; // iColumn is 0
843 }
844 }
845 }
846 }
847 for (k = toOne_[jColumn]; k < toZero_[jColumn + 1]; k++) {
848 if (sequenceInCliqueEntry(fixEntry_[k]) == (iColumn)) {
849 if (oneFixesInCliqueEntry(fixEntry_[j])) {
850 if (oneFixesInCliqueEntry(fixEntry_[k])) {
851 if (printit && !iPass)
852 printf("%d to one implies %d to one and %d to one implies %d to one\n",
853 iColumn, jColumn, jColumn, iColumn);
854 nFix++; // iColumn == jColumn
855 } else {
856 if (printit && !iPass)
857 printf("%d to one implies %d to one and %d to one implies %d to zero\n",
858 iColumn, jColumn, jColumn, iColumn);
859 nFix++; // iColumn is 0
860 }
861 } else {
862 if (oneFixesInCliqueEntry(fixEntry_[k])) {
863 if (printit && !iPass)
864 printf("%d to one implies %d to zero and %d to one implies %d to one\n",
865 iColumn, jColumn, jColumn, iColumn);
866 nFix++; // jColumn is 0
867 } else {
868 if (printit && !iPass)
869 printf("%d to one implies %d to zero and %d to one implies %d to zero\n",
870 iColumn, jColumn, jColumn, iColumn);
871 // 1-1 illegal
872 if (iPass) {
873 CliqueEntry temp;
874 setOneFixesInCliqueEntry(temp, true);
875 setSequenceInCliqueEntry(temp, iColumn);
876 entry[numberEntries] = temp;
877 }
878 numberEntries++;
879 if (iPass) {
880 CliqueEntry temp;
881 setOneFixesInCliqueEntry(temp, true);
882 setSequenceInCliqueEntry(temp, jColumn);
883 entry[numberEntries] = temp;
884 }
885 numberEntries++;
886 if (iPass) {
887 // slack
888 cliqueType[numberCliques] = 'S';
889 cliqueStart[numberCliques + 1] = numberEntries;
890 }
891 numberCliques++;
892 }
893 }
894 }
895 }
896 }
897 }
898 }
899 if (!iPass)
900 printf("%d cliques and %d fixed (%d already from matrix))\n",
901 numberCliques - numberMatrixCliques, nFix, numberMatrixCliques);
902 }
903 int iClique;
904 outDupsEtc(numberIntegers_, numberCliques, numberMatrixCliques,
905 cliqueStart, cliqueType, entry, -1, printit ? 2 : 1);
906 printf("%d matrix cliques and %d found by probing\n", numberMatrixCliques, numberCliques - numberMatrixCliques);
907 int *zeroStart = new int[numberIntegers_ + 1];
908 int *oneStart = new int[numberIntegers_];
909 int *zeroCount = new int[numberIntegers_];
910 int *oneCount = new int[numberIntegers_];
911 char *mark = new char[numberIntegers_];
912 memset(mark, 0, numberIntegers_);
913 int nStrengthen = -1;
914 int iPass = 0;
915 while (nStrengthen && iPass < 20) {
916 iPass++;
917 int numberLastTime = numberCliques;
918 int *count = new int[numberCliques];
919 int i, iColumn;
920 for (i = 0; i < numberCliques; i++) {
921 count[i] = 0;
922 }
923 int *whichCount = new int[numberCliques];
924 CoinZeroN(zeroCount, numberIntegers_);
925 CoinZeroN(oneCount, numberIntegers_);
926 for (iClique = 0; iClique < numberCliques; iClique++) {
927 for (CoinBigIndex j = cliqueStart[iClique]; j < cliqueStart[iClique + 1]; j++) {
928 int iColumn = static_cast< int >(sequenceInCliqueEntry(entry[j]));
929 if (oneFixesInCliqueEntry(entry[j])) {
930 oneCount[iColumn]++;
931 } else {
932 zeroCount[iColumn]++;
933 }
934 }
935 }
936 int j;
937 zeroStart[0] = 0;
938 cliqueStart[0] = 0;
939 for (j = 0; j < numberIntegers_; j++) {
940 int n;
941 n = zeroCount[j];
942 zeroCount[j] = 0;
943 oneStart[j] = zeroStart[j] + n;
944 n = oneCount[j];
945 oneCount[j] = 0;
946 zeroStart[j + 1] = oneStart[j] + n;
947 }
948 for (iClique = 0; iClique < numberCliques; iClique++) {
949 for (CoinBigIndex j = cliqueStart[iClique]; j < cliqueStart[iClique + 1]; j++) {
950 int iColumn = static_cast< int >(sequenceInCliqueEntry(entry[j]));
951 if (oneFixesInCliqueEntry(entry[j])) {
952 int k = oneCount[iColumn];
953 oneCount[iColumn]++;
954 int put = oneStart[iColumn] + k;
955 whichClique[put] = iClique;
956 } else {
957 int k = zeroCount[iColumn];
958 zeroCount[iColumn]++;
959 int put = zeroStart[iColumn] + k;
960 whichClique[put] = iClique;
961 }
962 }
963 }
964 nStrengthen = 0;
965 CoinBigIndex numberEntries = cliqueStart[numberCliques];
966 CoinBigIndex maximumEntries = numberEntries;
967 int maximumCliques = numberCliques;
968 for (iColumn = 0; iColumn < numberIntegers_; iColumn++) {
969 int i;
970 int n;
971 int nCount = 0;
972 n = 0;
973 for (i = zeroStart[iColumn]; i < oneStart[iColumn]; i++) {
974 int jClique = whichClique[i];
975 //if (jClique<numberMatrixCliques)
976 //continue;
977 CoinBigIndex j = cliqueStart[jClique];
978 //assert (cliqueStart[jClique+1]==j+2);
979 for (; j < cliqueStart[jClique + 1]; j++) {
980 CliqueEntry eJ = entry[j];
981 int jColumn = sequenceInCliqueEntry(eJ);
982 if (jColumn > iColumn && !mark[jColumn]) {
983 mark[jColumn] = 1;
984 whichP[n++] = jColumn;
985 assert(n < numberIntegers_);
986 if (oneFixesInCliqueEntry(eJ)) {
987 for (int k = oneStart[jColumn]; k < zeroStart[jColumn + 1]; k++) {
988 int jClique = whichClique[k];
989 if (!count[jClique])
990 whichCount[nCount++] = jClique;
991 count[jClique]++;
992 }
993 } else {
994 for (int k = zeroStart[jColumn]; k < oneStart[jColumn]; k++) {
995 int jClique = whichClique[k];
996 if (!count[jClique])
997 whichCount[nCount++] = jClique;
998 count[jClique]++;
999 }
1000 }
1001 }
1002 }
1003 }
1004 std::sort(whichP, whichP + n);
1005 for (i = 0; i < nCount; i++) {
1006 int jClique = whichCount[i];
1007 int jCount = count[jClique];
1008 count[jClique] = 0;
1009 if (jCount == cliqueStart[jClique + 1] - cliqueStart[jClique]) {
1010 printf("Zero can extend %d [ ", jClique);
1011 for (CoinBigIndex i = cliqueStart[jClique]; i < cliqueStart[jClique + 1]; i++)
1012 printf("%d(%d) ", sequenceInCliqueEntry(entry[i]), oneFixesInCliqueEntry(entry[i]));
1013 printf("] by %d(0)\n", iColumn);
1014 nStrengthen++;
1015 if (numberEntries + jCount + 1 > maximumEntries) {
1016 maximumEntries = CoinMax(numberEntries + jCount + 1, (maximumEntries * 12) / 10 + 100);
1017 CliqueEntry *temp = new CliqueEntry[maximumEntries];
1018 memcpy(temp, entry, numberEntries * sizeof(CliqueEntry));
1019 delete[] entry;
1020 entry = temp;
1021 int *tempI = new int[maximumEntries];
1022 memcpy(tempI, whichClique, numberEntries * sizeof(int));
1023 delete[] whichClique;
1024 whichClique = tempI;
1025 }
1026 if (numberCliques == maximumCliques) {
1027 maximumCliques = (maximumCliques * 12) / 10 + 100;
1028 CoinBigIndex *temp = new CoinBigIndex[maximumCliques + 1];
1029 memcpy(temp, cliqueStart, (numberCliques + 1) * sizeof(CoinBigIndex));
1030 delete[] cliqueStart;
1031 cliqueStart = temp;
1032 char *tempT = new char[maximumCliques];
1033 memcpy(tempT, cliqueType, numberCliques);
1034 delete[] cliqueType;
1035 cliqueType = tempT;
1036 }
1037 CliqueEntry eI;
1038 eI.fixes = 0;
1039 setSequenceInCliqueEntry(eI, iColumn);
1040 setOneFixesInCliqueEntry(eI, false);
1041 entry[numberEntries++] = eI;
1042 whichM[0] = iColumn;
1043 int n = 1;
1044 for (CoinBigIndex i = cliqueStart[jClique]; i < cliqueStart[jClique + 1]; i++) {
1045 entry[numberEntries++] = entry[i];
1046 whichM[n++] = sequenceInCliqueEntry(entry[i]);
1047 }
1048 CoinSort_2(whichM, whichM + n, (reinterpret_cast< int * >(entry)) + numberEntries - n);
1049 cliqueType[numberCliques] = 'S';
1050 numberCliques++;
1051 cliqueStart[numberCliques] = numberEntries;
1052 }
1053 }
1054 for (i = 0; i < n; i++)
1055 mark[whichP[i]] = 0;
1056 nCount = 0;
1057 n = 0;
1058 for (i = oneStart[iColumn]; i < zeroStart[iColumn + 1]; i++) {
1059 int jClique = whichClique[i];
1060 //if (jClique<numberMatrixCliques)
1061 //continue;
1062 CoinBigIndex j = cliqueStart[jClique];
1063 //assert (cliqueStart[jClique+1]==j+2);
1064 for (; j < cliqueStart[jClique + 1]; j++) {
1065 CliqueEntry eJ = entry[j];
1066 int jColumn = sequenceInCliqueEntry(eJ);
1067 if (jColumn > iColumn && !mark[jColumn]) {
1068 mark[jColumn] = 1;
1069 whichP[n++] = jColumn;
1070 assert(n < numberIntegers_);
1071 if (oneFixesInCliqueEntry(eJ)) {
1072 for (int k = oneStart[jColumn]; k < zeroStart[jColumn + 1]; k++) {
1073 int jClique = whichClique[k];
1074 if (!count[jClique])
1075 whichCount[nCount++] = jClique;
1076 count[jClique]++;
1077 }
1078 } else {
1079 for (int k = zeroStart[jColumn]; k < oneStart[jColumn]; k++) {
1080 int jClique = whichClique[k];
1081 if (!count[jClique])
1082 whichCount[nCount++] = jClique;
1083 count[jClique]++;
1084 }
1085 }
1086 }
1087 }
1088 }
1089 std::sort(whichP, whichP + n);
1090 for (i = 0; i < nCount; i++) {
1091 int jClique = whichCount[i];
1092 int jCount = count[jClique];
1093 count[jClique] = 0;
1094 if (jCount == cliqueStart[jClique + 1] - cliqueStart[jClique]) {
1095 #if 1
1096 if (printit) {
1097 printf("One can extend %d [ ", jClique);
1098 for (CoinBigIndex i = cliqueStart[jClique]; i < cliqueStart[jClique + 1]; i++)
1099 printf("%d(%d) ", sequenceInCliqueEntry(entry[i]), oneFixesInCliqueEntry(entry[i]));
1100 printf("] by %d(1)\n", iColumn);
1101 }
1102 #endif
1103 nStrengthen++;
1104 if (numberEntries + jCount + 1 > maximumEntries) {
1105 maximumEntries = CoinMax(numberEntries + jCount + 1, (maximumEntries * 12) / 10 + 100);
1106 CliqueEntry *temp = new CliqueEntry[maximumEntries];
1107 memcpy(temp, entry, numberEntries * sizeof(CliqueEntry));
1108 delete[] entry;
1109 entry = temp;
1110 int *tempI = new int[maximumEntries];
1111 memcpy(tempI, whichClique, numberEntries * sizeof(int));
1112 delete[] whichClique;
1113 whichClique = tempI;
1114 }
1115 if (numberCliques == maximumCliques) {
1116 maximumCliques = (maximumCliques * 12) / 10 + 100;
1117 CoinBigIndex *temp = new CoinBigIndex[maximumCliques + 1];
1118 memcpy(temp, cliqueStart, (numberCliques + 1) * sizeof(CoinBigIndex));
1119 delete[] cliqueStart;
1120 cliqueStart = temp;
1121 char *tempT = new char[maximumCliques];
1122 memcpy(tempT, cliqueType, numberCliques);
1123 delete[] cliqueType;
1124 cliqueType = tempT;
1125 }
1126 CliqueEntry eI;
1127 eI.fixes = 0;
1128 setSequenceInCliqueEntry(eI, iColumn);
1129 setOneFixesInCliqueEntry(eI, true);
1130 entry[numberEntries++] = eI;
1131 whichM[0] = iColumn;
1132 int n = 1;
1133 for (CoinBigIndex i = cliqueStart[jClique]; i < cliqueStart[jClique + 1]; i++) {
1134 entry[numberEntries++] = entry[i];
1135 whichM[n++] = sequenceInCliqueEntry(entry[i]);
1136 }
1137 CoinSort_2(whichM, whichM + n, (reinterpret_cast< int * >(entry)) + numberEntries - n);
1138 cliqueType[numberCliques] = 'S';
1139 numberCliques++;
1140 cliqueStart[numberCliques] = numberEntries;
1141 }
1142 }
1143 for (i = 0; i < n; i++)
1144 mark[whichP[i]] = 0;
1145 }
1146 if (nStrengthen) {
1147 int numberDeleted = outDupsEtc(numberIntegers_, numberCliques, numberMatrixCliques,
1148 cliqueStart, cliqueType, entry, numberLastTime, printit ? 2 : 1);
1149 if (numberDeleted < 0 || (iPass > 1 && numberCliques - numberDeleted > 5000))
1150 nStrengthen = 0;
1151 }
1152 delete[] count;
1153 delete[] whichCount;
1154 }
1155 #if 0
1156 if (numberCliques>numberMatrixCliques) {
1157 // should keep as cliques and also use in branching ??
1158 double * element = new double [numberIntegers_];
1159 for (iClique=numberMatrixCliques;iClique<numberCliques;iClique++) {
1160 int n=0;
1161 double rhs=1.0;
1162 for (int i=cliqueStart[iClique];i<cliqueStart[iClique+1];i++) {
1163 CliqueEntry eI=entry[i];
1164 int iColumn = integerVariable_[sequenceInCliqueEntry(eI)];
1165 whichP[n]=iColumn;
1166 if (oneFixesInCliqueEntry(eI)) {
1167 element[n++]=1.0;
1168 } else {
1169 element[n++]=-1.0;
1170 rhs -= 1.0;
1171 }
1172 }
1173 stored->addCut(-COIN_DBL_MAX,rhs,n,whichP,element);
1174 }
1175 delete [] element;
1176 }
1177 #endif
1178 OsiSolverInterface *newSolver = NULL;
1179 if (numberCliques > numberMatrixCliques || alwaysDo) {
1180 newSolver = si.clone();
1181 // Delete all rows
1182 CoinBigIndex *start = new CoinBigIndex[CoinMax(numberRows, numberCliques + 1)];
1183 int i;
1184 int *start2 = reinterpret_cast< int * >(start);
1185 for (i = 0; i < numberRows; i++)
1186 start2[i] = i;
1187 newSolver->deleteRows(numberRows, start2);
1188 start[0] = 0;
1189 CoinBigIndex numberElements = cliqueStart[numberCliques];
1190 int *column = new int[numberElements];
1191 double *element = new double[numberElements];
1192 double *lower = new double[numberCliques];
1193 double *upper = new double[numberCliques];
1194 numberElements = 0;
1195 for (iClique = 0; iClique < numberCliques; iClique++) {
1196 double rhs = 1.0;
1197 for (CoinBigIndex i = cliqueStart[iClique]; i < cliqueStart[iClique + 1]; i++) {
1198 CliqueEntry eI = entry[i];
1199 int iColumn = integerVariable_[sequenceInCliqueEntry(eI)];
1200 column[numberElements] = iColumn;
1201 if (oneFixesInCliqueEntry(eI)) {
1202 element[numberElements++] = 1.0;
1203 } else {
1204 element[numberElements++] = -1.0;
1205 rhs -= 1.0;
1206 }
1207 }
1208 start[iClique + 1] = numberElements;
1209 assert(cliqueType[iClique] == 'S' || cliqueType[iClique] == 'E');
1210 if (cliqueType[iClique] == 'S')
1211 lower[iClique] = -COIN_DBL_MAX;
1212 else
1213 lower[iClique] = rhs;
1214 upper[iClique] = rhs;
1215 }
1216 newSolver->addRows(numberCliques, start, column, element, lower, upper);
1217 delete[] start;
1218 delete[] column;
1219 delete[] element;
1220 delete[] lower;
1221 delete[] upper;
1222 }
1223 delete[] mark;
1224 delete[] whichP;
1225 delete[] whichM;
1226 delete[] cliqueStart;
1227 delete[] entry;
1228 delete[] cliqueType;
1229 delete[] zeroStart;
1230 delete[] oneStart;
1231 delete[] zeroCount;
1232 delete[] oneCount;
1233 delete[] whichClique;
1234 return newSolver;
1235 }
1236 // Take action if cut generator can fix a variable (toValue -1 for down, +1 for up)
fixes(int variable,int toValue,int fixedVariable,bool fixedToLower)1237 bool CglTreeProbingInfo::fixes(int variable, int toValue, int fixedVariable, bool fixedToLower)
1238 {
1239 //printf("%d going to %d fixes %d at %g\n",variable,toValue,fixedVariable,fixedToValue);
1240 int intVariable = backward_[variable];
1241 if (intVariable < 0) // off as no longer in order FIX
1242 return true; // not 0-1 (well wasn't when constructor was called)
1243 int intFix = backward_[fixedVariable];
1244 if (intFix < 0)
1245 intFix = numberIntegers_ + fixedVariable; // not 0-1
1246 int fixedTo = fixedToLower ? 0 : 1;
1247 if (numberEntries_ == maximumEntries_) {
1248 // See if taking too much memory
1249 if (maximumEntries_ >= CoinMax(1000000, 10 * numberIntegers_))
1250 return false;
1251 maximumEntries_ += 100 + maximumEntries_ / 2;
1252 CliqueEntry *temp1 = new CliqueEntry[maximumEntries_];
1253 memcpy(temp1, fixEntry_, numberEntries_ * sizeof(CliqueEntry));
1254 delete[] fixEntry_;
1255 fixEntry_ = temp1;
1256 int *temp2 = new int[maximumEntries_];
1257 memcpy(temp2, fixingEntry_, numberEntries_ * sizeof(int));
1258 delete[] fixingEntry_;
1259 fixingEntry_ = temp2;
1260 }
1261 CliqueEntry entry1;
1262 entry1.fixes = 0;
1263 setOneFixesInCliqueEntry(entry1, fixedTo != 0);
1264 setSequenceInCliqueEntry(entry1, intFix);
1265 fixEntry_[numberEntries_] = entry1;
1266 assert(toValue == -1 || toValue == 1);
1267 assert(fixedTo == 0 || fixedTo == 1);
1268 if (toValue < 0)
1269 fixingEntry_[numberEntries_++] = intVariable << 1;
1270 else
1271 fixingEntry_[numberEntries_++] = (intVariable << 1) | 1;
1272 return true;
1273 }
1274 // Initalizes fixing arrays etc - returns true if we want to save info
initializeFixing(const OsiSolverInterface * model)1275 int CglTreeProbingInfo::initializeFixing(const OsiSolverInterface *model)
1276 {
1277 if (numberEntries_ >= 0)
1278 return 2; // already got arrays
1279 else if (numberEntries_ == -2)
1280 return numberEntries_;
1281 delete[] fixEntry_;
1282 delete[] toZero_;
1283 delete[] toOne_;
1284 delete[] integerVariable_;
1285 delete[] backward_;
1286 delete[] fixingEntry_;
1287 numberVariables_ = model->getNumCols();
1288 // Too many ... but
1289 integerVariable_ = new int[numberVariables_];
1290 backward_ = new int[numberVariables_];
1291 numberIntegers_ = 0;
1292 int i;
1293 // Get integer types
1294 const char *columnType = model->getColType(true);
1295 for (i = 0; i < numberVariables_; i++) {
1296 backward_[i] = -1;
1297 if (columnType[i]) {
1298 if (columnType[i] == 1) {
1299 backward_[i] = numberIntegers_;
1300 integerVariable_[numberIntegers_++] = i;
1301 } else {
1302 backward_[i] = -2;
1303 }
1304 }
1305 }
1306 toZero_ = NULL;
1307 toOne_ = NULL;
1308 fixEntry_ = NULL;
1309 fixingEntry_ = NULL;
1310 maximumEntries_ = 0;
1311 numberEntries_ = 0;
1312 return 1;
1313 }
1314 // Converts to ordered and takes out duplicates
convert()1315 void CglTreeProbingInfo::convert()
1316 {
1317 if (numberEntries_ >= 0) {
1318 CoinSort_2(fixingEntry_, fixingEntry_ + numberEntries_, fixEntry_);
1319 assert(!toZero_);
1320 toZero_ = new int[numberIntegers_ + 1];
1321 toOne_ = new int[numberIntegers_];
1322 toZero_[0] = 0;
1323 int n = 0;
1324 int put = 0;
1325 for (int intVariable = 0; intVariable < numberIntegers_; intVariable++) {
1326 int last = n;
1327 for (; n < numberEntries_; n++) {
1328 int value = fixingEntry_[n];
1329 int iVar = value >> 1;
1330 int way = value & 1;
1331 if (intVariable != iVar || way)
1332 break;
1333 }
1334 if (n > last) {
1335 // sort
1336 assert(sizeof(int) == 4);
1337 std::sort(reinterpret_cast< unsigned int * >(fixEntry_) + last,
1338 reinterpret_cast< unsigned int * >(fixEntry_) + n);
1339 CliqueEntry temp2;
1340 temp2.fixes = 0;
1341 setSequenceInCliqueEntry(temp2, numberVariables_ + 1);
1342 for (int i = last; i < n; i++) {
1343 if (sequenceInCliqueEntry(temp2) != sequenceInCliqueEntry(fixEntry_[i]) || oneFixesInCliqueEntry(temp2) || oneFixesInCliqueEntry(fixEntry_[i])) {
1344 temp2 = fixEntry_[i];
1345 fixEntry_[put++] = temp2;
1346 }
1347 }
1348 }
1349 toOne_[intVariable] = put;
1350 last = n;
1351 for (; n < numberEntries_; n++) {
1352 int value = fixingEntry_[n];
1353 int iVar = value >> 1;
1354 if (intVariable != iVar)
1355 break;
1356 }
1357 if (n > last) {
1358 // sort
1359 assert(sizeof(int) == 4);
1360 std::sort(reinterpret_cast< unsigned int * >(fixEntry_) + last,
1361 reinterpret_cast< unsigned int * >(fixEntry_) + n);
1362 CliqueEntry temp2;
1363 temp2.fixes = 0;
1364 setSequenceInCliqueEntry(temp2, numberVariables_ + 1);
1365 for (int i = last; i < n; i++) {
1366 if (sequenceInCliqueEntry(temp2) != sequenceInCliqueEntry(fixEntry_[i]) || oneFixesInCliqueEntry(temp2) || oneFixesInCliqueEntry(fixEntry_[i])) {
1367 temp2 = fixEntry_[i];
1368 fixEntry_[put++] = temp2;
1369 }
1370 }
1371 last = n;
1372 }
1373 toZero_[intVariable + 1] = put;
1374 }
1375 delete[] fixingEntry_;
1376 fixingEntry_ = NULL;
1377 numberEntries_ = -2;
1378 }
1379 }
1380 // Fix entries in a solver using implications
fixColumns(OsiSolverInterface & si) const1381 int CglTreeProbingInfo::fixColumns(OsiSolverInterface &si) const
1382 {
1383 int nFix = 0;
1384 const double *lower = si.getColLower();
1385 const double *upper = si.getColUpper();
1386 bool feasible = true;
1387 for (int jColumn = 0; jColumn < static_cast< int >(numberIntegers_); jColumn++) {
1388 int iColumn = integerVariable_[jColumn];
1389 if (upper[iColumn] == 0.0) {
1390 int j;
1391 for (j = toZero_[jColumn]; j < toOne_[jColumn]; j++) {
1392 int kColumn = sequenceInCliqueEntry(fixEntry_[j]);
1393 kColumn = integerVariable_[kColumn];
1394 bool fixToOne = oneFixesInCliqueEntry(fixEntry_[j]);
1395 if (fixToOne) {
1396 if (lower[kColumn] == 0.0) {
1397 if (upper[kColumn] == 1.0) {
1398 si.setColLower(kColumn, 1.0);
1399 nFix++;
1400 } else {
1401 // infeasible!
1402 feasible = false;
1403 }
1404 }
1405 } else {
1406 if (upper[kColumn] == 1.0) {
1407 if (lower[kColumn] == 0.0) {
1408 si.setColUpper(kColumn, 0.0);
1409 nFix++;
1410 } else {
1411 // infeasible!
1412 feasible = false;
1413 }
1414 }
1415 }
1416 }
1417 } else if (lower[iColumn] == 1.0) {
1418 int j;
1419 for (j = toOne_[jColumn]; j < toZero_[jColumn + 1]; j++) {
1420 int kColumn = sequenceInCliqueEntry(fixEntry_[j]);
1421 kColumn = integerVariable_[kColumn];
1422 bool fixToOne = oneFixesInCliqueEntry(fixEntry_[j]);
1423 if (fixToOne) {
1424 if (lower[kColumn] == 0.0) {
1425 if (upper[kColumn] == 1.0) {
1426 si.setColLower(kColumn, 1.0);
1427 nFix++;
1428 } else {
1429 // infeasible!
1430 feasible = false;
1431 }
1432 }
1433 } else {
1434 if (upper[kColumn] == 1.0) {
1435 if (lower[kColumn] == 0.0) {
1436 si.setColUpper(kColumn, 0.0);
1437 nFix++;
1438 } else {
1439 // infeasible!
1440 feasible = false;
1441 }
1442 }
1443 }
1444 }
1445 }
1446 }
1447 if (!feasible) {
1448 #ifdef COIN_DEVELOP
1449 printf("treeprobing says infeasible!\n");
1450 #endif
1451 nFix = -1;
1452 }
1453 return nFix;
1454 }
1455 // Fix entries in a solver using implications for one variable
fixColumns(int iColumn,int value,OsiSolverInterface & si) const1456 int CglTreeProbingInfo::fixColumns(int iColumn, int value, OsiSolverInterface &si) const
1457 {
1458 assert(value == 0 || value == 1);
1459 int nFix = 0;
1460 const double *lower = si.getColLower();
1461 const double *upper = si.getColUpper();
1462 bool feasible = true;
1463 int jColumn = backward_[iColumn];
1464 if (jColumn < 0 || !toZero_)
1465 return 0;
1466 if (!value) {
1467 int j;
1468 for (j = toZero_[jColumn]; j < toOne_[jColumn]; j++) {
1469 int kColumn = sequenceInCliqueEntry(fixEntry_[j]);
1470 kColumn = integerVariable_[kColumn];
1471 bool fixToOne = oneFixesInCliqueEntry(fixEntry_[j]);
1472 if (fixToOne) {
1473 if (lower[kColumn] == 0.0) {
1474 if (upper[kColumn] == 1.0) {
1475 si.setColLower(kColumn, 1.0);
1476 nFix++;
1477 } else {
1478 // infeasible!
1479 feasible = false;
1480 }
1481 }
1482 } else {
1483 if (upper[kColumn] == 1.0) {
1484 if (lower[kColumn] == 0.0) {
1485 si.setColUpper(kColumn, 0.0);
1486 nFix++;
1487 } else {
1488 // infeasible!
1489 feasible = false;
1490 }
1491 }
1492 }
1493 }
1494 } else {
1495 int j;
1496 for (j = toOne_[jColumn]; j < toZero_[jColumn + 1]; j++) {
1497 int kColumn = sequenceInCliqueEntry(fixEntry_[j]);
1498 kColumn = integerVariable_[kColumn];
1499 bool fixToOne = oneFixesInCliqueEntry(fixEntry_[j]);
1500 if (fixToOne) {
1501 if (lower[kColumn] == 0.0) {
1502 if (upper[kColumn] == 1.0) {
1503 si.setColLower(kColumn, 1.0);
1504 nFix++;
1505 } else {
1506 // infeasible!
1507 feasible = false;
1508 }
1509 }
1510 } else {
1511 if (upper[kColumn] == 1.0) {
1512 if (lower[kColumn] == 0.0) {
1513 si.setColUpper(kColumn, 0.0);
1514 nFix++;
1515 } else {
1516 // infeasible!
1517 feasible = false;
1518 }
1519 }
1520 }
1521 }
1522 }
1523 if (!feasible) {
1524 #ifdef COIN_DEVELOP
1525 printf("treeprobing says infeasible!\n");
1526 #endif
1527 nFix = -1;
1528 }
1529 return nFix;
1530 }
1531 // Packs down entries
packDown()1532 int CglTreeProbingInfo::packDown()
1533 {
1534 convert();
1535 int iPut = 0;
1536 int iLast = 0;
1537 for (int jColumn = 0; jColumn < static_cast< int >(numberIntegers_); jColumn++) {
1538 int j;
1539 for (j = iLast; j < toOne_[jColumn]; j++) {
1540 int kColumn = sequenceInCliqueEntry(fixEntry_[j]);
1541 if (kColumn < numberIntegers_)
1542 fixEntry_[iPut++] = fixEntry_[j];
1543 }
1544 iLast = toOne_[jColumn];
1545 toOne_[jColumn] = iPut;
1546 for (j = iLast; j < toZero_[jColumn + 1]; j++) {
1547 int kColumn = sequenceInCliqueEntry(fixEntry_[j]);
1548 if (kColumn < numberIntegers_)
1549 fixEntry_[iPut++] = fixEntry_[j];
1550 }
1551 iLast = toZero_[jColumn + 1];
1552 toZero_[jColumn + 1] = iPut;
1553 }
1554 return iPut;
1555 }
generateCuts(const OsiSolverInterface & si,OsiCuts & cs,const CglTreeInfo) const1556 void CglTreeProbingInfo::generateCuts(const OsiSolverInterface &si, OsiCuts &cs,
1557 const CglTreeInfo /*info*/) const
1558 {
1559 const double *lower = si.getColLower();
1560 const double *upper = si.getColUpper();
1561 const double *colsol = si.getColSolution();
1562 CoinPackedVector lbs(false);
1563 CoinPackedVector ubs(false);
1564 int numberFixed = 0;
1565 char *fixed = NULL;
1566 for (int jColumn = 0; jColumn < static_cast< int >(numberIntegers_); jColumn++) {
1567 int iColumn = integerVariable_[jColumn];
1568 assert(iColumn >= 0 && iColumn < si.getNumCols());
1569 if (lower[iColumn] == 0.0 && upper[iColumn] == 1.0) {
1570 double value1 = colsol[iColumn];
1571 int j;
1572 for (j = toZero_[jColumn]; j < toOne_[jColumn]; j++) {
1573 int kColumn = sequenceInCliqueEntry(fixEntry_[j]);
1574 kColumn = integerVariable_[kColumn];
1575 assert(kColumn >= 0 && kColumn < si.getNumCols());
1576 assert(kColumn != iColumn);
1577 if (lower[kColumn] == 0.0 && upper[kColumn] == 1.0) {
1578 double value2 = colsol[kColumn];
1579 bool fixToOne = oneFixesInCliqueEntry(fixEntry_[j]);
1580 if (fixToOne) {
1581 if (value1 + value2 < 0.99999) {
1582 OsiRowCut rc;
1583 int index[2];
1584 double element[2];
1585 index[0] = iColumn;
1586 element[0] = 1.0;
1587 index[1] = kColumn;
1588 element[1] = 1.0;
1589 rc.setLb(1.0);
1590 rc.setUb(COIN_DBL_MAX);
1591 rc.setRow(2, index, element, false);
1592 cs.insertIfNotDuplicate(rc);
1593 }
1594 } else {
1595 if (value1 - value2 < -0.00001) {
1596 OsiRowCut rc;
1597 int index[2];
1598 double element[2];
1599 index[0] = iColumn;
1600 element[0] = 1.0;
1601 index[1] = kColumn;
1602 element[1] = -1.0;
1603 rc.setLb(0.0);
1604 rc.setUb(COIN_DBL_MAX);
1605 rc.setRow(2, index, element, false);
1606 cs.insertIfNotDuplicate(rc);
1607 }
1608 }
1609 }
1610 }
1611 for (j = toOne_[jColumn]; j < toZero_[jColumn + 1]; j++) {
1612 int kColumn = sequenceInCliqueEntry(fixEntry_[j]);
1613 kColumn = integerVariable_[kColumn];
1614 assert(kColumn >= 0 && kColumn < si.getNumCols());
1615 assert(kColumn != iColumn);
1616 if (lower[kColumn] == 0.0 && upper[kColumn] == 1.0) {
1617 double value2 = colsol[kColumn];
1618 bool fixToOne = oneFixesInCliqueEntry(fixEntry_[j]);
1619 if (fixToOne) {
1620 if (value1 - value2 > 0.00001) {
1621 OsiRowCut rc;
1622 int index[2];
1623 double element[2];
1624 index[0] = iColumn;
1625 element[0] = 1.0;
1626 index[1] = kColumn;
1627 element[1] = -1.0;
1628 rc.setLb(-COIN_DBL_MAX);
1629 rc.setUb(0.0);
1630 rc.setRow(2, index, element, false);
1631 cs.insertIfNotDuplicate(rc);
1632 }
1633 } else {
1634 if (value1 + value2 > 1.00001) {
1635 OsiRowCut rc;
1636 int index[2];
1637 double element[2];
1638 index[0] = iColumn;
1639 element[0] = 1.0;
1640 index[1] = kColumn;
1641 element[1] = 1.0;
1642 rc.setLb(-COIN_DBL_MAX);
1643 rc.setUb(1.0);
1644 rc.setRow(2, index, element, false);
1645 cs.insertIfNotDuplicate(rc);
1646 }
1647 }
1648 }
1649 }
1650 } else if (upper[iColumn] == 0.0) {
1651 for (int j = toZero_[jColumn]; j < toOne_[jColumn]; j++) {
1652 int kColumn01 = sequenceInCliqueEntry(fixEntry_[j]);
1653 int kColumn = integerVariable_[kColumn01];
1654 assert(kColumn >= 0 && kColumn < si.getNumCols());
1655 bool fixToOne = oneFixesInCliqueEntry(fixEntry_[j]);
1656 if (lower[kColumn] == 0.0 && upper[kColumn] == 1.0) {
1657 if (!fixed) {
1658 fixed = new char[numberIntegers_];
1659 memset(fixed, 0, numberIntegers_);
1660 }
1661 if (fixToOne) {
1662 if ((fixed[kColumn01] & 1) == 0) {
1663 fixed[kColumn01] |= 1;
1664 lbs.insert(kColumn, 1.0);
1665 }
1666 } else {
1667 if ((fixed[kColumn01] & 2) == 0) {
1668 fixed[kColumn01] |= 2;
1669 ubs.insert(kColumn, 0.0);
1670 }
1671 }
1672 numberFixed++;
1673 } else if ((fixToOne && upper[kColumn] == 0.0) || (!fixToOne && lower[kColumn] == 1.0)) {
1674 // infeasible!
1675 // generate infeasible cut and return
1676 OsiRowCut rc;
1677 rc.setLb(COIN_DBL_MAX);
1678 rc.setUb(0.0);
1679 cs.insertIfNotDuplicate(rc);
1680 //printf("IMPINFEAS!\n");
1681 return;
1682 }
1683 }
1684 } else {
1685 for (int j = toOne_[jColumn]; j < toZero_[jColumn + 1]; j++) {
1686 int kColumn01 = sequenceInCliqueEntry(fixEntry_[j]);
1687 int kColumn = integerVariable_[kColumn01];
1688 assert(kColumn >= 0 && kColumn < si.getNumCols());
1689 bool fixToOne = oneFixesInCliqueEntry(fixEntry_[j]);
1690 if (lower[kColumn] == 0.0 && upper[kColumn] == 1.0) {
1691 if (!fixed) {
1692 fixed = new char[numberIntegers_];
1693 memset(fixed, 0, numberIntegers_);
1694 }
1695 if (fixToOne) {
1696 if ((fixed[kColumn01] & 1) == 0) {
1697 fixed[kColumn01] |= 1;
1698 lbs.insert(kColumn, 1.0);
1699 }
1700 } else {
1701 if ((fixed[kColumn01] & 2) == 0) {
1702 fixed[kColumn01] |= 2;
1703 ubs.insert(kColumn, 0.0);
1704 }
1705 }
1706 numberFixed++;
1707 } else if ((fixToOne && upper[kColumn] == 0.0) || (!fixToOne && lower[kColumn] == 1.0)) {
1708 // infeasible!
1709 // generate infeasible cut and return
1710 OsiRowCut rc;
1711 rc.setLb(COIN_DBL_MAX);
1712 rc.setUb(0.0);
1713 cs.insertIfNotDuplicate(rc);
1714 //printf("IMPINFEAS!\n");
1715 return;
1716 }
1717 }
1718 }
1719 }
1720 if (numberFixed) {
1721 int feasible = true;
1722 for (int i = 0; i < numberIntegers_; i++) {
1723 if (fixed[i] == 3) {
1724 feasible = false;
1725 break;
1726 }
1727 }
1728 delete[] fixed;
1729 fixed = NULL;
1730 if (feasible) {
1731 //printf("IMP fixed %d\n",numberFixed);
1732 OsiColCut cc;
1733 cc.setUbs(ubs);
1734 cc.setLbs(lbs);
1735 cc.setEffectiveness(1.0);
1736 cs.insert(cc);
1737 } else {
1738 // infeasible!
1739 // generate infeasible cut
1740 OsiRowCut rc;
1741 rc.setLb(COIN_DBL_MAX);
1742 rc.setUb(0.0);
1743 cs.insertIfNotDuplicate(rc);
1744 //printf("IMPINFEAS!\n");
1745 }
1746 }
1747
1748 if (fixed)
1749 delete[] fixed;
1750 }
1751
1752 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
1753 */
1754