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

Last change on this file since 1898 was 1898, checked in by stefan, 5 years ago

fixup examples

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.1 KB
Line 
1// $Id: CbcSolver3.cpp 1898 2013-04-09 18:06:04Z stefan $
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    bool takeHint;
249    OsiHintStrength strength;
250    getHintParam(OsiDoInBranchAndCut,takeHint,strength);
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      modelPtr_->setSpecialOptions(saveOptions|(64|1024));
257    } else {
258      if((specialOptions_&4)==0) 
259        modelPtr_->setSpecialOptions(saveOptions|(64|128|512|1024|4096));
260      else
261        modelPtr_->setSpecialOptions(saveOptions|(64|128|512|1024|2048|4096));
262    }
263    //printf("thin options %d size %d\n",modelPtr_->specialOptions(),modelPtr_->numberColumns());
264    setBasis(basis_,modelPtr_);
265    //modelPtr_->setLogLevel(1);
266    printf("model has %d rows\n",modelPtr_->numberRows());
267    modelPtr_->dual(0,0);
268    basis_ = getBasis(modelPtr_);
269    modelPtr_->setSpecialOptions(saveOptions);
270    if (modelPtr_->status()==0) {
271      count_++;
272      double * solution = modelPtr_->primalColumnSolution();
273      int i;
274      for (i=0;i<numberColumns;i++) {
275        if (solution[i]>1.0e-6||modelPtr_->getStatus(i)==ClpSimplex::basic) {
276          node_[i]=CoinMax(count_,node_[i]);
277          howMany_[i]++;
278        }
279      }
280    } else {
281      if (!algorithm_==2)
282        printf("infeasible early on\n");
283    }
284  } else {
285    // use counts
286    int i;
287    const double * lower = modelPtr_->columnLower();
288    const double * upper = modelPtr_->columnUpper();
289    setBasis(basis_,modelPtr_);
290    ClpSimplex * temp = new ClpSimplex(modelPtr_,nNewRow,whichRow,nNewCol,whichColumn);
291    //temp->setLogLevel(2);
292    //printf("small has %d rows and %d columns\n",nNewRow,nNewCol);
293    temp->setSpecialOptions(128+512);
294    temp->setDualObjectiveLimit(1.0e50);
295    printf("model has %d rows\n",nNewRow);
296    temp->dual();
297    if (temp->status()) {
298      // In some cases we know that it must be infeasible
299      if (believeInfeasible_||algorithm_==1) {
300        modelPtr_->setProblemStatus(1);
301        printf("assuming infeasible!\n");
302        //modelPtr_->writeMps("infeas.mps");
303        //temp->writeMps("infeas2.mps");
304        //abort();
305        delete temp;
306        delete [] whichRow;
307        delete [] whichColumn;
308        return;
309      }
310    }
311    double * solution = modelPtr_->primalColumnSolution();
312    if (!temp->status()) {
313      const double * solution2 = temp->primalColumnSolution();
314      memset(solution,0,numberColumns*sizeof(double));
315      for (i=0;i<nNewCol;i++) {
316        int iColumn = whichColumn[i];
317        solution[iColumn]=solution2[i];
318        modelPtr_->setStatus(iColumn,temp->getStatus(i));
319      }
320      double * rowSolution = modelPtr_->primalRowSolution();
321      const double * rowSolution2 = temp->primalRowSolution();
322      double * dual = modelPtr_->dualRowSolution();
323      const double * dual2 = temp->dualRowSolution();
324      memset(dual,0,numberRows*sizeof(double));
325      for (i=0;i<nNewRow;i++) {
326        int iRow=whichRow[i];
327        modelPtr_->setRowStatus(iRow,temp->getRowStatus(i));
328        rowSolution[iRow]=rowSolution2[i];
329        dual[iRow]=dual2[i];
330      }
331      // See if optimal
332      double * dj = modelPtr_->dualColumnSolution();
333      // get reduced cost for large problem
334      // this assumes minimization
335      memcpy(dj,modelPtr_->objective(),numberColumns*sizeof(double));
336      modelPtr_->transposeTimes(-1.0,dual,dj);
337      modelPtr_->setObjectiveValue(temp->objectiveValue());
338      modelPtr_->setProblemStatus(0);
339      int nBad=0;
340     
341      for (i=0;i<numberColumns;i++) {
342        if (modelPtr_->getStatus(i)==ClpSimplex::atLowerBound
343            &&upper[i]>lower[i]&&dj[i]<-1.0e-5)
344          nBad++;
345      }
346      //modelPtr_->writeMps("bada.mps");
347      //temp->writeMps("badb.mps");
348      if (nBad) {
349        assert (algorithm_==2);
350        //printf("%d bad\n",nBad);
351        timesBad_++;
352        modelPtr_->primal();
353      }
354    } else {
355      // infeasible - do all
356      modelPtr_->setSpecialOptions(64+128+512);
357      setBasis(basis_,modelPtr_);
358      //modelPtr_->setLogLevel(1);
359      modelPtr_->dual(0,0);
360      basis_ = getBasis(modelPtr_);
361      modelPtr_->setSpecialOptions(0);
362      if (modelPtr_->status()) {
363        printf("really infeasible!\n");
364        delete temp;
365        delete [] whichRow;
366        delete [] whichColumn;
367        return;
368      } else {
369        printf("initially infeasible\n");
370      }
371    }
372    delete temp;
373    delete [] whichRow;
374    delete [] whichColumn;
375    basis_ = getBasis(modelPtr_);
376    modelPtr_->setSpecialOptions(0);
377    count_++;
378    if ((count_%100)==0&&algorithm_==2)
379      printf("count %d, bad %d\n",count_,timesBad_);
380    for (i=0;i<numberColumns;i++) {
381      if (solution[i]>1.0e-6||modelPtr_->getStatus(i)==ClpSimplex::basic) {
382        node_[i]=CoinMax(count_,node_[i]);
383        howMany_[i]++;
384      }
385    }
386    if (modelPtr_->objectiveValue()>=modelPtr_->dualObjectiveLimit())
387      modelPtr_->setProblemStatus(1);
388  }
389}
390
391//#############################################################################
392// Constructors, destructors clone and assignment
393//#############################################################################
394
395//-------------------------------------------------------------------
396// Default Constructor
397//-------------------------------------------------------------------
398CbcSolver3::CbcSolver3 ()
399  : OsiClpSolverInterface()
400{
401  node_=NULL;
402  howMany_=NULL;
403  count_=0;
404  model_ = NULL;
405  memory_=300;
406  believeInfeasible_=false;
407  nestedSearch_ = 1.0;
408  algorithm_=0;
409}
410
411//-------------------------------------------------------------------
412// Clone
413//-------------------------------------------------------------------
414OsiSolverInterface * 
415CbcSolver3::clone(bool CopyData) const
416{
417  if (CopyData) {
418    return new CbcSolver3(*this);
419  } else {
420    printf("warning CbcSolveUser clone with copyData false\n");
421    return new CbcSolver3();
422  }
423}
424
425
426//-------------------------------------------------------------------
427// Copy constructor
428//-------------------------------------------------------------------
429CbcSolver3::CbcSolver3 (
430                  const CbcSolver3 & rhs)
431  : OsiClpSolverInterface(rhs)
432{
433  model_ = rhs.model_;
434  int numberColumns = modelPtr_->numberColumns();
435  node_=CoinCopyOfArray(rhs.node_,numberColumns);
436  howMany_=CoinCopyOfArray(rhs.howMany_,numberColumns);
437  count_=rhs.count_;
438  memory_=rhs.memory_;
439  believeInfeasible_ = rhs.believeInfeasible_;
440  nestedSearch_ = rhs.nestedSearch_;
441  algorithm_=rhs.algorithm_;
442}
443
444//-------------------------------------------------------------------
445// Destructor
446//-------------------------------------------------------------------
447CbcSolver3::~CbcSolver3 ()
448{
449  delete [] node_;
450  delete [] howMany_;
451}
452
453//-------------------------------------------------------------------
454// Assignment operator
455//-------------------------------------------------------------------
456CbcSolver3 &
457CbcSolver3::operator=(const CbcSolver3& rhs)
458{
459  if (this != &rhs) { 
460    OsiClpSolverInterface::operator=(rhs);
461    delete [] node_;
462    delete [] howMany_;
463    model_ = rhs.model_;
464    int numberColumns = modelPtr_->numberColumns();
465    node_=CoinCopyOfArray(rhs.node_,numberColumns);
466    howMany_=CoinCopyOfArray(rhs.howMany_,numberColumns);
467    count_=rhs.count_;
468    memory_=rhs.memory_;
469    believeInfeasible_ = rhs.believeInfeasible_;
470    nestedSearch_ = rhs.nestedSearch_;
471    algorithm_=rhs.algorithm_;
472  }
473  return *this;
474}
475//-------------------------------------------------------------------
476// Real initializer
477//-------------------------------------------------------------------
478void
479CbcSolver3::initialize (CbcModel * model, const char * keep)
480{
481  model_=model;
482  int numberColumns = modelPtr_->numberColumns();
483  if (numberColumns) {
484    node_ = new int[numberColumns];
485    howMany_ = new int[numberColumns];
486    for (int i=0;i<numberColumns;i++) {
487      if (keep&&keep[i])
488        node_[i]=COIN_INT_MAX;
489      else
490        node_[i]=0;
491      howMany_[i]=0;
492    }
493  } else {
494    node_=NULL;
495    howMany_=NULL;
496  }
497}
Note: See TracBrowser for help on using the repository browser.