source: branches/devel/Cbc/examples/CbcSolver3.cpp @ 425

Last change on this file since 425 was 174, checked in by forrest, 14 years ago

new includes

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