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

Last change on this file since 1574 was 1574, checked in by lou, 8 years ago

Change to EPL license notice.

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