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

Last change on this file since 1464 was 1464, checked in by stefan, 9 years ago

merge split branch into trunk; fix some examples

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