1 /*===========================================================================*
2 * This file is part of the Abstract Library for Parallel Search (ALPS). *
3 * *
4 * ALPS is distributed under the Eclipse Public License as part of the *
5 * COIN-OR repository (http://www.coin-or.org). *
6 * *
7 * Authors: *
8 * *
9 * Yan Xu, Lehigh University *
10 * Aykut Bulut, Lehigh University *
11 * Ted Ralphs, Lehigh University *
12 * *
13 * Conceptual Design: *
14 * *
15 * Yan Xu, Lehigh University *
16 * Ted Ralphs, Lehigh University *
17 * Laszlo Ladanyi, IBM T.J. Watson Research Center *
18 * Matthew Saltzman, Clemson University *
19 * *
20 * *
21 * Copyright (C) 2001-2019, Lehigh University, Yan Xu, Aykut Bulut, and *
22 * Ted Ralphs. *
23 * All Rights Reserved. *
24 *===========================================================================*/
25
26
27 //#############################################################################
28 // This file is modified from SbbModel.cpp
29 //#############################################################################
30
31 #include <iostream>
32
33 #include "OsiRowCut.hpp"
34 #include "OsiColCut.hpp"
35 #include "OsiRowCutDebugger.hpp"
36
37 #include "AlpsTreeNode.h"
38
39 #include "AbcCutGenerator.h"
40 #include "AbcHeuristic.h"
41 #include "AbcMessage.h"
42 #include "AbcModel.h"
43 #include "AbcTreeNode.h"
44 #include "AbcNodeDesc.h"
45
46
47 //#############################################################################
48 //#############################################################################
49
50 void
initialSolve()51 AbcModel::initialSolve()
52 {
53 assert (solver_);
54 solver_->initialSolve();
55 }
56
57 //#############################################################################
58 // Parameters:
59 // cuts: (o) all cuts generated in this round of cut generation
60 // numberTries: (i) the maximum number of iterations for this round of cut
61 // generation; no a priori limit if 0
62 // whichGenerator: (i/o) whichGenerator[i] is loaded with the index of the
63 // generator that produced cuts[i]; reallocated as
64 // required
65 // numberOldActiveCuts: (o) the number of active cuts at this node from
66 // previous rounds of cut generation
67 // numberNewCuts: (o) the number of cuts produced in this round of cut
68 // generation
69 // maximumWhich: (i/o) capacity of whichGenerator; may be updated if
70 // whichGenerator grows.
71 // cutDuringRampup: (i) Whether generating cuts during rampup
72 // found: (o) great than 0 means that heuristics found solutions;
73 // otherwise not.
74 bool
solveWithCuts(OsiCuts & cuts,int numberTries,AbcTreeNode * node,int & numberOldActiveCuts,int & numberNewCuts,int & maximumWhich,int * & whichGenerator,const bool cutDuringRampup,int & found)75 AbcModel::solveWithCuts(OsiCuts & cuts, int numberTries,
76 AbcTreeNode * node, int & numberOldActiveCuts,
77 int & numberNewCuts, int & maximumWhich,
78 int *& whichGenerator, const bool cutDuringRampup,
79 int & found)
80 {
81 found = -10;
82 bool feasible;
83 int lastNumberCuts = 0;
84 double lastObjective = -1.0e100 ;
85 int violated = 0;
86 int numberRowsAtStart = solver_->getNumRows();
87 int numberColumns = solver_->getNumCols();
88
89 numberOldActiveCuts = numberRowsAtStart - numberRowsAtContinuous_;
90 numberNewCuts = 0;
91
92 feasible = resolve();
93 if(!feasible) {
94 return false; // If lost feasibility, bail out right now
95 }
96
97 reducedCostFix();
98 const double *lower = solver_->getColLower();
99 const double *upper = solver_->getColUpper();
100 const double *solution = solver_->getColSolution();
101
102 double minimumDrop = minimumDrop_;
103 if (numberTries < 0) {
104 numberTries = -numberTries;
105 minimumDrop = -1.0;
106 }
107
108 //-------------------------------------------------------------------------
109 // Is it time to scan the cuts in order to remove redundant cuts? If so,
110 // set up to do it.
111 # define SCANCUTS 100
112 int *countColumnCuts = NULL;
113 int *countRowCuts = NULL;
114 bool fullScan = false;
115 if ((numberNodes_ % SCANCUTS) == 0) {
116 fullScan = true;
117 countColumnCuts = new int[numberCutGenerators_ + numberHeuristics_];
118 countRowCuts = new int[numberCutGenerators_ + numberHeuristics_];
119 memset(countColumnCuts, 0,
120 (numberCutGenerators_ + numberHeuristics_) * sizeof(int));
121 memset(countRowCuts, 0,
122 (numberCutGenerators_ + numberHeuristics_) * sizeof(int));
123 }
124
125 double direction = solver_->getObjSense();
126 double startObjective = solver_->getObjValue() * direction;
127
128 int numberPasses = 0;
129 double primalTolerance = 1.0e-7;
130
131 //-------------------------------------------------------------------------
132 // Start cut generation loop
133 do {
134 numberPasses++;
135 numberTries--;
136 OsiCuts theseCuts;
137
138 // First check if there are cuts violated in global cut pool
139 if (numberPasses == 1 && howOftenGlobalScan_ > 0 &&
140 (numberNodes_ % howOftenGlobalScan_) == 0) {
141 int numberCuts = globalCuts_.sizeColCuts();
142 int i;
143 for ( i = 0; i < numberCuts; ++i) {
144 const OsiColCut *thisCut = globalCuts_.colCutPtr(i);
145 if (thisCut->violated(solution) > primalTolerance) {
146 printf("Global cut added - violation %g\n",
147 thisCut->violated(solution));
148 theseCuts.insert(*thisCut);
149 }
150 }
151 numberCuts = globalCuts_.sizeRowCuts();
152 for ( i = 0; i < numberCuts; ++i) {
153 const OsiRowCut * thisCut = globalCuts_.rowCutPtr(i);
154 if (thisCut->violated(solution) > primalTolerance) {
155 printf("Global cut added - violation %g\n",
156 thisCut->violated(solution));
157 theseCuts.insert(*thisCut);
158 }
159 }
160 }
161
162 //---------------------------------------------------------------------
163 // Generate new cuts (global and/or local) and/or apply heuristics
164 // NOTE: Make sure CglProbing is added FIRST
165 double * newSolution = new double [numberColumns];
166 double heuristicValue = getCutoff();
167
168 #if defined(ABC_DEBUG_MORE)
169 std::cout << "numberCutGenerators_ = " << numberCutGenerators_
170 << "numberHeuristics_ = " << numberHeuristics_
171 << std::endl;
172 #endif
173 for (int i = 0; i < numberCutGenerators_ + numberHeuristics_; ++i) {
174 int numberRowCutsBefore = theseCuts.sizeRowCuts();
175 int numberColumnCutsBefore = theseCuts.sizeColCuts();
176 if (i < numberCutGenerators_) {
177 if (cutDuringRampup) {
178 bool mustResolve =
179 generator_[i]->generateCuts(theseCuts, fullScan);
180 if (mustResolve) {
181 feasible = resolve();
182 if (!feasible)
183 break;
184 }
185 }
186 }
187 else {
188 double saveValue = heuristicValue;
189 int ifSol = heuristic_[i-numberCutGenerators_]->
190 solution(heuristicValue, newSolution);
191 // solution(heuristicValue, newSolution, theseCuts);
192
193 if (ifSol > 0) {
194 found = i;
195 }
196 else if (ifSol < 0) {
197 heuristicValue = saveValue;
198 }
199 }
200 int numberRowCutsAfter = theseCuts.sizeRowCuts();
201 int numberColumnCutsAfter = theseCuts.sizeColCuts();
202 int numberBefore =
203 numberRowCutsBefore + numberColumnCutsBefore + lastNumberCuts;
204 int numberAfter =
205 numberRowCutsAfter + numberColumnCutsAfter + lastNumberCuts;
206 if (numberAfter > maximumWhich) {
207 maximumWhich = std::max(maximumWhich * 2 + 100, numberAfter);
208 int * temp = new int[2 * maximumWhich];
209 memcpy(temp, whichGenerator, numberBefore * sizeof(int));
210 delete [] whichGenerator;
211 whichGenerator = temp;
212 }
213 int j;
214 if (fullScan) {
215 countRowCuts[i] += numberRowCutsAfter -
216 numberRowCutsBefore;
217 countColumnCuts[i] += numberColumnCutsAfter -
218 numberColumnCutsBefore;
219 }
220 for (j = numberRowCutsBefore; j < numberRowCutsAfter; ++j) {
221 whichGenerator[numberBefore++] = i;
222 const OsiRowCut * thisCut = theseCuts.rowCutPtr(j);
223 if (thisCut->globallyValid()) {
224 globalCuts_.insert(*thisCut);
225 }
226 }
227 for (j = numberColumnCutsBefore; j < numberColumnCutsAfter; ++j) {
228 whichGenerator[numberBefore++] = i;
229 const OsiColCut * thisCut = theseCuts.colCutPtr(j);
230 if (thisCut->globallyValid()) {
231 globalCuts_.insert(*thisCut);
232 }
233 }
234 }
235
236 //---------------------------------------------------------------------
237 // If found a solution, Record it before we free the vector
238 if (found >= 0) {
239 bool better =
240 setBestSolution(ABC_ROUNDING, heuristicValue, newSolution);
241 // if (!better){
242 // found = -1;
243 //}
244 //std::cout << "better = " << better
245 // << "; found = " << found << std::endl;
246 }
247 if(newSolution != 0) delete [] newSolution;
248
249 int numberColumnCuts = theseCuts.sizeColCuts();
250 int numberRowCuts = theseCuts.sizeRowCuts();
251 violated = numberRowCuts + numberColumnCuts;
252
253 //---------------------------------------------------------------------
254 // Apply column cuts
255 if (numberColumnCuts) {
256 double integerTolerance = getDblParam(AbcIntegerTolerance);
257 for (int i = 0; i < numberColumnCuts; ++i) {
258 const OsiColCut * thisCut = theseCuts.colCutPtr(i);
259 const CoinPackedVector & lbs = thisCut->lbs();
260 const CoinPackedVector & ubs = thisCut->ubs();
261 int j;
262 int n;
263 const int * which;
264 const double * values;
265 n = lbs.getNumElements();
266 which = lbs.getIndices();
267 values = lbs.getElements();
268 for (j = 0; j < n; ++j){
269 int iColumn = which[j];
270 double value = solution[iColumn];
271 solver_->setColLower(iColumn, values[j]);
272 if (value < values[j] - integerTolerance)
273 violated = -1; // violated, TODO: when happen?
274 if (values[j] > upper[iColumn] + integerTolerance) {
275 violated = -2; // infeasible
276 break;
277 }
278 }
279 n = ubs.getNumElements();
280 which = ubs.getIndices();
281 values = ubs.getElements();
282 for (j = 0; j < n; ++j) {
283 int iColumn = which[j];
284 double value = solution[iColumn];
285 solver_->setColUpper(iColumn, values[j]);
286 if (value > values[j] + integerTolerance)
287 violated = -1;
288 if (values[j] < lower[iColumn] - integerTolerance) {
289 violated = -2; // infeasible
290 break;
291 }
292 }
293 }
294 }
295
296 if (violated == -2) {
297 feasible = false ;
298 break ; // break the cut generation loop
299 }
300
301 //---------------------------------------------------------------------
302 // Now apply the row (constraint) cuts.
303 int numberRowsNow = solver_->getNumRows();
304 assert(numberRowsNow == numberRowsAtStart + lastNumberCuts);
305 int numberToAdd = theseCuts.sizeRowCuts();
306 numberNewCuts = lastNumberCuts + numberToAdd;
307
308 // Get a basis by asking the solver for warm start information.
309 // Resize it (retaining the basis) so it can accommodate the cuts.
310 delete basis_;
311 basis_ = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart());
312 assert(basis_ != NULL); // make sure not volume
313 basis_->resize(numberRowsAtStart + numberNewCuts, numberColumns);
314
315 // Now actually add the row cuts and reoptimise.
316 if (numberRowCuts > 0 || numberColumnCuts > 0) {
317 if (numberToAdd > 0) {
318 int i;
319 OsiRowCut * addCuts = new OsiRowCut [numberToAdd];
320 for (i = 0; i < numberToAdd; ++i) {
321 addCuts[i] = theseCuts.rowCut(i);
322 }
323 solver_->applyRowCuts(numberToAdd, addCuts);
324 // AJK this caused a memory fault on Win32
325 delete [] addCuts;
326 for (i = 0; i < numberToAdd; ++i) {
327 cuts.insert(theseCuts.rowCut(i));
328 }
329 for (i = 0; i < numberToAdd; ++i) {
330 basis_->setArtifStatus(numberRowsNow + i,
331 CoinWarmStartBasis::basic);
332 }
333 if (solver_->setWarmStart(basis_) == false) {
334 throw CoinError("Fail setWarmStart() after cut install.",
335 "solveWithCuts", "SbbModel");
336 }
337 }
338 feasible = resolve() ;
339 }
340 else {
341 numberTries = 0;
342 }
343
344 //---------------------------------------------------------------------
345 if (feasible) {
346 int cutIterations = solver_->getIterationCount();
347 //takeOffCuts(cuts, whichGenerator,
348 // numberOldActiveCuts, numberNewCuts, true);
349 if (solver_->isDualObjectiveLimitReached()) {
350 feasible = false;
351 #ifdef ABC_DEBUG
352 double z = solver_->getObjValue();
353 double cut = getCutoff();
354 // printf("Lost feasibility by %g in takeOffCuts; z = %g, cutoff = %g\n",
355 // z - cut, z, cut);
356 #endif
357 }
358 if (feasible) {
359 numberRowsAtStart = numberOldActiveCuts +
360 numberRowsAtContinuous_;
361 lastNumberCuts = numberNewCuts;
362 if ((direction * solver_->getObjValue() <
363 lastObjective + minimumDrop) && (numberPasses >= 3)) {
364 numberTries = 0;
365 }
366 if (numberRowCuts+numberColumnCuts == 0 || cutIterations == 0)
367 { break; }
368 if (numberTries > 0) {
369 reducedCostFix();
370 lastObjective = direction * solver_->getObjValue();
371 lower = solver_->getColLower();
372 upper = solver_->getColUpper();
373 solution = solver_->getColSolution();
374 }
375 }
376 }
377
378 // We've lost feasibility
379 if (!feasible) {
380 numberTries = 0;
381 }
382 } while (numberTries);
383 // END OF GENERATING CUTS
384
385 //------------------------------------------------------------------------
386 // Adjust the frequency of use for any of the cut generator
387 double thisObjective = solver_->getObjValue() * direction;
388 if (feasible && fullScan && numberCutGenerators_) {
389 double totalCuts = 0.0;
390 int i;
391 for (int i = 0; i < numberCutGenerators_; ++i)
392 totalCuts += countRowCuts[i] + 5.0 * countColumnCuts[i];
393 // Root node or every so often - see what to turn off
394 if (!numberNodes_)
395 handler_->message(ABC_ROOT, messages_)
396 << numberNewCuts
397 << startObjective << thisObjective
398 << numberPasses
399 << CoinMessageEol;
400 int * count = new int[numberCutGenerators_];
401 memset(count, 0, numberCutGenerators_ * sizeof(int));
402 for (i = 0; i < numberNewCuts; ++i)
403 count[whichGenerator[i]]++;
404 double small = (0.5 * totalCuts) / ((double) numberCutGenerators_);
405 for (i = 0; i < numberCutGenerators_; ++i) {
406 int howOften = generator_[i]->howOften();
407 if (howOften < -99)
408 continue;
409 if (howOften < 0 || howOften >= 1000000) {
410 // If small number switch mostly off
411 double thisCuts = countRowCuts[i] + 5.0 * countColumnCuts[i];
412 if (!thisCuts || howOften == -99) {
413 if (howOften == -99)
414 howOften = -100;
415 else
416 howOften = 1000000 + SCANCUTS; // wait until next time
417 }
418 else if (thisCuts < small) {
419 int k = (int) sqrt(small / thisCuts);
420 howOften = k + 1000000;
421 }
422 else {
423 howOften = 1 + 1000000;
424 }
425 }
426 generator_[i]->setHowOften(howOften);
427 int newFrequency = generator_[i]->howOften() % 1000000;
428 // if (handler_->logLevel() > 1 || !numberNodes_)
429 if (!numberNodes_)
430 handler_->message(ABC_GENERATOR, messages_)
431 << i
432 << generator_[i]->cutGeneratorName()
433 << countRowCuts[i]
434 << countRowCuts[i] //<<count[i]
435 << countColumnCuts[i]
436 << newFrequency
437 << CoinMessageEol;
438 }
439 delete [] count;
440 }
441
442 delete [] countRowCuts;
443 delete [] countColumnCuts;
444
445 #ifdef CHECK_CUT_COUNTS
446 if (feasible) {
447 delete basis_;
448 basis_ = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart());
449 printf("solveWithCuts: Number of rows at end (only active cuts) %d\n",
450 numberRowsAtContinuous_+numberNewCuts+numberOldActiveCuts);
451 basis_->print();
452 }
453 if (numberNodes_ % 1000 == 0) {
454 messageHandler()->message(ABC_CUTS, messages_)
455 << numberNodes_
456 << numberNewCuts
457 << startObjective
458 << thisObjective
459 << numberPasses
460 << CoinMessageEol;
461 }
462 #endif
463
464
465 //takeOffCuts(cuts, whichGenerator, numberOldActiveCuts,
466 // numberNewCuts, true);
467 incrementNodeCount();
468
469 return feasible;
470 }
471
472 //#############################################################################
473
474 bool
resolve()475 AbcModel::resolve()
476 {
477 int iRow;
478 int numberRows = solver_->getNumRows();
479 const double * rowLower = solver_->getRowLower();
480 const double * rowUpper = solver_->getRowUpper();
481 bool feasible = true;
482 for (iRow = numberRowsAtContinuous_; iRow < numberRows; ++iRow) {
483 if (rowLower[iRow] > rowUpper[iRow] + 1.0e-8)
484 feasible = false;
485 }
486
487 // Reoptimize. Consider the possibility that we should fathom on bounds.
488 // But be careful --- where the objective takes on integral values, we may
489 // want to keep a solution where the objective is right on the cutoff.
490 if (feasible) {
491 solver_->resolve();
492 numberIterations_ += getIterationCount();
493 feasible = (solver_->isProvenOptimal() &&
494 !solver_->isDualObjectiveLimitReached());
495 }
496
497 return feasible;
498 }
499
500 //#############################################################################
501
502 double
checkSolution(double cutoff,const double * solution,bool fixVariables)503 AbcModel::checkSolution (double cutoff,
504 const double *solution,
505 bool fixVariables)
506 {
507 return 0.0;
508 }
509
510 //#############################################################################
511
512 bool
setBestSolution(ABC_Message how,double & objectiveValue,const double * solution,bool fixVariables)513 AbcModel::setBestSolution(ABC_Message how,
514 double & objectiveValue,
515 const double * solution,
516 bool fixVariables)
517 {
518 double cutoff = getCutoff();
519 // Double check the solution to catch pretenders.
520 if (objectiveValue >= cutoff) { // Bad news
521 if (objectiveValue > 1.0e30)
522 handler_->message(ABC_NOTFEAS1, messages_) << CoinMessageEol;
523 else
524 handler_->message(ABC_NOTFEAS2, messages_)
525 << objectiveValue << cutoff << CoinMessageEol;
526 return false;
527 }
528 else { // Better solution
529 bestObjective_ = objectiveValue;
530 int numberColumns = solver_->getNumCols();
531 if (bestSolution_ == 0) {
532 bestSolution_ = new double[numberColumns];
533 }
534
535 memcpy(bestSolution_, solution, numberColumns*sizeof(double));
536 cutoff = bestObjective_ - dblParam_[AbcCutoffIncrement];
537 setCutoff(cutoff);
538
539 if (how == ABC_ROUNDING)
540 numberHeuristicSolutions_++;
541 numberSolutions_++;
542 // std::cout << "cutoff = " << getCutoff()
543 // << "; objVal = " << bestObjective_
544 // << "; cutoffInc = " << dblParam_[AbcCutoffIncrement]
545 // << std::endl;
546
547 handler_->message(how, messages_)
548 << bestObjective_ << numberIterations_
549 << numberNodes_
550 << CoinMessageEol;
551
552 return true;
553 }
554 }
555
556 //#############################################################################
557
558 bool
feasibleSolution(int & numberIntegerInfeasibilities)559 AbcModel::feasibleSolution(int & numberIntegerInfeasibilities)
560 {
561 bool feasible = true;
562 numberIntegerInfeasibilities = 0;
563 int i = -1;
564 const int numCols = getNumCols();
565
566 if (currentSolution_ != 0) {
567 delete [] currentSolution_;
568 currentSolution_ = 0;
569 }
570
571 currentSolution_ = new double [numCols];
572 memcpy(currentSolution_, solver_->getColSolution(),sizeof(double)*numCols);
573
574 for (i = 0; i < numberIntegers_; ++i) {
575 if ( ! checkInteger(currentSolution_[integerVariable_[i]]) ) {
576 ++numberIntegerInfeasibilities;
577 feasible = false;
578 }
579 }
580
581 return feasible;
582 }
583
584 //#############################################################################
585
586 void
findIntegers(bool startAgain)587 AbcModel::findIntegers(bool startAgain)
588 {
589 assert(solver_);
590
591 int iColumn;
592 int numberColumns = getNumCols();
593 const double *objCoeffs = getObjCoefficients();
594
595 if (numberStrong_ == 0) { // set up pseudocost list
596 pseudoList_ = new AbcPseudocost* [getNumCols()];
597 pseudoIndices_ = new int [getNumCols()];
598 for (iColumn = 0; iColumn < numberColumns; ++iColumn) {
599 pseudoList_[iColumn] = NULL;
600 pseudoIndices_[iColumn] = -1;
601 }
602 }
603
604 if (numberIntegers_ && !startAgain) return;
605 delete [] integerVariable_;
606 numberIntegers_ = 0;
607
608 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
609 if (isInteger(iColumn)) numberIntegers_++;
610 }
611
612 if (numberIntegers_) {
613 integerVariable_ = new int [numberIntegers_];
614 numberIntegers_=0;
615 for (iColumn = 0; iColumn < numberColumns; ++iColumn) {
616 if(isInteger(iColumn)) {
617 integerVariable_[numberIntegers_++] = iColumn;
618 if (numberStrong_ == 0) {
619 double obj = fabs(objCoeffs[iColumn]);
620 AbcPseudocost *pcost = new AbcPseudocost(iColumn,
621 obj,
622 0,
623 obj,
624 0);
625 pseudoList_[iColumn] = pcost;
626 //printf("numberIntegers_ = %d\n", numberIntegers_);
627 }
628 //printf("out\n");
629 }
630 }
631 }
632 else {
633 handler_->message(ABC_NOINT, messages_) << CoinMessageEol ;
634 }
635
636
637
638 }
639
640 //#############################################################################
641 // Add one generator
642 void
addCutGenerator(CglCutGenerator * generator,int howOften,const char * name,bool normal,bool atSolution,bool whenInfeasible)643 AbcModel::addCutGenerator(CglCutGenerator * generator,
644 int howOften, const char * name,
645 bool normal, bool atSolution,
646 bool whenInfeasible)
647 {
648 AbcCutGenerator ** temp = generator_;
649 generator_ = new AbcCutGenerator * [numberCutGenerators_ + 1];
650 memcpy(generator_, temp, numberCutGenerators_*sizeof(AbcCutGenerator *));
651 delete[] temp;
652 generator_[numberCutGenerators_++]=
653 new AbcCutGenerator(this, generator, howOften, name,
654 normal, atSolution, whenInfeasible);
655
656 }
657
658 //#############################################################################
659 // Add one heuristic
660 void
addHeuristic(AbcHeuristic * generator)661 AbcModel::addHeuristic(AbcHeuristic* generator)
662 {
663 AbcHeuristic ** temp = heuristic_;
664 heuristic_ = new AbcHeuristic* [numberHeuristics_ + 1];
665 memcpy(heuristic_, temp, numberHeuristics_ * sizeof(AbcHeuristic *));
666 delete [] temp;
667 heuristic_[numberHeuristics_++] = generator;
668 }
669
670 //#############################################################################
671 // Perform reduced cost fixing on integer variables. The variables in
672 // question are already nonbasic at bound. We're just nailing down the
673 // current situation.
reducedCostFix()674 void AbcModel::reducedCostFix ()
675 {
676 double cutoff = getCutoff();
677 double direction = solver_->getObjSense();
678 double gap = cutoff - solver_->getObjValue()*direction;
679 double integerTolerance = getDblParam(AbcIntegerTolerance);
680
681 const double* lower = solver_->getColLower();
682 const double* upper = solver_->getColUpper();
683 const double* solution = solver_->getColSolution();
684 const double* reducedCost = solver_->getReducedCost();
685
686 int numberFixed = 0 ;
687 for (int i = 0; i < numberIntegers_; i++) {
688 int iColumn = integerVariable_[i];
689 double djValue = direction * reducedCost[iColumn];
690 if (upper[iColumn] - lower[iColumn] > integerTolerance) {
691 if (solution[iColumn] < lower[iColumn] + integerTolerance &&
692 djValue > gap) {
693 solver_->setColUpper(iColumn, lower[iColumn]);
694 numberFixed++;
695 }
696 else if (solution[iColumn] > upper[iColumn] - integerTolerance &&
697 -djValue > gap) {
698 solver_->setColLower(iColumn, upper[iColumn]);
699 numberFixed++;
700 }
701 }
702 }
703 }
704
705 //#############################################################################
706 #if 0
707 void
708 AbcModel::takeOffCuts(OsiCuts &newCuts, int *whichGenerator,
709 int &numberOldActiveCuts, int &numberNewCuts,
710 bool allowResolve)
711 {
712 int firstOldCut = numberRowsAtContinuous_;
713 int totalNumberCuts = numberNewCuts + numberOldActiveCuts;
714 int *solverCutIndices = new int[totalNumberCuts];
715 int *newCutIndices = new int[numberNewCuts];
716 const CoinWarmStartBasis* ws;
717 CoinWarmStartBasis::Status status;
718 bool needPurge = true;
719
720 // The outer loop allows repetition of purge in the event that
721 // reoptimisation changes the basis. To start an iteration, clear the
722 // deletion counts and grab the current basis.
723
724 while (needPurge) {
725 int numberNewToDelete = 0;
726 int numberOldToDelete = 0;
727 int i;
728 ws = dynamic_cast<const CoinWarmStartBasis*>(solver_->getWarmStart());
729
730 // Scan the basis entries of the old cuts generated prior to this
731 // round of cut generation. Loose cuts are `removed'.
732 for (i = 0; i < numberOldActiveCuts; ++i) {
733 status = ws->getArtifStatus(i + firstOldCut);
734 if (status == CoinWarmStartBasis::basic) {
735 solverCutIndices[numberOldToDelete++] = i + firstOldCut;
736 }
737 }
738
739 // Scan the basis entries of the new cuts generated with this round
740 // of cut generation. At this point, newCuts is the only record of
741 // the new cuts, so when we delete loose cuts from newCuts, they're
742 // really gone. newCuts is a vector, so it's most efficient to
743 // compress it (eraseRowCut) from back to front.
744 int firstNewCut = firstOldCut + numberOldActiveCuts;
745 int k = 0;
746 for (i = 0; i < numberNewCuts; ++i) {
747 status = ws->getArtifStatus(i + firstNewCut);
748 if (status == CoinWarmStartBasis::basic) {
749 solverCutIndices[numberNewToDelete + numberOldToDelete] =
750 i + firstNewCut ;
751 newCutIndices[numberNewToDelete++] = i;
752 }
753 else { // save which generator did it
754 whichGenerator[k++] = whichGenerator[i];
755 }
756 }
757 for (i = numberNewToDelete - 1 ; i >= 0 ; i--) {
758 int iCut = newCutIndices[i];
759 newCuts.eraseRowCut(iCut);
760 }
761
762 // Did we delete anything? If so, delete the cuts from the constraint
763 // system held in the solver and reoptimise unless we're forbidden
764 // to do so. If the call to resolve() results in pivots, there's the
765 // possibility we again have basic slacks. Repeat the purging loop.
766
767 if (numberNewToDelete + numberOldToDelete > 0) {
768 solver_->deleteRows(numberNewToDelete + numberOldToDelete,
769 solverCutIndices);
770 numberNewCuts -= numberNewToDelete;
771 numberOldActiveCuts -= numberOldToDelete;
772 # ifdef ABC_DEBUG
773 std::cout << "takeOffCuts: purged " << numberOldToDelete << "+"
774 << numberNewToDelete << " cuts." << std::endl;
775 # endif
776 if (allowResolve) {
777 solver_->resolve();
778 if (solver_->getIterationCount() == 0) {
779 needPurge = false;
780 }
781 # ifdef ABC_DEBUG
782 else {
783 std::cout << "Repeating purging loop. "
784 << solver_->getIterationCount() << " iters."
785 << std::endl;
786 }
787 # endif
788 }
789 else {
790 needPurge = false;
791 }
792 }
793 else {
794 needPurge = false;
795 }
796 }
797
798 delete ws;
799 delete [] solverCutIndices;
800 delete [] newCutIndices;
801 }
802 #endif
803
804 //#############################################################################
805 #if 1
806 //void
807 //AbcModel::takeOffCuts(OsiCuts &newCuts, int *whichGenerator,
808 // int &numberOldActiveCuts, int &numberNewCuts,
809 // bool allowResolve)
810 void
takeOffCuts()811 AbcModel::takeOffCuts()
812 {
813 // assert(!numberOldActiveCuts);
814 int totalNumberCuts = solver()->getNumRows() - numberRowsAtContinuous_;
815 int *solverCutIndices = new int[totalNumberCuts];
816 // const CoinWarmStartBasis* ws;
817
818 for (int i = 0; i < totalNumberCuts; ++i) {
819 solverCutIndices[i] = i + numberRowsAtContinuous_;
820 }
821
822 // Delete all new cuts
823 solver_->deleteRows(totalNumberCuts, solverCutIndices);
824 solver_->setWarmStart(sharedBasis_);
825 //numberOldActiveCuts = numberNewCuts = 0;
826 delete [] solverCutIndices;
827 }
828 #endif
829
830 //#############################################################################
831 /**
832 This routine sets the objective cutoff value used for fathoming and
833 determining monotonic variables.
834
835 If the fathoming discipline is strict, a small tolerance is added to the
836 new cutoff. This avoids problems due to roundoff when the target value
837 is exact. The common example would be an IP with only integer variables in
838 the objective. If the target is set to the exact value z of the optimum,
839 it's possible to end up fathoming an ancestor of the solution because the
840 solver returns z+epsilon.
841
842 Determining if strict fathoming is needed is best done by analysis.
843 In sbb, that's analyseObjective. The default is false.
844
845 In sbb we always minimize so add epsilon
846 */
setCutoff(double value)847 void AbcModel::setCutoff (double value)
848 {
849 double tol = 0;
850 int fathomStrict = getIntParam(AbcFathomDiscipline);
851 double direction = solver_->getObjSense();
852 if (fathomStrict == 1) {
853 solver_->getDblParam(OsiDualTolerance, tol);
854 tol = tol * (1 + fabs(value));
855 value += tol;
856 }
857
858 // Solvers know about direction
859 solver_->setDblParam(OsiDualObjectiveLimit, value * direction);
860 }
861
862
863
createRoot()864 AlpsTreeNode * AbcModel::createRoot() {
865 return new AbcTreeNode();
866 }
867
868
869
870 //#############################################################################
871 // Initial solve and find integers
872 bool
setupSelf()873 AbcModel::setupSelf()
874 {
875 bool feasible = true;
876 solver_->messageHandler()->setLogLevel(0);
877 initialSolve();
878 sharedBasis_ = dynamic_cast<CoinWarmStartBasis*>
879 (solver_->getWarmStart());
880
881 # ifdef ABC_DEBUG_MORE
882 std::string problemName;
883 solver_->getStrParam(OsiProbName, problemName);
884 printf("Problem name - %s\n", problemName.c_str());
885 solver_->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0);
886 # endif
887
888 status_ = 0;
889
890 findIntegers(true);
891
892 bestObjective_ = 1.0e50;
893 double direction = solver_->getObjSense();
894
895 int numberColumns = getNumCols();
896 if (!currentSolution_)
897 currentSolution_ = new double[numberColumns];
898
899 //continuousSolver_ = solver_->clone();
900 numberRowsAtContinuous_ = getNumRows();
901
902 maximumNumberCuts_ = 0;
903 currentNumberCuts_ = 0;
904
905 // FIXME:
906
907 return feasible;
908 }
909
910 //#############################################################################
911 // Send model and root so that initial solve
encode(AlpsEncoded * encoded) const912 AlpsReturnStatus AbcModel::encode(AlpsEncoded * encoded) const {
913 AlpsReturnStatus status = AlpsReturnStatusOk;
914
915 //------------------------------------------------------
916 // Encode Alps part.
917 //------------------------------------------------------
918 status = AlpsModel::encode(encoded);
919
920 //------------------------------------------------------
921 // Encode Abc part.
922 //------------------------------------------------------
923
924 // Write the model data into representation_
925 const CoinPackedMatrix* matrixByCol = solver_->getMatrixByCol();
926 int numRows = getNumRows();
927 encoded->writeRep(numRows);
928 int numCols = getNumCols();
929 encoded->writeRep(numCols);
930 #if defined(ABC_DEBUG_MORE)
931 std::cout << "AbcModel::encode()-- numRows="<< numRows << "; numCols="
932 << numCols << std::endl;
933 #endif
934
935 const double* collb = solver_->getColLower();
936 encoded->writeRep(collb, numCols);
937 const double* colub = solver_->getColUpper();
938 encoded->writeRep(colub, numCols);
939 const double* obj = solver_->getObjCoefficients();
940 encoded->writeRep(obj, numCols);
941 const double objSense = solver_->getObjSense();
942 encoded->writeRep(objSense);
943 const double* rowlb = solver_->getRowLower();
944 encoded->writeRep(rowlb, numRows);
945 const double* rowub = solver_->getRowUpper();
946 encoded->writeRep(rowub, numRows);
947 int numElements = solver_->getNumElements();
948 encoded->writeRep(numElements);
949 const double* elementValue = matrixByCol->getElements();
950 encoded->writeRep(elementValue, numElements);
951 const CoinBigIndex* colStart = matrixByCol->getVectorStarts();
952 int numStart = numCols + 1;
953 encoded->writeRep(colStart, numStart);
954 const int* index = matrixByCol->getIndices();
955 encoded->writeRep(index, numElements);
956 encoded->writeRep(numberIntegers_);
957 encoded->writeRep(integerVariable_, numberIntegers_);
958 #if defined(ABC_DEBUG_MORE)
959 std::cout << "AbcModel::encode()-- objSense="<< objSense
960 << "; numElements="<< numElements
961 << "; numberIntegers_=" << numberIntegers_
962 << "; numStart = " << numStart <<std::endl;
963 #endif
964 #if defined(ABC_DEBUG_MORE)
965 std::cout << "rowub=";
966 for (int i = 0; i < numRows; ++i){
967 std::cout <<rowub[i]<<" ";
968 }
969 std::cout << std::endl;
970 std::cout << "elementValue=";
971 for (int j = 0; j < numElements; ++j) {
972 std::cout << elementValue[j] << " ";
973 }
974 std::cout << std::endl;
975 #endif
976
977 return status;
978 }
979
980 //#############################################################################
981 // Decode and load model data to LP solver.
decodeToSelf(AlpsEncoded & encoded)982 AlpsReturnStatus AbcModel::decodeToSelf(AlpsEncoded & encoded) {
983 AlpsReturnStatus status = AlpsReturnStatusOk;
984
985 //------------------------------------------------------
986 // Decode Alps part.
987 //------------------------------------------------------
988
989 status = AlpsModel::decodeToSelf(encoded);
990
991 //------------------------------------------------------
992 // Decode Abc part.
993 //------------------------------------------------------
994
995 int numRows;
996 encoded.readRep(numRows);
997 int numCols;
998 encoded.readRep(numCols);
999 #if defined(ABC_DEBUG_MORE)
1000 std::cout << "AbcModel::decode()-- numRows="<< numRows << "; numCols="
1001 << numCols << std::endl;
1002 #endif
1003 double* collb;
1004 encoded.readRep(collb, numCols);
1005 double* colub;
1006 encoded.readRep(colub, numCols);
1007 double* obj;
1008 encoded.readRep(obj, numCols);
1009 double objSense;
1010 encoded.readRep(objSense);
1011 double* rowlb;
1012 encoded.readRep(rowlb, numRows);
1013 double* rowub;
1014 encoded.readRep(rowub, numRows);
1015 int numElements;
1016 encoded.readRep(numElements);
1017 double* elementValue;
1018 encoded.readRep(elementValue, numElements);
1019 CoinBigIndex* colStart;
1020 int numStart = numCols + 1;
1021 encoded.readRep(colStart, numStart);
1022 int* index;
1023 encoded.readRep(index, numElements);
1024 encoded.readRep(numberIntegers_);
1025 encoded.readRep(integerVariable_, numberIntegers_);
1026 #if defined(ABC_DEBUG_MORE)
1027 std::cout << "AbcModel::decode()-- objSense="<< objSense
1028 << "; numElements="<< numElements
1029 << "; numberIntegers_=" << numberIntegers_
1030 << "; numStart = " << numStart <<std::endl;
1031 #endif
1032 #if defined(ABC_DEBUG_MORE)
1033 std::cout << "rowub=";
1034 for (int i = 0; i < numRows; ++i){
1035 std::cout <<rowub[i]<<" ";
1036 }
1037 std::cout << std::endl;
1038 std::cout << "elementValue=";
1039 for (int j = 0; j < numElements; ++j) {
1040 std::cout << elementValue[j] << " ";
1041 }
1042 std::cout << std::endl;
1043 std::cout << "index=";
1044 for (int j = 0; j < numElements; ++j) {
1045 std::cout << index[j] << " ";
1046 }
1047 std::cout << std::endl;
1048 std::cout << "colStart=";
1049 for (int j = 0; j < numElements+1; ++j) {
1050 std::cout << colStart[j] << " ";
1051 }
1052 std::cout << std::endl;
1053 #endif
1054
1055 // Check if solver_ is declared in main
1056 assert(solver_);
1057
1058 //-------------------------------------------------------------------------
1059 // load the standardized problem into stdSi
1060 #if 0 // load matrix doesn't work. Don't know why.
1061 CoinPackedMatrix * matrixByCol =
1062 new CoinPackedMatrix(true, numCols, numRows, numElements,
1063 elementValue, index, colStart, 0);
1064 solver_->loadProblem(*matrixByCol,
1065 collb, colub,
1066 obj,
1067 rowlb, rowub);
1068 #endif
1069 //-------------------------------------------------------------------------
1070 solver_->loadProblem(numCols, numRows,
1071 colStart, index, elementValue,
1072 collb, colub,
1073 obj,
1074 rowlb, rowub);
1075
1076 solver_->setObjSense(objSense);
1077 solver_->setInteger(integerVariable_, numberIntegers_);
1078
1079 delete [] collb;
1080 collb = NULL;
1081 delete [] colub;
1082 colub = NULL;
1083 delete [] obj;
1084 obj = NULL;
1085 delete [] rowlb;
1086 rowlb = NULL;
1087 delete [] rowub;
1088 rowub = NULL;
1089 delete [] elementValue;
1090 elementValue = NULL;
1091 delete [] colStart;
1092 colStart = NULL;
1093 delete [] index;
1094 index = NULL;
1095 return status;
1096 }
1097
1098 /// Abc does not need this for now.
decode(AlpsEncoded & encoded) const1099 AlpsKnowledge * AbcModel::decode(AlpsEncoded & encoded) const {
1100 std::cerr << "Not implemented!" << std::endl;
1101 throw std::exception();
1102 }
1103