source: trunk/Cbc/examples/CbcSolver3.cpp @ 2469

Last change on this file since 2469 was 2469, checked in by unxusr, 9 months ago

formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.7 KB
Line 
1// $Id: CbcSolver3.cpp 2469 2019-01-06 23:17:46Z unxusr $
2// Copyright (C) 2005, 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 <cassert>
7
8#include "CoinTime.hpp"
9
10#include "CoinHelperFunctions.hpp"
11#include "CoinIndexedVector.hpp"
12#include "ClpFactorization.hpp"
13#include "ClpObjective.hpp"
14#include "ClpSimplex.hpp"
15#include "CbcSolver3.hpp"
16#include "CbcModel.hpp"
17#include "ClpPresolve.hpp"
18#include "CbcHeuristicGreedy.hpp"
19#include "CbcBranchActual.hpp"
20#include "CbcCutGenerator.hpp"
21#include "CbcCompareUser.hpp"
22#include "CbcCompareActual.hpp"
23#include "CbcCompareObjective.hpp"
24// Cuts
25
26#include "CglGomory.hpp"
27#include "CglProbing.hpp"
28#include "CglKnapsackCover.hpp"
29#include "CglOddHole.hpp"
30#include "CglClique.hpp"
31#include "CglFlowCover.hpp"
32#include "CglMixedIntegerRounding.hpp"
33#include "CglTwomir.hpp"
34
35static int timesBad_ = 0;
36//#############################################################################
37// Solve methods
38//#############################################################################
39void CbcSolver3::initialSolve()
40{
41  modelPtr_->scaling(0);
42  setBasis(basis_, modelPtr_);
43  // Do long thin by sprint
44  ClpSolve options;
45  options.setSolveType(ClpSolve::usePrimalorSprint);
46  options.setPresolveType(ClpSolve::presolveOff);
47  options.setSpecialOption(1, 3, 30);
48  modelPtr_->initialSolve(options);
49  basis_ = getBasis(modelPtr_);
50  modelPtr_->setLogLevel(0);
51}
52
53//-----------------------------------------------------------------------------
54void CbcSolver3::resolve()
55{
56  int *whichRow = NULL;
57  int *whichColumn = NULL;
58  // problem may be small enough to do nested search
59  const double *colLower = modelPtr_->columnLower();
60  const double *colUpper = modelPtr_->columnUpper();
61
62  int numberIntegers = model_->numberIntegers();
63  const int *integerVariable = model_->integerVariable();
64  int numberRows = modelPtr_->numberRows();
65  int numberColumns = modelPtr_->numberColumns();
66
67  int i;
68  int nFix = 0;
69  int nNewRow = 0;
70  int nNewCol = 0;
71  int smallOriginalNumberRows = 0;
72  if (algorithm_ == 0) {
73    for (i = 0; i < numberIntegers; i++) {
74      int iColumn = integerVariable[i];
75      if (colLower[iColumn] == colUpper[iColumn])
76        nFix++;
77    }
78  } else {
79    whichRow = new int[numberRows];
80    whichColumn = new int[numberColumns];
81    for (i = 0; i < numberRows; i++)
82      whichRow[i] = i;
83    nNewRow = numberRows;
84    smallOriginalNumberRows = nNewRow;
85    for (i = 0; i < numberIntegers; i++) {
86      int iColumn = integerVariable[i];
87      if (colLower[iColumn] == colUpper[iColumn])
88        nFix++;
89      bool choose;
90      if (algorithm_ == 1)
91        choose = true;
92      else
93        choose = (node_[i] > count_ - memory_ && node_[i] > 0);
94      if ((choose && colUpper[i])
95        || (modelPtr_->getStatus(i) != ClpSimplex::atLowerBound && modelPtr_->getStatus(i) != ClpSimplex::isFixed)
96        || colLower[i] > 0.0)
97        whichColumn[nNewCol++] = i;
98    }
99  }
100  if (nestedSearch_ < 1.0 && model_ && model_->phase() == 2) {
101    if (nFix > nestedSearch_ * numberIntegers) {
102      // Do nested search
103      // back to original number of rows
104      nNewRow = smallOriginalNumberRows;
105      // and get rid of any basics
106      int nNewCol = 0;
107      for (i = 0; i < numberColumns; i++) {
108        if (colUpper[i] || colLower[i] > 0.0)
109          whichColumn[nNewCol++] = i;
110      }
111      // start again very simply
112      ClpSimplex temp(modelPtr_, nNewRow, whichRow, nNewCol, whichColumn);
113      int returnCode;
114      OsiClpSolverInterface temp2(&temp);
115      temp2.setupForRepeatedUse(2);
116      int numberColumns2 = temp.numberColumns();
117      const double *colUpper2 = temp2.getColUpper();
118      const double *colLower2 = temp2.getColLower();
119      const double *solution2 = temp.getColSolution();
120      double *cleanSolution2 = new double[numberColumns2];
121      for (i = 0; i < numberColumns2; i++) {
122        temp2.setInteger(i);
123        double value = solution2[i];
124        value = CoinMin(CoinMax(value, colLower2[i]), colUpper2[i]);
125        cleanSolution2[i] = value;
126      }
127      temp2.setColSolution(cleanSolution2);
128      delete[] cleanSolution2;
129      CbcModel modelSmall(temp2);
130      modelSmall.setNumberStrong(0);
131      CglProbing generator1;
132      generator1.setUsingObjective(true);
133      generator1.setMaxPass(3);
134      generator1.setMaxProbe(100);
135      generator1.setMaxLook(50);
136      generator1.setRowCuts(3);
137
138      CglGomory generator2;
139      // try larger limit
140      generator2.setLimit(3000);
141
142      CglKnapsackCover generator3;
143
144      CglOddHole generator4;
145      generator4.setMinimumViolation(0.005);
146      generator4.setMinimumViolationPer(0.00002);
147      // try larger limit
148      generator4.setMaximumEntries(200);
149
150      CglClique generator5;
151      generator5.setStarCliqueReport(false);
152      generator5.setRowCliqueReport(false);
153
154      CglMixedIntegerRounding mixedGen;
155      CglFlowCover flowGen;
156      // Add in generators (just at root)
157      modelSmall.addCutGenerator(&generator1, -99, "Probing", true, false, false, -1);
158      modelSmall.addCutGenerator(&generator2, -99, "Gomory", true, false, false, -99);
159      modelSmall.addCutGenerator(&generator3, -99, "Knapsack", true, false, false, -99);
160      modelSmall.addCutGenerator(&generator4, -99, "OddHole", true, false, false, -99);
161      modelSmall.addCutGenerator(&generator5, -99, "Clique", true, false, false, -99);
162      modelSmall.addCutGenerator(&flowGen, -99, "FlowCover", true, false, false, -99);
163      modelSmall.addCutGenerator(&mixedGen, -99, "MixedIntegerRounding", true, false, false, -100);
164#if 1
165      const CoinPackedMatrix *matrix = temp2.getMatrixByCol();
166      const int *columnLength = matrix->getVectorLengths();
167      int *priority = new int[numberColumns2 + 1];
168      // do pseudo costs and priorities - take a reasonable guess
169      CbcObject **objects = new CbcObject *[numberColumns2 + 1];
170      int n = 0;
171      const double *objective = modelSmall.getObjCoefficients();
172      for (i = 0; i < numberColumns2; i++) {
173        CbcSimpleIntegerPseudoCost *newObject = new CbcSimpleIntegerPseudoCost(&modelSmall, n, i, objective[i], 0.5 * objective[i]);
174        newObject->setMethod(3);
175        objects[n] = newObject;
176        priority[n++] = 10000 - columnLength[i];
177      }
178      modelSmall.addObjects(n, objects);
179      for (i = 0; i < n; i++)
180        delete objects[i];
181      delete[] objects;
182      modelSmall.passInPriorities(priority, false);
183      delete[] priority;
184#endif
185      modelSmall.setCutoff(model_->getCutoff());
186      modelSmall.messageHandler()->setLogLevel(1);
187      modelSmall.solver()->messageHandler()->setLogLevel(0);
188      modelSmall.messagesPointer()->setDetailMessage(3, 9);
189      modelSmall.messagesPointer()->setDetailMessage(3, 6);
190      modelSmall.messagesPointer()->setDetailMessage(3, 4);
191      modelSmall.messagesPointer()->setDetailMessage(3, 13);
192      modelSmall.messagesPointer()->setDetailMessage(3, 14);
193      modelSmall.messagesPointer()->setDetailMessage(3, 1);
194      modelSmall.messagesPointer()->setDetailMessage(3, 3007);
195      // Definition of node choice
196      CbcCompareObjective compare;
197      modelSmall.setNodeComparison(compare);
198      // And Greedy heuristic
199      CbcHeuristicGreedyCover heuristic2(modelSmall);
200      // Use original upper and perturb more
201      heuristic2.setAlgorithm(11);
202      modelSmall.addHeuristic(&heuristic2);
203      modelSmall.branchAndBound();
204      temp2.releaseClp();
205      if (modelSmall.bestSolution()) {
206        double objValue = 0.0;
207        const double *solution2 = modelSmall.bestSolution();
208        double *solution = modelPtr_->primalColumnSolution();
209        const double *objective = modelPtr_->objective();
210        for (i = 0; i < numberColumns; i++)
211          solution[i] = colLower[i];
212        for (i = 0; i < nNewCol; i++) {
213          int iColumn = whichColumn[i];
214          solution[iColumn] = solution2[i];
215        }
216        for (i = 0; i < numberColumns; i++)
217          objValue += solution[i] * objective[i];
218        assert(objValue < model_->getCutoff());
219        if (objValue < model_->getCutoff()) {
220          //printf("good solution \n");
221          model_->setBestSolution(CBC_TREE_SOL, objValue, solution);
222          returnCode = 0;
223        } else {
224          returnCode = 2;
225        }
226      } else {
227        returnCode = 2;
228      }
229      if (returnCode != 0 && returnCode != 2) {
230        printf("pretending entire search done\n");
231        returnCode = 0;
232      }
233      if (returnCode == 0 || returnCode == 2) {
234        modelPtr_->setProblemStatus(1);
235        delete[] whichRow;
236        delete[] whichColumn;
237        return;
238      }
239    }
240  }
241  if ((count_ < 100 && algorithm_ == 2) || !algorithm_) {
242    delete[] whichRow;
243    delete[] whichColumn;
244    assert(!modelPtr_->specialOptions());
245    int saveOptions = modelPtr_->specialOptions();
246    bool takeHint;
247    OsiHintStrength strength;
248    getHintParam(OsiDoInBranchAndCut, takeHint, strength);
249    if (strength != OsiHintIgnore && takeHint) {
250      // could do something - think about it
251      //printf("thin hint %d %c\n",strength,takeHint ? 'T' :'F');
252    }
253    if ((specialOptions_ & 1) == 0) {
254      modelPtr_->setSpecialOptions(saveOptions | (64 | 1024));
255    } else {
256      if ((specialOptions_ & 4) == 0)
257        modelPtr_->setSpecialOptions(saveOptions | (64 | 128 | 512 | 1024 | 4096));
258      else
259        modelPtr_->setSpecialOptions(saveOptions | (64 | 128 | 512 | 1024 | 2048 | 4096));
260    }
261    //printf("thin options %d size %d\n",modelPtr_->specialOptions(),modelPtr_->numberColumns());
262    setBasis(basis_, modelPtr_);
263    //modelPtr_->setLogLevel(1);
264    printf("model has %d rows\n", modelPtr_->numberRows());
265    modelPtr_->dual(0, 0);
266    basis_ = getBasis(modelPtr_);
267    modelPtr_->setSpecialOptions(saveOptions);
268    if (modelPtr_->status() == 0) {
269      count_++;
270      double *solution = modelPtr_->primalColumnSolution();
271      int i;
272      for (i = 0; i < numberColumns; i++) {
273        if (solution[i] > 1.0e-6 || modelPtr_->getStatus(i) == ClpSimplex::basic) {
274          node_[i] = CoinMax(count_, node_[i]);
275          howMany_[i]++;
276        }
277      }
278    } else {
279      if (!algorithm_ == 2)
280        printf("infeasible early on\n");
281    }
282  } else {
283    // use counts
284    int i;
285    const double *lower = modelPtr_->columnLower();
286    const double *upper = modelPtr_->columnUpper();
287    setBasis(basis_, modelPtr_);
288    ClpSimplex *temp = new ClpSimplex(modelPtr_, nNewRow, whichRow, nNewCol, whichColumn);
289    //temp->setLogLevel(2);
290    //printf("small has %d rows and %d columns\n",nNewRow,nNewCol);
291    temp->setSpecialOptions(128 + 512);
292    temp->setDualObjectiveLimit(1.0e50);
293    printf("model has %d rows\n", nNewRow);
294    temp->dual();
295    if (temp->status()) {
296      // In some cases we know that it must be infeasible
297      if (believeInfeasible_ || algorithm_ == 1) {
298        modelPtr_->setProblemStatus(1);
299        printf("assuming infeasible!\n");
300        //modelPtr_->writeMps("infeas.mps");
301        //temp->writeMps("infeas2.mps");
302        //abort();
303        delete temp;
304        delete[] whichRow;
305        delete[] whichColumn;
306        return;
307      }
308    }
309    double *solution = modelPtr_->primalColumnSolution();
310    if (!temp->status()) {
311      const double *solution2 = temp->primalColumnSolution();
312      memset(solution, 0, numberColumns * sizeof(double));
313      for (i = 0; i < nNewCol; i++) {
314        int iColumn = whichColumn[i];
315        solution[iColumn] = solution2[i];
316        modelPtr_->setStatus(iColumn, temp->getStatus(i));
317      }
318      double *rowSolution = modelPtr_->primalRowSolution();
319      const double *rowSolution2 = temp->primalRowSolution();
320      double *dual = modelPtr_->dualRowSolution();
321      const double *dual2 = temp->dualRowSolution();
322      memset(dual, 0, numberRows * sizeof(double));
323      for (i = 0; i < nNewRow; i++) {
324        int iRow = whichRow[i];
325        modelPtr_->setRowStatus(iRow, temp->getRowStatus(i));
326        rowSolution[iRow] = rowSolution2[i];
327        dual[iRow] = dual2[i];
328      }
329      // See if optimal
330      double *dj = modelPtr_->dualColumnSolution();
331      // get reduced cost for large problem
332      // this assumes minimization
333      memcpy(dj, modelPtr_->objective(), numberColumns * sizeof(double));
334      modelPtr_->transposeTimes(-1.0, dual, dj);
335      modelPtr_->setObjectiveValue(temp->objectiveValue());
336      modelPtr_->setProblemStatus(0);
337      int nBad = 0;
338
339      for (i = 0; i < numberColumns; i++) {
340        if (modelPtr_->getStatus(i) == ClpSimplex::atLowerBound
341          && upper[i] > lower[i] && dj[i] < -1.0e-5)
342          nBad++;
343      }
344      //modelPtr_->writeMps("bada.mps");
345      //temp->writeMps("badb.mps");
346      if (nBad) {
347        assert(algorithm_ == 2);
348        //printf("%d bad\n",nBad);
349        timesBad_++;
350        modelPtr_->primal();
351      }
352    } else {
353      // infeasible - do all
354      modelPtr_->setSpecialOptions(64 + 128 + 512);
355      setBasis(basis_, modelPtr_);
356      //modelPtr_->setLogLevel(1);
357      modelPtr_->dual(0, 0);
358      basis_ = getBasis(modelPtr_);
359      modelPtr_->setSpecialOptions(0);
360      if (modelPtr_->status()) {
361        printf("really infeasible!\n");
362        delete temp;
363        delete[] whichRow;
364        delete[] whichColumn;
365        return;
366      } else {
367        printf("initially infeasible\n");
368      }
369    }
370    delete temp;
371    delete[] whichRow;
372    delete[] whichColumn;
373    basis_ = getBasis(modelPtr_);
374    modelPtr_->setSpecialOptions(0);
375    count_++;
376    if ((count_ % 100) == 0 && algorithm_ == 2)
377      printf("count %d, bad %d\n", count_, timesBad_);
378    for (i = 0; i < numberColumns; i++) {
379      if (solution[i] > 1.0e-6 || modelPtr_->getStatus(i) == ClpSimplex::basic) {
380        node_[i] = CoinMax(count_, node_[i]);
381        howMany_[i]++;
382      }
383    }
384    if (modelPtr_->objectiveValue() >= modelPtr_->dualObjectiveLimit())
385      modelPtr_->setProblemStatus(1);
386  }
387}
388
389//#############################################################################
390// Constructors, destructors clone and assignment
391//#############################################################################
392
393//-------------------------------------------------------------------
394// Default Constructor
395//-------------------------------------------------------------------
396CbcSolver3::CbcSolver3()
397  : OsiClpSolverInterface()
398{
399  node_ = NULL;
400  howMany_ = NULL;
401  count_ = 0;
402  model_ = NULL;
403  memory_ = 300;
404  believeInfeasible_ = false;
405  nestedSearch_ = 1.0;
406  algorithm_ = 0;
407}
408
409//-------------------------------------------------------------------
410// Clone
411//-------------------------------------------------------------------
412OsiSolverInterface *
413CbcSolver3::clone(bool CopyData) const
414{
415  if (CopyData) {
416    return new CbcSolver3(*this);
417  } else {
418    printf("warning CbcSolveUser clone with copyData false\n");
419    return new CbcSolver3();
420  }
421}
422
423//-------------------------------------------------------------------
424// Copy constructor
425//-------------------------------------------------------------------
426CbcSolver3::CbcSolver3(
427  const CbcSolver3 &rhs)
428  : OsiClpSolverInterface(rhs)
429{
430  model_ = rhs.model_;
431  int numberColumns = modelPtr_->numberColumns();
432  node_ = CoinCopyOfArray(rhs.node_, numberColumns);
433  howMany_ = CoinCopyOfArray(rhs.howMany_, numberColumns);
434  count_ = rhs.count_;
435  memory_ = rhs.memory_;
436  believeInfeasible_ = rhs.believeInfeasible_;
437  nestedSearch_ = rhs.nestedSearch_;
438  algorithm_ = rhs.algorithm_;
439}
440
441//-------------------------------------------------------------------
442// Destructor
443//-------------------------------------------------------------------
444CbcSolver3::~CbcSolver3()
445{
446  delete[] node_;
447  delete[] howMany_;
448}
449
450//-------------------------------------------------------------------
451// Assignment operator
452//-------------------------------------------------------------------
453CbcSolver3 &
454CbcSolver3::operator=(const CbcSolver3 &rhs)
455{
456  if (this != &rhs) {
457    OsiClpSolverInterface::operator=(rhs);
458    delete[] node_;
459    delete[] howMany_;
460    model_ = rhs.model_;
461    int numberColumns = modelPtr_->numberColumns();
462    node_ = CoinCopyOfArray(rhs.node_, numberColumns);
463    howMany_ = CoinCopyOfArray(rhs.howMany_, numberColumns);
464    count_ = rhs.count_;
465    memory_ = rhs.memory_;
466    believeInfeasible_ = rhs.believeInfeasible_;
467    nestedSearch_ = rhs.nestedSearch_;
468    algorithm_ = rhs.algorithm_;
469  }
470  return *this;
471}
472//-------------------------------------------------------------------
473// Real initializer
474//-------------------------------------------------------------------
475void CbcSolver3::initialize(CbcModel *model, const char *keep)
476{
477  model_ = model;
478  int numberColumns = modelPtr_->numberColumns();
479  if (numberColumns) {
480    node_ = new int[numberColumns];
481    howMany_ = new int[numberColumns];
482    for (int i = 0; i < numberColumns; i++) {
483      if (keep && keep[i])
484        node_[i] = COIN_INT_MAX;
485      else
486        node_[i] = 0;
487      howMany_[i] = 0;
488    }
489  } else {
490    node_ = NULL;
491    howMany_ = NULL;
492  }
493}
Note: See TracBrowser for help on using the repository browser.